MAGNETIC FLUX CONSERVATION IN THE HELIOSHEATH INCLUDING SOLAR CYCLE VARIATIONS OF MAGNETIC FIELD INTENSITY

Voyager 1 (V1) and Voyager 2 (V2) measurements in the heliosheath (HS), the region between the termination shock (TS) and the interstellar medium (ISM), are challenging the standard theories and models. In particular, the flows within the HS remain a puzzle. After crossing the TS, V1 observed a steady, almost linear, decrease in the inferred radial velocity of the solar wind plasma, reaching zero and possibly negative values in what was termed the "stagnation region" (Krimigis et al. 2011), while the average V2 speed did not change but remained around 145 km s⁻¹ (Richardson et al. 2009). Another puzzling measurement is the magnetic flux. For a radial flow, the magnetic flux, V_RBR, is expected to be a conserved quantity (Parker 1963). While this is true along V2's trajectory, at V1 the magnetic flux decreased by over an order of magnitude through the HS (Richardson et al. 2013). In order to conserve magnetic flux, the magnetic field strength should have increased to compensate for the decrease in radial velocity along V1, however, the magnetic field stayed between 0.1 and 0.2 nT (Burlaga & Ness 2012).

Models currently can explain neither the direction, the magnitude of the flows, nor the drop of magnetic flux observed by V1. Current steady state heliospheric models predict that the radial velocity component will asymptotically reduce to zero only at the heliopause (HP) and cannot explain the relatively constant speed observed by V2.

There are several models suggested to explain the stagnation region. Opher et al. (2012) suggest that reconnection in the sector region, the region where the solar magnetic field switches polarity due to crossing the heliospheric current sheet, as a possible mechanism. Pogorelov et al. (2013) suggest that the stagnation region can be formed by solar cycle variations of the magnetic axis tilt and extent of the sector region, which causes time variable magnetic barriers at the edge of the HS. These time variable magnetic barriers cause the radial velocity to decrease within this region (Pogorelov et al. 2009). Finally, Lallement et al. (2014) suggest that, as the flow speed decreases to 10 km s⁻¹, the charge-exchange rate increases considerably, leading rapid momentum loss and further deceleration in the HS.

The magnetic flux observations at V1 could be due to a non-ideal MHD effect or temporal variations of the magnetic flux at the solar source. The non-conservation of magnetic flux at V1 could be explained if reconnection were taking place within the sector region (Drake et al. 2010; Opher et al. 2011) since V1 was immersed within this region throughout the HS while V2 went in and out. Lallement et al. (2014) also note that in the stagnation region, where a large fraction of the plasma is undergoing charge-exchange, magnetic flux may not be conserved as well.

While neither of these suggested effects are included in this present work, temporal effects due to the solar wind plasma and magnetic field strength are. The 11 yr solar cycle has been shown to affect the location of the TS and the HP, while greatly affecting the dynamics and structure of the HS (Karmesin et al. 1995; Scherer & Fahr 2003; Zank & Müller 2003; Izmodenov et al. 2008). In their recent time-dependent model, Washimi et al. (2011) used V2 data from 2002 to 2008 as boundary conditions to reproduce V1 and V2 TS crossings. Solar cycle variability of the solar wind density and speed affect the flows within the HS. Pogorelov et al. (2012) suggest that the negative radial velocity observed at V1 could possibly be due to solar cycle variations of the solar wind. Provornikova et al. (2014) shows that solar cycle variations of the solar wind velocity and density can possibly explain the constant radial velocity along V2, but does not explain the difference between the flows along V1 and V2.

The magnetic field also plays an important role in the HS. Burlaga & Ness (2012) observe a decrease in the magnetic field strength directly after the TS and then an increase in 2010 in the region where V1 measured ${{V}_{R}}\sim 0$. Our current steady state model does not match this behavior; it predicts an increase in magnetic field strength after the TS as the flow speed decreases toward the HP. As the solar wind approaches the HP and slows down, the magnetic field piles up, increasing the magnetic field strength if no dissipation processes occur. Time-dependent models suffer a similar problem. The time-dependent model of Provornikova et al. (2014) predicts a magnetic field intensity in the HS over twice as large as measured by V1 and V2.

In addition to solar cycle plasma and magnetic field variations observed from Ulysses, the Pogorelov et al. (2013) model includes variations in the tilt of the Sun's magnetic field derived from data from the Wilcox Solar Observatory (WSO). Washimi et al. (2012) included photospheric magnetic field distributions from WSO in conjunction with interplanetary scintillations (IPS) data to model the expansion of the oscillating heliospheric current sheet into the heliosphere. However, no previous models have included solar cycle variations of the magnetic field strength, which is vital to modeling the plasma environment within the HS.

In this work, we include time-dependent boundary conditions for the solar wind magnetic field into a 3D time-dependent model of the solar wind interaction with the ISM. We discuss the effects of the solar cycle variations of the magnetic field on the plasma parameters in the HS and show that this model can accurately predict the environment observed by V2, while the inclusion of other physical processes are needed to match the flows at V1.

Our model is a global 3D multi-fluid MHD simulation of the outer heliosphere based on the BATS-R-US code that describes the plasma and four neutral hydrogen species (Opher et al. 2003, 2009; Tóth et al. 2012). The model solves the ideal MHD equations for the plasma and a separate set of Eulers equations for the different populations of neutral atoms, which are coupled to the plasma through charge exchange with source terms (Alouani-Bibi et al. 2011). Since we do not describe the pickup ions as a separate component, after the neutrals undergo charge-exchange, we assume that the newly created ions are immediately picked up and take on the characteristics of the local plasma. Due to the various plasma environments in the interaction region between the solar wind and the ISM—the supersonic solar wind inside the TS, the subsonic solar wind in the HS, the subsonic ISM in between the HP and the slow bow shock (Zieger et al. 2013), and the pristine ISM—we define a neutral fluid corresponding to each region. A kinetic description of the hydrogen atoms is needed since the charge exchange mean free path is of the order of the size of the heliosphere (Izmodenov et al. 2000). However, the plasma solutions for the density, velocity, and temperature from multi-fluid MHD and kinetic-MHD models differ only by roughly 5%, due mostly to the different TS and HP locations found from the two models (Alexashov & Izmodenov 2005). This consistency allows us to compare our model to plasma observations by the Voyager spacecrafts.

The model is a cartesian grid that is Sun-centered with dimensions x ± 1500 AU, y ± 1500 AU, and z ± 1500 AU with a spherical inner boundary at 30 AU. We use a non-uniform numerical grid, achieving 0.5 AU resolution within a region in the HS containing both Voyager trajectories. The boundary conditions at 30 AU for the solar wind density and velocity are taken from Provornikova et al. (2014). Provornikova et al. (2014) used solar wind density derived from SOHO/SWAN Lyα emission intensity maps and velocity from IPS data between the years 1991 and 2012 to develop time and latitudinally dependent solar wind density and velocity boundary conditions at 30 AU.

To include solar cycle variations of the magnetic field intensity into our model we fit 27 day averages of the field magnitude average of the magnetic field strength at 1 AU from OMNIWeb⁵ provided by NASA/GSFC and the Space Physics Data Facility (SPDF). We model only the large scale trends in the magnetic field strength from 1991 to 2012, so we use the fit to the overall trend observed in Figure 1, shown by the blue curve. We assume a Parker spiral direction for the solar wind magnetic field and decompose the fit of the field magnitude average of the magnetic field into its components in a spherical coordinate system, radial B_R and azimuthal B $_{\phi }$, using a Parker spiral angle, $\theta ={{{\rm tan} }^{-1}}({{B}_{\phi }}/{{B}_{R}})$, of 45°. Over the course of the solar cycle, the solar wind speed changes, causing the Parker spiral angle to vary on short timescales. Since we are modeling the large-scale structure of the magnetic field intensity variation, the short timescale variation of the Parker spiral angle does not impact our results. We propagate B_R and B $_{\phi }$ derived from the fit of the magnetic field strength at 1 AU to the inner boundary of the model, assuming the Parker solution for the interplanetary magnetic field, ${{{\boldsymbol{B}} }_{{\rm IMF}}}={{B}_{S}}{{\left( \frac{{{R}_{S}}}{r} \right)}^{2}}({\boldsymbol{\hat{r}}} -\frac{{\Omega }(r-{{R}_{s}}){\rm sin} (\theta )}{{{u}_{{\rm sw}}}}\;\hat{\phi })$, where R_S is the radius of the source with a radial magnetic field B_S, Ω is the angular frequency of the Sun, θ is the polar angle, the angle between the Sun's rotational axis and solar latitude, u_sw is the solar wind speed, and ${\boldsymbol{\hat{r}}}$ and $\hat{\phi }$ are unit vectors (Parker 1958). The time-dependent solar wind velocity dependence of ${{B}_{\phi }}$ is included to produce the magnetic field boundary conditions at 30 AU. In order to assess the model boundary conditions, the red curve in Figure 1 compares observations to the boundary conditions of the magnetic field strength at 30 AU propagated back to 1 AU.

Zoom In Zoom Out Reset image size

We assume that the magnetic field intensity has an identical temporal dependence at all latitudes and that the magnitude is modulated by the polar angle in the Parker solution. We justify this assumption a posteriori since the model accurately predicts the magnetic field intensity at V1 and V2, which are at heliolatitudes 35° and −31°, as shown in Figure 2. Figure 3 details the full boundary conditions, latitudinal and temporal, for the solar wind density, speed, temperature, and magnetic field intensity at 30 AU used in this model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, we set the tilt of the solar magnetic axis with respect to the rotation axis to zero and use a monopolar field in the same direction as the ISM field to minimize numerical reconnection effects at the nose of the HP. The magnetic field azimuthal angle, λ, of the solar wind and the ISM is set to 90° which corresponds to a negative magnetic field polarity. V1 and V2 have different magnetic field polarities since they are in different heliospheric hemispheres; however, the results of our model remain the same if the polarity of the ISM and solar wind magnetic field is inverted to have a magnetic field azimuthal angle of 270°. There is no heliospheric current sheet in the monopole model, causing the polarity of the solar wind magnetic field to remain constant within the heliosphere contrary to observations by Burlaga & Ness (2012) of the sector region. This configuration of the solar wind magnetic field also reduces any numerical magnetic dissipation effects in the HS due to the removal of the heliospheric current sheet, allowing us to study only the effects of the time-varying solar wind conditions on the HS. The parameters used for the ISM are n$_{{{H}_{{\rm ISM}}}}=0.18$ cm⁻³, ${{{\rm n}}_{{{p}_{{\rm ISM}}}}}=0.06$ cm⁻³, V$_{{\rm ISM}}=26.4$ km s⁻¹, T$_{{\rm ISM}}=6519$ K, and B $_{{\rm ISM}}=4.37\mu$G. ${{{\boldsymbol{V}} }_{{\rm ISM}}}$ and ${{{\boldsymbol{B}} }_{{\rm ISM}}}$ are offset by 20° while the angle between the plane containing ${{{\boldsymbol{V}} }_{{\rm ISM}}}$ and ${{{\boldsymbol{B}} }_{{\rm ISM}}}$ and the plane of the solar equator was set to 60° since these values provide good agreement with the TS asymmetries observed by V1 and V2 (Opher et al. 2009). We initialize the time dependent run with a longitudinally symmetric steady state solution which includes a latitudinal speed and density variation (see Provornikova et al. 2014). We ran the time-dependent solution over seven iterations, corresponding to 14 solar cycles, until the solutions converged.

Figures 4 and 5 (top panels) compare the HS magnetic field strengths predicted by the models with and without the time-dependent magnetic field boundary conditions and the observed daily averaged magnetic field at V1 and V2, respectively. For the magnetic field of the time-independent magnetic field model, we take a Parker solution with a radial component of $7.17\times {{10}^{-3}}$ nT and an azimuthal component of 0.21 nT in the ecliptic plane at 30 AU. These values were derived using a magnetic field strength of $6.5\times {{10}^{-5}}$ G at 1 AU. Both models predict similar magnetic field intensities in 1991 at both V1 and V2, however, the magnetic field strength in the time-dependent magnetic field model decreases faster than the nominal 1/r dependence of ${{B}_{\phi }}$ in the time-independent model due to the observed decrease in field strength from 1991 to 2012 seen in Figure 1.

Zoom In Zoom Out Reset image size

The realistic time-dependent magnetic field boundary conditions significantly reduced the HS magnetic field intensity and better matches the intensity of the observed magnetic fields. The time-dependent magnetic field model predicts a minor increase in the magnetic field intensity in the HS along both Voyager trajectories. This increase is qualitatively different from the observations, but the magnetic field strength in the HS does match. The gradual decrease in the large scale magnetic field intensity over the course of the solar cycle causes the magnetic field strength to decrease in the outer heliosphere, lowering the magnetic pressure in the supersonic solar wind. While this decrease causes negligible change in the total pressure of the solar wind, since it is dominated by ram pressure, it changes the plasma beta and causes the magnetosonic speed to decrease, raising the magnetosonic Mach number of the TS. With a higher Mach number, the TS in the model is located further from the Sun and closer to the HP and is 6 and 4 AU further than observed at V1 and V2, respectively.

The radial speeds predicted by the model at V1 and V2 are relatively flat and very similar, as found by Provornikova et al. (2014). The model predicts very well the radial flow along V2, matching the flow speed almost exactly and resulting in an improvement from Provornikova et al. (2014). The inclusion of solar cycle variations of the magnetic field strength cause the radial velocity at V2 to be lower in the HS than the time-independent magnetic field model by ∼20 km s⁻¹, matching V2 observations more closely. The decrease in radial velocity to ${{V}_{R}}\sim 0$ along V1 is not predicted by the model. The model predicts a decrease of 75 km s⁻¹ in the first 5 AU after entering the HS, but then predicts a roughly constant speed of 100 km s⁻¹ along V1, very different from observations. The HP location in the model modulates around 150 AU, whereas V1 crossed the HP at 122 AU (Gurnett et al. 2013), making the predicted HS thickness twice as thick as observed. The model does predict a decrease in radial velocity as the Voyager spacecraft move toward the HP, therefore a direct comparison of the flows should be cautioned since the spacecraft is still 40 AU from the model HP.

The time-dependent magnetic field model predicts that magnetic flux is conserved along both Voyager trajectories. For V2, the flow predicted is mostly radial within the HS and V_RBR is constant. Along V1, V_N and V_R are comparable, therefore we calculated the magnetic flux as ${{V}_{\bot }}BL$, where ${{V}_{\bot }}=\sqrt{(V_{R}^{2}+V_{N}^{2})}$ and L is the distance separating the streamlines (Richardson et al. 2013). In both cases we normalized the magnetic field using the magnetic field at 30 AU. The comparison between the magnetic flux predicted by the models and the data along V1 and V2 is shown in the bottom panels of Figures 4 and 5. Similar to the magnetic field intensity and the radial velocity, the time-dependent magnetic field model accurately predicts the magnetic flux at V2, but not the order of magnitude decrease in magnetic flux at V1. We conclude that the drop in magnetic flux along V1 is not due to time-dependence of the solar wind source.

Zoom In Zoom Out Reset image size

While the time-dependent magnetic field model is an improvement to the time-independent magnetic field model, it predicts similar HS conditions at both V1 and V2 contrary to observations. The total plasma speed needed to conserve the magnetic flux in the HS is very different than the flow speeds measured by V1 and remains roughly constant, similar to the behavior at V2. The observations suggest that the two spacecrafts are traveling through very different physical environments. Our model shows that these differences in environment are unlikely caused by solar cycle effects but, rather, are due to missing physical processes in the model that causes the velocity to decrease faster closer to the nose of the heliosphere while the magnetic field remains roughly constant. Some form of non-ideal MHD effect is needed in order to disconnect the radial velocity and the magnetic field strength required for magnetic flux conservation. Therefore, a solution for the flows along V1 warrants a non-ideal MHD approach. Magnetic reconnection is one such form of dissipation. While V2 has gone in and out of the sector region throughout the HS, V1 remained within the sector region where Drake et al. (2010) and Opher et al. (2011) suggest that magnetic reconnection could be ubiquitous.

The inclusion of a solar cycle varying magnetic field into our 3D time-dependent model using realistic solar wind boundary conditions allows us to accurately match the intensity of the magnetic field both at V1 and V2. Our model predicts the Voyager spacecrafts will observe similar plasma parameters within the HS. While this model accurately predicts the observations at V2, it does not reproduce the decrease in radial velocity or drop in magnetic flux observed by V1. We conclude that the change in solar wind magnetic field observed at 1 AU does not cause the reduction of magnetic flux observed by Richardson et al. (2013) and that a non-ideal MHD approach is required to understand this region.

This work was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program—grant NNX14AO14H. M.O. and A.M. acknowledge the support of NASA Grand Challenge NNX14AIB0G and NASA award NNX13AE04G. The calculations were performed at NASA AMES Pleiades Supercomputer. J.D.R. was supported under NASA contract 959203 from the Jet Propulsion Laboratory to the Massachusetts Institute of Technology.