Propagation of large-scale solar wind events in the outer heliosphere from a numerical MHD simulation

1.   Introduction

The heliosphere is a product of the interaction between expanding solar wind and inflowing interstellar medium. The Sun is centered in the heliosphere and there are two discontinuities, the termination shock and heliopause, that have been explored by the two Voyager probes in the outer heliosphere (e.g. Burlaga et al., 2008; Krimigis et al., 2013). When solar wind transient events propagate into the outer heliosphere, some consequent phenomena, e.g., plasma waves or pulses, may be detected by the spacecraft (e.g. Whang YC and Burlaga, 1985; Burlaga et al., 1997; Wang C and Richardson, 2002). From the year 2012 to 2015, Voyager 1 (V1) observed two forward shocks at 2012.92 and 2014.66 and a possible reverse shock at 2013.35 from magnetic field data in the interstellar medium (Burlaga and Ness, 2016), which showed a sudden enhancement and drop of magnetic field strength, respectively. These shocks were further associated with the radio emissions observed from the plasma wave instrument on V1 (Gurnett et al., 2015). Some attempts have been made to investigate the origins of these shocks in the inner heliosphere. Liu YD et al., (2014) first used a one-dimensional MHD simulation to investigate the origin of a shock and radio emission event in interstellar space, suggesting that the observed radio emissions and the associated shock are the consequence of a series of interplanetary coronal mass ejections (ICMEs) in March 2012. However, due to the lack of inclusion of the inflowing interstellar plasma, there are no termination shock and heliopause in the model. They used an analog of the Earth's magnetosheath to take into account the propagation time in the heliosheath. A more recent complex global MHD model was used to investigate the consequences of solar wind transients in the outer heliosphere (Fermo et al., 2015). In this model, 1-hour averaged OMNI data were used at a spherical boundary at 1 AU, with interstellar shocks reproduced beyond the heliopause ~30 AU larger than the observation by V1 and Voyager 2 (V2). Kim et al. (2017) further extended the model, with the solar wind at 1 AU consisting of OMNI data at low latitude, and the prescribed high-speed solar wind at high latitude based on empirical models (Pogorelov et al., 2013). Their results indicated that multiple co-rotating interaction regions (CIRs) played an important role in the formation of the first forward shock, and determined the second forward shock in conjunction with the ICMEs.

The solar wind origins of interstellar shocks could be traced through observations by another spacecraft in the heliosphere. For example, V2 was still in the inner heliosheath when the three shocks were observed by V1. There was a time lag between the two spacecraft as they detected similar transient events. A recent work based on observations indicated that the pressure pulses measured by V2 may directly correspond with the shocks in the interstellar medium (Richardson et al., 2017). They proposed that five possible merged interaction regions (MIRs) observed by V2 may drive the transients observed by V1 in interstellar space. However, details of the linkage between the pressure pulses in the inner heliosheath and the shock events in the interstellar medium need to be investigated further, which can be carried out by means of numerical simulations.

In this work, we will investigate the propagation of large-scale solar wind events (e.g., MIRs or shocks) in the outer heliosphere by means of a revised one-dimensional MHD model. This model is physically similar to the well-developed 1D MHD model (Wang C et al., 2000; Liu YD et al., 2014) except that an inherent termination shock exists by imposing a constant pressure at the outer boundary, so that the propagation of shocks in the heliosheath is more realistic than their model. The simulation results will be compared with in situ observations from three spacecraft, New Horizons (NH), V2, and V1, as their heliocentric distances increase. The possible relation between the pressure pulses in the inner heliosheath and the shocks in the interstellar medium will be discussed based on the simulation results.

2.   Numerical Model

We use a simplified approach to model the solar wind, assuming a spherically symmetric flow near the equatorial plane and thus no side impact for the solar wind during its radial propagation. Under such an approximation, the normalized MHD equations can be expressed as follows:

where $\rho$, ${{u}}$, ${{B}}$, and $E$ are the plasma density, velocity, magnetic field, and energy density, respectively. The energy density is $E = p/(\gamma-1)+\rho{u^2}/2+B^2/2$, and the total pressure $p^T = p+B^2/2$, where $p$ is the thermal pressure. The charge exchange between plasma and interstellar neutrals is considered in the model, including the mass $Q_N$, momentum $Q_M$ and energy $Q_E$ terms (Pauls et al., 1995). The solar gravity effect appears in the equation, where $G$ is the solar surface gravity, and $M_{\rm s}$ is the mass of the Sun. An MUSCL numerical scheme, integrated with the extended HLLC Riemann solver (Guo XC, 2015), is used to implement the above equations using the finite volume method (Florinski et al., 2013). Overall, the simulation code has a second order accuracy both for the spatial reconstruction and time evolution.

The simulation domain ranges from 1−200 AU, with a total of 10,000 grid points along the radial direction. Using dense grids is beneficial to resolve discontinuities in the simulation; here the grid spacing $\Delta r$ is ~0.0062 AU near the inner boundary, ~0.024 AU at 90 AU, and ~0.03 AU at 120 AU. The typical solar wind conditions are initially set at the inner boundary of 1 AU as follows: number density is 5 ${\rm{cm}}^{-3}$, velocity $u_r$ = 400 ${\rm{km/s}}$, temperature $T$ = $1.5\times 10^5$ ${\rm{K}}$, and magnetic field $B_r$ = 2.8 ${\rm{nT}}$. The interstellar neutrals are simplified as hydrogen, which is fixed along the radial direction as a constant inflow fluid with a velocity of 26.2 ${\rm{km/s}}$ and temperature of 6,300 ${\rm{K}}$. The hydrogen density has a distribution along the radial distance $r$ as the function $n_{\rm H} = n_{\rm H0}\exp(-r_0/r)$ (Axford, 1972), where $n_{\rm H0}$ = 0.15 ${\rm{cm}}^{-3}$ is the number density of undisturbed interstellar hydrogen, and the initial penetration depth $r_0$ = 4 ${\rm{au}}$. The interstellar plasma is not introduced because of the limitation of the symmetric flow approximation, thus no heliopause (HP) is generated in the simulation. A two or three dimensional MHD simulation is necessary to overcome this deficiency. However, similar to the previous approach (Florinksi et al., 2004), the termination shock (TS) can be obtained by applying a boundary condition at 200 AU as $p = 7.1p_{\rm{ISM}}$, where $p_{\rm{ISM}}$ = 0.07 ${\rm{pdyn}}/{\rm{cm}}^2$ is the thermal pressure of the interstellar medium. Compared with the previous non-TS simulation (Liu YD et al., 2014), the introduction of TS in our model is more suitable to capture the solar wind evolution in the inner heliosheath because the previous observation by V2 showed that the solar wind was compressed and heated, and the bulk speed decreased after the TS crossing (Richardson et al., 2008), such that the propagation of solar wind events in the subsonic solar wind are expected to be much different from the supersonic solar wind.

The simulation was run for ~5 years and a steady state of solar wind along the radial direction was obtained, with the results shown as the dotted curves in Figure 1. A TS can be easily identified by physical jumps at the distance of ~84 AU, a value intentionally obtained by choosing a certain $p$ at the outer boundary to agree with the observation by V2 (Richardson et al., 2008). After that, the constant solar wind was replaced with the time-dependent data at 1 AU, which was measured from the spacecraft STEREO A and B (STA and STB), and the virtual spacecraft OMNI, with a duration from 2010.0−2017.0. Unfortunately, the data from STA is only partially available from September 2014 to November 2015, and STB lost contact after October 1, 2014. To overcome this difficulty, the unavailable data will be replaced as constant solar wind input in the initial setting. As we will show later, the main events discussed in this paper will be limited to the time before 2015, so that the negative impact of the unavailable data will be minimized. In Figure 1, the solid curves show the radial profiles of the simulated number density, velocity, magnetic field strength, and temperature of the solar wind from 1−200 AU at 2013.76. The results show large-scale variation of the four physical quantities due to the time-dependent data input at 1 AU from OMNI observations. The TS moves further to ~86 AU, across which the solar wind is compressed and heated, with the speed suddenly decreasing from ~300 to 100 km/s. The magnetic field strength increases after the TS, and the sector structure tends to be closely spaced near the TS due to the slowing of the solar wind. The separation between magnetic sectors becomes larger as radial distance increases, which is similar to the observation by V1 (Burlaga et al., 2014). The positions of NH, V1 and V2 are plotted as the vertical lines in the top panel. Two large amplitude shocks are propagating outward beyond the HP, and are expected to have been observed by the virtual spacecraft at V1 and V2. Details about the shocks will be discussed in the next section.

Figure  1.  The simulated radial profiles of the solar wind number density (N), velocity (V), magnetic field strength ($|B|$) and temperature (T), with the OMNI data input at 1 AU. The dotted curves show the initial state of the simulation at time 2010.0, and the solid curves are the simulation results at time 2014.03, with the positions of NH, V1 and V2 plotted as well. The HP position of 121.6 AU from V1 observation in 2012 (e.g. Krimigis et al., 2013) is shown as the vertical dotted line.

DownLoad: Full-Size Img PowerPoint

More knowledge about the trajectories of the three spacecraft (NH, V1 and V2) needs to be known in order to compare their observations with the simulation. Figure 2 shows the trajectory projections of the spacecraft in the equatorial (xy) and meridional (xz) planes using heliographic inertial (HGI) coordinates. During the period from 2010.0 to 2017.0, STA, STB, and OMNI rotated around the Sun seven times, with only a small change in the longitude for the three distant spacecraft, e.g., 173.6°−174.8° for V1, 216.8°−218.0° for V2, and 199.0°−208.9° for NH, which are shown in panel A as the two adjacent dashed line pairs in three different colors. The trajectories of STA, STB, and OMNI from 2013.0−2013.1 are plotted in the same panel for comparison. When a spacecraft at 1 AU is located at a similar longitude as the other one, it is expected that they will detect the same solar wind events if their separation is relatively small, e.g., within 1 AU. As for NH, V1 and V2, their longitudinal similarity to the three inner ones is not an important factor since the solar wind events will extend to a larger latitude and longitude when propagating to the outer heliosphere (e.g. Intriligator et al., 2005). In panels B and C, the trajectories of V1, V2, and NH are shown as the three arrows. During this period, NH traveled outward at a very low latitude from 4.1° to 5.6°, while V1 was at a latitude from 34.4° to 34.8°, and V2 was from −32.3° to −28.8°. The radial distance for V1 was from 112.1 to 144.3 AU, from the inner heliosheath to the outer heliosheath. For V2, the distance was from 91 to 119.5 AU, located within the inner heliosheath. NH was still measuring the supersonic solar wind from 15.4 to 43.3 AU. As a consequence, the outward solar wind measured by STA, STB, and OMNI is expected to be observed in turn by NH, V2 and V1 in the distant heliosphere.

Figure  2.  Trajectory projections of the six spacecraft, STEREO A and B, OMNI, Voyager 1 and 2, and New Horizons are shown as the solid curves in HGI coordinates, with their directions given by the arrows. From left to right: (a) $xy$ plane near 1 AU, the red, green and black dashed lines indicate the directions of V1, V2 and NH, respectively; the yellow circle shows the orbit of Earth. (b) $xy$ and (c) $xz$ planes from 0 to 150 AU, the dashed and dotted circles represent the locations of the HP and TS, at the radii of 120 and 94 AU, respectively. The yellow circle is the orbit of Jupiter at ~5 AU.

DownLoad: Full-Size Img PowerPoint

3.   Numerical Results

In the numerical model, the solar wind is simplified as a spherically symmetric flow near the equatorial plane. As shown in Figure 2, the trajectories of the three spacecraft in the outer heliosphere indicate that the model is best approximated for NH because of its very low latitude. For the two Voyagers, the latitudes are much higher at ~29°−35°, so the comparison between simulation and observations will be carried out according to the corresponding radial distances; however, the latitude difference cannot be ignored and will be discussed during each comparison.

3.1   Comparison with New Horizons

Figure 3 shows the comparison between the simulation results and the observations by NH from the year 2011.0 to 2015.0, including the density, velocity, and temperature of the solar wind. NH did not measure the magnetic field because of the absence of a magnetometer, thus the simulated magnetic field strength is plotted as a prediction. The trajectory of NH was between 15.45 and 40.3 AU from the Sun. NH simultaneously measured the thermal plasma and pickup ions using the Solar Wind Around Pluto (SWAP) instrument (McComas et al., 2008) during this period, shown as red and green dots in the figure, respectively. The number density slightly decreased along the radial direction both for the simulation and the observations, with a large fluctuation around 0.01 ${\rm{cm}}^{-3}$. However, the fluctuation from SWAP data is much larger than all the three simulations beyond 20 AU, and the observed density is markedly higher than the simulated values at most distances. This discrepancy is similar to the previous global MHD simulation in which OMNI data were applied at the inner boundary (Kim et al., 2016).

Figure  3.  From top to bottom: number density (N), velocity (V), temperature (T) and magnetic field strength ($|B|$) from simulation results (blue) and observations by NH (red: thermal plasma; green: Pickup ions) from 2011–2018. The three origins of the simulation data are listed above from left to right: OMNI, STA, and STB.

DownLoad: Full-Size Img PowerPoint

Solar wind streams with different velocities gradually merge into a larger structure when propagating outward to a larger heliocentric distance (Gazis, 2000). The bulk velocity of solar wind varies steeply in a small time scale of several days, and gradually in a larger time scale from 300−500 ${\rm{km/s}}$. Although STA, STB, and OMNI measured solar wind data for three different longitudes at 1 AU, the corresponding simulated velocities seem to be alike, showing a similar large-scale variation of ~1 year period to the observation. This periodic variation was occasionally found in the solar wind at larger heliocentric distances based on the observations from Pioneer 10 and 11 (Gazis, 1987), and from spacecraft near the Earth as well (Bolton, 1990). Richardson et al. (1994) showed a very strong modulation in solar wind velocities with a period of ~1.3 years from IMP-8 and V2 observations from 1987−1993, speculating that the periodicity was associated with the topology of coronal holes and the formation of open magnetic structures. It seems that these phenomenon have a long-term variation at the solar wind source, in contrast with smaller-scale structures such as CIRs, MIRs and Global MIRs (GMIRs), which are dynamic in origin (Gazis, 1996).

In the numerical model, pickup ions (PUIs) are treated as part of the thermal plasma shown in Equation (1)−(6). More precise description of PUIs as a separate fluid in MHD simulations can be found in the literature (e.g. Wang C and Richardson, 2001; Usmanov et al., 2012). The number density of PUIs is from 0.0002−0.001 ${\rm{cm}}^{-3}$, relatively small compared to that of the thermal plasma shown in the top panels of Figure 3. In the three bottom panels of the figure, the simulated temperature is much larger than that of the observed thermal plasma, and smaller than that of the PUIs. This value reflects a temperature mixture between solar wind and PUIs since the energy gain for ions due to charge exchange is entirely assimilated by the solar wind, thus no individual PUIs are considered in the simulation. In addition, the simulated |B| was shown as a prediction, being mainly around 0.1 nT at these distances. This result is consistent with the previous estimation based on V2 observations (Richardson et al., 2003; Bagenal et al., 2015).

3.2   Comparison with Voyager 2

V2 crossed the TS at a distance of 84 AU in 2007 (Burlaga et al., 2008), and stayed in the heliosheath until the end of 2018 when it crossed the HP and entered the interstellar medium (e.g. Richardson et al., 2019). Here we present the comparison between the simulated and observational dynamic pressure of the solar wind from 2011.0−2017.0, when V2 was located from 94.17 to 113.14 AU, as shown in Figure 4. The six letters A–F and the vertical dotted lines mark the maximum dynamic pressure of MIRs recorded by V2, which are believed to be associated with transient events observed in the interstellar medium by V1 (Richardson et al., 2017). The green thick lines indicate that plasma oscillation events occurred in the interstellar medium, which were detected by the plasma wave instruments (PWS) on V1 (Gurnett and Kurth, 2019). The green vertical solid lines correspond to three shock crossings by identifying jumps in the magnetic field strength from the V1 magnetometer (Burlaga and Ness, 2016). It is clear that the first three shock events are well-correlated with the plasma oscillation events because the latter were caused by low-energy electron beams propagating ahead of interstellar shocks originating from the Sun, such as type Ⅱ solar radio bursts (Bale et al., 1999). The fourth plasma oscillation event in 2015 did not have a shock signature in the magnetic field, and was believed to be produced by electron beams from distant shocks that were not detected by the magnetometer on V1.

Figure  4.  Comparison of dynamic pressure between the observations by V2 (red) and simulations (blue) with three different origins at 1 AU from OMNI, STA, and STB, respectively. The horizontal green thick lines show the four radiation events recorded by the plasma wave instrument of V1. The four short vertical solid lines in the top panel indicate the observed shock crossings by V1. The six vertical dotted lines roughly correspond to the maxima of observed dynamic pressure.

DownLoad: Full-Size Img PowerPoint

Overall the simulated $P_{\rm{dyn}}$ (blue curve) is comparable with the observation by V2 (red curve), except that the simulated $P_{\rm{dyn}}$ has a much lower value than the observed one at a particular time. For example, the simulated $P_{\rm{dyn}}$ approached zero at 2013.1, and kept a low level from 2013.9−2014.3 in the top panel. The observed $P_{\rm{dyn}}$ was larger than ~$1.5\times 10^{-5}$ nPa during the period. This inconsistency is possibly caused by inward propagating waves bouncing from the outer boundary, which may be replaced by the heliopause in reality.

Additionally, this discrepancy may be caused by the difference between the two solar wind flows at latitudes of 34° (for V2) and ~0° (in this model). From April 2010, V1 continuously detected a flow transition region with nearly zero radial speed in the HP upstream, inferred from the anisotropies of convected energetic ion distributions (Krimigis et al., 2011). However, as depicted by the red curve in Figure 4, V2 did not detect a similar signature to V1 until after the latest HP crossing at the end of 2018 (Richardson et al., 2019), which is not shown here. Note that the simulation result reflects the solar wind propagation near the equatorial plane, which locates at a latitude roughly between the two Voyager probes. It is possible for the simulation to detect similar signatures as V1 and V2, especially for large transient events, e.g., MIRs. Because MIRs have large-scale extension both in latitude and longitude when propagating in the outer heliosphere (e.g. Guo X and Florinski, 2014), the simulation results are expected to have coincidence with V2 for some transient events regardless of their different latitudes.

Apparently, the three blue curves illustrate different solar wind evolution in dynamic pressure because of their different inner sources at 1 AU. It seems that the result originated from STA is the most consistent with the V2 observation among the three. For each of the MIRs (A–F) observed by V2, there are corresponding signatures with large enhancement of dynamic pressure (e.g. shocks) from the simulation. For example, V2 detected the largest MIR (panel F) with a maximum dynamic pressure of $8.0\times 10^{-5}$ nPa in the second half of 2015; similarly the most intense simulated pressure pulse with similar magnitude in dynamic pressure occurred at 2015.9 and lasted for several months. The exact timings of the observed and simulated events may differ to some extent, which might be due to the latitude difference. A global simulation may be better to capture the MIRs since the high-speed solar wind from the high latitude region of the solar surface will be taken into account (Kim et al., 2017).

Richardson et al. (2017) noted that the observed pressure pulses C, D and E by V2 were believed to be associated with the shock/plasma oscillation events, shown as the first, third and forth thick green lines in the top panel, respectively; however, the second plasma oscillation event did register a corresponding pressure pulse by V2. Here we find some simulated pressure pulses, for example, at 2012.70 and 2013.2 for the OMNI source, at 2012.80 and 2013.1 for the STA source, or at 2012.78 for the STB source, that may be responsible for the second plasma oscillation event that occurred at ~2013.3. Detailed analysis will be presented in the next section.

3.3   Comparison with Voyager 1

V1 crossed the HP at a heliocentric distance of 121.6 AU in August, 2013 (Stone et al., 2013), after which it continuously measured plasma waves (Gurnett et al., 2015), galactic cosmic rays (Cummings et al., 2016) and magnetic fields (Burlaga and Ness, 2016) in interstellar space. Unfortunately, the plasma instrument of V1 has been non-functional since 1980, thus a direct plasma comparison between simulation and observations is impossible. Figure 5 shows the simulated number density, velocity, and dynamic pressure at V1 locations from 2012.0−2016.0, with three solar wind sources at 1 AU from OMNI, STA, and STB. Similar to Figure 4, the plasma oscillation events are shown as four thick green lines. The interstellar shock events are plotted as vertical lines, namely two solid lines for the foreshocks at 2012.90 (FS1) and 2014.64 (FS2), and one dotted at 2013.36 for the possible reverse shock (RS) (Burlaga and Ness, 2016).

Figure  5.  The simulated plasma number density $N$, radial speed $V_r$, dynamic pressure at V1 locations from 2012 to 2016, with three different solar wind origins at 1 AU from OMNI, STA, and STB, respectively.

DownLoad: Full-Size Img PowerPoint

Due to the absence of interstellar plasma in the simulation, no HP forms beyond the TS in the simulation so the shock propagation in interstellar space is uncertain. As an approximation, the simulated shock in the solar wind beyond the observed HP at 121.6 AU (shown as the vertical dashed line in Figure 5) will be compared with the observation by V1 in interstellar space. At first glance, the first encountered weak shock FS1 was well-reproduced from the two simulation cases of the sources OMNI and STA, in which jumps in density, velocity, and dynamic pressure are clearly identified around 2012.90 at 122.49 AU. Taking the OMNI case for example, the averaged shock speed is estimated to be 314 km/s if we trace the shock back to the history shown in Figure 4, where the shock front locates to 2012.55 and 99.06 AU. This result agrees with the rough estimation of the shock speed in the inner heliosheath from the shock propagation in the Earth's magnetosphere (Richardson et al., 2017). However, the fast-mode speed of the plasma was estimated to be ~40 km/s in interstellar space (Burlaga et al., 2013), thus the shock speed will decrease a small amount after the HP crossing. The time delay $\Delta t$ is estimated by:

where $V_{\rm s0} = 40$ km/s is the assumed shock speed in interstellar space, $V_{\rm s1}$ is the averaged shock speed, and $r_{\rm s}$ and $r_{\rm{HP}}$ indicate the heliocentric distances of the shock and HP, respectively. The simulated shock at 2012.90 will be ~0.09 years later, to be detected after the correction for the shock speed. Therefore, the simulated shock at ~2019.81, which is apparently seen in the cases from OMNI and STA, will arrive at ~2019.87 with a delay of 0.06 years after the shock speed correction beyond HP, being preferably consistent with the observed FS1. In Figure 4, the shock respectively originates at time ~2012.41 and 2012.45 for the sources OMNI and STA, and is between the observed pressure pulses B and C by V2. Based on the numerical result, we propose that both B and C pressure pulses by V2 are potentially linked with the first plasma oscillation event and the corresponding shock (FS1).

A similar correction is made to the simulated shocks at ~2013.15 (distance 123.37 AU) in the cases from OMNI and STB, and at ~2013.20 (distance 123.56 AU) in the case from STA. The shocks are expected to arrive with a delay of 0.18 years from Equation (7), and approach the observed possible RS at 2013.36. This coincidence indicates that the simulated shock may be the driver of the second plasma oscillation event (2013.28−2013.42) and the observed possible RS, although there was no clear corresponding signature observed by V2 (Richardson et al., 2017). Moreover, the shock is traced back to 2012.70 in the case from the OMNI source in Figure 4, where the shock can be grouped as pressure pulse C. This result indicates that a branch of pressure pulse C in the inner heliosheath could be the driver of the second plasma oscillation event beyond HP, and further result in the observed shock by V1.

For the two simulated shocks at 2013.60 and 2013.92 in the OMNI case, the corresponding shock speeds are estimated to be 284 km/s and 300 km/s if measured from the shock events detected at 2013.20 and 2013.54, respectively, shown in the top panel of Figure 4. When similar correction of the shock speed is applied, the two shocks will arrive at the same place about 126 and 171 days later, at around 2013.94 and 2014.39, respectively. There is a large overlap between the duration of the two shocks and the third plasma wave event (2014.1−2014.87), so that we may conclude that the simulated MIRs from 2013.60 to ~2014.64 shown in the case from OMNI may lead to the third plasma oscillation event and the corresponding shock (FS2). For the cases from the STA and STB sources, the profile of the shocks is similar to those from OMNI, except that the first arrival time is earlier for STA at 2013.52, and later for STB at 2013.77. In Figure 4, the origin of these shocks is similar to the observed pressure pulse D, which has been connected with the third plasma wave event by V1 (Richardson et al., 2017). Here we show the numerical evidence that they are connected with FS2 and the associated plasma wave events.

The shock at 2014.88 in the OMNI case has an origin at 2014.45 in Figure 4. The calculated shock speed is about 270 km/s, which is much larger than the fast-mode speed of the plasma beyond the HP. After the correction, it will arrive at 2015.69, which is larger than the start time of the fourth plasma wave event that occurred from 2015.55. Similar shock structures are shown in the case from STA at 2014.83, and STB at 2014.99. They are all expected to be associated with the fourth plasma wave event. Moreover, as Figure 4 shows, these shocks can be traced back to the inner heliosheath at 2014.44, 2014.42, and 2014.55 for the three sources OMNI, STA, and STB, respectively, which are closer to the pressure pulse E than the others.

In summary, our numerical result confirms most of the previous conclusions (Richardson et al., 2017) that the pressure pulses C, D, and E by V2 in the inner heliosheath (shown in Figure 4) are linked with the first, third and fourth plasma wave events by V1 in the interstellar medium, and the possible corresponding shocks shown as the thick green lines in Figure 5. In addition, the simulation result indicates that the pressure pulse B may contribute to the first plasma wave event, and a branch of pressure pulse C is believed to determine the formation of the second plasma wave event and the shock detected by V1.

4.   Summary

In this work, a simplified MHD simulation is conducted to model solar wind propagation near the equatorial plane in the outer heliosphere, with inputs from OMNI, STA, and STB at 1 AU from the years 2010.0−2017.0. The numerical results are respectively compared with the distant observations by NH, V1 and V2, regardless of their latitude differences; this shows that the long-term variations observed by NH are roughly reproduced by the model except for the temperature because the modeled value reflects a mixture of solar wind and PUIs. In the inner heliosheath, the modeled long-term pressure pulses of MIRs are similar to those observed by V2 and believed to be associated with the transient events recorded by V1. There is no HP in the simulation due to the absence of interstellar plasma, thus the simulated shock speed beyond a prescribed HP is reset according to an estimation from V1 observations. Then, the arrival times of the simulated shocks with three different inner sources are compared with V1 observations in the interstellar medium. It is found that each plasma oscillation event, even the shocks identified by the magnetic field strength at V1, agrees with those from the simulations. The possible linkage between the shocks at V1 in interstellar space and the pressure pulses at V2 in the inner heliosheath is established by the simulation.

Acknowledgments

The work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences Grant No. XDB41000000, XDA15052500,XDA17010301, QYZDJ-SSW-JSC028, NNSFC grants 41874171, 41674146, 41731070, and 41774179, in part by the Specialized Research Fund for State Key Laboratories of China and NSSC research fund for key development directions, and the Civil Aerospace "13th Five-Year Plan" PreliminaryResearch in space science (project No: D020301, D030202).We acknowledge the use of NASA GSFC's Space Physics Data Facility's CDAWeb service, and OMNI data.