Modeling Shocks Detected by Voyager 1 in the Local Interstellar Medium

The Voyager 1 (V1) spacecraft crossed the heliopause into the local interstellar medium (LISM) at 122 au in 2012 August. This major milestone was preceded by a two-step increase in galactic cosmic-ray flux accompanied by a significant drop in anomalous cosmic-ray intensities leading to the event and was confirmed by the high electron densities observed by the plasma wave science (PWS) instrument in 2013 April–May (Burlaga et al. 2013; Gurnett et al. 2013; Stone et al. 2013). The magnetometer (MAG) on board V1 has observed consistently large magnetic field strength above 0.4 nT in the very local interstellar medium (VLISM) near the heliopause, which was considerably stronger than any previously measured values in the inner heliosheath (Burlaga et al. 2014; Burlaga & Ness 2014, 2016). The interstellar magnetic field (ISMF) exhibited compressible, weakly turbulent fluctuations, while the azimuthal angle λ and the elevation angle δ were observed to increase/decrease linearly. The direction of the ISMF at V1 is significantly different from the Parker spiral direction, indicating a draped ISMF in the VLISM (Burlaga et al. 2014, 2015; Burlaga & Ness 2014, 2016).

Since 2012.65 (DOY 238), V1 MAG recorded two significant jumps in ISMF strength around 2012.92 (DOY 336) and 2014.65 (DOY 236). Both of these disturbances were preceded by electron plasma oscillation events that are strong indicators of shocks (Gurnett et al. 2015) and were followed by relatively quiet periods that ended with abrupt decreases in ISMF strength in 2013.35 (DOY 130) and 2015.37 (DOY 136; Burlaga & Ness 2016). Additionally, V1 PWS observed an intense plasma oscillation event in 2013 April–May, though MAG did not measure a meaningful jump in ISMF strength following the event unlike in 2012 and 2014. It is clear from these observations that V1 is still immersed in a region disturbed by shocks and pressure pulses of solar origin. First suggested by Gurnett et al. (1993), the idea of heliospheric shocks propagating across the heliopause into the LISM has been supported by a number of models (Whang & Burlaga 1994; Zank & Müller 2003; Pogorelov & Zank 2005; Washimi et al. 2011, 2015; Pogorelov et al. 2012).

There have been attempts to model shocks beyond the heliopause using near-Earth solar wind (SW) data. Liu et al. (2014) investigated one-dimensional (1D) propagation of interplanetary coronal mass ejections (ICMEs) and associated shocks from 1 to 120 au using a magnetohydrodynamic (MHD) model and the Wind spacecraft data. The 1D MHD model showed formation of a large merged interaction region (MIR) from a series of ICMEs encountered at 1 au by Wind in 2012 March. The model included the effects of pickup ions, but neglected transition across the termination shock to the inner heliosheath, and also across the heliopause to the LISM. To account for the large uncertainties of shock propagation through the inner heliosheath and VLISM that were not included in the model, Liu et al. (2014) had to make ad hoc adjustments to the shock arrival time at 120 au using shock passage through the Earth's bow shock and the magnetosheath as an analog.

More recently, Fermo et al. (2015) used a fully three-dimensional (3D) multi-fluid MHD–neutral model to simulate shocks in the LISM, which had the advantage of including heliospheric structures lacking in 1D models. The Fermo et al. (2015) model assumed spherical symmetry at 1 au where hourly averaged OMNI data (spacecraft-interspersed, near-Earth SW data available at https://omniweb.gsfc.nasa.gov/) were used as inner boundary conditions. Although the model showed magnetic field and density enhancements related to global MIRs (GMIRs) propagating across the heliopause, the location of the heliopause was incorrect by 20–30 au, which made it difficult to directly assess each individual event observed by V1. The main objectives of our study are as follows: (1) to attribute the shocks observed by V1 in the LISM to SW observations in the OMNI database at 1 au and (2) to identify the relative contributions of ICMEs and corotating interaction regions (CIRs) to the modeled shocks.

We devised a 3D multi-fluid model within the framework of Multi-scale Fluid-kinetic Simulation Suite to simulate the interaction between the SW and the partially ionized LISM (Borovikov et al. 2013; Pogorelov et al. 2014). The model consists of five separate fluids: one plasma fluid and four populations of neutral hydrogen atoms originating in different regions—i.e., in the undisturbed LISM, VLISM around the heliopause, inner heliosheath, and the super-Alvénic SW. While the plasma flow is governed by ideal MHD equations, the flows of neutral hydrogen atoms are described hydrodynamically by means of multi-component Euler equations (Zank et al. 1996; Pogorelov et al. 2006). We make a simplifying assumption that pickup ions resulting from the charge-exchange process between ions and neutral atoms are immediately equilibrated with thermal ions to form an isotropic mixture.

For computational efficiency, we divide the time-dependent simulation into two parts. In the first part, the SW is propagated from 1 to 12 au using a base grid of 256 × 128 × 64 cells in a spherical coordinate system (r, ϕ, θ). Subsequently, we use the time-dependent solutions saved at 12 au as inner boundary conditions for the second part where we employ a base grid of 640 × 128 × 64 cells with the outer boundary defined at 1000 au. In both parts, the radial grid size varies with heliocentric distance such that Δr is approximately 0.03 au at 1 au, 0.13 au at 12 au, 0.4 au at 120 au, and 20 au at 1000 au, for example. In the non-radial directions, the base grid ${\rm{\Delta }}\phi$ and ${\rm{\Delta }}\theta$ are both ∼28. The base grid is too coarse to resolve shocks at large distances, so adaptive mesh refinement (AMR) technique is used to selectively refine the grid in the inner heliosheath and VLISM, particularly around the heliopause. With AMR, the radial grid size is reduced to Δr = 0.036 au at 80 au, 0.043 au at 100 au, 0.025 au at 120 au, and 0.045 au at 140 au, whereas the non-radial grid size becomes as little as 018 at 120 au.

We used OMNI daily averaged plasma and interplanetary magnetic field data to introduce time-dependent effects at 1 au. The inner boundary frame is divided into different regions filled with OMNI data and idealized polar coronal hole (PCH) values whose latitudinal extents vary with time as illustrated in the top panel of Figure 1. In the equatorial OMNI region of each daily frame at 1 au, we filled the 360° longitudinal space centered around Earth using 27 days of data (up to 13 days from the past/future) with a simplifying assumption that the SW propagated radially outward with a spiral magnetic field at 1 au. We estimated the SW parameters in PCHs using empirical correlations (Pogorelov et al. 2013) to best fit the Ulysses data at high heliographic latitudes. The interface between OMNI and PCH regions is 20°–30° wide and are linearly interpolated over. The procedure is described in more detail by Kim et al. (2016), who used the same boundary conditions to reproduce large-scale fluctuations of the SW properties at Ulysses, Voyager, and New Horizons at various distances and latitudes between 1 and 80 au.

Zoom In Zoom Out Reset image size

While we estimated the magnetic field components at 1 au from OMNI $| {\boldsymbol{B}}|$ data in the form of Parker spiral field as done by Kim et al. (2016), we further introduced a heliospheric current sheet (HCS) in the form of a tilted circle (Pogorelov et al. 2007) whose tilt with respect to the Sun's rotation axis changed as a function of time according to the average HCS tilt provided by Wilcox Solar Observatory (WSO; see the bottom panel of Figure 1). Furthermore, the polarity of the magnetic field above and below the HCS at 1 au was instantaneously and simultaneously reversed during solar maximum in 2000 and 2013. In reality, polarity reversal is more complicated and occurs gradually over several months.

At the outer boundary at 1000 au, we set the inflow speed, direction, and temperature of interstellar hydrogen to 25.4 km s⁻¹, 757 ecliptic inflow longitude, −51 ecliptic inflow latitude, and 7500 K, respectively, as suggested by McComas et al. (2015) based on observations by the Interstellar Boundary Explorer (IBEX). We also use interstellar proton and hydrogen atom densities of 0.09 cm⁻³ and 0.154 cm⁻³ as well as ISMF strength and direction of 3 μG, 22699 ecliptic longitude, and 3482 ecliptic latitude that produced the best model fit to the "ribbon" of intense energetic neutral atom emissions observed by IBEX (Zirnstein et al. 2016).

The top two panels of Figure 2 show the model $| {\boldsymbol{B}}|$, the azimuthal angle λ, and the elevation angle δ compared with V1 MAG data in RTN coordinates between 2012.652 and 2016.652, though the model is shown extended out to 2020.0 until the last shock generated by the time-dependent boundary conditions reaches V1. Fluctuations of the model interstellar proton number density are also shown along with a spectrogram of the wideband electric field spectral densities measured by the V1 PWS instrument (Gurnett et al. 2015) in the two bottom panels of Figure 2. The plasma science instrument on V1 is not functioning, but it is still possible to derive density from electron plasma oscillation events detected by PWS assuming the observed electron plasma frequency ${f}_{p}=8980\sqrt{{n}_{e}}$ Hz, where n_e is in cm⁻³ (Gurnett et al. 2015).

Zoom In Zoom Out Reset image size

There is a difference of 0.06–0.11 nT between the model $| {\boldsymbol{B}}|$ and V1 MAG data immediately after the heliopause crossing from 2012.652 to 2012.880. The observed $| {\boldsymbol{B}}|$ decreased from 0.44 to 0.39 nT during this interval, but the model $| {\boldsymbol{B}}|$ remained steady around 0.330 nT, while the model density smoothly increased from 0.0172 to 0.0360 cm⁻³. The azimuthal and elevation angles λ and δ of the model $| {\boldsymbol{B}}|$ changed from 297° to 292° and from 36° to 28°, respectively, whereas the observed λ and δ increased from 286° to 295° and from 16° to 20°. The initial discrepancy of 11° (20°) between the model and observed λ (δ) decreased to 3° (8°) at 2012.880 and remained within 16° (8°) until 2016.650. The observed decrease in $| {\boldsymbol{B}}|$ may be explained by a heliospheric boundary layer resulting from draping of the ISMF around the heliopause (Pogorelov et al. 2017). The heliospheric boundary layer is characterized by increasing density on the LISM side of the heliopause, which is clearly reproduced by this model. The initial discrepancy between the model and observed $| {\boldsymbol{B}}|$, λ, and δ in the vicinity of the heliopause may be associated with uncertainties in the heliopause motion due to time-dependent effects and plasma instability at the heliopause (e.g., Borovikov & Pogorelov 2014).

The first modeled shock arrived at V1 around 2012.890 marked by step-like increases in $| {\boldsymbol{B}}|$ from 0.330 to 0.537 nT and density from 0.0360 to 0.0607 cm⁻³. The ratio of the model $| {\boldsymbol{B}}|$ across the shock ${B}_{2}/{B}_{1}$ is 1.63 which is 11% larger than the observed ratio of 1.47 (Burlaga & Ness 2016). Considering the uncertainties of the observations, these ratios agree favorably. Following the shock passage, the model $| {\boldsymbol{B}}|$ decreased smoothly until the middle of 2013, in general agreement with observations. However, the model did not reproduce the sharp decrease in the observed $| {\boldsymbol{B}}|$ at 2013.353 while showing a moderate bump around 2013.600 that was not observed by V1 MAG. Burlaga & Ness (2016) suggest two possibilities for the abrupt decrease in $| {\boldsymbol{B}}|$: stationary current sheets embedded in the LISM plasma or reverse shock/pressure waves. In the first case, we do not expect our model to reproduce this structure because such structures are not included in the model, but it may be possible in the second case with a more refined grid. We point out that the abrupt decrease at 2013.353 occurred over 2–3 days during which V1 moved away from the Sun by 0.02–0.03 au. However, the radial grid size of the model at V1 at that time was 0.0256 au, which might not be sufficiently small to resolve such a narrow, small-scale structure (i.e., ${B}_{2}/{B}_{1}=1.07$).

The abrupt decrease in the observed $| {\boldsymbol{B}}|$ was followed by a relatively quiet interval from 2013.362 to 2014.641 characterized by small amplitude fluctuations in $| {\boldsymbol{B}}|$ having a Kolmogorov spectrum (Burlaga et al. 2015). In contrast, the model $| {\boldsymbol{B}}|$ increased by ∼0.03 nT around 2013.600 and decreased almost linearly down to 0.358 nT until mid-2014. While we did not attempt to track shocks driven by each individual ICME in the model, we speculate that the moderate increase in the model $| {\boldsymbol{B}}|$ at 2013.600 may have been triggered by a pressure wave delivered by a structure formed primarily by ICMEs in the OMNI data with large longitudinal separation from V1. We further point out the systematically lower model $| {\boldsymbol{B}}|$ and density in comparison to observations during this quiet interval. The pattern of lower-than-observed model values that persisted until at least mid-2016 may be associated with a large rarefaction region trailing the GMIR that drove the 2012 shock in the model. When the rarefaction region reached the heliopause in the model, it caused the heliopause to oscillate, and a significant rarefaction region developed behind the heliopause and propagated into the LISM. Observations do not suggest any dramatic movement of the heliopause or such a large decrease of the ISMF after the 2012 shock. The modeled shocks behind the heliopause are consistently stronger than observed, so we suspect that GMIRs and associated structures in the model may have been somewhat exaggerated in their scale and influence on the heliopause.

The end of the quiet interval was marked by the second modeled shock that arrived at V1 around 2014.665 when the model $| {\boldsymbol{B}}|$ and density jumped from 0.358 to 0.455 nT and from 0.0690 to 0.0891 cm⁻³, respectively. The arrival time of the modeled shock closely matches that of the observed shock at 2014.648, but the ratios of the model $| {\boldsymbol{B}}|$ (1.27) and density (1.29) are 12% and 16% higher than the observed values of 1.13 and 1.11 (Gurnett et al. 2015; Burlaga & Ness 2016), though we estimate the uncertainty of the observed ratios to be ∼5%. Shortly after the passage of the second shock, the model showed another step-like enhancement in $| {\boldsymbol{B}}|$ and density at ∼2014.978. The model $| {\boldsymbol{B}}|$ steadily decreased during the relatively undisturbed period afterwards, in agreement with observations.

Burlaga & Ness (2016) reported a linear decrease in the observed $| {\boldsymbol{B}}|$ from 2014.833 to 2015.370 ending with an abrupt decrease from 0.494 to 0.456. Burlaga & Ness (2016) also pointed out a 28 day quasi-period oscillation of $| {\boldsymbol{B}}|$ during this interval, which was followed by a quiet period of steady decline in the observed $| {\boldsymbol{B}}|$. The model reproduced the steady decrease in $| {\boldsymbol{B}}|$ from late-2014 to mid-2016. We offer comparison between the model and observations until 2016.652 because V1 MAG data are only available up to that date at the moment. However, we extend the model out to 2020.0 until every shock generated by the time-varying boundary conditions propagates out to V1.

The third modeled shock arrives at V1 around 2017.395 when $| {\boldsymbol{B}}|$ and density increase from 0.357 to 0.435 nT and from 0.0913 to 0.113 cm⁻³. The ratios of the model $| {\boldsymbol{B}}|$ and density across the shock are 1.22 and 1.24, respectively. The final shock generated by the time-dependent boundary conditions arrives at V1 around 2019.502 when $| {\boldsymbol{B}}|$ and density increase from 0.369 to 0.474 nT and from 0.107 to 0.139 cm⁻³ with $| {\boldsymbol{B}}|$ and density ratios of 1.28 and 1.30 across the shock. Similar to the mid-2014 shock, the model shows another step-like increase in $| {\boldsymbol{B}}|$ from 0.470 to 0.545 nT and density from 0.138 to 0.161 cm⁻³ around 2019.754 shortly after the shock passage. The widths of these shocks appear broader than the previous shocks due to the relatively large grid size at larger distances. The first two shocks passed V1 at 122.44 and 128.80 au where the radial grid sizes are 0.0252 and 0.0262au, respectively. When the latter two shocks reach V1 at 138.55 and 146.07 au, the radial grid sizes are 0.0427 and 0.0528 au, respectively.

It is interesting to see how much ICMEs affect the shocks modeled at V1 beyond the heliopause. To estimate the contribution of ICMEs, we used the ICME lists for ACE (http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm) and Wind (https://wind.nasa.gov/ICMEindex.php) as reference to identify and remove ICMEs at 1 au, and the resulting data gaps were linearly interpolated. Thus, the boundary conditions at 1 au in this case would only consist of ambient and corotating streams. We followed the same procedure described in the previous section to perform a time-dependent simulation with these boundary conditions. The results are labeled as Model 2 in Figure 3, where Model 1 refers to the original results shown in Figure 2.

Zoom In Zoom Out Reset image size

The arrival time of the late-2012 shock in Model 2 is around 2012.916, and the $| {\boldsymbol{B}}|$ and density ratios across this shock are 1.62 and 1.68, respectively. These values are nearly the same as in Model 1, suggesting that this shock may have been driven by an MIR consisting primarily of CIRs in both models. Contrastingly, the second shock that reached V1 at 2014.665 in Model 1 arrived significantly later in Model 2 at 2015.090. The $| {\boldsymbol{B}}|$ and density ratios across this shock are 1.38 and 1.40, respectively, which are slightly larger than the ratios in Model 1 by 8.7% and 8.5%. Comparison of the models with observations suggests that the mid-2014 shock may have been driven by an MIR consisting of multiple ICMEs and CIRs and that the shock was considerably accelerated by ICMEs. Similarly, the arrival times for the third and the fourth shocks in Model 2 are also delayed by ∼200 days compared to Model 1. We also note that the modest increase in $| {\boldsymbol{B}}|$ of Model 1 around 2013.600 is not present in Model 2, supporting our view that it was affected by inclusion of all ICMEs in the OMNI data, some of which were directed far away from V1 in reality.

Using daily averaged SW parameters from OMNI data as time-varying boundary conditions, we performed a global 3D time-dependent simulation to reproduce shocks propagating beyond the heliopause. The modeled shock arrival times closely match those of the late-2012 and mid-2014 shocks observed by V1, though the changes in the model $| {\boldsymbol{B}}|$ and density across the shocks are slightly larger than observed, considering the uncertainties in the measurements. Furthermore, the model predicts shock arrivals at V1 around 2017.395 and 2019.502. A variant of the model that excludes ICMEs from OMNI data suggests that ICMEs may accelerate some of the shocks significantly.

Although we employed a reasonably high spatial resolution using multiple levels of AMR, the model did not reproduce the relatively steep drops in $| {\boldsymbol{B}}|$ at 2013.353 and 2015.372, or the quasi-periodic oscillations of ISMF after the shock passage in mid-2014. We note that the radial grid of this model is small enough to reproduce daily fluctuations associated with solar activity well into the inner heliosheath. However, the non-radial grid size becomes too large to resolve small-scale fluctuations (e.g., ∼28 days) deeper in the inner heliosheath where the flow develops significant non-radial components as the SW is diverted toward the tail, even with multiple levels of AMR. The computational cost for resolving these fluctuations would have been too excessive. A more detailed investigation of this phenomenon using sufficiently fine Cartesian AMR grids will follow.

The authors acknowledge use of the SPDF COHOWeb database and WSO data (wso.stanford.edu). This work is supported by the NSF PRAC award OCI-1144120 and related computer resources from the Blue Waters sustained-petascale computing project. Supercomputer time allocations were also provided on SGI Pleiades by NASA High-End Computing Program award SMD-15-5860 and on Stampede and Comet by NSF XSEDE project MCA07S033. T.K.K. and N.V.P. acknowledge support from the NSF SHINE project AGS-1358386, SAO subcontract SV4-84017, and NASA contract NNX14AF41G. L.F.B. was supported by NASA contract NNG14PN24P.