Stability of two-fluid partially-ionised slow-mode shock fronts

A magnetohydrodynamic (MHD) shock front can be unstable to the corrugation instability, which causes a perturbed shock front to become increasingly corrugated with time. An ideal MHD parallel shock (where the velocity and magnetic fields are aligned) is unconditionally unstable to the corrugation instability, whereas the ideal hydrodynamic (HD) counterpart is unconditionally stable. For a partially ionised medium (for example the solar chromosphere), both hydrodynamic and magnetohydrodynamic species coexist and the stability of the system has not been studied. In this paper, we perform numerical simulations of the corrugation instability in two-fluid partially-ionised shock fronts to investigate the stability conditions, and compare the results to HD and MHD simulations. Our simulations consist of an initially steady 2D parallel shock encountering a localised upstream density perturbation. In MHD, this perturbation results in an unstable shock front and the corrugation grows with time. We find that for the two-fluid simulation, the neutral species can act to stabilise the shock front. A parameter study is performed to analyse the conditions under which the shock front is stable and unstable. We find that for very weakly coupled or very strongly coupled partially-ionised system the shock front is unstable, as the system tends towards MHD. However, for a finite coupling, we find that the neutrals can stabilise the shock front, and produce new features including shock channels in the neutral species. We derive an equation that relates the stable wavelength range to the ion-neutral and neutral-ion coupling frequencies and the Mach number. Applying this relation to umbral flashes give an estimated range of stable wavelengths between 0.6 and 56 km.


INTRODUCTION
Magnetohydrodynamic (MHD) shock waves are a fundamental feature of astrophysical systems, occurring across a range of physical systems, for example, the solar atmosphere (Beckers & Tallant 1969;Hollweg 1982) and molecular clouds (Draine et al. 1983). The formation mechanism behind astrophysical shocks is as varied and results as a consequence of a number of physical processes, such as magnetic reconnection (Yamada et al. 2010;Petschek 1964) and wave steepening (Suematsu et al. 1982). For a hydrodynamic (HD) system, sonic shocks can exist as a transition from above to below the sound speed. In magnetohydrodynamics (MHD) there are three characteristic wave speeds (slow, Alfvén, fast) which leads to a wealth of possible shock transitions (Tidman & Krall 1971). Here we study the stability of slow mode shock fronts in partially ionised systems.
For a steady-state shock, the MHD shock jumps can be described in terms of the conservative quantities sufficiently upstream and downstream of the shock. It can be shown that MHD systems support a variety of shock transitions (broadly categorised as slow, intermediate, and fast) as well as contact and tangential discontinuies. If a non-zero component of magnetic field exists normal to the steady shock ( ≠ 0) then jumps in the tangential magnetic field imply jumps in the tangential velocity. As such, vorticity (∇ × v) can exist at an MHD shock front. For an MHD contact discontinuity, where there is zero mass flux across the interface, only jumps in the density are permitted. As such, there is no magnetic field jump to support a ★ E-mail: b.snow@exeter.ac.uk jump in tangential velocity hence vorticity cannot exist at a steadystate MHD contact discontinuity. In a hydrodynamic system, the magnetic field is zero and hence the shock jumps are easier to compute. With regards to vorticity, the HD system has the opposite behaviour to an MHD system; vorticity can be exist across a HD contact discontinuity, but not a HD shock (since a transverse velocity is no longer supported by a magnetic field, Hayes 1957). This has consequences for instabilities such as the Richtmyer-Meshkov instability which results in a hydrodynamic contact discontinuity being unstable, however the instability is suppressed in MHD (Wheatley et al. 2009). For the corrugation instability (where a shock front is perturbed), the opposite can be true where the instability grows in MHD but is suppressed in HD.
The stability of shocks to the corrugation instability depends on the type of shock. Fast MHD shocks are categorically stable provided the adiabatic index < 3 (Gardner & Kruskal 1964), whereas parallel MHD shocks (where the velocity and magnetic fields are parallel to the shock front) are always unstable to the corrugation instability (Stone & Edelman 1995). The stability of other types of shocks (e.g., switch-off) depends on the angle of the magnetic field (relative to the shock front) and the Alfvén Mach number (Stone & Edelman 1995;Édel'Man 1989;Lessen & Deshpande 1967). Here we focus on a slow-mode (parallel) shock that is categorically unstable for MHD, and stable for HD.
The stability of shock fronts to perturbations is important for compressible turbulent astrophysical systems, such as the interstellar medium (Elmegreen & Scalo 2004) and solar atmosphere (Reardon et al. 2008), where slow-mode shocks can encounter a non-uniform This may make slow-mode shocks more difficult to identify than fast mode shocks due to the corrugated shock front (Park & Ryu 2019). Similar results are present in reconnection studies that feature turbulence (Zenitani & Miyoshi 2011. In addition to MHD and HD systems, shocks regularly occur in warm plasma, where the medium is only partially ionised, such as the solar chromosphere and molecular clouds. These partially ionised mediums represent a challenging area of study, specifically for instabilities, where the general behaviour is neither MHD-like or HD-like. Partially-ionised mediums can be modelled using two-fluid (neutral, ion+electron) equations that allow for coupling and decoupling of the two species, and non-MHD-like behaviour. Two-fluid interactions in interface instabilities (where the initial magnetic field is parallel to the jump) can suppress small-scale features (Popescu Braileanu et al. 2020) and allow cross-field transport of momentum and energy (Hillier 2019).
The stability of two-fluid shocks to the corrugation instability has not been studied and likely depends on the level of coupling of the plasma and neutral species. A consequence of partial ionisation on shocks is that the shock wave has a finite-width that is determined by the physical parameters of the system. Hillier et al. (2016) extensively studied the two-fluid effects in a switch-off slow-mode shock. Within the finite width of the shock additional shock transitions can occur due to the species decoupling (Snow & Hillier 2019). One can hypothesise that under the extreme conditions of the finite width being zero or infinity, the plasma should reduce to MHD-like behaviour. However, in the more realistic finitely-coupled regimes the behaviour is undetermined, which forms the basis of the study presented here.
In this paper, we study the stability of two-fluid slow-mode shock fronts using 2D numerical simulations. We investigate the consequences and stability of the shock front for different levels of collisional coupling and different ionisation fractions. We find that the shock-front can be stable or unstable depending on the wavelength of the perturbation relative to the finite width of the shock. This may have consequences for observations by giving an expected range at which partially-ionised slow-mode shocks can be stable.

CORRUGATION INSTABILITY SCHEMATIC
In this paper, we focus on specifically the parallel slow-mode shock, where the flow is aligned with the magnetic field either side of the shock. This type of shock is unconditionally unstable in the MHD model (Stone & Edelman 1995), and stable in the hydrodynamic model, to the corrugation instability. In this section, we present the basic schematic for the HD stability and MHD instability for a 2D shock front encountering an upstream density perturbation, and make conjectures about the potential stability or instability expected for a partially ionised system.

MHD
An ideal magnetohydrodynamic (MHD) shock is unstable to the corrugation instability under certain conditions. For a parallel MHD shock, where the velocity and magnetic fields are parallel to the shock front either side of the shock (i.e., v = [ , 0, 0], B = [ , 0, 0] in the shock frame), the shock front is always unstable to the corrugation instability (Stone & Edelman 1995). The steady-state shock jump conditions can be written as where the above notation relates conserved quantities upstream (superscript ) and downstream (superscript ) of the shock front as: for any conserved quantity . From these equations (particularly Equation 3) it can be seen that a jump in the tangential component of the velocity can be supported across a steady-state shock front by a jump in the tangential magnetic field component. When the shock front corrugates it generates a tangential magnetic field component and hence a tangential velocity. This velocity is continuous across the interface and hence the flow is directed behind the peaks, see Figure1. The flow pattern increases the pressure behind the peaks and decreases it behind the troughs. As such, the peaks are pushed up and the troughs are pushed down, increasing the magnitude of the corrugation and leading to growth of the instability. This is shown by the schematic in Figure 1.

Hydrodynamics
An ideal hydrodynamic (HD) shock is stable to the corrugation instability. In a HD system, tangential velocity is continuous across the shock front. Therefore, at the distorted shock front the flow is refracted away from the 'peaks'. The flow then enhances the pressure behind the 'troughs', and decreases the pressure behind the 'peaks', see Figure 2. As such, the pressure provides a restoring force whereby the peaks are pushed down and the troughs are pushed up. The pressure restoring force therefore acts to mitigate the corrugation of the shock front and the system stabilises, returning to a flat shock front. A schematic of this is shown in Figure 2 where a high pressure region forms behind the troughs.

Interpolating these results to partially-ionised plasma (PIP)
In the two fluid partially ionised case, the schematic is not as straightforward. The ionised species behaves like an MHD system, whereas the neutral species is HD. As such, one may consider that the isolated neutral species tends towards stability, whereas the isolated plasma species tends towards being unstable. The stability of the bulk (plasma+neutral) PIP shock depends on the balance of these forces, and the timescales on which the species interact.
The coupling between the ionised and neutral species is governed by the collisional coefficient (Hillier et al. 2016) which is discussed further in Section 3. In the the extreme of = 0, the system is fully decoupled, and the ionised species will be unstable. At the other extreme of = ∞ the system is fully coupled and MHD-like, hence one would also expect the system to be unstable. In the finitely coupled regime, interactions can allow the system to be unstable or stable. The finitely coupled regime is the focus of our paper.
Similarly, one has to consider the perturbation wavelength relative to the finite width of the shock. Consider a fixed size perturbation of wavelength that is parallel to the shock front, at different levels of coupling. For a infinitely coupled system, the finite width of the shock approaches zero and hence the perturbation wavelength is much larger than the finite width of the shock. As such, one may expect that the system will behave similar to a bulk-MHD model. At the other extreme of coupling tending to zero, the finite width of the shock becomes large. The upstream perturbation in this case will encounter the plasma shock front first and therefore also behave like an MHD simulation. Between these two extremes, we expect the system to be finitely coupled. As such, the stability of a partiallyionised shock can be connected to the wavelength of the perturbation relative to the finite shock width.

METHODS
The two-fluid equations governing the behaviour of a neutral species (subscript n) and a charge-neutral electron+ion species (subscript p) are: for density , energy , pressure , velocity v and magnetic field B. The adiabatic index = 5/3 and is constant. The plasma and neutral species are coupled through thermal collisions which is controlled through the coupling coefficient , where , are the neutral and plasma temperatures respectively. The constant 0 can be specified and governs the coupling between the two species. We study a hydrogen-only system.
We use a 4th-order central difference solver with artificial viscosity to limit the numerical oscillations around the shock front. It was found this solver is far less susceptible to the numerical Carbunkle instability (which forms due to a curved feature forming on a Cartesian grid, Elling 2009) than the first order HLLD solver used in previous shock studies with the (PIP) code (Snow & Hillier 2019, 2021. The 4th order scheme in the (PIP) code has been used successfully in prior studies of partially ionised plasma dynamics (Hillier 2019;Murtas et al. 2021).

Initial conditions
The initial conditions used in this paper are in the shock frame, where the shock is stationary at = 0, with the inflow coming from the upstream conditions ( > 0), and the post-shock region downstream ( < 0). This has the advantage of allowing us to increase the resolution of the system since no space is required for the shock to propagate into.
A sonic Mach 2 parallel shock is specified analytically, see Appendix A, where 'parallel' refers to the initial flow being aligned with the magnetic field. Here, the shock is in the -direction and = = 0 initially and corresponds to a slow-mode shock. The pressure and density are specified such that we have a thermal equilibrium and the initial plasma-= 0.1 in the upstream medium.
The corrugation is initiated by placing a small pseudo-random density perturbation upstream of the shock for the form where 0 , 1 , 0 , 1 give the horizontal and vertical locations for the perturbation. 10 wavelengths are imposed parallel to the shock front (i.e., in the direction), with random amplitude and phase (determined by ). The perturbation is placed sufficiently upstream of the shock that the shock front has time to relax to its numerical equilibrium before encountering the density enhancement. Note that the initial shock front is uniform in the −direction so evolves like a 1D system before it encounters the density perturbation. We initially set the perturbation between 0 = 10, 1 = 20 and 0 = 0, 1 = 10 covering the full vertical extent of our box and horizontally located upstream of the shock front. The same random seed (and hence the same perturbation) is used in all simulations. The initial conditions are shown in Figure 3 along with a 1D plot of the density perturbation. The perturbation is placed sufficiently far upstream to allow the shock front to relax slightly to its numerical equilibrium before the shock front encounters the density enhancement. For the two-fluid case considered here (where the species are thermally coupled), the plasma conditions sufficiently upstream and downstream of the shock and described by the MHD shock jump equations (Snow & Hillier 2019, 2021. As such, the bulk (plasma+neutral) density, pressure, velocity and magnetic field jumps are the same as the MHD case, see Appendix A. The partial pressures and densities are calculated using the neutral fraction and the bulk pressure and density. The neutral and plasma parameters are therefore defined as: where a subscript denotes the MHD value and and are the neutral and plasma fractions respectively. In this formation, the MHD MHD (a=¥) simulation is effectively the infinitely coupled case where the system behaves like a bulk (plasma+neutral) fluid. At the other extreme of = 0 the fluids are completely decoupled and the plasma will behave MHD-like, free of any interactions with the neutral species. The shock front is allowed to numerically stabilise before encountering the upstream density perturbation.

Boundary conditions
Upper and lower ( = 0, 10) boundary conditions are set to periodic. The left and right boundaries ( = −25, 25) are specified as damping layers using a similar formulation as Felipe et al. (2010). Several waves are generated from the initial shock wave relaxing and the interactions with the density perturbation. This damping layer effectively removes these unwanted waves from the system.
where U is the instantaneous evolved vector of conserved variables and U 0 is the equilibrium quantity. The initial quantities U 0 are determined analytically using the shock jump conditions. The perturbation to be damped is then U − U 0 which is gradually damped over = 20 grid cells.
is the value at the end of the damping region. The factor 0 scales the maximum damping to be 0.3. We find that this boundary condition is sufficient to remove unwanted perturbations from system and prevents reflections from the boundary. The boundaries are also placed far from the shock front to further minimise their impact on the results.

MHD simulation
For the MHD system, the shock front encounters the density perturbation, corrugates, then the shock front perturbation grows with time, see Figure 4. The initial perturbation of the shock front allows a non-zero field to exist at the shock interface leading to growth of the instability, as described in Section 2.
In this paper, we quantify the corrugation using the integrated enstrophy (vorticity squared) around the shock front as: where is the location of the shock front. The integral limits capture the full extend of the domain in the -direction, and the horizontal  . Initial shock for the MHD (black) and PIP (red neutrals, blue plasma) at = 6. The snapshot is taken before the shock front collides with the density enhancement. Small numerical oscillations exist at the shock front as a consequence of the numerical scheme.
window is chosen to capture the full extent of the corrugation in the -direction. The vorticity is non-zero at the distorted shock front, and zero elsewhere. As such, the integrated enstrophy allows us to determine if the corrugation is growing or decaying with time.
For the MHD model, Figure 5, the integrated enstrophy around the shock is close to zero initially, as it should be for a uniform shock front. When the shock front encounters the density enhancement, there is a sudden rise in the enstrophy as the shock front corrugates. As time increase, the corrugation continues to grow, as seen in Figure 4, and hence the integrated enstrophy grows. As time increases towards infinity the growth of the instability will saturate, as seen in Stone & Edelman (1995), however this has not yet occurred in time frame studied here.
A 1D shock front can be extracted from the 2D simulation by taking the displacement at each height , with the location of the shock determined by the maximum density gradient in the -direction. Plotting the shock front through time shows the growth of the corrugation, as shown in Figure 6. The initial disturbance warps the shock front. As time increases, the displacement grows across all frequencies. This shows the expected result that the MHD simulation is unstable to the corrugation instability.

PIP simulation
Here we present an initial two-fluid PIP simulation where the neutral fraction is set as = 0.9 and a coupling coefficient of 0 = 1, giving a mostly neutral medium that is relatively weakly coupled on the time scales considered here. This choice of coupling coefficent and initial conditions means that on average, upstream collisions occur on a time scale of unity. The bulk (neutral+plasma) density and pressure are the same as in the MHD case, see Appendix A. For a two-fluid system, the bulk (neutral+plasma) shock jumps reduce to the MHD shock jumps sufficiently upstream and downstream of the shock (Snow & Hillier 2019). However, within the finite-width of the shock, drift velocities form between the two species and substructure exists. The role of this substructure on the stability of the shock front is investigated here. Figure 4 shows the plasma and neutral densities through time for this fiducial PIP simulation. Initially the shock front is uniform in the −direction, as with the MHD case, however the two-fluid effects produce a finite width to the shock which is highlighted in Figure 7 (the finite width of two-fluid shocks is discussed in detail in Hillier et al. (2016)). When the shock front encounters the density enhancement, corrugations occur in both the plasma and neutral shock fronts. These corrugations initially grow in the plasma and decay in the neutrals. Finally, the distortions decay in both species and the shock front stabilises. This can be seen in Figure 5, where there is a rise in corrugation, followed by a decay towards a stable shock front.
For this simulation, the plasma shock is at the leading edge, followed by the neutral shock, see Figure 7. As such, the plasma species encounters the density perturbation first, and the instability grows in an MHD-like manner initially. At this time, there is also a slight corrugation of the neutral shock front from the encountered perturbation. However, as time advances, the displacement of the shock front decreases in both the neutral and plasma species and the shock front stabilises.
The schematic for flow entering the initially distorted shock is given in Figure 8 and can be compared to the schematics for MHD ( Figure 1) and HD (Figure 2). The flow entering the plasma shock is slightly corrugated in the same way as the MHD case. The magnitude is far less in the PIP case than the MHD case due to the resistance provided by the neutral species. In the neutral species, the flow is directed towards the troughs, creating a high pressure region that suppresses the instability. In this regime with a mostly neutral medium ( = 0.9) and a weakly coupled system ( 0 = 1) the interactions between the plasma and neutral species are sufficient to suppress the corrugation instability and stabilise the shock front. This is further shown in the evolution of the shock front displacement with time, Figure 6, where the plasma shock front distorts initially, then stabilises. The neutral shock front distorts slightly after encountering the density enhancement but is quickly stabilised. Physically, the neutral species is trying to stabilise the shock front, whereas the plasma species is trying to go unstable. Here the overall physics is dominated neutral species since medium is mostly neutral (hence the plasma-neutral coupling is prevalent) and the perturbation wavelengths in the −direction are equally affected by the coupling (this is discussed further in Section 5).
The finite width of the shock front here is ≈ 0.6, which is smaller than parallel perturbation wavelengths that are in the range of = 1 − 10 nondimensional units. The ratio of the shock width to the perturbation wavelength should govern the behaviour of the system. For a perturbation much larger than the finite width, the system should behave like a fully coupled MHD simulation since the finite-shock width becomes negligible and overall behaviour is determined by bulk properties. For a perturbation much smaller than the finite width, the simulation should behave like a fully decoupled simulation, also MHD like, since on short time scales the perturbation will only encounter the plasma shock front. Between these two limits, the finite coupling leads to non-MHD-like behaviour.

Coupling coefficient
Changing the coupling coefficient changes the finite width of the shock, with larger coupling leading to a narrower finite width (Hillier et al. 2016). For the two-fluid = 2 shock investigated in this paper, with an upstream plasma-= 0.1, the finite width of the shock is a function of the coupling coefficient, as shown in Figure 9 from a series of 1D numerical simulations. An analytic estimate is also shown, which is calculated as the difference in propagation speeds of the isolated fluids, divided by the coupling frequency, i.e., for the plasma and neutral sound speeds and . The analytical and numerical results for the shock width pair up very well, and scale well in the coupling regime investigated here. Also, changing the neutral fraction does not greatly affect the shock width and this estimate remains a good fit to the 1D simulation data.
At the limit of = ∞ and = 0, the system is MHD with the fluid acting as a bulk (neutral+plasma) system or a completely decoupled system respectively. Between these limits, we expect twofluid effects to become important, as seen in Section 4.2. Here we investigate the consequences of changing the coupling coefficient on the stability of the shock front. Figure 10 shows a time evolution of the shock front for different levels of coupling, and two MHD simulations that represent the infinitely coupled and fully decoupled extremes. Figure 11 shows the corresponding integrated entrophy around the shock front through MHD (a=0) time for these cases. One can see that the MHD simulations, the corrugation grows with time and persist throughout the simulation. The growth of the = 0 case saturates near the end of the simulation, whereas the = ∞ case continues to grow. In both cases however, the shock front has become unstable and the corrugation is a consistent feature of the system. In the PIP simulations, the initial perturbation distorts the shock front (as with the MHD cases) however the perturbation is seen to stabilise to varying degrees depending on the coupling coefficient. For the = 1 case, all frequencies are damped at roughly the same rate, with the shock front tending to a flat shock front. For = 10 the moderate coupling has a severe effect in rapidly damping the perturbation, quickly leading to a stable/flat shock front. At = 100 the perturbation grows and has saturated, however the smaller frequencies have been damped from the system and only larger frequencies exist.
The integrated enstrophy around the shock is shown in Figure 11a where the time zero corresponds to the shock encountering the perturbation. The corresponding time derivative at time = 10 is shown in Figure 11b. The MHD ( = 0) saturates during the simulation and hence the time derivative is very close to zero, and in fact slightly negative. It can be seen from Figure 10 that the perturbed shock front remains distorted with time and the growth rate being close to zero is due to saturation of the instability. For a weakly coupled system ( = 0.1) the system is unstable and growing, however the rate is very small implying that it is approaching saturation. At the end time, the simulation is comparable to the MHD ( = 0) case. The perturbation width for = 0.1 is approximately the same size as the finite width of the shock. As the coupling increases ( = 1, 10), the system becomes stable as the neutrals start to influence the plasma. This is seen clearly in Figure 10 where these cases approach an undisturbed state at late times. At strong coupling ( = 100), the system begins to act like an ensemble bulk MHD ( = ∞) system and the perturbation persists through time. From Figure 11b it can be seen that the growth rate here is very close to zero and likely saturated for = 100. The structure of the perturbed shock front in this case is also slightly different to the MHD = ∞ case; only the low frequency components are present in the perturbation, and it is much smoother than the MHD case.
The wavelengths that persist in the shock front vary for different coupling coefficients, see Figure 12 which shows the normalised shock front displacement. For weak coupling ( = 1) the wavelengths are all damped fairly uniformly, as seen in Figure 6. For moderate coupling ( = 10) the system stabilises very rapidly and only the longest wavelengths persists. This is shown in Figure 12 where only the fundamental mode is present. For strong coupling ( = 100), the shorter scales are damped, and several longer scales persist. It is known that there are several scales on which the neutral and plasma species interact (Hillier 2019). On the largest scales, the neutrals and plasma species are coupled. On intermediate length scales the neutrals are decoupled from the plasma, and on the smallest length scales both species are decoupled from each other. These scales are determined using the collision coupling frequencies: where and are the ion-neutral and neutral-ion collisional frequencies respectively. As increases, the coupling frequencies increase. One can separate the coupling with wavelength into three frequency bands. For low frequencies, the neutrals and plasma are well coupled. For intermediate frequencies, the neutrals decouple from the plasma (and hence the magnetic field), however the plasma is still coupled to the neutrals. For high frequencies, both species decouple. For the simulations performed here, increasing the coupling coefficient effectively decreases the wavelength at which the neutrals decouple from the plasma. This can be used to explain the damping of short wavelengths in the = 100 case. The approximate coupling frequency at which the system becomes unstable can be estimated by equating the coupling frequency to the frequency of the wavelength. For the lower limit, the frequency at which the plasma decouples from the neutrals for the largest wavelength is ,min Similarly, the upper limit is provided by the frequency at which the neutral decoupled from the plasma using the smallest wavelength: The velocity jump in our simulation is | − | ≈ 1.1, the upstream densities are = 0.1, = 0.9, and our wavelengths are in the range ∈ Z : ∈ [1, 10], with max = 10 and min = 1. As such, we estimate that the upper and lower limits on stability of our system are ,min ≈ 0.754 and ,max ≈ 69.1, which are shown by the dashed lines on Figure 11. One can see that these approximations are a reasonably good fit to the simulation results, and all simulations between ,min < < ,max are stable to the corrugation instability.
We note that a time scale argument exists for the weakly coupled systems. On the time scales considered here the corrugation evolves and reaches saturation, however, as time advances, the plasma will experience more collisions from the neutral species. As such, as time tends to infinity, the saturated, corrugated PIP shocks may further evolve due to a gradual influence of the neutral species.

Ionisation fraction
Simulations thus far used a neutral fraction of = 0.9 meaning that the medium is mostly neutral particles. Here we investigate changing the neutral fraction by looking at = 0.99, 0.5, 0.1. The = 0.1 case is mostly ionised. The = 0.99 case is dominated by the neutral species. The = 0.5 is an interesting case where the system is equally ionised and neutral. A time series showing the plasma and neutral densities at different times is shown in Figure 13. The = 0.99 case behaves fairly similar to the = 0.9 case studied previously. Both the plasma and neutral shock fronts corrugate slightly when they encounters the upstream density perturbation. As time advances, the shock front stabilises. Since this simulation is dominated by the neutral species, the stabilising forces from the neutral species are stronger than the destabilising forces of the plasma. The corrugation of the plasma shock front increases between times = 2, 5. However the neutral coupling then leads to a stabilisation of the shock front. This is comparable to the results of the = 0.9 case in Section 4.2.
For a neutral fraction of = 0.5 the system is equally neutral and ionised. The evolution of the shock front through time is shown in Figure 13. The shock front is unstable to the corrugation instability here and the shock front distortion grows with time. Interestingly, this x n =0.99 simulation forms far longer fingers than has been seen in the MHD or PIP simulations in Section 4.2 Between the elongated fingers, the neutral species also develops shock channels, as seen by the sharp density structures.
In the MHD case, the flow is directed behind the peaks (see Figure  1). For the PIP = 0.5 case, the neutral species resists this pattern and the plasma flow is directed towards the troughs instead, see Figure 14. As time advances, this forces the troughs to extend further back and creates the elongated fingers seen in the Figure 13. In the neutral species, similar behaviour occurs where the flow is directed towards the troughs, resulting in an increased neutral pressure and the formation of shocks due to the converging flow, Figure 14. When the plasma can exert a strong influence on the neutrals, these shock channels can form as flow pattern of the neutral species is dominated by the plasma motion.

The
= 0.1 case is dominated by the plasma species, and, as such, the corrugation of the plasma shock front grows, similar to the MHD model. The behaviour of the neutral species is mostly determined through the coupling with the plasma species. The result is that the neutral flow is channelled towards the troughs, as with the = 0.5 case. However, here the neutral species constitutes far less of the bulk medium and hence there is very little feedback from the neutrals to the plasma.

Perturbation length scale
One can re-frame the effect of different collisional coefficients in terms of the length scale relative to finite width of the shock. The initial case studied in Section 4.2 has a finite shock width of approximately = 0.68 and an initial parallel perturbation wavelengths in the range = [1, 10], therefore the wavelength of the perturbation is larger than the finite width of the shock. The perturbation wavelength relative to the finite width of the shock is an important parameter. The simulations in Section 5 showed that when the perturbation wavelength was between 10 and 1000 times larger than the finite width, the neutral interactions stabilised the shock front.
The finite-width of the parallel slow shock considered here can be estimated using Equation 29. Rewriting this in terms of bulk fluid properties: The bulk sound speed in the chromosphere is approximately 8 km/s. The ion-neutral collisional frequency in the chromosphere varies in the range of 10 3 − 10 6 s −1 (Popescu Braileanu et al. 2019). This leads to a finite-width of approximately ≈ 4 m. In Section 5 it is shown that perturbations 1000 times larger than finite width were unstable. Therefore, we expect that in the partially ionised solar chromosphere, slow-mode shock fronts will be stable to perturbations on the order of approximately 4 − 4000 m. Physically, this means that the perturbations that we are able to observe should result in unstable shock fronts. Below our observational limits, the partially ionised shocks should stabilise. The analytical approximation for the stability range of a partially ionised parallel shock (given by Equation 32-33) is a reasonable approximation to the simulation results. This formulation can be written for a parallel slow-mode shock of an arbitrary sonic Mach number as: ,max for compressible ratio = / , where > 1 for a shock. Using the relationship between and for a parallel slow-mode shock: we can rewrite Equations 39,40 as ,min ,max Therefore, given the upstream plasma and neutral densities, the upstream sound speed, the shock Mach number, and the range of perturbation wavelengths upstream of the shock, one can estimate the coupling frequencies at which the the shock may become unstable. These coupling frequencies can then be compared to the expected range of frequencies for the medium (e.g., Popescu Braileanu et al.

2019)
The results can be used to approximate the stability of shocks in a solar sunspot, where the magnetic field is fairly straight and propagating shocks are regularly observed in the form of umbral flashes (Beckers & Tallant 1969), which have average lifetimes of around 44.2 seconds (Nelson et al. 2017). Umbral flashes have Mach numbers between = 1 and = 1.7 (Anan et al. 2019) and using the coupling frequencies in the upper chromosphere (e.g., see Figure  2 in Popescu Braileanu et al. 2019) of = 3 × 10 2 , = 3 s −1 with a typical sound speed of = 8 km/s gives a stable range of wavelengths between approximately 0.6 and 56 km (for M=1.7).
For the parallel slow-mode shock investigated here, the shock width has very little dependence on the ionisation fraction, see Figure  9. However, it is known that for switch-off shocks, the finite width has a dependence on the neutral fraction and plasma− (Hillier et al. 2016). Also, the physical width of the shocks studied in this paper is much smaller than the > 300 km shocks in Snow & Hillier (2020) As such, one would expect that the stability range changes for different shock types. Further work is needed to calculate this stable range for different types of partially-ionised shocks.

Slow-mode shocks in turbulence
In MHD simulations of turbulence, fast-mode shocks are identified more readily than slow-mode shocks (Park & Ryu 2019). The potential reason for this is that the slow-mode shocks are unstable to the corrugation instability, hence the shock front deforms and is difficult to identify. However, analysing the MHD simulation of the corrugation instability by using an automated shock detection method (based on the SHOCKFIND algorithm developed by Lehmann et al. 2016) we see that the slow-mode shocks are still detected over the vast majority of the shock front, see Figure 15. In face, due to the increased length of the shock, approximately 10 times more points are detected in the corrugated shock front at time = 30 than in the initial stable shock front. This implies that the corrugation instability should lead to a greater detection of slow-mode shocks, contrary to the claim of  Figure 15. Detected pixels that satisfy a slow-mode shock transition in the reference MHD simulation at time = 30. Note that the detected shock pixels cover the majority of the shock front. Park & Ryu (2019). Further study is needed to determine the effect of the corrugation instability on the presence of shocks for different levels of corrugation and Mach numbers.
If the corrugation instability is reducing the detection of slow mode shocks (as suggested by Park & Ryu 2019), then in a partially ionised system, we may detect statistically more slow-mode shocks in a two-fluid turbulence than MHD turbulence. Here we show that for a partially ionised medium, the neutral interactions can stabilise slow-mode shocks. As such, the corrugation of the shock front is a transient feature and shocks should damp the perturbation. Stone & Edelman (1995) show that a wider range of shock types are unstable in 3D MHD, since perturbation are allowed to grow in the out-of-plane direction and MHD shocks that are stable in 2D can be unstable in 3D. Here we see that for an unstable 2D shock in a partially-ionised plasma, the neutral species can act to stabilise the system. In 3D we would expect that the neutral species could act to stabilise the corrugation growth in the out-of-plane direction, resulting in a stable 3D partially ionised shock. A future study will investigate the stability of 3D partially ionised shocks.

CONCLUSIONS
In this paper we have investigated the stability of slow-mode shocks in partially ionised plasma to the corrugation instability. As the shock front encounters a density perturbation, the shock front distorts leading to the corrugation instability. For an MHD system, a parallel slow-mode shock is always unstable, whereas a HD shock is usually stable. For the partially-ionised cases presented here, we show that the instability can be stable or unstable depending on the system parameters.
In the reference case ( = 1, = 0.9), the corrugation grows in the plasma species initially, then the coupling with the neutral species leads to stabilisation of the shock front. The stability of the shock front is dependent on the coupling coefficient. At large coupling ( >> 100) the system is unstable to the corrugation instability and behaves like a bulk MHD system. At the other extreme ( < 0.1), the system is again unstable and MHD-like, behaving like a fully decoupled plasma. Between these limits (0.1 < < 100) the neutral interactions with the plasma species lead to a stabilisation of the shock front to the corrugation instability. The stability range can be rewritten in therms of the perturbation width relative to the finite-width of the shock, and determines the stability of the shock front. When the perturbation width is on the same order (or smaller) than the finite-width of the shock, the system is unstable and MHD-like. Similarly, when the perturbation is 1000 times larger than the finite-width, the system is again unstable and MHD-like. Between these limits, the two-fluid interactions become important and the neutral species can stabilise the shock front.
In conclusion, two-fluid interactions can lead to stabilisation of MHD parallel slow-mode shock fronts to the corrugation instability. Further work is needed to determine the stability criteria of different shock types.

APPENDIX A: PARALLEL SHOCKS
A parallel shock is where both the velocity and magnetic field vectors are parallel shock front in both the upstream and downstream states, i.e., v = ( , 0, 0), B = ( , 0, 0). Using these, the shock equations reduce to: [ ] = 0 (A1) The compression across the shock can be defined as / = . This leads to the relation that / = 1/ (from equation A1). Introducing Finally, can be determined using Equation A3 and some algebra as: Therefore, for given upstream properties, the downstream conditions can be easily calculated. In these simulations, the system is normalised such that the upstream sound speed is one by setting = 1, = 1/ . In the shock frame, our upstream and downstream states are given in table A1. This paper has been typeset from a T E X/L A T E X file prepared by the author.