In the single-fluid MHD approximation, the distinction between the different species that compose the plasma, namely ions, electrons, and neutrals, is lost. The governing equations are written in terms of total or averaged quantities that represent the plasma as a whole (see, e.g., Goedbloed et al., 2019). Therefore, the single-fluid approximation is applicable to situations in which the temporal and spatial scales of interest are much larger than the corresponding scales of the interactions between the individual plasma components. In the case of PIP, the relevant temporal scale is the ion-neutral (or neutral-ion) collision time, i.e., the average time required for an ion (neutral) to encounter a neutral (ion). The corresponding length scale is the ion-neutral (neutral-ion) collision length, i.e., the average distance that an ion (neutral) needs to travel to encounter a neutral (ion). Although the details of the ion-neutral interactions are lost in the single-fluid approximation, effective remnants of those interactions remain in the form of nonideal terms in the MHD equations. Thus, the usual expressions of the single-fluid MHD equations for a PIP are D ρ D t = − ρ ∇ ⋅ v , (1) ρ D v D t = − ∇ p + 1 μ ∇ × B × B + ρ g − ∇ ⋅ π ̂ , (2) ∂ B ∂ t = ∇ × v × B − ∇ × η ∇ × B + ∇ × η A ∇ × B × B × B − ∇ × η H ∇ × B × B , (3) D p D t = − γ p ∇ ⋅ v + γ − 1 L , (4) p = ρ R T μ ̃ . (5)

In these equations, D/Dt = ∂/∂t + v ⋅∇ is the total or Lagrangian derivative, ρ is the plasma total density, p is the plasma total gas pressure, v is the center-of-mass velocity, B is the magnetic field, μ is the magnetic permeability, g is the acceleration of gravity, π ̂ is the viscosity tensor, η is Ohm’s diffusion coefficient, η _A is the ambipolar diffusion coefficient, η _H is Hall’s coefficient, γ is the adiabatic coefficient, L is the total heat-loss function of the plasma, R is the ideal gas constant, T is the plasma common temperature, and μ ̃ is the mean molecular weight. Expressions of the coefficients that appear in these equations are not given here for the sake of simplicity and can be found elsewhere (see, e.g., Khomenko et al., 2014a; Ballester et al., 2018)

Eq. 1 is the usual continuity equation for the whole fluid. Eq. 2 is the momentum equation, which includes the pressure force, Lorentz’s force, gravity, and viscosity. The full viscosity tensor, π ̂ , should include contributions from ion, electron, and neutral viscosities. The full treatment of viscosity in a collisional partially ionized plasma is rather cumbersome (see Vranjes, 2014). A common approximation is considering the classic Braginskii (1965) viscosity for ions and electrons, together with an isotropic viscosity for neutrals (see, e.g., Leake et al., 2012; Soler et al., 2015a). In the lower solar atmosphere, neutral viscosity is several orders of magnitude larger than ion viscosity (see Soler et al., 2015b). Eq. 3 is the induction equation, which includes the ideal inductive term followed by several nonideal terms: Ohm’s diffusion term, ambipolar diffusion term, and Hall’s term. Ohm’s diffusion is generally caused by collisions between charge-carriers and neutrals. In the solar atmosphere, Ohm’s diffusion is essentially caused by electron-ion collisions and is greatly enhanced by electron-neutral collisions in the weakly ionized solar chromosphere (see, e.g., Soler et al., 2019). Ambipolar difussion is caused by ion-neutral collisions and typically is the dominant dissipation mechanism in a PIP with strong magnetic fields (Zweibel, 2015). Hall’s term is not a dissipative but a dispersive term that arises because of the imperfect magnetization of ions. In a PIP, ion-neutral collisions enhance this mechanism by further decoupling ions from the magnetic field (see Pandey and Wardle, 2008). The importance of Ohm’s, ambipolar, and Hall’s terms in the lower solar atmosphere was explored by Khomenko and Collados (2012, see their Figure 1), who showed that Ohm’s and Hall’s terms are relevant at lower heights, while the importance of the ambipolar term grows with height and becomes the dominant term in the middle and higher chromosphere. Other terms that could be included in the induction equation are the diamagnetic current term and Biermann’s battery term, which couple pressure gradients with the magnetic field (see, e.g., Khomenko et al., 2014a; Ballester et al., 2018). These terms are omitted here because they are generally less relevant than the included terms (see also a discussion in Cowling, 1956).

On the other hand, Eq. 4 is the internal energy equation, which is written here as an equation for the pressure. All sources and sinks of energy are enclosed within the heat-loss function, namely L = ∇ ⋅ κ ̂ ⋅ ∇ T − Λ + J ⋅ E * − π ̂ : ∇ v + H other , (6) where the various terms, in order from left to right, represent thermal conduction, radiative cooling, Joule heating (with J and E* the current density and the effective electric field, respectively), viscous heating, and other unspecified sources of heating. In Eq. 6, κ ̂ is the thermal conductivity tensor that, as happens for the viscosity tensor, should include contributions from ion, electron, and neutral thermal conductivities. Again, a common approximation is considering the classic Braginskii (1965) thermal conductivity for electrons (ion thermal conductivity is usually negligible), but a simple isotropic thermal conductivity for neutrals (see, e.g., Leake et al., 2012; Soler et al., 2015a). As in the case of viscosity, neutral thermal conductivity is also larger than electron thermal conductivity in the lower solar atmosphere (see Soler et al., 2015b). The cooling function, Λ, accounts for radiative losses. Assuming optically-thin plasma, several parametrizations of the cooling function are available for solar atmospheric plasmas (see, e.g., Hildner, 1974; Athay, 1986; Klimchuk and Cargill, 2001; Schure et al., 2009; Soler et al., 2012a). The Joule heating can be split in the following way J ⋅ E * = μ η J ∥ 2 + μ η C J ⊥ 2 , (7) where η C = η + B 2 η A is the so-called Cowling’s total coefficient (Cowling, 1956) and J _∥ and J _⊥ denote the components of the current density parallel and perpendicular to the magnetic field direction. Thus, in a PIP the dissipation of perpendicular currents is more efficient than that of parallel currents because of the effect of ambipolar diffusion (see also Cowling, 1956; Khomenko and Collados, 2012). Finally, the ideal gas law is considered as equation of state (Eq. 5) to close the system, where the mean molecular weight, μ ̃ , depends on the ionization degree. In a pure hydrogen plasma, it varies from μ ̃ = 0.5 for a fully ionized plasma to μ ̃ = 1 for a neutral gas.

The single-fluid MHD approximation breaks down when the spatial and temporal scales of interest approach the characteristic ion-neutral or neutral-ion scales. In that scenario, the plasma cannot longer be treated as a single fluid and a multi-fluid treatment is needed (see, e.g., Martínez-Gómez et al., 2016, 2017). A specific version of the multi-fluid theory is the two-fluid approximation, in which neutrals are treated as a separate fluid while ions and electrons are still assumed to remain strongly coupled and to form another fluid (see, e.g., Zaqarashvili et al., 2011; Leake et al., 2012; Maneva et al., 2017; Popescu Braileanu et al., 2019). The neutral fluid and the ion-electron fluid interact with each other by means of ion-neutral and electron-neutral collisions, along with ionization and recombination. The two-fluid model is justified by the fact that the frequency of collisions between ions and electrons is typically much higher than the frequency of collisions between ions and neutrals. This is the realistic situation in many astrophysical plasmas, including the solar atmosphere (see, e.g., Figure 11 of Ballester et al., 2018). The basic equations in the two-fluid approximation are D ρ c D t = − ρ c ∇ ⋅ v c − S , (8) D ρ n D t = − ρ n ∇ ⋅ v n + S , (9) ρ c D v c D t = − ∇ p c + 1 μ ∇ × B × B + ρ c g − ∇ ⋅ π ̂ c − R , (10) ρ n D v n D t = − ∇ p n + ρ n g − ∇ ⋅ π ̂ n + R , (11) ∂ B ∂ t = ∇ × v c × B − ∇ × η ∇ × B − ∇ × η H ∇ × B × B − ∇ × χ v c − v n , (12) D p c D t = − γ p c ∇ ⋅ v c + γ − 1 L c , (13) D p n D t = − γ p n ∇ ⋅ v n + γ − 1 L n , (14) p c = n c k B T c , (15) p n = n n k B T n , (16) where all symbols have the same meaning as before, with the subscripts “c” and “n” explicitly denoting quantities related to the charged (ion-electron) fluid and the neutral fluid, respectively. Quantities that have not been defined before are the number densities of charges and neutrals, namely n _c and n _n, and Boltzmann’s constant, k _B. We note that because of the small electron mass, ρ _c ≈ ρ _i, with ρ _i the ion density.

There are some terms in the two-fluid equations that were absent from the equivalent single-fluid equations. The term S in the two continuity equations (Eqs 8, 9) represents ionization/recombination, namely S = ρ c Γ rec − ρ n Γ ion , (17) with Γ^ion/rec the ionization/recombination rates. Explicit expressions for these rates in terms of the plasma properties can be found in, e.g., Meier and Shumlak (2012); Leake et al. (2012); Maneva et al. (2017); Popescu Braileanu et al. (2019) and references therein. The term R in the two momentum equations (Eqs. 10, 11) represents collisional momentum transfer and can be cast as R = α cn v c − v n + ρ c Γ rec v c − ρ n Γ ion v n , (18) where the first term is the momentum transfer due to collisions (that may include both elastic collisions and charge exchange collisions) and the last two terms are associated to the loss or gain of momentum during ionization/recombination processes. We note that Eq. 18 assumes that the velocity drift, v _c − v _n, is much smaller than the thermal speed (Draine, 1986). In Eq. 18, α _cn is the total friction coefficient, which is computed as the sum of the friction coefficients for ion-neutral, α _in, and electron-neutral, α _en, collisions, namely α cn = α in + α en . (19)

For simplicity, expressions for the friction coefficients are not given here and can be found elsewhere (see, e.g., Braginskii, 1965; Draine, 1986, among others). The collision frequency of two species β and β′ is computed from the corresponding friction coefficient as ν β β ′ = α β β ′ ρ β , (20) where ρ _β is the mass density of species β. Although the friction coefficients are symmetric, i.e., α _ββ′ = α _β′β , in general ν _ββ′ ≠ ν _β′β because of the different densities that the two colliding species may have. These collision frequencies introduce the relevant time scales for the interactions between fluids.

The frictional term (Eq. 18) is one of the major forces behind the district dynamics that a PIP displays compared to a fully ionized plasma. In the two-fluid formalism, the frictional term appears explicitly in the momentum equations of charges and neutrals, while in the single-fluid equations this frictional effect was hidden behind the ambipolar diffusion term. Many effects that are caused by ambipolar diffusion of the magnetic field in the single-fluid approximation become more physically transparent in the two-fluid case, where the coupling between charges and neutrals is mediated by a more plain frictional force. For instance, the damping of MHD waves that propagate in a PIP is caused by this friction (see, e.g., Kulsrud and Pearce, 1969; Balsara, 1996; Soler et al., 2013a,b), while such a damping is attributed to ambipolar diffusion of the magnetic field perturbations in the single-fluid model. An analysis of the differences between the two-fluid and single-fluid models using analytical results of MHD waves was performed by Zaqarashvili et al. (2011). The ion-neutral friction is also fundamental to understand the properties of shocks in astrophysical PIP like molecular clouds. In this regard, the formation of C-shocks in which all hydrodynamic variables are continuous is a good example (see, e.g., Draine, 1980; Draine et al., 1983). In connection to the topic on this review, the reader should not be surprised that the frictional term is also relevant for the development of instabilities.

To summarize, the two-fluid induction equation (Eq. 12) contains three differences with respect to its single-fluid equivalent (Eq. 13). Firstly, in the ideal term the whole fluid centre-of-mass velocity, v, is replaced by the velocity of charges, v _c. Secondly, as a consequence of the previous difference, the ambipolar term is absent from the two-fluid induction equation, but this does not mean that the physical effect of ambipolar diffusion is missing form the equations. Indeed, as explained before, the effect of ambipolar difusion consistently remains in the form of the collisional terms in the two momentum equations (Eqs 10, 11). Thirdly, a new term that accounts for electron-neutral collisions is present. This term depends on the coefficient χ given in Ballester et al. (2018) and is equivalent to the coefficient η _D given in Popescu Braileanu et al. (2019). However, this additional term is usually negligible. We note again that some terms that couple pressure gradients with the magnetic field evolution have been omitted from the induction equation because of their minor importance (see, e.g., Popescu Braileanu et al., 2019, for the full expression of the induction equation).

The two energy equations (Eqs 13, 14) now contain separate heat-loss functions for charges and for neutrals, namely L c = Q c c, n + ∇ ⋅ κ ̂ c ⋅ ∇ T c − Λ c + μ η J 2 − π ̂ c : ∇ v c + H c, other , (21) L n = Q n n, c + ∇ ⋅ κ ̂ n ⋅ ∇ T n − Λ n − π ̂ n : ∇ v n + H n, other , (22) where the meaning of the symbols is the same as in Eq. 6. Importantly, we see that the Joule heating term is only included in the charges heat-loss function (see the fourth term on the right-hand side of Eq. 21), and this term is now isotropic because of Ohm’s diffusion. In addition, both L c and L n contain an additional term that accounts for heating and heat exchange due to the interactions with the other fluid, i.e., collisions and ionization/recombination. These terms can be cast as (see, e.g., Martínez-Gómez et al., 2016, 2017; Popescu Braileanu et al., 2019) Q c c, n = 2 α cn m i + m n 1 2 m n v n − v c 2 + 1 γ − 1 k B T n − T c + 1 2 ρ n Γ ion v n − v c 2 + 1 γ − 1 2 k B m i + m n ρ n T n Γ ion − ρ c T c Γ rec , (23) Q n n, c = 2 α cn m i + m n 1 2 m i v c − v n 2 + 1 γ − 1 k B T c − T n + 1 2 ρ c Γ rec v c − v n 2 + 1 γ − 1 2 k B m i + m n ρ c T c Γ rec − ρ n T n Γ ion , (24) where m _i and m _n are the masses of ions and neutrals, respectively. The terms proportional to α _cn in Eqs 23, 24 are related with collisions, while the terms proportional to Γ^ion/rec are related with ionization/recombination. Among these terms, those involving the temperatures, T _c and T _n, represent heat exchange between the fluids and their role is to equalize both temperatures, whereas those involving the square of the velocity drift, v n − v c 2 , produce a net heating of the two fluids and so a global increase of the temperature. In particular, the terms proportional to both α _cn and v n − v c 2 are the two-fluid version of the ambipolar heating term, μ B 2 η A J ⊥ 2 , that appears in the single-fluid approximation. We again note that Eqs 23, 24 are obtained in the limit that the velocity drift is much smaller than the thermal speed (Draine, 1986).

The numerical implementation of the PIP equations is a challenging task because of the presence of diffusive (parabolic) terms and the dispersive (hyperbolic) Hall’s term, which impose extremely small time steps if the temporal evolution is done with standard explicit schemes. For this reason, different strategies are adopted to allow high-resolution numerical simulations to be practical. For instance, a convenient strategy is that of operator splitting, in which the ideal and nonideal terms are evolved separately. Ideal terms can be evolved with regular explicit schemes as in ideal MHD simulations, in which the time step is essentially determined by the largest value of the Alfvén/sound velocity in a low-β/high-β plasma. Then, different methods are used to deal with the nonideal terms.

On the other hand, diffusive terms can be advanced with the super time stepping technique (Alexiades et al., 1996; O’Sullivan and Downes, 2006, 2007), which is able to accelerate the explicit computations by using a mixture of large (unstable) time steps and short (stable) time steps that are determined by Chebyshev polynomials. Globally, the super time stepping scheme remains stable. Alternatively, the super time stepping method can also be implemented with Legendre polynomials, which might offer better stability properties (see Meyer et al., 2014). Other methods to deal with diffusive terms are, e.g., the heavy ion approximation (e.g., Li et al., 2006; McKee et al., 2010; Li et al., 2012) and the use of implicit or mixed implicit-explicit schemes (e.g., Falle, 2003).

Concerning the dispersive Hall’s term, the treatment of this term is even more problematic. The Hall term imposes explicit time steps that tend to zero in situations where it is dominant over diffusive terms (Falle, 2003). To solve this problem, the hyper-diffusivity method (see, e.g., Tóth et al., 2008) introduces a sort of artificial diffusion (or hyper-diffusion) with the aim of stabilizing the Hall term, which is advanced with standard explicit schemes. On the contrary, O’Sullivan and Downes (2006, 2007) proposed the Hall diffusion scheme, in which no arbitrary hyper-diffusivity is required. The Hall diffusion scheme relies on the fact that, because of Hall’s term, the instantaneous rate of change of any one component of the magnetic field depends only on the spatial gradients of the other two components. With this in mind, one can advance one component of the magnetic field explicitly, followed by an implicit-like discretization of the other components. However, if the different magnetic field components are updated in a particular order, then the difference equations are explicit in the sense that no matrix inversions, or approximations of matrix inversions, are required. The Hall diffusion scheme is does not lead to instability unless the Courant condition is not satisfied.

Readers interested in the benefits and limitations of these strategies and methods are referred to Section 3 of Ballester et al. (2018) and references therein for detailed explanations.

In the spirit of the classic treatment of Chandrasekhar (1961b), the simplest situation in which the onset of the KHi can be studied corresponds to the case of two fully ionized unbounded plasmas with uniform densities ρ ₁ and ρ ₂ and separated by an abrupt interface. Both plasmas are permeated by a straight and constant magnetic field, which is orientated parallel to the interface. The magnetic field strength at the two sides of the interface is denoted by B ₁ and B ₂, respectively. In addition, the two plasmas are assumed to be in relative motion, with v ₁ and v ₂ the constant flow velocities at the two sides of the interface. The flow is assumed to be along the background magnetic field direction and the velocity shear across the interface is Δ v = v 1 − v 2 . We note that this set-up composed of an infinitely thin shear layer is ill-defined from a numerical point of view, but it is very convenient from an analytic perspective. A stability analysis of the interface can be done using the linearized single-fluid ideal MHD equations. Considering incompressible perturbations, a temporal dependence of the form exp − i ω t , where ω is the frequency, and a Fourier-analysis in space, the following dispersion relation of the surface waves propagating on the plasma interface can be derived, namely ρ 1 ω − k ∥ v 1 2 + ρ 2 ω − k ∥ v 2 2 − k ∥ 2 B 1 2 + B 2 2 μ = 0 , (25) where k _∥ is the component of the wavevector parallel to the magnetic field and the flow. The analytic solution for the frequency is ω = ρ 1 v 1 + ρ 2 v 2 ρ 1 + ρ 2 k ∥ ± B 1 2 + B 2 2 μ ρ 1 + ρ 2 − Δ v 2 ρ 1 ρ 2 ρ 1 + ρ 2 2 1 / 2 k ∥ . (26)

The + and − signs in Eq. 26 correspond to parallel propagating (forward) waves and anti-parallel propagating (backward) waves, respectively, with respect to the magnetic field direction. There is a critical velocity shear, Δv _crit., denoted by Δ v crit. = ρ 1 + ρ 2 ρ 1 ρ 2 B 1 2 + B 2 2 μ 1 / 2 , (27) which determines the stability of the wave modes. When Δv < Δv _crit. the frequency is real and the modes are stable regardless of the presence of the velocity shear. This is so because of the stabilizing influence of magnetic tension. However, when Δv > Δv _crit. magnetic tension is not able to stabilize the shear flow any more. The physical reason is that the kinetic energy associated with the shear flow becomes larger than the background magnetic energy. Then, the frequency becomes complex. The imaginary part of the frequency of the mode with the + sign is positive and, according to the prescribed dependence exp − i ω t , the amplitude of this mode grows with time. This corresponds to the onset and exponential growth of the KHi. The growth rate of the instability is given by γ KHi = k ∥ ρ 1 ρ 2 ρ 1 + ρ 2 Δ v 2 − Δ v crit. 2 1 / 2 . (28)

The growth rate is proportional to k _∥, which means that small scales grow faster than large scales. In addition, for Δv ≫Δv _crit. the growth rate is approximately linear in the velocity shear. Eq. 27 can be rewritten as Δ v crit. = ρ 1 + ρ 2 ρ 1 ρ 2 ρ 1 v A, 1 2 + ρ 2 v A, 2 2 1 / 2 , (29) where v A = B / μ ρ is the Alfvén velocity. Eq. 29 evidences that the critical shear flow is a sort of root-mean-square average of the Alfvén velocities at the two sides of the interface. This points out that the shear flow velocities needed for the onset of the KHi are necessarily super-Alfvénic. We note that the condition that the densities at the two sides of the interface are different is not actually needed for the triggering of the KHi; only a shear flow is needed. We also note that if the fluid is unmagnetized so that B ₁ = B ₂ = 0, this results in Δv _crit. = 0. In other words, there is no shear flow threshold for the triggering of the KHi in the absence of the stabilizing effect of the magnetic field. It would be the case of a fluid exclusively composed of neutrals. As seen later, this has a very important implication in the case of PIP, where only the charges are directly affected by the magnetic field.

The simple configuration explored above suffices for our later purpose to illustrate the effect of partial ionization. Additional ingredients in the background configuration as, e.g., compressibility, dissipation, flows not strictly parallel to the magnetic field, time-dependent flows, nonplanar geometry, a smooth transition, etc., cause refinements to the theory described above at the expense of a more involved mathematical analysis. It is not our goal to analyze all these additional effects, whose discussion can be found elsewhere. For instance, the effect of dissipation has been discussed in, e.g., Ryutova (2015); Ballai et al. (2015). The effect of replacing the abrupt interface by a smooth transition of thickness δ was explored in, e.g., Miura and Pritchett (1982); Berlok and Pfrommer (2019), who found that modes with k _∥ δ ≳ 2 become stable because of the presence of the inhomogeneity. The KHi driven by time-dependent flow patterns was investigated in, e.g., Browning and Priest (1984). The role of compressibility is worth to mention in more detail because of its relevant influence on the stability of the shear flow. For a flow parallel to the magnetic field, compressibility has a destabilizing effect when the velocity shear is small and the ratio of densities is large, but a strong stabilizing influence for large velocity shear, so that the KHi can even be completely suppressed for sufficiently large Δv (see, e.g., Fejer, 1964; Sen, 1964; Gerwin, 1968; Miura and Pritchett, 1982; Ferrari and Trussoni, 1983; Rae, 1983; Soler et al., 2012b). Importantly, Δv _crit. is not affected by compressibility, so that Eq. 27 obtained in the incompressible limit remains valid. The value of Δv needed for the compressible stabilization depends upon the plasma β and the orientation of the wavevector with respect to the magnetic field direction, with the most stable situation being that of a low-β plasma with the wavevector parallel to the magnetic field (Miura and Pritchett, 1982; Soler et al., 2012b). In that case, compressible stabilization may happen for a velocity shear smaller than Δv _crit., so that the KHi cannot grow. However, a small perpendicular component of the wavevector is enough to allow the development of the KHi for a restricted range of shear flow velocities between Δv _crit. and the value for the compressible stabilization (see Figure 1). More details about the effect of compressibility can be checked in, e.g., Miura and Pritchett (1982); Soler et al. (2012b).

The stability of a shear flow in the case of a magnetized plasma composed by a mixture of charges and neutrals poses a fundamental problem. Considering again the same simple scenario discussed in the previous section, we anticipate that charges and neutrals would have different stability properties. On the one hand, charges are subjected to the stabilizing effect of the magnetic field, so that the KHi would only be possible for super-Alfvénic shear flow velocities larger than Δv _crit.. On the other hand, neutrals do not feel the magnetic force, so that there should not be a threshold velocity shear for the KHi. In the paradigmatic case that charges and neutrals do not interact, a sub-Alfvénic velocity shear may be unstable for the neutrals but stable for the charges. Therefore, the KHi may be triggered in the plasma for a sub-Alfvénic shear because of the presence of a neutral component. However, in a real PIP charges and neutrals interact through collisions. Then, the question arises: Do neutrals remain unstable for Δv < Δv _crit. even when they are collisionally coupled with charges? What is the effect of collisions between charges and neutrals on the onset of the KHi?

The question of whether the KHi can be triggered in a PIP for sub-Alfvénic shear is not only relevant from the theoretical point of view, but it also has important observational implications. Again we restrict ourselves to the solar context. Observations of turbulent flow and vortex structures in solar prominences (see Berger et al., 2008, 2010; Ryutova et al., 2010; Berger et al., 2017) already suggested the presence of the KHi, which was confirmed by recent direct observations (Li D. et al., 2018; Hillier and Polito, 2018; Yang et al., 2018). The typical flow velocities measured in quiescent prominences are in the range 10–70 km s⁻¹ (see Mackay et al., 2010), which are typically sub-Alfvénic. The value of the Alfvén velocity for typical values of the magnetic field strength, ∼ 10 G, and the mass density, ∼ 5 × 1 0 − 11 kg m⁻³, in quiescent prominences is ∼ 126 km s⁻¹. Assuming that the flows are along the magnetic field direction, the observed flow velocities should be stable according to the classic theory for a fully ionized plasma. It has been proposed that partial ionization effects are responsible for the development of the KHi in solar prominence conditions when slow flows are present (Soler et al., 2012b; Martínez-Gómez et al., 2015).

The onset and initial exponential growth of the KHi in a PIP has theoretically been investigated in a number of different astrophysical contexts by, e.g., Hans (1968); Ershkovich et al. (1986); Prialnik et al. (1986); Chhajlani and Vyas (1990); Birk and Wiechen (2002); Watson et al. (2004); Michikoshi and Inutsuka (2006); Shadmehri and Downes (2007, 2008); Kunz (2008); Soler et al. (2012b); Martínez-Gómez et al. (2015). For this problem, the two-fluid approximation offers a more transparent physical picture because the different stability properties of charges and neutrals can be better described than with the single-fluid model. In particular, the work by Watson et al. (2004) is relevant for our present discussion because the case of field-aligned shear flows at a sharp interface is explored. See, e.g., Kunz (2008) for a detailed study in the single-fluid case. In the analysis of Watson et al. (2004), collisions between charges and neutrals are considered in the two-fluid momentum equations, but ionization and recombination are neglected. Also, diffusion mechanisms other than collisions are omitted. In the incompressible limit, Watson et al. (2004) derived a dispersion relation in the two-fluid approximation, which was reobtained by Soler et al. (2012b) with a slightly different notation. The dispersion relation was also recovered in some limits by Shadmehri and Downes (2007) and Martínez-Gómez et al. (2015), who considered a plasma layer and a cylindrical tube, respectively. In the notation of Soler et al. (2012b), the dispersion relation can conveniently be cast as D c D n + D cn 2 = 0 , (30) with D c = ρ c,  1 ω − k ∥ v 1 ω − k ∥ v 1 + i ν cn,  1 + ρ c,  2 ω − k ∥ v 2 ω − k ∥ v 2 + i ν cn,  2 − k ∥ 2 B 1 2 + B 2 2 μ (31) D n = ρ n,  1 ω − k ∥ v 1 ω − k ∥ v 1 + i ν nc,  1 + ρ n,  2 ω − k ∥ v 2 ω − k ∥ v 2 + i ν nc,  2 (32) D cn = ρ n,  1 ω − k ∥ v 1 ν nc,  1 + ρ n,  2 ω − k ∥ v 2 ν nc,  2 . (33)

Here, D c is the dispersion relation associated with the charged fluid (see Eq. 25), which is modified by collisions between charges and neutrals through the terms with ν _{cn, 1} and ν _{cn, 2}. In turn, D n is the dispersion relation associated with the neutral fluid (note the absence of the magnetic field term), which is modified by collisions between neutrals and charges through the terms with ν _{nc, 1} and ν _{nc, 2}. Finally, D cn is a coupling term. Eq. 30 is a forth-order polynomial, so it has four different solutions: two associated with the charges and other two with the neutrals. The role of the coupling term D cn is to entangle the modes associated with charges and neutrals, so that they become global modes of the whole PIP. Numerical investigation of the solutions of Eq. 30 reveals that collisions are not able to completely stabilize the neutral fluid unstable mode for sub-Alfvénic shear (Watson et al., 2004; Soler et al., 2012b; Martínez-Gómez et al., 2015). Collisions reduce the growth rate compared with the collisionless case but a instability is always present for small, sub-Alfvénic shear. Under the approximation of immobile charges and in the case that the densities at the two sides of the interface are the same, Watson et al. (2004) found an approximate dependence of the neutral unstable mode growth rate with the shear flow velocity as ∼ Δ v 2 instead of a linear dependence with Δv as in the hydrodynamic collisionless case that would result from Eq. 28 with Δv _crit. = 0. Martínez-Gómez et al. (2015) obtained a similar quadratic dependence of the growth rate with the flow velocity by considering a perturbation analysis in the full Eq. 30 and the particular cases of v ₂ = 0, v ₁ ≪ v _A,1, and a strong coupling. The approximate growth rate derived by Martínez-Gómez et al. (2015) is γ KHi ≈ 2 ρ n,  1 ρ n,  2 ρ n,  1 + ρ n,  2 k ∥ 2 v 1 2 ρ n,  1 ν nc,  1 + ρ n,  2 ν nc,  2 . (34)

In addition to the quadratic dependence with the flow velocity already found by Watson et al. (2004), Eq. 34 also evidences that the growth rate is inversely proportional to the neutral-charge collision frequencies. Collisions tend to stabilize the KHi in the neutral fluid, but they never manage to suppress it completely. A general investigation for arbitrary coupling strength, showing the transition from the linear dependence to the quadratic dependence in the velocity shear, can be checked in Soler et al. (2012b); Martínez-Gómez et al. (2015), but still in the incompressible case. Martínez-Gómez et al. (2015) applied Eq. 34 to the case of thin threads in solar prominences and computed KHi growth rates compatible with the observations (see Figure 2).

The derivation of a dispersion relation in the partially ionized case becomes much more cumbersome when compressibility is taken into account. A simple, approximate expression for the growth rate is no longer easy to obtain. Soler et al. (2012b) obtained a dispersion relation in the compressible case, whose solutions were numerically investigated. Soler et al. (2012b) found that the domain of instability becomes strongly dependent on the plasma parameters, especially the values of the collision frequencies and the density jump across the interface (see their Figures 6, 7). Soler et al. (2012b) found that, in general, compressibility tends to stabilize the neutral component of the plasma. If charges and neutrals are strongly coupled and the density jump is large, the results from Soler et al. (2012b) suggest that the threshold velocity slightly increases when the density jump increases (see their Figure 7B). For realistic physical properties in solar prominences, the threshold velocity shear for the onset of the KHi remains sub-Alfvénic under certain conditions. In particular, modes with a wavevector forming a sufficiently oblique angle with the magnetic field remain unstable and behave, approximately, as in the incompressible limit (see Figure 3).

Compressible effects in the presence of partial ionization were previously investigated by Ershkovich et al. (1986) and Prialnik et al. (1986), but in both studies the flow of neutrals was assumed to be perpendicular to the interface and not parallel to it as in our exemplary configuration. In their case, compressibility was also found to have either a stabilizing or destabilizing effect on the KHi growth rate depending on the parameters considered.

The case investigated by Watson et al. (2004) and Soler et al. (2012b) neglects all nonideal terms in the induction equation. Pandey (2018) and Pandey and Vladimirov (2019) have explored the triggering of the KHi in the presence of Hall’s term in the induction equation using the single-fluid equations. In the case of fully ionized plasmas, the effect of Hall’s term has been explored by, e.g., Zhelyazkov et al. (2020). Pandey (2018) showed that Hall’s effect opens up a new way through which the fluid can become Kelvin–Helmholtz unstable sub-Alfvénic flows. Pandey (2018) have explored the applicability of this mechanism in different astrophysical contexts, including the solar atmosphere, protoplanetary discs, molecular clouds, and Earth’s ionosphere. As Hall’s effect is enhanced by the presence of neutrals (Pandey and Wardle, 2008), this alternative mechanism may be of relevance in PIP. Readers are referred to Pandey (2018) for details. Recently, Martínez-Gómez et al. (2021) have numerically studied the KHi evolution due to shear flow in an initially unmagnetized plasma including the Biermann battery term in the induction equation. The presence of this term does not seem to strongly affect the linear growth of the KHi because similar growth rates than those predicted by the ideal theory are found in the simulations of Martínez-Gómez et al. (2021). However, the Biermann battery term heavily impacts on the nonlinear evolution by generating a magnetic field (see more details in Section 3.3).

While linear theory is very useful to understand the onset of the KHi, full nonlinear numerical simulations are requiered to study its later evolution. In fully ionized astrophysical plasmas, the nonlinear phase of the KHi has been investigated in detail in both 2D and 3D simulations (see, e.g., Frank et al., 1996; Malagoli et al., 1996; Keppens et al., 1999; Keppens and Tóth, 1999; Ryu et al., 2000; Baty et al., 2003; Baty and Keppens, 2006; Matsumoto and Hoshino, 2004; Matsumoto and Seki, 2010, just to name a few representative papers among many other works). Important results obtained from the simulations are that the nonlinear development of the KHi leads to the formation of vortices, mixing of plasmas, and eventual generation of turbulence. Particularly important is the transition from the laminar to the turbulent regime associated with the later evolution of the plasma mixing (Matsumoto and Hoshino, 2004; Matsumoto and Seki, 2010) and its implications concerning the energy cascade towards the dissipative scales. The magnetic field plays an important role in this process. For moderate shear flow velocities, the magnetic field can nonlinearly inhibit the formation of vortices, thus making the system nonlinearly stable although being linearly unstable (Ryu et al., 2000; Hillier, 2019). During the vortex formation, the winding of the magnetic field can cause field amplification, formation of current sheets, and field line reconnection. Another common feature of the KHi nonlinear evolution obtained from simulations is the formation of magnetic islands (see, e.g., Keppens et al., 1999; Nakamura et al., 2008). Readers are referred to the papers cited above, and references therein, for further details.

Considering nonlinear simulations of the KHi in multi-fluid plasmas, the effects of the interactions between the various species have been explored in the last decade in a relatively small number of publications. Most of the studies (see, e.g., Birk and Wiechen, 2002; Johansen et al., 2006; Wiechen, 2006; Barranco, 2009; Hendrix and Keppens, 2014) deal with the case of partially ionized dusty plasmas and focus on analyzing the effect of dust. Typically, the dust grains are treated as particles with a certain mass and charge (they may be neutral) that collide with ions, electrons, and neutrals that may also be present in the plasma. As in the single-fluid fully ionized simulations, these works show that once unstable modes are excited they evolve relatively quickly nonlinearly and result in the formation of vortices, current sheets, and turbulence. Wiechen (2006) concludes that the mass and charge of the dust particles affect the stability, so that a higher mass of the dust grains has a stabilizing effect, while a higher charge number of the dust has a destabilizing effect. With hydrodynamic simulations, Hendrix and Keppens (2014) also found a stabilizing effect of dust when the dust-to-gas ratio is high. Also, Hendrix and Keppens (2014) showed that during the formation of vortices, filamentary high density dust structures can be formed, which may be relevant in molecular clouds. However, none of these studies analyse in detail the role to collisions between charges and neutrals.

The works by Jones and Downes (2011, 2012) are more related with the present discussion. Jones and Downes (2011) performed 2.5D simulations of the KHi in weakly ionized plasmas composed of ions, electrons, and neutrals. The governing equations solved by Jones and Downes (2011) consider the neutral fluid velocity as the centre-of-mass velocity, neglect the inertia of charges, and use a generalized Ohm’s law. Jones and Downes (2011) investigated the ambipolar-dominated and Hall-dominated regimes. They find in their simulations that partial ionization effects do not significantly influence the growth rate of the instability in the linear phase, which seem to be in contradiction with the linear results of the two-fluid model of Watson et al. (2004) and may be a consequence of the different models used by Watson et al. (2004) and Jones and Downes (2011). On the other hand, Jones and Downes (2011) obtain that ambipolar diffusion, which is caused by charge-neutral collisions, significantly decreases the growth of magnetic energy during the nonlinear phase of the instability compared with the evolution in ideal simulations, while the Hall effect, when dominant, causes the system to fail to settle to a quasi-steady state after saturation of the instability. Subsequently, Jones and Downes (2012) performed a specific application to molecular cloud conditions and also added the presence of dust. As in their previous paper, ambipolar diffusion leads to less amplification of the magnetic energy as the instability develops and, therefore, a stronger wind-up of the neutral fluid is possible during vortex formation. However, in molecular cloud conditions Hall’s term does not appear to have an important effect. As explained by Jones and Downes (2012), ambipolar diffusion largely suppresses the influence of Hall’s term. While the simulations of Jones and Downes (2011, 2012) are relevant in molecular clouds, their results cannot easily translated to the case of the partially ionized solar plasma, as solar prominences, where the weakly-ionized version of the governing equations is not applicable.

A thoughtful 2D numerical study of the nonlinear evolution of the KHi in a partially ionized medium has been performed by Hillier (2019). The considered configuration was the classic set-up: a sharp interface with a velocity shear between two partially ionized plasmas with different densities. A straight magnetic field was assumed parallel to the interface. A white noise perturbation was imposed initially. Hillier (2019) evolved in time the nonlinear two-fluid equations, explicitly solving the dynamics of both charges (ions) and neutrals, which were collisionaly coupled in the equations. By using a very high resolution, Hillier (2019) was able to study the different scales where the strength of the ion-neutral coupling varies. Hillier (2019) finds that at the small scales where neutrals are decoupled from ions, fully formed neutral vortices are present, whereas the ion velocity remains laminar in nature. However, and despite the very different velocity patters, similar density structures are found in both neutral and ion fluids because of the role of frictional heating and heat transfer between fluids (see Figure 4). At larger scales, neutrals become coupled with ions and the magnetic field, although secondary, smaller vortices can remain decoupled. The effect of the magnetic field may cause the nonlinear suppression of the instability even when the shear flow is linearly unstable (the same effect but in ideal MHD simulations was already discussed by Ryu et al., 2000). Hillier (2019) concludes that the turbulent energy cascade should probably show a transition from MHD turbulence at the coupled scales to neutral fluid HD turbulence at the smallest scales. This should be confirmed by future 2D and 3D high-resolution simulations.

Recently, 2D numerical simulations of the KHi in a partially ionized two-fluid plasma have been presented by Martínez-Gómez et al. (2021). These authors considered an initially unmagnetized medium composed of a dense partially ionized plasma slab embedded in a lighter environment. A constant flow longitudinal to the slab was assumed, with different flow velocities inside and outside the slab and a continuous transition between the two regions. The Biermann battery term was included in the Ohm’s Law. The simulations showed that the KHi grows after applying a small perturbation on the slab boundaries. During the KHi nonlinear development, the collisional interaction between the charges and neutrals drives the generation of a magnetic field through the Biermann battery mechanism (see, e.g., Kulsrud et al., 1997). The magnetic field is predominantly generated in the direction perpendicular to the flow, so it has no important influence on the development of the KHi. It is found that the strength of the generated magnetic field increases when the ionisation degree of the plasma decreases. Variations in the charges temperature are produced because of collisions with the neutrals. This increases the charges baroclinic term and enhances the Biermann battery effect (see detailed explainations in Martínez-Gómez et al., 2021). The results of the numerical experiments by Martínez-Gómez et al. (2021) are a good example of how interesting physical effects are driven by partial ionization.

Chandrasekhar (1961b) presents the simplest configuration in which the onset of the RTi can be studied (see also Sharp, 1984). Let us consider the case of two fully ionized unbounded plasmas with uniform densities laying one above the other in the presence of gravity. The two plasmas are separated by an abrupt interface and the plasma with density ρ ₁ is on top of the plasma with density ρ ₂, with ρ ₁ > ρ ₂. A straight and constant magnetic field that is orientated parallel to the interface, i.e., in the horizontal direction, permeates both plasmas. The direction of the magnetic field is assumed to be the same at both sides of the interface. The magnetic field strength at the two sides of the interface is denoted by B ₁ and B ₂, respectively. Using the linearized single-fluid ideal MHD equations and incompressible perturbations, we can perform a stability analysis of the interface by assuming a temporal dependence of the form   exp( − iωt), where ω is the frequency, and a Fourier-analysis in space. A dispersion relation for the surface waves propagating on the interface can be derived, whose analytic solution for the square of frequency yields, ω 2 = B 1 2 + B 2 2 μ ρ 1 + ρ 2 k ∥ 2 − ρ 1 − ρ 2 ρ 1 + ρ 2 g k , (35) where k _∥ is the component of the wavevector along the background magnetic field, k is the modulus of the wavevector, g is the acceleration of gravity, and the factor ρ 1 − ρ 2 / ρ 1 + ρ 2 is the Atwood number. Eq. 35 neglects fluid surface tension. The first term on the right-hand side of Eq. 35 is the square of the frequency of the incompressible surface Alfvén wave. The second term contains the effect of gravity and is related to the so-called interchange mode (see, e.g., Hillier, 2016). The RTi appears when ω ² < 0. This is always the case in the absence of a magnetic field, i.e., the HD case, so that for B ₁ = B ₂ = 0 the first term on the right-hand side of Eq. 35 vanishes. However, in the presence of a magnetic field, magnetic tension plays a stabilizing influence on the interface through the presence of a surface Alfvén wave, which is also called the undular mode in the literature. Perturbations whose wavevector is mostly parallel to the direction of the magnetic field are stable. Conversely, perturbations propagating with a sufficiently oblique angle remain unstable. The critical longitudinal wavenumber for stabilization due to magnetic tension is k ∥ , c r i t . = μ ρ 1 − ρ 2 B 1 2 + B 2 2 g k . (36)

We note that Eq. 36 depends on the modulus of the full wavenumber, k, so that the perpendicular lengthscale of the perturbations also influence the longitudinal lengthscales that are stable or unstable, making the stability properties of the interface to be rather complex. Readers are referred to Hillier (2016), who provides a more detailed account of the properties of the ideal RTi in its linear regime beyond the simple discussion provided here.

Terradas et al. (2012) investigated the triggering of the RTi in a model where the interface was replaced by a slab with a finite thickness (see Figure 5). The slab was meant to represent a thin thread of a solar prominence. They found that two different modes appear in the slab configuration, owing the presence of two separate interfaces. One mode is always stable and its character varies from being localised at the upper interface of the slab when the magnetic field is weak, to having a global nature and resembling the transverse kink mode of the whole slab when the magnetic field is strong. On the other hand, there is another mode that is unstable and localised at the lower interface when the magnetic field is weak, but it becomes a stable sausage magnetic mode when the magnetic field is increased. The criterion to know whether the magnetic field is weak or strong comes, again, from the comparison of the gravity force with the magnetic tension force. Ruderman et al. (2014) investigated the RTi in the presence of a sheared magnetic field (see also Hillier, 2016). They considered both the interface and slab scenarios. They showed that magnetic shear can have a strong effect on the growth rates of the instability. For small shear angles the RTi growth rate is linearly proportional to the shear angle, while in the limit of large angles the growth rate becomes independent of the shear angle.

As in the case of the KHi, the effect of compressibility on the onset of the RTi has been long under debate in both HD and MHD (see. e.g., Vandervoort, 1961; Shivamoggi, 1982; Bernstein and Book, 1983; Ribeyre et al., 2004; Livescu, 2004; Liberatore and Bouquet, 2008, among others). The complexity of the dispersion relation in the compressible case makes it necessary to resort to numerical solutions unless some heavy simplifications are made (see, e.g., Shivamoggi, 2012). Studies by, e.g., Ribeyre et al. (2004); Livescu (2004); Liberatore et al. (2009); Díaz et al. (2012); Ruderman (2017) point out that compressibility has a complex influence, so that the linear growing rates in the compressible case are smaller or larger than those obtained in incompressible approximation depending on the wavenumber range and the plasma β. However, the critical longitudinal wavenumber for magnetic tension stabilization (Eq. 36) appears to be unaltered by compressibility. Detailed studies of the RTi triggering and evolution in more realistic configurations necessarily require the use of numerical simulations (see Section 4.3).

As in the case of the KHi, the presence of a neutral species may have a relevant effect on the RTi. Neutrals are not affected by the stabilizing influence of magnetic tension. Therefore, a PIP composed of both charges and neutrals should be, in principle, more unstable regarding the RTi than a fully ionized plasma. Of course, ion-neutral collisions couple neutrals and charges, so that neutrals do feel the magnetic field influence in an indirect way. Again, the relevant question is whether collisions are able to fully stabilize the neutral fluid when the classic stability threshold is verified (Eq. 36) or, on the contrary, neutrals remain unstable.

The work by Hans (1968) is probably the first study of the role of ion-neutral collisions on the onset of the RTi in a PIP. Hans (1968) considered the classic setup of a denser fluid located on top of a lighter fluid and separated by an abrupt interface in the presence of a horizontal magnetic field. Both fluids are partially ionized. This author considered an incompressible two-fluid model that included collisions between charges and neutrals and finite Larmor radius effects expressed through a gyro-viscosity term for the charges. However, the gas pressure term was absent from the neutrals momentum equation and the induction equation only contained the ideal term. After performing a linear stability analysis, Hans (1968) concludes that the resulting dispersion relation has always an unstable root, whose growth rate decreases as the collision frequency increases. Ogbonna and Bhatia (1984) revisited the same configuration of Hans (1968) but performed a more in-depth paramenter study. In agreement with the previous results, Ogbonna and Bhatia (1984) conclude that friction between charges and neutrals decrease the growth rate of the RTi compared to the purely HD case, but an instability always remains. The problem was explored once again by Chhajlani and Vaghela (1989), who also included surface tension in their analysis (see Chandrasekhar, 1961b). They found that the conditions for the occurrence of the RTi, i.e., the instability threshold because of surface tension, remain unaffected by the presence of neutrals, but the growth rate is reduced by collisions.

All the works cited in the above paragraph ignored the role of neutral gas pressure and only included the ion-neutral coupling term in the neutrals momentum equation. Such a simplification forces the dynamics of neutrals to be entirely dependent to that of charges. The first study where the dynamics of neutrals was consistently described was done by Díaz et al. (2012). They considered a two-fluid model with collisions between charges and neutrals and included the gas pressure term in the neutrals momentum equation, but neglected all nonideal terms in the induction equation for simplicity. Again, the simple interface configuration was adopted. Díaz et al. (2012) considered both the incompressible and compressible cases. Shadmehri et al. (2013) performed a similar study to that of Díaz et al. (2012) but only in the incompressible limit. The dispersion relation derived in the incompressible case (Eq. 47 of Díaz et al. (2012) and Eq. 23 of Shadmehri et al. (2013)) can be written with a similar structure to that of the KHi dispersion relation (see Eq. 30), i.e., a dispersion relation associated with the neutral fluid multiplied by a dispersion relation associated with the charges, plus an additional term that couples both relations. The coupling term owes its existence to collisions. Unfortunately, unlike in the KHi case, the full dispersion relation is a rather complicated expression even in the incompressible limit, so that no simple analytic approximation for the growth rate was obtained by Díaz et al. (2012) and a numerical study of the solutions was necessary. Díaz et al. (2012) found two unstable modes that are related to the RTi in the neutral and charged fluids separately, with the neutral unstable mode having a larger growth rate. Díaz et al. (2012) computed their results as function of their parameter ϒ = ν _in L/c _s , where ν _in is ion-neutral collision frequency, L is a length scale, and c _s is the sound speed. Hence, large values of ϒ correspond to strong coupling. Díaz et al. (2012) found that collisions between ions and neutrals decrease the growth rates as ϒ increases (see Figure 6). The charges unstable mode can be stabilized by magnetic tension for a sufficiently large longitudinal wavenumber that is independent from the collision frequency (Eq. 36). However, the neutrals-related solution remains always unstable for finite ϒ. So, the results of the more complete model of Díaz et al. (2012) confirmed the previous tentative conclusions obtained in more simplified works (Hans, 1968; Ogbonna and Bhatia, 1984; Chhajlani and Vaghela, 1989): collisional coupling with the charges is not able to fully suppress the RTi of the neutral fluid even when the longitudinal wavenumber is larger than the critical one.

Díaz et al. (2012) extended their results to the compressible case and showed that the stability thresholds are not modified by compressibility: the mode related with the neutral fluid is always unstable. Compressibility appears to have a largely stabilizing effect, since smaller growth rates are obtained in the compressible case compared to those in the incompressible limit. However, the exact influence of compressibility depends upon the values of some parameters as, e.g., the density contrast and the ratio of the perpendicular to longitudinal wavenumbers. Therefore, although compressibility does not seem to have a relevant influence on the actual triggering of the RTi, it may be relevant for the subsequent growth of the perturbations towards the nonlinear development of the instability.

Díaz et al. (2014) revisited the same problem studied in Díaz et al. (2012) but considering the single-fluid approximation instead of the two-fluid approach. In the single-fluid model, the role of the ion-neutral collisions remains in the form of the ambipolar diffusion term of the induction equation. Regarding the properties of the RTi, the results of the single-fluid model consistently agreed to those previously obtained in the two-fluid formalism. The main difference between the single-fluid and two-fluid cases is that in the single-fluid case it is not possible to disentangle the nature of the two unstable modes, since the single-fluid approximation already assumes a strong coupling between charges and neutrals, so that only one unstable “global” mode is found. Díaz et al. (2014) also considered other nonideal terms in the induction equation as Ohm’s diffusion, Hall’s term, and the battery term that were ignored in Díaz et al. (2012). However, for typical conditions of solar prominences, Díaz et al. (2014) found that the ambipolar term is by far the dominant term in the induction equation, while the other terms are largely negligible, thus confirming the appropriateness of the ideal induction equation used in the two-fluid model of Díaz et al. (2012). Later, Ruderman et al. (2018) expanded the work of Díaz et al. (2014) and performed a detailed mathematical study of the effect of magnetic shear. Ruderman et al. (2018) found that the larger the shear angle, the smaller the maximum RTi growth rate. In addition, Ruderman et al. (2018) concluded that ambipolar diffusion only affects the growth rate when the plasma β is small.

Astrophysical applications of the linear theory discussed above include solar prominences and local clouds. Díaz et al. (2012) argue that the obtained growth rates for physical conditions in partially ionized prominences may explain the existence of fine structures with lifetimes of the order of 30 min, while the timescales derived from the classical theory for fully ionized plasma are about one order of magnitude shorter and incompatible with the observed lifetimes. In addition, Díaz et al. (2014) conclude that their partially ionized model provides an instability timescale comparable to observed lifetimes of RTi plumes in prominences. On the other hand, Shadmehri et al. (2013) explained that, owing to two-fluid effects, the RTi may operate less effectively in local clouds than previously though according to classical theory.

Valuable but limited information is provided by the analysis of the linear regime of the RTi. The linear theoretical analysis must necessarily be complemented with nonlinear numerical simulations to understand the later evolution of the instability. An extensive literature on the nonlinear evolution of the RTi exists. Here we only discuss a few representative works.

Early attempts to simulate with low resolutions the development of the RTi in ideal MHD include, e.g., the works by Wang and Robertson (1985) and Jun et al. (1995), among other relevant papers. Wang and Robertson (1985) considered 2D compressible simulations with the goal to study the mixing process occurring at later stages of the instability. Their simulations show how the interface is deformed by the formation of mushroom-like structures (fingers and plumes) that ultimately lead to the presence of swirling motions and plasma mixing. An energy cascade towards small scales occurs during the nonlinear phase. An inverse cascade towards large scales is found later in the evolution owing to the merging of smaller structures. Jun et al. (1995) performed quasi-incompressible simulations in 2D and 3D. The quasi-incompressible regime was achieved by considering a high gas pressure. They found that turbulence associated with the nonlinear evolution of the RTi amplifies the magnetic field, preferentially on small scales. This effect is more important in 3D than in 2D. The growth of fingers and plumes display the generation of vortex structures associated with secondary KHi. The fully developed structures appear to be sensitive to the initial magnetic field strength and orientation. The growth of instabilities across the magnetic field is more prominent, a result already anticipated by the theoretical linear analysis.

Subsequent 3D simulations with higher resolutions provided more in-depth information of the small-scale evolution of the perturbations. Isobe et al. (2005) performed 3D simulations of magnetic flux emergence in the solar atmosphere in which the RTi is shown to develop. The instability leads to the formation of filamentary structure and small-scale current sheets. Dissipation of the current sheets produces plasma heating, which may me relevant in the solar corona. The simulations also indicate that magnetic reconnection is initiated locally by the RTi, while in turn the reconnection process affects the growth of the instability, resulting in a spatially intermittent reconnection. Thus, the simulations of Isobe et al. (2005) clearly show that reconnection is an essential ingredient of the nonlinear evolution of the RTi.

Stone and Gardiner (2007) performed a rather general numerical study of the RTi in 3D by considering different initial configurations for the magnetic field. They show that uniform magnetic fields cannot suppress the instability in 3D. As linear theory predicts, interchange modes perpendicular to the field can grow at the same rate as in HD, while the magnetic field only stabilizes modes with large longitudinal wavenumbers. In the nonlinear simulations, this results in a highly anisotropic structure as the RTi evolves, as already anticipated by, e.g., Jun et al. (1995). The magnetic tension can inhibit secondary instabilities and reduce the growth of small scales and the mixing of plasmas. As a consequence of the restrained turbulent mixing, the fingers and plumes associated with the primary RTi can grow faster. However, a sheared magnetic field can significantly delay the instability and modify the structures that are formed in the full nonlinear regime. Later results by Ruderman et al. (2014) in the linear regime seem to support the importance of magnetic shear for the RTi growth rate. Stone and Gardiner (2007) used their results to explain the morphology of the optical filaments observed in the Crab nebula. However, more recent adaptive mesh refinement simulations by Porth et al. (2014) indicate that with very high resolutions, the filamentary structure driven by the RTi becomes less similar to the one observed in the Crab nebula.

Due to its observational importance (see Berger et al., 2008; Ryutova et al., 2010; Berger et al., 2010, 2011, 2017), the nonlinear development of the RTi in solar prominences has been intensely studied in the last decade. Also in the solar context, Moschou et al. (2015) showed indications of the RTi or interchange instability during numerical simulations or coronal rain formation and downfall. From here on, we shall put the focus of our discussion to the case of simulations in solar prominences. Again, we refer readers to Hillier (2018) for a comprehensive review on RTi simulations in prominences. What follows is a brief summary of some relevant results.

Hillier et al. (2011) performed the first attempt to simulate the RTi in a quiescent prominence. They considered a modification of the classic Kippenhahn-Schlüter model (Kippenhahn and Schlüter, 1957) in which a high-temperature, low-density tube was placed in the center of the model. The RTi drives upflows that interact with each other and create larger plumes through an inverse cascade process. The upflows advect the magnetic field lines through the structure, but the field line curvature is not heavily affected. The dynamics of the simulations compares qualitatively well with the observations, but the obtained upflow velocities in the model are lower than those observed.

Hillier et al. (2012a,b) revisited the same Kippenhahn-Schlüter configuration and performed an improved analysis of the simulations. Hillier et al. (2012a) carried out a detailed parameter study. They showed that the instability creates low-density filaments inside the prominence that are aligned with the direction of the magnetic field. This implies that a 3D RTi mode, i.e., a mode with both parallel and perpendicular components of the wavenumber, grows inside the prominence. For some range of parameters and initial conditions, the velocity and width of the simulated upflows can match the observed values. Subsequently, Hillier et al. (2012b) investigated a particular aspect of the simulations related to the process of magnetic reconnection. The RTi fingers and plumes and their associated shear flows resulted in the formation of current sheets that can reconnect. This reconnection allows the formation of downflowing blobs that may be related to the observed knots in prominences (see, e.g., Chae, 2010).

The Kippenhahn-Schlüter model used in Hillier et al. (2011, 2012a,b) can be understood as a local model of the prominence. Other works have studied the RTi in prominences considering more global models. Using a 3D arcade model for the prominence magnetic field, Terradas et al. (2015a) showed that the RTi can also grow in fully detached prominences suspended above the photosphere. The chromosphere was not included in the model. The magnetic field structure was anchored at the photosphere by means of the line-tying condition. The downflows developed in the model of Terradas et al. (2015a) may or may not reach the photosphere depending on the strength of the magnetic field in the prominence. For sufficiently intense magnetic fields, an oscillatory behavior of the flow is obtained, suggesting that the magnetic field is nonlinearly suppressing the RTi. The photospheric line-tying may also play a role. Magnetic shear is able to reduce or even to suppress the RTi completely for the considered spatial resolution. Indeed, further simulations by Terradas et al. (2016) using a twisted flux rope model display a much more stable behaviour than the arcade model regarding the RTi.

The development of the RTi in a whole-prominence model was also explored in the numerical simualtions by Keppens et al. (2015). They performed ideal MHD simulations at high resolution with an adaptive mesh refinement code. A horizontal background magnetic field nonuniform in the vertical direction was assumed, which introduces a local magnetic shear inside the prominence. The considered vertical profiles for the density and temperature aimed to represent a suspended prominence mass above the photosphere and chromosphere. The simulations show the quick generation of nonlinear magnetoconvective motions. The downflows associated with the RTi impact and reflect on the chromosphere, so that chromospheric plasma gets mixed with the prominence plasma suspended high in the corona. The results of Keppens et al. (2015) suggest that the nonlinear RTi may be intimately involved in the mass cycle of prominences and the interchange of material between prominences and the chromosphere.

An extension of the work by Keppens et al. (2015) was done in Xia and Keppens (2016). The main difference with the previous work was that the prominence was assumed to be formed by two parallel slab-like layers instead of a single monolitic mass. The simulations show that the two layers of the prominence evolve coherently due to their magnetic connectivity, since the mainly horizontal magnetic field transversely crosses both slabs. The RTi similarly evolves in the two layers following an equivalent dynamics to that discussed in Keppens et al. (2015). A detail of the strongly nonlinear evolution can be seen in Figure 7. Interestingly, the vertical density structures formed during the nonlinear evolution of the RTi may appear as horizontal thread-like structures when seen from the top of the prominence, which may explain the different structures observed in prominences above the solar limb and filaments on the solar disk.

All the papers cited above performed MHD simulations and none considered the effect of partial ionization. Hillier et al. (2010) investigated the evolution of the Kippenhahn-Schlüter prominence model under the presence of Cowling’s diffusion. However, they were only concerned with the diffusion of the magnetic field and their simplified model neglected all coupling with the plasma dynamics. A pioneering work where the role of partial ionization effects on RTi simulations was explored was done by Khomenko et al. (2014b), who considered the single-fluid equations for a PIP including the ambipolar diffusion term. Khomenko et al. (2014b) performed 2.5D simulations of the RTi initiated at the corona–prominence interface with a constant magnetic field perpendicular or almost perpendicular to the plane. They found that the configuration is always unstable, which agrees with the expected behavior based on linear theory (Díaz et al., 2012; Díaz et al., 2014). The growth rate of the small-scale modes in the non-linear regime is up to 50% larger than that obtained in equivalent ideal MHD simulations. Significant ion-neutral drift occurs at the corona–prominence interface (drift differences between ions and neutrals have been detected in prominences, see Khomenko et al., 2016). A faster downward motion of the neutral component with respect to the ionized component is obtained (see also Terradas et al., 2015b). The differences in temperature of the RTi bubbles between the ideal and ambipolar cases can be as large as 30% because of the additional heating associated with the dissipation of perpendicular currents by ambipolar diffusion (see Figure 8). The results of Khomenko et al. (2014b) clearly show that partial ionization effects have a measurable influence of the RTi onset and evolution in prominences, particularly in the small scales.

Recently, Popescu Braileanu et al. (2021a), Popescu Braileanu et al. (2021b) have extended the single-fluid simulations of Khomenko et al. (2014b) to the two-fluid case. Although their simulations remained 2.5D as in Khomenko et al. (2014b), the background model was improved by considering a smooth interface instead of a sharp transition and by including magnetic shear. The two-fluid simulations took into account viscosity, thermal conduction, ionization/recombination, and energy and momentum transfer through collisions between neutrals and charges. In Popescu Braileanu et al. (2021a) the study focused on assessing the effects of a smooth interface and magnetic shear. They showed that magnetic shear reduces or even suppresses the instability growth rate, as previous theoretical studies in the single-fluid approximation already anticipated (see, e.g., Ruderman et al., 2014). In turn, the inclusion of a continuous transition affects the length scales of the perturbations that develop due to the RTi, especially for wavelengths comparable to the density gradient length scale. In Popescu Braileanu et al. (2021b) the emphasis was put in investigating the effects of collisions within the framework of the two-fluid model. For prominence conditions, ionization and recombination do not significantly influence the development of the RTi main structures. Secondary structures formed during the later nonlinear development and mixing seem to be more affected by ionization and recombination processes. Ion-neutral collisions play a role in determining the evolution and dissipation of small structures during the nonlinear stage of the RTi. Ion-neutral decoupling affects smaller scales than viscosity. The nonlinear development of the RTi drives decoupling between neutrals and charges. The decoupling is more pronounced on small spatial scales and at locations of strong gradients in density and/or magnetic field. The ion-neutral flow decoupling is more pronounced in the horizontal direction. So, compared with the single-fluid simulations of Khomenko et al. (2014b), the two-fluid simulations of Popescu Braileanu et al. (2021a), Popescu Braileanu et al. (2021b) display a more complex and rich behaviour at small scales, where two-fluid effects become of relevance and the ion-neutral decoupling takes place (see Figure 9). Another important difference between the models of Khomenko et al. (2014b) and Popescu Braileanu et al. (2021a), Popescu Braileanu et al. (2021b) is that the latter papers considered a sheared magnetic configuration. For sufficiently large shear, the cutoff in the linear growth rate disappears, so that the dissipative effect of the collisions can be better observed at small scales.

In the solar context, and related with prominences, thermal modes were studied in detail in non-homogeneous and fully ionized plasmas (e.g., van der Linden and Goossens, 1991b,a; van der Linden, 1993; Ireland et al., 1995; Carbonell et al., 2004; Soler et al., 2011, 2012a, and references therein). Furthermore, Hildner (1974); Oran et al. (1982); Dahlburg and Mariska (1988); Karpen et al. (1989b,a); Cargill and Hood (1989); Carbonell et al. (2004) investigated both linear and nonlinear thermal instabilities since this instability represents a key guide to understand the formation of prominences and coronal rain in the solar corona. Effects like thermal conduction across magnetic field lines and resistivity, could play an important role in the development of prominence fine structure (van der Linden, 1993; Ireland et al., 1998). Field’s analysis has been recently revised in a thorough study made by Waters and Proga (2019) who, together with a discussion about the stability of acoustic and condensation modes between the two wavelengths corresponding to isobaric and isochoric instabilities, included a numerical analysis of the non-linear evolution of condensation modes with applications to active galactic nuclei (AGNs) and giant molecular clouds (GCM). In solar corona conditions, Claes and Keppens (2019) have also revisited Field’s treatment and have combined it with numerical simulations in order to check the predicted growth rates and to study the non-linear regime. Claes et al. (2020) have investigated the formation of threads by non-linear thermal instability, obtaining that the threads are misaligned with the underlying magnetic structure. This study reveals intricate multidimensional processes that occur through in situ condensations in a low plasma β regime representative of the coronal medium. Falle et al. (2020) have also done another re-analysis of Field’s results using a new approach based on Whitham’s theory of wave hierarchies (Whitham, 1974), which simplifies calculations and establishes relationships between dispersion relations and the involved physical processes. On the other hand, nonlinear thermal instabilities in magnetized solar plasmas have also been studied by Dahlburg et al. (1987); Karpen et al. (1988, 1989b) and the nonlinear dynamics of radiative condensations in the case of optically thin plasmas was reviewed by Meerson (1996). Readers are referred to the cited papers for details on the physics of the instability.

In the case of fully ionized astrophysical plasmas beyond the solar case, there is an extense literature about the role of thermal instability in the interstellar medium and star formation. As a matter of example (Gomez-Pelaez and Moreno-Insertis, 2002), considered the linear thermal stability of an expanding and cooling medium with thermal conduction and self-gravity included, and the obtained results were applied to the case of a hot optically thin interstellar medium such as supernova remnants. Wareing et al. (2016); Ji et al. (2018); Fragile et al. (2018); Choi and Stone (2012) showed that the thermal instability mechanism can be responsible for the formation of molecular loops in the galactic central region, the appearance of higher-density filaments in star-forming clouds, the cold gas found in galactic haloes, the vertical collapse of the accretion disk around stellar-mass black holes, and that, together with anisotropic thermal conduction, thermal instability can significantly affect the shapes and sizes of cold clouds of the interstellar medium. For a thorough review on thermal instability in the interstellar medium, see Inutsuka et al. (2005) and references therein. These examples point out the relevance of the thermal instability in the astrophysical context. The importance of all these studies rests in the fact that thermal instability can explain the formation of dense and cool localized structures in astrophysical and laboratory plasmas, when their masses are less than those required for gravitational contraction i.e. smaller than the Jean’s mass. (See Section 8).

The study of thermal instability in partially ionized plasmas is a topic of great interest in astrophysics, however, the first studies were performed with laboratory plasmas. The first demonstration of thermal instability in fluids giving place to thermal convection was done by Benard (1900). The theoretical interpretation was made by Rayleigh who established that for any Rayleigh number, R, greater than a critical one, R _c , thermal instability sets in. Chandrasekhar (1961a) studied the thermal instability of a layer of fluid in which an adverse temperature gradient is maintained by heating the layer from below, and different situations such as the Bénard problem, the presence of rotation and/or magnetic field were considered. Following the approach by Chandrasekhar (1961a), Sharma (1976) studied the collisional effects between ions and neutrals on the thermal instability of a horizontal layer of finite thickness made of a partially ionized plasma, assumed incompressible, and permeated by a vertical magnetic field. This layer is heated from below in such a way that a steady temperature gradient is maintained. The boundaries of the layer were assumed to be free and the limiting fluid was considered non-conducting. After linearizing a simplified version of the two-fluid equations, including thermal conduction, viscosity, resistivity, gravity, and collisions between neutrals and charges, a dispersion relation was obtained that was expressed in terms of the Rayleigh number. For stationary convection, it was shown that collisions do not have any effect on thermal instability. The possibility that the instability appeared as an overstability was examined and conditions for the non-existence of overstability were obtained. Furthermore, the effect of the magnetic field on the thermal instability was studied by considering the variation of the Rayleigh number with the Chandrasekhar number, which includes the magnetic field. It was found that the magnetic field has an stabilizing effect on thermal instability. Later, Sharma and Sharma (1978), Sharma and Sharma (1985) extended the previous study to the case in which the fluid is rotating and arrived at the same conclusions for the case of stationary convection. Regarding overstability, they studied the variation of the Rayleigh number with respect to the Taylor number, which includes rotation, and with respect to the Chandrasekhar number, finding that rotation and magnetic field have an stabilizing effect on thermal instability. Next, rotation was removed and the Hall term was taken into account in the calculations. As before, for stationary convection collisions have no effect on thermal instability in a partially ionized Hall plasma, however, they found that the Hall currents have a destabilizing effect on thermal hydromagnetic instability when the variation of the Rayleigh number with respect to the Hall parameter was studied. In summary, the obtained results are similar to those of Chandrasekhar (1961a) since collisions between ions and neutrals do not have any influence on stationary convection. Finally, Sharma and Misra (1986) considered the same configuration as before but assuming that the fluid is compressible, and studied the effect of compressibility on thermal instability. Following a linear analysis, they found that the main effect of compressibility is to delay the onset of the instability, therefore, compressibility has an stabilizing influence.

On the other hand, observations have pointed out that small-scale structures pervade the interstellar medium (ISM). In order to understand the formation of these small-scale structures in the weakly ionized cold HI and molecular clouds of the ISM, Fukue and Kamaya (2007) made a thorough study of thermal instability in a magnetized and weakly ionized plasma, based on the idea that thermal instability could be a mechanism for further fragmentation in molecular clouds. The reason is that even in systems stable against gravitational instability, the systems can become thermally unstable because the critical wavelength for thermal instability is smaller than that for gravitational (Jeans) instability, which means that the size of the critical length is smaller than Jean’s length. Using the two-fluid approach (see Subsection 2.2), Fukue and Kamaya (2007) studied the effect of ion-neutral friction on the growth of thermal instability, focusing in the condensation mode rather than in the oscillatory modes. After linearizing the two-fluid equations, a sixth order dispersion relation was obtained, which describes two sets of three modes: one of the modes in each set is the condensation mode and the other two are oscillatory modes. From the dispersion relation it can be seen that the critical wavelength for the condensation mode, which is obtained by setting the growth rate equal to zero, is not affected by ion-neutral friction. In order to study the importance of the condensation mode for structure formation, it is necessary to solve numerically the full dispersion relation, and this was done for different values of the strength of magnetic field and friction. In absence of magnetic field and friction, two independent pure modes for ion and neutrals are present [see thin curves in Figure 2 of Fukue and Kamaya (2007)], however, when friction is allowed pure ion or neutral modes do not exist and, for both components, the growth rate at larger scales decreases because friction is more efficient at these scales [see thick curves in Figure 2 of Fukue and Kamaya (2007)]. This figure also shows that, when friction is allowed, the curves representing the behaviour of the modes originate from the curves corresponding to the case in which magnetic field and friction are not considered [see thin curves in Figure 2 of Fukue and Kamaya (2007)]. When magnetic field is taken into account, the critical wavelength for the ion mode is enlarged, as the criterion for thermal instability in a fully ionized magnetized plasma indicates. However, the critical wavelength for the neutral mode is not affected [see Figure 3 of Fukue and Kamaya (2007)]. When the strength of the magnetic field is strongly increased, the magnetic field completely stabilizes the ion mode, and only the neutral mode appears in the thermal instability diagram [see Figure 5 of Fukue and Kamaya (2007)]. When the strength of ion-neutral friction increases, ion mode remains stabilized for sufficiently intense magnetic fields, while the neutral mode remains unstable but with a reduced growth rate [see Figure 6 of Fukue and Kamaya (2007)]. Baruah et al. (2010) investigated thermal instability in a weakly ionized plasma with ionization and recombination such as in the envelopes of planetary nebulae that surround red giant stars. The presence of ionization and recombination affects the instability behaviour and limits the size of the small structures which can be formed by this instability.

As it is well-known, instability of thermal modes play an important role in the solar context, in particular in solar prominences or coronal rain, since it provides a mechanism for plasma condensation in the solar corona. Following the instability criteria introduced by Field (1965), it is straightforward to conclude that the radiative loss function is of great importance to determine the stability of thermal modes. However, an accurate description of radiative losses in cool plasmas, like prominences, is not an easy task because solutions of the radiative transfer problem in non-local thermodynamic equilibrium are needed. Taking into account partial ionization effects, Soler et al. (2012a) studied the stability of thermal modes in unbounded and uniform plasmas with physical properties akin to those of solar prominences. The single-fluid approximation was used. They considered three different parametrizations for the radiative loss function: the function proposed by Hildner (1974), the Raymond-Klimchuk function (Klimchuk and Cargill, 2001), and a radiative loss function derived from the CHIANTI atomic database (Parenti et al., 2006; Parenti and Vial, 2007). Starting from single-fluid equations (see Subsection 2.1) with ambipolar diffusion and non-adiabatic terms included, a linear analysis was performed in order to obtain a dispersion relation that describes two damped slow modes, two damped fast modes, and one thermal mode. From the full dispersion relation, and after performing a first-order expansion for a low-β plasma, an approximate dispersion relation for slow and thermal modes can be obtained. Since for the thermal mode in a static plasma the real part of the frequency is equal to zero, we can write ω = is, therefore, due to the temporal dependence exp( − iωt) assumed by Soler et al. (2012a), thermal mode perturbations are proportional to exp(is), where s is the growth rate. An approximate growth rate was obtained, which depends on the non-adiabatic terms and Cowling’s diffusion. When Cowling’s diffusion is set to zero, the growth rate of the thermal mode in the case of a fully ionized plasma is recovered. Next, the full dispersion relation was solved numerically. Figure 10 shows a plot of the growth rate versus wavelength for the case of fully and partially ionized plasma and different temperatures and ionization degrees. In both plots, the approximate and the full numerical solution of the dispersion relation are compared. For the considered temperatures, it is found that the approximate and full numerical solutions are in good agreement which points out that magnetic diffusion plays a negligible role in the value of the growth rate. On the other hand, Figure 10 also shows that for fully or partially ionized plasma the thermal mode is stabilized at short wavelengths which is due to the fact that thermal conduction is stronger for short wavelengths. On the opposite, for long wavelengths the growth rate saturates becoming independent from the wavelength. Furthermore, the growth rate decreases with the temperature. When Cowling’s diffusion and thermal conduction by neutrals are taken into account together, we can observe that the growth rate increases when the ionization degree decreases. This is due to an increase of the effective density of the plasma because, in the constant pressure calculations of Soler et al. (2012a), the number density of neutrals is increased when the ionization degree decreases. In addition, when the ionization degree decreases the critical wavelength increases due to the effect of thermal conduction by neutrals. As in the case of fully ionized plasma, the effect of magnetic diffusion on the growth rate is negligible.

Finally, the effect of thermal instability on the onset of the gravitational or Jeans instability, described in subsect 8, has been also considered.