© ESO 2021

1. Introduction

A fundamental feature of the solar atmosphere (and astrophysical systems in general) is shock waves. Shocks can form in a number of ways, for example the natural steepening of upwardly propagating waves (for example umbral flashes) or from reconnection events. These events all occur in a magnetised plasma and as such the solar atmosphere can support several different types of shock. For magnetohydrodynamic (MHD) systems, shocks can be broadly defined as slow, fast, or intermediate depending on the velocity on either side of the shock (Tidman & Krall 1971). Here we focus on the slow-mode shock, which is a fundamental feature of magnetic reconnection (Petschek 1964). The results presented here are also applicable to more general shocks in partially ionised plasmas across astrophysical systems.

The MHD model works well for studying shocks in fully ionised media. However, the lower solar atmosphere consists of both ionised and neutral particles, and hence two-fluid effects become important in this region (see the review of Khomenko 2017). Events such as umbral flashes or Ellerman bombs occur in the lower solar atmosphere and are intrinsically linked to shocks. Partially ionised shocks are important for other astrophysical systems, such as molecular clouds (Chernoff 1987; Lehmann & Wardle 2016; Godard et al. 2019).

For partially ionised two-fluid (ion+electron, neutral) systems, shocks become more complicated. There is a larger set of characteristic wave speeds, including the neutral sound speed, the slow, Alfvén, and fast plasma speeds, and the slow, Alfvén, and fast bulk (plasma plus neutral) speeds (Snow & Hillier 2020). An effective set of wave speeds also exists that involves using the collisionallity of the system to determine the slow, Alfvén, and fast speeds (Soler et al. 2013). In a two-fluid shock, the species decouple and recouple, leading to a finite width and large drift velocities in the shock (Hillier et al. 2016). Within this finite shock width, additional shock transitions can form (Snow & Hillier 2019).

The rapid changes in local properties in a shock alter the ionisation and recombination rates of the system (for example Carlsson & Stein 2002). Due to the temperature changes in and across a shock, the ionisation and recombination rates change, which has the effect of changing the composition of the medium, increasing or decreasing the neutral fraction. The most common way to model the rates of collisional ionisation and recombination is to use either empirical fits to lab data (Smirnov 2003; Voronov 1997) or a multi-level atom model (Jefferies 1968). The empirical formulas assume a statistical equilibrium, where the hydrogen species ionises as an ensemble single fluid, and do not account for different levels of excitation within the species itself. A more complete multi-level model includes non-local thermal equilibrium where the ionisation and recombination rates of individual hydrogen levels are modelled separately (for example Carlsson & Stein 2002).

During collisional ionisation, kinetic energy of a free electron is used to release an electron that is trapped in the electric field of a nucleus. The macroscopic effect of the work done to ionise an atom is the reduction in internal energy (and with it the temperature) of the plasma fluid. However, as photons can be released as part of the process of collisional recombination or de-excitation, the recombination-ionisation cycle of an atom does not necessarily conserve the energy of the fluid, potentially resulting in cooling.

In this paper, we focus on the solar atmosphere however, partially ionised shocks are also a common feature of the interstellar medium. When radiative losses are considered, the energy loss within the finite-width of the shock acts to limit the maximum temperature increase (Draine & McKee 1993). Radiative losses then cool the post-shock region as further energy is lost, despite an increase in both density and pressure across the shock.

In this paper, we study the role of collisional ionisation and recombination in a reconnection-driven slow mode shock for a two-fluid partially ionised hydrogen plasma. Our previous studies neglected ionisation and recombination and focused on the fundamental two-fluid effects (Hillier et al. 2016; Snow & Hillier 2019). Here we expand the model further towards reality, accounting for collisional ionisation, recombination and ionisation potential. We present a semi-analytical method for determining the possible equilibrium shock solutions for a radiating fluid. We then perform and analyse 1D numerical simulations. The resultant shocks have different behaviour to the MHD solutions due to the non-conservative energy equation. A key feature is that the ionisation potential acts as an energy loss term and cools the system, despite an increase in both pressure and density across the shock interface.

The results of this paper have application to a wide range of partially ionised mediums, including the solar chromosphere and molecular clouds, where the multi-fluid effects lead to substantial deviation from single-fluid MHD results. This work focuses on reconnection-driven switch-off slow-mode shocks however the results can easily be extended to other types of shock, including steepening events.

2. Methodology

In this paper, we study ionisation, recombination, ionsiation potential, and arbitrary heating in the framework of two-fluid (neutral, and ion+electron) model for hydrogen using a semi-analytical approach and numerical simulations. Specifically, in this study, we use the following set of equations:

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}})={\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}},\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial }{\partial t}({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}})+\mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}+P_{\mathrm{n}}\mathbf{I})=-{\alpha }_c{\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}({\mathit{v}}_{\mathrm{n}}-{\mathit{v}}_{\mathrm{p}})\hfill \\ \hfill & +{\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial e_{\mathrm{n}}}{\partial t}+\mathrm{\nabla }·\left[{\mathit{v}}_{\mathrm{n}}(e_{\mathrm{n}}+P_{\mathrm{n}})\right]=-{\alpha }_c{\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}[\frac{1}{2}({\mathit{v}}_{\mathrm{n}}^2-{\mathit{v}}_{\mathrm{p}}^2)+\frac{3}{2}(\frac{P_{\mathrm{n}}}{{\rho }_{\mathrm{n}}}-\frac{1}{2}\frac{P_{\mathrm{p}}}{{\rho }_{\mathrm{p}}})]\hfill \\ \hfill & +\frac{1}{2}({\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}^2-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}^2)\hfill \\ \hfill & +\frac{1}{(\gamma -1)}(\frac{1}{2}{\mathrm{\Gamma }}_{\mathrm{rec}}{P}_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{P}_{\mathrm{n}}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& e_{\mathrm{n}}=\frac{P_{\mathrm{n}}}{\gamma -1}+\frac{1}{2}{\rho }_{\mathrm{n}}{v}_{\mathrm{n}}^2,\hfill \end{array}$$(4)

$$\begin{array}{cc}& \frac{\partial {\rho }_{\mathrm{p}}}{\partial t}+\mathrm{\nabla }·({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}})=-{\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}+{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}\hfill \end{array}$$(5)

$$\begin{array}{cc}& \frac{\partial }{\partial t}({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}})+\mathrm{\nabla }·({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}+P_{\mathrm{p}}\mathbf{I}-\mathit{BB}+\frac{{\mathit{B}}^2}{2}\mathbf{I})={\alpha }_c{\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}({\mathit{v}}_{\mathrm{n}}-{\mathit{v}}_{\mathrm{p}})\hfill \\ \hfill & -{\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}+{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}},\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{\partial }{\partial t}(e_{\mathrm{p}}+\frac{{\mathit{B}}^2}{2})+\mathrm{\nabla }·[{\mathit{v}}_{\mathrm{p}}(e_{\mathrm{p}}+P_{\mathrm{p}})-({\mathit{v}}_{\mathrm{p}}\times \mathit{B})\times \mathit{B}]\hfill \\ \hfill & ={\alpha }_c{\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}[\frac{1}{2}({\mathit{v}}_{\mathrm{n}}^2-{\mathit{v}}_{\mathrm{p}}^2)+\frac{3}{2}(\frac{P_{\mathrm{n}}}{{\rho }_{\mathrm{n}}}-\frac{1}{2}\frac{P_{\mathrm{p}}}{{\rho }_{\mathrm{p}}})]\hfill \\ \hfill & -\frac{1}{2}({\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}^2-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}^2)-{\varphi }_I+A_{\mathrm{heat}}\hfill \\ \hfill & -\frac{1}{(\gamma -1)}(\frac{1}{2}{\mathrm{\Gamma }}_{\mathrm{rec}}{P}_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{P}_{\mathrm{n}}),\hfill \end{array}$$(7)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}-\mathrm{\nabla }\times ({\mathit{v}}_{\mathrm{p}}\times \mathit{B})=0,\hfill \end{array}$$(8)

$$\begin{array}{cc}& e_{\mathrm{p}}=\frac{P_{\mathrm{p}}}{\gamma -1}+\frac{1}{2}{\rho }_{\mathrm{p}}{v}_{\mathrm{p}}^2,\hfill \end{array}$$(9)

$$\begin{array}{cc}& \mathrm{\nabla }·\mathit{B}=0,\hfill \end{array}$$(10)

for a charge neutral plasma (subscript p) and neutral (subscript n) species. The fluid properties are given by density ρ, pressure P, velocity v, magnetic field B and energy e. Both species follow ideal gas laws for temperature T, namely T_n = γP_n/ρ_n and $T_{\mathrm{p}}=\frac{1}{2}\gamma P_{\mathrm{p}}/{\rho }_{\mathrm{p}}$.

The species are thermally coupled through the collisional coefficient α_c which is calculated as:

$$\begin{array}{c}\hfill {\alpha }_c={\alpha }_0\sqrt{\frac{T_{\mathrm{p}}+T_{\mathrm{n}}}{2}}\sqrt{\frac{1}{T_{\mathrm{init}}}}.\end{array}$$(11)

The factor of $\sqrt{\frac{1}{T_{\mathrm{init}}}}$ is to normalise the collisional coefficient using the initial temperature T_init such that α_c(t = 0) = α₀.

2.1. Collisional ionisation and recombination

The collisional (three-body) ionisation (Γ_ion) and recombination (Γ_rec) rates for a hydrogen atom are given by the empirical forms from Voronov (1997) and Smirnov (2003). In normalised form, these equations are:

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{rec}}=\frac{{\rho }_{\mathrm{p}}}{\sqrt{T_{\mathrm{p}}}}\frac{\sqrt{T_f}}{{\xi }_{\mathrm{p}0}}{\tau }_{\mathrm{IR}}=F(T){\rho }_{\mathrm{p}},\hfill \end{array}$$(12)

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{ion}}={\rho }_{\mathrm{p}}\frac{{\mathrm{e}}^{-\chi }{\chi }^{0.39}}{0.232+\chi }\frac{\widehat{R}}{{\xi }_{\mathrm{p}0}}{\tau }_{\mathrm{IR}}=G(T){\rho }_{\mathrm{p}},\hfill \end{array}$$(13)

$$\begin{array}{cc}& \chi =13.6\frac{T_f}{T_{e0}{T}_{\mathrm{p}}},\hfill \end{array}$$(14)

$$\begin{array}{cc}& \widehat{R}=\frac{2.91\times {10}^{-14}}{2.6\times {10}^{-19}}\sqrt{T_{e0}},\hfill \end{array}$$(15)

where T_f is a normalisation factor to ensure the simulation temperature is the reference temperature T₀ in Kelvin, and T_e0 is the reference temperature converted to electron volts. A free parameter exists (τ_IR) that governs the rate of ionisation and recombination relative to the collisional time. A detailed derivation is included in Appendix A.

In Eq. (7) for the plasma energy, the recombination term has a factor of 1/2, which ensures that the ionisation and recombination terms are in equilibrium at T_n = T_p, consistent with Leake et al. (2012), Popescu Braileanu et al. (2019). We note that the equations used in Singh et al. (2019) are missing this factor, implying that their ionisation and recombination equilibrium exists at T_n = 2T_p, whereas their collisional equilibrium occurs at T_n = T_p, leading to non-trivial equilibrium conditions.

The energy equation for the plasma species includes an additional term to approximate the ionisation potential and the energy removed from the system due to ionisation, denoted in Eq. (7) as ϕ_I. We make several assumptions for the ionisation potential:

– On recombination only the proton thermal energy is transferred to the neutral fluid. The thermal energy previously held by the recombined electron remains in the electron fluid (and with that the plasma as a whole) as it is transferred to the free electron involved in the three-body recombination.

– The electron is assumed to recombine to an arbitrary level of the atom.

– We assume all the other energy involved in the recombination process (coming from the work done on the recombining electron by the proton electric field) is released as a photon.

– The medium is assumed to be optically thin.

– The electron in any neutral cascades down from higher levels to the ground state by releasing photons.

– To ionise the neutral from the ground state through collisions with another electron, work has to be done by the free electron to release the bound electron, resulting in an energy loss. As we assume this happens to electrons in the ground state, the work done is 13.6 eV.

All these assumptions allow us to avoid the modelling of the internal structure of the atom and have a detailed understanding of the energy balance between thermal and photon energy involved in the recombination process.

In dimensional form, the energy removed from the plasma due to the ionisation potential is calculated as:

$$\begin{array}{c}\hfill {\varphi }_{I,\dim }=13.6{\mathrm{\Gamma }}_{\mathrm{ion},\dim }{n}_{\mathrm{n}},\end{array}$$(16)

where n_n is the neutral number density, 13.6 is the ionisation energy of hydrogen (in electron volts), and Γ_ion, dim is the dimensional ionisation rate.

The non-dimensional form used in our simulations is therefore calculated as

$$\begin{array}{c}\hfill {\varphi }_I={\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}\widehat{\varphi },\end{array}$$(17)

where $\widehat{\varphi }$ is a constant factor such that this ionisation potential is consistent with the rest of the normalisation. If the system is normalised to c_s = 1, ρ_tot = 1 then $\widehat{\varphi }=\frac{13.6}{K_{\mathrm{B}}\gamma T_0}$. If the system is normalised to V_A = 1, then $\widehat{\varphi }=\frac{13.6\beta }{2K_{\mathrm{B}}{T}_0}$.

An arbitrary heating term A_heat is also included to obtain initial equilibrium. The form of this is equal to the ionisation potential energy loss using the equilibrium quantities, namely:

$$\begin{array}{c}\hfill A_{\mathrm{heat}}={\mathrm{\Gamma }}_{\mathrm{ion}}(t=0){\rho }_{\mathrm{n}}(t=0)\widehat{\varphi }.\end{array}$$(18)

This heating term is constant throughout the simulation.

2.2. Ionisation equilibrium conditions

When ionisation and recombination are studied (with ionisation potential neglected), the ionisation equilibrium can be easily determined from the continuity equation as:

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}={\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}},\hfill \end{array}$$(19)

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{ion}}/{\mathrm{\Gamma }}_{\mathrm{rec}}={\rho }_{\mathrm{p}}/{\rho }_{\mathrm{n}}.\hfill \end{array}$$(20)

Therefore, the only constraint for ionisation equilibrium is that the ratio of ionisation to recombination rates equals the ratio of plasma to neutral densities. As such, ionisation equilibrium relies on the ratio of densities, and not the specific density values.

Including the ionisation potential has an interesting effect on the equilibrium states that can be achieved by the system. For an equilibrium, we require that the ionisation potential and heating terms balance:

$$\begin{array}{c}\hfill \widehat{\varphi }{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}=\widehat{\varphi }{\mathrm{\Gamma }}_{\mathrm{ion}}(t=0){\rho }_{\mathrm{n}}(t=0).\end{array}$$(21)

Assuming a constant ionisation energy $\widehat{\varphi }$, combined with the equilibrium continuity equation we find

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}={\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}={\mathrm{\Gamma }}_{\mathrm{ion}}(t=0){\rho }_{\mathrm{n}}(t=0)=\mathrm{const}.\end{array}$$(22)

Now, the ionisation and recombination rates can be expressed as a function of temperature and a plasma density:

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mathrm{ion}}& =G(T){\rho }_{\mathrm{p}},\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\mathrm{rec}}& =F(T){\rho }_{\mathrm{p}}.\hfill \end{array}$$(24)

Hence Eq. (22) becomes:

$$\begin{array}{c}\hfill G(T){\rho }_{\mathrm{p}}{\rho }_{\mathrm{n}}=F(T){\rho }_{\mathrm{p}}^2=\mathrm{const}.\end{array}$$(25)

This is an interesting result since it states that for a given temperature, an ionisation equilibrium only exists at specific density values (rather than density ratios when the ionisation potential term is neglected).

In terms of bulk density ρ, Eq. (25) becomes

$$\begin{array}{c}\hfill F(T){\xi }_i^2{\rho }^2=\mathrm{const}.,\end{array}$$(26)

where the ionisation fraction ξ_i can be expressed as a function of temperature only using the conservation of mass equation:

$$\begin{array}{c}\hfill {\xi }_i=\frac{1}{\frac{F(T)}{G(T)}+1}.\end{array}$$(27)

The behaviour of $F(T){\xi }_i^2$ and ξ_i with respect to temperature are shown in Fig. 1. The shape of $F(T){\xi }_i^2$ implies that there are at most two potential temperatures that satisfy Eq. (25). In the range 5000 K <  T <  20 000 K the function is monotonically increasing meaning that a compression of the system requires a decrease in temperature. As such, for a compressible shock one may expect the downstream temperature to be less than the upstream temperature, in contrast to the MHD result where temperature increases across a shock. Above this temperature range, $F(T){\xi }_i^2$ is monotonically decreasing, meaning that a compressible shock has an increase in temperature as is expected for MHD. The change in gradient of the function occurs when the ionisation fraction approaches 1. As such, this equation is consistent with standard MHD theory for a fully ionised plasma, however introduces new behaviour for the temperatures lower than T ≈ 20 000 K when the plasma is partially ionised.

Fig. 1. $F(T){\xi }_i^2$ as a function of temperature in dimensional form. The red line shows the ionisation fraction ξi for the associated temperature.

In this paper, we study the partially ionised range for T <  20 000 K. Above this limit, we encounter limitations to our model. A more realistic ionisation energy model would be required to account for the increased likelihood of exciting a higher level when the temperature is higher. Limitations also arise from the numerical model when the neutral density and pressures become too low and the treatment of the neutral species as a fluid breaks down.

3. Analytic solutions to shock jump equations

3.1. MHD

For the MHD equations, the shock jump conditions can be derived by considering the equations upstream (^u) and downstream (^d) of the shock in the deHoffmann-Teller shock frame (zero electric field either side of the shock). The conservative nature of the MHD equations means that for a steady-state shock in the deHoffmann-Teller frame, mass, momentum and energy can be equated sufficiently upstream and downstream. The MHD equations have no source terms and hence can be integrated across the shock to become:

$$\begin{array}{cc}& {\left[\rho v_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(28)

$$\begin{array}{cc}& {[\rho v_x^2+P+\frac{B_y^2}{2}]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(29)

$$\begin{array}{cc}& {[\rho v_x{v}_y-B_x{B}_y]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(30)

$$\begin{array}{cc}& {\left[v_x(\frac{\gamma }{\gamma -1}P+\frac{1}{2}\rho v^2)\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(31)

$$\begin{array}{cc}& {\left[B_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(32)

$$\begin{array}{cc}& {[v_x{B}_y-v_y{B}_x]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(33)

where this notation means that

$$\begin{array}{c}\hfill {\left[Q\right]}_{\mathrm{d}}^{\mathrm{u}}\equiv Q^{\mathrm{u}}-Q^{\mathrm{d}}\end{array}$$(34)

for any quantity Q.

Rearranging these equations, it is possible to derive a single equation for possible shock transitions for a given upstream plasma-β and upstream angle of magnetic field θ, that relates the upstream and downstream Alfvén Mach numbers (Hau & Sonnerup 1989):

$$\begin{array}{cc}\hfill A_x^{\mathrm{u}2}& =[A_x^{\mathrm{d}2}(\frac{\gamma -1}{\gamma }(\frac{\gamma +1}{\gamma -1}-{tan}^2\theta ){(A_x^{\mathrm{d}2}-1)}^2\hfill \\ \hfill & +{tan}^2\theta (\frac{\gamma -1}{\gamma }{A}_x^{\mathrm{d}2}-1)(A_x^{\mathrm{d}2}-2))-\frac{\beta }{{cos}^2\theta }{(A_x^{\mathrm{d}2}-1)}^2]\hfill \\ \hfill & /[\frac{\gamma -1}{\gamma }\frac{{(A_x^{\mathrm{d}2}-1)}^2}{{cos}^2\theta }-A_x^{\mathrm{d}2}{tan}^2\theta (\frac{\gamma -1}{\gamma }{A}_x^{\mathrm{d}2}-1)].\hfill \end{array}$$(35)

This is shown in Fig. 2 for a choice of β = 0.1, θ = π/4 for the upstream plasma. A trivial solution exists where the upstream and downstream Alfvén Mach numbers are equal. The non-trivial solution intersects the trivial solution at three points where the upstream velocity is equal to the sound speed, the Alfvén speed, and the fast speed. Solutions below the trivial solution (where the where the velocity is higher downstream than upstream) are non-physical because these require that the entropy decreases across the shock. Additional details regarding this MHD solution can be found in Hau & Sonnerup (1989), Snow & Hillier (2019).

Fig. 2. Solutions to the shock jump equations for the IRIP (red) and MHD (black) equations relating the upstream (u) and downstream (d) Alfvén Mach numbers Ax. The trivial solution ($A_x^{\text{d}}=A_x^{\text{u}}$) exists for both sets of equations. The parameters used for this plot are T0 = 10 000 K, β = 0.1, θ = π/4, and γ = 5/3.

In the MHD model, there is a limitation on the compressionality of the system that arises from the energy equation, namely

$$\begin{array}{c}\hfill 1\le r\le \frac{\gamma +1}{\gamma -1}.\end{array}$$(36)

For γ = 5/3, this leads to a maximum compression ratio of r = 4.

3.2. Collisional ionisation and recombination (no ionisation potential – IR model)

When the two-fluid equations are used with ionisation and recombination (IR) in the absence of ionisation potential terms, one can again find an analytical solution for the jumps across the shock. For the deHoffmann-Teller shock frame, we set the time derivatives to be zero and assume the electric field is zero across the shock.

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}})=S_{\mathrm{mass}},\hfill \end{array}$$(37)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}+P_{\mathrm{n}}\mathbf{I})=S_{\mathrm{mom}},\hfill \end{array}$$(38)

$$\begin{array}{cc}& \mathrm{\nabla }·\left[{\mathit{v}}_{\mathrm{n}}(e_{\mathrm{n}}+P_{\mathrm{n}})\right]=S_{\mathrm{eng}},\hfill \end{array}$$(39)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}})=-S_{\mathrm{mass}},\hfill \end{array}$$(40)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}+P_{\mathrm{p}}\mathbf{I}-\mathit{BB}+\frac{{\mathit{B}}^2}{2}\mathbf{I})=-S_{\mathrm{mom}},\hfill \end{array}$$(41)

$$\begin{array}{cc}& \mathrm{\nabla }·[{\mathit{v}}_{\mathrm{p}}(e_{\mathrm{p}}+P_{\mathrm{p}})-({\mathit{v}}_{\mathrm{p}}\times \mathit{B})\times \mathit{B}]=-S_{\mathrm{eng}},\hfill \end{array}$$(42)

$$\begin{array}{cc}& \mathrm{\nabla }\times ({\mathit{v}}_{\mathrm{p}}\times \mathit{B})=0,\hfill \end{array}$$(43)

$$\begin{array}{cc}& S_{\mathrm{mass}}={\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}},\hfill \end{array}$$(44)

$$\begin{array}{cc}& S_{\mathrm{mom}}=-{\alpha }_c{\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}({\mathit{v}}_{\mathrm{n}}-{\mathit{v}}_{\mathrm{p}})+{\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}},\hfill \end{array}$$(45)

$$\begin{array}{cc}& S_{\mathrm{eng}}=-{\alpha }_c{\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}[\frac{1}{2}({\mathit{v}}_{\mathrm{n}}^2-{\mathit{v}}_{\mathrm{p}}^2)+\frac{3}{2}(\frac{P_{\mathrm{n}}}{{\rho }_{\mathrm{n}}}-\frac{1}{2}\frac{P_{\mathrm{p}}}{{\rho }_{\mathrm{p}}})]\hfill \\ \hfill & +\frac{1}{2}({\mathrm{\Gamma }}_{\mathrm{rec}}{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}^2-{\mathrm{\Gamma }}_{\mathrm{ion}}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}^2)\hfill \\ \hfill & +\frac{1}{(\gamma -1)}(\frac{1}{2}{\mathrm{\Gamma }}_{\mathrm{rec}}{P}_{\mathrm{p}}-{\mathrm{\Gamma }}_{\mathrm{ion}}{P}_{\mathrm{n}}).\hfill \end{array}$$(46)

However, these equations have a non-zero right hand side due to the source terms and therefore cannot be integrated in their current form.

Since the IR equations are conservative, the source terms are equal and opposite for the plasma and neutral equations, that is, momentum, energy, and mass that is lost or gained by the neutrals is gained or lost by the plasma. As such, the equations can be added together to remove the source terms and integrated. The partial densities and pressures can also be expressed in terms of the neutral fraction (ξ_n), bulk density (ρ_B) and bulk pressure (P_B).

$$\begin{array}{cc}& {[{\xi }_{\mathrm{n}}{\rho }_{\mathrm{B}}{v}_{\mathit{nx}}+(1-{\xi }_{\mathrm{n}}){\rho }_{\mathrm{B}}{v}_{\mathit{px}}]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(47)

$$\begin{array}{cc}& [{\rho }_{\mathrm{B}}{v}_{\mathit{nx}}^2+\frac{{\xi }_{\mathrm{n}}}{2-{\xi }_{\mathrm{n}}}{P}_{\mathrm{B}}+(1-{\xi }_{\mathrm{n}}){\rho }_{\mathrm{B}}{v}_{\mathit{px}}^2\hfill \\ \hfill & {+\frac{2(1-{\xi }_{\mathrm{n}})}{2-{\xi }_{\mathrm{n}}}{P}_{\mathrm{B}}+\frac{B_y^2}{2}]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(48)

$$\begin{array}{cc}& {[{\xi }_{\mathrm{n}}{\rho }_{\mathrm{B}}{v}_{\mathit{nx}}{v}_{\mathit{ny}}+(1-{\xi }_{\mathrm{n}}){\rho }_{\mathrm{B}}{v}_{\mathit{px}}{v}_{\mathit{py}}-B_x{B}_y]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(49)

$$\begin{array}{cc}& [v_{\mathit{nx}}(\frac{\gamma }{\gamma -1}\frac{{\xi }_{\mathrm{n}}}{2-{\xi }_{\mathrm{n}}}{P}_{\mathrm{B}}+\frac{1}{2}{\xi }_{\mathrm{n}}{\rho }_{\mathrm{B}}{v}_n^2)\hfill \\ \hfill & {+v_{\mathit{px}}(\frac{\gamma }{\gamma -1}\frac{2(1-{\xi }_{\mathrm{n}})}{2-{\xi }_{\mathrm{n}}}{P}_{\mathrm{B}}+\frac{1}{2}(1-{\xi }_{\mathrm{n}}){\rho }_{\mathrm{B}}{v}_p^2)]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(50)

$$\begin{array}{cc}& {\left[B_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(51)

$$\begin{array}{cc}& {[v_{\mathit{px}}{B}_y-v_{\mathit{py}}{B}_x]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(52)

for upstream (superscript u) and downstream (superscript d) states.

To simplify the equations further one can make some assumptions regarding the medium far away from the shock. Sufficiently upstream and downstream of the shock, one would expect that the drift velocity is zero (v_nx = v_px, v_ny = v_py). This greatly simplifies the equations:

$$\begin{array}{cc}& {\left[{\rho }_{\mathrm{B}}{v}_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(53)

$$\begin{array}{cc}& {[{\rho }_{\mathrm{B}}{v}_x^2+P_{\mathrm{B}}+\frac{B_y^2}{2}]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(54)

$$\begin{array}{cc}& {[{\rho }_{\mathrm{B}}{v}_x{v}_y-B_x{B}_y]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(55)

$$\begin{array}{cc}& {\left[v_x(\frac{\gamma }{\gamma -1}{P}_{\mathrm{B}}+\frac{1}{2}{\rho }_{\mathrm{B}}{v}^2)\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(56)

$$\begin{array}{cc}& {\left[B_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(57)

$$\begin{array}{cc}& {[v_x{B}_y-v_y{B}_x]}_{\mathrm{d}}^{\mathrm{u}}=0.\hfill \end{array}$$(58)

This set of equations is identical to the MHD shock equations meaning that the analytic solution is identical to the MHD solution. This only holds sufficiently upstream and downstream of the shock since two-fluid effects lead to a finite shock width and substructure (Snow & Hillier 2019).

3.3. Collisional ionisation, recombination, and ionisation potential (IRIP model)

Including ionisation potential and arbitrary heating for the IRIP model introduces two new terms into the plasma energy equation:

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}})=S_{\mathrm{mass}},\hfill \end{array}$$(59)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}+P_{\mathrm{n}}\mathbf{I})=S_{\mathrm{mom}},\hfill \end{array}$$(60)

$$\begin{array}{cc}& \mathrm{\nabla }·\left[{\mathit{v}}_{\mathrm{n}}(e_{\mathrm{n}}+P_{\mathrm{n}})\right]=S_{\mathrm{eng}},\hfill \end{array}$$(61)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}})=-S_{\mathrm{mass}},\hfill \end{array}$$(62)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}+P_{\mathrm{p}}\mathbf{I}-\mathit{BB}+\frac{{\mathit{B}}^2}{2}\mathbf{I})=-S_{\mathrm{mom}},\hfill \end{array}$$(63)

$$\begin{array}{cc}& \mathrm{\nabla }·[{\mathit{v}}_{\mathrm{p}}(e_{\mathrm{p}}+P_{\mathrm{p}})-({\mathit{v}}_{\mathrm{p}}\times \mathit{B})\times \mathit{B}]=-S_{\mathrm{eng}}-{\varphi }_I+A_{\mathrm{heat}},\hfill \end{array}$$(64)

$$\begin{array}{cc}& \mathrm{\nabla }\times ({\mathit{v}}_{\mathrm{p}}\times \mathit{B})=0,\hfill \end{array}$$(65)

with S_mass, S_mom, S_eng given by Eqs. (44)–(46). Since the ionisation potential ϕ_I and the arbitrary heating A_heat only exist in the plasma energy equation, adding together the equations for neutral and plasma species results in:

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}+{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{n}})=0,\hfill \end{array}$$(66)

$$\begin{array}{cc}& \mathrm{\nabla }·({\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}+P_{\mathrm{n}}\mathbf{I}+{\rho }_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}{\mathit{v}}_{\mathrm{p}}+P_{\mathrm{p}}\mathbf{I})=0,\hfill \end{array}$$(67)

$$\begin{array}{cc}& \mathrm{\nabla }·[{\mathit{v}}_{\mathrm{n}}(e_{\mathrm{n}}+P_{\mathrm{n}})+{\mathit{v}}_{\mathrm{p}}(e_{\mathrm{p}}+P_{\mathrm{p}})-({\mathit{v}}_{\mathrm{p}}\times \mathit{B})\times \mathit{B}]=-{\varphi }_I+A_{\mathrm{heat}},\hfill \end{array}$$(68)

$$\begin{array}{cc}& \mathrm{\nabla }\times ({\mathit{v}}_{\mathrm{p}}\times \mathit{B})=0,\hfill \end{array}$$(69)

$$\begin{array}{cc}& \mathrm{\nabla }·\mathit{B}=0.\hfill \end{array}$$(70)

Now it is fairly clear that all except the energy Eq. (68) can be trivially integrated across the shock. Since the energy equation is non-conservative, we cannot use this equation to derive any shock jump conditions. A semi-analytical solution can still be found by replacing the energy equation with the relation that the ionisation potential and arbitrary heating terms must balance sufficiently upstream and downstream of the shock, Eq. (26). The assumption of zero drift sufficiently upstream and downstream is also applied to get a solution for the bulk fluid. Our shock jump equation for the IRIP model are then:

$$\begin{array}{cc}& {\left[{\rho }_{\mathrm{B}}{v}_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(71)

$$\begin{array}{cc}& {[{\rho }_{\mathrm{B}}{v}_x^2+P_{\mathrm{B}}+\frac{B_y^2}{2}]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(72)

$$\begin{array}{cc}& {[{\rho }_{\mathrm{B}}{v}_x{v}_y-B_x{B}_y]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(73)

$$\begin{array}{cc}& {\left[B_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(74)

$$\begin{array}{cc}& {[v_x{B}_y-v_y{B}_x]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(75)

$$\begin{array}{cc}& {\left[F(T){\xi }_i^2{\rho }_{\mathrm{B}}^2\right]}_{\mathrm{d}}^{\mathrm{u}}=0.\hfill \end{array}$$(76)

The ideal gas law can be used to relate the temperature and pressure and close the system. Therefore, all possible shock transitions can be obtained by solving the following equations:

$$\begin{array}{cc}& \frac{T^{\mathrm{d}}}{T^{\mathrm{u}}}=\frac{A_x^{d2}}{A_x^{u2}}[1+\frac{2}{\beta (1+{tan}^2(\theta ))}\times \hfill \\ \hfill & [A_x^{u2}-A_x^{d2}+\frac{{tan}^2(\theta )}{2}(1-{\left(\frac{A_x^{u2}-1}{A_x^{d2}-1}\right)}^2)]],\hfill \end{array}$$(77)

$$\begin{array}{cc}& \frac{F(T^{\mathrm{d}})}{F(T^{\mathrm{u}})}{\left(\frac{F(T^{\mathrm{u}})/G(T^{\mathrm{u}})+1}{F(T^{\mathrm{d}})/G(T^{\mathrm{d}})+1}\right)}^2=\frac{A_x^{d4}}{A_x^{u4}}.\hfill \end{array}$$(78)

A derivation of these equations is included in Appendix B. Equation (77) gives a temperature jump across the shock (using the ideal gas law and Eqs. (71)–(75)) and can be solved numerically in conjunction with Eq. (78) to get the possible shock transitions for a given upstream plasma-β, angle of magnetic field, and reference temperature. A solution for T^u = 10 000 K is given in Fig. 2. The trivial solution (A^u = A^d) still satisfies the IRIP jump conditions and has three intersections with the non-trivial solution.

Both the IRIP and MHD models have an intersection when the upstream and downstream Alfvén Mach numbers are unity which is the rotational discontinuity. For A^d <  1 there are slow (A^u <  1) and intermediate (A^d >  1, A^u <  1) shock solutions. Interestingly, these solution imply much greater compression than the MHD solution since the upstream Alfvén Mach number for the IRIP case is much larger relative to the downstream Alfvén Mach number. Unlike the MHD model, there is no restriction on the compressionality of the system.

In this paper, the focus is on a switch-off slow-mode shock which occurs when the upstream Alfvén Mach number is unity. From Fig. 2 it can easily be seen that the switch-off shock in the IRIP model will have a much smaller downstream shock velocity than in the MHD case. It will also have a much greater compressional ratio and will have a decrease in temperature across the shock.

4. Numerical model

Numerical simulations are performed using the (PIP) code (Hillier et al. 2016) that solves for two fluids (neutral, and ion+electron) in normalised units. The (PIP) code has been developed since earlier editions and now features terms to account for ionisation, recombination, and ionisation potential energy. Specifically, Eqs. (1)–(10) are solved numerically using a first order HLLD scheme. The source terms for collisions and ionisation/recombination are integrated explicitly. The system is resolved by 256 000 grid cells.

Initial conditions. The initial conditions are designed to mimic a slow-mode shock produced by fast magnetic reconnection. Specifically, our initial conditions consist of a thermal equilibrium and a discontinuity in magnetic field across the origin:

$$\begin{array}{cc}& B_x=0.1,\hfill \end{array}$$(79)

$$\begin{array}{cc}& B_y=-1.0(x>0),1.0(x<0),\hfill \end{array}$$(80)

$$\begin{array}{cc}& {\rho }_{\mathrm{n}}={\xi }_{\mathrm{n}}{\rho }_{\mathrm{t}},\hfill \end{array}$$(81)

$$\begin{array}{cc}& {\rho }_{\mathrm{p}}={\xi }_{\mathrm{i}}{\rho }_{\mathrm{t}}=(1-{\xi }_{\mathrm{n}}){\rho }_{\mathrm{t}},\hfill \end{array}$$(82)

$$\begin{array}{cc}& P_{\mathrm{n}}=\frac{{\xi }_{\mathrm{n}}}{{\xi }_{\mathrm{n}}+2{\xi }_{\mathrm{i}}}{P}_{\mathrm{t}}=\frac{{\xi }_{\mathrm{n}}}{{\xi }_{\mathrm{n}}+2{\xi }_{\mathrm{i}}}\beta \frac{B_0^2}{2},\hfill \end{array}$$(83)

$$\begin{array}{cc}& P_{\mathrm{p}}=\frac{2{\xi }_i}{{\xi }_{\mathrm{n}}+2{\xi }_{\mathrm{i}}}{P}_{\mathrm{t}}=\frac{2{\xi }_{\mathrm{i}}}{{\xi }_{\mathrm{n}}+2{\xi }_{\mathrm{i}}}\beta \frac{B_0^2}{2},\hfill \end{array}$$(84)

where ξ_n and ξ_i are the neutral and ion fractions respectively. These initial conditions have been studied in the absence of ionisation and recombination in Hillier et al. (2016), Snow & Hillier (2019).

In normalised form, the equilibrium neutral fraction depends on the reference temperature T₀. Taking the continuity of mass (Eq. (1)) and rearranging, we get Γ_rec/Γ_ion = ρ_n/ρ_p. Therefore, the equilibrium neutral fraction ξ_n can be calculated using

$$\begin{array}{c}\hfill {\xi }_{\mathrm{n}}=\frac{{\mathrm{\Gamma }}_{\mathrm{rec}}/{\mathrm{\Gamma }}_{\mathrm{ion}}}{({\mathrm{\Gamma }}_{\mathrm{rec}}/{\mathrm{\Gamma }}_{\mathrm{ion}}+1)}.\end{array}$$(85)

The initial neutral fraction is controlled by varying the reference temperature T₀.

We consider three different cases in this work. Firstly, for reference, the MHD model is considered where the plasma is fully ionised. Secondly, we have the IRIP model for a two-fluid system where ionisation, recombination, ionisation potential, and arbitrary heating is included. Finally, again for reference, an IR model is used which is multi-fluid with ionisation and recombination but neglecting the ionisation potential and heating terms. The main focus of this paper is the IRIP model.

5. Results

For the initial simulations, we used T₀ = 10 000 K, which leads to an equilibrium neutral fraction of ξ_n ≈ 0.99734. We set the recombination rate to be Γ_rec = 10⁻³ and the collisional frequency was set to unity. Physically, this means that collisions occur on a timescale of unity, and recombination occurs on a timescale of 1000. The ionisation timescale is calculated from the recombination timescale and the normalisation temperature. Using T₀ = 10 000 K, the initial ionisation rate is Γ_ion = 2.66 × 10⁻⁶. The recombination and ionisation rates depend on the instantaneous temperature of the system which varies in time and hence the rates also vary.

5.1. Equilibrium solution (t = 40 000)

After 40 000 collisional times, the system is sufficiently steady. The rarefaction wave and shock have separated, and variables between these features have very small gradients. As such, an (approximate) equilibrium state has been reached. Figure 3 shows this pseduo-steady state for the IRIP model compared to the MHD solution. For comparison, an IR simulation is also presented which does not include the ionisation potential and heating terms.

Fig. 3. Context figure for the equilibrium state in the MHD (black), IRIP (solid) and IR (dashed) cases. The blue (red) line is for the plasma (neutral) species for the two-fluid cases.

The first result to state is that the IRIP result is very different to the MHD and IR solutions. This is not unexpected for the reasons given in previous sections, namely that the energy equation is non-conserving, equilibrium only exists at certain densities, and we expect greater compression across the shock. In particular, the inflow to the shock front is much larger in the IRIP case, and the propagation speed of the IRIP shock is much lower than the MHD case.

5.1.1. Inflow

To understand the inflow to the shock, the physics of the rarefaction wave must be understood. Figure 4 shows the energy terms across the system. For the MHD rarefaction wave, magnetic and thermal energy are converted to kinetic energy driving an inflow towards the slow-mode shock.

Fig. 4. Energies across the system for the MHD (black) and IRIP cases (red neutrals, blue plasma). The quantity Aheat − ϕI shows the balance between the heating and loss terms in the energy equation, where red denotes net energy addition (heating) and blue is net energy loss (cooling).

For the IRIP case, there is a similar thermal energy reduction, and a much larger loss of magnetic energy, compared to MHD, see Fig. 4. The thermal energy reduction is mostly in the neutral species which constitutes the majority of the bulk fluid. The rarefaction wave is expansive, meaning that there is a decrease in total density across the rarefaction wave. The density reduction is in the neutral species and plasma density is roughly the same in the static medium and inflow regions, see Fig. 3.

In the IRIP case, the total density decrease across the rarefaction wave results in a slight increase in temperature due to the requirement that for a steady (equilibrium) solution to exist,

$$\begin{array}{c}\hfill G(T){\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}=F(T){\rho }_{\mathrm{p}}^2={\mathrm{\Gamma }}_{\mathrm{ion}}(t=0){\rho }_n(t=0)=\mathrm{const}.\end{array}$$(86)

The decreased total density necessitates an increased ionisation rate which can only be obtained through a temperature increase of the plasma species (since both the ionisation and recombination rates rely on the plasma temperature). The temperature increase is provided through the arbitrary heating term, which, as the system undergoes expansion, is larger than the energy lost through ionisation potential in the rarefaction wave.

The leading edge of the rarefaction wave propagates at the same speed in the MHD and IRIP models. However, the IRIP rarefaction wave has a far greater expansion than the MHD case. The rarefaction wave in all simulations drives velocity towards the slow-mode shock at the Alfvén speed.

Across the rarefaction wave, the IR jumps in bulk quantities match the MHD solution (see Fig. 3). This matches the theoretical result from Sect. 3.2 that, since the IR equation reduce to the MHD equations sufficiently upstream and downstream of a feature, the MHD result should be obtained.

5.1.2. Slow-mode shock

The inflow conditions for the slow-mode shock in the IRIP case are different to the MHD case, and, as such, the jump conditions are expected to be different. Figure 5 shows the velocity (top left), temperature (top right), density (lower left), and ionisation and recombination rates (lower right) across the slow mode shock. Since the system has been evolved for 40 000 collisional times the slow-mode shock is fairly localised on the x/t axis but is well resolved by our grid resolution.

Fig. 5. Close-up of the slow-mode shock for the IRIP model showing vx velocity (top left), temperature (top right), and density (lower left) for plasma (blue) and neutral (red) species. Lower right panel: ionisation (orange) and recombination (green) rates.

The rarefaction wave drives inflow towards the slow-mode shock. In the IRIP case, the v_x velocity is far greater than in the MHD case, however, in the shock frame, both systems have an inflow Alfvén Mach number of one, see Fig. 6. The specific shock type here is a switch-off slow-mode shock where the perpendicular magnetic field reduced to zero across the shock. Switch-off shocks are the most extreme slow-mode shocks and occur when the inflow Alfvén Mach number is unity.

Fig. 6. Alfven Mach numbers for the MHD (black) and plasma (blue) based on the bulk density. The red line shows the neutral sonic Mach number.

Fig. 7. Time series for the IRIP case after 1 (top), 10 (middle) and 100 (lower) collisional times showing the vx velocity (left), temperature (centre), and ionisation and recombination rates (right).

The velocity structure across the shock (Fig. 5) displays similar structures to the PIP case in Hillier et al. (2016), Snow & Hillier (2019) namely that the neutral species shocks first, followed by a gradual coupling to the plasma resulting in a finite width to the slow-mode shock. Two-fluid effects result in a finite shock width where the species separate and interesting physics occurs.

In the absence of ionisation and recombination, a key feature of two-fluid interaction in shocks is an overshoot of the neutral velocity (Hillier et al. 2016; Snow & Hillier 2019). The large drift between the ions and neutrals inside the shock results in frictional heating leading to a Sedov-Taylor-like expansion of the neutral species. The neutral overshoot is not present in the IRIP model presented here (see Fig. 5). The initial discontinuity drives the plasma only, which then couples to the neutral species. At early times in the simulation (t = 10) an overshoot in the neutrals exists but ionisation prevents it from being a sustained feature of the system (the time evolution of the system is discussed further in Sect. 5.2).

At the leading edge of the shock, a sonic shock occurs in the neutral species, see Fig. 6. The sonic shock features a temperature increase (due to adiabatic heating) in the neutral species. Collisional coupling transfers this heat from the neutrals to the plasma and hence leads to a local enhancement of ionisation. Behind this, the plasma produces a slow-mode shock, with another localised temperature increase and ionisation rate.

The post-shock region has a lower temperature than the pre-shock region. This is related to the relation:

$$\begin{array}{c}\hfill G(T){\rho }_{\mathrm{n}}{\rho }_{\mathrm{p}}=F(T){\rho }_{\mathrm{p}}^2={\mathrm{\Gamma }}_{\mathrm{ion}}(t=0){\rho }_{\mathrm{n}}(t=0).\end{array}$$(87)

Both the slow-mode and the sonic shock that forms within it are compressional, thus the increase in density leads to a decrease in temperature from the ionisation potential energy losses. An interesting corollary of the temperature reduction across the shocks is that observational signatures of chromospheric shocks may manifest as a decrease in intensity due to the temperature drop. We note that the slow-mode shock features an increase in both density and pressure across the shock.

The Mach numbers for the system in the shock frame are shown in Fig. 6. The MHD case is a slow-mode switch-off shock, which is the strongest possible slow-mode shock where the inflow velocity is the Alfvén speed. In the IRIP case, despite the very different inflow velocity and magnetic field, we still have an Alfvénic inflow and a super- to sub-slow transition across the shock indicating the same switch-off slow-mode shock as in the MHD case. In the two-fluid case, the species decouple and recouple around the shock front. At the leading edge of the shock the species separate and a sonic shock forms in the neutral species, as seen in Hillier et al. (2016). Furthermore, the plasma exceeds the Alfvén speed inside the finite width of the shock, resulting in an intermediate shock. The magnitude of the magnetic field reversal associated with the intermediate shock decreases with time due to the balance between the neutral pressure and the Lorentz force, as described in Snow & Hillier (2019). After 40 000 collisional times, the magnetic field reversal is very small.

Using the upstream properties from the numerical simulation in the analytical solution, the inflow Alfvén Mach number pairs remarkably well (see Fig. B.1). This gives us confidence in the applicability of the semi-analytical solution described in Sect. 3.3.

5.2. Early times

Now that the equilibrium state is understood, the system can be considered on transient timescales. For this simulation, the initial recombination rate was set to Γ_rec = 10⁻³ of the collisional time. As such, one expects the system to reach a collisional equilibrium long before an ionisation equilibrium is reached.

After 1 collisional time the neutrals have experienced very few collisions. The initial conditions feature a discontinuity of magnetic field across the origin, triggering a fast-mode wave and a slow-mode shock in the plasma species only. At such an early time, the fast-mode and slow-mode are still interacting and have not decoupled. The slow-mode shock features a temperature increase and hence increased local collisions so the neutral species has a slight response. The increased temperature boosts the ionisation and recombination rates by approximately nine and two orders of magnitude respectively. As such, these effects become important even on such small timescales. The large increase in ionisation rate results in large energy losses through the ionisation potential in these regions.

After 10 collisional times, the neutral species has a clear response to the plasma motion due to collisional coupling. The profiles are still very distinct and far from equilibrium state. The ionisation and recombination effects are starting to have an effect on the system here, as can be seen by the post-shock cooling in the plasma temperature.

By 100 collisional times, similar structure begins to exist in the plasma and neutral species. A sonic shock can be seen in the neutral species. The post shock temperature is very similar to the case after 40 000 collisional times in Fig. 5. The pre-shock region is still far from equilibrium, with large drift velocities existing between the two species.

5.3. Steady solution

As time advances, the shock tends towards a steady solution, where the finite-width of the shock is determined by the diffusive terms of the system. Figure 8 shows the drift velocity throughout the slow-mode shock at different times, and shock width can be determined by the departure from zero drift (that is, a coupled fluid). The peak drift velocity is approximately 15% of the bulk Alfvén speed. After 10 000 collisional times, the shock width changes very little and the system can be considered to be a steady solution. The finite width of the shock is determined based on the physical diffusion mechanisms, here the ion-neutral collisions and ionisation, recombination and ionisation potential. 10 000 collisional times corresponds to 10 recombination times based on the background recombination rate of τ_IR = 10⁻³. However, inside the shock the local temperature increases result in a recombination rate a few orders of magnitude higher than the background value, see Fig. 5. Hence recombination (and ionisation) occurs on a far faster timescale inside and around the shock than in implied by background quantities. Taking the peak recombination rate inside the shock of 0.1, we can estimate that 10 000 collisional times corresponds to approximately 1000 recombination times. As such, the system has had plenty of time to reach ionisation/recombination equilibrium inside the shock.

Fig. 8. Drift velocity at different times for the IRIP simulation. The shocks are aligned using xs, which is the location of the minimum gradient of the plasma velocity, that is, the shock terminus.

The finite width of the shock is much smaller in the IRIP model than when ionisation and recombination is neglected. In the IRIP simulation, the finite width is approximately L = 35. For a simulation using the same initial neutral fraction but without ionisation and recombination, the shock width is approximately L = 414. The smaller finite width can be explained by the increased ion-fraction inside the shock shown in Fig. 5 (this was shown in Hillier et al. 2016, where a higher ion-fraction leads to a smaller shock). A consequence of the reduced shock width is that a package of fluid moving though the shock will lose less energy than it would through a wider shock. As such, cooling inside the shock becomes less important for a narrower shock. The post-shock cooling is still very important in removing energy from the system.

6. Parameter study

6.1. Recombination timescale

In the previous simulations the initial recombination rate was scaled to τ_IR = 10⁻³ of the collisional time. The initial ionisation rate is then a factor of the recombination time based on the reference temperature such that the system is in equilibrium. In this section we investigate changing the reference recombination rate τ_IR, which also changes the ionisation rate. In Sect. 5.2 it was found that the ionisation rate increases by approximately nine orders of magnitude at the slow-mode shock after 1 collisional time. Therefore, the rates are increased due to the initial temperature enhancements, resulting in the ionisation and recombination forces balancing before the collisional components equalise. As such, smaller reference recombination times are a primary focus in this section. All other system parameters are kept the same (T₀ = 10 000 K, β = 0.1).

Decreasing the recombination rate causes the shock to propagate more slowly and leads to a greater shock width, Fig. 9. Similar velocity structure exists in all simulations, that is, a neutral shock at the leading edge, with a sharp jump in plasma velocity at the rear of the shock. In the shock frame, the inflow conditions for all simulations is the Alfvén speed, and as such all simulations can be classified as having a switch-off slow-mode shock. Plotting the shock width as a function of the initial recombination rate, Fig. 10, one can see the shock width decreases as the recombination rate increases. For τ_IR <  10⁻³ the system is still evolving and the shock is thinning. The τ_IR = 10⁻² case has roughly the same shock thickness as the τ_IR = 10⁻³ case because the system has reached a steady state.

Fig. 9. Plasma (blue) and neutral (red) vx velocities for recombination timescales of 10−3, 10−5, 10−6, 10−7, from top to bottom.

Fig. 10. Finite width of the shock (black line) as a function of the initial recombination rates. Integrated cooling for a parcel of fluid travelling through the shock is shown by the red line.

The initial conditions produce a large temperature from the adiabatic heating of the slow-mode shock in the plasma (see Hillier et al. 2016). The temperature increase boosts the ionisation rate, which in turn increases the energy lost through ionisation potential. Hence decreasing the rate of recombination (and with it, ionisation) results in energy being removed from the system slower and the medium undergoes more gradual cooling, see Fig. 11. All simulations have the same restraint that the ionisation potential losses and arbitrary heating terms must balance for equilibrium and hence as distance from the shock front tends to infinity, the upstream temperature should be the same across all simulations. It can be concluded that all simulations are in collisional equilibrium (v_p = v_n sufficiently upstream and downstream) however they are not in ionisation-recombination equilibrium (Γ_ionρ_n ≠ Γ_recρ_p). The low ionisation rates for τ_IR = 10⁻⁷ result in a very gradual cooling and one would expect the equilibrium to be obtained on far longer timescales than studied here (end time t = 10 000). The decrease in cooling also affects the shock temperature since energy is being removed at a much lower rate. For τ = 10⁻⁷ the shock temperature approaches the downstream MHD temperature, see Fig. 11.

Fig. 11. Plasma (solid) and neutral (dashed) temperatures for different initial recombination rates. The orange dashed line shows the downstream MHD temperature.

The interior velocity structure is similar for the different initial recombination rates (Fig. 9) however the neutral fraction is very different, as shown in Fig. 12. In all cases, the neutral shock at the leading edge results in a decrease in ionisation fraction due to the increased losses from the adiabatic heating. Following this, there is a gradual rise in the ionisation fraction throughout the shock, reaching an apex at the shock terminus. Post-shock, the neutral fraction is much larger for smaller τ_IR due to the decreased cooling rate for these simulations, compared to the reference value. The snapshots were taken after 10 000 collisioal times however it is known that the τ_IR = 10⁻³ case only just reaches a steady state around this time. One would expect that decreasing the recombination rate by an order of magnitude increases the time required for the steady solution by an order of magnitude. Therefore, it would require 10⁸ collisional times for the τ_IR = 10⁻⁷ case to reach its steady solution.

Fig. 12. Ionisation fraction across the shock for different recombination rates.

A parcel of fluid travelling through the shock experiences an energy loss due to the ionisation potential loss term being larger than the arbitrary heating. The integrated cooling across the finite width of the shock can be calculated using:

$$\begin{array}{c}\hfill \mathrm{Cooling}=-{\mathrm{\Sigma }}_a^b(A_{\mathrm{heat}}-{\varphi }_I)\frac{\mathrm{\Delta }x}{v_{\mathit{xp}}-v_{\mathrm{s}}},\end{array}$$(88)

where a, b are the physical limits of the shocks finite width, Δx is the grid cell size, and v_s is the shock velocity. The integrated cooling through the shock is plotted in Fig. 10 for different initial recombination rates. Decreasing the initial recombination rate lowers the ionisation potential and hence results in less cooling through the shock, and a higher temperature is obtained for smaller reference recombination rates (Fig. 11). The shock also increases in width for smaller reference recombination rates. As such, the cooling of a parcel of fluid travelling across the shock is not as severe as one may expect; reducing the reference recombination time from 10⁻² to 10⁻⁷ only reduces the cooling across the shock by one order of magnitude.

To put this parameter study in the context of the solar atmosphere, in dimensional units, ionisation/recombination rates are typically of the order of 10⁻³ − 10⁻⁵ per second and are locally enhanced in shocks (Carlsson & Stein 2002). For a VALC chromosphere, the ion-neutral (ν_in) and neutral-ion (ν_ni) collisional frequencies are approximately in the ranges of 10³ <  ν_in <  10⁶ and 10⁻¹ <  ν_ni <  10² per second varying with height due to stratification (from Fig. 8 in Popescu Braileanu et al. 2019). We therefore expect the ionisation/recombination rates to be between 2 and 7 orders of magnitude smaller than the neutral-ion collisional frequency. As such, this parameter study covers the expected range of recombination timescales expected in the lower solar atmosphere.

6.2. Reference temperature

The reference temperature is a fairly important parameter since it governs the initial neutral fraction and relative ionisation and recombination rates. Previously, the temperature was set to T₀ = 10 000 K resulting in a predominantly neutral medium with a neutral fraction of ξ_n ≈ 0.99734. In this section, the reference temperature is modified resulting in different neutral fractions, as shown in Fig. 13. Other initial system parameters are kept as plasma-β = 0.1 and recombination rate as 10⁻³.

Fig. 13. vx plasma velocity for different reference temperatures showing the equilibrium state after 40 000 collisional times.

The system studied produces a rarefaction wave that drives inflow towards a switch-off slow-mode shock at the Alfvén speed. At the equilibrium state, both the ionisation potential and heating terms must balance resulting in unique solutions for the density jump. All solutions show a similar equilibrium state to the one discussed in detail in Sect. 5, see Fig. 13. Increasing the reference temperature results in a slower propagating slow-mode shock, which is also predicted by the semi-analytical solution. Other than that, there are relatively few differences for different reference temperatures. We note that the simulations performed are in the partially ionised regime. As one tends towards the fully ionised solution (that is, MHD), one would expect greater differences to occur.

7. Discussion

In the interstellar medium, molecules were discovered downstream of shocks that should have been disassociated due to the MHD predicted temperature inside the shock. However, including radiative losses in the system reduces the maximum obtained temperature in the shock, allowing these molecules to survive into the downstream medium (Draine & McKee 1993). The results presented in this paper show a similar effect where there is a reduction in the maximum temperature obtained in the shock due to radiative losses (see Fig. 11). A parcel of fluid moving through the shock experiences significant cooling and therefore lower temperature spectral lines may be present after the shock, compared with an MHD prediction. This result may be relevant to the IRIS result of De Pontieu et al. (2015) where the non-thermal line broadening in shocks was found to be greater than expected due to non-equilibrium ionisation. Non-equilibrium ionisation does not explain the full extend of the observational non-thermal broadening and hence the two-fluid effects discussed in this paper may be a way of increasing the non-thermal line width. Further study is needed to investigate the non-thermal line broadening due to two-fluid effects.

Within the finite-width of the shock, drift velocities reach up to 15% of the bulk Alfvén speed, as shown in Fig. 8. In terms of dimensional units, for plasma-β = 0.1, and a chromospheric sound speed of 8 km s⁻¹, this drift velocity equates to an Alfvén speed of approximately 20 km s⁻¹ and hence a drift velocity of approximately 3 km s⁻¹. This is fairly substantial and lends credibility to a previously predicted signatures of two-fluid effects in shocks, namely the electric field felt by the neutral species (Anan et al. 2017; Snow & Hillier 2019).

Including the ionisation, recombination, and ionisation potential terms narrows the finite width of the shock. However, the finite-width is far from discontinuous. A previously proposed effect of two-fluid propagating shocks is a large finite-width due to ion-neutral drag (Snow & Hillier 2020). This remains a potential observable for two-fluid shocks despite the narrowing of the finite-width, due to the large drift velocities present in our simulations.

8. Conclusions

In this paper we analyse the consequences of collisional ionisation, recombination, and ionisation potential on a partially ionised switch-off slow-mode shock. The ionisation potential term leads to a non-conservative energy equation that results in a equilibrium shock structure very different from the MHD case.

For an equilibrium to exist, the heating and cooling terms must balance. As such, a simple equation can be derived that must be satisfied for equilibrium that relates the upstream temperature and density to the arbitrary heating constant. In the limit of T <  20 000 K, the function is monotonically increasing with temperature, implying that a compressional shock results in upstream cooling, as opposed to the temperature increase obtained from MHD results. It also implies that for a given temperature, equilibrium only exists at specific density values.

Shock jump equations can be derived for the IRIP model (two-fluid equations with collsional ionisation, recombination, ionisation potential, and arbitrary heating) and a semi-analytic solution for possible shock transitions can be derived. The analytic solution for the shock Alfvén Mach number pairs well with the numerical solution for the slow-mode switch-off shock studied here. An empirical form of ionisation and recombination was used in this paper however the methodology could easily be applied to more advanced models.

The semi-analytic jump conditions derived state that for a compressible, steady-state, partially ionised, shock with ionisation potential included, the temperature should decrease across the shock, see Eq. (25). For an equilibrium to exist, the ionisation potential energy losses must balance the arbitrary heating term. In the partially ionised regime (T <  20 000 K), the equilibrium constraint necessitates that a steady-state compressible shock has a decrease in temperature across the shock. Our numerical simulation show this result, where sufficiently downstream of the shock, the temperature is lower than upstream, as shown by Fig. 5. This is contrary to the MHD and IR result where a compressible shock requires a temperature increase across the interface. Furthermore, including ionisation potential leads to cooling of the plasma within the finite-width of the shock and lowers the maximum temperature obtained, as shown in Fig. 11. As such, including ionisation potential fundamentally changes the behaviour of shocks in partially ionised plasmas.

An intermediate transition exists within the finite shock-width, as seen in the two-fluid simulations of Snow & Hillier (2019) without ionisation and recombination. This is a promising result that provides further evidence that additional shock transitions may occur within two-fluid shocks in the solar atmosphere.

Similar results are obtained when the recombination rate or reference temperature are altered. Decreasing the recombination rate results in the system taking longer to reach a steady state. Increasing the reference temperature reduces the ionisation fraction and results in a slower propagating shock. The slower propagation can be predicted from the semi-analytical solution.

In conclusion, non-conservative energy effects can have significant consequences for shocks in the partially ionised systems. In particular, for the equations studied here, a compressional shock will result in downstream cooling of the system if the medium is in the partially ionised regime. A semi-analytical solution for the equilibrium shock jumps can be derived which pairs well with the numerical simulations. The results of this paper are framed in the context of the lower solar atmosphere however are applicable to other partially ionised systems, for example, molecular clouds.

Acknowledgments

BS and AH are supported by STFC research grant ST/R000891/1. AH also acknowledges support by STFC Ernest Rutherford Fellowship grant number ST/L00397X/2 We would like to acknowledge and thank the useful discussions with Elena Khomenko, Nikola Vitas, Malcolm Druett and Giulio Del Zanna.

References

Appendix A: Ionisation and recombination rates

In dimensional units (of s⁻¹), the recombination and ionisation rates are given by the empirical rates in Voronov (1997), Smirnov (2003):

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{rec},\dim }=\frac{n_0{\rho }_{\mathrm{p}}}{\sqrt{T_{e0}{T}_{\mathrm{p}}/T_f}}2.6\times {10}^{-19},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{ion},\dim }=n_0{\rho }_{\mathrm{p}}\frac{2.91\times {10}^{-14}}{0.232+\chi }{\chi }^{0.39}{e}^{\chi },\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& \chi =13.6\frac{T_f}{T_{e0}{T}_{\mathrm{p}}},\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& T_f=\frac{1}{4}\beta \gamma \frac{2{\xi }_{p0}}{{\xi }_{n0}+2{\xi }_{p0}},\hfill \end{array}$$(A.4)

where n₀ is a reference dimensional density, and T_e0 is a reference dimensional temperature in electron volts. The instantaneous plasma temperature T_p is modified by T_f due to the normalisation in this paper being a Alfvén speed of unity. At time t = 0, T_p/T_f = 1.

From the conservation of mass equation, the initial (dimensionless) neutral fraction can be calculated from the ratio of the dimensional rates as:

$$\begin{array}{c}\hfill {\xi }_{\mathrm{n}}={\mathrm{\Gamma }}_{\mathrm{rec},\dim }/{\mathrm{\Gamma }}_{\mathrm{ion},\dim }/({\mathrm{\Gamma }}_{\mathrm{rec},\dim }/{\mathrm{\Gamma }}_{\mathrm{ion},\dim }+1),\end{array}$$(A.5)

which is a function of normalisation temperature only. Therefore, for a given reference (dimensional) temperature T₀, the initial neutral ξ_n0 and plasma ξ_p0 fractions can be uniquely determined.

We normalise the ionisation and recombination rates by the recombination rate using the normalisation parameters (also in units of s⁻¹):

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{rec},0}={\xi }_{\mathrm{p}0}\frac{n_0}{\sqrt{T_{e0}}}2.6\times {10}^{-19}.\\ \hfill \end{array}$$(A.6)

The normalised form of the rates is then given by:

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{rec}}=\frac{{\mathrm{\Gamma }}_{\mathrm{rec},\dim }}{{\mathrm{\Gamma }}_{\mathrm{rec},0}}\tau =\frac{{\rho }_{\mathrm{p}}}{\sqrt{T_{\mathrm{p}}}}\frac{\sqrt{T_f}}{{\xi }_{\mathrm{p}0}}{\tau }_{\mathrm{IR}},\hfill \end{array}$$(A.7)

$$\begin{array}{cc}& {\mathrm{\Gamma }}_{\mathrm{ion}}=\frac{{\mathrm{\Gamma }}_{\mathrm{ion},\dim }}{{\mathrm{\Gamma }}_{\mathrm{rec},0}}\tau ={\rho }_{\mathrm{p}}\frac{{\mathrm{e}}^{-\chi }{\chi }^{0.39}}{0.232+\chi }\frac{\widehat{R}}{{\xi }_{\mathrm{p}0}}{\tau }_{\mathrm{IR}},\hfill \end{array}$$(A.8)

$$\begin{array}{cc}& \widehat{R}=\frac{2.91\times {10}^{-14}}{2.6\times {10}^{-19}}\sqrt{T_{\mathrm{e}0}},\hfill \end{array}$$(A.9)

where τ = τ_IR/(α(0)ρ_t) is the collisional time multiplied by a factor τ_IR that governs the rate of ionisation and recombination relative to the collisional time. These non-dimensional ionisation and recombination rates are independent of the reference density n₀ and therefore the reference temperature T₀ alone determines the initial rates and equilibrium.

Appendix B: Shock jump derivation

The IRIP equations studied in this paper involve a non-conservative energy equation, where the ionisation potential and heating terms provide energy losses and gains from the system. As such, analysing the jump conditions across the shock requires a slightly different set of equations since the upstream energy is not necessarily equal to the downstream energy. We note that momentum and mass across the shock are conserved so these equations have the same jump equations as the MHD case. For the IRIP case we replace the energy equation with Eq. (22) for the balance of the ionisation potential and heating terms in a slightly different form such that the bulk density ρ = ρ_p + ρ_n is used. Choosing a point sufficiently upstream and downstream of the system such that the system is in thermal equilibrium (T = T_n = T_p), ionisation equilibrium (Γ_recρ_p = Γ_ionρ_n), with zero drift velocity (v = v_n = v_p), and a ionisation potential balance (ϕ_I = A_heat), the jump equations in the deHoffman-Teller frame for the IRIP case are as follows:

$$\begin{array}{cc}& {\left[\rho v_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(B.1)

$$\begin{array}{cc}& {[\rho v_x^2+P+\frac{B_y^2}{2}]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(B.2)

$$\begin{array}{cc}& {[\rho v_x{v}_y-B_x{B}_y]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(B.3)

$$\begin{array}{cc}& {\left[B_x\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(B.4)

$$\begin{array}{cc}& {[v_x{B}_y-v_y{B}_x]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(B.5)

$$\begin{array}{cc}& {\left[F(T){\rho }^2{\xi }_i^2\right]}_{\mathrm{d}}^{\mathrm{u}}=0,\hfill \end{array}$$(B.6)

for upstream (superscript u) and downstream (superscript d) states.

It can easily be seen from Eq. (B.4) that B_x is constant across the shock, that is,

$$\begin{array}{c}\hfill \frac{B_x^{\mathrm{d}}}{B_x^{\mathrm{u}}}=1.\end{array}$$(B.7)

By defining a compressional ratio r as

$$\begin{array}{c}\hfill \frac{{\rho }^{\mathrm{d}}}{{\rho }^{\mathrm{u}}}=r,\end{array}$$(B.8)

the mass conservation (Eq. (B.1)) gives

$$\begin{array}{c}\hfill \frac{v_x^{\mathrm{d}}}{v_x^{\mathrm{u}}}=\frac{1}{r}·\end{array}$$(B.9)

The xy-momentum (Eq. (B.3)) can be written as

$$\begin{array}{cc}& {\rho }^{\mathrm{d}}{v}_x^{\mathrm{d}}{v}_y^{\mathrm{d}}-B_x{B}_y^{\mathrm{d}}={\rho }^{\mathrm{u}}{v}_x^{\mathrm{u}}{v}_y^{\mathrm{u}}-B_x{B}_y^{\mathrm{u}},\hfill \end{array}$$(B.10)

$$\begin{array}{cc}& \frac{{\rho }^{\mathrm{d}}}{{\rho }^{\mathrm{u}}}{v}_x^{\mathrm{d}}{v}_y^{\mathrm{d}}-\frac{B_x{B}_y^{\mathrm{d}}}{{\rho }^{\mathrm{u}}}=v_x^{\mathrm{u}}{v}_y^{\mathrm{u}}-\frac{B_x{B}_y^{\mathrm{u}}}{{\rho }^{\mathrm{u}}},\hfill \end{array}$$(B.11)

$$\begin{array}{cc}& r{v}_x^{\mathrm{d}}{v}_y^{\mathrm{d}}-\frac{B_x}{{\rho }^{\mathrm{u}}}{B}_y^{\mathrm{d}}={\rho }^{\mathrm{u}}{v}_x^{\mathrm{u}}{v}_y^{\mathrm{u}}-\frac{B_x}{{\rho }^{\mathrm{u}}}{B}_y^{\mathrm{u}},\hfill \end{array}$$(B.12)

$$\begin{array}{cc}& v_x^{\mathrm{u}}{v}_y^{\mathrm{d}}-\frac{B_x}{{\rho }^{\mathrm{u}}}{B}_y^{\mathrm{d}}={\rho }^{\mathrm{u}}{v}_x^{\mathrm{u}}{v}_y^{\mathrm{u}}-\frac{B_x}{{\rho }^{\mathrm{u}}}{B}_y^{\mathrm{u}}.\hfill \end{array}$$(B.13)

Now, let B_y = tan(θ) and define the Alfvén Mach number as $V_{\mathit{Ax}}^2=\frac{B_x^2}{{\rho }^{\mathrm{u}}}$ to give

$$\begin{array}{cc}& v_x^{\mathrm{u}}{v}_y^{\mathrm{d}}-V_{\mathit{Ax}}^{2}tan({\theta }^{\mathrm{d}})=v_x^{\mathrm{u}}{v}_y^{\mathrm{u}}-V_{\mathit{Ax}}^{2}tan({\theta }^{\mathrm{u}}),\hfill \end{array}$$(B.14)

$$\begin{array}{cc}& v_x^{\mathrm{u}}(v_y^{\mathrm{d}}-v_y^{\mathrm{u}})=V_{\mathit{Ax}}^2(tan({\theta }^{\mathrm{d}})-tan({\theta }^{\mathrm{u}})).\hfill \end{array}$$(B.15)

In the Hoffman-Teller frame, the electric field is zero, hence v_y = v_xtan(θ^u). Additionally, v_x/v_y = B_x/B_y = 1/tan(θ).

$$\begin{array}{cc}& v_x^{\mathrm{u}}(v_x^{\mathrm{d}}tan({\theta }^{\mathrm{d}})-v_x^{\mathrm{u}}tan({\theta }^{\mathrm{u}}))=V_{\mathit{Ax}}^2(tan({\theta }^{\mathrm{d}})-tan({\theta }^{\mathrm{u}})),\hfill \end{array}$$(B.16)

$$\begin{array}{cc}& v_x^{\mathrm{u}}(v_x^{\mathrm{u}}tan({\theta }^{\mathrm{d}})/r-v_x^{\mathrm{u}}tan({\theta }^{\mathrm{u}}))=V_{\mathit{Ax}}^2(tan({\theta }^{\mathrm{d}})-tan({\theta }^{\mathrm{u}})),\hfill \end{array}$$(B.17)

$$\begin{array}{cc}& v_x^{2u}(tan({\theta }^{\mathrm{d}})/r-tan({\theta }^{\mathrm{u}}))=V_{\mathit{Ax}}^2(tan({\theta }^{\mathrm{d}})-tan({\theta }^{\mathrm{u}})),\hfill \end{array}$$(B.18)

$$\begin{array}{cc}& tan({\theta }^{\mathrm{d}})(\frac{1}{r}-\frac{V_{\mathit{Ax}}^2}{v_x^{2u}})=tan({\theta }^{\mathrm{u}})(1-\frac{V_{\mathit{Ax}}^2}{v_x^{2u}}),\hfill \end{array}$$(B.19)

$$\begin{array}{cc}& \frac{tan({\theta }^{\mathrm{d}})}{tan({\theta }^{\mathrm{u}})}=\frac{B_y^{\mathrm{d}}}{B_y^{\mathrm{u}}}=r\frac{A_x^{u2}-1}{A_x^{u2}-r}.\hfill \end{array}$$(B.20)

From electric field requirements: In the deHoffmann-Teller frame, the electric field is zero on either side of the shock:

$$\begin{array}{cc}& B_x=\frac{v_x^{\mathrm{d}}{B}_y^{\mathrm{d}}}{v_y^{\mathrm{d}}}=\frac{v_x^{\mathrm{u}}{B}_y^{\mathrm{u}}}{v_y^{\mathrm{u}}},\hfill \end{array}$$(B.21)

$$\begin{array}{cc}& \frac{v_y^{\mathrm{d}}}{v_y^{\mathrm{u}}}=\frac{v_x^{\mathrm{d}}}{v_x^{\mathrm{u}}}\frac{B_y^{\mathrm{d}}}{B_y^{\mathrm{u}}},\hfill \end{array}$$(B.22)

$$\begin{array}{cc}& \frac{v_y^{\mathrm{d}}}{v_y^{\mathrm{u}}}=\frac{A_x^{u2}-1}{A_x^{u2}-r}.\hfill \end{array}$$(B.23)

From the xx-momentum equation:

$$\begin{array}{cc}& {\rho }^{\mathrm{d}}{v}_x^{d2}+P^{\mathrm{d}}+\frac{B_y^{d2}}{2}={\rho }^{\mathrm{d}}{v}_x^{u2}+P^{\mathrm{u}}+\frac{B_y^{u2}}{2},\hfill \end{array}$$(B.24)

$$\begin{array}{cc}& \frac{P^{\mathrm{d}}}{P^{\mathrm{u}}}=1+\frac{1}{P^{\mathrm{u}}}({\rho }^{\mathrm{u}}{v}_x^{u2}-{\rho }^{\mathrm{d}}{v}_x^{d2}+\frac{B_y^{u2}}{2}-\frac{B_y^{d2}}{2})\hfill \end{array}$$(B.25)

$$\begin{array}{cc}& =1+\frac{1}{P^{\mathrm{u}}}({\rho }^{\mathrm{u}}(v_x^{u2}-\frac{1}{r}{v}_x^{u2})\hfill \\ \hfill & +\frac{B_y^{2u}}{2}(1-r^2{\left(\frac{A_x^{2u}-1}{A_x^{2u}-r}\right)}^2)).\hfill \end{array}$$(B.26)

The plasma-β can be defined as

$$\begin{array}{cc}\hfill \beta & =\frac{P}{B^2/2}=\frac{2P}{B_x^2+B_y^2}=\frac{2P}{B_x^2(1+{tan}^2(\theta ))},\hfill \end{array}$$(B.27)

$$\begin{array}{cc}\hfill P& =\frac{B_x^2\beta (1+tan(\theta ))}{2}.\hfill \end{array}$$(B.28)

Therefore,

$$\begin{array}{cc}\hfill \frac{P^{\mathrm{d}}}{P^{\mathrm{u}}}& =1+\frac{2}{\beta (1+{tan}^2(\theta ))}\frac{A_x^{u2}}{r}[r-1+\frac{r{tan}^2(\theta )}{2A_x^{u2}}\times \hfill \end{array}$$(B.29)

$$\begin{array}{cc}& (1-r^2{\left(\frac{A_x^{2u}-1}{A_x^{2u}-r}\right)}^2)].\hfill \end{array}$$(B.30)