© ESO 2020

1. Introduction

A compressional wave propagating upwards in the solar atmosphere naturally steepens due to the stratification of the atmosphere and can readily develop non-linearities and shock (e.g. Suematsu et al. 1982). This process occurs throughout the solar atmosphere, for example, it is believed to be the cause of umbral flashes (Beckers & Tallant 1969; Houston et al. 2018; Anan et al. 2019). These waves are magnetohydrodynamic (MHD) in nature and thus depend on the inclination of the magnetic field.

Mangetohydrodynamic (MHD) systems are capable of supporting numerous types of shock transitions due to the three characteristic speeds (slow, Alfvén, fast; e.g. Delmont & Keppens 2011). Slow-mode shocks feature a transition from superslow to subslow and are of particular importance for magnetic reconnection events (Petschek 1964; Innocenti et al. 2015; Shibayama et al. 2015). Fast-mode shocks propagate perpendicular to the magnetic field and are also important for magnetic reconnection (Yokoyama & Shibata 1996). There in also a class of intermediate shocks which have a transition from above to below the Alfvén speed, which can arise due to two-fluid effects in slow-mode shocks (Snow & Hillier 2019).

It is widely known that MHD waves can undergo mode conversion in the solar atmosphere (e.g. Schunker & Cally 2006; Cally & Hansen 2011; Khomenko & Cally 2019), and the same is true of shocks. A recent paper by Pennicott & Cally (2019) studied the propagation of an MHD shock in a stratified, isothermal atmosphere and the separation into slow- and fast-mode shocks above the β_v point, which is defined as the height in the atmosphere where the sound and Alfvén speeds are equal. Their results showed that as a shock passes through the β_v point, it can can undergo mode conversion. In cases where the angle between the magnetic field and the vertical is large, the slow component of the shock formed after mode conversion becomes smoothed.

Many events involving shocks can occur in the lower solar atmosphere where both neutral and ionised species are present and the interactions between these species are thought to play a key role in the underlying dynamics. The lower solar atmosphere is partially ionised, being mostly neutral fluid in the photosphere and chromosphere, and mostly ionised above the transition region. Partial ionisation is known to have several physical consequences in this region including additional heating through dissipation (e.g. Khomenko & Collados 2012) and enhanced reconnection rates (Leake et al. 2013). Studies of two-fluid waves in a stratified atmosphere show that the two-fluid effects lead to additional damping (Popescu Braileanu et al. 2019). Reviews of partial ionisation can be found in Khomenko (2017, 2020), Ballester et al. (2018).

Shocks in partially-ionised plasma (PIP) are of particular interest due to the decoupling and recoupling of species around the shock front (Hillier et al. 2016). As such, shocks are a promising area for the two-fluid effects to have a significant impact on the physics of the system. Previous work has shown that additional shocks can form within the finite width of a reconnection-driven slow-mode shock (Snow & Hillier 2019).

In this paper, we investigate the effect of gravitational stratification and partial ionisation on a propagating shock initiallised through wave steepening. In particular, we investigate the role of two-fluid effects on mode conversion of the shock across the β_v point and the consequences of two-fluid interactions on the structure of the propagating shocks. The collisional coefficient is critical in determining the mode conversion point, the type of shocks present, and the width of these shocks. The fast-mode shock width can exceed the pressure scale height for a broad range of coupling coefficients, which we present as a new potential observable of two-fluid effects in the lower solar atmosphere.

2. Numerical methods

The numerical simulations are performed using the (PIP) code (Hillier et al. 2016) to solve the evolution of a two fluid (neutral and ion+electron) hydrogen plasma. The non-dimensional equations for the evolution of the neutral (subscript n) and charge-neutral plasma (subscript p) used in our study are given by:

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

$$\begin{array}{l}\frac{\partial }{\partial t}({\rho }_{\text{n}}{\mathit{v}}_{\text{n}})+\nabla \cdot ({\rho }_{\text{n}}{\mathit{v}}_{\text{n}}{\mathit{v}}_{\text{n}}+P_{\text{n}}\mathit{I})\\ ={\rho }_{\text{n}}\mathit{g}-{\alpha }_{\text{c}}{\rho }_{\text{n}}{\rho }_{\text{p}}({\mathit{v}}_{\text{n}}-{\mathit{v}}_{\text{p}}),\end{array}$$(2)

$$\begin{array}{l}\frac{\partial e_{\text{n}}}{\partial t}+\nabla \cdot \left[{\mathit{v}}_{\text{n}}(e_{\text{n}}+P_{\text{n}})\right]\\ ={\rho }_{\text{n}}\mathit{g}\cdot {\mathit{v}}_{\text{n}}-\frac{{\alpha }_{\text{c}}{\rho }_{\text{n}}{\rho }_{\text{p}}}{2}\left[({\mathit{v}}_{\text{n}}^2-{\mathit{v}}_{\text{p}}^2)+3\left(\frac{P_{\text{n}}}{{\rho }_{\text{n}}}-\frac{1}{2}\frac{P_{\text{p}}}{{\rho }_{\text{p}}}\right)\right],\end{array}$$(3)

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

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

$$\begin{array}{l}\frac{\partial }{\partial t}({\rho }_{\text{p}}{\mathit{v}}_{\text{p}})+\nabla \cdot \left({\rho }_{\text{p}}{\mathit{v}}_{\text{p}}{\mathit{v}}_{\text{p}}+P_{\text{p}}\mathit{I}-\mathit{B}B+\frac{{\mathit{B}}^2}{2}\mathit{I}\right)\\ ={\rho }_{\text{p}}\mathit{g}+{\alpha }_{\text{c}}{\rho }_{\text{n}}{\rho }_{\text{p}}({\mathit{v}}_{\text{n}}-{\mathit{v}}_{\text{p}}),\end{array}$$(6)

$$\begin{array}{l}\frac{\partial }{\partial t}\left(e_{\text{p}}+\frac{{\mathit{B}}^2}{2}\right)+\nabla \cdot \left[{\mathit{v}}_{\text{p}}(e_{\text{p}}+P_{\text{p}})-({\mathit{v}}_{\text{p}}\times \mathit{B})\times \mathit{B}\right]\\ ={\rho }_{\text{p}}\mathit{g}\cdot {\mathit{v}}_{\text{p}}+\frac{{\alpha }_{\text{c}}{\rho }_{\text{n}}{\rho }_{\text{p}}}{2}\left[({\mathit{v}}_{\text{n}}^2-{\mathit{v}}_{\text{p}}^2)+3\left(\frac{P_{\text{n}}}{{\rho }_{\text{n}}}-\frac{1}{2}\frac{P_{\text{p}}}{{\rho }_{\text{p}}}\right)\right],\end{array}$$(7)

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

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

$$\nabla \cdot \mathit{B}=0,$$(10)

for density ρ, velocity v, pressure P, internal energy e, magnetic field B, the collisional coupling term α_c, and gravity g. We note that the thermal component of the collisional term in Eqs. (3) and (7) is different from in previous papers using the (PIP) code by a factor of 1/2 (Hillier et al. 2016; Snow & Hillier 2019). The corrected term was found to have no significant changes for the reconnection driven shock considered in previous works. The thermal coupling term is now consistent with other two-fluid numerical codes, for example, that of Popescu Braileanu et al. (2019).

This system of equations has been non-dimensionalised in the following way: the velocity v has been non-dimensionalised using the sound speed of the plasma fluid c_s, the density ρ by a reference plasma density ρ_p, 0, and the lengthscale by the characteristic pressure scale height of the plasma Λ_p, 0. This implies that the timescales of the system are nondimensionalised by ${\mathrm{\Lambda }}_{\mathrm{p},0}{c}_{\mathrm{s}}^{-1}$, the pressure P by $c_{\text{s}}^2{\rho }_{\text{p},\text{0}}$ and the magnetic field B by $B_0/\sqrt{4\pi }=c_{\mathrm{s}}\sqrt{{\rho }_{\mathrm{p},0}}$. Here the subscript 0 is used to represent a reference value of a quantity (where these values have been chosen to be those from the x = 0 point) and γ is the adiabatic index. Non-dimenstionalisation using the plasma components was selected as it means that the results of reference MHD simulations would be exactly reproduced in the case of α_c = 0. We assume that both fluids follow the ideal gas law which for the two fluids in non-dimensional form become $T_{\mathrm{n}}=\frac{P_{\mathrm{n}}\gamma }{{\rho }_{\mathrm{n}}}$ and $T_{\mathrm{p}}=\frac{P_{\mathrm{p}}\gamma }{2{\rho }_{\mathrm{p}}}$, respectively.

Here we use the first order HLLD solver in the (PIP) code to accurately capture shocks. The HLLD scheme is described in detail in Miyoshi & Kusano (2005) and can accurately model shocks. This scheme has been used previously with the (PIP) code to model MHD and two-fluid shocks across a range of plasma-β values and ionisation fractions (Hillier et al. 2016; Snow & Hillier 2019). A discussion of the numerical diffusion of the scheme is included in Appendix B. The domain is spatially resolved by 8000 grid cells. The equations are solved for 1.5D evolution in that there is variation in the x direction, but velocities and magnetic fields can exist in the y direction.

2.1. Initial conditions

In this paper, we consider an isothermal atmosphere to study the propagation of two-fluid shocks in an simplified system. The initial density profile for the MHD and PIP simulations is shown in Fig. 1. This initial set up is similar to that of Pennicott & Cally (2019) except here we have the addition of two-fluid effects. Future work will include studying shock behaviour in more realistic atmospheres.

Fig. 1. Density for the initial equilibrium. MHD density is the black-solid line. The PIP simulation has two densities: plasma (blue dashed) and neutral (red dashed). The normalisation leads to the MHD density being the same as the PIP plasma density. The total density for the PIP case is shown in green-dashed. Gravity is constant between the two vertical black lines.

At the simulation boundaries, gravity rolls down to zero. This is purely to stabilise the boundary and simplify the driver formulation. This has no consequence on the propagation of the waves. Gravity is applied in the x direction only (g = [g_x, 0, 0]), where g_x is specified as a piecewise function of height as:

$$\begin{array}{c}\hfill g_x=\{\begin{array}{cc}0\hfill & x\le g_1;\hfill \\ g_0(1+\cos \left(\pi \frac{x-g_2}{g_1-g_2}\right))\hfill & g_1\le x\le g_2;\hfill \\ g_0\hfill & g_2\le x\le g_3;\hfill \\ g_0(1+\cos \left(\pi \frac{x-g_3}{g_4-g_3}\right))\hfill & g_3\le x\le g_4;\hfill \\ 0\hfill & g_4\le x,\hfill \end{array}\end{array}$$

where g₀ = 1/γ from the normalisation. g₁ = −3.7,  g₄ = 3.7 are the upper and lower grid limits, such that gravity is zero in our ghost cells. g₂ = −3.15 and g₃ = 3.15 and between these limits the gravity is constant. The region where gravity is constant is marked on Fig. 1 by the vertical black lines.

Our equilibrium pressure and density profiles are calculated by integrating upwards using the pressure scale heights for the two species, with the density calculated from the ideal gas:

$$\begin{array}{cc}& P_{\mathrm{p};\mathrm{n}}(x)=P_{\mathrm{p};\mathrm{n}}(0)exp{\displaystyle (-{\int }_{g_1}^x\frac{1}{{\mathrm{\Lambda }}_{\mathrm{p};\mathrm{n}}}\mathrm{d}{x}^{\prime }),}\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\mathrm{n}}=\frac{T_{\mathrm{n}}}{\gamma g_x},\hfill \end{array}$$(12)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\mathrm{p}}=\frac{2T_{\mathrm{p}}}{\gamma g_x}·\hfill \end{array}$$(13)

Normalisation is chosen such that the plasma density is unity at x = 0. Initially T_p = T_n = 1/2. All simulations have the same plasma density profile, which results in the PIP simulations having greater total mass, see Fig. 1. We set the magnetic field strength as |B| = 1 such that the β_v point is at x = 0.

The non-dimensional collisional terms are calculated using the average temperature:

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{c}}={\alpha }_0\sqrt{\frac{T_{\mathrm{n}}+T_{\mathrm{p}}}{2}}·\end{array}$$(14)

The collisional coupling coefficient is set to α₀ = 100 unless otherwise stated.

2.2. Boundary conditions

A pulse of one half period is applied to the system at the lower boundary of the form:

$$\begin{array}{cc}& v_{x,\mathrm{p};\mathrm{n}}=Asin(2\pi t/\tau ),\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\rho }_{\mathrm{p};\mathrm{n}}=A\frac{{\rho }_{\mathrm{p};\mathrm{n}}(-3.7)}{c_{\mathrm{s}}}sin(2\pi t/\tau ),\hfill \end{array}$$(16)

$$\begin{array}{cc}& P_{\mathrm{p};\mathrm{n}}=A\frac{\gamma P_{\mathrm{p};\mathrm{n}}(-3.7)}{c_{\mathrm{s}}}sin(2\pi t/\tau ),\hfill \end{array}$$(17)

for period τ, initial density ρ_−3.7, initial pressure P_−3.7 at the base of the domain, and initial sound speed c_s. This form of driver has the advantage of preventing additional modes being driven, compared to exciting velocity only, and is equivalent to the drivers used in Khomenko & Cally (2012) and Felipe (2019). For details on the numerical implementation see Appendix A.

The period is chosen as τ = 0.5 in normalised time units, and the normalised wave amplitude A = 0.1. The choice of amplitude is such that the wave steepens into a shock below the β_v point (where c_s = v_A). The acoustic cutoff frequency can be calculated using ${\omega }_{\mathrm{c}}=\frac{c_{\mathrm{s}}}{2H}\sqrt{1-2\frac{\mathrm{d}H}{\mathrm{d}x}}$, where $H=\frac{c_{\mathrm{s}}^2}{\gamma g}$. Where gravity is constant (−3.15 <  x <  3.15), ω_c = 0.5 which is smaller than the driver frequency of 2π/τ = 4π. Due to the gravity variations near our boundaries the acoustic cutoff varies near the edges of our domain however the maximum value is max(ω_c)≈1.22 which is still smaller than the driver frequency.

In the two-fluid case, the velocity is applied to both the neutral and plasma species. As such, more total kinetic energy is applied to the two-fluid case than the MHD case since the PIP simulation have higher total density, see Fig. 1.

The boundary at the upper end of the domain is taken to be an open boundary. Simulations are stopped before the fast-mode component reaches this boundary.

3. Fundamentals of two-fluid waves

3.1. Characteristic speeds in two-fluid plasma

For two-fluid plasma, one can define the wave speeds based on the total density, or the plasma density alone:

$$\begin{array}{cc}\hfill c_{\mathrm{s}}^2& =\gamma P_{\mathrm{p}}/{\rho }_{\mathrm{p}},\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill c_{\mathrm{sn}}^2& =\gamma P_{\mathrm{n}}/{\rho }_{\mathrm{n}},\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill c_{\mathrm{st}}^2& =\gamma (P_{\mathrm{n}}+P_{\mathrm{p}})/({\rho }_{\mathrm{n}}+{\rho }_{\mathrm{p}}),\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill v_{\mathrm{A}}^2& =(B_x^2+B_y^2)/\rho ,\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill v_{\mathrm{At}}^2& =(B_x^2+B_y^2)/({\rho }_{\mathrm{p}}+{\rho }_{\mathrm{n}}),\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill {\widehat{v}}_A& =v_{\mathrm{At}}cos(\theta ),\hfill \end{array}$$(23)

$$\begin{array}{cc}\hfill {\theta }_B& =arctan(B_y/B_x),\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill V_{\mathrm{s}}^2& =\frac{1}{2}(c_{\mathrm{s}}^2+v_{\mathrm{A}}^2-\sqrt{{(v_{\mathrm{A}}^2+c_{\mathrm{s}}^2)}^2-4v_{\mathrm{A}}^2{c}_{\mathrm{s}}^{2}cos{({\theta }_B)}^2)}),\hfill \end{array}$$(25)

$$\begin{array}{cc}\hfill V_{\mathrm{st}}^2& =\frac{1}{2}(c_{\mathrm{st}}^2+v_{\mathrm{At}}^2-\sqrt{{(v_{\mathrm{At}}^2+c_{\mathrm{st}}^2)}^2-4v_{\mathrm{At}}^2{c}_{\mathrm{st}}^{2}cos{({\theta }_B)}^2)}),\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill V_{\mathrm{f}}^2& =\frac{1}{2}(c_{\mathrm{s}}^2+v_{\mathrm{A}}^2+\sqrt{{(v_{\mathrm{A}}^2+c_{\mathrm{s}}^2)}^2-4v_{\mathrm{A}}^2{c}_{\mathrm{s}}^{2}cos{({\theta }_B)}^2)}),\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill V_{\mathrm{ft}}^2& =\frac{1}{2}(c_{\mathrm{st}}^2+v_{\mathrm{At}}^2+\sqrt{{(v_{\mathrm{At}}^2+c_{\mathrm{st}}^2)}^2-4v_{\mathrm{At}}^2{c}_{\mathrm{st}}^{2}cos{({\theta }_B)}^2)}),\hfill \end{array}$$(28)

for the sound c_s, Alfvén v_A, slow V_s and fast V_f wave speeds. θ_B is the angle of the magnetic field. ${\widehat{v}}_{\mathrm{A}}$ is the component of the Alfvén velocity parallel to the x-axis and is necessary to calculate the shock jump conditions for intermediate shocks (see Sect. 3.3.1 for details).

For an entirely coupled system (α_c → ∞), the wave speeds depend on the total density (ρ_n + ρ_p). For a entirely decoupled system (α_c = 0) the wave speeds in the plasma depend on only the plasma density (ρ_p) and those in the neutrals depend on only the neutral density (ρ_n). Here we are interested in the region where a propagating wave undergoes a finite number of collisions as it crosses a pressure scale height. In this situation, the majority of the velocity drifts are expected to exist in the shock front.

A dispersion relation exists that can be numerically solved to calculate the effective wave speeds for a finitely coupled plasma Soler et al. (2013). We do not take this approach in this paper and instead we consider the wave speeds at the extremes of coupling defined by the above equations.

3.2. Wave steepening

Waves naturally steepen as they propagate through a stratified atmosphere, provided that the non-linear terms are greater than the dissipative effects. The rate at which this occurs depends on the scale height of the system. Figure 2 shows a wave steepening in our atmosphere for the MHD (solid) and plasma species in the PIP (dashed) cases at different snapshots. Here the PIP case has a neutral fraction of ξ_n = 0.9 at the lower boundary and the same driver amplitude, hence more applied kinetic energy since there is also velocity applied to the neutral species.

Fig. 2. Snapshots showing the wave steepening for the MHD (solid) and the plasma component of the PIP (dashed) cases with a vertical magnetic field at different times.

One can see from Fig. 2 that the amplitude of the waves is approximately equal after 0.5t. As time advances, the amplitude increases much more rapidly for the PIP case, with the velocity being around 20% larger at 4.5t. The increased rate of steepening can be explained using by the pressure scale heights. In the MHD case, the pressure scale height is constant away from the boundaries (between −3 ≤ x ≤ 3). For the PIP case, the pressure scale height for a coupled wave changes as it propagates from a mostly neutral to a mostly plasma medium.

It is also clear in Fig. 2 that the PIP shock propagates slower than the MHD shock. For a vertical magnetic field, the propagation speed of an isolated plasma wave is the sound speed (c_s) below the β_v point, and the slow speed (which aligned with the magnetic field is the sound speed) above the β_v point. An isolated neutral wave will propagate at the neutral sound speed (c_sn) through all heights in the atmosphere, which is a factor of $\sqrt{2}$ smaller than the plasma sound speed, see Fig. 3. Here, the coupling between the neutral and plasma species results in a propagation speed slower than for an isolated plasma (i.e. MHD), and faster than an isolated neutral species. This is discussed in more detail in later sections.

Fig. 3. Wave speeds for a PIP simulation with a vertical magnetic field. The sound (black), Alfvén (red) and slow (blue) speeds are shown for the initial atmosphere using the total densities (solid) and plasma density only (dashed).

3.3. Shock fundamentals

3.3.1. Shock transitions

Shocks are classified by comparing the upstream and downstream normal flow velocity v_⊥ to the characteristic speeds (e.g. Delmont & Keppens 2011):

• (1) superfast: V_f <  |v_⊥|,

• (2) subfast: ${\widehat{V}}_{\mathrm{A}}<|v_{\perp }|<V_{\mathrm{f}}$,

• (3) superslow: $V_{\mathrm{s}}<|v_{\perp }|<{\widehat{V}}_{\mathrm{A}}$,

• (4) subslow: 0 <  |v_⊥|< V_s,

• (∞) static: v_⊥ = 0.

Defining the upstream condition u and downstream condition d, several shocks of the form u → d are possible:

• 1 → 2 fast shocks,

• 3 → 4 slow shocks,

• 1 → 2 = 3 switch-on,

• 2 = 3 → 4 switch-off,

• 1 → 3, 1 → 4, 2 → 3, 2 → 4 intermediate shocks.

In the neutral species, there is only one characteristic wave speed, the sound speed c_s, and hence only a sonic shock can exist in this fluid. A purely hydrodynamic shock in our setup will only exist in the parallel component of velocity.

In partially ionised plasma, both neutral and plasma species exist and can interact through collisional coupling. Previous studies have shown that this can have several interesting effects, including the formation of shock substructure and additional shock transitions (Hillier et al. 2016; Snow & Hillier 2019).

3.3.2. Shock frame

To calculate the shock transitions accurately, the shock must be moved to the shock frame, where it is stationary in time. The shock frame here is defined as:

$$\begin{array}{c}\hfill x_{\mathrm{s}}=(x-x_{\mathrm{shock}})/t,\end{array}$$(29)

where x_shock is the location of the shock determined by the maximum absolute velocity gradient (if discontinuous), or the maximum absolute velocity (if smooth). The shock velocity is calculated using a quadratic derivative of the shock location in time.

4. Results

In this section we present the analysis of the compressional waves as they propagate through the stratified atmosphere with an inclined magnetic field. The velocity from our simulations is decomposed into parallel and perpendicular components relative to the instantaneous magnetic field. This is approximately a decomposition into slow (parallel) and fast (perpendicular) wave modes above the β_v point (i.e. x >  0 in our model). Below the β_v point (x <  0) this approximation does not hold and only an acoustic shock exists.

4.1. Vertical magnetic field: B_x = 1.0, B_y = 0.0

4.1.1. MHD

For a vertical magnetic field (B_x = 1,  B_y = 0) the wave steepens into a shock below the β_v point. A vertical magnetic field acts as a reference case since the propagation is restricted to be along the magnetic field direction and v_⊥ = 0. The wave steepens with height, forming a shock below the β_v point, see Fig. 4. The shock front is characterised by the jump in the parallel velocity and density.

Fig. 4. MHD (left) and PIP (right) time series for a vertical field (Bx = 1,  By = 0), for the plasma (black) and neutral (red) species. Velocity parallel (solid) and perpendicular (dashed) is plotted and corresponds to the left axis. The density increase normalised by the initial density profile is shown as the dot-dash line and corresponds to the right axis.

Figure 5 shows the shock transitions corresponding to the last frame of the time series (Fig. 4), where the horizontal axis has been shifted into the shock frame (as described in Sect. 3.3). The shock transitions show that this is a slow-mode shock (i.e. a transition from above to below the slow speed). This is as expected for this atmosphere, as above the β_v point only the slow mode transition parallel to the field is permitted.

Fig. 5. MHD shock transitions for the vertical magnetic field (Bx = 1,  By = 0) at t = 5.68 (left). Horizontal axis has been transferred to the shock frame xs. PIP shock transitions for the vertical magnetic field (B − x = 1,  By = 0) for the parallel shock transitions for the plasma (black) and neutral (red) species at t = 5.90 (right).

4.1.2. PIP

For the PIP case, the overall shock behaviour is similar to the MHD case however there are interesting features localised around the shock front, Fig. 4. The wave steepens into a shock below the β_v line showing sharp discontinuities in both the plasma (black) and neutral (red) species. The perpendicular component of both species is zero, as expected. A sharp rise in density is present at the shock front in both species. The neutral density increases by a larger relative amount than the plasma, however the neutral density is very small here (see Fig. 1).

There is a slight lag between the maximum velocities in the two species, with the neutral fluid slightly behind the plasma. Physically, the plasma sound speed c_s is larger than the neutral sound speed c_sn, see Fig. 3. For isolated species, a plasma shock would be expected to propagate faster in this atmosphere (see Fig. 2). High in the atmosphere, the neutral fraction is very small and hence the collisional coupling is weaker. As such, the neutral species is dragged by the plasma via collisional coupling.

The shock transitions (Fig. 5) for the PIP case show a show mode shock in the plasma and a sonic shock in the neutral species. We note that the shock transitions are only applicable around x_s  =  0 and the transition at x_s  ≈   − 0.017 is an artefact from the limitations of the shock fit and not a real shock transition.

Figure 6 shows the temperature across the domain (a) and the heating terms (b), where the frictional heating is defined as α_c(T_n, T_p)ρ_nρ_p(v_n − v_p)²/2 and thermal damping is 3α_c(T_n, T_p)ρ_nρ_p(P_n/ρ_n − P_p/(2ρ_p))/2. The shock causes a temperature increase in both plasma and neutral species, however the temperature in the neutrals increase comparatively more. Two reasons exist for the higher neutral temperature: the neutrals undergo greater compression because the wave is dragged faster than it would travel in a single fluid simulation, and the lower neutral density at this point which allows for a greater temperature increase from the frictional heating.

Fig. 6. Top: temperature in the PIP simulation for the plasma (black) and neutral (red) species. The black dashed line shows the MHD temperature at the same time. Lower: frictional heating (black) and thermal damping (red and blue). Red line indicates ions losing heat to neutrals. Vice versa for the blue line.

The MHD temperature at the same time is also shown in Fig. 6 (black dashed line). The PIP simulation has a higher temperature than the MHD case. This is likely to be as a result of the increased non-linearity of the wave created by the smaller pressure scale height of the neutrals.

The drift velocities (v_n − v_p) are shown for a few different times in Fig. 7. As the shock propagates, the drift velocities increase. At time τ = 5.50 the drift velocity is fairly close to the equilibrium neutral sound speed ($1/\sqrt{2}$). It should be noted that such large drift velocities are approaching the limits of applicability of our drift velocity model (Eq. (14)), because once the drift velocity becomes much greater than the sound speed it is the drift and not thermal motions that drive the collisions between the fluids (Draine 1986).

Fig. 7. Drift velocities at different times for the vertical magnetic field (Bx = 1.0,  By = 0.0) at different times, corresponding to the snapshots in Fig. 4.

4.2. Inclined field: B_x = 0.8, B_y = 0.6

When the initial conditions feature an angle in the magnetic field, energy can be exchanged between the different wave modes at the β_v point and the shock fronts can separate into their fast (perpendicular) and slow (parallel) components (studied in MHD by Pennicott & Cally 2019). Here we present results from a inclined magnetic field with B_x = 0.8,  B_y = 0.6, preserving the total magnetic field strength as unity. The initial pressure and density profiles are identical to the vertical field (see Fig. 1), as is the driver.

4.2.1. MHD case (B_x = 0.8, B_y = 0.6)

In the MHD case (see left column of Fig. 8), a notable difference (compared to the vertical field) is the non-zero component of velocity perpendicular to the magnetic field (i.e. the fast-mode component). Below the β_v point, there is a single shock front that consists of both perpendicular and parallel components of velocity, see Fig. 8. As the shock propagates through the β_v point, the energy is exchanged between modes and the fast- and slow-mode separate. Notably, the parallel (slow) component of the velocity becomes smoothed: below x = 0 there is a steep discontinuity, above x = 0 the wavefront becomes smoothed and no longer satisfies the conditions to be called a shock. The perpendicular (fast) component retains its steep shock structure and as a consequence of being a shock there is a steep change in the density relative to the initial density (Fig. 8). A small non-zero parallel velocity component exists at the location of the fast-mode shock which is an artefact from the assumption of decomposing into parallel and perpendicular velocities, (see, for example, Pennicott & Cally 2019).

Fig. 8. MHD (left) and PIP (right) time series for an inclined field with Bx = 0.8,  By = 0.6 for the plasma (black) and neutral (red) species. Solid (dashed) line shows the parallel (perpendicular) component of velocity relative to the instantaneous magnetic field. Density normalised by the density at t = 0 is plotted as the dot-dashed line and corresponds to the right axis.

The shock transitions have been calculated for the wide inclined field, B_x = 0.8,  B_y = 0.6, and are shown in Fig. 9. In this case, there is no shock transition associated with the parallel component of velocity. In the perpendicular component of velocity, there is a clear transition from superfast to subfast, that is, a fast-mode shock. We note that shock frame here was calculated using the maximum absolute gradient for the perpendicular velocity and the maximum absolute velocity for the parallel component, see Sect. 3.3.2.

Fig. 9. Shock transitions for the MHD (left panel) and PIP (right panel) shocks with an inclined field (Bx = 0.8,  By = 0.6). Parallel (solid) and perpendicular (dashed) components are plotted in their relative shock frames. The neutral transitions (red) correspond to the right axis.

4.2.2. PIP case (B_x = 0.8, B_y = 0.6)

Below the β_v point, the plasma and neutral species are well coupled, as are the parallel and perpendicular components of each species, see Fig. 8. Changes arise above the β_v point. The plasma components separates into its parallel (slow) and perpendicular (fast) components, with the parallel components becoming smoothed and the perpendicular component remaining sharp (as seen in the MHD case). The shock transitions (Fig. 9) reveal a fast-mode shock in the perpendicular component of the plasma species, with no shock transitions present in the parallel plasma velocity or the neutral species.

The neutral species exhibits a non-zero perpendicular velocity that propagates with (but behind) the plasma perpendicular velocity. As the neutral species does not have an associated fast speed, the perpendicular component of the neutral species arises from collisional coupling with the plasma species. This results in a reasonably large increase in the neutral temperature (compared to the plasma temperature) around the fast shock, Fig. 10. Here, the increase in neutral temperature is at the location of the fast-mode shock only. The slow component has been smoothed and the neutral and plasma species are reasonably well coupled, unlike the vertical field where the species decouple at the shock front and there is an associated increase in frictional heating of the neutrals (see Figs. 4 and 6).

Fig. 10. Inclined field: (left) temperature in the PIP simulation for the plasma (black) and neutral (red) species, and the MHD case (black dashed) Note that the MHD is at τ = 4.31 and the PIP simulation is at time τ = 4.80. Right: frictional heating (black) and thermal damping (red and blue). Red line indicates ions losing heat to neutrals. Vice versa for the blue line.

The drift velocities (v_n − v_p) are shown in Fig. 11 for a few different times. The drift velocity is fairly large at the location of the fast-mode shock. The drift velocity associated with the parallel component of velocity is relatively small. For a linear wave, coupling is much more effective than a non-linear wave. The parallel component is smoothed and no longer satisfies the conditions for a shock above the β_v point, and the plasma and neutral species are reasonably well-coupled.

Fig. 11. Parallel (solid) and perpendicular (dashed) drift velocities for the wide inclined field (Bx = 0.8,  By = 0.6) at different times, corresponding to the snapshots in Fig. 8.

5. Effect of different neutral fractions and collisional coupling

The cases considered thus far are using a neutral fraction of ξ_n = 0.9 at the lower boundary with the collisional coefficient taken to be α₀ = 100. This profile becomes plasma-dominated at x ≈ −1, which is below the β_v point. As such, the neutrals effect on the plasma is relatively small near the top of this domain. By increasing the initial neutral fraction, the neutral species has a larger impact on the plasma through increased collisions. A large neutral fraction, along with larger values of α₀, is a regime that is closer to that of the solar chromosphere making it an important regime to study. In this section, larger neutral fractions and a range of α₀ values are investigated, but first we discuss subsequent peculiarities characterising this set-up.

5.1. Different β_v locations

In single-fluid MHD, the β_v point is singularly defined since there is only one sound speed, and one Alfvén speed. In the two-fluid case, one can consider the β_v point using isolated plasma speeds, or the coupled speeds (see Sect. 3.1) here referred to as β_vt (where c_st = v_At). As the neutral fraction is increased, the relative locations of the β_v and β_vt points change.

Figure 12 shows the initial atmospheres for different neutral fractions specified at the base of the domain, specifically ξ_n = 0.9, 0.99, 0.999. The normalisation means that the plasma density is equal to the MHD density in all atmospheres and a higher neutral fraction adds more mass to the system. The solid vertical lines show the β_v and β_vt points for all atmospheres. The black line uses the plasma speeds only (c_s = v_A) and hence is the same for all atmospheres here, whereas the coloured lines use the total speeds (c_st = v_At) and varies with the neutral fraction. For ξ_n = 0.9, the two lines are fairly close together (red for coupled, black for MHD). For larger ξ_n values, the coupled β_v occurs much higher in the atmosphere than the plasma β_v.

Fig. 12. Plasma (solid black) and neutral (dashed) densities with height for ξn0 = 0.9 (red), 0.99 (blue), and 0.999 (green). The MHD density is given by the black line and the plasma component of the two-fluid simulation is identical to the MHD density profile. The vertical line indicate the βv point (where cs = vA) for the plasma species (and MHD) only. The coloured vertical lines show the βvt points using the total densities (ρn + ρp), that is, cst = vAt.

The previous sections of this paper use a neutral fraction of ξ_n = 0.9 and there is very little difference between the β_v and β_vt points (Fig. 12). The obvious question here is: does the mode conversion occur at the β_v (c_s = v_A) or the β_vt (c_st = v_At) point? This is investigated by varying the collisional coefficient and the neutral fraction.

Similarly, care must be taken when defining the wave speeds; for a weakly coupled case, the plasma wave speeds (c_s, v_A, V_s, V_f) are likely most appropriate, whereas for a strongly coupled case the combined wave speeds (c_st, v_At, V_st, V_ft) should be more suitable, see Sect. 3.1. For finite coupling, it is not immediately obvious which wave speeds are most appropriate. For the case where ξ_n = 0.9 at the base of the domain, the neutral density and pressure are relatively small compared to the plasma values above the β_v point, hence there is not too much of a difference. Choosing the appropriate wave speeds becomes more important when the neutral fraction is higher, see Fig. 12 where the neutral and plasma densities are approximately equal near the top of the domain.

5.2. Effect of different collisional coefficients

For the simulations thus far, the collisional coefficient is set as α₀ = 100. Changing this value modifies the number of collisions per unit time and hence determines the level of coupling between the neutral and plasma species. Taking α₀ = 0 would yield the same result as the reference MHD calculation in the plasma since the two species propagate independently. At the other extreme, taking α₀ = ∞ would also yield an MHD result but the wave speeds and the rate at which non-linearities develop would all depend on the total density (ρ_p + ρ_n). Between these extremes, the coupling is finite and one would expect two-fluid effects to be most important. One may also expect separation of wave modes to occur between β_v and β_vt.

Simulations have been performed using B_x = 0.8,  B_y = 0.6 (corresponding to Sect. 4.2), using different values of collisional coupling, and a neutral fraction of ξ_n = 0.999. We extend the domain to allow sufficient height above the β_vt point. The gravity in this extended region is zero. This prevents the fast mode wave from accelerating above x ≈ 3.7. The new atmosphere is shown in Fig. 13.

Fig. 13. Long atmosphere used for the collisional coefficient simulations. Black line shows the plasma (and MHD) density. Red dashed line shows the neutral density for an initial neutral fraction of ξn = 0.999. The βv and βvt points are shown by the black and green vertical lines.

Figure 14 shows the PIP simulations using different collisional coefficients, and the MHD case for reference. The corresponding shock transitions are included in Appendix C for reference. The MHD results shows a sharp perpendicular (fast-mode) shock, and a smoothed parallel (slow-mode) component, as discussed in Sect. 4.2.1. For the weakly collsional case (α₀ = 0.001), the plasma component of the PIP simulation is very similar to the MHD simulations, in terms of both features and amplitudes. The plasma component is very weakly influenced by the neutral species and therefore they propagate fairly independently and the plasma wave modes separate at the β_v point. A fast transition occurs in the perpendicular plasma velocity using the characteristic wave speeds of the isolated plasma. There is no parallel plasma shock transition since this has been smoothed (as seen in the MHD case). In the neutral species, a sonic shock exists in both the parallel and perpendicular components of velocity. The neutral parallel and perpendicular shocks are located at the same height since the neutrals are behaving independently of the plasma, and hence independently of the magnetic field.

Fig. 14. Snapshots of the parallel (solid) and perpendicular (dashed) velocities for the plasma (black) and the neutral (red) species for different collisional coefficients.

As the collisional coefficient increases (α₀ = 0.1) changes start to appear in the plasma component of the PIP simulation, compared to the MHD case; the amplitude of the fast- and slow- mode shocks have increased, and the propagation time has slowed. The increase in shock amplitude and propagation time was described in Sect. 3.2 and is due to the different pressure scale heights of the neutral and plasma species. The same shock transitions exist as for α₀ = 0.001 in the neutral and perpendicular plasma velocities. A new slow-mode shock exists in the parallel plasma velocity. The increased coupling of the system has lead to a steepening of the slow-mode component through the increased compressionallity of the system.

At α₀ = 1, we start to see similarities between the plasma and neutral velocities parallel to the instantaneous magnetic field. The parallel plasma velocity corresponds to the slow component and propagates at a comparable speed to the neutral parallel velocity, indicating an increased coupling over the lower α₀ cases. The plasma perpendicular velocity is clearly responding to the neutral shock even for this low value of α_c, as evidenced from the change in shape of the shock. This is partially a result of the increased neutral density and temperature post shock increasing the coupling of the plasma to the neutrals through increased collisions. The perpendicular velocity corresponds to the fast-mode shock and propagates much faster than the parallel velocity. As a result, the neutrals feel less collisions from the plasma as this wave propagates past them, resulting in a much weaker coupling between the species in this wave front compared to the parallel (slow) component. The perpendicular neutral velocity appears unaffected by the plasma. The shock transitions for α₀ = 1 show a sonic shock in the neutral species and a fast-mode transition in the perpendicular plasma velocity, as seen for α₀ = 0.1 (see Appendix C). However, the parallel plasma velocity shows two different transitions depending on the wave speeds chosen: a slow-mode shock (3 → 4) in the plasma wave speeds, and an intermediate transition (2 → 3) in the total wave speeds. Physically, we can deduce that the appropriate choice here is the slow-mode shock using the plasma properties since there is no reversal of the magnetic field (which is a requirement for an intermediate transition). This highlights the importance of choosing the appropriate wave speeds in two-fluid plasma.

At α₀ = 100, the parallel neutral and plasma velocities are well-coupled and the plasma perpendicular velocity creates a drag on the neutrals driving a neutral perpendicular velocity. This is similar to the results in Sect. 4.2.2 for a neutral fraction of ξ_n = 0.9. Here, the plasma perpendicular shock is fairly broad and smoothed (compared to the discontinuous jump in the MHD case). The coupling is sufficient that no shock transitions exist in the parallel neutral or plasma velocity. The slow-mode component in the plasma is smoothed and the coupling transfers this effect to the neutral species, smoothing the parallel sonic shock. A sonic shock exists in the perpendicular neutral velocity, located at an overshoot in the perpendicular neutral velocity, see Fig. 14. As the coupling increases, the tail end of the fast-mode shock begins to couple with the neutrals, leading to a large frictional heating of the neutral species and a pressurised expansion in the neural species, creating the overshoot in perpendicular neutral velocity. The perpendicular plasma velocity shows a shock transition when the plasma wave speeds are used but not when the total wave speeds are used. The shock front is decoupled from the neutrals so one would expect that the plasma wave speeds are most appropriate here. This level of coupling presents an interesting case where the slow-mode plasma velocity is strongly coupled to the neutral species, however the fast-mode plasma velocity is only partially coupled to the neutrals. The fast-mode component travels much faster than the slow mode and hence experiences less collisions per unit distance.

Increasing the collisional coefficient further (α₀ = 1000) encourages more coupling between the neutral and plasma species. The perpendicular drift velocity between the two species is small and the plasma component is broader than the MHD case. This shows a significant increase in the level of coupling since both the fast- and slow-mode plasma velocities are strongly coupled to the neutral species. Here, no shock transitions exist in the neutral species. In the perpendicular plasma velocity two possible shocks exist depending on the wave speeds taken; an intermediate transition (2 → 3) for the plasma wave speeds, and a fast-mode shock (1 → 2) using the total wave speeds. We can rule out the intermediate shock since there is no reversal of magnetic field and deduce that here the total wave speeds are most appropriate and a fast-mode shock exists. This indicated that the perpendicular plasma velocity is strongly coupled to the neutral species and the system is now behaving like an MHD case.

Figure 15 shows the perpendicular shock width for the plasma (calculated as the distance between the locations of the maximum velocity and the maximum gradient of velocity) as a function of the collisional coefficient. The shock width is calculated when the perpendicular velocity has reached x ≈ 6.5, corresponding to the snapshots in Fig. 14. For α₀ <  0.1 (and the MHD case), the shock is effectively discontinuous. A small numerical shock width exists due to numerical dissipation of approximately 0.02 pressure scale heights, see the discussion in Appendix B. For collisional coefficients in this range, the coupling between plasma and neutral species is weak for the fast-mode component and the perpendicular plasma velocity in the PIP simulation is similar to the MHD simulation, see Fig. 14. Between 0.1 <  α₀ ≤ 30 the perpendicular shock becomes increasingly broad as the fast-mode plasma shock begins to couple with the neutral species. The neutral drag causes a change in the shape of the plasma fast-mode shock.

Fig. 15. Perpendicular shock widths for different collisional coefficients. Corresponds to the snapshots in Fig. 14. Shock width is calculated as the distance between the maximum velocity, and the maximum absolute gradient of velocity.

A sharp increase in shock width occurs for 30 ≤ α₀ ≤ 100 with a maximum occurring around α₀ = 100. This occurs because the fast-mode shock has not fully separated from the slow-mode component. Performing the α₀ = 100 simulation in an atmosphere that is stratified up to x = 6 allows the parallel and perpendicular velocities to separate and the new shock width is 1.19 indicating that even with the reduced total collisions in a fully stratified atmosphere, one can expect shock widths of larger than the pressure scale height for finitely coupled systems.

As the collisional coefficient increases further (α₀ ≥ 1000) the shock width decreases and the simulation tends towards an MHD simulation with the dynamics determined by the total density. It should be noted that while these shocks are very broad, they still satisfy the condition for a fast-mode shock to exist, namely a transition from superfast to subfast across the leading edge of the shock. Finite shock widths are a fundamental feature of two-fluid shocks and arise due to the separation of ionised and neutral species at the shock front (e.g. Hillier et al. 2016; Snow & Hillier 2019). At the leading edge of the shock, the plasma separates from the neutral species and has a superfast to subfast transition. Following this, the collisional coupling with the neutral species creates a drag that shifts the maximum velocity to far behind the leading edge, resulting in a finite width of the shock that can exceed the pressure scale height.

The different collisional coefficients change the behaviour of the shock propagation, from MHD-like when the collisional coefficient is small (α₀ <  1), to a partially-coupled regime (1 <  α₀ ⪅ 1000), and to a strongly coupled system (1000 ⪅ α₀) which is again MHD-like but the properties depend on the total density (ρ_p + ρ_n). One would expect the mode separation location to also depend on the collisionality of the system. To investigate the height at which the wave modes separate, the location of the maximum velocity is plotted through time for the parallel and perpendicular components for different collisional coefficients, see Fig. 16. Here the neutral fraction at the base of the domain is ξ_n = 0.999, i.e., the atmosphere in Fig. 13.

Fig. 16. Shock locations through time for different collisional coupling coefficients. Parallel (solid) and perpendicular (dashed) velocities and shown for the plasma (black) and neutral (red) species. The βv and βvt locations are shown as the green and blue line respectively. The graph is plotted until the perpendicular shock approaches the top of the domain.

In the MHD case, the wave modes separate at the β_v point. Similarly, for a weakly coupled system (α ⪅ 3) the coupling is sufficiently weak that the physics is mostly determined by the plasma component such that the separation occurs at approximately the β_v point.

For the highly collisional case (α = 1000) the parallel and perpendicular separate above the β_vt point indicating that the system is strongly coupled and the total densities are determining the physics. This can be considered as MHD-like since there is such strong coupling between the neutral and plasma species. The perpendicular velocity in the neutral species is strongly coupled to the perpendicular plasma velocity.

Between these limits, the separation happens between the β_v and β_vt points. At α₀ = 10, the separation occurs just above the β_v point. At α₀ = 100 the separation occurs just below the β_vt point. The finite coupling between the neutral and plasma species causes the two-fluid effects to become important and the point at which wave mode separate is less trivial to determine.

As with Pennicott & Cally (2019), we find that the mode conversion height is advected with slow-mode shock for a short time after the mode separation, as shown in Fig. 16. A larger parameter study would be required to analyse the consequences of this since it is likely related to the velocity amplitude.

5.3. Effect of neutral fraction

In this section, simulations are now performed using three initial neutral fraction of ξ_n = 0.9,  0.99,  0.999 at the base of the domain and the collisional coefficient is set to α₀ = 100. The domain is extended for all atmosphere with gravity g = 0 in the extended region, see Fig. 13 for the ξ_n = 0.999 case. Two angles are investigated: vertical field (B_x = 1.0,  B_y = 0.0) and inclined field (B_x = 0.8,  B_y = 0.6). Snapshots of the data are taken when the slow shock reaches x = 3 for the vertical field, and when the fast-mode shock reaches x = 6.5 for the wide inclined field. The parallel and perpendicular velocities for the plasma and neutral species are shown in Fig. 17.

Fig. 17. Parallel (solid) and perpendicular (dashed) velocities in the plasma (black) and neutral (red) species for the vertical magnetic field (top row) and inclined field (bottom row) for different neutral fractions (ξn = 0.9, 0.99, 0.999 across the columns). The vertical lines indicate the βv (black) and βvt (green) points.

For the vertical field, increasing the neutral fraction at the base does not significantly change the results. There is increased compression of the of the shock due to the increased coupled pressure scale height since, for this setup, a higher neutral fraction has more total mass.

For the inclined field (B_x  =  0.8,  B_y  =  0.6) the results are more different. Again, there is an increased wave amplitude due to the increased compression. For the parallel (slow) plasma velocity, the main difference in the increased neutral fraction is an increased amplitude. The parallel plasma velocity couples well with the parallel neutral velocity. The perpendicular plasma velocity shows increased coupling with perpendicular neutral velocity for higher neutral fractions, resulting in a broader shock front. This shock broadening is the same result as seen for increased collisional coefficient (see Sect. 5.2). The momentum and energy exchange between the plasma and neutral species depend on the density and collisional coefficient. Here the total density increases as the neutral fraction increases so the coupling terms become larger, hence increasing the coupling coefficient and the neutral fraction have similar effects in this set-up.

6. Observational consequences

The results of this paper indicate that the widths of fast-mode shocks in two-fluid stratified atmospheres can regularly exceed the pressure scale height. Relating this study to the solar chromosphere, the temperature is roughly constant, and the pressure scale height is L₀ ≈ 300 km, and one would expect fast-mode shocks in a finitely-coupled two-fluid medium to be of approximately this width. This is very different from the intuitive image of shocks as a steep discontinuity and provides a potential observable of two-fluid effects in the lower solar atmosphere.

To relate these results to the solar atmosphere, our results can be dimensionalised using the pressure scale height of the chromosphere (L₀ ≈ 300 km), and the sound speed (c_s0 ≈ 8 km s⁻¹), leading to a time scale t₀ ≈ 40 s. Sampling the Doppler velocity along the direction of stratification (v_xp) at a fixed height above the β_vt point, we find that the fast-mode shock should appear as a gradual increase in Doppler velocity over approximately 6 s. As such, with high temporal cadence and a narrow spectral line, one would expect this feature to be observationally resolvable.

In addition, previous studies have shown that another potential observable of two-fluid effects is the electric field felt by the neutrals (Anan et al. 2014; Snow & Hillier 2019). The large drift velocities lead to polarisation of the neutrals as they flow past the magnetic field, and hence have an associated electric field term of the form E_n = (v_n − v_p)×B. The large drift velocities (up to 80% of the plasma sound speed) present over the large shock width in the fast-mode shock provide an excellent scenario for observing the electric field felt by the neutrals.

7. Conclusions

In this paper, we investigate the propagation and mode conversion of two-fluid shocks in an isothermal, stratified atmosphere. The collisionallity of the system is found to have a strong effect on determining the behaviour of the system. Two angles of magnetic field are investigated: vertical (B_x = 1,  B_y = 0), and inclined (B_x = 0.8,  B_y = 0.6).

There are several quirks that characterise two-fluid simulations that are important in this study. Firstly, the neutral and ionised species have different pressure scale heights, meaning that a coupled two-fluid wave propagates slower and steepens more than an isolated MHD wave. Secondly, the characteristic wave speeds of a two-fluid plasma are multiply defined depending on the level of coupling, see Sect. 3.1. Finally, the point at which mode conversion occurs is multiply defined: β_v (c_s = v_A) using the MHD (or plasma) properties, and β_vt (c_st = v_At) using the properties of the coupled system (plasma + neutrals).

For a vertical magnetic field (B_x = 1,  B_y = 0) the initial wave steepens into a slow-mode shock. When the shock passes the β_v point in the MHD simulation, there is no significant changes in the shock since the propagation is restricted to along the magnetic field direction. In the PIP case (for ξ_n = 0.9), there is a separation of the neutral and plasma species that appears at the leading edge of the shock above the β_v point. The plasma sound speed is larger than the neutral sound speed by a factor of $\sqrt{2}$ hence the plasma shock propagates ahead of the neutral shock. The frictional coupling between the two species forces the neutral species to propagate only slightly behind the plasma, and there is an associated increased temperature of the neutral species (Fig. 6). Increasing the neutral fraction increases the amplitude of the shock due to the increased compressionality of the system. The propagation speed of the shock remains the same since the sound speed does not change by increasing the neutral fraction.

When the magnetic field is highly inclined (B_x = 0.8,  B_y = 0.6) the initially coupled wave separates into fast- and slow-mode components above the β_v point. In the MHD case, this results in a smoothing of the parallel (slow) component of the velocity (see Fig. 4), as analysed by Pennicott & Cally (2019). For the two-fluid PIP case (for ξ_n = 0.9,  α₀ = 100), the smoothing of the parallel (slow) shock component results in a strong coupling between the neutral and plasma species due to the smoothing and does not possess a shock transition in either the neutral or plasma species. The perpendicular (fast) component of velocity exhibits a reasonably strong drift velocity near the shock front. The neutral species does not possess an MHD fast wave speed (since it is hydrodynamic and only supports sonic shocks). The perpendicular component of the neutral species arises from the coupling with plasma fast shock. As such, there is a lag in the peak velocities of the plasma and neutral perpendicular velocities (Fig. 8) and strong increase in neutral temperature due to the frictional heating (Fig. 10).

Increasing the neutral fraction and collisional coefficient drastically changes the physics of the system, see Sect. 5. For ξ_n = 0.999, the collisional coefficient results in a MHD-like decoupled system (α₀ ≤ 0.1), a finitely-coupled system (0.1 ≤ α₀ ≤ 100), and finally a MHD-like coupled system (100 ≤ α₀). In the extremes of low α₀ ≤ 0.1, the plasma and neutral species propagate fairly independently, with the plasma fast-mode virtually unaffected by the neutrals, and the plasma slow-mode only weakly affected. The fast-mode propagates much more quickly that the slow-mode and hence experiences fewer collisions per unit distance. As α₀ increases, the parallel plasma and neutral species become more coupled. The perpendicular plasma shock becomes increasingly broad as it experiences more collisions with the plasma and results in a shock width of larger than the pressure scale height. For α₀ = 1000, both the parallel and perpendicular plasma velocity modes are well coupled to the neutral species and the system becomes MHD-like.

The location at which mode-conversion occurs also depends on the collisionallity of the system. In the MHD case and PIP α₀ ≤ 1, the fast- and slow-modes separate at the β_v point since the dominant separation is determined by the isolated plasma. For a α₀ = 1000 the modes separate above the β_vt point since the system is strongly coupled and MHD-like. between these two limits, the system is finitely coupled and the mode conversion occurs between the β_v and β_vt points.

Relating our results to the solar chromosphere, one would expect fast-mode shocks in finitely-coupled, two-fluid mediums to have a width comparable to the pressure scale height (≈300 km in the chromosphere). Such broad shocks could be a indirect way of observing two-fluid effects in the lower solar atmosphere.

In summary, two-fluid shocks in an isothermal stratified atmosphere feature several peculiarities compared to MHD shocks. The mode-conversion point, shock width, and the shock transitions present all depend on the collisionality of the system. The results presented in this paper aid in the understanding of shocks in the lower solar atmosphere.

Acknowledgments

BS and AH are supported by STFC research grant ST/R000891/1. AH acknowledges the support of STFC Ernest Rutherford Fellowship grant number ST/L00397X/2. The authors wish to acknowledge scientific discussions with the Waves in the Lower Solar Atmosphere (WaLSA; www.WaLSA.team) team, which is supported by the Research Council of Norway (project number 262622).

References

Appendix A: Numerical implementation of velocity driver

The PIP code evolves conserved quantities. As such, the numerical implementation of the driver described by Eqs. (15)–(17) in normalised form is:

$$\begin{array}{cc}& c_{\mathrm{s}0}=\sqrt{\gamma P_0/{\rho }_0},\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& v_{\mathrm{pert}}=Asin(2\pi t/\tau ),\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& v_{\mathrm{pert}}^{\prime }=\frac{2\pi A}{\tau }cos(2\pi t/\tau ),\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& P_{\mathrm{pert}}^{\prime }=\frac{P_0\gamma }{c_{\mathrm{s}0}}{v}_{\mathrm{pert}}^{\prime },\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& {\rho }_{\mathrm{pert}}=\frac{{\rho }_0}{c_{\mathrm{s}0}}{v}_{\mathrm{pert}},\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& {\rho }_{\mathrm{pert}}^{\prime }=\frac{{\rho }_0}{c_{\mathrm{s}0}}{v}_{\mathrm{pert}}^{\prime },\hfill \end{array}$$(A.6)

$$\begin{array}{cc}& \rho =\rho +{\rho }_{\mathrm{pert}}^{\prime }\mathrm{d}t,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}& \rho v_x=\rho v_x+({\rho }_{\mathrm{pert}}^{\prime }{v}_{\mathrm{pert}}+({\rho }_{\mathrm{pert}}+{\rho }_0)v_{\mathrm{pert}}^{\prime })\mathrm{d}t,\hfill \end{array}$$(A.8)

$$\begin{array}{cc}& e=e+(\frac{1}{2}{v}_{\mathrm{pert}}^2{\rho }_{\mathrm{pert}}^{\prime }+({\rho }_{\mathrm{pert}}+{\rho }_0)v_{\mathrm{pert}}{v}_{\mathrm{pert}}^{\prime }+\frac{P_{\mathrm{pert}}^{\prime }}{(\gamma -1)})\mathrm{d}t,\hfill \end{array}$$(A.9)

where dt is the simulation time step, and ρ,  ρv_x,  e are the evolved simulation density, x momentum and energy on the driven boundary.

Appendix B: Numerical diffusion

We use a first order HLLD numerical scheme, as described by Miyoshi & Kusano (2005). The (PIP) solver have been used with this scheme previously to model MHD shocks in a non-stratified atmosphere (Hillier et al. 2016; Snow & Hillier 2019) and the numerical diffusion of the scheme was found to be small, with the shock existing across a few grid cells, as shown in Fig. B.1. In the stratified atmosphere used in this paper, the shock is found to be more diffuse since the densities either side of the shock are constantly changing and hence a true steady-state solution does not exist. This is demonstrated in Fig. B.1 which shows the MHD slow-mode shock for a vertical magnetic field and corresponds to the MHD simulation in Fig. 4 at time t = 5.68. In the stratified case, the physical shock width is approximately 0.02 pressure scale heights and hence we consider the numerical diffusion a negligible effect.

Fig. B.1. Numerical diffusion tests for a slow-mode shock in non-stratified (top) and stratified (bottom) atmospheres. Top panel: MHD slow shock used in Hillier et al. (2016); bottom panel: MHD simulation in Fig. 4 at time t = 5.68. The numerical shock width in the stratified atmosphere is less then 0.02 pressure scale heights and considered negligible.

Appendix C: Shock transitions for the collisional tests

Figure C.1 shows the shock transitions for the different collisional fraction in Sect. 5.

Fig. C.1. Shock transitions for different collisional coefficients. Parallel (solid) and perpendicular (dashed) lines are shown for the plasma (black and green) and neutral (red) species. The black lines show the transitions relative to the isolated plasma, whereas the green lines show the transitions relative to the fully coupled system.

All Figures

	Fig. 1. Density for the initial equilibrium. MHD density is the black-solid line. The PIP simulation has two densities: plasma (blue dashed) and neutral (red dashed). The normalisation leads to the MHD density being the same as the PIP plasma density. The total density for the PIP case is shown in green-dashed. Gravity is constant between the two vertical black lines.
In the text

	Fig. 2. Snapshots showing the wave steepening for the MHD (solid) and the plasma component of the PIP (dashed) cases with a vertical magnetic field at different times.
In the text

	Fig. 3. Wave speeds for a PIP simulation with a vertical magnetic field. The sound (black), Alfvén (red) and slow (blue) speeds are shown for the initial atmosphere using the total densities (solid) and plasma density only (dashed).
In the text

	Fig. 4. MHD (left) and PIP (right) time series for a vertical field (Bx = 1,  By = 0), for the plasma (black) and neutral (red) species. Velocity parallel (solid) and perpendicular (dashed) is plotted and corresponds to the left axis. The density increase normalised by the initial density profile is shown as the dot-dash line and corresponds to the right axis.
In the text

	Fig. 5. MHD shock transitions for the vertical magnetic field (Bx = 1,  By = 0) at t = 5.68 (left). Horizontal axis has been transferred to the shock frame xs. PIP shock transitions for the vertical magnetic field (B − x = 1,  By = 0) for the parallel shock transitions for the plasma (black) and neutral (red) species at t = 5.90 (right).
In the text

	Fig. 6. Top: temperature in the PIP simulation for the plasma (black) and neutral (red) species. The black dashed line shows the MHD temperature at the same time. Lower: frictional heating (black) and thermal damping (red and blue). Red line indicates ions losing heat to neutrals. Vice versa for the blue line.
In the text

	Fig. 7. Drift velocities at different times for the vertical magnetic field (Bx = 1.0,  By = 0.0) at different times, corresponding to the snapshots in Fig. 4.
In the text

	Fig. 8. MHD (left) and PIP (right) time series for an inclined field with Bx = 0.8,  By = 0.6 for the plasma (black) and neutral (red) species. Solid (dashed) line shows the parallel (perpendicular) component of velocity relative to the instantaneous magnetic field. Density normalised by the density at t = 0 is plotted as the dot-dashed line and corresponds to the right axis.
In the text

	Fig. 9. Shock transitions for the MHD (left panel) and PIP (right panel) shocks with an inclined field (Bx = 0.8,  By = 0.6). Parallel (solid) and perpendicular (dashed) components are plotted in their relative shock frames. The neutral transitions (red) correspond to the right axis.
In the text

	Fig. 10. Inclined field: (left) temperature in the PIP simulation for the plasma (black) and neutral (red) species, and the MHD case (black dashed) Note that the MHD is at τ = 4.31 and the PIP simulation is at time τ = 4.80. Right: frictional heating (black) and thermal damping (red and blue). Red line indicates ions losing heat to neutrals. Vice versa for the blue line.
In the text

	Fig. 11. Parallel (solid) and perpendicular (dashed) drift velocities for the wide inclined field (Bx = 0.8,  By = 0.6) at different times, corresponding to the snapshots in Fig. 8.
In the text

Fig. 12. Plasma (solid black) and neutral (dashed) densities with height for ξn0 = 0.9 (red), 0.99 (blue), and 0.999 (green). The MHD density is given by the black line and the plasma component of the two-fluid simulation is identical to the MHD density profile. The vertical line indicate the βv point (where cs = vA) for the plasma species (and MHD) only. The coloured vertical lines show the βvt points using the total densities (ρn + ρp), that is, cst = vAt.

In the text

	Fig. 13. Long atmosphere used for the collisional coefficient simulations. Black line shows the plasma (and MHD) density. Red dashed line shows the neutral density for an initial neutral fraction of ξn = 0.999. The βv and βvt points are shown by the black and green vertical lines.
In the text

	Fig. 14. Snapshots of the parallel (solid) and perpendicular (dashed) velocities for the plasma (black) and the neutral (red) species for different collisional coefficients.
In the text

	Fig. 15. Perpendicular shock widths for different collisional coefficients. Corresponds to the snapshots in Fig. 14. Shock width is calculated as the distance between the maximum velocity, and the maximum absolute gradient of velocity.
In the text

	Fig. 16. Shock locations through time for different collisional coupling coefficients. Parallel (solid) and perpendicular (dashed) velocities and shown for the plasma (black) and neutral (red) species. The βv and βvt locations are shown as the green and blue line respectively. The graph is plotted until the perpendicular shock approaches the top of the domain.
In the text

	Fig. 17. Parallel (solid) and perpendicular (dashed) velocities in the plasma (black) and neutral (red) species for the vertical magnetic field (top row) and inclined field (bottom row) for different neutral fractions (ξn = 0.9, 0.99, 0.999 across the columns). The vertical lines indicate the βv (black) and βvt (green) points.
In the text

	Fig. B.1. Numerical diffusion tests for a slow-mode shock in non-stratified (top) and stratified (bottom) atmospheres. Top panel: MHD slow shock used in Hillier et al. (2016); bottom panel: MHD simulation in Fig. 4 at time t = 5.68. The numerical shock width in the stratified atmosphere is less then 0.02 pressure scale heights and considered negligible.
In the text

	Fig. C.1. Shock transitions for different collisional coefficients. Parallel (solid) and perpendicular (dashed) lines are shown for the plasma (black and green) and neutral (red) species. The black lines show the transitions relative to the isolated plasma, whereas the green lines show the transitions relative to the fully coupled system.
In the text