A POROUS, LAYERED HELIOPAUSE

The Voyager 1 and 2 spacecrafts have been mapping the structure of the outer heliosphere as they leave the solar system. In 2005, Voyager 1 (V1) crossed the termination shock (Stone et al. 2005; Burlaga et al. 2005; Decker et al. 2005), where the supersonic solar wind becomes subsonic, and has since been traversing the heliosheath (HS). The heliopause (HP), whose location and structure are unknown, separates the magnetic field and plasma associated with the Sun from that of the local interstellar medium (LISM; Parker 1963; Baranov et al. 1979). The magnetic field in the HS has remained dominantly in the azimuthal (east–west) direction given by the Parker spiral but could rotate and acquire measurable north–south and radial components upon crossing the HP. In ideal (non-dissipative) models of the heliosphere, the local magnetic field is transverse to the boundary and the HP is a tangential discontinuity (Parker 1963; Baranov et al. 1979). However, whether the HP is a smooth interface or breaks up due to instabilities has been the subject of substantial discussion (Fahr et al. 1986; Baranov et al. 1992; Liewer et al. 1996; Zank et al. 1996; Swisdak et al. 2010). The structure of the HP, and in particular whether or not the boundary is porous to some classes of particles, is of great importance because of its impact on the transport of energetic particles into and out of the heliosphere.

Beginning on day 210 of 2012, V1 measured a series of dropouts in the intensities of energetic particles produced in the heliosphere: the anomalous cosmic rays (ACRs) and the lower-energy termination shock particles (TSPs; Webber & McDonald 2013; Stone et al. 2013; Krimigis et al. 2013). Simultaneous with the dropouts were abrupt increases in the Galactic cosmic ray (GCR) electrons and protons and an increase in the magnetic field intensity (Burlaga et al. 2013). Finally, on around day 238, the heliosphere-produced particles dropped to noise levels and the GCRs underwent a final increase. Both have since exhibited no significant variations, which suggests that V1 crossed the HP, with the repeated dropouts and increases perhaps due to radial fluctuations caused by changes in the solar wind dynamic pressure. However, during this time, the direction of the magnetic field remained dominantly azimuthal (Burlaga et al. 2013), consistent with the spacecraft remaining in the HS. While magnetohydrodynamic (MHD) models of the heliosphere suggested that the rotation of the field across the HP at the location of V1 would be small (Opher et al. 2009a, 2009b), the lack of any significant change in direction across the final transition on day 238 suggested that V1 remained within the HS's magnetic domain.

Here we present a model of the magnetic structure of the HP at V1's location that produces particle and magnetic signatures consistent with the observations. By pairing a global MHD simulation with a local particle-in-cell (PIC) simulation, we show that magnetic reconnection can produce a complex, nested set of magnetic islands at the HP. Tongues of LISM plasma penetrate into the HS along reconnected field lines. These tongues trace local depletions of HS plasma and enhancements in the magnetic pressure. A key result of the simulations is that sharp anti-correlated jumps in the HS and LISM number density can occur across the separatrices emanating from reconnection sites while the magnetic field undergoes essentially no rotation. Such behavior undercuts the primary argument, suggesting that V1 has not crossed the HP—that no field rotation was seen on day 238 when the final ACR drop occurred (Burlaga et al. 2013). We therefore suggest that V1 actually crossed the HP on day 209, the time of the last reversal in the azimuthal magnetic field B_T, and that the steady values of the normal, B_N ∼ 0.12 nT, and azimuthal B_T ∼ −0.40 nT, fields since then measure the draped interstellar field just outside of the HP.

To establish local conditions at the HP, we first explore the heliosphere's large-scale structure with a global MHD simulation that includes both neutral and ionized components (and both thermal and pick-up ions in the solar wind; Zieger et al. 2013). The LISM field, B_ISM, has a magnitude of 0.44 nT and a direction defined by α_Bv = 159 and β_Bv = 515, where α_Bv and β_Bv are the angles between B_ISM and the velocity of the interstellar wind v_ISM and the B_ISM–v_ISM plane and the solar equator (for further discussion, see Opher et al. 2009b). The Z-axis is along the solar rotation axis and the X-axis is chosen so that v_ISM lies in the X–Z plane. The MHD simulation did not include the sector zone (where the Parker spiral field periodically reverses polarity due to the tilt between the solar magnetic and rotation axes) since this leads to field reversals that cannot be numerically resolved upstream of the HP, and therefore produces incorrect values of B = |B| (Opher et al. 2011; Borovikov et al. 2011). The solar field polarity corresponds to solar cycle 24, with the azimuthal angle λ (between the radial and T directions in heliospheric coordinates) 90° in the north and 270° in the south.

In Figure 1, B from the simulation reveals the solar wind compression at the termination shock, the downstream HS, and the HP. Profiles (solid curves in Figure 2) along the V1 trajectory of the pick-up (n_pui) and thermal (n_th) ion densities and the azimuthal (B_T) and normal (B_N) magnetic fields near the HP are inputs for the PIC simulations. n_pui decreases from ≈7 × 10⁻⁴ cm⁻³ in the HS to zero in the LISM while n_th rises from 0.003 cm⁻³ to ≈0.08 cm⁻³. B_N (Figure 2(C)) is small at V1's latitude in the LISM. B_T flips direction across the HP, but remains the dominant component on both sides of the boundary (Figure 2(D)). The polar angle δ (the angle between B_N and the equatorial field) in the simulations approaches 14° just outside of the HP, consistent with the steady values seen in the V1 data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The MHD simulation does not match Voyager observations in several respects—the sign of the HS azimuthal magnetic field orientation, the strength of the flows in the HS, and the characteristic scale length of the HP—none of which is essential for calculating initial conditions for the PIC simulations. First, because V1 continued to measure sector boundaries in the HS during 2012, and therefore probably remained in the sector zone, the sign of B_T in the HS in the MHD model is irrelevant since a "correct" model should include the reversals associated with the sector region. Second, in contrast to the simulation, indirect measurements by the V1 LECP instrument indicate little to no normal flow in the HS (Decker et al. 2012). No published global model has explained the observed flows, although simulations including the sectored field (e.g., Opher et al. 2012) are closer to the observations than those presented here (see Pogorelov et al. 2012 for an alternative explanation). Finally, since the MHD model does not include the physics necessary to describe the structure of the HP, the scale length of this transition is not physical, but is instead a numerical artifact. On the other hand, what is essential for input into the PIC simulations is the strength of the HS field and the strength and orientation of the field in the LISM.

The initial profiles of the magnetic field, density, and temperature for the two-dimensional PIC simulations (dotted lines in Figure 2; right-hand scale) were constructed with input from the MHD profiles although, in keeping with the V1 observations, there are no initial flows. The PIC code is written in normalized units based on a field strength B₀ and density n₀ (lengths normalized to the ion inertial length d_i = c/ω_pi, with ω_pi the ion plasma frequency, times to the ion cyclotron time $\Omega _{i0}^{-1}$ and velocities to the Alfvén speed c_A0). In the HS, thermal (n_th = 0.25n₀, $T_{{\rm th}}=0.25 m_ic_{A0}^2$) and pick-up ions (n_pui = 0.01n₀, $T_{{\rm pui}}=15.0m_ic_{A0}^2$) were included as independent species while the LISM only included a thermal component (n_th = 2.0n₀, $T_{{\rm th}}=0.2m_ic_{A0}^2$). The simulations were performed in a domain with dimensions (L_T, L_R) = (409.6d_i, 204.8d_i). The ion-to-electron mass ratio was 25 and the velocity of light was 15c_A0. Not shown in Figure 2 are the three current sheets, of initial half-width 0.5d_i, that produce the sectored HS field. This scale reflects satellite measurements at Earth's magnetopause indicating that such boundaries collapse to kinetic scales (Sonnerup et al. 1981). Pressure balance across each reversal is achieved by adjusting the out-of-plane component B_N (Smith 2001). For HS-appropriate values, n = 10⁻² cm⁻³ and B = 0.3 nT, d_i ≈ 2 × 10⁻⁵ AU, $\Omega _{i0}^{-1} \approx 30 \,{\rm s}$, and c_A ≈ 100  km s⁻¹. Resolving kinetic scales forces the simulation domain to be much smaller than the actual system.

The simulations are evolved with no initially imposed perturbations. Because of the lower density, which leads to a locally higher c_A and effectively thinner current sheets (when normalized to the local d_i), magnetic reconnection first starts in the sectored HS. Small magnetic islands grow on individual current layers in the HS and merge to form larger islands until they are comparable in size to the sector spacing (Figure 3(A)). At the HP, reconnection of the HS and LISM fields produces small islands that also grow, merge, and are eventually compressed by HS islands pushing against the HP (Figure 3(A)). Employing periodic boundary conditions limits the available magnetic flux and so we conclude the simulation when reconnection exhausts the supply. At the real HP, the solar wind continually replenishes reconnected magnetic flux.

Zoom In Zoom Out Reset image size

We define the HP as the location, on a cut in the R direction, of the final flip in the direction of B_T, ΔR/d_i ≈ 25 in Figure 3(D). This feature is quite sharp, ≲ 3d_i, in both the simulation and the Voyager data. Other signatures of the HP (e.g., plasma composition) can be less distinct, as is also the case at the terrestrial magnetopause.

We can independently track every particle in the PIC model and therefore can explore the mixing of the LISM and HS plasmas. The final state is a highly structured distribution of the densities of the LISM (n_LISM) and HS (n_HS), with each experiencing sharp jumps across the separatrices bounding the outflows ejected from reconnection sites (Figures 3(A)–(C)). Particles initially in the LISM continue to dominate the density on the unreconnected LISM field lines, have mixed with HS particles in the nested islands formed from HS–LISM reconnection, and are largely excluded from islands formed from reconnection of the HS sectored field (Figure 3(B)). Particles initially in the HS dominate islands resulting from reconnection of the sectored field, are mixed with LISM particles in the HP islands, and are nearly excluded from unreconnected regions of the LISM (Figure 3(C)).

Radial cuts through the simulation reveal that the increases and decreases in n_LISM and n_HS are typically anti-correlated (Figure 3(G)). Moving from a pure HS magnetic island into an island or outflow jet where LISM and HS plasmas have mixed reduces n_HS and increases n_LISM. Along the cut, the first drop in n_HS occurs downstream of a magnetic separatrix, where HS particles have an open path to the LISM along open field lines (ΔR/d_i ∼ 8 in Figure 3(G)). Similar behavior has been documented in satellite measurements at Earth's magnetopause (Sonnerup et al. 1981) and echoes V1's observations of the anti-correlated variations in the fluxes of ACRs/TSPs and galactic electrons/GCRs (Stone et al. 2013; Krimigis et al. 2013; Webber & McDonald 2013). Most importantly, the cuts reveal that when crossing the final separatrix on the LISM side of the HP and entering pristine LISM plasma, the sharp decrease in n_HS and increase in n_LISM (in the interval ΔR/d_i = 38–50 in Figure 3(G)) are accompanied by no directional change in the magnetic field (Figures 3(D) and (E)). The absence of a directional change in B at locations with strong variations in the particle densities is consistent with one of the most significant of the V1 observations (Burlaga et al. 2013). The fact that this occurs in our simulation on the LISM side of the HP therefore suggests that it may be incorrect to conclude that V1 has not crossed the HP.

The simulation cuts also reveal that local decreases in the HS density typically correspond to increases in the local magnetic field (Figures 3(F) and (G)). The total pressure across the HP is balanced. While the dominant pressure in the HS is from the plasma, the dominant pressure in the LISM is magnetic. Thus, when reconnection opens a path for HS plasma to escape into the LISM and mix with the lower-pressure LISM plasma, there is nothing to balance the total pressure and so the region compresses to increase the magnetic field amplitude. This behavior is primarily seen at separatrix crossings remote from where reconnection locally reduces the magnetic field strength (the interval ΔR/d_i = 38–50 in Figures 3(F) and (G)). In the V1 data, the magnetic field strength is also observed to increase where the local flux of HS plasma decreases (Stone et al. 2013; Krimigis et al. 2013; Webber & McDonald 2013; Burlaga et al. 2013).

Thus, based on our simulations, we suggest that V1's observations of simultaneous drops (increases) in HS (LISM) particle fluxes occur at a series of separatrix crossings outside that are associated with a pair of magnetic islands (Figure 4). At such crossings, the magnetic field direction does not change significantly, while, as seen in the simulation data, particle fluxes can change sharply. Three active reconnection sites at the HP, and associated separatrices with two nested islands, are sufficient to explain the Voyager observations. On day 166, the spacecraft crossed a current layer, on day 190 the flux of HS electrons began dropping, on day 209 another current layer was crossed, and on days 210, 222, and 238 three successive drops (increases) in the HS (LISM) particle fluxes occurred. The day 190 drop in the HS electrons suggests that after this time the magnetic field was no longer laminar and that electrons, with their small Larmor radii and large velocities, could leak into the LISM. Note that the schematic displays the simplest HP structure consistent with the V1 data and does not correspond exactly to the PIC simulation.

Zoom In Zoom Out Reset image size

Islands and x-lines flowing away from an active x-line (e.g., the rightmost x-line in Figure 4) correspond to reconnection sites that developed earlier in time. The separatrices connect to x-lines and, for two reasons, act as bottlenecks to particle transport across the HP. First, at the x-line, the magnetic field turns into the N direction since the R and T components vanish. Thus, the x-line halts the field-aligned streaming of particles across the HP. Second, to the extent that B_N is weak compared with B_T, energetic particles scatter near x-lines further limiting transport. In contrast, downstream of separatrices (to the left in Figure 4) particles can freely stream across the HP. Thus, the day 210 drop in ACRs (rise in GCRs) occurred at the separatrix corresponding to the leftmost x-line of Figure 4 which blocked the transport of ACRs (GCRs). The ACR (GCR) intensity rose (dropped) as the spacecraft crossed field lines that formed a corridor across the HP to the left of the middle x-line in Figure 4. In this region, the flow of GCRs into the HS acts as a sink. The second ACR drop, on day 222, and subsequent recovery are similar. The final drop of the HS particle fluxes, on day 238, occurred at the separatrix of the rightmost x-line. LECP measurements of ACR anisotropies (Krimigis et al. 2013) show that particles propagating parallel to the field dropped more rapidly than those with perpendicular pitch angles. Such behavior is consistent with our schematic—parallel-propagating particles can more easily escape on field lines open to the LISM.

A seeming inconsistency between our PIC simulations and the observations concerns the spatial region where sharp jumps in the ACRs and GCRs occur. In the simulations, the anti-correlated jumps are found on both sides of the HP but, if our illustration is correct, they occur only on the LISM side in the V1 data. The contradiction is possibly due to the simulations being of a two-dimensional system. A magnetic field in a real three-dimensional system will likely be at least mildly stochastic so that wandering field lines will smooth the variability of ACR and GCR intensities far from the HP boundary. More challenging three-dimensional simulations are necessary to explore this issue.

A second issue is the angle δ of the field outside of HP which, for the present MHD simulation, is 14°. B_ISM for this simulation was chosen to match heliospheric asymmetries (Opher et al. 2009a, 2009b). Other MHD simulations based on fitting the IBEX ribbon yield δ ∼ 30° (N. V. Pogorelov 2013, private communication). Nevertheless, our conclusion that significant variations in the ACR and GCR densities can occur in regions with essentially no variation in the magnetic orientation is not sensitive to the value of δ. Any rotation in the field outside of the HP propagates at the local Alfvén speed, which is well below the velocities of the particles of interest. Thus, outside of the HP, separatrices should retain their original LISM orientation in locations where there are significant variations in particle intensity.

The conclusions in this work depend upon the validity of scaling the PIC results to realistic sizes. By repeatedly doubling the size of the computational domain, Schoeffler et al. (2012) demonstrated that magnetic islands in HS current layers grow to the scale of the sector spacing at rates independent of the box size. Combining a typical aspect ratio of reconnection-produced islands (0.1), the typical time between dropouts (10 days), and the speed of V1 with respect to the HS plasma (≈20 km s⁻¹), yields the approximate size of the islands in Figure 4 as 1 AU. This is consistent with the HS sector spacing and supports the extrapolation from the PIC simulation to the schematic.

If nested magnetic islands do occur at the HP, then the dropouts in the HS particle fluxes occurred on the LISM side of the HP on field lines with a LISM source. Thus, according to this picture, V1 crossed the HP on day 208 and has been measuring LISM fields since that time. Our results thus suggest that B_T in the LISM is negative (i.e., has a polarity of 270°).

The authors acknowledge the support of the NSF grant AGS-1202330 to the University of Maryland, and the NSF grant ATM-0747654 and the NASA grant NNX07AH20G to Boston University. The PIC simulations were performed at the National Energy Research Scientific Computing Center and the MHD simulations at the Ames NASA Supercomputer Center. We acknowledge fruitful discussions with L. F. Burlaga, R. B. Decker, M. E. Hill, and E. C. Stone on the Voyager observations. This research benefited greatly from discussions held at the meetings of the Heliopause International Team at the International Space Science Institute.