郭孝城，周昱成，王 赤，李 暉

（1.中國科學(xué)院 國家空間科學(xué)中心 空間天氣學(xué)國家重點實驗室，北京 100190；2.中國科學(xué)院大學(xué) 地球和行星科學(xué)學(xué)院，北京 100049）

引 言

源于太陽表面的超聲速太陽風(fēng)等離子體向外徑向膨脹時會與包含星際等離子體和中性原子的鄰近星際介質(zhì)相互作用形成一泡狀結(jié)構(gòu)即日球?qū)樱鼈兊慕唤缑婢褪侨涨驅(qū)禹擺1]。在經(jīng)歷漫長的旅程后，“旅行者1號”（ Voyager 1）和“旅行者2號”（ Voyager 2）飛船分別在距太陽94 AU和84 AU左右穿越終止激波，最終在120 AU左右穿越日球?qū)禹斶M(jìn)入星際空間[2-3]。太陽風(fēng)事件，包括共轉(zhuǎn)太陽風(fēng)作用區(qū)（Corotating Interaction Regions，CIRs）[4]和行星際日冕拋射事件（Interplanetary Coronal Mass Ejections，ICMEs）[5]，在日球?qū)觾?nèi)傳輸時會引起系列反應(yīng)，比如飛船發(fā)現(xiàn)在幾十AU處太陽風(fēng)事件傳輸時會逐漸匯聚合并形成更大尺度的等離子體脈沖或激波事件[6-8]

“旅行者1號”在2012年8月穿越日球?qū)禹敽?，磁強計? 012.92（十進(jìn)制年份）和2014.66年（十進(jìn)制年份）分別探測到兩個前向激波，在2013.35年探測到一個疑似后向激波[9]。這些激波被認(rèn)為與“旅行者1號”的電場波動儀器探測到的射電爆信號相關(guān)[10]。

學(xué)者嘗試對這些星際空間激波的來源做分析。Liu等[11]首先采用了一維磁流體力學(xué)數(shù)值模擬研究了太陽風(fēng)事件在外日球?qū)訁^(qū)的傳播和演化，提出“旅行者1號”在2013.35年觀測到的射電爆和激波信號可能源自Wind飛船在2012年3月探測到的一系列的行星際日冕物質(zhì)拋射事件。然而，他們的模擬沒有考慮星際等離子體的入流，所以結(jié)果沒有終止激波和日球?qū)忧蕝^(qū)的形成，太陽風(fēng)始終為超聲速；在計算激波傳播時類比地球弓激波上下游的情況并對預(yù)定終止激波以上區(qū)域的激波速度進(jìn)行人為校正。Fermo等[12]采用更復(fù)雜的全球磁流體力學(xué)數(shù)值模擬研究了太陽風(fēng)事件在外日球?qū)拥膫鞑ズ脱莼诰嗳? AU處引入OMNI（https://omniweb.gsfc.nasa.gov）的太陽風(fēng)觀測數(shù)據(jù)，發(fā)現(xiàn)產(chǎn)生的日球?qū)禹敽图げㄎ恢帽葘嶋H觀測到的數(shù)據(jù)要遠(yuǎn)約30 AU。后來，Kim等[13]擴展了該工作，在1 AU的低緯度地區(qū)采用OMNI數(shù)據(jù)，但高緯度區(qū)采用經(jīng)驗給出的高速太陽風(fēng)數(shù)值[14]。結(jié)果顯示多重的太陽風(fēng)共轉(zhuǎn)作用區(qū)對第一個前向激波的形成具有重要的作用，而第二個前向激波可能與太陽風(fēng)共轉(zhuǎn)作用區(qū)和行星際日冕拋射事件的共同作用有關(guān)。

當(dāng)“旅行者1號”在星際空間探測到這3個激波時，“旅行者2號”此時在內(nèi)日球?qū)忧蕝^(qū)的太陽風(fēng)里。Richandson等[15]通過“旅行者2號”的觀測數(shù)據(jù)發(fā)現(xiàn)鞘區(qū)內(nèi)的太陽風(fēng)大尺度壓力脈沖與星際空間的激波具有一定的關(guān)聯(lián)性，經(jīng)過對波速的計算認(rèn)為5個合并太陽風(fēng)作用區(qū)可能最終驅(qū)動了相應(yīng)激波在星際空間的形成。然而，由于“旅行者1號”和“旅行者2號”并不沿著相同的方向飛行，這些壓力脈沖事件與星際激波的因果關(guān)聯(lián)仍有待確認(rèn)。本文將通過一個一維磁流體力學(xué)數(shù)值模擬對這個因果關(guān)系進(jìn)行分析和討論。

1 數(shù)值模型

首先假設(shè)太陽風(fēng)在赤道面附近為沿著徑向方向運動的球?qū)ΨQ流，這里忽略了太陽風(fēng)在傳播過程中可能受到的來自側(cè)面方向太陽風(fēng)的影響。此時理想磁流體力學(xué)方程組在球坐標(biāo)系下可以簡化為沿著徑向 r方向的一維問題為

其中：ρ、u、B和E分別是太陽風(fēng)等離子體的密度、速度、磁場強度和能量密度。

總壓 pT和能量密度 E的具體形式為

其中：p是等離子熱壓。

考慮了太陽風(fēng)等離子體離子與星際中性原子的電荷交換，方程組右邊的 QN、 QMr、 QM?和 QE分別表示因電荷交換引起的質(zhì)量、徑向動量、方位角動量和能量的變化[16-17]。太陽的重力效應(yīng)也被考慮，其中G 是萬有引力常數(shù)，Ms是太陽質(zhì)量。采用MUSCL數(shù)值格式和有限體積方法數(shù)值求解上述方程組[18]，并運用拓展型HLLC黎曼算子完成對方程通量的計算[19]，在時間推進(jìn)上采用龍格-庫塔格式，最終數(shù)值格式在空間和時間上均具有二級精度。

數(shù)值模擬區(qū)域范圍為r方向從1～200 AU，共1萬個網(wǎng)格點。采用非均勻網(wǎng)格，設(shè)N為總格點數(shù)，則格面值ri(i=3,···,N ?1)由代數(shù)式為

其余四個格面值 r1、r2、rN、rN+1由插值給出；而格心值為兩相鄰格面值的中點；δ 參數(shù)表征最大網(wǎng)格與最小網(wǎng)格的差異程度，δ →0時網(wǎng)格趨向均勻網(wǎng)格，這里取 δ = 2.0。這種網(wǎng)格規(guī)模能夠使得模擬獲得較高的分辨率，比如在內(nèi)邊界網(wǎng)格尺寸約0.006 2 AU，在終止激波（～90 AU）附近約0.024 AU，在日球?qū)禹敚ā?20 AU）約為0.03 AU。初始時刻，內(nèi)邊界處太陽風(fēng)取值為：數(shù)密度5 / cm3，速度ur=400 km/s，溫度1.5×105K，磁場強度 Br=2.8 nT。星際介質(zhì)方面，這里只考慮中性原子，其分布形式為[20]

其中：nH0=0.15/cm3為臨近星際空間的中性原子密度初始值； r0=4 AU為星際中性原子進(jìn)入日球?qū)拥纳疃龋镁嚯x日心的距離表征。

簡化計，模擬沒有考慮星際等離子體，因此太陽風(fēng)與星際介質(zhì)相互作用的形式主要以電荷交換形式出現(xiàn)。

從初態(tài)開始，太陽風(fēng)經(jīng)過5年的演化達(dá)到準(zhǔn)穩(wěn)態(tài)，如圖1的點線所示。注意到，由于沒有星際等離子體的作用，終止激波和日球?qū)禹斁鶝]有出現(xiàn)。由于太陽風(fēng)離子與中性原子的電荷交換作用導(dǎo)致太陽風(fēng)被加熱，在外日球?qū)訁^(qū)域其溫度沿著徑向并沒有出現(xiàn)絕熱冷卻的現(xiàn)象，而是在5 AU以上逐漸升高[21]。進(jìn)一步以此為初態(tài)，在1 AU的內(nèi)邊界處代之以2010—2017年間的OMNI、STEREO A和B（以下簡寫STA和STB）等飛船的觀測數(shù)據(jù)。因為STA在2014年9月—2015年11月間部分?jǐn)?shù)據(jù)缺失，而STB自2014年10月起已經(jīng)失去聯(lián)系，所以缺失的數(shù)據(jù)這里以上述恒定的初態(tài)條件代替。后面討論的外日球?qū)蛹げㄊ录?015年之前，假設(shè)激波速度500 km/s，那么從1～120 AU以上區(qū)域激波傳播需要1.14 a左右，因此可以認(rèn)為這種代替的負(fù)面影響是小的。圖1中的實線表示以O(shè)MNI數(shù)據(jù)為輸入，模擬數(shù)據(jù)在2013.42年時刻的物理量徑向分布，可見在太陽風(fēng)事件沿著徑向方向向外傳播，出現(xiàn)的激波信號將會在120 AU以上的星際空間被探測到。

圖1 沿著徑向方向的太陽風(fēng)參數(shù)Fig.1 The profiles of solar wind plasma along radial direction

2 模擬結(jié)果

需要說明的是，在模擬中太陽風(fēng)被簡化為接近赤道面的球?qū)ΨQ流?！奥眯姓?號”在2010—2017年間在HGI（Heliopsheric Inertial Coordinate）坐標(biāo)系下的緯度變化為34.4°N～ 34.8°N，“旅行者2號”從32.3°N～28.8°N；經(jīng)度變化上，“旅行者1號”是173.6°S～175.1°S，“旅行者2號”是216.8°S～218.3°S；可見兩艘飛船均離赤道面有一定的距離。但由于太陽風(fēng)結(jié)構(gòu)在外日球?qū)訁^(qū)域傳播存在大區(qū)域傳播的特性，在經(jīng)度和緯度上覆蓋面廣[22]。這里假設(shè)模擬的結(jié)果能夠一定程度上反映兩艘飛船觀測得到的等離子體大尺度結(jié)構(gòu)的特征。

2.1 與“旅行者2號”的比較

“旅行者2號”最終在2007年距日84 AU處穿越終止激波[2]，之后一直在內(nèi)日球?qū)忧蕝^(qū)直至2018年底被確認(rèn)穿越日球?qū)禹斶M(jìn)入星際空間[23]。給出2011—2017年間的模擬與觀測的動壓對比，在此期間“旅行者2號”從距日94.17 AU處運動至113.14 AU，如圖2所示。綠色橫線段表示發(fā)生在星際空間的等離子體振蕩事件，由“旅行者1號”攜帶的等離子體電場波動儀器測量得到[24]。綠色的豎線對應(yīng)于“旅行者1號”磁強計測量得到的磁場跳變信號，對應(yīng)于星際空間傳播的激波[9]。

一般認(rèn)為圖2中的前3個激波信號與等離子體振蕩時間關(guān)聯(lián)得較好，后者被認(rèn)為來自星際激波前方的低能電子束流，類似于II型的太陽射電爆發(fā)信號[25]。第4個等離子體振蕩信號由于沒有對應(yīng)的磁場跳變，一般認(rèn)為它可能來源于星際其他地方的激波且未被“旅行者1號”的磁強計探測得到。6個字母A～F和對應(yīng)的垂直虛線表示“旅行者2號”觀測到的太陽風(fēng)合并作用區(qū)動壓的局地最高點，這些壓力脈沖被認(rèn)為與“旅行者1號”探測到的等離子體爆發(fā)事件信號相關(guān)，其中C、D和E被認(rèn)為對應(yīng)于第1、3和4個等離子體振蕩[15]。

圖2 2011—2017年采用3種不同內(nèi)邊界數(shù)據(jù)源的太陽風(fēng)動壓對比Fig.2 The comparision of solar wind dynamical pressure between the observation（red）and the simulations（blue）during the years 2011—2017

由于沒有終止激波，模擬得到的動壓相對觀測值較大；出于方便考慮，將觀測的動壓結(jié)果乘以1.5后與模擬結(jié)果做對比。圖2中可見，很難將模擬結(jié)果與“旅行者2號”的觀測做精確的逐一比較，一方面是模擬的有效區(qū)域在近赤道面，與“旅行者2號”的緯度存在較大差異，另外球?qū)ΨQ流不存在與徑向方向相切的流動，實際的太陽風(fēng)傳播要復(fù)雜得多，此外終止激波的缺失也是一方面，以及其他的物理因素如湍流、拾起離子的影響等。

大體上，3個不同的內(nèi)邊界輸入源會導(dǎo)致不同的模擬結(jié)果，但趨勢大致相近，存在近似1年左右的長周期變化[26]。觀測數(shù)據(jù)發(fā)現(xiàn)，激波穿越地球弓激波后速度為原先速度的0.7～1倍左右[27]。因此，設(shè)終止激波在84 AU左右，激波速度在～500 km/s，假設(shè)其在圖2中距離（94～113 AU）上穿越終止激波后以350 km/s運動10～30 AU左右，預(yù)計圖2中的模擬激波信號比實際傳播快0.04～0.12年左右，這種差異在圖中的刻度范圍內(nèi)并不顯著。從圖2中看，OMNI源的模擬結(jié)果對A、B、C和D動壓脈沖具有一定的表現(xiàn)，但對后期兩個壓力脈沖E和F并不明顯；而STA源的模擬結(jié)果對B、D和F具有一定的表現(xiàn)，對A、C和E表現(xiàn)不明顯；對STB源的模擬結(jié)果來說，前期在動壓的變化趨勢上比較一致，但對后期的脈沖比如F表現(xiàn)不明顯，這種情況很可能是因為STB數(shù)據(jù)在自2014年10月后的缺失引起。

模擬結(jié)果表明，對一維模擬來說，1 AU處不同經(jīng)度上的內(nèi)邊界輸入在外日球?qū)訁^(qū)雖然在大尺度結(jié)構(gòu)上具有一定的一致性，但細(xì)節(jié)方面仍具有差異，原因在于OMNI，STA和STB在1 AU軌道上的存在經(jīng)度差異，它們對太陽風(fēng)事件的測量并不總是一致，比如存在有些太陽風(fēng)事件只被某一經(jīng)度上的飛船探測到的情況[28]。因此在與“旅行者2號”的觀測數(shù)據(jù)對比時，需要綜合考慮3個不同內(nèi)邊界輸入的情況。以往人們在模擬外日球?qū)犹栵L(fēng)傳播時側(cè)重于使用OMNI數(shù)據(jù)[13]，該模擬算例表明在更高維度上如何處理好內(nèi)邊界問題是一個很具有挑戰(zhàn)性的難題。

2.2 與“旅行者1號”的比較

“旅行者1號”在2013年8月于距日121.6 AU處穿越了日球?qū)禹擺29]，自那以后人們開始能夠持續(xù)測量星際空間的等離子體波動[10]，銀河宇宙線[30]和磁場信號[9]。然而，由于“旅行者1號”的等離子體儀器自1980年開始就已經(jīng)失效，無法將模擬得到的等離子體參數(shù)與觀測數(shù)據(jù)直接比較。

圖3 2012—2016年用3種不同內(nèi)邊界數(shù)據(jù)源，模擬得到的“旅行者1號”可能探測的動壓數(shù)據(jù)Fig.3 The simulated dynamical pressures that could be encountered by Voyager 1 during the years 2012—2016

圖3分別給出了采用OMNI，STA和STB源作為內(nèi)邊界輸入模擬得到“旅行者1號”自2012—2016年間測量得到的等離子體動壓變化。黑色豎虛線表示“旅行者1號”在2 012.65時刻觀測到的穿越日球?qū)禹數(shù)奈恢茫硗?個豎線表征時刻2 012.904、2 013.356和2 014.644，分別對應(yīng)于“旅行者1號”通過磁場強度跳變觀測得到的前向激波FS1，后向激波RS和前向激波FS2[9]?？梢钥闯觯M得到3個激波附近均具有持續(xù)時間較長的壓力脈沖結(jié)構(gòu)甚至激波組合，時間范圍大致包括2 012.4—2 013.0、2 013.2—2 013.8和2 014.6—2 015.2。

由于模擬并未包含終止激波和日球?qū)禹?，實際上這些壓力脈沖在與終止激波和日球?qū)禹斪饔煤竽軌蜻M(jìn)一步形成激波，進(jìn)而在星際空間被探測到。以O(shè)MNI數(shù)據(jù)的第2個動壓脈沖結(jié)構(gòu)中的2 013.36處的激波為例，注意到在模擬后的第1 086 d到達(dá)85 AU，而在第1 212 d到達(dá)120 AU左右（本文未列出），可以估算激波速度為482 km/s。如果激波在85 AU左右穿越終止激波后速度下降為70%，則在120 AU處激波會遲54 d左右到達(dá)。而在星際空間里，磁場觀測顯示快磁聲波速大概為40 km/s[31]，因此激波穿越日球?qū)禹敽笏俣葧抵猎瓉淼?/10左右，這樣越往后模擬到的動壓脈沖實際對應(yīng)的激波信號延遲越大。Richardson等[15]提出經(jīng)過修正后，“旅行者2號”觀測到的動壓脈沖C、D、E分別對應(yīng)第1、3和4個等離子體振蕩信號，而第2個振蕩信號來源仍然未知。比較圖2和圖3可見，圖2中的C，D和E基本對應(yīng)圖3中的三個動壓脈沖組合，即模擬得到的3個動壓脈沖也同樣可能對應(yīng)于實際觀測到的第1、3和4個等離子體振蕩信號。此外，從模擬結(jié)果看，仍然無法對第2個振蕩信號的來源做出判斷。

以往的三維數(shù)值模擬結(jié)果認(rèn)為第1前向激波FS1來自于太陽風(fēng)的共轉(zhuǎn)作用區(qū)在外日球?qū)訁^(qū)的合并發(fā)展，第2個前向激波是由共轉(zhuǎn)作用區(qū)與行星際日冕物質(zhì)拋射共同作用的結(jié)果[13]。可能的后向激波RS最初被認(rèn)為與2012年3月爆發(fā)的大規(guī)模日冕物質(zhì)拋射事件有關(guān)[11]。

由于模型的限制，這里沒有對模擬得到的3個動壓脈沖結(jié)構(gòu)的驅(qū)動源做分析，后續(xù)將采用更為全面的日球?qū)尤驍?shù)值模擬開展相關(guān)研究。

3 結(jié) 論

本文通過簡化的一維太陽風(fēng)?星際中性原子相互作用的磁流體力學(xué)數(shù)值模型，對太陽風(fēng)在外日球?qū)訁^(qū)的演化進(jìn)行了數(shù)值模擬，其中內(nèi)邊界設(shè)在距日1 AU處，外邊界在200 AU的星際空間。采用1 AU軌道上不同經(jīng)度的OMNI，STA和STB這3個不同觀測數(shù)據(jù)作為輸入條件，時間范圍從2010—2017年。模擬結(jié)果表明，采用不同內(nèi)邊界輸入后，雖然整體上看太陽風(fēng)在外日球?qū)訁^(qū)的大尺度結(jié)構(gòu)比較類似，但細(xì)節(jié)表現(xiàn)會有不同，因此在做數(shù)值模擬時需要做綜合考慮。模擬了“旅行者1號”在星際空間可能探測到的大尺度動壓脈沖結(jié)構(gòu)，這些結(jié)構(gòu)被認(rèn)為能夠與日球?qū)禹斚嗷プ饔貌⑦M(jìn)一步產(chǎn)生激波，從而觸發(fā)“旅行者1號”觀測到的等離子電場振蕩信號，因此可以認(rèn)為這類結(jié)構(gòu)也與觀測到的星際激波信號相關(guān)。