© ESO 2020

1. Introduction

Following the Voyager-1 and -2 crossings of the solar wind termination shock, the scientific community somehow expects to be informed about the plasma conditions at and beyond this shock. However, it has become clear meanwhile that the Voyager plasma instruments are unfortunately unable to detect the essential features of the heliosheath plasma, such as “middle-energetic” pick-up protons and low-energetic solar wind electrons. Therefore, even at the present time there is little observational information available for example about shock-processed pick-up protons and electrons, if not otherwise derivable from independent data. This means that our current understanding of the particular characteristics of this multi-fluid magnetohydrodynamics (MHD) shock, the thermodynamic essentials, and its internal entropy dissipation is still based on purely theoretical guesses and principles (see early papers by Decker et al. 2008; Izmodenov & Baranov 2006; Fahr & Chalov 2008; Wu et al. 2009; Zank et al. 2010; Chalov & Fahr 2011, and more recent papers by Fahr et al. 2015, 2017; Chalov 2019a,b; Fahr & Dutta-Roy 2019).

To nevertheless find out something more conclusive about the thermodynamic conditions of that poorly understood heliosheath plasma, we probe this plasma by studying the propagation properties of traveling shock disturbances inside it, assuming homogeneous plasma conditions. We use a kind of seismological procedure to investigate the properties of traveling shock waves. However, the assumption of homogeneous heliosheath conditions at first glance appears inappropriate in view of theoretically calculated heliosheath pressure maps published by Fahr & Siewert (2013) which show strongly pronounced 3D structures of the downstream plasma conditions. Such a simplifying assumption could nevertheless turn out to be surprisingly reasonable, since, in a broader view, at least the upwind heliosheath plasma is characterized as a system of low-Mach-number plasma flows. This fact approximately enforces plasma incompressibility and guarantees nearly isobaric conditions. This also points to nearly isothermal conditions, as indeed we show below by deriving theoretical expressions for the relevant pressure carriers in the heliosheath and by developing a method to determine the local shock velocity from in-situ heliosheath proton data.

2. Fluid pressures downstream of the multi-fluid solar wind termination shock

The problem of describing the solar wind plasma at its passage over the solar wind termination shock has been approached in the past by the application of Rankine-Hugoniot relations regulating the continuity of conserved momentum fluxes from the upstream to the downstream region. During the last decade these relations have been applied in the form of a systematically growing generality and complication, according to currently available information. This means these relations were formulated and solved in a growing refinement process, ascending from two-fluid approaches (see Chalov & Fahr 1994, 1995), via three-fluid approaches (see Chalov & Fahr 1996, 1997), finally to five-fluid approaches (Fahr et al. 2000) by systematically taking into account more and more kinetically independent, but dynamically coupled plasma fluids like solar protons, pick-up protons, electrons, anomalous cosmic rays (ACRs) and galactic cosmic rays (GCRs).

More recently, the particular role of electrons at the plasma passage over the shock has been investigated in more complexity (Chalov & Fahr 2013; Fahr & Dutta-Roy 2018, 2019; Chalov 2019a,b) and from these studies it became evident that, concerning their pressure, downstream electrons are comparable with downstream protons, even when the pressures of solar and pick-up protons are taken together.

The most interesting property in the context of the present paper is perhaps the local magnetic configuration of the shock, expressed in terms of the tilt angle α with cos(α) = cos(B, n) of the local magnetic field B with respect to the local normal n of the shock surface. Including the solar wind magnetic field B, the termination shock must be described by the more complicated magneto-hydrodynamic Rankine-Hugoniot shock conditions. As shown in Fahr (2007), Fahr & Chalov (2008), Fahr et al. (2012), two kinetic quantities ${\mathrm{\Gamma }}_{\Vert }=\frac{{\mathit{v}}_{\Vert }B}{\rho }$ (differential field parallel driftmotion) and ${\mathrm{\Gamma }}_{\perp }=\frac{{\mathit{v}}_{\perp }^2}{B}$ (magnetic particle moment) have to be conserved as particle invariants at the passage of the particles over the shock where v_∥ and v_⊥ are particle velocity components parallel and perpendicular to the field B. Under conservation of these kinetic invariants it can be shown (Fahr & Siewert 2013) that the following pressure transformations are valid at the termination shock passage (subindex 1 means upstream, subindex 2 means downstream):

$$\begin{array}{cc}& P_{2,\mathrm{p}}=\frac{s}{3}(2A(\alpha )+B(\alpha ))P_{1,\mathrm{p}},\hfill \end{array}$$(1)

$$\begin{array}{cc}& P_{2,\mathrm{pui}}=\frac{s}{3}(2A(\alpha )+B(\alpha ))P_{1,\mathrm{pui}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& P_{2,\mathrm{e}}=\frac{M}{m}\frac{s^2-1}{s^2}\frac{U_1^2}{c_{1,\mathrm{e}}^2}[A(\alpha ){sin}^2\alpha +B(\alpha ){cos}^2\alpha ]P_{1,\mathrm{p}},\hfill \end{array}$$(3)

where the suffixes p, pui, and e denote the respective quantities for protons, pick-up ions, and electrons, respectively. The tilt angle functions A(α) and B(α) are given by:

$$\begin{array}{c}\hfill A(\alpha )=\sqrt{{cos}^2\alpha +s^2{sin}^2\alpha },\end{array}$$(4)

and

$$\begin{array}{c}\hfill B(\alpha )=\frac{s^2}{A^2(\alpha )},\end{array}$$(5)

where U₁ denotes the upstream proton bulk velocity, c_1, e is the upstream electron thermal velocity, and s = n₂/n₁ denotes the shock compression ratio. M and m are the masses of protons and electrons, respectively.

To apply the above Eqs. (1)–(3) for the calculation of the plasma pressures resulting downstream of the termination shock, i.e., of protons, P_2, p, pick-up protons, P_2, pui, and electrons, P_2, e, we require knowledge of the magnetic tilt angle α (see Fig. 1) at the shock as a function of the shock position coordinates r_S and θ. For this purpose one can use the solar wind magnetic field model given in the form of the Parker spiral (see back to Parker 1958) and evaluate this field B at the shock surface which is approximated as a three-dimensional prolate ellipsoid with the sun in one of the ellipsoid focus points (see the model given by Scherer & Fahr 2009). Figures 2 and 3 show Mollweide projection plots in galactic coordinates (see e.g. McComas et al. 2009) showing the tilt angle α = α(r_S, θ) and the total plasma pressure downstream of the termination shock, that is, P₂(θ) = P_2, p(θ)+P_2, pui(θ)+P_2, e(θ). As evident from Fig. 3 the total pressure of the downstream plasma has a pronounced 3D structure, and consequently the heliosheath plasma leaving the termination shock also has a 3D temperature inhomogeneity.

Fig. 1. Schematic view of the MHD-termination shock configuration (upwind hemisphere).

Figures 2 and 3 show the tilt angle α(r_s(θ),θ) and the total pressure P₂(r_s(θ),θ) for the shock surface which in our case here approximates solar maximum conditions (i.e., high-speed solar wind all over). For solar minimum conditions with low-speed winds at small ecliptic latitudes and high-speed winds at moderate and high ecliptic latitudes the shock surface needs to be approximated by a three-axis ellipsoid, slightly complicating the above results correspondingly. Also, the variation of the neutral magnetic sheath plays an important role for the actual tilt angle α as was explicitly studied by Scherer & Fahr (2009). This further complicates the above mentioned results in view of actual plasma data. Therefore, for our purposes here, in our present endeavours we are satisfied with Fig. 3, showing one example for the course of the solar activity cycle and pointing to the fact that the termination shock surface is the origin of a heliosheath plasma with pronounced inhomogeneities of pressure and temperature in galactic longitude and latitude.

Fig. 2. Magnetic tilt angle α at the shock surface (solar maximum configuration). Data are given in a Mollweide plot in Galactic coordinates. Copyright by Fahr & Siewert (2013).

Fig. 3. Total downstream pressure in dynes cm−2, displayed in a Mollweide plot in galactic coordinates. Copyright by Fahr & Siewert (2013) (upwind nose is in the center; solar maximum configuration). The dashed red curve gives the sound velocity vHS (see left scale) at an ecliptic cut.

3. The initial fluid pressures downstream of the termination shock

Here we follow calculations carried out recently by Fahr & Dutta-Roy (2019). For the initial pressures of the two fluids, i.e., electrons and protons, immediately downstream of the termination shock, assuming kappa functions $f_{\mathrm{e}}^{\kappa }(s_0,\mathit{v})$ and $f_{\mathrm{i}}^{\kappa }(s_0,\mathit{v})$ as initial distribution functions at the initial streamline coordinate s = s₀, for example at the shock-crossing-point of Voyager-1 (see Fahr et al. 2017), we find the following expressions for the ion and electron kappa pressure:

$$\begin{array}{c}\hfill P_{\mathrm{i}}^{\kappa }(s_0)=\frac{M}{2}{n}_{\mathrm{i}0}{\mathrm{\Theta }}_{\mathrm{i}0}^2\frac{{\kappa }_{\mathrm{i}}(s_0)}{{\kappa }_{\mathrm{i}}(s_0)-3/2},\end{array}$$(6)

and for electrons,

$$\begin{array}{c}\hfill P_{\mathrm{e}}^{\kappa }(s_0)=\frac{m}{2}{n}_{\mathrm{e}0}{\mathrm{\Theta }}_{\mathrm{e}0}^2\frac{{\kappa }_{\mathrm{e}}(s_0)}{{\kappa }_{\mathrm{e}}(s_0)-3/2},\end{array}$$(7)

where the Θ_i0, κ_i0, and Θ_e0, κ_e, 0 are parameters of the initial kappa functions $f_{\mathrm{e}}^{\kappa }(s_0,\mathit{v})$ and $f_{\mathrm{i}}^{\kappa }(s_0,\mathit{v})$. With the reasonable assumption that due to charge neutrality requirements, the densities of electrons and ions downstream from the shock, after passage through the space charge layer and after flow relaxation, need to be identical, i.e., n_i0 = n_e0, and with kappa parameters calculated in Fahr & Dutta-Roy (2019) this would give the following pressure ratio:

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{\mathrm{i}/\mathrm{e}}(s_0)=\frac{P_{\mathrm{i}}^{\kappa }(s_0)}{P_{\mathrm{e}}^{\kappa }(s_0)}=\frac{M{\mathrm{\Theta }}_{\mathrm{i}0}^2\frac{{\kappa }_{\mathrm{i}}(s_0)}{{\kappa }_{\mathrm{i}}(s_0)-3/2}}{m{\mathrm{\Theta }}_{\mathrm{e}0}^2\frac{{\kappa }_{\mathrm{e}}(s_0)}{{\kappa }_{\mathrm{e}}(s_0)-3/2}}=\frac{M}{m}\frac{{(0.4U_0)}^2}{\frac{3}{m}M{U}_1^2}\frac{\frac{1.6}{1.6-1.5}}{\frac{1.517}{1.517-1.5}}·\end{array}$$(8)

With a compression ratio at the termination shock of Γ_1, 2 = n₀/n₁ = U₁/U₀ = 2.5 (see Richardson et al. 2008) one then finds the surprising result:

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{\mathrm{i}/\mathrm{e}}(s)=\frac{P_{\mathrm{i}}^{\kappa }(s_0)}{P_{\mathrm{e}}^{\kappa }(s_0)}\simeq 0.512,\end{array}$$(9)

that is, astonishingly, the electron pressure dominates the ion pressure by a factor of two – and will, as is easily proven, continue doing so over the next year of travel time down from the termination shock.

Moreover, this reveals an important and until-now unidentified fact: the effective streaming configuration, i.e., the flux tube geometry, may essentially be determined both by the electron and the ion pressure, and eventually the ion population simply behaves as a co-flowing ghost fluid embedded in an electron fluid. This is completely different from the modeling approach carried out by Heerikhuisen et al. (2019), for example, in which the electron pressure near the termination shock was adopted to only amount to 10% of the ion pressure with an unmotivated, fictionally linear approach to the ion pressure over the next 800 AU of travel path in the heliosheath.

In contrast, Chalov (2019a), based on a preceding work by Chalov & Fahr (2013), began to take an alternative approach, trying to fit the Voyager-2 electron flux data with the assumption of kappa-distributed electrons, and found that the electron pressure downstream from the termination shock is larger than the ion pressure by a factor of 6.7, while in the earlier approach Chalov & Fahr (2013) even found a factor of 10. The heliosheath plasma is therefore a matter of intense debate, and would perhaps benefit from new additional information arising from a theoretical study of kappa distributions. This is discussed in detail in a recent publication by Scherer et al. (2019), where the authors analyze the fact that velocity moments of kappa functions, constructed as integrals over the velocity space, especially the pressure moment, when using an unlimited upper integration border for the velocity, would illegally contain contributions from particles with superluminal velocities (v >  c). These latter particles do not exist in reality and therefore cannot contribute to the physical pressure. Since such moment contributions do not make sense in such a calculation, one has to try to remove their contributions from the calculated values of the pressures. Depending on the kappa parameters of the electron distribution function shown by Scherer et al. (2019), the calculated pressure values have to be reduced in order to be free from these unreasonable, artefactual contributions by superluminal particles. The question then only remains as to whether or not the complete removal of contributions by particles with v ≥ c is recommended here, or whether or not a better, more physically justified calculation is needed with some compensation for the lost population. These questions have been addressed in more detail in the most recent paper by Fahr & Dutta-Roy (2019), showing that electron pressures relative to proton pressures are more reduced by this circumstance.

4. Cooperation of the driver and the ghost fluid downstream of the termination shock

The general physics of a twin-fluid plasma has been treated by Goedbloed & Poedts (2004). However, here a very special twin-fluid situation becomes apparent that requires special attention.

As outlined in this section, it became clear that, when studied as decoupled fluids, it may turn out that the electrons in the heliosheath represent the pressure-dominating fluid, while the protons as fluid when leaving the termination shock only contribute a minor part of the total fluid pressure. This might mean that the electrons by their dominant pressure also essentially determine the plasma streaming velocity U and the flow geometry in the heliosheath, while at least at first glance, protons might appear to rather behave like a co-convected ghost fluid (tracer gas).

However, it must also be taken into account that while the electrons perhaps dominate the pressure, the protons clearly dominate the plasma mass and momentum flows. This brings up the unusual condition that wherever the ion flow is accelerated or decelerated, the electron pressure gradient might be responsible for these momentum changes and has to do work correspondingly. As thermodynamic feedback to this energy expense, naturally the electron pressure has to react accordingly.

Therefore, looking back into the recent study on electron pressures in the heliosheath by Fahr & Dutta-Roy (2019) one can identify the need for an additional term that in view of this new situation must perhaps be included in their pressure transport Eq. (3) to more correctly describe the plasma dynamics, namely a term that describes the change of the electron pressure due to energy that has to be expended by the electrons in order to cause changes in the plasma flow velocity, that is, the ion bulk velocity U(s) along the streamline.

At first glance, it may be suggested that this term in the solar rest frame consequently should be formulated in the following form:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{P}_{\mathrm{e}}}{\mathrm{d}s}=-\frac{\mathrm{d}}{\mathrm{d}s}\left[n\frac{M}{2}{U}^2\right],\end{array}$$(10)

and as such needs to be added as an additional term on the right side of Eq. (3) of Fahr & Dutta-Roy (2019). In view of the above equation and the expected decrease of the bulk velocity U along streamlines downstream from the termination shock (see dynamic models used in Fahr et al. 2016) under the present view given with a dominating electron pressure, would have the counter-intuitive effect of even increasing the dominant electron pressure to cause a decrease of the plasma bulk velocity U. The latter decrease is expected when deriving the plasma bulk velocity and streaming geometry, as is usually done with the help of a streaming potential Φ(r) (see Parker 1963).

This potential flow approximation is permitted and appropriate in cases where the plasma can be assumed to be incompressible, i.e., dn/ds = 0. However, this would then lead to a counter-intuitive situation whereby in the heliosheath regions downstream from the upwind termination shock, where bulk velocities according to standard models decrease by about 30% (e.g., see Fig. 3 in Fahr et al. 2016), the electron pressures should, in reaction to that, even continue to increase in these regions.

However, this counterintuitive situation may have been provoked by the exaggerating assumption that exclusively the electron pressure dictates the plasma flow. Perhaps in addition one should also face the problem here, that under the discussed situation of a plasma determined in its dynamics by the pressures of two independent plasma fluids, i.e., the electron fluid and the ion fluid, single-fluid solutions derived from one streaming potential Φ may no longer be applicable. What to do under these conditions?

This unexpected phenomenon suggested by the above equation may certainly appear a slightly controversial, and therefore for more clarification of this point, we look at the situation from a slightly different view, following the standard thermodynamical procedure which states that the work done by the pressure at a change of the co-moving plasma volume ΔV is reflected by an associated change of the internal energy ϵ of that volume. This requires that in the solar rest frame (SRF) the following equation has to be valid

$$\begin{array}{c}\hfill -(P_{\mathrm{e}}+P_{\mathrm{i}})\frac{\mathrm{d}\mathrm{\Delta }V}{\mathrm{d}s}=\frac{\mathrm{d}}{\mathrm{d}s}[({ϵ}_{\mathrm{i}}+{ϵ}_{\mathrm{e}})\mathrm{\Delta }V],\end{array}$$(11)

where ΔV, as explained in Fahr & Dutta-Roy (2019), denotes the co-moving plasma volume on the streamline, i.e., a fluid volume that locally moves with the plasma bulk velocity U. Hereby the indices “i, e” indicate ion- or electron-related quantities as pressures and internal energies, respectively.

Now taking into account the fact that in the SRF the ion energy density is ϵ_i = nMU²/2 + (3/2π)P_i while the electron energy density only is given by ϵ_e = (3/2π)P_e (i.e., strongly subsonic electron flow). Furthermore assuming, in order to begin from some basis, that the electron pressure dominates over the ion pressure, i.e., P_e ≫ P_i, will then bring us to the following net equation:

$$\begin{array}{c}\hfill -P_{\mathrm{e}}\frac{\mathrm{d}\mathrm{\Delta }V}{\mathrm{d}s}=\frac{\mathrm{d}}{\mathrm{d}s}[(nM{U}^2/2+\frac{3}{2\pi }{P}_{\mathrm{e}})·\mathrm{\Delta }V].\end{array}$$(12)

When additionally recognizing here that for an incompressible flow, as given in the case of a strongly subsonic flow, the comoving plasma volume is given by the following relation: ΔV = ΔV₀ ⋅ (U₀/U) (see Fahr & Dutta-Roy 2019). Then the above equation simplifies into:

$$\begin{array}{cc}& -P_{\mathrm{e}}\frac{\mathrm{d}\frac{1}{U}}{\mathrm{d}s}=\frac{\mathrm{d}}{\mathrm{d}s}\left[\frac{1}{2}nMU\right]+\frac{3}{2\pi }\frac{\mathrm{d}}{\mathrm{d}s}\left(P_{\mathrm{e}}\frac{1}{U}\right),\hfill \end{array}$$(13)

$$\begin{array}{cc}& \frac{2\pi +3}{2\pi }{P}_{\mathrm{e}}\frac{\mathrm{d}U}{\mathrm{d}s}-\frac{3}{2\pi }U\frac{\mathrm{d}{P}_{\mathrm{e}}}{\mathrm{d}s}=\frac{1}{2}nM{U}^2\frac{\mathrm{d}U}{\mathrm{d}s}·\hfill \end{array}$$(14)

This equation requires a solution of the form P = P(U), and with dU/ds ≠ 0 leads to

$$\begin{array}{c}\hfill \frac{2\pi +3}{2\pi }{P}_{\mathrm{e}}-\frac{3}{2\pi }U\frac{\mathrm{d}{P}_{\mathrm{e}}}{\mathrm{d}U}=\frac{1}{2}nM{U}^2.\end{array}$$(15)

Obviously the solution of the upper differential equation requires P_e to be a function of U, tentatively by the following representation

$$\begin{array}{c}\hfill P=P_0·{(U/U_0)}^2.\end{array}$$(16)

When inserting this into the upper differential equation, one then finds the requirement

$$\begin{array}{c}\hfill \frac{2\pi -3}{2\pi }{P}_o\frac{U^2}{U_0^2}=\frac{1}{2}nM{U}^2,\end{array}$$(17)

or yielding the initial electron pressure in the form:

$$\begin{array}{c}\hfill P_o=\frac{2\pi }{2\pi -3}\frac{1}{2}nM{U}_0^2=0.96·nM{U}_0^2.\end{array}$$(18)

This solution now expresses the fact that the electron pressure falls off with the plasma bulk velocity like U², and that at a decrease of the bulk velocity with increasing line element s in the heliosheath this pressure should also decrease accordingly. In addition it would mean that with the electrons moving downstream from the termination shock along streamlines with decreasing bulk velocity (see e.g., Fahr et al. 2016), the electron pressure should also decrease accordingly, so as to consequently bring the electron and ion pressures closer together at the motion downstream from the termination shock and thus, in principle, to re-establish in this way a kind of thermodynamic quasi-equilibrium plasma condition further down from the shock.

With that result, coming now back to the fact that the electron pressure performs thermodynamical work, when driving down the streamline the electron plasma, one must conclude that without any interaction of ions and electrons, this energy, which has to be thermodynamically expended, has to be taken from the internal thermal energy ϵ_e of the electrons themselves. This leads to the following term describing the decrease of the electron thermal energy

$$\begin{array}{c}\hfill -P_{\mathrm{e}}\frac{\mathrm{d}\mathrm{\Delta }V}{\mathrm{d}s}=\frac{\mathrm{d}}{\mathrm{d}s}\left[{ϵ}_{\mathrm{e}}\mathrm{\Delta }V\right]=\frac{\mathrm{d}}{\mathrm{d}s}\left[\frac{3}{2\pi }{P}_{\mathrm{e}}\mathrm{\Delta }V\right].\end{array}$$(19)

Together with the relation for the comoving fluid volume in incompressible flows, ΔV = ΔV₀(U₀/U), this leads to the following expression:

$$\begin{array}{c}\hfill -P_{\mathrm{e}}\frac{\mathrm{d}}{\mathrm{d}s}\frac{1}{U}=\frac{\mathrm{d}}{\mathrm{d}s}\left[\frac{3}{2\pi }{P}_{\mathrm{e}}\frac{1}{U}\right]=\frac{3}{2\pi }[P_{\mathrm{e}}\frac{\mathrm{d}}{\mathrm{d}s}\frac{1}{U}+\frac{1}{U}\frac{\mathrm{d}{P}_{\mathrm{e}}}{\mathrm{d}s}],\end{array}$$(20)

which can be simplified to

$$\begin{array}{c}\hfill -\frac{2\pi }{3}(1+\frac{3}{2\pi })U\frac{\mathrm{d}}{\mathrm{d}s}\frac{1}{U}=\frac{1}{P_{\mathrm{e}}}\frac{\mathrm{d}{P}_{\mathrm{e}}}{\mathrm{d}s},\end{array}$$(21)

and consequently yielding a pressure change due to volume work as given by:

$$\begin{array}{c}\hfill \frac{\mathrm{d}ln{P}_{\mathrm{e}}}{\mathrm{d}s}=\frac{2\pi +3}{3}\frac{\mathrm{d}lnU}{\mathrm{d}s}·\end{array}$$(22)

When this is combined with the other terms in the pressure transport equation (see Eq. (14) in Fahr & Dutta-Roy 2019), it then leads to the following completed form:

$$\begin{array}{c}\hfill \frac{\mathrm{d}ln{P}_{\mathrm{e}}}{\mathrm{d}s}=\frac{4}{3}\frac{\mathrm{d}lnB}{\mathrm{d}s}+\frac{10D_0}{U}+\frac{2\pi +3}{3}\frac{\mathrm{d}lnU}{\mathrm{d}s}·\end{array}$$(23)

This differential equation can be integrated and leads to the following completed solution for the electron pressure:

$$\begin{array}{c}\hfill P_{\mathrm{e}}(s)=P_{\mathrm{e}0}{\left(\frac{B}{B_0}\right)}^{4/3}{\left(\frac{U}{U_0}\right)}^{\frac{2\pi +3}{3}}exp\left[10D_0{\displaystyle {\int }_{s_0}^s}\frac{\mathrm{d}s}{U}\right]·\end{array}$$(24)

This solution shows that the electron pressure decreases with the plasma bulk velocity proportional to $U{(s)}^{\frac{2\pi +3}{3}}=U{(s)}^{3.09}$, but also that, in addition, frozen-in magnetic fields enforcing the conservation of magnetic particle moments and wave-electron diffusion may independently and additionally modify the electron pressure.

5. Plasma sound velocity and incompressibility

To bring the description of electron and ion flow dynamics into a more consistent form, one might check on two things which up to now were tacitly assumed in the derivation: The first is the assumption of plasma incompressibility, i.e., n_e(s ≥ s₀) = n_i(s ≥ s₀) = const. The second is that a joint streaming potential Φ(r) = Φ_e, i(r) can be found for this incompressible plasma from which the joint streamlines of this twin-fluid plasma can be derived with the help of a relation like U(r) = grad[Φ_e, i(r)] ≃ grad[Φ_e(r)+Φ_i(r)].

The problem of incompressibility is closely related to the effective sound velocity $c_{\mathrm{s}}^{\mathrm{eff}}$ of this twin-fluid plasma which is defined by

$$\begin{array}{c}\hfill {\left(c_{\mathrm{s}}^{\mathrm{eff}}\right)}^2=\frac{\mathrm{d}P}{\mathrm{d}\rho }=\frac{1}{M}\frac{\mathrm{d}}{\mathrm{d}n}[P_{\mathrm{e}}+P_{\mathrm{i}}].\end{array}$$(25)

Assuming, despite the prevailing nonequilibrium conditions, the validity of an effective polytropic relation in the form c_γ = P/ρ^γeff = (P_e + P_i)/ρ^γeff, with γ^eff ≤ 5/3 denoting the effective polytrope index, one then finds:

$$\begin{array}{c}\hfill c_{\gamma }{\rho }^{\gamma }=P_{\mathrm{e}}+P_{\mathrm{i}},\end{array}$$(26)

meaning

$$\begin{array}{c}\hfill \frac{\mathrm{d}(P_{\mathrm{e}}+P_{\mathrm{i}})}{\mathrm{d}\rho }=\gamma c_{\gamma }{\rho }^{\gamma -1}=\gamma c_{\gamma }\frac{{\rho }^{\gamma }}{\rho },\end{array}$$(27)

and yielding:

$$\begin{array}{c}\hfill \frac{\mathrm{d}(P_{\mathrm{e}}+P_{\mathrm{i}})}{\mathrm{d}\rho }=\gamma \frac{P_{\mathrm{e}}+P_{\mathrm{i}}}{\rho },\end{array}$$(28)

and the effective shock propagation velocity given by:

$$\begin{array}{c}\hfill c_{\mathrm{s}}^{\mathrm{eff}}=\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho }[C_{\gamma }{\rho }^{{\gamma }^{\mathrm{eff}}}]}=\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho }}\sqrt{[P_{\mathrm{e}}+P_{\mathrm{i}}]},\end{array}$$(29)

already indicating that the prevailing effective sound velocity, under approximately equal pressure conditions with P_e ≃ P_i, will be greater by a factor of $\sqrt{2}$ compared to pure ion sound speed. Furthermore, given the values for the electron and ion pressures downstream of the termination shock as mentioned above, this speed turns out to be higher than the bulk velocity U by at least a factor of ten, yielding an effective Mach number of this twin-fluid plasma flow by $M^{\mathrm{eff}}=U/c_{\mathrm{s}}^{\mathrm{eff}}$ ≤0.1.

Concerning the incompressibility, now the following argument can be made: Assuming irrotational heliosheath flows, then along the streamlines the Bernoulli constant Be = (1/2)ρU² + P holds. For a streamline extending from the termination shock to the stagnation point, connecting regions of maximum pressure differences, this then easily tells us that the pressure at the shock P_shock = P(s = s₀) is nearly equal to the pressure at the stagnation point P_stag with (P_stag − P_shock)/P_stag ≤ 10⁻², and hence with the polytropic law C = P/ρ^γeff one can also derive:

$$\begin{array}{c}\hfill n_{\mathrm{stag}}/n_{\mathrm{shock}}={(P_{\mathrm{shock}}/P_{\mathrm{stag}})}^{{\gamma }^{\mathrm{eff}}}\simeq 1,\end{array}$$(30)

that is, incompressibility can be assumed under the given conditions and the theoretical concept of potential flows can in principle also be used for the description of this twin-fluid plasma; the question rather is what the twin-fluid potential may look like.

First we want to draw some conclusions from the above-derived result for the effective sound velocity and start by considering the conditions in the upwind hemisphere of the heliosheath. There we find strongly wound Archimedian B-fields, and therefore the Alfvénic shock modes are propagating mostly in the Archimedian field direction and essentially do not count for the shock wave propagation in radial directions. In addition, the magnitude of the Alfvén velocity is estimated to amount to 39 km s⁻¹ (see Rankin et al. 2019). Hence one can safely neglect these Alfvénic contributions to the radial shock propagation speed.

This means that we remain with an effective sound speed of:

$$\begin{array}{c}\hfill c_{\mathrm{s}}^{\mathrm{eff}}=\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho }}\sqrt{[P_{\mathrm{e}}+P_{\mathrm{i}}]}\end{array}$$(31)

downstream of the termination shock. As the geometry of the solar maximum termination shock is approximated by a quasi-spherical ellipsoid, one can assume that, in the upwind hemisphere, the shock radius r_s is approximately constant with r_s(θ) ≃ r_s ≃ 90 AU, which leads to the result that the downstream density ρ(θ) at the shock point r_s(θ) is given by

$$\begin{array}{c}\hfill \rho (\theta )={\rho }_{\mathrm{E}}·{(r_{\mathrm{E}}/r_{\mathrm{s}})}^2·s(\theta ).\end{array}$$(32)

Taking the result displayed in Fig. 8a of Fahr & Siewert (2013) one can see that the shock compression ratios s(θ) in the upwind hemisphere are given by values 3.2 ≤ s(θ) ≤ 3.4 and thus allow us to approximate the effective sound velocity by:

$$\begin{array}{cc}\hfill c_{\mathrm{s}}^{\mathrm{eff}}(\theta )=& \sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho (\theta )}}\sqrt{[P_{\mathrm{e}}(\theta )+P_{\mathrm{i}}(\theta )]}\hfill \\ \hfill & =\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{{\rho }_{\mathrm{E}}·{(r_{\mathrm{E}}/r_{\mathrm{s}})}^2}}\sqrt{\frac{1}{s(\theta )}}\sqrt{[P_{\mathrm{e}}(\theta )+P_{\mathrm{i}}(\theta )]}\hfill \\ \hfill & \simeq \sqrt{\frac{{\gamma }^{\mathrm{eff}}}{3.3·{\rho }_{\mathrm{E}}{(r_{\mathrm{E}}/r_{\mathrm{s}})}^2}}\sqrt{P(\theta )},\hfill \end{array}$$(33)

where P(θ) = P_e(θ)+P_i(θ) denotes the total pressure downstream of the termination shock displayed in Fig. 3.

From Fig. 3, it can now easily be deduced that over the upwind hemisphere of the heliosheath (i.e., within angles θ = 0° up to θ = 90°) the effective sound velocity should vary by about a factor of $\sqrt{P({90}^{°})/P(0^{°})}=\sqrt{8.4/4.4}=1.38$.

This expresses the fact that the shock wave near the upwind axis at θ = 0° is propagating with a shock speed of about a factor of 1.4 higher than the shock wave leaving the termination shock at angles θ ≃ 90° which all together leads us to predict that the global shock front leaving the termination shock and propagating through the heliosheath, instead of keeping its radial shape, would become strongly deformed, creating a strong bulge at the upwind axis – the stronger the longer it propagates through the heliosheath – provided the pressure asymmetry conditions are enduring. Adopting the observationally derived shock speed of about 300 km s⁻¹ obtained by Rankin et al. (2019), with a shock velocity asymmetry factor of 1.4, we arrive at a shock front distance difference of Δr(θ ≃ 0°, Θ ≃ 90° ) ≃ (121−90) AU/300 km s⁻¹ ⋅ [(1.4−1)300 km s⁻¹] = (0.4 ⋅ 31) AU = 12.4 AU. This characterizes the relative nonsphericity of the great merging interaction region (GMIR) shock front by a value of δr = Δr/r = 12.4/121 = 0.103, which does not look too impressive. The nonsphericity of the GMIR shock may be put further into doubt by the fact that such a shock already in the inner heliosphere is not spherical (see e.g., simulations carried out by Opher et al. 2015).

Another reason why pressure inhomogeneity does not strongly complicate further considerations is that streamlines originating from the regions of the termination shock with the highest pressures will diverge more strongly compared to those originating from moderate latitudes with lower pressures. This stronger divergence of the upwind streamlines, due to a stronger cooling by expansion work along the flux tubes, leads to faster cooling of the heliosheath plasma at radial propagation (see an illustration of the principle behind this phenomenon in Fig. 4). The effect of this can easily be predicted in qualitative terms: At its radial propagation the hotter upwind heliosheath plasma will cool faster compared to the cooler plasma propagating from regions at higher inclinations Θ. The effect in qualitative terms then is that at radial propagation the originally different pressures of neighboring plasma flows will assimilate each other creating approximately isobaric and isothermal heliosheath plasma conditions at moderate distances from the termination shock.

Fig. 4. Example for streamlines given in cylindrical coordinates z and 𝜚, given in units of the shock radius, and generated from the scalar flow potential $\mathrm{\Phi }(z,\varrho )=U_0(z+\frac{2}{\varrho })$. The inner circle approximates the termination shock. Copyright by Fahr et al. (2016).

Perhaps one additional point to be inspected here is the fact that the GMIR shock cannot be assumed to propagate essentially with the sound speed of the protons. In the case of the heliosheath plasma where protons and electrons contribute approximately equal partial pressures, the dominant factor is the effective sound speed $c_{\mathrm{s}}^{\mathrm{eff}}$ given by

$$\begin{array}{c}\hfill c_{\mathrm{s}}^{\mathrm{eff}}=\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho }}\sqrt{[P_{\mathrm{e}}+P_{\mathrm{i}}]}.\end{array}$$(34)

Based on the results of the above sections, there are good reasons to assume that electron and proton pressures are of the same magnitude, i.e., P_e ≃ P_i = P/2, thus leading to

$$\begin{array}{c}\hfill c_{\mathrm{s}}^{\mathrm{eff}}(P)=\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho }}\sqrt{2P_{\mathrm{i}}},\end{array}$$(35)

and the proton pressure is derived from that as

$$\begin{array}{c}\hfill P_{\mathrm{i}}=\frac{1}{2}(\rho /{\gamma }^{\mathrm{eff}})c_{\mathrm{s}}^{{\mathrm{eff}}^2},\end{array}$$(36)

yielding the total pressure P = P_i + P_e = 2P_i as

$$\begin{array}{c}\hfill P=2P_{\mathrm{i}}=(\rho /{\gamma }^{\mathrm{eff}})c_{\mathrm{s}}^{{\mathrm{eff}}^2}.\end{array}$$(37)

6. How to deduce the shock velocity from available data?

We now address the question of how the magnitude of the above shock velocity could be extracted from in-situ data that are available for example from the Voyager-1/-2 measurements, even though, as mentioned earlier, neither the pick-up ions contributing the majority of the ion pressure nor the solar wind electrons can be distinctly measured. Nevertheless, the following offers a good opportunity to extract the relevant information on the total ion pressure P_i from Voyager-2 data, even though most of the ions contributing to this pressure P_i cannot be registered directly by Voyager-2.

This is connected with the fact that a traveling shock wave passing over Voyager-2 will give a kick to the whole local plasma and consequently in the shock ramp, for a certain time period τ_X, the local density n_p(r, t) of the local low-energy protons will also increase, showing a typical density enhancement signature. However, the local proton density increase is not directly related to the total pressure P_i, but can serve as a helpful signature or as a useful timer of the shock wave passage, by indicating its typical propagation velocity by means of the associated signature of the associated proton density enhancement as function of time, provided that this temporal enhancement signature can be clearly analyzed in the data.

In a paper by Fahr et al. (2014), we looked into the physics of such traveling shock waves and found with the help of Burger’s nonlinear shock velocity relation (see e.g., Treumann & Baumjohann 1997) that a shock velocity profile of the following form must be expected for the wave:

$$\begin{array}{c}\hfill \delta U=\mathrm{\Delta }U[1-tanh\left(\frac{\mathrm{\Delta }U}{2\xi }y\right)],\end{array}$$(38)

where y = x − ΔU ⋅ t, with ΔU ≃ c_s denoting the total bulk velocity step, x being the local space coordinate in the direction of the shock velocity ΔU, and t the local time coordinate. Here, ξ is a nonlinear wave dissipation or spatial diffusion coefficient with dimension [cm² s⁻¹] whose value can be estimated by the typical timescale of density wave steepening or dissolution at the traveling wave front given in the form of the following solution:

$$\begin{array}{c}\hfill \delta U(x,t)=c_{\mathrm{s}}\frac{2D}{\sqrt{4\pi \xi t}}exp[-\frac{x^2}{4\xi t}]·\end{array}$$(39)

The resulting shock profile stabilized by nonlinear wave dissipation then evaluates to a typical length dimension of X = 4ξ/c_s (see also illustration in Fig. 5). The typical time period τ_X of the passage of this shock-induced density structure is therefore given by:

$$\begin{array}{c}\hfill {\tau }_X=\frac{X}{c_{\mathrm{s}}}=\frac{4\xi }{c_{\mathrm{s}}^2}·\end{array}$$(40)

Fig. 5. Schematic illustration of the traveling shock front showing the length X of the proton density enhancement.

Assuming now that this characteristic time period τ_X is extractable from Voyager-2 low-energy proton data (see e.g., Fig. 2 of Richardson & Decker 2014, or Richardson et al. 2009, 2008; Li et al. 2008) one can directly derive from it the total plasma pressure P = P_i + P_e at the place of Voyager-2 by use of the simple relation:

$$\begin{array}{c}\hfill \frac{{\gamma }^{\mathrm{eff}}}{\rho }(P_{\mathrm{i}}+P_{\mathrm{e}})=c_{\mathrm{s}}^2=\frac{4\xi }{{\tau }_X}·\end{array}$$(41)

The relevant spatial diffusion coefficient ξ based on space diffusion of protons in the inner heliosheath may be given with a value of ξ = 10¹⁹ cm² s⁻¹ (see e.g., Chalov et al. 2003). For this method, the only data that can be used are those with sufficiently small running averaging periods of τ_a ≪ τ_X ≤ 10⁵ s.

To also arrive at some essential conclusions concerning the state of the heliosheath plasma pressure, Rankin et al. (2019) followed a somewhat similar procedure to that described above. These latter authors used a special mono-causal solar-physical shock-wave event that touched the heliosheath plasma during the years 2012−2013 and, as assumed by the authors, passed over both the positions of the Voyager space probes with a time delay of 130 ± 5 days. This causal event was associated with a traveling shock wave (GMIR) that entered the heliosheath and finally, after passing over the heliopause, even propagated into the ambient very local interstellar medium (VLISM) plasma. The passage of the shock wave over the two Voyagers in both cases was followed by a remarkable depression of the local GCR particle fluxes continuously monitored by the onboard Voyager detectors. By the time delay between the two GCR depression events, seen at the two Voyagers, the authors, assuming homogeneous density and temperature conditions in the heliosheath, determined the average propagation velocity of this traveling shock between the two Voyagers to be about 300 km s⁻¹ and calculated the associated average proton pressure from it. Their result regarding the associated mean heliosheath proton pressure lead to a pioneering first estimate of the heliosheath plasma properties.

7. Conclusions

In the present paper, we show that the total plasma pressure at the solar wind termination shock is strongly inhomogeneous, varying with the inclination angle Θ with respect to the upwind axis (see Fig. 3). We also show that the heliosheath plasma is a twin-fluid plasma with electrons and protons contributing about equal magnitudes to the total plasma pressure. Therefore, the typical shock propagation speed is not only inhomogeneous with respect to the angle Θ, but is also characterized by the so called “effective” sound velocity given by

$$\begin{array}{c}\hfill c_{\mathrm{s}}^{\mathrm{eff}}=\sqrt{\frac{{\gamma }^{\mathrm{eff}}}{\rho }}\sqrt{[P_{\mathrm{e}}+P_{\mathrm{i}}]}.\end{array}$$(42)

Therefore, a look at our Fig. 3, which displays the total pressure P = P_e + P_i, reveals that the traveling shock wave will propagate with a speed of $c_{\mathrm{s}}^{\mathrm{eff}}=c_{\mathrm{s}}^{\mathrm{eff}}(\mathrm{\Theta })$ with maximum values at Θ = 0° higher by a factor of 1.4 compared to propagation speeds at higher inclinations of Θ ≃ 90°. This induces a nonspherical bulge-like deformation of the shock propagating downstream into the heliosheath as we discuss in the above sections. However, this deformation will not persist or even grow during the shock propagation over the whole heliosheath because the typical plasma flow geometry shown in Fig. 4 works against the pressure and temperature inhomogeneity. This is because near the upwind axis, the flow lines have the strongest divergence; the connected plasma has to do thermodynamical work and consequently cools down much more efficiently than in regions with streamlines of less divergence. It is therefore expected that, after leaving the inhomogeneous pressure conditions at the termination shock, the plasma downstream of the shock will systematically smooth the temperature and pressure inhomogeneities. Therefore, at larger heliosheath distances, the assumption of pressure homogeneity might apply relatively well to the actual heliosheath plasma conditions.

We also derive a relation that allows us to deduce the effective sound velocity, $c_{\mathrm{s}}^{\mathrm{eff}}$, i.e., the total plasma pressure, from a special density signature in the low-energy proton data measured at Voyager-2 when the traveling shocks pass over. As we show, the low-energy protons react in the form of a temporal density enhancement when a traveling shock is passing over, and, as was shown in the previous section, using the characteristic time period ${\tau }_X=4\xi /c_{\mathrm{s}}^2$ of this density enhancement it is possible to locally measure the effective sound velocity $c_{\mathrm{s}}^{\mathrm{eff}}$.

In order to also obtain some information on the LISM counterpart pressure along a similar procedure, Rankin et al. (2019) used the short part of the Voyager-1 trajectory outside of the heliopause (i.e., 0.9 AU). This is very short compared to the distance of Voyager-2 to the expected heliopause at 121 AU radial distance and thus strongly reduces the given interstellar sensitivity to the total propagation time. This means the influence of the heliosheath plasma on the GMIR shock propagation time is vastly greater than that of the short LISM shock passage beyond the heliopause. This situation will improve in the future, as Voyager-1 moves further out into the outer heliosphere beyond the heliopause.

References

All Figures

	Fig. 1. Schematic view of the MHD-termination shock configuration (upwind hemisphere).
In the text

	Fig. 2. Magnetic tilt angle α at the shock surface (solar maximum configuration). Data are given in a Mollweide plot in Galactic coordinates. Copyright by Fahr & Siewert (2013).
In the text

	Fig. 3. Total downstream pressure in dynes cm−2, displayed in a Mollweide plot in galactic coordinates. Copyright by Fahr & Siewert (2013) (upwind nose is in the center; solar maximum configuration). The dashed red curve gives the sound velocity vHS (see left scale) at an ecliptic cut.
In the text

	Fig. 4. Example for streamlines given in cylindrical coordinates z and 𝜚, given in units of the shock radius, and generated from the scalar flow potential $\mathrm{\Phi }(z,\varrho )=U_0(z+\frac{2}{\varrho })$. The inner circle approximates the termination shock. Copyright by Fahr et al. (2016).
In the text

	Fig. 5. Schematic illustration of the traveling shock front showing the length X of the proton density enhancement.
In the text