A Turbulent Heliosheath Driven by Rayleigh Taylor Instability

The heliosphere is the bubble formed by the solar wind as it interacts with the interstellar medium (ISM). Studies show that the solar magnetic field funnels the heliosheath solar wind (the shocked solar wind at the edge of the heliosphere) into two jet-like structures 1-2 . Magnetohydrodynamic simulations show that these heliospheric jets become unstable as they move down the heliotail 1,3 and drive large-scale turbulence. However, the mechanism that produces of this turbulence had not been identified. Here we show that the driver of the turbulence is the Rayleigh-Taylor (RT) instability caused by the interaction of neutral H atoms streaming from the ISM with the ionized matter in the heliosheath (HS). The drag

Astrophysical jets are prone to instabilities such as kink or pinch modes 5 . Global simulations of the heliosphere have revealed that the heliospheric jets also drive instabilities, allowing the jets to evolve into a complex turbulent outflow with macroscales that are 100's of au (astronomical units) in size. The scale of this turbulence dwarfs that which develops elsewhere in the heliosphere. While it was suggested that the kink instability was the driver of the turbulence 5 initiated by plasma jets, there has been no definitive exploration of jet instabilities in partially ionized flows nor an exploration of the impact of this turbulence on particle energization. The goal of the present paper is to establish the driver of heliospheric jet instabilities. We carry out a variety of simulations, including several with zero flow velocity in the interstellar medium (ISM), to identify the driving mechanism. We first establish that the kink instability is not the driving mechanism, finding in fact that the heliospheric jets are stable to the kink instability.
The computational model has a single ionized component and 4 neutral components as in ref 1. The cold thermal solar wind and hot pickup ions (PUIs) are treated as a single species 1,6 .
The neutral H component is captured with a four-fluid approximation 6,7 . In all simulations we use the same inner boundary conditions (see Methods). We set the interstellar magnetic field (BISM) in the z (for the coordinate system see Methods) direction to suppress the Kelvin-Helmholtz instability. We cover all regions of interest with a high grid resolution and run timedependent fashion cases for hundreds of years of simulated time. For comparison, it takes one year for the solar wind to reach the termination shock (TS) (~100 au) at a velocity of 400 km s −1 .
In the heliosheath (HS), the flow speed is ~50 km s −1 and it takes ~40 yrs to traverse the heliosphere, ~400 au. The details are given in Methods.
In the absence of neutral H streaming, the heliospheric jets are stable (Methods - Figure   1), so the kink instability is not active. Even if the ionized ISM is streaming, if there were no neutrals then the jets are stable (see Figure 3A). When streaming of neutral H is included ( Figure   1), an instability develops along the axis of the heliospheric jets. Due to charge exchange between the neutral H atoms and the plasma, even if there were no initial plasma motion in the ISM, a wind forms in the neutral flow direction that bends the heliospheric magnetic jets tailward. The instability is shown by the blue color ( Figure 3a) where the magnetic field intensity is close to zero. This instability is the equivalent of the Rayleigh-Taylor (RT) instability, first proposed to occur along the heliopause (HP) 8-10 due to the effective gravity produced by neutral atoms charge exchanging and therefore dragging the plasma ions (Methods - Figure 3B). Later works [11][12] expanded the analysis to include the stabilizing effect of hot HS neutrals. The predicted RT growth timescale was ~3 yrs at the HP for typical values of the ISM and HS plasma when no interstellar magnetic field (BISM) was present. When BISM is included in the models the RT at the HP is less unstable and the growth timescale is 12 to 50-100 yrs with wavelength scales ~ 40 au.
In Methods we revisit the RT instability, examining two classical cases: a single interface between high and low density fluids that is the analog of the HP; and a high density fluid slab surrounded by low density fluid that is the analog of the situation along the axis of the heliospheric jets. The RT instability in fluids can be derived analytically considering a sharp interface between two fluids with densities " and # (ref. 13) (see Methods - Figure 3A). At the HP, the jump in density is "~ 10 # ; " being the ISM density and # the density in the HS (see Methods). For the case with no magnetic field all wavelengths are unstable. The growth rate 8,10 , which is the inverse of the growth timescale, is = * * ( " − # ) ( # + " ) ⁄ , where k is the wave number and g * is the effective gravity (see Methods Eq. (7)). At the HP, the growth rate The confinement of the heliospheric flow by the solar magnetic field creates a profile of the plasma thermal pressure and density that drops from near the jet axis towards the edges of the heliospheric jets 2 .
To understand the RT instability for the heliospheric jets, one needs to generalize the RT instability for a case where the increase in density occurs between two surrounding regions of lower density (Methods- Figure 3D). The growth rate for a slab of width A with high density " , between two surrounding regions of lower density # (Methods-Eq. (18)) is # = * 3 4 5 (6 7 86 9 ) 9 (6 7 :6 9 4 5 )(6 7 4 5 :6 9 ) where < = tanh <A # . This growth rate is smaller than the growth rate for a single density jump because the jumps on the stable and unstable sides compete with each other.
To ascertain that the instability is in fact an RT instability, we ran the same model as in Figure 1 but varied the density of the neutrals in the pristine ISM. As nH is increased from 0.045 to 0.36 cm -3 the instability becomes stronger (Figure 2). Along the axis of the heliospheric jets the fluid is heavier than the surrounding fluids with a jump of "~ 2 # ( Figure 2E, Methods-   The model used in this paper was based on a single fluid plasma approximation (i.e., the cold thermal plasma and hot pick-up ions (PUIs) were considered as one fluid). More advanced models such as ref. 16 that separate the cold solar wind from the hot PUI ions and treat the neutrals in a kinetic fashion 15 are much more computationally challenging and will be explored in the future. A two-fluid model treatment might affect the evolution of the instability, although the instability was also seen using this more advanced model 16 .
The turbulence triggered by the RT instability is macroscopic and has potentially important implications for particle acceleration. The expected acceleration of ions driven by the turbulence will be calculated in future studies by the extension of models such as kglobal 17-18 that include self-consistently the interaction between electrons and ions with large-scale dynamics of plasmas turbulence.
Until a few years ago the shape of the heliosphere was thought to be comet-like with a tail extending to 1000's of au 19,20 . Currently there is an active debate, propelled by recent modeling and observations, about the shape and structure of the heliospheric tail 19-23;1,16 . The RT instability along the jets seems to be the responsible for opening the heliospheric tail, bringing the northern lobe towards the southern lobe, and triggering reconnection that further opens the heliospheric tail. This result therefore explains why models 24 that maintain a perfectly ideal heliopause and do not allow exchange of solar plasma with that of the interstellar medium show confinement by the solar magnetic field with a long comet-like tail.
The mixing of solar material with the ISM via this mechanism is important, it may provide a general process whereby stellar material may be mixed into a partially ionized local ISM.
The signatures of large-scale RT turbulence in Energetic Neutral Atom (ENA) maps should be explored as well. Figure 4 shows that high density nodules form near the axis of the heliospheric jets when the turbulence is fully developed. The nodules are large in size and will vary with time as the RT instability evolves. The formation and evolution of the nodules will be a subject of a future paper. We expect that high resolution ENA maps from future missions such as IMAP 4 would be able to detect these nodules.

Data availability
Our  Figure 3 at t=396 years. The heliopause is shown as a yellow surface defined by the temperature isosurface at T=85,000K that captures the closed solar field lines 5 . The white lines are ISM streamlines that pierce the heliospheric tail. The red structures are nodules of high density formed along the heliospheric jet axes as a result of the RT instability. They are defined as iso-surfaces of density of 0.019 cm -3 .

Description of the Rayleigh-Taylor Instability at a single density interface: applicable to the heliopause:
The Rayleigh-Taylor (RT) instability in fluids can be derived considering an interface between two fluids with uniform densities " and # (see Figure M1). Suppose a high-density fluid " is adjacent to a lower density fluid # in a "gravitational" field with acceleration g * , the equation for a sharp discontinuity with the gravity pointing in the x direction is Considering a single discontinuity at x=0; at ≠ 0 Eq. (1) is written as Making sure that g as a function of s is continuous; E f fg g v 8w w = * < 9 u 9 ( # − " ) g (6) And the growth rate is # = 6 7 86 9 6 7 :6 9 * (7) The neutral-ion drag acts as an effective gravity * xxxx⃗ = OP ( ⃗ P − ⃗ T ) where OP is the charge-exchange collision frequency, ⃗ P is the bulk neutral velocity and ⃗ T is the bulk ion velocity. We will consider here only the dominant species of neutrals that come from the ISM. Ref. 10 explored the effect of including the HS neutrals in the RT at the HP. Our numerical MHD model presented here also includes the HS neutrals. The ion bulk velocity ⃗ T must flow parallel to the density interface, so the perpendicular component of the gravity that induces the RT instability simplifies to * = OP P . The collision frequency OP = S * , where is the charge exchange cross section, S is the neutral density and * = 3h c C k ( " # + # # ) + (∆ ) # is the effective velocity difference that includes ∆ = " − # is the relative speed between species 1 and 2 and the thermal speeds " and # for the two species. At the HP " is the ISM density, # is the HS plasma density and "~ 10 # . The growth rate from Eq. (7) is ~3 { "" * ~* * .

Rayleigh-Taylor instability for two density jumps: applicable to the heliospheric jets
For the case of the heliospheric jets, one needs to generalize the RT instability for a case where the increase in density occurs between two surrounding regions of lower density (Methods- Figure 3D). Based on the simple RT calculation, one expects that one of the interfaces is stable ("L") and the other one is unstable ("R" ), however we need to investigate the stability of the slab as a whole.
Consider a slab of high density " of a width A (Methods- Figure 3D) between two regions of lower density # . The equation for a sharp boundary Eq. (1) For x> A/2 the perturbation; g = : 8<hg8 } 9 k ; and for x<A/2 g = 8 <hg: } 9 k and for A/2<x<A/2, we take g = O cosh( ) + € sinh( ). One needs to make sure that g is continuous across = ± /2, so Defining a new variable ds=dx/ and taking # = − # as the square of growth rate; Making sure the g as a function of s is continuous; at "R" edge; = + /2, integrating: Numerical Model. The coordinate system is such that the z axis is parallel to the solar rotation axis, the x axis is oriented in the direction of the interstellar flow (which points 5° upward in the x-z plane) and the y axis completes the right-handed coordinate system. Inner boundary: For all models presented in the paper the parameters of the solar wind at the inner boundary at 30 au are the same as used in our previous works 2 : vSW = 417 km/s, nSW = 8.74 x 10 -3 cm -3 , TSW = 1.087 x 10 5 K (OMNI solar data; http://omniweb.gsfc.nasa.gov/). The magnetic field is given by the Parker spiral magnetic field with BSW = 7.17x10 -3 nT at the equator. We use a monopole configuration for the solar magnetic field. This description, while capturing the topology of the field lines, does not capture the change of polarity with solar cycle or across the heliospheric current sheet. This choice, however, minimizes artificial reconnection effects, especially in the heliospheric current sheet. We assume that the magnetic axis is aligned with the solar rotation axis. The solar wind flow at the inner boundary is assumed to be spherically symmetric. Figure 3 Cartoons of (A) ideal situation of a heavier fluid "1" sitting atop of a lighter fluid "2" with "gravity g" directed perpendicular to the interface; (B) RT at the HP: the streaming of neutrals from the ISM acts as an effective gravity g * between the dense ISM and the lighter fluid in the HS; (C) RT along the axis of the heliospheric jets in the situation where there is no motion through the ISM. The streaming of neutrals perpendicular to the axis of the heliospheric jets acts as an effective gravity between the denser fluid at the center of the jet (in gray), in the HS compared to the lighter fluid surrounding it. (D) Cartoon of RT in slab geometry. In isolation, the left density jump "L" would be stable to RT, while the right interface "R" would be unstable.
ISM conditions: Methods- Figure 1: vISM =0; , nISM = 0.06 cm −3 and TISM = 6,519 K and no neutral streaming. The strength of the BISM in the model is 4.4 G=0.44nT in the z direction, parallel to the axis of rotation of the sun to stabilize Kelvin-Helmholtz instabilities. The simulation was run to 526 years. Figure 1: same ISM conditions as in Methods- Figure 1 but including neutral hydrogen atoms in the interstellar medium with number density nH = 0.18 cm −3 ; the neutral temperature is the same as for the interstellar plasma. The speed of the neutral H atoms was set to 26.4 km/s in the x direction. Figure 2: All cases use the same ISM conditions as in Methods - Figure 1 but with different number densities of hydrogen atoms in the interstellar medium. The temperature of the neutral H atoms is the same as for the interstellar plasma. The speed of the neutral H atoms was set to 26.4 km/s in the x direction in all cases. For panels A and E: nH = 0.36 cm −3 , panels B and F same as in Figure 2; panels C and G nH = 0.09 cm −3 and panels D and H nH = 0.045 cm −3 . Figure 3: This case uses realistic ISM conditions; vISM = 26.4 km s −1 , nISM = 0.09 cm −3 (the higher than the usual density is used to keep ISM pressure the as when we include a interstellar magnetic field BISM) and TISM = 6,519 K. The number density of hydrogen atoms in the interstellar medium is nH = 0.18 cm −3 and the velocity and temperature are the same as for the interstellar plasma. BISM was not included to avoid reconnection. The model was run to 872 years with no charge exchange between the neutral H and the ionized component. At 872 years charge exchange was turned on and the model was run to 1307 years.
Grid resolution: For the model and cases presented in Figures 1-2 and in Methods- Figure 1; the computational domain extends to ±3000 au, in the x, y and z directions. The smallest grid cell size is 1.47 au inside x=y= ±105 au and z=± 608 au and the coarsest grid has 187.5au resolution near the outer boundaries. For Figure 3: the computational domain extends to ±1500 au in the x, y and z directions. The smallest grid is 0.37 au near the inner boundary and z=± 608 au and a coarser grid of 93.75 au is used near the outer boundaries. This case resolved the entire heliosphere with high grid resolution, including the heliospheric tail (with resolutions of 3 au throughout the tail for case A and 2 au for case B until 400 au and 4 au until 600 au).