ADVECTIVE TRANSPORT OF INTERSTELLAR PLASMA INTO THE HELIOSPHERE ACROSS THE RECONNECTING HELIOPAUSE

The relative motion of the Sun with respect to the local interstellar medium (LISM) leads to the formation of a cavity in the ambient interstellar medium called the heliosphere. The solar wind (SW) and the LISM plasma flow are assumed to be separated by the heliopause (HP), located between a termination shock (TS) in the SW and possibly a bow shock (BS) in the LISM. The inner heliosheath (IHS) is defined to be located between the TS and HP, while the outer heliosheath (OHS; located between the HP and BS, if the BS exists) is a region where significant modification of the LISM flow occurs.

Recent measurements of the Voyager 1 (V1 hereafter) spacecraft provided puzzling observational data that resulted in controversy concerning their interpretation. The V1 spacecraft observed two partial depletions in anomalous cosmic ray (ACR) fluxes, followed by a decrease to the instrumental background ∼122 AU from the Sun (Krimigis et al. 2013). The variations in the ACRs were anticorrelated with changes in the galactic cosmic ray (GCR) fluxes that significantly increased (Webber & McDonald 2013; Krimigis et al. 2013; Stone et al. 2013). At the same time sudden enhancements in the magnetic pressure were observed, but the lack of significant change in the direction of the magnetic field vector has led to doubts whether the observed changes were associated with the HP (Burlaga et al. 2013). However, theoretical studies of the magnetic reconnection at the HP based on numerical simulations have suggested that the V1 observations are consistent with the crossing of a structure related to the HP modified by magnetic reconnection (Swisdak et al. 2013; Strumik et al. 2013). These doubts were finally dispelled with the detection of local plasma oscillations by the V1 plasma wave instrument, where the deduced value of the number density of the surrounding plasma clearly indicated that V1 has entered the interstellar medium (Gurnett et al. 2013). In view of the recent V1 observations the question of the structure of the HP and fine-scale phenomena around the discontinuity surface has become an interesting and timely problem.

We present results of numerical simulations that aim to give a detailed explanation of the V1 measurements of the magnetic field vector changes associated with the ACR boundary crossing. Starting our simulation from a simple configuration of two laminar current sheets (one of them representing the HP), we calculate the time evolution of the plasma parameters and the magnetic field in an area of linear size 4 AU in the normal direction to the HP. Since the linear scale of our computational problem is large in comparison with the ion inertial length λ_i = V_A/ω_ci ≲ 10⁻⁵ AU, we use a magnetohydrodynamical (MHD) approach. Length scales considered in our Letter are large in comparison with the tens or hundreds of λ_i typical for kinetic (particle-in-cell, PIC) simulations of the magnetic reconnection in the heliosphere (Drake et al. 2010; Swisdak et al. 2010, 2013). The PIC simulation results can be linked to larger scales by scaling arguments (Schoeffler et al. 2012; Swisdak et al. 2013), but a fully consistent description of processes on larger (several AU) scales is feasible at present by using the MHD approach. Note also that the scales considered in our Letter are much smaller than the characteristic scales of the HP instabilities caused by charge exchange (Liewer et al. 1996; Zank et al. 1996; Florinski et al. 2005; Borovikov et al. 2008).

It has been suggested that (for the LISM magnetic field of the order of a fraction of nT) reconnection processes at the HP may play an important role in mixing between heliospheric and interstellar plasmas (Macek & Grzedzielski 1985; Fahr et al. 1986). We discuss in this context properties of advective transport of the LISM plasma into the heliosphere that is revealed by our simulations.

MHD equations are solved numerically in 2.5 dimensional geometry (vanishing out-of-plane derivatives) using the high-resolution MUSCL scheme (Kurganov & Tadmor 2000). Resistive and viscous effects are not included explicitly but result from small numerical diffusion of the numerical scheme. A divergence-free magnetic field is ensured by applying a flux-constrained (staggered mesh) approach (Balsara & Spicer 1999). Since phenomena related to the interaction of plasma and neutrals are typically associated with much larger spatial scales, no neutral particle background is included in our model. Periodic boundary conditions are applied in the x-direction and open boundary conditions in the y-direction. Velocity is normalized to the Alfvén speed V_{A, IHS} in the IHS. The MHD equations are scale-less, which allows us to choose the length unit to be 1 AU; the resulting time unit is then T₀ =AU/V_{A, IHS} ≈ 17.5 days. The resolution of the simulation grid is 3840 × 768 points and the computational domain size is L_x × L_y = 20 AU × 4 AU. The density and magnetic field in the simulation are normalized to their averages in the IHS, N₀ and B₀.

The measurements of the magnetic field vector variations by the V1 spacecraft are conventionally presented in the RTN frame (see, e.g., Burlaga & Ness 2012 and references therein). The orientation of the magnetic field vector is given by two angles. The λ angle describes deviation of the magnetic field vector from the radial (R) direction in the radial–tangential (R–T) plane, while the δ angle specifies the deviation of the magnetic field vector from the R–T plane. Figure 1 shows the V1 measurements of the magnetic field vector variations based on data that we digitally extracted from Figure 2 in Burlaga et al. (2013). In the extraction process the data are obtained on a daily timescale, neglecting the details on the fine (48 s) timescale. In Figure 1 two crossings of current sheets observed by V1 are indicated by the λ angle changing from ∼270° to ∼90° (days 162–170 of 2012, V1CS1 hereafter) and then from ∼90° to ∼270° (day 208 of 2012, V1CS2), and are associated with variations of the δ angle and the magnetic field strength B. Note that for the time interval shown in Figure 1, V1 was ∼34° north of the solar equator and the magnetic polarity of the northern hemisphere of the Sun was negative (λ ≈ 90°). Minimum variance analysis (see, e.g., Sonnerup & Scheible 1998 for a description of the method) suggests that the normal vector to the current sheet V1CS2 has the direction (0.7923, 0.0551, −0.6076) in the RTN frame, which is close to the direction of the unperturbed LISM flow. To include this particular physical configuration in our simulations, we assume that the simulation frame xyz is related to the RTN frame through the following transformation:

$$\begin{equation} \left(\begin{array}{@{}c@{}}\hat{\mathbf {e}}_R \\ \hat{\mathbf {e}}_T \\ \hat{\mathbf {e}}_N \end{array} \right) = \left(\begin{array}{@{}rrr@{}}0.0527 & 0.7923 & 0.6079 \\ -0.9983 & 0.0551 & 0.0147 \\ -0.0218 & -0.6076 & 0.7939 \end{array} \right) \left(\begin{array}{@{}c@{}}\hat{\mathbf {e}}_x \\ \hat{\mathbf {e}}_y \\ \hat{\mathbf {e}}_z \end{array} \right). \end{equation} \tag{ 1 }$$

Figure 2 shows the relative orientation of the two frames with respect to the current sheets and magnetic field lines arranged in the initial condition of our simulation. Note that the transformation matrix implies x ≈ −T, the y-axis is the normal vector to the current sheets, and the angle between y- and R-axes is ∼37°. One should also note that the 2.5 dimensional geometry assumption means that the x–y plane is the simulation domain and ∂/∂z = 0 but the magnetic field and velocity vectors may have three non-zero components in general.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Length scales for our simulations are roughly of the same order of magnitude as the separation of the current sheets V1CS1 and V1CS2 observed by V1. The V1 speed (∼17 km s⁻¹) and the time interval (∼43 days) between the current sheet observations yield a separation of ∼0.42 AU with a large uncertainty resulting from neglecting an unknown mean plasma flow. As illustrated in Figure 1 (solid line), in the initial condition for our simulations we set up two current sheets SCS1 (at y = 2 AU) and SCS2 (at y = 2.75 AU), corresponding to the current sheets V1CS1 and V1CS2 observed by the V1 spacecraft. We assume that the normal component of the magnetic field B_y is initially zero everywhere and that B_x/B₀ = 1 and B_z/B₀ = 0 for y < 2 AU. Appropriate rotations of the magnetic field vector in the B_x–B_z plane (by 180° at the SCS1 and 157° at the SCS2) are applied to obtain the variations of the angles λ and δ shown by the solid line in Figure 1. We interpret the region y < 2.75 to be the IHS and the remaining region to be the OHS. Therefore, the current sheet SCS2 can be associated with an initially laminar HP and SCS1 with an occurrence of the heliospheric current sheet (HCS). For SCS2 we arrange the magnetic strength jump B_OHS/B_IHS = 2 and the number density jump N_OHS/N_IHS = 20, similar to our previous work (Strumik et al. 2013). The magnetic strength jump is consistent with the V1 observations (Burlaga et al. 2013). The magnitude of the number density jump is chosen to agree with recent measurements in the IHS by V2 that give N_IHS = 0.0025 cm⁻³ (Richardson & Wang 2012) and with the OHS value N_OHS = 0.05 cm⁻³ determined from the outset of the upward drifting radio emission at frequency ∼1.8 kHz as reported by Gurnett et al. (2013). To ensure total (thermal+magnetic) pressure equilibrium between the IHS and OHS regions we use Δβ = β_IHS − β_OHS = 3.3 at SCS2 where β = 2μ₀p/B² is the ratio of the kinetic and magnetic energy densities. No mean flow in the simulation box is imposed in the initial condition.

To initiate magnetic reconnection we locally impose small-amplitude random perturbations on the plasma pressure in the proximity of the SCS2 location in the simulation box. A locally increased level of noise is supposed to accelerate the growth of the tearing instability and the development of the reconnection sites.

The simulation starting from the initial condition described above results in a time evolution, where the development of the reconnection sites at random locations on SCS2 initiates the growth of magnetic islands at the current sheet (see, e.g., Strumik et al. 2013 for a more detailed discussion of a qualitatively similar scenario). Merging interacting magnetic islands is observed, leading to the growth of the islands' size, while the number of magnetic islands decreases in time. Fluctuations associated with the reconnection at SCS2 propagate throughout the simulation box and initiate reconnection at the current sheet SCS1. At a later time, magnetic islands at the two current sheets start to interact. These dynamical processes cause plasma density and magnetic field compressions and variations in the magnetic field and velocity vectors.

Figure 3 shows comparison of the V1 observations and similar profiles found in our simulations along a virtual spacecraft trajectory r = r₀ + v_SC(t − t₀), where in the simulation frame xyz the initial position is r₀ = (5.35, 1, 0) AU, the virtual spacecraft velocity vector is v_SC/V_{A, IHS} = (0.016, 0.161, 0.058), and t₀/T₀ = 14.5. The virtual spacecraft speed v_SC/V_{A, IHS} ≈ 0.172 in the simulation corresponds to the V1 speed v_V1 ≈ 17 km s⁻¹. The magnetic field vector is assumed constant in the z-direction in the simulation, and linear interpolation in the space domain (based on the simulation grid) and time domain (using frames separated by Δt = 0.1/T₀) is used to obtain the magnetic field vector components between nodes of the simulation grid. Note that the profiles presented in Figure 3 result from the motion of the virtual spacecraft in a spatially inhomogeneous environment accompanied by temporal changes of the spatial distribution of the magnetic field vector. As illustrated in Figure 3 the simulation properly reproduces the following features of the V1 measurements: the magnetic field enhancement of a factor 1.4 between the V1CS1 and V1CS2 crossings (marked by MC1) and strong pulsations of the magnetic field strength after the V1CS2 crossing (marked by MC2). The V1CS1 crossing is associated with δ angle pulsations in the negative range of values, whereas mostly positive values of δ were measured for the V1CS2 crossing. One can see that the simulated δ profiles correspond to the V1 data.

Zoom In Zoom Out Reset image size

Further inspection of the numerical solution allows us to identify the fine-scale physical effects responsible for the appearance of the observational features discussed above. Figure 4(a) shows the spatial distribution of B for t/T₀ = 23 (occurrence of MC1), and Figure 4(b) shows the same but for t/T₀ = 26.1 (occurrence of MC2). The ∼1.4 fold increase of the magnetic field strength (MC1 in Figure 3) between V1CS1 and V1CS2 is apparently related to a magnetic compression occurring when two magnetic islands from neighboring current sheets collide with each other, as illustrated in Figure 4(a). This is a rather typical feature in this physical configuration resulting from the independent motion of magnetic islands at neighboring current sheets during the early stage of development of the islands. When the islands at two current sheets have grown to a size comparable to the distance between the sheets, they start to interact, causing significant magnetic compressions. Even stronger magnetic compressions are seen after the V1CS2 crossing, marked by MC2 in Figure 3. As shown in Figure 4(b), these magnetic strength enhancements are associated with the internal structure of the magnetic islands, appearing during the merging process of two islands significantly different in size. It is again a rather typical situation that reconnection sites initiated at random locations on a current sheet produce magnetic islands that have generally different sizes. Therefore, the merging process is likely to include magnetic islands of different sizes.

Zoom In Zoom Out Reset image size

Due to the purely MHD nature of our simulations, the cosmic ray component is not included in our computations. Therefore, accurate tracking of high-energy particle fluxes is not possible. Heuristically we may argue that since the magnetic compression MC2 in the simulation is part of a magnetic island that has grown at the SCS2 discontinuity, it is likely to contain higher energy particles from the LISM and a deficit of cosmic rays of heliospheric origin. However, a detailed description of cosmic ray fluxes is beyond the scope of this Letter.

Another novel feature revealed by our simulation is advective transport of the relatively dense LISM plasma into the heliosphere through the reconnecting HP. Due to the magnetic strength and density jump across the current sheet SCS2 we obtain a bulge of magnetic separatrices and reconnection jets toward the IHS (small values of y), which is a typical feature of asymmetric reconnection (see, e.g., Swisdak et al. 2010). During the development of the reconnection sites in random locations on the current sheet, reconnection jets from different reconnection sites interact, causing local density and pressure enhancements. The plasma flow velocity component is V_y < 0 on average in these compressions, and thus plasma is transported away from the HP toward the IHS (see Figure 5(a), y < 2.75), which is the lower density region. At a later moment in time, the transport of plasma across the SCS2 occurs in several ducts (see Figure 5(b)) spontaneously appearing on the current sheet as a result of merging of smaller structures. When magnetic islands from the two current sheets in the simulation start to interact, plasma from the density intrusions is transported effectively deep into the IHS. Later, the magnetic sector initially set up between the SCS1 and SCS2 is annihilated, and the reconnection process and transport of the dense plasma across the SCS2 ceases. As illustrated in Figure 5(c), in the final state we observe dense plasma intrusions separated from the SCS2 and layered in the IHS region. Our results indicate that the IHS penetration depth for the intrusions is ∼1.5 AU, but the simulation does not include the IHS turbulence that could presumably transport the density intrusions even deeper into the IHS. The average flux inside the transport ducts estimated from our computations is ∼7.5 × 10⁵ cm⁻² s⁻¹, the typical linear size of the transport ducts is ∼0.5 AU, the separation distance between the ducts is ∼7 AU, and the transport event duration time is ∼1.5 yr.

Zoom In Zoom Out Reset image size

We present results of the modeling of the plasma dynamics in the proximity of the HP that shed light on its internal structure and the surrounding region. The model is validated by finding virtual spacecraft observations corresponding to the magnetic field vector measurements obtained during V1 crossing through the ACR boundary. This allows us to provide an interpretation of particular segments of the observations in terms of the fine-scale physical effects responsible for their appearance. Results of our simulations suggest that characteristic features of the V1 spacecraft measurements can be explained as resulting from dynamical processes driven by magnetic reconnection occurring at two closely separated current sheets, the HCS and the HP. Our simulation also reveals a transport mechanism of the relatively dense LISM plasma through the HP to the IHS region that may be responsible for the formation of density intrusions in the IHS filled by dense plasma of interstellar origin. The values of the basic parameters characterizing the transport process are reported in Section 4 of our Letter.

The advective transport of plasma has been obtained within an MHD approach that describes the dynamics of the thermal plasma component. The plasma experiment for V1 has not been in operation since 1981, and thus plasma density measurements are not available for 2012. The density estimates deduced by Gurnett et al. (2013) from the plasma wave instrument were obtained for much later times than the final jump in the magnetic field strength (seen on day 237 of 2012). The region of density intrusions is predicted by our computations to appear ahead of the jump, thus the measurements by Gurnett et al. (2013) do not allow for observational confirmation of the transport process. Mechanisms of transport of high-speed cosmic rays are generally different from the advective transport of thermal plasma described by the MHD approximation. Therefore, it is not clear at present how the density intrusions revealed by our simulations may be related to abrupt jumps of ACR or GCR fluxes observed by V1 (days 210–237 of 2012). Possible HP crossing by V2 in the future could provide interesting data in this respect.

Our results are interesting in the context of the so-called pile-up problem for the HCS. The HCS foldings are constantly produced by the changing magnetic polarity of the Sun and advected by the SW to the IHS region (Nerney et al. 1995; Czechowski et al. 2010). The recent V1 observations indicate that at least some of the HCS foldings and the magnetic sectors between them may survive the advective transport in turbulent SW and appear very close to the HP. This suggests that the annihilation of magnetic sectors in the direct heliospheric neighborhood of the HP could be a quite likely solution to the pile-up problem for the HCS. For this reason one may expect that the scenario of appearance of density intrusions in the IHS shown in our simulations repeats periodically in nature. Since the increase of the average density in the IHS may change the pressure balance between the IHS and OHS, the effects discussed in our Letter may have consequences for the global structure of the heliosphere. The density profile in the neighborhood of the HP is an interesting question in this context as well. We believe that the transport mechanism revealed by our simulations is important for a proper understanding of the interaction of the SW and interstellar plasmas.

M.S. and W.M. acknowledge support by the Polish National Science Center (N N307 0564 40). R.R. acknowledges support by the Institute of Aviation, HECOLS project, and ISSI.