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

（1.中国科学院 国家空间科学中心 空间天气学国家重点实验室，北京 100190；2.中国科学院大学 地球和行星科学学院，北京 100049）

引 言

源于太阳表面的超声速太阳风等离子体向外径向膨胀时会与包含星际等离子体和中性原子的邻近星际介质相互作用形成一泡状结构即日球层，它们的交界面就是日球层顶[1]。在经历漫长的旅程后，“旅行者1号”（ Voyager 1）和“旅行者2号”（ Voyager 2）飞船分别在距太阳94 AU和84 AU左右穿越终止激波，最终在120 AU左右穿越日球层顶进入星际空间[2-3]。太阳风事件，包括共转太阳风作用区（Corotating Interaction Regions，CIRs）[4]和行星际日冕抛射事件（Interplanetary Coronal Mass Ejections，ICMEs）[5]，在日球层内传输时会引起系列反应，比如飞船发现在几十AU处太阳风事件传输时会逐渐汇聚合并形成更大尺度的等离子体脉冲或激波事件[6-8]

“旅行者1号”在2012年8月穿越日球层顶后，磁强计在2 012.92（十进制年份）和2014.66年（十进制年份）分别探测到两个前向激波，在2013.35年探测到一个疑似后向激波[9]。这些激波被认为与“旅行者1号”的电场波动仪器探测到的射电爆信号相关[10]。

学者尝试对这些星际空间激波的来源做分析。Liu等[11]首先采用了一维磁流体力学数值模拟研究了太阳风事件在外日球层区的传播和演化，提出“旅行者1号”在2013.35年观测到的射电爆和激波信号可能源自Wind飞船在2012年3月探测到的一系列的行星际日冕物质抛射事件。然而，他们的模拟没有考虑星际等离子体的入流，所以结果没有终止激波和日球层鞘区的形成，太阳风始终为超声速；在计算激波传播时类比地球弓激波上下游的情况并对预定终止激波以上区域的激波速度进行人为校正。Fermo等[12]采用更复杂的全球磁流体力学数值模拟研究了太阳风事件在外日球层的传播和演化，在距日1 AU处引入OMNI（https://omniweb.gsfc.nasa.gov）的太阳风观测数据，发现产生的日球层顶和激波位置比实际观测到的数据要远约30 AU。后来，Kim等[13]扩展了该工作，在1 AU的低纬度地区采用OMNI数据，但高纬度区采用经验给出的高速太阳风数值[14]。结果显示多重的太阳风共转作用区对第一个前向激波的形成具有重要的作用，而第二个前向激波可能与太阳风共转作用区和行星际日冕抛射事件的共同作用有关。

当“旅行者1号”在星际空间探测到这3个激波时，“旅行者2号”此时在内日球层鞘区的太阳风里。Richandson等[15]通过“旅行者2号”的观测数据发现鞘区内的太阳风大尺度压力脉冲与星际空间的激波具有一定的关联性，经过对波速的计算认为5个合并太阳风作用区可能最终驱动了相应激波在星际空间的形成。然而，由于“旅行者1号”和“旅行者2号”并不沿着相同的方向飞行，这些压力脉冲事件与星际激波的因果关联仍有待确认。本文将通过一个一维磁流体力学数值模拟对这个因果关系进行分析和讨论。

1 数值模型

首先假设太阳风在赤道面附近为沿着径向方向运动的球对称流，这里忽略了太阳风在传播过程中可能受到的来自侧面方向太阳风的影响。此时理想磁流体力学方程组在球坐标系下可以简化为沿着径向 r方向的一维问题为

其中：ρ、u、B和E分别是太阳风等离子体的密度、速度、磁场强度和能量密度。

总压 pT和能量密度 E的具体形式为

其中：p是等离子热压。

考虑了太阳风等离子体离子与星际中性原子的电荷交换，方程组右边的 QN、 QMr、 QMϕ和 QE分别表示因电荷交换引起的质量、径向动量、方位角动量和能量的变化[16-17]。太阳的重力效应也被考虑，其中G 是万有引力常数，Ms是太阳质量。采用MUSCL数值格式和有限体积方法数值求解上述方程组[18]，并运用拓展型HLLC黎曼算子完成对方程通量的计算[19]，在时间推进上采用龙格-库塔格式，最终数值格式在空间和时间上均具有二级精度。

数值模拟区域范围为r方向从1～200 AU，共1万个网格点。采用非均匀网格，设N为总格点数，则格面值ri(i=3,···,N −1)由代数式为

其余四个格面值 r1、r2、rN、rN+1由插值给出；而格心值为两相邻格面值的中点；δ 参数表征最大网格与最小网格的差异程度，δ →0时网格趋向均匀网格，这里取 δ = 2.0。这种网格规模能够使得模拟获得较高的分辨率，比如在内边界网格尺寸约0.006 2 AU，在终止激波（～90 AU）附近约0.024 AU，在日球层顶（～120 AU）约为0.03 AU。初始时刻，内边界处太阳风取值为：数密度5 / cm3，速度ur=400 km/s，温度1.5×105K，磁场强度 Br=2.8 nT。星际介质方面，这里只考虑中性原子，其分布形式为[20]

其中：nH0=0.15/cm3为临近星际空间的中性原子密度初始值； r0=4 AU为星际中性原子进入日球层的深度，用距离日心的距离表征。

简化计，模拟没有考虑星际等离子体，因此太阳风与星际介质相互作用的形式主要以电荷交换形式出现。

从初态开始，太阳风经过5年的演化达到准稳态，如图1的点线所示。注意到，由于没有星际等离子体的作用，终止激波和日球层顶均没有出现。由于太阳风离子与中性原子的电荷交换作用导致太阳风被加热，在外日球层区域其温度沿着径向并没有出现绝热冷却的现象，而是在5 AU以上逐渐升高[21]。进一步以此为初态，在1 AU的内边界处代之以2010—2017年间的OMNI、STEREO A和B（以下简写STA和STB）等飞船的观测数据。因为STA在2014年9月—2015年11月间部分数据缺失，而STB自2014年10月起已经失去联系，所以缺失的数据这里以上述恒定的初态条件代替。后面讨论的外日球层激波事件在2015年之前，假设激波速度500 km/s，那么从1～120 AU以上区域激波传播需要1.14 a左右，因此可以认为这种代替的负面影响是小的。图1中的实线表示以OMNI数据为输入，模拟数据在2013.42年时刻的物理量径向分布，可见在太阳风事件沿着径向方向向外传播，出现的激波信号将会在120 AU以上的星际空间被探测到。

图1 沿着径向方向的太阳风参数Fig.1 The profiles of solar wind plasma along radial direction

2 模拟结果

需要说明的是，在模拟中太阳风被简化为接近赤道面的球对称流。“旅行者1号”在2010—2017年间在HGI（Heliopsheric Inertial Coordinate）坐标系下的纬度变化为34.4°N～ 34.8°N，“旅行者2号”从32.3°N～28.8°N；经度变化上，“旅行者1号”是173.6°S～175.1°S，“旅行者2号”是216.8°S～218.3°S；可见两艘飞船均离赤道面有一定的距离。但由于太阳风结构在外日球层区域传播存在大区域传播的特性，在经度和纬度上覆盖面广[22]。这里假设模拟的结果能够一定程度上反映两艘飞船观测得到的等离子体大尺度结构的特征。

2.1 与“旅行者2号”的比较

“旅行者2号”最终在2007年距日84 AU处穿越终止激波[2]，之后一直在内日球层鞘区直至2018年底被确认穿越日球层顶进入星际空间[23]。给出2011—2017年间的模拟与观测的动压对比，在此期间“旅行者2号”从距日94.17 AU处运动至113.14 AU，如图2所示。绿色横线段表示发生在星际空间的等离子体振荡事件，由“旅行者1号”携带的等离子体电场波动仪器测量得到[24]。绿色的竖线对应于“旅行者1号”磁强计测量得到的磁场跳变信号，对应于星际空间传播的激波[9]。

一般认为图2中的前3个激波信号与等离子体振荡时间关联得较好，后者被认为来自星际激波前方的低能电子束流，类似于II型的太阳射电爆发信号[25]。第4个等离子体振荡信号由于没有对应的磁场跳变，一般认为它可能来源于星际其他地方的激波且未被“旅行者1号”的磁强计探测得到。6个字母A～F和对应的垂直虚线表示“旅行者2号”观测到的太阳风合并作用区动压的局地最高点，这些压力脉冲被认为与“旅行者1号”探测到的等离子体爆发事件信号相关，其中C、D和E被认为对应于第1、3和4个等离子体振荡[15]。

图2 2011—2017年采用3种不同内边界数据源的太阳风动压对比Fig.2 The comparision of solar wind dynamical pressure between the observation（red）and the simulations（blue）during the years 2011—2017

由于没有终止激波，模拟得到的动压相对观测值较大；出于方便考虑，将观测的动压结果乘以1.5后与模拟结果做对比。图2中可见，很难将模拟结果与“旅行者2号”的观测做精确的逐一比较，一方面是模拟的有效区域在近赤道面，与“旅行者2号”的纬度存在较大差异，另外球对称流不存在与径向方向相切的流动，实际的太阳风传播要复杂得多，此外终止激波的缺失也是一方面，以及其他的物理因素如湍流、拾起离子的影响等。

大体上，3个不同的内边界输入源会导致不同的模拟结果，但趋势大致相近，存在近似1年左右的长周期变化[26]。观测数据发现，激波穿越地球弓激波后速度为原先速度的0.7～1倍左右[27]。因此，设终止激波在84 AU左右，激波速度在～500 km/s，假设其在图2中距离（94～113 AU）上穿越终止激波后以350 km/s运动10～30 AU左右，预计图2中的模拟激波信号比实际传播快0.04～0.12年左右，这种差异在图中的刻度范围内并不显著。从图2中看，OMNI源的模拟结果对A、B、C和D动压脉冲具有一定的表现，但对后期两个压力脉冲E和F并不明显；而STA源的模拟结果对B、D和F具有一定的表现，对A、C和E表现不明显；对STB源的模拟结果来说，前期在动压的变化趋势上比较一致，但对后期的脉冲比如F表现不明显，这种情况很可能是因为STB数据在自2014年10月后的缺失引起。

模拟结果表明，对一维模拟来说，1 AU处不同经度上的内边界输入在外日球层区虽然在大尺度结构上具有一定的一致性，但细节方面仍具有差异，原因在于OMNI，STA和STB在1 AU轨道上的存在经度差异，它们对太阳风事件的测量并不总是一致，比如存在有些太阳风事件只被某一经度上的飞船探测到的情况[28]。因此在与“旅行者2号”的观测数据对比时，需要综合考虑3个不同内边界输入的情况。以往人们在模拟外日球层太阳风传播时侧重于使用OMNI数据[13]，该模拟算例表明在更高维度上如何处理好内边界问题是一个很具有挑战性的难题。

2.2 与“旅行者1号”的比较

“旅行者1号”在2013年8月于距日121.6 AU处穿越了日球层顶[29]，自那以后人们开始能够持续测量星际空间的等离子体波动[10]，银河宇宙线[30]和磁场信号[9]。然而，由于“旅行者1号”的等离子体仪器自1980年开始就已经失效，无法将模拟得到的等离子体参数与观测数据直接比较。

图3 2012—2016年用3种不同内边界数据源，模拟得到的“旅行者1号”可能探测的动压数据Fig.3 The simulated dynamical pressures that could be encountered by Voyager 1 during the years 2012—2016

图3分别给出了采用OMNI，STA和STB源作为内边界输入模拟得到“旅行者1号”自2012—2016年间测量得到的等离子体动压变化。黑色竖虚线表示“旅行者1号”在2 012.65时刻观测到的穿越日球层顶的位置，另外3个竖线表征时刻2 012.904、2 013.356和2 014.644，分别对应于“旅行者1号”通过磁场强度跳变观测得到的前向激波FS1，后向激波RS和前向激波FS2[9]。可以看出，模拟得到3个激波附近均具有持续时间较长的压力脉冲结构甚至激波组合，时间范围大致包括2 012.4—2 013.0、2 013.2—2 013.8和2 014.6—2 015.2。

由于模拟并未包含终止激波和日球层顶，实际上这些压力脉冲在与终止激波和日球层顶作用后能够进一步形成激波，进而在星际空间被探测到。以OMNI数据的第2个动压脉冲结构中的2 013.36处的激波为例，注意到在模拟后的第1 086 d到达85 AU，而在第1 212 d到达120 AU左右（本文未列出），可以估算激波速度为482 km/s。如果激波在85 AU左右穿越终止激波后速度下降为70%，则在120 AU处激波会迟54 d左右到达。而在星际空间里，磁场观测显示快磁声波速大概为40 km/s[31]，因此激波穿越日球层顶后速度会降至原来的1/10左右，这样越往后模拟到的动压脉冲实际对应的激波信号延迟越大。Richardson等[15]提出经过修正后，“旅行者2号”观测到的动压脉冲C、D、E分别对应第1、3和4个等离子体振荡信号，而第2个振荡信号来源仍然未知。比较图2和图3可见，图2中的C，D和E基本对应图3中的三个动压脉冲组合，即模拟得到的3个动压脉冲也同样可能对应于实际观测到的第1、3和4个等离子体振荡信号。此外，从模拟结果看，仍然无法对第2个振荡信号的来源做出判断。

以往的三维数值模拟结果认为第1前向激波FS1来自于太阳风的共转作用区在外日球层区的合并发展，第2个前向激波是由共转作用区与行星际日冕物质抛射共同作用的结果[13]。可能的后向激波RS最初被认为与2012年3月爆发的大规模日冕物质抛射事件有关[11]。

由于模型的限制，这里没有对模拟得到的3个动压脉冲结构的驱动源做分析，后续将采用更为全面的日球层全球数值模拟开展相关研究。

3 结 论

本文通过简化的一维太阳风−星际中性原子相互作用的磁流体力学数值模型，对太阳风在外日球层区的演化进行了数值模拟，其中内边界设在距日1 AU处，外边界在200 AU的星际空间。采用1 AU轨道上不同经度的OMNI，STA和STB这3个不同观测数据作为输入条件，时间范围从2010—2017年。模拟结果表明，采用不同内边界输入后，虽然整体上看太阳风在外日球层区的大尺度结构比较类似，但细节表现会有不同，因此在做数值模拟时需要做综合考虑。模拟了“旅行者1号”在星际空间可能探测到的大尺度动压脉冲结构，这些结构被认为能够与日球层顶相互作用并进一步产生激波，从而触发“旅行者1号”观测到的等离子电场振荡信号，因此可以认为这类结构也与观测到的星际激波信号相关。