Three-dimensional Features of the Outer Heliosphere Due to Coupling between the Interstellar and Heliospheric Magnetic Field. V. The Bow Wave, Heliospheric Boundary Layer, Instabilities, and Magnetic Reconnection
N. V. Pogorelov, J. Heerikhuisen, V. Roytershteyn, L. F. Burlaga, D. A. Gurnett, W. S. Kurth
aa r X i v : . [ a s t r o - ph . S R ] J un Three-dimensional Features of the Outer Heliosphere Due toCoupling between the Interstellar and Heliospheric MagneticField. V.The Bow Wave, Heliospheric Boundary Layer, Instabilities, andMagnetic Reconnection
N. V. Pogorelov, , J. Heerikhuisen , , V. Roytershteyn, L. F. Burlaga, D. A. Gurnett, and W. S. Kurth Department of Space Science, The University of Alabama in Huntsville, AL 35805, USA Center for Space Plasma and Aeronomic Research, The University of Alabama inHuntsville, AL 35805, USA Space Science Institute, Boulder, CO 80301, USA NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA Department of Physics and Astronomy, The University of Iowa, Iowa City, IA 52242, USA [email protected] (corresponding author)
Received ; acceptedAstrophys. J., in press 2 –
ABSTRACT
The heliosphere is formed due to interaction between the solar wind (SW) andlocal interstellar medium (LISM). The shape and position of the heliosphericboundary, the heliopause, in space depend on the parameters of interactingplasma flows. The interplay between the asymmetrizing effect of the interstellarmagnetic field and charge exchange between ions and neutral atoms plays animportant role in the SW–LISM interaction. By performing three-dimensional,MHD plasma / kinetic neutral atom simulations, we determine the width of theouter heliosheath – the LISM plasma region affected by the presence of the helio-sphere – and analyze quantitatively the distributions in front of the heliopause.It is shown that charge exchange modifies the LISM plasma to such extent thatthe contribution of a shock transition to the total variation of plasma parametersbecomes small even if the LISM velocity exceeds the fast magnetosonic speed inthe unperturbed medium. By performing adaptive mesh refinement simulations,we show that a distinct boundary layer of decreased plasma density and enhancedmagnetic field should be observed on the interstellar side of the heliopause. Weshow that this behavior is in agreement with the plasma oscillations of increas-ing frequency observed by the plasma wave instrument onboard
Voyager
1. Wealso demonstrate that
Voyager observations in the inner heliosheath between theheliospheric termination shock and the heliopause are consistent with dissipationof the heliospheric magnetic field. The choice of LISM parameters in this anal-ysis is based on the simulations that fit observations of energetic neutral atomsperformed by
IBEX . Subject headings:
ISM: kinematics and dynamics — magnetic fields — solar wind —Sun: heliosphere
1. Introduction
The interaction of the solar wind (SW) with the local interstellar medium (LISM) isessentially the combination of a blunt-body and a supersonic jet flows. Head-on collisionof the SW and LISM plasma flows creates a tangential discontinuity (the heliopause, HP),which extends far into the wake region (see Fig. 1). The SW flow in the direction parallelto the Sun’s motion resembles a jet immersed into a medium with lower thermal pressure.The LISM plasma is decelerated at the HP, which may result, depending on the LISMparameters, in the formation of a so-called bow shock (BS). The SW flow, on the otherhand, is decelerated due to its interaction with the HP, charge exchange with interstellarneutral atoms, and by the LISM counter-pressure in the heliotail region. Since the neutralhydrogen (H) density in the LISM is greater that the proton density, resonant chargeexchange between ions and neutral ions plays a major role in the SW–LISM interactions(Blum & Fahr 1969; Holzer 1977; Wallis 1971, 1975). In particular, the SW in the tail isdecelerated and cooled down by charge exchange until the heliotail disappears at a fewtens of thousands of AU. Because of the large mean free path, charge exchange and, ingeneral, the transport of neutral atoms should be performed kinetically, by solving theBoltzmann equation. The first self-consistent simulation of this kind was performed byBaranov & Malama (1993) in an axially-symmetric statement of the problem neglectingthe effect of the heliospheric and interstellar magnetic fields (HMF and ISMF). This modelwas extended to time-dependent (Izmodenov et al. 2005b) and 3D flows (Izmodenov et al.2005a) much later. Another class of models assume that neutral atoms can be treated asa fluid, or rather a set of fluids, each of them describing the flow of neutral atoms born inthermodynamically different regions of the SW–LISM interaction. These are usually (i)the unperturbed LISM; (ii) the LISM region substantially modified by the presence of theheliosphere; (iii) the region between the TS and HP; and (iv) the supersonic SW region(Pauls et al. 1995; Zank et al. 1996; Fahr et al. 2000; Florinski et al. 2004; Pogorelov et al. 4 –2006). Such multi-fluid approaches are easily applicable to genuinely time-dependentproblems (see, e.g., Zank & M¨uller 2003; Sternal et al. 2008; Pogorelov et al. 2009a, 2013c),which are very expensive computationally when neutrals atoms are treated kinetically(Izmodenov et al. 2005b; Zirnstein et al. 2015b).Figure 1 shows a typical simulation result that takes into account solar cycle effects.The inner boundary conditions, corresponding to a nominal solar cycle with the radialvelocity, ion density, and temperature in the fast and slow SW, are specified at the Earthorbit ( R = 1 AU). It is assumed that the latitudinal extent of the slow wind varies withan 11-year period from θ = 28 ◦ at solar minima to 90 ◦ at solar maxima. Additionally,the tilt between the Sun’s magnetic and rotation axes varies from θ = 8 ◦ at solar minimato 90 ◦ at solar maxima and flips to the opposite hemisphere at each maximum. Thiscreates a sequence of regions possessing opposite HMF polarities in the heliotail. Weperform all simulations in a so-called heliospheric coordinate system, where the z -axisis aligned with the Sun’s rotation axis, the x -axis belongs to the plane formed by the z -axis and V ∞ , and directed upstream into the LISM. The y -axis completes the rightcoordinate system. The boundary conditions in the SW and LISM are taken from theexisting simulation (Borovikov & Pogorelov 2014) and are for illustration purposes only.In summary, the heliosphere is characterized by the presence of a very long heliotail,which extends to distances exceeding 5,000 au, and is compressed approximately in thedirection perpendicular to the BV -plane. The latter is defined by the LISM velocity andISMF vectors, V ∞ and B ∞ , in the unperturbed LISM. The width of the heliotail in theplane of its maximum flaring (the BV -plane) decreases with distance from the Sun. Thiscreates an illusion that the heliotail disappears when we look at the mutually perpendicularcross-sections shown in Fig. 1. Three-dimensional pictures of the heliosphere can be foundin Borovikov & Pogorelov (2014). We will discuss the boundary conditions in the LISM inSection 3. 5 –In principle, some sort of kinetic treatment of neutral atoms is preferred becausethe charge exchange mean free path is about 50–100 AU, depending on the region ofthe heliosphere and the origin of H atoms. This is especially important for simulationsaimed to provide input to calculations of energetic neutral atoms (ENAs) observed by the Interstellar Boundary Explorer ( IBEX ), see McComas et al. (2017) for a review of themission results over the past 7 years. In particular, the secondary H atoms born in the SWare of importance if we are interested in the distance to which the heliosphere may affect theLISM flow. It is known from theory and simulations (Gruntman 1982; Baranov & Malama1993; Zank et al. 1996) that secondary neutral atoms can travel far upwind where theymay experience charge exchange and affect the LISM flow. A number of the
IBEX ribbonmodels (Heerikhuisen et al. 2010; Chalov et al. 2010; Isenberg 2014; Giacalone & Jokipii2015) involve secondary neutral atoms. Radio emission in a 2–3 kHz range observed by V1 also relies upon global shock waves propagating outward due to various solar events and“the neutral SW” (H atoms born inside the TS) as a source of pickup ions (PUIs) thatinitiate a series of physical processes which ultimately result in the observed wave activity(Gurnett et al. 1993, 2006, 2013, 2015; Cairns & Zank 2002; Pogorelov et al. 2008, 2009b;Mitchell et al. 2009). It is believed that the ring-beam instability of PUIs born in the outerheliosheath (OHS), i.e., in the region of the LISM affected by the presence of the heliosphere,resonantly accelerate ambient electrons by lower hybrid waves. These pre-acceleratedelectrons are further accelerated by transient shocks creating the foreshock electron beams,plasma waves, and radio emission. Secondary neutral atoms are also important to establishthe geometrical size the OHS. The speed and temperature of the unperturbed LISM canbe derived from the properties of He atoms observed by such Earth-bound spacecraft as Ulysses and
IBEX (Witte 2004; Bzowski et al. 2015; McComas et al. 2015). It was shownthat the pristine LISM flow is supersonic and one would expect a bow shock to be formedin front of the HP. However, the LISM is magnetized, so its flow may become subfast 6 –
Model n ∞ , cm − n H ∞ , cm − B ∞ , µ G B ∞ /B ∞ V ∞ , km s − V ∞ /V ∞ T ∞ , K1 0.11 0.165 2 (0 . , − . , . − . , , . . , − . , . − . , , . . , − . , . − . , , . . , − . , . − . , , . . , − . , . − . , , . . , − . , . − . , , . Table 1: Model description for our MHD plasma / kinetic neutral atoms simulation of theSW–LISM interaction.magnetosonic (its speed being less than the fast magnetosonic speed), which will eliminatethe fast-mode bow shock. In principle, slow-mode shocks may still exist in front of the HP(Florinski et al. 2004; Pogorelov et al. 2006, 2011; Zieger et al. 2013) if the angle between V ∞ and B ∞ is small. In this paper, we will show that this is an unlikely scenario in thepresence of charge exchange of LISM ions with secondary H atoms because a fast-modeshock is not just disappearing when B ∞ reaches some threshold value. It is eroding, itsstrength decreasing until no shock is observed.The objective of this paper is to investigate the structure of the LISM region perturbedby the presence of the heliosphere as a function of LISM parameters. In particular, wedetermine the width of the LISM region perturbed by the heliosphere and the contributionof a shocked transition to the overall change of LISM properties across this region. Inaddition, we will consider some issues related to the formation of a boundary layer in theLISM plasma near the HP, development of instabilities and magnetic reconnection, andthe HMF distribution in the presence of the heliospheric current sheet (HCS). Differentproblems require different models for their solution. For this reason, we use an MHD-kineticmodel to investigate the bow shock behavior and the distribution of quantities in theLISM flowing around the heliopause. On the other hand, a multi-fluid approach is moreappropriate for modeling the HP instabilities. Comparisons between MHD-kinetic and 7 – Model ( λ V ∞ , β V ∞ ) ( ◦ ) ( λ B ∞ , β B ∞ ) ( ◦ )1 (255.7, 5.0) (233.20, 29.86)2 (255.7, 5.0) (229.61, 32.83)3 (255.7, 5.0) (228.34, 33.81)4 (255.7, 5.0) (226.99, 34.82)5 (255.7, 5.0) (224.46, 36.61)6 (255.7, 5.0) (222.31, 38.02) Table 2: The directions of the unperturbed LISM velocity and ISMF vectors from Table 1in the ecliptic J2000 coordinate system.multi-fluid models are presented by Alexashov & Izmodenov (2005); Heerikhuisen et al.(2006); M¨uller et al. (2008) and Pogorelov et al. (2009c).
2. Constraints on the LISM Properties
As discussed in the Introduction, the LISM velocity vector and temperature can bederived from the He observations in the inner heliosphere. The rest of quantities in theunperturbed LISM are not measured directly. Following Zirnstein et al. (2016), they can bechosen to satisfy a number of observational results:1. By analyzing the Ly- α backscattered emission in the Solar and HeliosphericObservatory ( SOHO ) solar win anisotropy (SWAN) experiment, Lallement et al.(2005, 2010) discovered a deflection ( ∼ ◦ ) of the neutral H atom flow in theinner heliosphere from its original direction, V ∞ . These two directions define aso-called hydrogen deflection plane (HDP). MHD-plasma/kinetic-neutrals simulations(Izmodenov et al. 2005a; Pogorelov et al. 2008, 2009b; Katushkina et al. 2015) showedthat the average deflection occurs predominantly parallel to the BV -plane.2. It is also possible to restrict LISM properties by fitting the IBEX ribbon of enhanced 8 –ENA flux (see Schwadron et al. 2009, where it was determined that the directionstowards the ribbon strongly correlate with the lines of sight perpendicular to theISMF lines draping around the HP). As shown by Pogorelov et al. (2010) andHeerikhuisen & Pogorelov (2011), the position of the ENA ribbon strongly dependson rotation of the BV -plane about the V ∞ vector. Kinetic ENA flux simulations ofHeerikhuisen et al. (2010, 2014); Heerikhuisen & Pogorelov (2011); Zirnstein et al.(2015a,b, 2016) reproduced the ribbon using the BV -plane consistent with the HDP.It is worth noticing here that the accuracy of SOHO
SWAN measurements allowssubstantial margins in the determination of the BV -plane (see, e.g., Pogorelov et al.2007). It is of interest from this viewpoint that the BV -plane from Zirnstein et al.(2016) lies in the middle of the range derived from the HDP analysis. Astrophysicalobservations restricting the ISMF properties (Frisch et al. 2015) are also consistentwith the above considerations. New IBEX results (McComas et al. 2017) also suggestsecondary ENA’s as accepted ribbon source.3. One-point-per-time, in situ measurements performed in the LISM by V1 , also providerestriction on the direction and strength of B ∞ . E.g., numerical simulations ofPogorelov et al. (2009b) provided B ∞ · R = 0 ( R is a unit vector in the radialdirection) directions consistent with the IBEX ribbon (McComas et al. 2009). Thesame choice of the LISM properties also reproduced the elevation and azimuthalangles in the ISMF beyond the HP (see Pogorelov et al. 2013c; Borovikov & Pogorelov2014). Simulations presented in Zirnstein et al. (2016) are also restricted by theHP position and magnetic field angles observed by V1 . It is certainly possible (see,e.g., Pogorelov et al. 2015) to shift the HP position to 122 AU where the HP wascrossed by V1 (Gurnett et al. 2013). However, the TS position in the V1 directionbecomes substantially smaller (by ∼
20 AU) than at the time of crossing. Thisshould not be discouraging, since there is no information about the TS position 9 –after
Voyagers crossed it. Numerical simulations based on
Ulysses observations(Pogorelov et al. 2013c) indeed show that the TS was moving inward between 2004and 2010. A decrease in the SW ram pressure from one solar cycle to another couldcontribute to the HP shift inward by a few AU. As shown in Malama et al. (2006)and Pogorelov et al. (2016), the IHS width also decreases when PUIs are treated as aseparate plasma component.4. The H density at the TS derived from PUI measurements (Bzowski et al. 2009) canalso be used to constrain the models.5. Modeled anisotropy in the 1–10 TeV galactic cosmic ray (GCR) flux and itscomparison with multiple air shower observations can also improve our knowledge ofthe LISM (Schwadron et al. 2014; Zhang et al. 2014; Zhang & Pogorelov 2016).
3. Bow Wave and Heliospheric Boundary Layer
McComas et al. (2012) and Zank et al. (2013) investigated conditions to be satisfiedfor the bow shock to exist. However, a definitive answer to this question depends on thedetails of the global SW–LISM simulation and is not readily available from a direct analysisof the LISM properties far away from the HP. The LISM near the HP is a weakly collisionalmedium (the mean free path with respect to Coulomb collisions is about 1 AU, Baranov& Ruderman, 2013), so the bow shock, if it exists, is rather well described by ideal MHDequations, which have the t -hyperbolic type. The source terms in these equations are dueto charge exchange and therefore contain no delta-functions. As a result, the Hugoniot-typeboundary conditions at a bow shock cannot be modified by such source terms. Thus, theonly “shock structure” to be expected is of numerical origin. However, charge exchangemodifies plasma and magnetic field both in front and behind the shock in a way unknown 10 –before the problem of the SW–LISM interaction is solved as a whole. It has been shown inthe gas dynamic plasma ( B = )/ kinetic neutrals simulations of Izmodenov (2000), theincrease in the LISM neutral H density results in a weaker bow shock. This conclusionholds in the presence of ISMF.Figure 2 shows the distributions of the plasma density and fast magnetosonic Machnumber M f = v/a f along the x -axis behind the heliopause. Figure 3 shows the plasmadensity distributions for the same set of parameters in the meridional plane. The LISMparameters are taken from our previous simulation in Zirnstein et al. (2016), and aresummarized in Table 1. For convenience, the V ∞ and B ∞ directions are also given inthe J2000 ecliptic coordinates (Table 2). The SW properties are somewhat different fromZirnstein et al. (2016), being closer to OMNI data averaged over a substantial period oftime 2120 days starting between 2010 DOY 1 and 2015 DOY 294) to obtain a sphericallysymmetric distribution at 1 AU: the plasma density is 5.924 cm − , velocity 409.8 km s − ,temperature 82,336 K, and radial HMF component 39 µ G.It is clear from this figure that the LISM properties are substantially modified bycharge exchange, the changes being stronger behind the fast magnetosonic shocks seenin Figs. 2a–2b. This is not surprising since (1) the density of secondary neutral atomsaffecting the LISM ions decreases with heliocentric distance and (2) the LISM plasmadensity increases as it approaches the heliospheric boundary layer (HBL) on the interstellarside of the HP. It can be seen from Fig. 2c that the shocked transition essentially disappearsalready at B ∞ = 2 . µ G: M f ≈ .
02 ahead of it. For B ∞ = 3 µ G, M f becomes smallerthan 1 smoothly. Further increase in B ∞ makes the density variation from the unperturbedLISM to the HP weaker, but wider. This effect is well-known (see, e.g. Izmodenov et al.2005a; Pogorelov et al. 2006; Borovikov et al. 2008a; Zank et al. 2010; Heerikhuisen et al.2015). However, the solutions presented here are for the first time obtained using adaptive 11 –mesh refinement (AMR) in the OHS region (see the patch edges in Fig. 2a), which madeit possible to identify shocks inside the OHS plasma and the HBL near the HP. In someof presented simulations, such shocks are situated inside regions of substantial variationof the LISM properties. From this viewpoint, the OHS itself can be called a “bow wave,”which forms in front of the HP due to the SW–LISM collision. This bow wave can eitherhave a shock inside it or not, depending on the full set of SW and LISM parameters. Thepresence of a shocked transition inside the bow wave clearly is not determined by thecondition M f > R E for the SW Alfv´en number M A = 8 and rapidly decreasing as M A increases. Whensimplistically scaled to the size of the outer heliosphere, a depletion layer at the HP wouldbe 1%–2% of the heliocentric distance of the latter, which gives us about 1.25–2.5 AU.The ions of the terrestrial magnetosheath are typically observed to have bi-Maxwellianvelocity distributions with T ⊥ /T k >
1, where the superscripts ⊥ and k denote directionsperpendicular and parallel to the background magnetic field. Anderson & Fuselier (1993)showed that the temperature anisotropy may be important for the PDL formation because 12 –it can launch an electromagnetic ion cyclotron instability, which makes scattered ionspropagate along the magnetic field and leave the draping region. In the more recent analysisof Fuselier & Cairns (2013), the PDL width is in the range 1–5 AU, which emphasizes thenecessity to incorporate microphysical processes of the magnetic field draping/depletionlayer formation into the PDL analysis (see also Gary et al. 1994; Denton & Lyon 1996). Itis interesting to notice, however, in this connection that a HBL exists in on the LISM sideof the HP in simulations without magnetic field (Baranov & Malama 1993), which meansthat charge exchange affects them somehow. To separate the density decrease from theLISM side to the SW side of the HP, it is necessary either to fit the HP surface ensuringthe satisfaction of the boundary conditions suitable for tangential discontinuities in MHD,or use AMR. Izmodenov & Alexashov (2015) show that such boundary layer exists also onthe SW side of the HP. It is interesting to note in this connection that Belov & Ruderman(2010) argue that the density jump across the HP may disappear at the LISM stagnationpoint, the density variation being smooth and occurring mostly in the boundary layers.It is seen from Figs. 2 that the LISM plasma density reaches its maximum at 50–100AU from the HP surface, which is of the order of 1–2 charge exchange free paths in theOHS. The maximum values are: (1) 0.2; (2) 0.175; (3) 0.155; (4) 0.14; (5) 0.115; and (6)0.095 cm − , respectively. Further on, the plasma density is only decreasing in the sunwarddirection, and the maximum gradient is at the HP surface. Since we determine the HPposition exactly, by solving a level-set equation for the boundary between the SW andthe LISM (Borovikov et al. 2011), we can look closer at the magnetic field and densityvariations in the HBL. In Fig. 4, we show the magnetic field vector magnitude, B , and itselevation and azimuthal angles, δ and λ , along the V1 trajectory in the simulation with B ∞ = 2 . µ G from Table 1. In agreement with observations, both B and λ are continuousacross the heliopause. The grid resolution near the HP is 1.2 au cubed. The HP position isshown with the vertical dashed line. The numerical smearing of quantities is about ± δ increases to about 28 ◦ . The observed values are: δ = 22 ◦ ± ◦ and λ = 291 ◦ ± ◦ (Burlaga & Ness 2014a). In addition, the simulation shows a consistent undraping of theISMF as V1 propagates deeper into the LISM. The gradient in δ and λ are consistentwith observations reported by Burlaga & Ness (2014b). However, these gradients becomesmaller if averaged from the crossing time to 2016 (Burlaga & Ness 2016). This leads us toa conclusion that time-dependent phenomena, such as described in Fermo et al. (2015) arelikely to affect the undraping.Figure 5 shows the plasma density distribution in the meridional plane ( right panel )and along the V1 trajectory ( left panel ) in a time-dependent simulation which particularlyfocuses on the heliopause resolution. The LISM parameters are taken from Table 1( B = 3 µ G), but a nominal (periodic with the 11-year period) solar cycle is taken intoaccount, similarly to Pogorelov et al. (2009a) and Borovikov & Pogorelov (2014). Thelength scale has been decreased by a factor of 1.1 on the left panel, to assure the HPposition in the observed point. The density increase with distance from the heliopauseshould result in the increase of the electron oscillation plasma frequency, in accordance with V1 observations with the Plasma Wave Instrument (PWS) (Gurnett et al. 2015) (see alsoFig. 6). Initially, the density increased from 0.06 cm − to 0.08 cm − from from Nov 2012 toApril-May 2013 (on the distance of ∼ V1 ).) In the figure, this distance issomewhat greater ( ∼ . V1 occurred inNovember 2014 and showed the density increase to about 0.09–0.11 cm − . The spacecrafttraveled approximately 7 AU between these events. The simulation shows the distance ofabout 10 AU. Two more plasma oscillation events were measured by PWS: on Sep–Nov2015, when the density measured on Day 298, 2015 was 0.115 cm − at a distance of about133 AU from the Sun, and on Day 201, 2016, when the density was determined to be0.113 cm − at a distance of 135.5 AU. We show the latter event in Fig. 5, where it occurred 14 – ∼
139 AU. In addition, the measured density change between the latest two events wasalmost negligible, which is not seen in our simulation. These discrepancies should not besurprising because our model is not detailed enough to identify MHD shocks propagatingthrough the OHS due to realistic time-dependent boundary conditions. The figure doesshow, however, the perturbations propagating through the LISM due to the solar cycle. Itcan be seen that there are intervals as large as 5–6 AU without any density increase. Thedensity gradient becomes smaller as V1 continues to traverse the LISM and it may take ∼
4. Instabilities and magnetic reconnection near the heliopause
As mentioned above, the HP instability may be responsible for the HP “structure”observed by V1 before it entered the LISM completely. The HP is also a likely venuefor magnetic reconnection. It is known from both theory and simulations that the HP isunstable both at its nose and on the flanks (see Ruderman & Belov 2010; Pogorelov et al.2014, 2017; Avinash et al. 2014, and references therein). Charge exchange between ionsand neutral atoms play an important role here through the action of the source terms inthe momentum and energy equations. Near the stagnation point, where the shear betweenthe SW and LISM flow is small, charge exchange results in a sort of Rayleigh–Taylor (RT)instability (Liewer et al. 1996; Zank 1999; Florinski et al. 2005; Borovikov et al. 2008b),which is known to take place when a heavier fluid lies upon a lighter one. Farther fromthe stagnation point, the Kelvin–Helmholtz (KH) and other instabilities may develop(Ruderman 2015). The instabilities of the HP nose are seen only at high numericalresolution in this region ( ∼ . B , the y -component of the magnetic field vector, B y , and plasma density.The boundary conditions are taken from Borovikov & Pogorelov (2014), but the resolutionis higher (0.22 au cubed). For better understanding of the solution behavior at Voyager spacecraft, the cross-cuts are made by the plane defined by the V1 (in the northernhemisphere) and V2 (in the southern hemisphere) trajectories. One can see that the solarcycle creates magnetic barriers of opposite polarity that propagate through the IHS towardsthe HP. The HMF polarity in such barriers changes every 11-years. As a barrier approachesthe HP, it becomes exceedingly thinner (see Fig. 8), creating the possibility of magneticreconnection across the HP if the orientations of B become suitable. The inspection ofthese figures demonstrates that this is especially true for the southern hemisphere, where 16 – V2 is approaching the HP. As seen in Figs. 8d–e, magnetic reconnection reveals itself as atearing mode (or plasmoid) instability. Similar features are seen in numerical modeling ofmagnetic reconnection during solar eruptions presented, e.g., in Pontin & Wyper (2015).It has been shown that plasmoid (tearing-mode) instability of extended currentsheets provides a mechanism for fast magnetic reconnection in large-scale systems.Within an MHD framework, the instability has a growth rate that increases with theLundquist number, while its nonlinear development results in a turbulent reconnectionlayer and average reconnection rates that are independent of or weakly dependent onresistivity (Shibata & Tanuma 2001; Loureiro et al. 2007, 2013; Bhattacharjee et al. 2009;Uzdensky et al. 2010; Higginson et al. 2016, 2017). The conditions for such instabilityare satisfied for high Lundquist numbers, which can be reached by increasing the gridresolution, and large aspect ratios of the reconnecting current layers. While we do notexplicitly include resistivity, magnetic diffusion proportional to ∆ enters our system dueto discretization. Henceforth, both conditions are clearly satisfied in our simulations.While plasmoid instability ensures that reconnection can remain fast at very largeLundquist numbers, it is by no means the only such mechanism. Simulations andtheoretical considerations demonstrate that in the presence of turbulence large-scalemagnetic reconnection can proceed with large, resistivity-independent rates (e.g.,Lazarian & Vishniac 1999; Eyink et al. 2011, and references therein). In fact, Beresnyak(2017) argues that the physical reason for the resistivity-independent reconnection rateis a consequence of turbulence locality, similar to models of reconnection due to ambientturbulence. Furthermore, under certain conditions, the plasmoid instability can directlytransition the system to a kinetic regime where local reconnection rates again becomeindependent of resistivity (see, e.g., Daughton et al. 2009; Daughton & Roytershteyn 2012;Ji & Daughton 2011). In addition, Zweibel et al. (2011) show that the presence of neutral 17 –atoms may modify the reconnection process. These effects, however, are beyond the scopeof this paper and will be addressed elsewhere.By tracing magnetic field lines that pass through one of the plasmoid regions shownin Fig. 8, we arrive at another conclusion: the actual magnetic reconnection occurs at adistance of a few AU away from this region. This conclusion has a far-reaching consequence,i.e., magnetic reconnection events have global, macroscopic consequences, which cannot beaddressed directly by kinetic simulations because of the length scale limitations intrinsicto them. A similar situations may be observed in solar flares (Liu et al. 2013), where thelength of a magnetic reconnection sheet is in excess of 10 ion inertial lengths.While more “reconnection” is seen in the southern hemisphere and at V2 , theconsequences of the HP instability are stronger at V1 . Magnetic field distributions inFigures 7–8 demonstrate the possibility that V1 could cross the regions belonging to theSW and LISM consecutively. This means that on the way out of the heliosphere it couldbe magnetically connected either to the HMF, and observe enhanced fluxes of anomalouscosmic rays (AMRs) and depressed GCR fluxes, or to the ISMF, where ACRs virtuallydisappear, while the GCR flux increases. This scenario requires that diffusion parallel tothe magnetic field should be substantially greater than perpendicular diffusion. Luo et al.(2015) show that the abrupt increase in the GCR flux observed by V1 when it crossed theHP is possible only if the ratio between the parallel and perpendicular diffusion coefficientsexceeds 10 . Our simulations provide a plausible explanation of the changes in the ACRand GCR fluxes before V1 entered the LISM permanently.To supplement the results shown in Figs. 7–9, we add animations of the same quantitiesto the on-line version of the paper. 18 –
5. Magnetic field dissipation in the IHS
Issues related to the magnetic field behavior at V1 and V2 are of great importancebecause both spacecraft provide us with appropriate measurements (Burlaga & Ness 2014a;Burlaga et al. 2014). In the idealized simulation considered in the previous section, theangle between the Sun’s magnetic and rotation axes is a periodic function of time. Theminimum tilt of 8 ◦ is attained at solar activity minima, whereas the maximum of 90 ◦ isreached at solar activity maxima, where the magnetic dipole flips from one hemisphereto another. This is, of course, a simplification. Pogorelov et al. (2013c) considered solarcycle effects with the tilt being a function of time from WSO data. As a result, themagnetic barriers described in the previous section had a layered structure, which was dueto local non-monotonicities in the tilt angle in the vicinity of the spacecraft latitude. Everynon-monotonicity of this kind creates an additional current sheet.Clearly, resolving the sectors of alternating magnetic field polarity in the IHS isimpossible, even for a simplified solar cycle. This is because the sector width is proportionalto the SW velocity, provided that the HCS propagates kinematically and exerts no backreaction onto the plasma surrounding it. Borovikov et al. (2011) proposed another approachto track the HCS surface. Other approaches to track the HMF polarity were used inCzechowski et al. (2010); Alexashov et al. (2016). In the approach of Borovikov et al.(2011), the HMF is assumed unipolar and a special, level-set equation is solved to propagatethe HCS surface from the inner boundary towards the HP. Once the HCS surface is known,it is easy to assign proper signs to the HMF vector components at a postprocessing stage.However, this turned out to be impossible even for the level-set approach because thedistances to be resolved near the HP become too small for any practically acceptable grid,so the HCS was accurately resolved only half way from the TS to the HP. In principle, forany chosen grid resolution ∆, this approach fails once ∆ /T > V , where T ≈
25 days is the 19 –period of the Sun’s rotation. E.g., for the radial SW velocity component of ∼
90 km s − ,which is currently observed at V2 the sector width in the solar equatorial plane is ∼ . ≈ .
13 AU to resolve the sector structure.Such resolutions are impossible except for over a very limited region.
Voyager 1 , on theother hand, had been observing negative radial velocity components for two years beforeit crossed the HP (Decker et al. 2012). The sector width should be negligible in this case.Moreover, the sector width decreases to zero at the HCS tips (see Fig. 14 in Pogorelov et al.2013c), which makes attempts to resolve the traditional HCS structure questionable. Wecall the HCS traditional if the sector structure is entirely due to the Sun’s rotation with afixed period.The intervals between HCS crossings depend on the relative velocity of the SW withrespect to the moving spacecraft. If V = 350 km s − in front of the TS and becomes150 km s − behind it, the maximum sector width was about 0 . × V AU, i.e, 4.9 AUin front of the TS and 2.1 AU behind it. The velocities of V1 and V2 are approximately16.6 km s − and 14.2 km s − , respectively. So V1 should have been crossing an idealizedHCS every 25.5 days before the TS and every 27 days after it, which is not the case (seeFig. 10). Voyager V1 (Fig. 10) and V2 (Fig. 11) spacecraft. It is especially interesting that there is no HCS crossings at V1 forat least 100 days after crossing the TS. The HCS is not crossed if a spacecraft moves withthe velocity of the ambient SW (17.1 km s − for V1 and 15.7 km s − for V2 ). As shownabove, the expected “nominal” decrease in the HCS crossing time immediately beyond theTS is small. This means that either the HMF strength is too small in front of the TS,so that polarity reversals are caused not only by the HCS crossings but also by turbulentfluctuations in the SW plasma, or time-dependent phenomena make the sector structureirregular. Notice that the sectored region of the SW plasma turns northward in the IHS, sothe spacecraft should remain in the polarity-reversal region.The sector boundaries, if not destroyed by turbulence should be piling up in front ofthe HP moving inward (at V1 before it crossed the HP) at 20 km s − , so the the numberof sector crossings should have increased dramatically, but it had not. Our numericalsimulations show that the absence of magnetic field polarity reversals observed by V1 nearthe HP may be due to its entering a magnetic barrier. The following few polarity reversalsmay be caused by the complicated structure of the HP caused by its instability.The extended periods with no polarity reversals are also seen along the V2 trajectory.Richardson et al. (2016) have investigated the effect of the magnetic axis tilt on the numberof HCS crossings and compared the observed and expected numbers. It has been reportedthat the number of HCS crossings substantially decreased two years after V1 and V2
22 –crossed the TS. However, V2 might have entered the unipolar region at that time. It wasultimately concluded that there are indications of magnetic field dissipation possibly dueto magnetic reconnection across the HCS. However, occasional deviations between theobserved and expected HCS crossings are to be expected also for the reasons of streaminteraction discussed above and should not be necessarily interpreted as clear evidencefor magnetic reconnection. On the other hand, as shown by Drake et al. (2017), V2 datareveal that fluctuations in the density and magnetic field strength are anticorrelated inthe sectored regions, as expected from their magnetic reconnection modeling, but not inunipolar regions. A possible annihilation of the HMF in such regions may also be anexplanation of a sharp reduction in the number of sectors, as seen from the V1 data.Richardson et al. (2016) assumed that the radial and latitudinal velocity componentsat the boundary between the unipolar and sectored regions are determined by spacecraftmeasurements, i.e., are independent of latitude. However, this is not quite true. Figures 12show that the variations in the latitudinal components can be substantial, which is notsurprising because the boundary of the sectored region should propagate from ∼ ◦ duringsolar minima to 90 ◦ at solar maxima in 11 years. Moreover, we remember that a layer ofthe sectored magnetic field never disappears on the inner side of the HP surface, at leastabove the equatorial plane. Only its width is a function of time. Thus, a more detailedanalysis of observational results may be required.While the extent to which the HMF dissipates in the IHS remains the subject ofinvestigation, numerical simulations allow us to find out what happens to B if the HMF isassumed to be unipolar (see the discussions, e.g., in Pogorelov et al. 2015, 2017). Figure 13shows the magnetic field magnitude, B , and the spherical components of B along the V1 trajectory for the simulation with B ∞ = 3 µ G from Table 1. It is clear from this simulationthat the calculated magnetic field strength is substantially overestimated (see Fig. 10). 23 –This behavior of the modeled HMF at small distances beyond the TS was noticed earlier byBurlaga et al. (2009), but attributed to possible transient effects. On the other hand, basedon
Ulysses measurements, solar-cycle simulations in Pogorelov et al. (2013c), which takeinto account the observed variations in the magnetic axis tilt with respect to the rotationaxis, although not resolving the HCS, are able to reproduce V1 observations relatively wellon the average. Thus, the possibility of an occasional HMF annihilation in the IHS shouldnot be disregarded.
6. Conclusions
In this paper we have addressed a variety of physical phenomena related to
Voyager observations. Some of these phenomena have clear physical explanation, whereas additionalinvestigations, both theoretical and numerical, are necessary to describe the others. Themodification of bow shocks by charge exchange between ions and neutral atoms is awell-known phenomenon, frequently referred to not only in the heliospheric bow wavecontext, but also in astrophysics (Chevalier & Raymond 1978; Chevalier et al. 1980;Blasi et al. 2012; Morlino et al. 2012, 2013; Morlino & Blasi 2016). Charge exchangecannot modify the Hugoniot relations at shocks (essentially the conservation laws of mass,momentum, and energy) propagating through plasma, but can substantially change theplasma properties ahead of and behind the shock. As a result, secondary H atoms of theSW origin propagate far upstream into the LISM and decrease the bow shock intensity ascompared with ideal MHD flows. The structure of the bow wave in front of the heliopauseis of importance for the interpretation of
IBEX and
Voyager observations. In addition,Schwadron et al. (2014); Zhang et al. (2014), and Zhang & Pogorelov (2016) demonstratethat it also affects the observed anisotropy of 1–10 TeV cosmic rays. In the situationrelevant to the heliospheric bow shock, the range of possible ISMF strengths, B ∞ , is 24 –such that the influence of charge exchange becomes dominant, i.e., the contribution ofthe shock compression to the total density (and magnetic field) enhancement on the HPsurface is small. We have shown that any attempt to predict the presence of a shock in thecompression region is impossible only by analyzing the properties of the LISM not affectedby the presence of the HP. In the absence of such shock, the LISM interaction with theheliosphere produces a rarefaction wave propagating outwards into the LISM.High-resolution simulations show the presence of a HBL (a region of depressed plasmadensity and increased magnetic field strength) on the LISM side of the HP. The identificationof the internal structure of such layers is a challenge for discontinuity-capturing numericalmethods because of a dramatic change of plasma density across the HP. Discontinuity-fittingmethods like that of Izmodenov & Baranov (2006) are more suitable for this purpose,unless substantial adaptive mesh refinement is applied. A drawback of HP-fitting methodsin the difficulty to address related physical instabilities. We demonstrated that thedensity increase with distance from the heliopause is consistent with the plasma wavefrequency observations at V1 (Gurnett et al. 2015). It is of interest that the boundarylayers seen in global simulations are not due to magnetic field effects, since they are alsopresent in simulations without magnetic field (Baranov & Malama 1993; Zank et al. 1996).Comparison of multi-fluid and ideal MHD simulations performed by Pogorelov & Zank(2005) suggest that the “depth” of the boundary layer increases with the LISM neutralH density. On the other hand, the width of the boundary layer seems to be comparablewith the charge exchange mean free path in the LISM (40–50 AU), which means that thisboundary layer is somehow related to change exchange. The relative contribution of theplasma pressure anisotropy on the HBL structure requires further investigation.While there is little doubt that HP instabilities and magnetic reconnection areintrinsic to the heliospheric interface, it remains a challenge to relate observational data to 25 –simulation results. This is because observations are limited to one point per time, whichmakes it difficult to put them into the context of large-scale, 3D phenomena occurringnear the HP. Indeed, the HP instability may result in substantial penetration of the LISMplasma into the SW. This, in turn, creates the possibility for a spacecraft like V1 tocross the LISM and SW plasmas several times consecutively, which may be a reason forthe changes in the observed ACR and GCR flux intensity while V1 was crossing the HP.While this scenario still requires confirmation from a direct simulation, it is clear thata 3D, data-driven model of the SW–LISM interaction is necessary to explain the ISMFbehavior along the V1 trajectory. A few attempts to create such model have been presentedrecently by Fermo et al. (2015) and Kim et al. (2016). In addition to the HP instability, thesimulations presented here predict that V2 should observe more consequences of magneticreconnection near the HP than V1 did. This is apparently the result of the HMF–ISMFcoupling at the HP. Further investigation is necessary to understand the physics of plasmawave generation and radio emission observed by V1 in the OHS. Such investigation shouldalso be data-driven because the ISMF undraping observed by V1 strongly suggest thattime-dependent phenomena are deeply involved in this process. Steady-state MHD-kineticsimulations presented here (see also Zirnstein et al. 2016) show that such undraping shouldbe monotone.The HCS, its behavior in the IHS, and possible influence on the magnetic field andplasma distributions has been one of the most controversial subjects of discussion in thepast few years (for a detailed analysis, see Pogorelov et al. 2017). It is clearly impossibleto resolve micro-scale phenomena related to the HCS in an ideal MHD model, especiallyclose to the HP. On the other hand, kinetic simulations of magnetic reconnection in theIHS and across the HP are too local to be able to describe the macroscopic effect of thisphenomenon. While the HCS is an inherent component of magnetic field distribution inthe IHS, numerical simulations allow us to perform a thought experiment where the HCS 26 –is excluded by assuming the HMF to be unipolar. As discussed above, one may try toassign the signs to the magnetic field components post factum , after a unipolar simulationis finished, provided that we can track the HCS surface as a discontinuity kinematicallypropagating towards the HP with the SW velocity. We demonstrated here that thisapproach is not well justified. If the HMF is assumed unipolar, the simulated magneticfield strength is considerably higher than it was measured by V1 and V2 . This is a possibleexplanation of the discrepancies in the simulations performed with the unipolar and dipolarHMF presented by Opher et al. (2015); Pogorelov et al. (2015, 2016) and further discussedin Pogorelov et al. (2017).In summary, our results imply that there is some dissipation of HMF in the IHS.There many reasons for such dissipations: (1) SW turbulence, which is especially enhancedbeyond the TS; (2) magnetic reconnections; (3) kinetic and MHD instabilities, etc. Afew evidences of such dissipation are discussed in this paper (stochastic destruction of theHCS in the IHS and tearing-mode instability destroying time-dependent magnetic barrierswhen theirs aspect ratio becomes small). Pogorelov et al. (2013b) showed that this mayresult in additional plasma heating and changes in the SW radial velocity componentgradients. However, as far as the plasma heating is concerned, it should be clear thatits analysis is impossible without proper treatment of PUIs and anomalous cosmic raysfrom the TS into the IHS. Pogorelov et al. (2016) demonstrated that specifically designedboundary conditions for PUIs at the TS may be able to describe the preferential heatingof PUIs as compared with thermal SW ions (Zank et al. 2010). Such boundary conditionsare of kinetic nature and therefore cannot be derived from any continuum mechanicsapproach. Approaches which are not based on the conservation-law principles and involvingstraightforward calculations of the PUI pressure derivatives across the TS (e.g., Usmanov2016) are mathematically flawed. Fahr & Siewert (2013) (see also references therein)proposed a number of theoretical approaches to derive the above-mentioned boundary 27 –conditions. Local particle simulations may also serve as tools to help derive the boundaryconditions that would be able to reproduce the ion acceleration at any point of the TS. Theion distribution function in the shock vicinity is highly anisotropic, which makes continuummechanics approaches not applicable. We note also a recently proposed generalized systemof such equations that take into account dissipative affects and heat flux terms (Zank et al.2014), but not taking into account possible reflection of PUIs and their further accelerationat the TS.This work was supported NASA grants NNX14AJ53G, NNX14AF43G, NNX15AN72G,and NNX16AG83G. The work done in the University of Iowa was supported by NASAthrough Contract 1279980 with the Jet Propulsion Laboratory. This work was alsopartially supported by the IBEX mission as a part of NASA’s Explorer program. Weacknowledge NSF PRAC award ACI-1144120 and related computer resources from the BlueWaters sustained-petascale computing project. Supercomputer time allocations were alsoprovided on SGI Pleiades by NASA High-End Computing Program award SMD-16-7570and Stampede by NSF XSEDE project MCA07S033.The authors are grateful to G. P. Zank for stimulating discussions. 28 –Fig. 1.— SW–LISM interaction pattern in the presence of solar cycle effects. The y -component (top panels) of the magnetic field vector and its magnitude (bottom panels)in the meridional (left panels) and equatorial (right panels) planes. Distances are given inAU and magnetic field in µ G. One can see the TS, HP, and bow wave. 29 –Fig. 2.— Distributions of the plasma density (black lines) and fast magnetosonic Machnumber (red lines) along the x -axis in the simulations from Table 1. The vertical bluedashed lines shows the HP position. 30 –Fig. 3.— Distributions of the plasma density in the meridional plane for the simulationsfrom Table 1. 31 –Fig. 4.— Distributions of the magnetic field vector magnitude B (black line) and its elevation(blue line) and azimuthal (red line) angles, δ and λ , along the Voyager 1 trajectory. 32 –Fig. 5.— The distribution of plasma density (left panel) along the V1 trajectory and itscomparison with the plasma waves events detected by the spacecraft beyond the heliopause,and ( right panel ) in the meridional plane.Fig. 6.— Voyager V1 and V2 trajec-tories in the simulation that takes into account solar cycle effects. Time is increasing fromthe left to the right and from the top to the bottom. 34 –Fig. 8.— The same as in Figure 7 but for B y . 35 –Fig. 9.— The same as in Figure 7 but for plasma density. 36 –Fig. 10.— Hourly averages of the magnetic field vector magnitude and azimuthal angle alongthe V1 trajectory ( Voyager data courtesy of CohoWeb). 37 –Fig. 11.— Hourly averages of the magnetic field vector magnitude and azimuthal angle alongthe V2 trajectory ( Voyager data courtesy of CohoWeb). 38 –Fig. 12.— Time-dependent distributions of the latitudinal component of the velocity vectorin the meridional plane in the simulation of Pogorelov et al. (2013c) demonstrate that thiscomponent strongly depends on latitude. 39 –Fig. 13.— The distribution of the magnetic field strength (black line) and its radial (blueline), θ - (red line), and φ -components (green line) along the V1 trajectory in an MHD-kinetic simulation where the heliospheric magnetic field is unipolar. Presented simulationsdemonstrate that the assumption of unipolar HMF results in a considerably overestimationof the magnetic field. 40 – REFERENCES