Intermediate shock substructures within a slow-mode shock occurring in partially ionised plasma

TL;DR

This study uses high-resolution simulations to reveal the formation and persistence of intermediate shocks within slow-mode shocks in partially ionised plasmas, relevant to solar and astrophysical phenomena.

Contribution

It demonstrates the formation of long-lived intermediate (Alfven) shocks within slow-mode shocks in a partially ionised plasma using a two-fluid numerical approach.

Findings

Intermediate shocks form within slow-mode shocks in partially ionised plasma.

Neutral fluid heating causes a Sedov-Taylor-like expansion and velocity overshoots.

The intermediate shock persists across various physical parameters and initial conditions.

Abstract

Slow-mode shocks are important in understanding fast magnetic reconnection, jet formation and heating in the solar atmosphere, and other astrophysical systems. The atmospheric conditions in the solar chromosphere allow both ionised and neutral particles to exist and interact. Under such conditions, fine substructures exist within slow-mode shocks due to the decoupling and recoupling of the plasma and neutral species. We study numerically the fine substructure within slow-mode shocks in a partially ionised plasma, in particular, analysing the formation of an intermediate transition within the slow-mode shock. High-resolution 1D numerical simulations are performed using the (P\underline{I}P) code using a two-fluid approach. We discover that long-lived intermediate (Alfven) shocks can form within the slow-mode shock, where there is a shock transition from above to below the Alfven speed and a reversal of the magnetic field across the shock front. The collisional coupling provides frictional heating to the neutral fluid, resulting in a Sedov-Taylor-like expansion with overshoots in the neutral velocity and neutral density. The increase in density results in a decrease of the Alfven speed and with this the plasma inflow is accelerated to above the Alfven speed within the finite width of the shock leading to the intermediate transition. This process occurs for a wide range of physical parameters and an intermediate shock is present for all investigated values of plasma-$\beta$, neutral fraction, and magnetic angle. As time advances the magnitude of the magnetic field reversal decreases since the neutral pressure cannot balance the Lorentz force. The intermediate shock is long-lived enough to be considered a physical structure, independent of the initial conditions.

Equations

$${\rho }^{\text{u}}{v}_{\perp }^{\text{u}}={\rho }^{\text{d}}{v}_{\perp }^{\text{d}},$$

$${\rho }^{\text{u}}{v}_{\perp }^{\text{u}}{v}_{\parallel }^{\text{u}}-\frac{1}{{\mu }_0}{B}_{\perp }^{\text{u}}{B}_{\parallel }^{\text{u}}={\rho }^{\text{d}}{v}_{\perp }^{\text{d}}{v}_{\parallel }^{\text{d}}-\frac{1}{{\mu }_0}{B}_{\perp }^{\text{d}}{B}_{\parallel }^{\text{d}},$$

$${\rho }^{\text{u}}{v}_{\perp }^{\text{u}}{v}_{\perp }^{\text{u}}+P^{\text{u}}+\frac{B^{\text{u}2}}{2{\mu }_0}={\rho }^{\text{d}}{v}_{\perp }^{\text{d}}{v}_{\perp }^{\text{d}}+P^{\text{d}}+\frac{B^{\text{d}2}}{2{\mu }_0},$$

$$\frac{\gamma }{\gamma -1}\frac{P^{\text{u}}}{{\rho }^{\text{u}}}+\frac{v^{\text{u}2}}{2}=\frac{\gamma }{\gamma -1}\frac{P^{\text{d}}}{{\rho }^{\text{d}}}+\frac{v^{\text{d}2}}{2},$$

$$v_{\perp }^{\text{u}}{B}_{\parallel }^{\text{u}}-v_{\parallel }^{\text{u}}{B}_{\perp }^{\text{u}}=v_{\perp }^{\text{d}}{B}_{\parallel }^{\text{d}}-v_{\parallel }^{\text{d}}{B}_{\perp }^{\text{d}}$$

$$B_{\perp }^{\text{u}}=B_{\perp }^{\text{d}}.$$

$$A_{\perp }^{\text{u}2}=[A_{\perp }^{\text{d}2}(\frac{\gamma -1}{\gamma }\left(\frac{\gamma +1}{\gamma -1}-{\tan }^2{\theta }^{\text{u}}\right){\left(A_{\perp }^{\text{d}2}-1\right)}^2$$

$$+{\tan }^2{\theta }^{\text{u}}\left(\frac{\gamma -1}{\gamma }{A}_{\perp }^{\text{d}2}-1\right)\left(A_{\perp }^{\text{d}2}-2\right))-\frac{{\beta }^{\text{u}}}{{\cos }^2{\theta }^{\text{u}}}{\left(A_{\perp }^{\text{d}2}-1\right)}^2]$$

$$/\left[\frac{\gamma -1}{\gamma }\frac{{\left(A_{\perp }^{\text{d}2}-1\right)}^2}{{\cos }^2{\theta }^{\text{u}}}-A_{\perp }^{\text{d}2}{\tan }^2{\theta }^{\text{u}}\left(\frac{\gamma -1}{\gamma }{A}_{\perp }^{\text{d}2}-1\right)\right].$$

$$\frac{\partial {\rho }_{\text{n}}}{\partial t}+\nabla \cdot ({\rho }_{\text{n}}{\text{v}}_{\text{n}})=0,$$

$$\frac{\partial }{\partial t}({\rho }_{\text{n}}{\text{v}}_{\text{n}})+\nabla \cdot ({\rho }_{\text{n}}{\text{v}}_{\text{n}}{\text{v}}_{\text{n}}+P_{\text{n}}\text{I})$$

$$=-{\alpha }_c{\rho }_{\text{n}}{\rho }_{\text{p}}({\text{v}}_{\text{n}}-{\text{v}}_{\text{p}}),$$

$$\frac{\partial e_{\text{n}}}{\partial t}+\nabla \cdot \left[{\text{v}}_{\text{n}}(e_{\text{n}}+P_{\text{n}})\right]$$

$$=-{\alpha }_c{\rho }_{\text{n}}{\rho }_{\text{p}}\left[\frac{1}{2}({\text{v}}_{\text{n}}^2-{\text{v}}_{\text{p}}^2)+3\left(\frac{P_n}{{\rho }_n}-\frac{1}{2}\frac{P_p}{{\rho }_p}\right)\right],$$

$$e_{\text{n}}=\frac{P_{\text{n}}}{\gamma -1}+\frac{1}{2}{\rho }_{\text{n}}{v}_{\text{n}}^2,$$

$$\frac{\partial {\rho }_{\text{p}}}{\partial t}+\nabla \cdot ({\rho }_{\text{p}}{\text{v}}_{\text{p}})=0$$

$$\frac{\partial }{\partial t}({\rho }_{\text{p}}{\text{v}}_{\text{p}})+\nabla \cdot \left({\rho }_{\text{p}}{\text{v}}_{\text{p}}{\text{v}}_{\text{p}}+P_{\text{p}}\text{I}-\text{B B}+\frac{{\text{B}}^2}{2}\text{I}\right)$$

$$={\alpha }_c{\rho }_{\text{n}}{\rho }_{\text{p}}({\text{v}}_{\text{n}}-{\text{v}}_{\text{p}}),$$

$$\frac{\partial }{\partial t}\left(e_{\text{p}}+\frac{{\text{B}}^2}{2}\right)+\nabla \cdot \left[{\text{v}}_{\text{p}}(e_{\text{p}}+P_{\text{p}})-({\text{v}}_p\times \text{B})\times \text{B}\right]$$

$$={\alpha }_c{\rho }_{\text{n}}{\rho }_{\text{p}}\left[\frac{1}{2}({\text{v}}_{\text{n}}^2-{\text{v}}_{\text{p}}^2)+3\left(\frac{P_n}{{\rho }_n}-\frac{1}{2}\frac{P_p}{{\rho }_p}\right)\right],$$

$$\frac{\partial \text{B}}{\partial t}-\nabla \times ({\text{v}}_{\text{p}}\times \text{B})=0,$$

$$e_{\text{p}}=\frac{P_{\text{p}}}{\gamma -1}+\frac{1}{2}{\rho }_{\text{p}}{v}_{\text{p}}^2,$$

$$\nabla \cdot \text{B}=0,$$

$${\alpha }_c={\alpha }_0\sqrt{\frac{T_n+T_p}{2T_0}}$$

$$B_x=0.1$$

$$B_y=-1.0(x>0),1.0(x<0)$$

$${\rho }_{\text{n}}={\xi }_{\text{n}}{\rho }_{\text{t}}$$

$${\rho }_{\text{p}}={\xi }_i{\rho }_{\text{t}}\text{=}(1-{\xi }_{\text{n}}){\rho }_{\text{t}}$$

$$P_{\text{n}}=\frac{{\xi }_{\text{n}}}{{\xi }_{\text{n}}+2{\xi }_i}{P}_{\text{t}}=\frac{{\xi }_{\text{n}}}{{\xi }_{\text{n}}+2{\xi }_i}\beta \frac{B_0^2}{2}$$

$$P_{\text{p}}=\frac{2{\xi }_i}{{\xi }_{\text{n}}+2{\xi }_i}{P}_{\text{t}}=\frac{2{\xi }_i}{{\xi }_{\text{n}}+2{\xi }_i}\beta \frac{B_0^2}{2}$$

$${\rho }_{\text{n}}{v}_{\perp \text{n}}=\text{mbox}const,$$

$${\rho }_{\text{n}}{v}_{\perp \text{n}}{v}_{\parallel \text{n}}=-I_1+\text{mbox}const,$$

$${\rho }_{\text{n}}{v}_{\perp \text{n}}{v}_{\perp \text{n}}+P_{\text{n}}=-I_2+\text{mbox}const,$$

$$v_{\perp \text{n}}\left(\frac{\gamma }{\gamma -1}{P}_{\text{n}}+\frac{1}{2}{\rho }_{\text{n}}{v}_{\text{n}}^2\right)=-I_3+\text{mbox}const,$$

$${\rho }_{\text{p}}{v}_{\perp \text{p}}=\text{mbox}const,$$

$${\rho }_{\text{p}}{v}_{\perp \text{p}}{v}_{\parallel \text{p}}-\frac{1}{{\mu }_0}{B}_{\perp }{B}_{\parallel }=I_1+\text{mbox}const,$$

$${\rho }_{\text{p}}{v}_{\perp \text{p}}{v}_{\perp \text{p}}+P_{\text{p}}+\frac{B^2}{2{\mu }_0}=I_2+\text{mbox}const,$$

$$v_{\perp \text{p}}\left(\frac{\gamma }{\gamma -1}{P}_{\text{p}}+\frac{1}{2}{\rho }_{\text{p}}{v}_{\text{p}}^2\right)=I_3+\text{mbox}const,$$

$$v_{\perp \text{p}}{B}_{\parallel }-v_{\parallel \text{p}}{B}_{\perp }=0,$$

$$B_{\perp }=\text{mbox}const,$$

$$I_1=\int {\alpha }_c(T_{\text{n}},T_{\text{p}}){\rho }_{\text{n}}{\rho }_{\text{p}}(v_{\parallel \text{n}}-v_{\parallel \text{p}})\mathrm{d}\perp ,$$

$$I_2=\int {\alpha }_c(T_{\text{n}},T_{\text{p}}){\rho }_{\text{n}}{\rho }_{\text{p}}(v_{\perp \text{n}}-v_{\perp \text{p}})\mathrm{d}\perp ,$$

$$I_3=\int {\alpha }_c(T_{\text{n}},T_{\text{p}}){\rho }_{\text{n}}{\rho }_{\text{p}}$$

$$\times \left[\frac{1}{2}(v_{\text{n}}^2-v_{\text{p}}^2)+3R_g(T_{\text{n}}-T_{\text{p}})\right]\mathrm{d}\perp .$$

$${\rho }_{\text{n}}{v}_{\perp \text{n}}+{\rho }_{\text{p}}{v}_{\perp \text{p}}=\text{mbox}const,$$

$${\rho }_{\text{n}}{v}_{\perp \text{n}}{v}_{\parallel \text{n}}+{\rho }_{\text{p}}{v}_{\perp \text{p}}{v}_{\parallel \text{p}}-\frac{1}{{\mu }_0}{B}_{\perp }{B}_{\parallel }=\text{mbox}const,$$

$${\rho }_{\text{n}}{v}_{\perp \text{n}}{v}_{\perp \text{n}}+P_{\text{n}}+{\rho }_{\text{p}}{v}_{\perp \text{p}}{v}_{\perp \text{p}}+P_{\text{p}}+\frac{B^2}{2{\mu }_0}=\text{mbox}const,$$

$$v_{\perp \text{n}}\left(\frac{\gamma }{\gamma -1}{P}_{\text{n}}+\frac{1}{2}{\rho }_{\text{n}}{v}_{\text{n}}^2\right)+v_{\perp \text{p}}\left(\frac{\gamma }{\gamma -1}{P}_{\text{p}}+\frac{1}{2}{\rho }_{\text{p}}{v}_{\text{p}}^2\right)$$

$$=\text{mbox}const,$$

$$v_{\perp \text{p}}{B}_{\parallel }-v_{\parallel \text{p}}{B}_{\perp }=0,$$

$$B_{\perp }=\text{mbox}const.$$

$${\rho }_{\text{n}}$$

Full text

*Aims. *We study numerically the fine substructure within slow-mode shocks in a partially ionised plasma, in particular, analysing the formation of an intermediate transition within the slow-mode shock.

*Methods. *High-resolution 1D numerical simulations are performed using the (PIP) code using a two-fluid approach.

*Results. *We discover that long-lived intermediate (Alfvén) shocks can form within the slow-mode shock, where there is a shock transition from above to below the Alfvén speed and a reversal of the magnetic field across the shock front. The collisional coupling provides frictional heating to the neutral fluid, resulting in a Sedov-Taylor-like expansion with overshoots in the neutral velocity and neutral density. The increase in density results in a decrease of the Alfvén speed and with this the plasma inflow is accelerated to above the Alfvén speed within the finite width of the shock leading to the intermediate transition. This process occurs for a wide range of physical parameters and an intermediate shock is present for all investigated values of plasma-$\beta$, neutral fraction, and magnetic angle. As time advances the magnitude of the magnetic field reversal decreases since the neutral pressure cannot balance the Lorentz force. The intermediate shock is long-lived enough to be considered a physical structure, independent of the initial conditions.

*Conclusions. *Intermediate shocks are a physical feature that can exist as shock substructure for long periods of time in partially ionised plasma due to collisional coupling between species.

1 Introduction

Shocks occur readily in the lower solar atmosphere, driven by wave steepening, e.g., umbral flashes (Beckers & Tallant, 1969; Houston et al., 2018), or magnetic reconnection driven events, such as Ellerman bombs (Ellerman, 1917; Nelson et al., 2013). The listed phenomena occur in regions of the sun where partially ionised effects are thought to play a key role in the underlying physics. As such, to understand phenomena occurring in the lower solar atmosphere, we must also understand the role of partial ionisation in shocks.

In magnetohydrodynamics (MHD), there are three characteristic wave speeds (slow, Alfvén and fast) which leads to a multitude of potential shock transitions (see for example Delmont & Keppens, 2011). The slow-mode shock (transition from super-slow to sub-slow flow speeds) is of particular importance due to its role in magnetic reconnection, where a change in connectivity of the magnetic field results in a release of stored magnetic energy. The energy can be released instantaneously via Joule heating, or post-reconnection by the influence of the magnetic field on the plasma. The classical 2D schematic for fast reconnection features slow-mode shocks (Petschek, 1964). Recent work suggests more complicated reconnection configurations occur in the solar atmosphere, e.g., plasmoid instability, 3D topology, extended MHD (see reviews by Yamada et al., 2010; Pontin, 2011; Loureiro & Uzdensky, 2016; Cassak et al., 2017). In the Petschek model, the presence of slow-mode shocks in the outflow region allow efficient transport of energy away from the reconnection region, increasing the reconnection rate. Slow-mode shocks are also found to form in alternative models for fast reconnection (Liu et al., 2012; Innocenti et al., 2015; Shibayama et al., 2015).

The intermediate shock (transition from above to below the Alfvén velocity) has been shown to exist in single-fluid resistive MHD representative of the Earth’s magnetopause (Karimabadi, 1995). The physical mechanism being that the non-ideal region around the shock allows a separation of the magnetic field and fluid, and the formation of an intermediate shock between the two ideal regions. The stability of intermediate shocks was proved by Wu & Hada (1991) who studied numerically the formation of these shocks due to wave steepening. The relationship between intermediate shocks and slow-mode shocks was shown analytically by Hau & Sonnerup (1989) for resistive MHD shocks, finding a relation between the downstream slow-mode speed and the strength of the intermediate transition.

In partially ionised plasma, such as in the solar chromosphere or prominances, shocks become more complex. Hillier et al. (2016) studied the formation of slow-mode shocks in partially ionised plasmas. Around the shock front, the plasma and neutral species decouple and recouple resulting in a finite width slow-mode shock. The Lorentz force indirectly affects neutrals through collisions and hence the drift between ion and neutral species becomes an important parameter in partially-ionised plasma. However, observing the ion-neutral drift has been found to be difficult due to line-of-sight effects (Anan et al., 2017). A review of partially-ionised modelling can be found in Khomenko (2017).

In this paper, we utilise high-resolution two-fluid numerical simulations to investigate the substructure within slow-mode shocks. We discover that long-lived intermediate shocks can form within the finite width of the shock as a result of the fluid coupling and decoupling around the shock front. This result has applications to the solar chromosphere, as well as interplanetary and interstellar partially ionised plasma. We study the formation of compound intermediate shocks in partially ionised plasma, and the conditions in which such shock structures can form.

The outline of this paper is as follows. First we define the shock transitions and classifications to be used in this paper, and the analytical solution to the magnetohydrodynamic (MHD) and partially ionised plasma (PIP) equations in the Hoffman-Teller shock frame. Next, the numerical methods and initial conditions under consideration are introduced. Following this, a reference case is presented where we compare an MHD and a PIP simulation, and identify the key differences, investigating the physical evolution that produces the intermediate shock. Finally, we consider a parameter study and investigate the effect of the plasma properties on the lifetime and magnitude of the shock. We find that for our initial conditions, an intermediate shock will always form in a partially-ionised plasma. However the lifetime and magnitude depend heavily on the plasma properties.

2 Methodology

2.1 Shock classifications

MHD waves have three characteristic speeds: Alfvén ($V_{A}$), slow ($V_{s}$), and fast ($V_{f}$). As such, multiple shock transitions are possible, depending on the magnitude of the velocity, relative to the characteristic speeds in the pre- and post-shock regions.

Following the approach of Delmont & Keppens (2011), we classify shock transitions using the relationship between the normal flow velocity $v_{\perp}$ and the characteristic speeds:

• •

(1) superfast: $V_{f}<\left\lvert v_{\perp}\right\rvert$,

• •

(2) subfast: $V_{A}<\left\lvert v_{\perp}\right\rvert<V_{f}$,

• •

(3) superslow: $V_{s}<\left\lvert v_{\perp}\right\rvert<V_{A}$,

• •

(4) subslow: $0<\left\lvert v_{\perp}\right\rvert<V_{s}$,

• •

($\infty$) static: $v_{\perp}=0$.

Defining the upstream condition u and downstream condition d, several shocks of the form $\text{u}\rightarrow\text{d}$ are possible:

• •

$1\rightarrow 2$ fast shocks,

• •

$3\rightarrow 4$ slow shocks,

• •

$1\rightarrow 2=3$ switch-on,

• •

$2=3\rightarrow 4$ switch-off,

• •

$1\rightarrow 3,1\rightarrow 4,2\rightarrow 3,2\rightarrow 4$ intermediate shocks.

2.2 Hoffman-Teller Equations

The Hoffman-Teller frame allows for jump relations to be derived from the MHD equations. In this choice of rest frame, the velocity and magnetic field vectors are in the same plane either side of the shock, i.e., the electric field across the shock is zero. This reduces the MHD equations to a 2-dimensional problem for variables perpendicular ($\perp$) and parallel ($\parallel$) to the shock front.

2.2.1 MHD Solution

In the Hoffmann-Teller frame, the MHD equations can be integrated and the upstream (${}^{\text{u}}$) and downstream (${}^{\text{d}}$) conditions can be equated:

[TABLE]

We can also define the following: Alfvén Mach number as $A_{\perp}^{2}=v_{\perp}^{2}\frac{\mu_{0}\rho}{B_{x}^{2}}$, plasma-$\beta$ as $\beta=\frac{2\mu_{0}P}{B^{2}}$, and the angle between the velocity and magnetic field $\theta$.

A solution to Equations (1-6) can be found in Hau & Sonnerup (1989) relating the upstream Alfvén velocity to the upstream plasma-$\beta$, upstream $\theta$, and downstream Alfvén velocity, i.e.,

[TABLE]

This is shown graphically in Figure 1 for a choice of plasma-$\beta^{\text{u}}=0.1$ and $\theta^{\text{u}}=\pi/4$. The trivial solution is that the upstream and downstream velocities are identical, i.e., no shock transition. The shaded region shows impossible solutions, where the velocity is higher downstream than upstream. The non-trivial solutions show the possible shock transitions for given plasma-$\beta$ and $\theta$ values. This curve intersects the $A^{\text{u}}=A^{\text{d}}$ line at three points. The intersect labelled $s^{\text{u}}$ denotes the upstream slow speed. At the point $A_{\perp}^{\text{u}}=A_{\perp}^{\text{d}}=1$, a rotational discontinuity is possible, where there is a change in the angle of the magnetic field but plasma properties remain the same. Slow-mode transitions (blue line) are bounded by $A_{\perp}^{\text{u}}=1$. The downstream slow speed (point labelled on Figure 1 as $s^{\text{d}}$) corresponds to the critical value that separates strong ($2\rightarrow 4$) and weak ($2\rightarrow 3$) intermediate transitions. The strong transition ($2\rightarrow 4$) is closely linked to the slow-mode transition ($3\rightarrow 4$) through its relation to the downstream slow speed.

2.2.2 PIP Solution

In the PIP equations, collisional terms are taken into account and therefore, the equations have a non-zero right-hand side, see Equations (8-17). As such, when the equations are integrated in the Hoffmann-Teller frame, there are integral terms on the right hand side, see Equations (25-37), as opposed to constants in the MHD Equations. However, by adding the neutral and plasma species together, the integral terms can be eliminated, and by choosing a point upstream and downstream such that the drift velocity equals zero, one recovers the MHD equations in the Hoffman-Teller frame, where the Alfvén speed and plasma-$\beta$ values depend on the total plasma density ($\rho_{\text{n}}+\rho_{\text{p}}$), see Appendix A. These equations are also independent from the neutral fraction and hence, in the steady state solution, the values either side of the shock are governed by the MHD solution, when points are chosen such that the drift velocity is zero. However, within the shock, the species decouple and recouple hence in partially-ionised plasmas, it is of interest to study the substructure that occurs within the finite-width of the shock.

2.3 Numerical methods

Two-fluid numerical simulations are performed using the (PIP) code (Hillier et al., 2016) which solves the interactions of a neutral fluid, and a coupled ion electron plasma. The simulations are 1D and use high resolution to resolve the substructure present within the shock. A first order HLLD (Harten-Lax-van Leer-Discontiunities) scheme is used to prevent spurious oscillations from occurring around the shock interface. The normalised PIP equations are given below:

[TABLE]

for neutral (subscript n) and plasma (subscript p) species. The neutral equations (8-11) are independent of the magnetic field. The plasma equations (12-17) are similar to the MHD equations however include collisional terms that couple the plasma to the neutral fluid. The collisional coefficient $\alpha_{c}$ is defined as

[TABLE]

where $T_{0}$ is the normalisation temperature such that the sound speed $c_{s}^{2}=V_{A}^{2}$ and, in the normalised form, $\alpha_{0}=1$. The magnetic field is normalised using $\textbf{B}=\textbf{B}/\sqrt{\mu_{0}}$. Details of the equations and their implementation in the (PIP) code can be found in Hillier et al. (2016).

Our simulation data is translated into the shock frame by calculating the shock propagation speed in the MHD solution ($v_{s}$). The position and velocity variables are remapped to the shock frame, i.e., $v_{x{\text{n}}s,x{\text{p}}s}=v_{x{\text{n}},x{\text{p}}}-v_{s}$, $x_{s}=x/t-v_{s}$.

2.4 Initial conditions

The initial conditions in this paper are an extension of the model for slow-mode shock created in reconnection proposed by Petschek (1964). Initially, a discontinuous magnetic field is specified across the boundary. The plasma and neutral fluids are assumed to be in thermal equilibrium at the start of the simulation.

[TABLE]

where $\xi_{n}$ and $\xi_{i}$ are the neutral and ion fractions respectively.

The $x=0$ boundary is treated as reflective, but such that magnetic field can penetrate the boundary. The system is normalised to have an Alfvén velocity of unity. The collisional time as determined by the bulk fluid density is calculated as $\tau=1/(\alpha(T_{0})\rho_{\text{t}})=1$.

These initial and boundary conditions were used in the work of Hillier et al. (2016) to investigate the slow-mode shock. Here we investigate the formation and lifetimes of a feature that was present in Figures 4a and 5a of Hillier et al. (2016) but not discussed, namely the intermediate shock transition, which features a reversal of magnetic field across the shock interface.

3 Results

In this section, we present an MHD simulation and a PIP simulation and highlight the key differences. A snapshot of the results is shown in Figure 2. Both simulations use 256,000 grid cells to resolve the spatial dimension. These two simulations use the following parameters: $\beta=1.0$, $B_{x}=0.1$, $\xi_{\text{n}}=0.9$. Note that the neutral fraction is only used in the PIP case and the MHD simulation is fully ionised. The effect of different parameters is investigated in Section 4.2.

3.1 Reference MHD solution

In the MHD case (Figure 2, black lines), the initial conditions produce a rarefaction wave that drives fluid at the local Alfvén speed towards a slow-mode shock. The solution is expanding in time, however, is steady-state when plotted on a time-normalised axis ($x/t$), hence can be described as pseudo-steady.

It is important to note that in the MHD solution, the slow-mode shock possesses no complex substructure, has a discontinuous jump (no finite width), and is purely a $3\rightarrow 4$ transition. There is no reversal in magnetic field across the shock (see Figure 2). These features do not hold when the model is extended to include additional effects (e.g. partial ionisation) where substructure can form in the shock front affecting the shock dynamics and resultant heating and energy transport.

3.2 Reference PIP solution

The red and blue lines in Figure 2 show the neutral and ion fluid respectively for the PIP simulation, using the same parameters as the MHD case. The snapshot is taken after 2500 collisional times when the system can be assumed to be in a quasi-self-similar state. In the PIP case, the ion and neutral species decouple and recouple around the shock front, resulting in a finite width slow-mode shock where shock sub-structure can occur (compared to the discontinuous jump in the MHD case).

The initial conditions are a discontinuity in the magnetic field only, therefore directly affect the plasma only. The neutral fluid can only be influenced indirectly by collisions with the ionised plasma. As such, the system tends towards a pseudo-steady-state as time tends towards infinity. Hillier et al. (2016) presented the dynamic evolution for a similar case:

• •

Initialisation The initial conditions are a discontinuity in the magnetic field and hence directly affect the plasma only, and the neutral fluid is affected via collisions. At $\tau=1$ the system is highly decoupled. A Sedov-Taylor-like expansion occurs in the neutral fluid due to frictional heating between the two species.

• •

Weak coupling A fast-mode rarefaction wave forms and drives fluid towards a slow-mode shock. At this time, the rarefaction wave and slow-mode shock front are coupled and interact.

• •

Strong coupling In this phase, the rarefaction wave and slow-mode shock decouple. The slow-mode shock can be considered to be independently evolving.

• •

Quasi-self-similar state System tends towards a self-similar state as time tend towards infinity. This is plotted in Figure 2. The rarefaction wave is separated from the slow-mode shock, Figure 2a. The two fluids are reasonably well coupled except within the slow-mode shock where large drift velocities are present, Figure 2a. The neutral overshoot is present at late times indicating that it is a stable feature of the system (Figure 2a).

A feature present but not discussed in the work of Hillier et al. (2016) is the reversal of the magnetic field across the shock front (Figure 2b). This is a long-lived (but transient) feature of the system that exists long after the rarefaction wave and slow-mode shock have separated. The magnetic field reversal is a signature that an intermediate transition exists within the slow-mode shock.

4 Analysis

4.1 Intermediate Shock

The reversal in magnetic field in Figure 2b is a key indicator of an intermediate shock, where the plasma transitions from above the local Alfvén speed, to below it. An intermediate shock is present in all partially ionised simulations as a transient feature of the system, however is not present in the MHD cases. The feature arises due to the interactions of ion and neutral particles. To analyse this feature, a new frame of reference is chosen such that the shock is stationary, i.e., $x_{s}=x/t-v_{s}$ where $v_{s}$ is the propagation speed of the MHD shock. In this frame, the velocity is also adjusted to account for the shock propagation speed, i.e., $v_{x}-v_{s}$. The shock frame allows us to correctly compute the transitions across the shock front and compare with the analytical results. The evolution of this structure is as follows:

(1) Sedov-Taylor-like expansion of the neutral fluid:

The initial conditions are a discontinuity in the magnetic field that affects the plasma only, resulting in a separation of neutral and plasma species in the shock front (Figure 3a). The large drift velocity creates a rapid heating of the neutral fluid. Subsequently, a hydrodynamic, Sedov-Taylor-like expansion occurs in the neutral fluid, resulting in an overshoot of neutral velocity (Figure 3a), and a sudden increase in neutral density (Figure 3b).

(2) Acceleration of plasma:

The rarefaction wave drives velocity towards the slow-mode shock at the Alfvén speed (Figure 3c). This is also true in the MHD case and is a feature of the shock problem being studied. Within the slow-mode shock front, the bulk Alfvén speed decreases due to the increase in neutral density (Figure 3d) and the reversal in the magnetic field (Figure 2b). The plasma inside the shock is accelerated to above the Alfvén speed (Figure 3c), resulting in an intermediate transition (Figure 3e).

(3) Decay with time:

The intermediate shock is sustained by the neutral pressure (Figure 3f), which, because of the finite coupling, cannot balance the Lorentz force. As time tends to infinity, the neutral pressure equalises and the magnetic field reversal becomes decreases in magnitude, see Figure 4. This is an indicator that the intermediate shock disappears with time. For these parameters, the magnetic field reversal is still prominent after 10000 collisional times, indicating that whilst this is a transient feature of the system, it is sufficiently long-lived to be physical and independent of the initial conditions.

The frictional heating ($\alpha_{c}(T_{\text{n}},T_{\text{p}})\rho_{\text{n}}\rho_{\text{p}}(v_{\text{n}}-v_{\text{p}})^{2}$) and thermal damping ($3\alpha_{c}(T_{\text{n}},T_{\text{p}})\rho_{\text{n}}\rho_{\text{p}}(P_{\text{n}}/\rho_{\text{n}}-P_{\text{p}}/(2\rho_{\text{p}}))$) between the two species are shown in Figure 5 through the finite-width of the shock. Both the frictional heating and thermal damping reach their maximum values at the intermediate transition. Here, there is the largest drift velocity (Figure 3a). Either side of the shock, the species are reasonably well coupled and hence there is minimal frictional heating.

4.2 Parameter Study

4.2.1 Changing ionisation fraction

The simulations in this section use $\beta=0.1,B_{x}=0.3$ with the ionisation fraction $\xi_{i}$ varying. The results can be scaled by a rough estimate of the time scale changes by multiplying by the ionisation fraction, i.e., $t\xi_{i}$. As the neutral fraction increases, the influence of the neutral pressure becomes more important. The initial evolution stages generate an overshoot in the neutral pressure due to the collisional coupling. There is then an interplay between the neutral pressure and the magnetic tension. When $\xi_{\text{n}}\leq 0.9$, the neutral pressure increases with time to balance the magnetic tension. When $\xi_{\text{n}}>0.9$ the neutral pressure overshoot decreases with time. This is shown in Figure 6 with the $\xi_{\text{n}}=0.9$ (black dashed) and the $\xi_{\text{n}}=0.99999$ (magenta dashed) lines. Whilst the two results look very different, the same mechanism is occurring whereby the system seeks an equilibrium; the main difference is the magnitude of the neutral pressure overshoot due to the initial heating. All cases tend towards a constant maximum neutral pressure.

All neutral fractions investigated feature the intermediate transition as a substructure within the slow mode shock. The normalised time scale does not map the behaviour perfectly and there are differences in magnitude and gradients of the magnetic field across the different neutral fractions. The neutral pressure is the key variable due to its role in the equalising the Lorentz force. As the neutral fraction tends to unity, the maximum neutral pressure is reasonably similar across simulations.

4.2.2 Changing Beta

Figure 7 shows the neutral pressure and magnetic field reversal for the low-beta regime. All three cases plotted contain roughly the same gradient in magnetic field over time but have different equilibrium neutral pressures. As the plasma-$\beta$ increases, the maximum equilibrium neutral pressure increases drastically. This results in a larger restoring force and the magnetic field decreases more rapidly.

For large plasma-$\beta$ values, a sonic shock can occur in the neutral fluid, as seen in Hillier et al. (2016). In these simulations, there can be two shocks occurring within the finite-width shock region: intermediate shock in the plasma, and sonic shock in the neutrals.

4.2.3 Changing magnetic angle

The initial $B_{x}$ component of magnetic field is modified to alter the angle of the magnetic field. Note that the initial $B_{y}$ component remains the same and hence the total $B$ strength varies across simulations. Interestingly, the pre-shock region has $|B|\approx 1$ in all simulations as a consequence of the investigated system.

Changing the $B_{x}$ value changes the propagation speeds of waves and hence phenomena occurs on different time scales. From Figure 8 it is clear that a larger $B_{x}$ results in faster equalisation of the neutral pressure and hence the magnetic field reversal and the intermediate shock decay faster, compare to low $B_{x}$ values. Magnetic tension is the equalising force for the intermediate shock. The magnetic tension increases for larger $B_{x}$ values and hence the intermediate shock decays faster.

5 Comments on potential observations

In the MHD case, the fluid and the magnetic field are frozen together, so in the rest frame of the fluid the electric field becomes $\textbf{E}\propto\textbf{v}\times\textbf{B}=0$. However, for the two fluid case, there is a separation of the species within the finite-width of the shock. Within the finite-width, even though there is no electric field felt by the plasma fluid there is therefore an electric field felt by the neutrals. This can be calculated by looking at how fast the neutrals move across the magnetic field, which is the the drift velocity, i.e., $\textbf{E}_{\rm n}\propto(\textbf{v}_{\rm n}-\textbf{v}_{\rm p})\times\textbf{B}$. This metric is zero either side of the shock since the species are fully coupled. Inside the shock, the ions and neutrals decouple and hence there is a localised increase in the neutral electric field inside the shock substructure. The electric field felt by the neutrals is a potential observable for measuring the neutral-ion drift (e.g., Anan et al., 2014).

The maximum electric field felt by the neutrals within the shock is plotted through time in Figure 9 for the reference PIP simulation ($\beta=1$, $B_{x}=0.1$, $\xi=0.9$). As time advances, this tends towards a constant value indicating that there will always be a localised increase in neutral electric filed within the shock. The presence of this enhanced neutral electric field could be a potential observable for the effects of partial ionisation in shock waves.

6 Summary

This paper has demonstrated that intermediate shocks can occur as substructure inside slow-mode shocks due to partial ionisation. High-resolution 1D numerical simulations were performed to fully resolve the substructure that occurs within the finite width of the slow-mode shock. The physical process involved in forming and dissipating the intermediate shock substructure is as follows:

Collisional coupling results in overshoots in neutral velocity and density The initial conditions drive the plasma only. As such, there is a large drift velocity between the two species. The neutral fluids response is to create a Sedov-Taylor-like expansion which creates localised increases in the neutral density and velocity at the interface. 2. 2.

Acceleration of the plasma Our system has velocity driven towards the slow-mode shock at the Alfvén speed. Within the finite-width of the shock, the Alfvén velocity decreases due to the increased neutral pressure and reversal of the magnetic field, and hence the plasma velocity is accelerated inside the shock to above the Alfvén speed. 3. 3.

Decay with time The neutral pressure is not sufficient to balance the Lorentz force and hence there is a gradual evolution whereby the magnetic field reversal tends to zero as time tends towards infinity. The intermediate shock is however sufficiently long-lived to be independent of the initial conditions and considered a physical feature of partially ionised shocks.

The intermediate shock was present as substructure for all tested parameter regimes but the magnitude of the magnetic field reversal was parameter dependent. The larger the equilibrium neutral pressure, the larger the magnetic field reversal. As such, the intermediate shock is strongest for high plasma-$\beta$ values. We would therefore expect this feature to be most significant in the lower atmosphere, e.g., Ellerman bombs.

The work in this paper is analogous to the magnetospheric intermediate shocks discussed in Karimabadi (1995) where the resistivity creates a dispersive region accelerating the plasma towards the shock region and resulting in an intermediate transition. Here we have a similar effect except it is driven by two-fluid interaction in the absence of resistivity.

We have shown that an intermediate shock exists for a wide range of parameters, with varying degrees of magnitude and lifetime. Theoretically, this has the potential to be present across a wide range of phenomena in the solar atmosphere where partial-ionisation effects are important, from wave-steepening events (e.g., umbral flashes), to magnetic reconnection (Ellerman bombs, spicules). Future work will be to analyse the implications and observability of intermediate shocks in the lower solar atmosphere.

In summary, there are four main conclusions in this paper:

Ideal MHD heating across the shock We have shown analytically that the two-fluid equations reduce to the ideal MHD shock equations when the species are coupled either side of the shock (see Appendix A). The consequence of this is that one would expect to obtain ideal MHD-like heating across the shock, with no additional heating from the collisions. 2. 2.

Shocks as substructure in PIP case Within the finite width of a partially-ionied shock, interactions between the two species can lead to the formation of intermediate transitions within the larger shock structure. Intermediate shocks are a transition from above to below the Alfvén speed (here super-Alfvén to sub-slow) and feature a reversal in the magnetic field (see Figure 2). This feature forms due to the collisional effects between the two species inside the finite-width shock. An intermediate transition was present for all tested parameter regimes. 3. 3.

Potential for large currents The formation of an intermediate shock features a sharp reversal in the magnetic field across the intermediate shock front, hence there is the potential for large currents to form inside the large-scale slow-mode transition. Large currents are known to play a role in particle acceleration. Hence the formation of intermediate shocks in partially ionised plasma may lead to an additional particle acceleration mechanism. 4. 4.

Localised electric field experienced by the neutrals Within the finite width of the shock, there is a localised increase in neutral electric field. This tends towards a constant value as time tends to infinity, hence may be a potential observable of ion-neutral interactions within shock fronts.

Acknowledgements

BS and AH are supported by STFC research grant ST/R000891/1. AH is also supported by STFC Ernest Rutherford Fellowship grant number ST/L00397X/2.

Appendix A Hoffman-Teller PIP Equations

In the Hoffman-Teller frame, the 2-fluid PIP equations are:

[TABLE]

The integral terms ($I_{1},I_{2},I_{3}$) can be removed by adding the neutral and ion equations together, i.e.,

[TABLE]

The partial pressure and density can be expressed in terms of a total value using the neutral fraction $\xi_{n}$:

[TABLE]

Substituting these into Equations (38-43) gives:

[TABLE]

Furthermore, we can impose an additional constraint such that either side of the shock the drift velocity equals zero ($v_{\parallel\text{p}}=v_{\parallel\text{n}}=v_{\parallel}$ and $v_{\perp\text{p}}=v_{\perp\text{n}}=v_{\perp}$):

[TABLE]

Equations (54-59) are identical to the MHD equations except here the density is the total density $\rho_{\text{t}}$ and pressure $P_{\text{t}}$. Therefore, the solution to these equations is identical to the Hau & Sonnerup (1989) solution for MHD (Equation 7), independent of the neutral fraction. It should be noted that this is only true over the larger shock structure and inside the shock there is substructure that is highly dependent on the collisional effects.