Shocks and instabilities in the partially ionised solar atmosphere

The low solar atmosphere is composed of mostly neutral particles, but the importance of the magnetic field for understanding observed dynamics means that interactions between charged and neutral particles play a very important role in controlling the macroscopic fluid motions. As the exchange of momentum between fluids, essential for the neutral fluid to effectively feel the Lorentz force, is through collisional interactions, the relative timescale of these interactions to the dynamic timescale determines whether a single-fluid model or, when the dynamic frequency is higher, the more detailed two-fluid model is the more appropriate. However, as many MHD phenomena fundamentally contain multi-time-scale processes, even large-scale, long-timescale motions can have an important physical contribution from two-fluid processes. In this review we will focus on two-fluid models, looking in detail at two areas where the multi-time-scale nature of the solar atmosphere means that two-fluid physics can easily develop: shock-waves and instabilities. We then connect these ideas to observations attempting to diagnose two-fluid behaviour in the solar atmosphere, suggesting some ways forward to bring observations and simulations closer together.


The partially ionised solar atmosphere
Observations of the cool material in the lower solar atmosphere show a wide range of both interesting (from a perspective of theoretically understanding them) and important (in terms of their role in mass transfer and heating of the solar atmosphere) dynamic phenomena. At the lowest level of the solar atmosphere, the solar photosphere, observations show cell-like structures known as granules. These are turbulent convective cells carrying heat from the solar interior into the atmosphere. As we move up into the solar chromosphere we find a change in the dynamics as we move from the fluid dominated layers into those where the magnetic field dominates, characterised by the plasma b (ratio of gas to magnetic pressure) dropping below 1 (Gary, 2001). Observations show many dynamic features in this layer of the atmosphere including jets (e.g. Nishizuka et al., 2011;Singh et al., 2011) and spicules (e.g. Pereira et al., 2014), which are driven by or excite a wide range of physical processes including waves (e.g. Morton et al., 2014;Okamoto and De Pontieu, 2011), shocks (e.g. Houston et al., 2018;Houston et al., 2020) and magnetic reconnection (e.g. Nishizuka et al., 2008). We can find chromospheric material even high up in the solar atmosphere, suspended in the solar corona as prominences/filaments. Observations of these fascinating structures reveal waves (e.g. Okamoto et al., 2007;Hillier et al., 2013), instabilities (Berger et al., 2008;Berger et al., 2010;Berger et al., 2011;Berger et al., 2017;Hillier and Polito, 2018), reconnection (Hillier and Polito, 2021) and turbulence (Leonardis et al., 2012;Hillier et al., 2017). Overall there is a wide range of dynamic phenomena, all connecting with fundamental MHD theoretical concepts. An example of a prominence and many of the chromospheric dynamics discussed above is presented in Fig. 1.
There are multiple aspects that physically connect the different dynamics observed in these different layers of the solar atmosphere. The similar temperature of the material, in the range of a few thousand to 10 4 K, allows all these layers to be treated as the cool phase of the solar atmosphere. And being cool has a very important consequence. The (relatively) low temperatures of the plasma found in these layers of the atmosphere, along with the local densities, is insufficient to fully ionise the material. This results in the plasma in these layers being classed as partially ionised plasma, as a significant proportion of the species that compose it are neutral (Khomenko et al., 2014a).
Another key aspect that connects these dynamics is the importance of magnetic fields for either driving the dynamics, e.g. in launching chromospheric jets, or as a conduit for energy to be transported from lower levels to the atmosphere to higher regions. However, neutral species don't naturally feel the Lorentz force, so the fact that these lower regions of the solar atmosphere are partially ionised means that there has to be a physical process coupling the neutral material to the magnetic field. We will present how this physically happens, and its consequences for some dynamics observed in the solar atmosphere in the subsequent sections.
In this review article, we will present some fundamental ideas behind the modelling of partially ionised plasmas, and that allow us to understand where this will lead to differences in dynamics from a fully ionised MHD approximation. We will focus purely on advances in two-fluid modelling (see the next section for more details) with a par-ticular focus on multi-scale dynamics. The particular examples of this we focus on are shocks (see Section 3) and and instabilities (see Section 4).

Dynamics in partially ionised plasmas
With the plasma found in the lower solar atmosphere being partially ionised, and the connection between this plasma and the magnetic field of high importance to understand a myriad of observed phenomena, the question is: how does a plasma that is predominantly neutral couple with the magnetic field. For the solar atmosphere, the answer to this question comes through collision-like interactions between charged species (which feel the Lorentz force of the magnetic field) and the neutral particles. These might be hard-sphere collisions or through charge exchange (e.g. Vranjes and Krstic, 2013;Meier and Shumlak, 2012). There are other processes that couple the fluids, with ionisation and recombination an important one for solar plasma due to the role they have in determining the ionisation degree of the different regions of the solar atmosphere. However, the strongest coupling comes through collision processes (e.g. Vranjes and Krstic, 2013;Nó brega-Siverio et al., 2020a).
A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 sion (Braginskii, 1965). As the charged plasma component of a partially ionised plasma would feel fundamentally different forces than that of the neutral fluid (i.e. whether the fluid feels the Lorentz force or not), one sensible way to model this system is to look at a two-fluid system where the charged species are modelled as a charge-neutral plasma and the neutral species as a separate fluid with the two fluids coupled through collisions and ionisation/recombination (e.g. Khomenko, 2020). The equations for this approximation of the system (assuming a purely hydrogen plasma and including only the collisional coupling terms) are given as @q n @t þ r Á ðq n v n Þ ¼ 0; ð1Þ @ @t ðq n v n Þ þ r Á ðq n v n v n þ P n IÞ À q n g ¼ Àa c q n q p ðv n À v p Þ; ð2Þ @en @t þ r Á v n ðe n þ P n Þ ½ Àq n v n Á g ¼ Àa c q n q p 1 2 ðv 2 n À v 2 p Þ þ 1 cÀ1 P n q n À 1 2 P p q p h i ; ð3Þ for the neutral fluid and @q p @t þ r Á ðq p v p Þ ¼ 0; ð5Þ @ @t ðq p v p Þ þ r Á q p v p v p þ P p I À BB þ B 2 2 I À q p g ¼ a c q n q p ðv n À v p Þ; ð6Þ @ @t e p þ B 2 2 þ r Á v p ðe p þ P p Þ À ðv p Â BÞ Â B Â Ã Àq p v p Á g ¼ a c q n q p 1 2 ðv 2 n À v 2 p Þ þ 1 cÀ1 P n q n À 1 2 P p q p h i ; ð7Þ @B @t À r Â ðv p Â BÞ ¼ 0; ð8Þ for the plasma. The variables q; P ; B; v; and e are the density, pressure, magnetic field, velocity field, and internal energy respectively, and c is the adiabatic index. Note we have performed the simplification of B= ffiffiffiffi ffi l 0 p ! B in these equations. Here the subscripts n and p refer to the neutral and plasma fluid respectively, and g is the gravitational acceleration. We have assumed that these are ideal gases with P n ¼ R g q n T n and P p ¼ 2R g q p T p , with R g the gas constant and T n and T p the neutral and plasma temperatures. General information of derivation of such equations can be found in, for example, Meier and Shumlak (2012); Khomenko et al. (2014a).
Here we have assumed that the different charges can be treated as a single fluid (a key assumption in formulating the MHD equations). However, there are situations where treating the different charges (positively charged ions and negatively charged electrons) as different fluids becomes appropriate. The key criteria that need to be satisfied to assume a single, charge-neutral plasma fluid are that the phenomena studied are sufficiently large scale (i.e. larger than the Larmor radius, Debye length and electron mean-free-path) and low frequency (frequency smaller than the plasma and electron gyro-frequency). We have also assumed they are slow (sub-relativistic to allow, for example, the displacement current to be neglected). In this article, we will assume that we are discussing dynamics that satisfy these conditions.
The terms on the RHS of the above equations are the collisional coupling terms. In the momentum equations these are terms that transport momentum between the fluids due to collisions and this is linearly proportional to the velocity drift (v D ¼ v n À v p ). In the energy equation there is a term that gives both the transfer of kinetic energy and frictional heating of the system, and a term in the temperature difference that forces the fluids to relax to a thermal equilibrium. These momentum and energy transfer terms have a prefactor including the parameter a c written so they take the form with m in and m ni are the collision frequencies of ions onto neutrals and neutrals onto ions, respectively. Note here there is the assumption that the important collisions in the transfer of momentum and energy between the fluids are those between ions and neutrals and not electrons and neutrals. The form of a c is approximately given by for both hard sphere collisions (Draine, 1986) and charge exchange (Zank et al., 2018), where with T n and T p the neutral and plasma temperatures, m n and m i the mass of the neutral and ion particles involved in the collisions and k B Boltzmann's constant. Though it is common to ignore the term that allows collisions to be driven by velocity drifts, it can become important in highly dynamic systems (e.g. Murtas et al., 2021). The collisional cross-section (R in ) is generally larger for charge exchange interactions than hard-sphere collisions but these crosssections can vary greatly with the collisional velocity (e.g. Vranjes and Krstic, 2013).
In this review article we will focus on the two-fluid approximation and results from these calculation. There have been a range of studies looking at solar dynamics A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 using this approximation. The study of drift velocities in prominences by Gilbert et al. (2002) predicts that there is a slow mass drainage of prominence material (with heavier elements draining faster). This was confirmed numerically by Terradas et al. (2015), who highlighted that the magnitude of the drift velocity is inversely proportional to the magnitude of the collision frequency. The observational study of Gilbert et al. (2007) seemed to confirm that the rate of draining of material from a filament was species dependent. There have also been a number of studies looking at the role of frictional heating by ion neutral drift for waves in the solar chromosphere (Kuźma et al., 2019;Wójcik et al., 2020). It is important to note that for the two-fluid approximation to hold, it is necessary that the particles in a fluid have significantly more interactions with particles in their own fluid than that of the other. If this is not the case, then assuming that the mean-free-path of a particle is controlled by interactions with particles of the same species is flawed and can lead to spurious results, e.g. possibly this is the cause of the overly large viscosity found for neutral species in the simulations of Braileanu et al. (2021a).

Connecting between the two-fluid and single fluid approximations
A big question that faces anyone trying to model partially ionised plasma dynamics in the solar atmosphere is what level of accuracy does one need to model the system, with the main choice being between using a single fluid or two-fluid approximation. The ultimate decision on which is more appropriate will boil down to a question of timescales. We can imagine that if the timescales associated with the collisional coupling are much shorter than the timescale of the dynamics of the system, then it would be natural to expect that the two fluids are strongly coupled such that they only have a small velocity drift and they move almost as a single fluid. Conversely, if the timescale of the dynamics is much shorter than the coupling timescale then the fluids will be decoupled on those timescales and move almost completely independently. There will also be a regime between these where the dynamic and coupling timescales are similar resulting in intermediate or partial coupling of the fluids. But can this be quantified in any way?.
As can be seen from the two-fluid partially ionised plasma equations, the key physical quantities that control the interaction between the two species are the velocity difference (or drift) and the temperature difference. Therefore it makes sense to develop a simple set of equations to perform a reduced analysis of the temporal evolution of these quantities. This can then be used to show how we expect these difference terms to vary with time, and on what timescales this occurs.
We start by looking at the two velocity equations (one for the neutral fluid and one for the plasma fluid): Here some important terms are missing from the equation, e.g. gravity, but to add these is relatively trivial and would not change the arguments we are about to present. By subtracting the second of these equations from the first, we can derive an equation for the time-evolution of the velocity drift, This equation shows there are two ways in which a velocity difference can occur at a given point: 1) The force terms (first three terms on the RHS) can create a velocity difference through driving different flows in the different fluids and 2) the transport of flow by the advection terms (on the LHS) can transport different velocity components from different parts of the flow creating a new velocity difference locally. We then have the coupling term, which through simple inspection we can see will result in the exponential damping of any velocity difference. When looking at this equation, we can see that there are three fundamental timescales (two of which will be discussed in more detail in Section 2.2) of this system, the timescale through which velocity differences are created, the timescale through which they are damped and the timescale through which the collisional coupling terms evolve (i.e. the timescale over which the collision frequencies change). As we are most interested in the consequence of the damping of the velocity drift, for the moment we will focus on the evolution of the velocity drift by investigating its evolution over a period of time that is sufficiently short to ''freeze" the dynamic terms and the collision frequencies as approximately constant. Over this short period, the equation we are analysing becomes with Note that C 1 is not treated as a function of time as these terms are assumed to be not changing over the time period we are considering. From this we can analyse the local evolution, over a short time period, of the velocity drift though a linear ODE (which is mathematically considerably simpler than the Eq. 19). The philosophy behind this set of assumptions comes from numerical simulations, where A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 time integrals are performed over a Dt where the quantities associated with the fluxes are treated as being fixed over that short time step (determined by the CFL condition). This then allows the equations to be integrated as difference equations. Here we treat it as though only the drift terms can evolve with time, fixing the other terms. As we are taking both a c and C 1 as being constant, this allows us to perform a separation of variables in Eq. 20 which leads us to the solution where v D ð0Þ is the velocity drift at t ¼ 0. This solution is of the exponential decay of the drift velocity to a constant value. Taking the limit of the applicability of the approximation to be at t ¼ t 0 , we can see that if a c ðq n þ q p Þt 0 ) 1 then the velocity difference at t 0 (v 0 D ) is given by i.e. a constant. This is known as the strong coupling limit of the velocity drift, and is applied (often in simplified forms, e.g. Braginskii, 1965;Khomenko, 2020) as a correction to the induction equation giving where v CM is the centre-of-mass velocity of the two fluids (e.g. Khomenko, 2020). This allows the mass, momentum and energy equations to be summed and solved as a single fluid approximation (e.g. Khomenko et al., 2014b;Nó brega-Siverio et al., 2020a). A further consequence of this results is that the existence of drift velocities for large (but not infinite) collisional frequencies means that there will be non-zero frictional heating. However, the magnitude of the heating will tend to zero as the coupling gets increasingly stronger (e.g. Hillier, 2019). Conversely, for the case wherea c ðq n þ q p Þt 0 ( 1 then which will be dominated by the initial velocity drift. This implies that over these timescales the fluids are effectively decoupled. Now that we have the easy part out the way, we can approach the trickier step of analysing the thermal coupling. Here we start with the two temperature equations for the two fluids with T D ¼ T n À T p . This allows us to formulate an equation for the temperature difference with C 3 ¼ Àv n Á rðT n Þ À ðc À 1ÞT n r Á v n þ v p Á rðT p Þ þðcÀ 1ÞT p r Á v p . Again we will be analysing the system over a sufficiently short time period such that the terms in C 3 and the collision frequencies can be treated as approximately constant. One thing that can be pointed out here is that C 3 could incorporate any heating term (e.g. Ohmic heating, viscous heating or even radiative losses) which can be considered approximately constant over the timescale of interest. The solution to this equation will have the form T D ¼ pðtÞ þ C expðÀa c ðq p l n þ q n l p ÞtÞ where pðtÞ is the particular integral and C is a constant we can determine once pðtÞ has been found. After some algebra this gives T D ðtÞ ¼ T D ð0Þ À D expðÀa c ðq p l n þ q n l p ÞtÞ ð 29Þ þ C 3 þ C 2 1 ðq p l n Àq n l p Þ acðq n þq p Þ 2 cÀ1 2c ða c ðq p l n þ q n l p ÞÞ À1 þA expðÀa c ðq n þ q p ÞtÞ þ B expðÀ2a c ðq n þ q p ÞtÞ A ¼ q p l n Àq n l p ðq p l n þq n l p Àq n Àq p Þ cÀ1 2c B ¼ q p l n À q n l p ðq p l n þ q n l p À 2q n À 2q p Þ There is one very interesting and often overlooked consequence here, in the limit of strong coupling the temperature difference does not reduce to zero (a standard assumption applied when deriving a single fluid model of partially ionised plasma dynamics, Braginskii, 1965). In fact, it becomes: A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 Here we can see that we can have a temperature difference that, like the velocity drift, is inversely proportional to the rate of the collisional coupling between the fluids. The temperature difference in any region is driven by advection of fluids of different temperatures into a region, different levels of compression of the fluids and the frictional heating created by the drift between the two species. Though it is likely that in many instances these effects are not so large, for example highly compressible dynamics like shocks are able to drive large temperature differences in two-fluid models but this is completely ignored in single fluid approximations of this phenomenon. This formulation may allow a more realistic temperature variance in single fluid calculations of the solar atmosphere where the correctly modelling the electron temperature is important to calculate collisional ionisation and recombination rates.
The alternate situation of very weak coupling is also meaningful to understand. Similarly to the velocity difference, the temperature difference is then determined by its initial value, i.e. T D ðtÞ % T D ð0Þ with a small correction that is linear in time. In this case the temperature difference does not vary greatly over the dynamic timescale, which implies that for fast changes in temperature the fluids are thermally decoupled.
An interesting point of note is this analytic formulation for v D and T D can be used directly in the numerical simulations of two-fluid phenomena by replacing the difference terms in the stiff terms on the RHS of the two-fluid equations presented in the two-fluid equations. This allows these terms to be directly integrated with time over a given timestep analytically due to their exponential form. Therefore, instead of considering using implicit methods, this form of exponential integrator may also prove to be a versatile tool for numerically modelling partially ionised plasma phenomena in a complex, partially ionised atmosphere. A simplified version of this form has been implemented in the (PIP) code (Hillier et al., 2016), but more work needs to be done to see if the it can be used as part of efficient, accurate two-fluid simulations of the Chromosphere.

The multiscale nature of partially ionised plasma dynamics
The previous subsection has shown that the ratio of the dynamic timescale to the collision timescale changes fundamentally how a system responds to any motion, with a key point of comparison being ratio of the characteristic timescales of the collisions to the dynamic timescale. This then leads us to the following question: What is the dynamic timescale of MHD phenomenon in the solar atmosphere?.
An analogy can be drawn between this ratio and the Deborah number, often used in rheology, which is given by the ratio of the relaxation time of a non-Newtonian fluid (which connects with our collisional time) to the timescale of the experiment (which connects to our dynamic timescale). The name for the Deborah number comes from the old testament, with the line from the song by the prophetess Deborah: ''The mountains flowed before the Lord". This name was used in the relaxation of fluids as given a sufficiently long relaxation time (e.g. geological timescales) even mountains can be observed to flow like a fluid. When performing experiments on non-Newtonian fluids it is important to compare data for experiments with similar Deborah numbers as the dynamics fundamentally changes as this number is varied. With Newtonian viscous flow-like behaviour at low Deborah number and non-Newtonian behaviour at high Deborah number.
As there are fundamentally many different timescales associated with dynamics in the solar atmosphere, considering the ratio of the collision frequency to the timescale dynamics of interest is of great importance to decide whether single or two-fluid approximations are most appropriate. However, it is of high importance to note that multiple timescales can all naturally exist in the same system at the same time, meaning the ratio of timescales can be both small and large, and it may also be that to correctly model the scale of greatest interest a wide variety of timescales must be considered. One could say that in solar physics it is important to simultaneously understand both solid and flowing mountains (or how high-frequency, two-fluid dynamics feeds back on larger-scale, slower dynamics). In the rest of this section we will provide some broadstrokes examples of the fundamentally multi-time-scale nature of solar dynamics.

Alfvén waves
Due to both their ubiquitous nature in the solar atmosphere, and the relative simplicity in providing analytical formula, it is good to start with waves. Focusing on the Alfvén wave (following Soler et al., 2013b), we can determine the frequency for an Alfvén wave in a two-fluid medium from the solutions of the following cubic equation where V A is the Alfvén speed in the plasma, e.g.
As this is a cubic equation its solutions give the frequencies for three different waves. These are the forward and backward propagating Alfvén waves, and a neutral shear wave. Due to the non-propagating nature of the neutral shear wave (and arguments on symmetry) we would not expect this extra wave to have any real frequency and analysis shows that this is purely damped (Soler et al., 2013b). For the two other solutions, generally these are the forward and backward propagating, damped Alfvén waves though there are situations when these both become purely damped. It is rather simple to look at the limits for Eq. 35. If a c ðq n þ q p Þ is small compared to x then we just have A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 i.e. one zero frequency neutral shear wave and two Alfvén waves, with the Alfvén wave frequency determined by the plasma Alfvén speed. For the photosphere, with its very high collision frequency this regime would occur when ReðxÞ is at least greater than 10 9 s À1 , though nearer the top of the chromosphere this becomes ReðxÞ'10 3 s À1 due to the drastically smaller collision frequency at this height. For the opposite extreme, a c q p is large, Eq. 35 reduces to here we have two waves, forward and backward propagating Alfvén waves where the wave frequency is determined by the bulk Alfvén speed, i.e. V At ¼ B= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q n þ q p p . Therefore, two-fluid Alfvén waves can behave like a single-fluid MHD Alfvén wave or like an Alfvén wave where the plasma doesn't know there is a neutral fluid (or somehwere in between). Either of these extremes can be easily excited. As a note, there has also been some detailed investigations into compressional MHD waves, and the reader is directed to Soler et al. (2013a) and Alharbi et al. (2022) for more details on these cases. Here we will look in greater detail at shock waves as a fundamentally multi-scale phenomenon (see Section 3).

Instabilities and magnetic reconnection
Moving on from waves it is natural to consider instabilities. When looking at the linear stability, similar arguments can be put forward as those for waves: if the growth rate is large enough then the two fluids will decouple during the growth of the instability, but for small growthrates they will behave like a single fluid. However, it can be difficult, sometimes impossible, to know a priori which of these regimes will be the most important for understanding the dynamics of interest. Nor is it clear a priori whether two-fluid effects will manifest themselves in an important contribution to the nonlinear dynamics of the instability. We will look at instabilities in greater detail in Section 4.
Another important dynamic phenomenon of the solar atmosphere is magnetic reconnection, with observations showing that fast, bursty reconnection can occur in the solar chromosphere . The fundamental timescale of reconnection is the Alfvén time based on the halflength of the current sheet as it gives a timescale for material to be ejected from the reconnection region. However, the dynamics in a reconnecting current sheet are more complex with plasmoids forming through resistive instabilities (Zweibel, 1989) and then interacting to drive multi-scale, potentially fractal reconnection (Singh et al., 2015;Singh et al., 2019). This naturally results in a multi-scale process that can scale the dynamic frequencies covered by singlefluid, two-fluid and into kinetic dynamics (Singh et al., 2015).

Turbulence
Finally, it is very interesting to look at turbulence. This is the archetypal example of multi-time-scale, multi-spatial-scale dynamics in a fluid system, and both MHD and partially ionised plasma turbulence also have this property. As such the existence of a turbulence cascade can take energy from large-scale, well-coupled flows and inject it in scales that are becoming more and more decoupled (e.g. Burkhart et al., 2015). This is easy to see from a simple argument, the frequency of motions at any scale in a steady-state turbulent cascade can be approximated by where is the energy cascade rate of the turbulent cascade and L is the dynamic scale under consideration. This frequency increases as the scale under consideration of the turbulent cascade decreases. Assuming that viscosity does not truncate the cascade too soon, there will exist scales L m ni and L m in such that the turbulent frequency surpasses the neutral ion and ion neutral collision frequencies, respectively. Therefore, during a turbulent cascade in partially ionised plasma it would be natural for highlycoupled dynamics to become decoupled requiring twofluid modelling.
In the following sections we will look in more detail at a few examples where the natural existence of multiple scales of interest or multi-scale dynamics make them key areas in which systems can naturally exhibit two-fluid partially ionised plasma effects. We will focus on shocks (Section 3) and and linear and non-linear instabilities (Section 4) in two-fluid systems.

The role of partially ionised plasma in shocks
Shocks are highly-nonlinear, highly-compressible waves characterised by very sharp transitions in density, velocity, pressure and, for MHD shocks, magnetic field around the shock front. This sharp transition in physical quantities results from the magnitude of the flow in a system abruptly changing from above a characteristic wave speed to below, with the three wave-speeds of MHD (slow, fast and Alfvén) leading to three different shock jumps (e.g. Delmont and Keppens, 2011). However, as we have seen for partially ionised plasmas when considering Alfvén waves (see Section 2.2 or Soler et al., 2013b) or for the slow and fast magnetoacoustic wave modes (e.g. Soler et al., 2013a;Alharbi et al., 2022) there are multiple version of the same wave speed. As such a partially ionised plasma system could possess multiple different flavours of the three types of shocks (potentially simultaneously).
The possibility of this multiple flavours of different shocks has to be understood in the context of their lifetime. Once a shock has propagated a sufficient distance, the upstream and downstream states should couple together making a system that on large scales behaves like strongly coupled MHD. The work of Hillier et al. (2016) studied the evolution of a switch-off shock (a type of slow-mode MHD shock) in a two-fluid system highlighting just such an effect.
A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 Initially a strong shock was created in the plasma, but as this propagated the neutrals began to respond, forcing this plasma shock to evolve towards one that obeyed a shock transition of a highly coupled system. This leads to fundamentally two timescales in the shock problem, the highfrequency timescales of the shock front, and the evolving timescales relating to the coupling between the fluids in the regions around the shock, this latter timescale is related to the length of time a shock has propagated. This makes shocks a good example where two-fluid dynamics may connect to the larger-scale, slower-evolving dynamic evolution of the solar atmosphere. Due to the importance shocks are believed to play in the heating of the lower solar atmosphere (e.g. Hollweg et al., 1982;Brady and Arber, 2016) and their clear observations (e.g. Chae et al., 2018;Houston et al., 2018;Houston et al., 2020) it is important to understand the nature of shocks in the partially ionised layers of the solar atmosphere. The study of partially ionised plasma shocks in relation to the solar atmosphere is a relatively new field, with a lot of the understanding coming as an extension of the work looking at shocks in the interstellar medium. For a review on this topic see, for example, Draine and McKee (1993). In this section we review some of the studies of particular importance for understanding two-fluid effects in shocks in the lower solar atmosphere.

The PIP shock jump
One of a theorists fundamental tools to analyse and understand a given shock are the shock jump conditions (e.g., Goedbloed et al., 2010). These are a set of 1D conditions on the mass, momentum, energy and magnetic fluxes to enforce conservation of these quantities across a shock front. The ideal MHD equations can be moved to the Hoffman-Teller frame (the shock frame with zero electric field both upstream and downstream of the shock) then integrated to yield: where this notation means that for any quantity Q. These can be reduced to a single equation that gives the possible stable shock jumps for a specified upstream plasma beta and upstream angle of magnetic field in terms of the upstream and downstream Alfvén Mach numbers (Hau and Sonnerup, 1989): for the Alfvén Mach number Bx evaluated in the upstream and downstream media. This is a very powerful tool as it allows us to calculate the types of shocks that can exist in a system from minimal details of the upstream medium. Note that this form is equivalent to the shock adiabatic including a linear (trivial) solution of A u This analysis can be extended to the two-fluid system given by Eqs. 1,4,5,10 (Snow and Hillier, 2019). For simplicity, the ion and neutral species are assumed to be coupled by thermal collisions only (with gravity, viscosity, magnetic diffusion etc. neglected). In the Hoffman-Teller frame, the two-fluid PIP equations become: Note that in these equations, integral terms (I 1 ; I 2 ; I 3 ) exist that govern the exchange of momentum and energy between the ion and neutral species. However, these integral terms can not be evaluated easily. Instead, due to the conservative nature of the system, the integral terms can be removed by adding the neutral and ion equations together, reducing the equations to: A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 q n v ?n þ q p v ?p ¼ const; The equations can be further simplified by expressing the partial pressures and densities as a total value using the neutral fraction n n : Substituting these into Eqs. 60,65 and furthermore imposing an additional constraint such that upstream and downstream of the shock the drift velocity equals zero Eqs. 70,75 are identical to the MHD jump equations (Eqs. 39,44) except here the density and pressure are given by their total quantities q t ; P t . Therefore, the solution to these equations is identical to the Hau and Sonnerup, 1989 solution for MHD (Eq. 46), independent of the neutral fraction. It should be noted that this is only true over the larger shock structure and inside the shock there is substructure that is highly dependent on the collisional effects. This result should not be of any surprise. The collisional coupling terms conserve both momentum and energy, so any model that looks at the overall conservation of these quantities over both fluids will reduce to the MHD solution sufficiently upstream and downstream of the shock.

Finite Width shocks and their substructure
In the previous subsection we have seen that at least for a steady-state shock front the inclusion of partially ionised plasma effects does not change the overall shock jump as this is determined by conservation of total mass, momen-tum and energy as well as magnetic flux. However, the inclusion of dissipative terms means that the discontinuous shock front that exists in ideal, inviscid MHD and inviscid hydrodynamics can no longer be sustained. Instead, the shock now has a finite width that is determined by a physical length scale of the dissipation terms (e.g. Hau and Sonnerup, 1989;Draine and McKee, 1993).
As the ionised and neutral species enter the shock they decouple before recoupling as they pass through the shock, see Fig. 2. This allows large drift velocities to develop in the shock, which results in frictional heating within the finite shock width. The thickness of the shock naturally scales with linearly with the inverse of the magnitude of the collisional coupling (e.g. Draine and McKee, 1993;Snow and Hillier, 2021b). However, for other parameters of the system, the correlation to the shock width is not so straight forward. For example, the scaling of the thickness of the two-fluid, switch-off slow-mode shock to the ionisation fraction (n i ) was found to be proportional to n À1:2 i (Hillier et al., 2016).
What happens inside this finite width, as the fluids decouple and recouple, can be very complex, though the basic steady-state solutions can be generally categorised as either continuous (C-) shocks or jump (J-) shocks (Draine and McKee, 1993). For the C-shocks, all the physical quantities vary smoothly over the shock jump. However for the J-shocks there exists a further shock jump inside the larger shock transition. Often this is a hydrodynamic shock in the neutral fluid (e.g. Draine and McKee, 1993;Hillier et al., 2016), which can form as a result of the supersonic velocity forming in the neutrals through being dragged by the plasma (Draine and McKee, 1993) or pushing ahead (Hillier et al., 2016) depending on the type of shock. Transient, though long lived, shocks can also form as the system evolves towards its steady state, with intermediate shocks found to occur in switch-off shocks (Snow and Hillier, 2019). These were found to be A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 created as the plasma gets accelerated in the shock becoming super-Alfvénic, allowing and intermediate shock transition to occur, accompanied by the tell-tale reversal of the magnetic field. This reversal of the magnetic field cannot be supported by the neutral pressure resulting in the slow decay of the shock, hence why this feature is transient.
There is an interesting analogy that can be drawn here between the finite width and substructure of partially ionised plasma shocks and similar physics that can be found in hydrodynamic shocks with thermal conduction. In a hydrodynamic shock with thermal conduction, the shock front becomes broad. However, for sufficiently strong shocks at low Prandtl number (ratio of viscosity to thermal conduction), the system develops an internal, isothermal shock (Landau and Lifshitz, 1987;Guidoni and Longcope, 2010). For this shock, the thickness of the global shock transition is fixed by thermal conduction, but the thickness of the internal sub-shock is fixed by viscosity. This is a very similar situation with partially ionised plasma shocks with an internal hydrodynamic shock. the global structure is determined by the coupling length of neutrals to the plasma (and the magnetic field), but the thickness of the internal shock is determined by the shorter lengthscale of coupling of the plasma to the neutrals.

Stratification and partially ionsied plasma shocks
Compressible waves efficiently steepen into shocks as they propagate upwards in the solar atmosphere due to the decreasing density with height creating increasing nonlinearity in the wave (e.g., Suematsu et al., 1982). These shocks can form in the lower solar atmosphere, where partial-ionisation plays a key role (for example, umbral flashes Beckers and Tallant, 1969). This leads to the obvious question: what influence does partial ionisation play on shock formation and evolution as a shock propagates upwards through the solar atmosphere?.
When stratification is included, an additional complexity arises. In the two-fluid description described in Section 2, the electrons contribute a pressure (that is equal to the ion pressure) but no density. As such the pressure scale height of the plasma (ions + electrons) is defined as K p ¼ 2T p cg , which differs from the neutral pressure scale height K n ¼ T n cg at equal temperatures. In the two-fluid description described in Section NUMBER, both the ions and the electrons contribute to the plasma pressure scale height, as such it differs from the neutral pressure scale height by a factor of 2, assuming equal temperatures. The different pressure scale heights of ionised and neutral species helps lead to a change in composition of the medium from the mostly neutral photosphere, to the almost entirely ionised transition region. Coupled with this is the reduction in density resulting in weaker coupling between the fluids as the shock propagates upwards. As such, a propagating shock is subject to different levels of coupling with height with different levels of importance for each fluid, which can greatly affect the shock heating of the med-ium (Niedziela et al., 2021;Zhang et al., 2021). The large neutral fraction in the photosphere results in a relatively strongly coupled media, however as waves or shock propagate further upwards, the coupling becomes weaker (Braileanu et al., 2019a;Braileanu et al., 2019b). This can lead to wave damping (and frictional heating) due to the ion-neutral interactions for both partially ionised waves and shocks. Generally it was found that both smaller period (which excite more two-fluid effects) and larger amplitude (which make shocks form earlier) wave drivers had the strongest damping (Braileanu et al., 2019b). It has been suggested that damping of upward propagating partially-ionised waves alone is a potential way to balance radiative losses from the lower solar atmosphere (Wó jcik et al., 2020).
Both waves and shocks can also undergo mode conversion, where one wave mode becomes another, as they propagate upwards in the solar atmosphere (e.g., Schunker and Cally, 2006;Khomenko and Cally, 2019). Since the mode conversion height (where the Alfvén speed equals the sound speed, i.e. V A ¼ C s ) occurs in lower solar atmosphere, it becomes of interest to consider how this process occurs for partially ionised plasmas. In a partially ionised plasma, the height at which mode conversion for waves and shocks occurs can be multiply defined depending on the level of coupling (Snow and Hillier, 2020). For a fully coupled system, the characteristic MHD wave speeds depend on the bulk parameters, i.e., V At ¼ B= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q p þ q n p ; , whereas for a completely decoupled systems the MHD wave speeds depend on the isolated plasma properties: V A ¼ B= ffiffiffiffi ffi q p p ; C s ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi cP p =q p q (as discussed in Section 2.2). This provides upper and lower limits limits for where mode-conversion can occurs with the exact height coming from the effective wave speeds calculated by solving a dispersion relation including the coupling frequencies (Soler et al., 2013a). For shocks in particular, as they pass through the C s =V A ¼ 1 height, they split into slow-and fastmagnetoacoustic modes. These can either be shocks or smoothed waves depending on the inclination angle of the magnetic field (Pennicott and Cally, 2019). For a partially ionised plasma, this process can be rather nontrivial as multiple waves and shocks can form depending on the level of coupling of the system (Snow and Hillier, 2020). Since the slow-magnetoacoustic plasma speed and the sonic plasma speed are relatively close, the slow component couples fairly readily to the neutral species for moderate levels of coupling. However, the relatively large difference between plasma fast-magnetoacoustic and neutral sonic can result in stark differences to the MHD model. Since the plasma fast magnetoacoustic mode speed is much larger than the neutral sound speed, the neutrals are dragged along by plasma, which can result in localised frictional heating, as well as increasing the finite width of fastmode shocks. In the study by Snow and Hillier (2020) the finite width of the fast-mode shock was found to exceed the A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 pressure scale height of the system for a wide range of coupling coefficients. Correspondingly, one may expect that fast-mode shocks in the lower solar atmosphere have a finite width of % 300 km. Note that despite the large shock width, a fast-mode shock transition in still satisfied at the location of the maximum velocity gradient.

Stability of Shock Fronts in partially ionised plasma
In hydrodynamics, shock fronts are considered to be very stable structures, at least in the absence of strong radiation (Laming and Grun, 2002;Grun et al., 1991). This can be understood by the looking at the baroclinic term in the hydrodynamic vorticity equation. As the pressure and density gradients of a shock both point in the same direction, perturbations to hydrodynamic shock fronts typically decay with time. Another way of stating this is the steady state shock solution does not allow any vorticity at the shock front (Zhou et al., 2021).
The inclusion of a magnetic field perpendicular to the shock front fundamentally changes this. The steady state shock front can have vorticity (in the form of jumps in the flow parallel to the shock front) but contact discontinuities cannot (Zhou et al., 2021). A consequence of this is MHD shock fronts may be unstable, but instability of contact discontinuities is suppressed. Fast mode shocks are found to be categorically stable to this instability (provided c < 3 Gardner and Kruskal, 1964), whereas the slow-mode shocks may develop the corrugation instability (which corrugates the shock front). The conditions for stability of these shock fronts depends on the Alfvén Mach number and the angle of the magnetic field relative to the shock front (Stone and Edelman, 1995;É del'Man, 1989;Lessen and Deshpande, 1967).
For a parallel shock, where the velocity and magnetic field are aligned, the corrugation instability grows unconditionally in MHD (Stone and Edelman, 1995), however the corresponding HD shock is unconditionally stable. For a partially ionised system, it is less easy to determine if the shock front will be stable or unstable to the corrugation instability. An argument of timescales tells us that both an instability that grows on timescales much shorter than the coupling time and much longer than the coupling time will experience different MHD limits, allowing them both to be unstable. However, a stability range exists whereby the coupling allows the neutral component to suppress the instability (Snow and Hillier, 2021b). Fig. 3 shows the evolution of the corrugation instability in a partially ionised plasma with different levels of coupling, and MHD simulations of the zero and infinite coupling cases. The MHD simulations are unstable, whereas the twofluid simulations can become stable for finite coupling. A caveat here is that this has only been studied in 2D parallel shocks. Shocks that are stable in 2D can become unstable  Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 to the corrugation instability in 3D due to wide varieties of angles that form between the velocity and magnetic field (Stone and Edelman, 1995;É del'Man, 1989). Propagating shocks in the solar atmosphere are almost guaranteed to encounter perturbations and hence corrugate due to the inherent inhomogeneity of the atmosphere. In the lower solar atmosphere, which is known to be both partially ionised and abundant with shocks, one can use characteristic properties of the media to estimate the stability range with respect to perturbation wavelength. Applying this to a sunspot, where the field is reasonably aligned with the direction of propagation and hence the analysis of partially ionsied parallel shocks in Snow and Hillier (2021b) becomes applicable, one can estimate that shocks are stable to perturbations between 0:6 and 56 km.

Shock ionisation and cooling
The rapid changes in density and temperature across a shock lead to locally enhanced ionisation and recombination rates (e.g., Carlsson and Stein, 2002). As such, these effects become critical in modelling partially ionised shocks and omitting ionisation and recombination can lead to a simultaneous over-prediction of shock heating and underprediction of energy dissipation (Zhang et al., 2021).
In its simplest form, ionisation and recombination appear in the two-fluid equations as an exchange of mass, momentum and energy between the neutral and ionised species such that these quantities are all conserved. Therefore, shock jump equations can be written including these terms and the same analysis can be performed as in Section 3.1 to show that the pre-and post-shock states can be modelled as MHD jumps (Snow and Hillier, 2021a). However, this is only part of the picture. A more realistic model for ionisation and recombination of hydrogen involves radiative losses. Ionisation as a three-body process involves a free electron expending energy to release a bound electron, hence energy is lost from the macroscopic plasma species. The energy lost is proportional to the ionisation rate multiplied by the energy required to release the bound electron (for example, 13.6 eV for ground state hydrogen). Including the ionisation loss modifies the plasma energy equation in the two-fluid model as: where A heat is an arbitrary heating term included to balance the system. Note that since ionisation potential loss terms are a energy sink in the plasma only, the neutral energy equation remains unchanged. A consequence of an energy loss term is that the energy equation is no longer conservative. As such, it cannot be used to construct shock jump relations. Instead, for specific types of cooling, requirements on the heating and cooling terms balancing sufficiently upstream and downstream of the shock can be used to create a semi-analytical description of the stable shock solutions (for example the empirical ionisation and recombination formulae used in Snow and Hillier, 2021a), as shown in Fig. 4. It is obvious that including the ionisation potential loss into the equations has resulted in a very different set of solutions to the MHD (or conservative two-fluid) solutions. This is in contrast to the conservative equations which reduce to MHD sufficiently far from the shock. The general shape of the solutions are both cubic with intersections with the linear (trivial) solution occurring near the slow, Alfvén and fast wave speeds.
As a result of their being cooling the shocks are much more compressible: in MHD the compressiblity limit is given by r ¼ q d =q u ¼ ðc þ 1Þ=ðc À 1Þ whereas here, the compression across the shock can be far greater. Taking for example the solution at A u x ¼ 1, (i.e., a switch-off slow mode shock), the compression across the shock can be calculated as r ¼ A u2 =A d2 % 1=0:05 ¼ 20, which is far greater than the MHD limit. This compression has been confirmed in the simulations of Snow and Hillier (2021a).
As a consequence of this cooling, it has also been found that the compressible shocks in this model will always be cooler downstream of the shock than upstream. This is because the requirements for equilibrium is that the losses are balanced by the background heating term. Equating the heating and the losses, one sees that a compression of the plasma (and the increased density this implies) can only be balanced by a post-shock reduction of temperature.
Until now, the physics we have been discussing would apply to a discontinuous shock as much as a partially ionised plasma shock with a finite width. However, one large difference will exist for these different scenarios: in a shock with a finite width there is the potential for cooling to happen inside the shock front stopping the plasma A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 reaching the maximum predicted temperature of the ideal shock jump. A manifestation of this effect is seen in molecular clouds where molecules that should be dissociated due to the temperature rise predicted by the theoretical MHD shock were found to still exist in the post shock media. The radiative losses inside the shock were found to be crucial in reducing the maximum temperature obtained within the shock and hence allow these molecules to survive the shock (Draine and McKee, 1993). For a hydrogen plasma more relevant for the solar atmosphere, the simulations of Snow and Hillier (2021a) showed that the radiative losses will reduce the maximum temperature obtained within the finite width of the shock. Fig. 5 shows the plasma and neutral temperatures across the finite-width of a switch-off shock for different recombination timescales. Within the shockfront, the ionisation and recombination rates are increased by several orders of magnitude such that even when the background recombination timescale is much longer than the evolution time of the system, the plasma within the shock has cooled to slightly below the MHD limit (orange line). For larger background rates, the cooling is more extreme. This could have important consequences for heating, ionisation and line formation around shocks in the chromosphere.

Instabilities in solar partially ionised plasma
Instabilities are another physical process where we can find the importance of partially ionised plasma effects. This might be through changes in the linear growth rate (e.g. Díaz et al., 2012;Soler et al., 2012) or important nonlinear dynamics driving large velocity drifts (e.g. Braileanu et al., 2021b). In this section we review some of the key results that have been found to date for instabilities in two-fluid modelling relating to the solar atmosphere. For details on a wider range of instabilities and their applications to different astrophysical systems the readers are referred to Soler and Ballester (2022).

Linear instabilities in a partially ionised plasma
A brief outline of the role of partial ionisation in linear stability theory was provided in Section 2.2. As with linear wave theory there is a simple argument of timescales to determine whether an instability develops independently in one of the fluids or as a coherent instability in both. If the timescales of the instability growth in one of the fluids is fast compared to the collisional coupling then a decoupled instability (with the eigenvector of the instability dominated by terms for only one of the fluids) will occur. If the timescales are slow, then the instability-driven motions will be relatively similar in both fluids. In this regime, there will be a leading fluid (that is most unstable) and a following fluid (that is being dragged along). The leading and following fluids can be characterised by the magnitudes of the different terms in the eigenvector. Soler et al. (2012) investigated the growth of the Kelvin-Helmholtz instability (an instability driven by shear flows) in a compressible two-fluid system. The growth rate (c KHI ) found in the incompressible limit of their calculations for a range of velocity differences (the driver of the instability) is shown in the top panel of Fig. 6 where instabilities dominated by both neutral species modes and charged species modes can be seen. As the Lorentz force works to suppress instability in the plasma, naturally these modes are less unstable than the neutral-dominated modes. Three different levels of coupling were investigated, showing that the growth rate reduces for both neutral-dominated and charged species-dominated modes as coupling increases. For the neutral-dominated modes this is a sign of the drag of the charge-species slowing the growth of the instability to pull it closer to that of a single instability mode based on the stability properties of the bulk partially ionised plasma. For the charged species, as separate neutraldominated and charged-species dominated modes do not exist in the fully coupled limit, the charged speciesdominated modes have their growth rate decrease towards stability. Extending these ideas to a prominence thread, Martínez-Gó mez et al. (2015) showed that even though magnetic fields might be able to stabilise the instability for a fully ionised plasma, the partial ionisation of the prominence plasma would allow for instability at observed flow speeds.
In a similar vein, the linear analysis of the Rayleigh-Taylor instability (an instability driven by baroclinicity, anti-alignment of density and pressure gradients, at a density inversion) in a partially ionised plasma has been investigated in both single-fluid (Díaz et al., 2014;Ruderman et al., 2018) and two-fluid frameworks . The growth rate of this instability (c RTI ) in the incompressible limit of the two-fluid instability against the magnitude of the driving buoyancy term is shown in the bottom panel of Fig. 6. As with the work of Soler et al. (2012), the Lorentz force is working to suppress the instability in the charged species-dominated modes, making the neutraldominated modes the most unstable. Similarly, again, to A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983Soler et al. (2012, as the coupling is increased the neutraldominated mode morphs into the coupled-fluids mode and the isolated charged species-dominated mode tends towards stability. The addition of compressibility in both the studies of Soler et al. (2012) and Díaz et al. (2012) leads to a further effect suppressing instability (energy has to be spent compressing the fluid instead of driving motions, e.g. Ruderman, 2017). This can switch the most unstable fluid from the neutral fluid to the plasma for these instabilities, or create situations where in a previously unstable system neither fluid is unstable at all, resulting in a complex picture of the stability (for example see Figure 6 of Soler et al., 2012).

Partial Ionisation in nonlinear instability dynamics
Linear stability theory is a powerful tool to understand whether a particular system can become unstable, however once that instability has grown and developed the subsequent behaviour is the realm of nonlinear theory. Due to the complexity of nonlinear dynamics we often require numerical simulations to make progress, the same is the case for instabilities in two-fluid partially ionised plasma models. For example, the work of Jones and Downes (2011) and Jones and Downes (2012) shows that for the Kelvin-Helmholtz instability in the nonlinear regime the decoupling of the bulk flow from the magnetic field results in a vast reduction in the magnetic energy of the nonlinear system. This instability has been observed to develop as a result of flows of prominence material (Berger et al., 2017;Yang et al., 2018;Hillier and Polito, 2018), as such merits understanding how partial ionisation will change its nonlinear behaviour.
The work of Hillier (2019) looks in detail at the nonlinear evolution of the Kelvin-Helmholtz instability in the case where multiple linear modes, spanning from coupled to decoupled, are excited in a system with both a velocity and density jump, and a stream-wise magnetic field. The natural response in this system is for neutral fluid to be unstable, but the magnetic field to suppress the development of the instability in the plasma leading to the drift of neutral material across the magnetic field. This can be seen clearly in Fig. 7 where the neutral vortices are moving the fluid across the magnetic field. These vortical structures were found to be where strong velocity drifts were present (and with this frictional heating). However, as the physical scale of the nonlinear instability layer became larger, the effective coupling became stronger reducing the velocity drift and with that the total heating from frictional heating in the layer (which for large layer widths scaled as the inverse of the layer width Hillier, 2019).
One of the key physical responses found in this system was in the thermal coupling. As neutral material moved across the magnetic field it interacted with plasma of a different temperature. When the ionisation fraction is small, the neutral fluid acts as a large source/sink of thermal energy, driving temperature changes in the local plasma. This leads to plasma pressure imbalances along the magnetic field, resulting in compressible motions developing as the plasma pressure works to create a pressure balance along the magnetic field. This turbulent transport of heat across the magnetic field by neutral drift has the potential to be very important when considering how prominence material or spicules interact with the coronal material that surrounds them, potentially aiding the mixing/cooling process proposed by Hillier and Arregui (2019). Martínez-Gó mez et al. (2021) studied the evolution of the Kelvin-Helmholtz in a partially ionised plasma, but with initially no magnetic field. They used Biermann battery term (e.g. Kulsrud and Zweibel, 2008) in induction equation which creates magnetic fields though the electron pressure and density distributions becoming out of alignment. In the simulations of Martínez-Gó mez et al. (2021) Kelvin-Helmholtz vortices drove this term to create a magnetic field perpendicular to the plane of the simulation. For A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 simulations where there was no collisional coupling between the fluids, there was (as would be expected) no variation in the results on the ionisation fraction. However, when the coupling terms between the fluids were switched on, it resulted in stronger magnetic fields being produced by a factor of $ 1:5. Even for these relatively simple settings, there are still many fundamental instability driven physical process where the role of partial ionisation is still to be quantified.
A key example of this would be the vortex disruption process investigated by Mak et al. (2017). In their study flow instabilities would wind up magnetic field until it became sufficiently stressed, underwent magnetic reconnection, ultimately destroying the vortex. How this behaviour will alter in partially ionised plasma, where we will see in Section 4.3 that partial ionisation can change reconnection dynamics, is still to be investigated.

Magnetic Rayleigh Taylor instability and partially ionised prominence dynamics
Moving beyond idealised models, we come to simulations of instabilities in partially ionised plasma that are set up to be directly relevant to understanding the dynamics of plasma of the solar atmosphere. In the studies of Braileanu et al. (2021a,b), the authors used two-fluid nonlinear simulations to investigate the development of the Rayleigh-Taylor instability in a prominence thread to understand how two-fluid effects may be important in the development of plumes in solar prominences (for a review see Hillier, 2018). The work of Braileanu et al. (2021a,b) follows on from the study presented in Leake et al. (2014) and in some sense extending the study performed using the single fluid approximation by Khomenko et al. (2014b).
By seeding the Rayleigh-Taylor instability at a smooth boundary between the predominantly neutral prominence thread and the ionised corona below, it is natural that the gravitational potential energy of theneutral fluid is the driver of the instability in this case . However, the actual linear stability of this system is quite complex. The fully-coupled limit (i.e. ideal MHD) was found to be unstable for any wavelength perturbation when there was no magnetic shear, with magnetic shear being necessary to create a cut-off wavelength (Braileanu et al., 2021a). However, the two-fluid simulations presented a cut-off wavelength below which the instability did not grow even without magnetic shear (Braileanu et al., 2021a). An explanation for this is that the initial conditions of both fluids (both out of equilibrium but held together through the collisional coupling) are baroclinically stable (i.e. the pressure and density gradients in each fluid are aligned). As it is the baroclinic term (rq Â rP ) in the vorticity equation that drives the Rayleigh-Taylor instability (e.g Zhou et al., 2021), it is only when the coupling is sufficient to force the fluids to work together that they become a single baroclinically unstable fluid and can produce instability.
In the nonlinear stages of the instability, large velocity drifts of the order of 1 km s À1 between the neutral fluid and plasma developed (Braileanu et al., 2021a). These are drift velocities driven by both the natural separation of the fluids in the development of Rayleigh-Taylor plumes, and through the interaction of plumes especially in the case where magnetic shear means that current sheets develop as the plumes interact. By investigating how the magnitude of the coupling changed the nonlinear dynamics, they could show that lower coupling leads to more diffuse structures in the neutral fluid, because the fluid was able to slip more efficiently across the magnetic field (Braileanu et al., 2021b), though the stability of the individual fluids (as discussed above) may also have played an important role in this process as well.
This work was extended to include ionisation and recombination in Braileanu et al. (2021b). In Fig. 8 the left panels show the neutral and plasma densities when ionisation and recombination are included, the right when those terms are switched off in the equations. Though the density evolution of the neutral fluid is not noticeably effected, the plasma density is vastly different, with the downward plumes significantly more visible when ionisation and recombination allow the ambient plasma to respond to the neutrals as they drift across the magnetic field. A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983

Instabilities and magnetic reconnection
In terms of instabilities relating to magnetic reconnection there are two key instabilities to consider: the tearing instability (e.g. Zweibel, 1989) that results in the development of plasmoids in a current sheet, and the coalescence instability (e.g. Tajima and Sakai, 1986) that drives the newly formed plasmoids to interact. In combination these instabilities are seen as key for driving turbulence in a reconnecting current sheet, which ultimately results in fast magnetic reconnection (e.g. Loureiro et al., 2007;Loureiro et al., 2012). This process could be important for explaining the presence of fast, bursty reconnection observed in chromospheric plasma Singh et al., 2012;Guo et al., 2020).
In a partially ionised plasma, as well as the natural collapse of gas-pressure supported current sheets (Brandenburg and Zweibel, 1994), the growth rate of the tearing mode can also greatly change. The analysis by Zweibel (1989) highlighted three regimes of the instability in a two-fluid model: the strongly coupled regime where both fluids evolve in the instability, the intermediate regime where the neutrals do not respond but exert a drag on the plasma and a decoupled regime where the instability grows in the plasma fluid. As these different regimes are associated with different levels of coupling, they can be connected with different scales of current sheets in the solar chromosphere (Singh et al., 2015). As a consequence, different physics controls the growth of plasmoids at different scales, modifying the potential fractal nature of plasmoids in a current sheet (Singh et al., 2015). See, for example, the review by Zweibel et al. (2011) for more details.
The development of plasmoids in a reconnecting current sheet in conditions relevant for chromospheric plasma has been studied through numerical simulations. As with fully ionised MHD simulations (e.g. Loureiro et al., 2012) once a sufficiently high Lundquist number (ratio of diffusion time to Alfvén time) is reached reconnecting current sheets in two-fluid partially ionised plasma simulations also become unstable to the formation of plasmoids (e.g. Leake et al., 2012;Ni et al., 2015). However, the change in instability physics with scale modifies the critical aspect ratio of a current sheet at which plasmoids can form (Pucci et al., 2020). The development of plasmoids in the current sheet produces highly time-dependent reconnection with average reconnection rates that are independent of the Lundquist number (Leake et al., 2012).
Once plasmoids have formed in the current sheet a second instability, the coalescence instability mentioned earlier, can play a crucial role in driving the interaction of neighbouring plasmoids. This instability is a current driven instability (developing as two currents flowing in the same direction attract each other). As the Lorentz force is the key driver of this instability it naturally develops on timescales of the order of the Alfvén time. Murtas et al. (2021) performed a detailed study of how the finite coupling between a plasma and neutral fluid changes the coalescence of plasmoids. As we can see from our discussion of Alfvén waves relating to Eq. 35 there are two Alfvén speeds, the bulk Alfvén speed and the plasma Alfvén speed, which results in two possible Alfvén times that could control the dynamics. The results of the parameter study by Murtas et al. (2021) show exactly this (see Fig. 9), weakly coupled systems coalescing on timescales determined by the plasma Alfvén time, and strongly coupled systems coalescing on timescales determined by the bulk Alfvén time. In between these two cluster points was an intermediate regime, covering a couple of decades of collision frequency, where the transition between coupled and decoupled occurs. the general consequence was that partially ionised plasma coalescence could occur on faster timescales than its fully ionised MHD counterpart. Related works by Smith and Sakai (2008) and Sakai and Smith (2009) focusing on plasmoid merging in the lower solar atmosphere to understand penumbral microjets Fig. 8. Neutral (top) and plasma (bottom) densities for the case of Rayleigh-Taylor plumes at the boundary between a prominence thread and the ambient corona form Braileanu et al. (2021b). Differences between the case where ionisation and recombination are included (left) and neglected (right) are shown. Credit: Popescu Braileanu, B., Lukin, V. S., Khomenko, E. et al., A&A, 650, A181, 2021b, reproduced with permission Ó ESO.
A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983. Their results also show that plasmoid merger is sped up in a two-fluid system. Also they highlight that initial large neutral inflow velocities can trigger the strongest reconnection. As plasmoids merge, the current sheet that develops between them can become unstable and form further plasmoids (as can be seen in Fig. 9). In general, it is expected that this subsequent plasmoid development is related to the timescale of the linear growth of the instability in the current sheet compared with the timescale for a plasmoid to be ejected (e.g. Pucci et al., 2020). This multiple plasmoid formation can result in the development of a fractal-like reconnection process (e.g. Singh et al., 2015), with the plasmoids potentially forming down to kinetic scales under the right conditions. In the simulations of Murtas et al. (2021) evidence was found of these tearing unstable plasmoids, i.e. the current sheet had reached a critical threshold and became unstable to the linear tearing instability, but also evidence of a separate plasmoid formation mechanism (called sub-critical plasmoid formation in their paper). In these simulations plasmoids were able to develop even though the current sheet did not reach any particular threshold for linear instability. In this case nonlinear two-fluid dynamics became important for driving secondary plasmoid formation, with the neutral drag changing the flow into the current sheet resulting in it pinching and forming a plasmoid as a result. Clearly the sum of all these effects implies that two-fluid physics can be very important for Chromospheric reconnection.

Discussion and looking towards the observations
In this review article we have looked at the general ideas behind using two-fluid modelling to investigate dynamics of the lower solar atmosphere. The key physical parameter that determines whether a fully ionised MHD or singlefluid partially ionised plasma or a two-fluid model (or on into kinetic modelling) is the most appropriate is the ratio of the dynamic frequency to the collision frequency. However, as has been argued here, many key dynamic phenomena of interest naturally scale multiple frequencies from those where single fluid models are applicable into the regimes where the physics becomes fundamentally twofluid. Shocks, instabilities, and the nonlinear turbulence they create, are important examples where this occurs.
This paper has focused heavily on theoretical studies, and associated arguments and ideas, to understand the current state of our studies of the partially ionised solar atmosphere. The main reason behind this is it has been possible to make a lot of progress in developing theory and models, but conclusive observational studies have proved to be very difficult. Fig. 9. Contour plots of field lines and current for the evolution of the coalescence instability in MHD (left column) and a two-fluid simulation (right column). Reproduced from Murtas, G., Hillier, A., and Snow, B. "Coalescence instability in chromospheric partially ionized plasmas." Physics ofPlasmas, 28(3), 032901, (2021), with the permission of AIP Publishing.
A. Hillier, B. Snow Advances in Space Research 71 (2023) 1962-1983 A number of attempts have been made to measure velocity differences between charged and neutral species in solar prominences. These seemed like the natural structure to make these observations as both the foreground and the background is the hot, tenuous, ionised corona which will not contribute to the emission in the cool spectral lines used to observe prominences. This then allows the Doppler velocity of different spectral lines to be measured, and by comparison of these velocities potentially identify if material of a different charge state is moving at a different velocity.
The results from these studies can be divided into two classes: those that claim an observation of ion-neutral drift (e.g. Khomenko et al., 2016;Wiehr et al., 2019;Wiehr et al., 2021) and those that see drift between all species (including between different neutral species, e.g. Anan et al., 2017). In the study of Khomenko et al. (2016) they measured velocity differences between neutral and ionised species up to a few km s À1 . Interestingly this is a similar magnitude as found in the simulations of Braileanu et al. (2021a).
The explanation given in the studies that find Doppler velocity difference between even different neutral species is that along any line of sight in a prominence the relative distribution of charged and neutral species, and even excitation states of a given neutral species, is not uniform. Therefore, as different packets of material that form our line-of-sight view of the prominence move, the Dopplershift observed in different spectral lines is different as it picks up the motion of different fluid packets and not because it is measuring the local drift between species (Anan et al., 2017). As this explanation can also explain the observations of drifts between neutral and charged species, it presents a challenge to really be able to identify genuine velocity drifts in prominences.
To remove these potential issues, Zapió r et al. (2022) applied strong constraint on the observed intensity of spectral lines based on radiative transfer calculations to remove any opacity effects in the calculation of prominence velocity drifts along with a number of spectral lines. This resulted in only a small area of the prominence being usable for the analysis, but with that area displaying systematic velocity differences of $ 1:7 km s À1 between charged and neutral species. As the systematic velocity differences between purely neutral species were one to two orders of magnitude smaller this could be the first categorical measurement of velocity differences between charges and neutral species in the solar atmosphere. It is interesting that the velocity difference is a systematic shift of $ 1:7 km s À1 as this is different to the velocity difference growing linearly with magnitude of the velocity of the prominence motions as found in Wiehr et al. (2019). Clearly more work, both observational and theoretical, is necessary to explain how these velocity drifts are developing.
Another interesting idea to measure velocity drifts was put forward by Anan et al. (2014), who proposed that the electric field felt by neutral particles as they drift across the magnetic field (i.e. the electric field given by Àv n Â B) could be measurable through the Stark effect. Simulations of shocks suggest this field could be strong inside the shock front (Snow and Hillier, 2019). However, to measure this would require significant polarimetric accuracy. The design of the European Solar Telescope does seem it would likely be suited for performing these observations.
In the field of studying solar partially ionised plasma, it is our belief that the most important future progress will be made through better connection between theory and observations. Further development of observational studies, as discussed above, will be possible now with DKIST online and EST promises greater potential in the future. To make interpretation of the observations possible, effort is also needed to make two-fluid numerical models more realistic (in terms of the excitation and ionisation states of the particles) to allow for a closer one-to-one comparison with observations.

Data availability
There is no data produced for this manuscript.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.