© ESO, 2014

1. Introduction

One of the well-known fluid instabilities that has been applied widely in many different astrophysical contexts is the Rayleigh-Taylor instability (RTI), which appears when a lighter fluid supports or accelerates a heavier one. We can cite as examples: the RTI in planetary nebulas (Bucciantini et al. 2004), supernova explosions (Fryxell et al. 1991), acretion disks (Wang & Nepveu 1983), relativistic jets (Matsumoto & Masada 2013), the evolution of the inner layers of red giants (Eggleton et al. 2006), the formation of hydrogen clouds in the local bubble (Breitschwerdt et al. 2000), the solar atmospheric flux tubes (Parker 1979), and prominences (Isobe et al. 2005). These topics are among the vast literature regarding this instability in astrophysical plasmas.

The starting point of these studies is the classical hydrodynamic instability with the addition of a magnetic field. To study the Rayleigh-Taylor instability in magnetohydrodynamics (MHD), it is better to use the divergence-free velocity condition (which is the type of movement more likely to produce instabilities, since the energy is not wasted in compressing the plasma). Hence, the dispersion relation for the two fluids with densities ρ₁ and ρ₂ that lay one below the other with a magnetic field parallel to the contact surface is (Chandrasekhar 1961; Priest 1982), $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{gk}\frac{{\mathit{\rho }}_{\mathrm{2}}\mathrm{-}{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{1}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}}{\mathit{\mu }\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}\mathrm{)}}\mathit{,}& & \end{array}$(1)where B₀ is the magnetic field strength, k is the modulus of the wavenumber, k_x is the component of the wavenumber along the magnetic field direction, k_y the component across the field, g is the acceleration of gravity, and μ is the magnetic permittivity. The hydromagnetic case is recovered if no magnetic field is considered, which means that the configuration is unstable if ρ₂>ρ₁. The magnetic field stabilizes perturbations down to a critical value of k_x along its direction, but the field cannot stabilize perturbations across the field (k_x = 0) no matter how strong it might be.

Other effects, such as compressibility, viscosity, tension forces, relativistic corrections, magnetic fields perpendicular to the surface or other types of stratification, and forces may introduce corrections to the classical stability criterion and linear growth rate of the RTI. Compressibility is normally the leading effect, which must be taken into account in most of the applications. There have been many studies of the effect of compressibility in the magnetic RTI (see for example Vandervoort 1961; Shivamoggi 1982; Bernstein & Book 1983; Ribeyre et al. 2004; Livescu 2004; Shivamoggi 2008; Liberatore & Bouquet 2008, and references therein). These studies show that the compressibility has mainly a stabilizing influence by lowering the linear growth rates, although the stability threshold that appears in Eq. (1) is not modified. On the other hand, by considering only this extra effect greatly complicates the solution, so simple relations, such as Eq. (1), are no longer obtained.

Prominences are a very likely candidate to display the RTI in the solar atmosphere, since they are composed of cool and dense plasma surrounded by much lighter coronal plasma (see the reviews Labrosse et al. 2010; Mackay et al. 2010). The magnetic field plays a key role in the structure and dynamics of prominences; hence, the RTI must be studied in the context of the MHD theory. Using computational techniques allow us to directly solve the MHD partial differential equations and study the non-linear regime of the RTI. These numerical studies of RTI (see for example Jun et al. 1995; Arber et al. 2007; Stone & Gardiner 2007) agree with the results of the linear theory and show some new interesting features, such as the growth enhancement of bubbles and fingers in the non-linear regime across the field, since it prevents secondary Kelvin-Helmholtz instabilities and mixing between the fluids. This non-linear phase agrees qualitatively with the observations of turbulent plumes and bubbles in prominences (Isobe et al. 2005; Berger et al. 2008, 2010, 2011; Heinzel et al. 2008; Ryutova et al. 2010) and has even been used to infer plasma properties from observational features (Hillier et al. 2012b). The presence of the RTI has also been confirmed numerically in more general prominence models with 3D geometry (Isobe et al. 2006; Hillier et al. 2011, 2012a,c).

Prominences also have another feature that can be relevant for the RTI because of their physical properties, namely partial ionization. Since prominences are relatively cool and dense objects, their plasma is expected to be partially ionized (Patsourakos & Vial 2002; Gilbert et al. 2007; Labrosse et al. 2010; Zaqarashvili et al. 2011b). The ionization degree of the prominence plasma has not been directly measured with accuracy but the plasma is neither in a completely ionized state or is a neutral gas with the typical physical parameters of this objects. This partially ionized prominence plasma can no longer be described with a one-fluid ideal MHD. Chhajlani & Vaghela (1989) studied the magnetic RTI in a two-fluid model without considering the full dynamics of the neutral fluid (which was only subject to collisions) and concluded that the stability threshold is not changed, but the growth rate is lowered. A two-fluid description with one neutral and one ion-electron fluid coupled only by ion-neutral collision (without electron collisions and diffusive terms in the induction equation) was considered in Soler et al. (2012b) for the Kelvin-Helmholtz instability due to shear flow at an interface between two partially ionized plasmas and for the RTI in a similar setup in Díaz et al. (2012). The conclusion in Díaz et al. (2012) was that the instability threshold of the RTI is modified, so the configuration with heavier fluid on the top of lighter fluid is always unstable. However, the linear growth rate is substantially reduced depending on the parameters. This approach has also been used by Shadmehri et al. (2013) in the context of the RTI in the borders of the local bubble, where partial ionization also plays a relevant role.

Our main objective in this paper is to study the stability threshold and the linear growth rate of the RTI in a plasma described in a one-fluid model. This takes into account the partial ionization effects in the form of a generalized induction equation and considers all the types of collisions, forces present, and the diffusive terms in the induction equation. We develop a general formulation, having in mind its global applicability to different astrophysical situations. After this general formulation has been set, we particularize it to the case of solar prominences and consider the importance of the different terms in the induction equation (many of them neglected in previous works) for the development of the RTI in this environment. We compare our results with those from the simple two-fluid approach in Díaz et al. (2012) and try to understand the linear phase of the RTI in a partially ionized plasma to use it as a comparison with more complicated computational models with complex geometry and the non-linear stages.

2. Multi-fluid equations for partially ionized plasma

2.1. Fluid equations for each species

The general transport equations for a multi-component plasma can be derived from the Boltzmann kinetic equation, taking into acount general properties of the collisions terms (Braginskii 1965; Bittencourt 1986; Balescu 1988). The most common form of the MHD theory can only be applied to totally ionized plasma, where the different species are completely coupled dynamically and thermally by collisions, so partial ionization cannot be described. The first extension from MHD in the presence of neutrals for strong collisional coupling is to only consider the modifications due to collisions in the generalized Ohm’s law and energy transport (Braginskii 1965; Khodachenko et al. 2004; Forteza et al. 2007), by assuming a strong thermal coupling and neglecting the transport coefficients. Another way of including partial ionization effects is to use a multi-fluid treatment, where ions and electrons are considered together as an ion-electron fluid (due to their strong electromagnetic coupling) and neutrals are separately considered with as many neutral species as one wishes to include (see, e.g., Zaqarashvili et al. 2011b,a; Soler et al. 2012a). The two-fluid approach (a hydrogen plasma only considered) was used in Díaz et al. (2012) for studying the RTI with the additional neglection of other diffusive terms in the induction equation. Here, we take the single fluid approach by combining the equations for each species in a single fluid equation and obtaining a generalized form of the MHD equations valid for PI plasmas.

We proceed with the derivation of the generalized one-fluid equations starting from the macroscopic equations for each species. In the following expressions, the subscripts i, n, and e stand for ions, neutrals, and electrons, respectively. The momentum equation for each species is $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{e}}\left(\frac{\mathit{\partial }{{v}}_{\mathrm{e}}}{\mathit{\partial t}}\mathrm{+}{{v}}_{\mathrm{e}}\mathrm{·}\mathrm{\nabla }{{v}}_{\mathrm{e}}\right)& \mathrm{=}& \mathrm{-}\mathrm{\nabla }{\mathit{p}}_{\mathrm{e}}\mathrm{-}\mathit{e}{\mathit{n}}_{\mathrm{e}}\mathrm{(}{E}\mathrm{+}{{v}}_{\mathrm{e}}\mathrm{\times }{{B}}^{\mathrm{)}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}{g}\mathrm{+}{{R}}_{\mathrm{e}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{i}}\left(\frac{\mathit{\partial }{{v}}_{\mathrm{i}}}{\mathit{\partial t}}\mathrm{+}{{v}}_{\mathrm{i}}\mathrm{·}\mathrm{\nabla }{{v}}_{\mathrm{i}}\right)& \mathrm{=}& \mathrm{-}\mathrm{\nabla }{\mathit{p}}_{\mathrm{i}}\mathrm{+}\mathit{e}{\mathit{n}}_{\mathrm{i}}\mathrm{(}{E}\mathrm{+}{{v}}_{\mathrm{i}}\mathrm{\times }{{B}}^{\mathrm{)}}\mathrm{+}{\mathit{\rho }}_{\mathrm{i}}{g}\mathrm{+}{{R}}_{\mathrm{i}}\mathit{,}\\ {\mathit{\rho }}_{\mathrm{n}}\left(\frac{\mathit{\partial }{{v}}_{\mathrm{n}}}{\mathit{\partial t}}\mathrm{+}{{v}}_{\mathrm{n}}\mathrm{·}\mathrm{\nabla }{{v}}_{\mathrm{n}}\right)& \mathrm{=}& \end{array}$(2)where v_i, v_e, and v_n are the velocity of the ion, electron, and neutral fluid, respectively; p_i, p_e, and p_n are the pressure of the ion, electron and the neutral fluid, respectively; ρ_i, ρ_e and ρ_n are the ion, electron and neutron densities, respectively; and R_i, R_e, and R_n are the momentum transfer terms due to collisions for ions, electrons, and neutrals, respectively. We also define the total density ρ = ρ_n + ρ_i, the neutral fraction ξ_n = ρ_n/ρ, and the ion fraction ξ_i = ρ_i/ρ, with ξ_n + ξ_i = 1. Hence, the parameter ξ_n indicates the ionization degree from ξ_n = 0 for a fully ionized plasma to ξ_n = 1 for a neutral gas. We have assumed the non-diagonal terms of the pressure tensors to be negligible and the diagonal terms to be equal, so the pressure tensor becomes isotropic and can be represented by the scalar pressure. The elastic collision term of each species is approximated as (Braginskii 1965; Bittencourt 1986) $\begin{array}{ccc}{{R}}_{\mathit{\alpha }}\mathrm{=}\mathrm{-}{\mathit{\rho }}_{\mathit{\alpha }}\sum _{\mathit{\beta }}{\mathit{\nu }}_{\mathit{\alpha \beta }}\mathrm{(}{{v}}_{\mathit{\alpha }}\mathrm{-}{{v}}_{\mathit{\beta }}\mathrm{)}\mathit{,}& & \end{array}$(3)where ν_αβ is the collisional frequency of species α with particles of species β. For our three species, they become $\begin{array}{ccc}{{R}}_{\mathrm{e}}& \mathrm{=}& \mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{ei}}\mathrm{\right(}{{v}}_{\mathrm{e}}\mathrm{-}{{v}}_{\mathrm{i}}\mathrm{)}\mathrm{+}{\mathit{\nu }}_{\mathrm{en}}\mathrm{(}{{v}}_{\mathrm{e}}\mathrm{-}{{v}}_{\mathrm{n}}\mathrm{\left)}\mathrm{\right)}\\ & \mathrm{=}& {\mathit{\rho }}_{\mathrm{e}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{J}\mathrm{(}{\mathit{\nu }}_{\mathrm{ei}}\mathrm{+}{\mathit{\nu }}_{\mathrm{en}}\mathrm{)}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}{\mathit{\nu }}_{\mathrm{en}}{w}\\ {{R}}_{\mathrm{i}}& \mathrm{=}& \mathrm{-}{\mathit{\rho }}_{\mathrm{i}}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{ie}}\mathrm{\right(}{{v}}_{\mathrm{i}}\mathrm{-}{{v}}_{\mathrm{e}}\mathrm{)}\mathrm{+}{\mathit{\nu }}_{\mathrm{in}}\mathrm{(}{{v}}_{\mathrm{i}}\mathrm{-}{{v}}_{\mathrm{n}}\mathrm{\left)}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{-}{\mathit{\rho }}_{\mathrm{i}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{\mathit{\nu }}_{\mathrm{ie}}{J}\mathrm{-}{\mathit{\rho }}_{\mathrm{i}}{\mathit{\nu }}_{\mathrm{in}}{w}\\ {{R}}_{\mathrm{n}}& \mathrm{=}& \mathrm{-}{\mathit{\rho }}_{\mathrm{n}}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{ni}}\mathrm{\right(}{{v}}_{\mathrm{n}}\mathrm{-}{{v}}_{\mathrm{i}}\mathrm{)}\mathrm{+}{\mathit{\nu }}_{\mathrm{ne}}\mathrm{(}{{v}}_{\mathrm{n}}\mathrm{-}{{v}}_{\mathrm{e}}\mathrm{\left)}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{-}{\mathit{\rho }}_{\mathrm{n}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{\mathit{\nu }}_{\mathrm{ne}}{J}\mathrm{-}{\mathit{\rho }}_{\mathrm{n}}\mathrm{(}{\mathit{\nu }}_{\mathrm{ni}}\mathrm{+}{\mathit{\nu }}_{\mathrm{ne}}\mathrm{)}{w}\mathit{,}\end{array}$(4)with w = v_i − v_n as the diffusion velocity of ions with respect to neutrals and J = en_e(v_i − v_e) as the total current density. This form of the friction momentum transfer does not alter the one-fluid equations of mass conservation and momentum conservation from their MHD counterparts when written in terms of the overall velocity of the plasma, v = ∑ _{α = i,e,i}(ρ_αv_α) /ρ ≈ ξ_iv_i + ξ_nv_n. However, the energy equation has to take into account the currents arising from these non-MHD terms, and an additional equation for the evolution of the magnetic field, namely a generalized induction equation, is also required to close the system.

2.2. Equation for the diffusion velocity

To obtain the momentum equation for the diffusion velocity w, we proceed as follows. The momentum equation for electrons and ions are added up, neglecting the electron inertial terms because of the small electron mass compared to the other species. We do not neglect the electron gravity compared to the ion gravity yet (ρ_eg compared to ρ_ig), $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{i}}(\frac{\mathit{\partial }{{v}}_{\mathrm{i}}}{\mathit{\partial t}}\mathrm{+}{{v}}_{\mathrm{i}}\mathrm{·}\mathrm{\nabla }{{v}}_{\mathrm{i}})& \mathrm{=}& {J}\mathrm{\times }{B}\mathrm{+}\mathrm{(}{\mathit{\rho }}_{\mathrm{i}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}{g}\mathrm{-}\mathrm{\nabla }\mathrm{(}{\mathit{p}}_{\mathrm{i}}\mathrm{+}{\mathit{p}}_{\mathrm{e}}\mathrm{)}\\ & & \mathrm{-}{\mathit{\alpha }}_{\mathrm{n}}{w}\mathrm{+}{\mathit{\rho }}_{\mathrm{n}}{\mathit{\nu }}_{\mathrm{ne}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{J}\\ {\mathit{\rho }}_{\mathrm{n}}(\frac{\mathit{\partial }{{v}}_{\mathrm{n}}}{\mathit{\partial t}}\mathrm{+}{{v}}_{\mathrm{n}}\mathrm{·}\mathrm{\nabla }{{v}}_{\mathrm{n}})& \mathrm{=}& {\mathit{\rho }}_{\mathrm{n}}{g}\mathrm{-}\mathrm{\nabla }{\mathit{p}}_{\mathrm{n}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{n}}{w}\mathrm{-}{\mathit{\rho }}_{\mathrm{n}}{\mathit{\nu }}_{\mathrm{ne}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{J}\mathit{,}\end{array}$(5)with the definition for the coefficient of friction between the plasma and the neutral gas as $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{n}}\mathrm{=}{\mathit{\rho }}_{\mathrm{n}}\mathrm{(}{\mathit{\nu }}_{\mathrm{ni}}\mathrm{+}{\mathit{\nu }}_{\mathrm{ne}}\mathrm{)}\mathrm{=}{\mathit{\rho }}_{\mathrm{i}}{\mathit{\nu }}_{\mathrm{in}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}{\mathit{\nu }}_{\mathrm{en}}\mathit{.}& & \end{array}$(6)Now, we add the upper equation multiplied by ξ_n and lower equation, multiplied by − ξ_i. The result is: $\begin{array}{ccc}{\mathit{\xi }}_{\mathrm{i}}{\mathit{\xi }}_{\mathrm{n}}\mathit{\rho }\left(\frac{{\mathrm{D}}_{\mathrm{i}}{{v}}_{\mathrm{i}}}{\mathrm{D}\mathit{t}}\mathrm{-}\frac{{\mathrm{D}}_{\mathrm{n}}{{v}}_{\mathrm{n}}}{\mathrm{D}\mathit{t}}\right)\mathrm{=}{\mathit{\xi }}_{\mathrm{n}}\left[{J}\mathrm{\times }{B}\right]\mathrm{-}{G}\mathrm{-}{\mathit{\alpha }}_{\mathrm{n}}{w}& & \\ \mathrm{+}{\mathit{\rho }}_{\mathrm{n}}{\mathit{\nu }}_{\mathrm{ne}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{J}\mathrm{+}{\mathit{\xi }}_{\mathrm{n}}{\mathit{\rho }}_{\mathrm{e}}{g}\mathit{,}& & \end{array}$(7)where $\begin{array}{ccc}\frac{{\mathrm{D}}_{\mathit{\alpha }}}{\mathrm{D}\mathit{t}}\mathrm{=}\frac{\mathit{\partial }}{{\mathit{\partial }}_{\mathit{t}}}\mathrm{+}{\mathit{v}}_{\mathit{\alpha }}\mathrm{·}\mathrm{\nabla }\mathit{.}& & \end{array}$(8)We have taken into account that ξ_i + ξ_n = 1, thus the gravity terms for ions and neutrals cancel out and the electron term might become relevant. This gravitional term is similar to those of the electron inertia that have been already neglected, but we keep it here to check its magnitude, especially since the RTI is driven by gravity. It is therefore important to assure its effect as much as possible. Regarding the pressure gradients, we have introduced the new PI pressure terms following Braginskii (1965) as $\begin{array}{ccc}{G}\mathrm{=}{\mathit{\xi }}_{\mathrm{n}}\mathrm{\nabla }\mathrm{(}{\mathit{p}}_{\mathrm{i}}\mathrm{+}{\mathit{p}}_{\mathrm{e}}\mathrm{)}\mathrm{-}{\mathit{\xi }}_{\mathrm{i}}\mathrm{\nabla }{\mathit{p}}_{\mathrm{n}}\mathit{.}& & \end{array}$(9)We still have the ion and neutral inertia terms in Eq. (7). These are neglected on the basis of the following argument. We can express the total derivative in terms of the diffusion velocity $\begin{array}{ccc}\left(\frac{{\mathrm{D}}_{\mathrm{i}}{{v}}_{\mathrm{i}}}{\mathrm{D}\mathit{t}}\mathrm{-}\frac{{\mathrm{D}}_{\mathrm{n}}{{v}}_{\mathrm{n}}}{\mathrm{D}\mathit{t}}\right)& \mathrm{=}& \frac{\mathit{\partial }{w}}{\mathit{\partial t}}\mathrm{+}\mathrm{(}{\mathit{\xi }}_{\mathrm{n}}\mathrm{-}{\mathit{\xi }}_{\mathrm{i}}\mathrm{)}{w}\mathrm{·}\mathrm{\nabla }{w}\\ & & \mathrm{+}\hspace{0.17em}{w}\mathrm{·}\mathrm{\nabla }{v}\mathrm{+}{v}\mathrm{·}\mathrm{\nabla }{w}\mathrm{\approx }\frac{\mathit{\partial }{w}}{\mathit{\partial t}}\mathrm{·}\end{array}$(10)In a linear regime, all the advection terms in this equation are second order effects, and the remaining time derivative of w can be neglected when compared with the friction terms, which are of the order of w/τ_col with τ_col ~ 1/ν_αβ being the characteristic timescale related to collisions.

Taking all the aforementioned simplifications into account, we obtain an expression for the diffusion velocity between ions and neutrals: $\begin{array}{ccc}{w}\mathrm{=}\frac{{\mathit{\xi }}_{\mathrm{n}}}{{\mathit{\alpha }}_{\mathrm{n}}}{J}\mathrm{\times }{B}\mathrm{-}\frac{{G}}{{\mathit{\alpha }}_{\mathrm{n}}}\mathrm{+}{\mathit{\rho }}_{\mathrm{n}}{\mathit{\nu }}_{\mathrm{ne}}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}{\mathit{\alpha }}_{\mathrm{n}}}{J}\mathrm{+}\frac{{\mathit{\xi }}_{\mathrm{n}}{\mathit{\rho }}_{\mathrm{e}}}{{\mathit{\alpha }}_{\mathrm{n}}}{g}\mathit{.}& & \end{array}$(11)From this expression, we can see that the ion and neutral fluids do not follow each other exactly, which raises additional dissipative effects. On the other hand, by neglecting the inertial terms, we have obtained an explicit expression for the diffusion velocity in terms of other variables, assuming that collisions lead to the terminal values of w given by Eq. (11). These are much faster than the evolution of the remaining variables. Thus, the velocities of the ion and neutral species do not longer appear in the equations and the diffusion velocity w can be computed from the single-fluid variables.

2.3. Induction equation

To proceed further, we need a generalized Ohm’s law and an induction equation to obtain an equation for the magnetic field evolution. These are obtained from the momentum equation for electrons in Eq. (2), again neglecting their inertial terms. After some algebra we obtain $\begin{array}{ccc}{E}\mathrm{+}{v}\mathrm{\times }{B}& \mathrm{=}& \mathrm{-}{\mathit{\xi }}_{\mathrm{n}}\hspace{0.17em}{w}\mathrm{\times }{B}\mathrm{+}\frac{\mathrm{1}}{{\mathit{n}}_{\mathrm{i}}\mathit{e}}{J}\mathrm{\times }{B}\\ & & \end{array}$(12)with the definition α_e = ρ_e(ν_ei + ν_en). We then substitute the expression for the diffusion velocity in Eq. (11) and insert the result in Faraday’s law, obtaining $\begin{array}{ccc}\frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& {\nabla }\mathrm{\times }[\mathrm{(}{v}\mathrm{\times }{B}\mathrm{)}\mathrm{-}\frac{{J}}{\mathit{\sigma }}\mathrm{-}\left(\frac{\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\epsilon }{\mathit{\xi }}_{\mathrm{n}}}{\mathit{e}{\mathit{n}}_{\mathrm{e}}}{J}\mathrm{\times }{B}\right)\\ & & \mathrm{+}\left(\frac{{\mathit{\xi }}_{\mathrm{n}}^{\mathrm{2}}}{{\mathit{\alpha }}_{\mathrm{n}}}\mathrm{(}{J}\mathrm{\times }{B}\mathrm{)}\mathrm{\times }{B}\right)\mathrm{-}\left(\frac{\mathit{\epsilon }{G}\mathrm{-}\mathrm{\nabla }{\mathit{p}}_{\mathrm{e}}}{\mathit{e}{\mathit{n}}_{\mathrm{e}}}\right)\mathrm{-}{\left(\frac{{\mathit{\xi }}_{\mathrm{n}}}{{\mathit{\alpha }}_{\mathrm{n}}}{G}\mathrm{\times }{B}\right)}^{}\\ & & \end{array}$(13)with the definitions of ε = ρ_eν_en/α_n (which is a small parameter) and the Ohmic conductivity σ = (en_e)²/(α_e − ε²α_n). In this equation, the terms on the right hand side are: ideal MHD induction term, Ohmic term, Hall term, ambipolar term, generalized battery term (which includes a part already present in fully ionized plasmas and a part depending on partial ionization by means of G), a G × B term of perpendicular currents caused by these pressure gradients (similar to the Hall term with currents), and gravity terms (again with a similar g × B part also included). We can define the coefficients in the different terms as $\begin{array}{ccc}\mathit{\eta }& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\sigma \mu }}\mathrm{=}\frac{{\mathit{\alpha }}_{\mathrm{e}}\mathrm{-}\mathrm{(}{\mathit{\rho }}_{\mathrm{e}}{\mathit{\nu }}_{\mathrm{en}}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{\alpha }}_{\mathrm{n}}}{\mathrm{(}\mathit{e}{\mathit{n}}_{\mathrm{e}}{\mathrm{)}}^{\mathrm{2}}\mathit{\mu }}\mathit{,}\\ {\mathit{\eta }}_{\mathrm{H}}& \mathrm{=}& \frac{\mathrm{1}\mathrm{-}\mathrm{2}\mathit{\epsilon }{\mathit{\xi }}_{\mathrm{n}}}{\mathit{e}{\mathit{n}}_{\mathrm{e}}\mathit{\mu }}{\mathit{B}}_{\mathrm{0}}\mathit{,}\\ {\mathit{\eta }}_{\mathrm{A}}& \mathrm{=}& \frac{{\mathit{\xi }}_{\mathrm{n}}^{\mathrm{2}}}{{\mathit{\alpha }}_{\mathrm{n}}\mathit{\mu }}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{,}\\ {\mathit{\chi }}_{\mathrm{p}}& \mathrm{=}& \frac{{\mathit{\xi }}_{\mathrm{n}}}{{\mathit{\alpha }}_{\mathrm{n}}}\mathit{,}\\ {\mathit{\chi }}_{\mathrm{g}}& \mathrm{=}& \end{array}$(14)with η, η_H, η_A, χ_p, and χ_g being the ohmic diffusivity, Hall diffusivity, ambipolar diffusivity, and coefficients related to the battery and gravity terms, respectively. Assuming that the acceleration of gravity is uniform, the curl of the second to last term in Eq. (13) vanishes, and the induction equation is finally written as $\begin{array}{ccc}\frac{\mathit{\partial }{B}}{\mathit{\partial t}}& \mathrm{=}& {\nabla }\mathrm{\times }[\mathrm{(}{v}\mathrm{\times }{B}\mathrm{)}\mathrm{-}\mathit{\eta }\mathrm{\nabla }\mathrm{\times }{B}\mathrm{-}{\mathit{\eta }}_{\mathrm{H}}\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\mathit{/}{\mathit{B}}_{\mathrm{0}}\\ & & \mathrm{+}{\mathit{\eta }}_{\mathrm{A}}\left\{\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\right\}\mathrm{\times }{B}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}\left(\mathit{\epsilon }{G}\mathrm{-}\mathrm{\nabla }{\mathit{p}}_{\mathrm{e}}\right)\mathit{/}\mathrm{(}\mathit{e}{\mathit{n}}_{\mathrm{e}}{\mathrm{)}}^{}\\ & & \end{array}$(15)in which we used Ampere’s law (neglecting Maxwell’s displacement current) to eliminate the current density in terms of the magnetic field, J = ∇ × B/μ. Equation (15) is a very general form of the induction equation in the one-fluid description of partial ionized plasmas, and it is a generalization of the well-known generalized induction equation in classical textbooks (see for example Braginskii 1965) with all the pressure gradient, gravity terms, and the expressions of the diffusion coefficients.

The formulation of the induction equation in Eq. (15) allows for very general analyses that can be useful in a broad context of astrophysical plasmas. Some of the terms are a priory expected to be smaller than others (as those related to the electron mass), but others cannot be ruled out just from general considerations. We discuss their importance for a case with parameter values appropriate for solar prominences below. We describe next the plasma and magnetic field configuration used to study the RTI in this environment in Sect. 3 and then explore the effect of the leading term under these conditions in Sect. 4, namely, with the ambipolar diffusion ${\mathit{\eta }}_{\mathrm{A}}\left\{\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\right\}\mathrm{\times }{B}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$. Then, we consider the full induction equation in Sect. 5 to test the magnitude of the remaining terms and finish discussing the results and drawing our conclusions.

3. Reference configuration

Since we aim to obtain some extensions to the well-known formula in Eq. (1), we restrict the analysis to the simple configuration of a contact surface following the classical analysis in Chandrasekhar (1961), Drazin & Reid (1981), Priest (1982), which is amenable to analytical solutions. We use this configuration to study the RTI in prominence threads, especially for choosing the values of the equilibrium and perturbation parameters, but the method developed in this paper is general and can be applied to other astrophysical situations, which involve the RTI in PI plasmas.

The reference configuration consists of two regions filled with uniform plasmas composed of ions, electrons, and neutrals separated by a contact surface at z = 0. We use Cartesian coordinates and denote the quantities in the plasma below the discontinuity (z< 0) with a subscript 1 and those in the plasma above the discontinuity (z> 0) with a subscript 2. The magnetic field permeating the plasma is uniform and tangent to the discontinuity, so . Since gravity is perpendicular to it, the gravity direction is g =  −gẑ. The whole configuration is invariant in the x and y-directions.

Fig. 1 Sketch of the equilibrium configuration used in the analysis of this work. The equilibrium state is a contact surface between two regions filled uniformly with plasma having different properties with the lower quantities labelled as “1” and the upper ones as “2”. The magnetic field is uniform and directed along the x-axis, while the whole configuration is invariant in the x- and y-directions.

In absence of hydrostatic pressure gradients or flows, the plasma described in the previous paragraph is not in equilibrium, since nothing counteracts the gravity force. We are not interested in the overall equilibrium, only in the local region where the instability is triggered. More precisely, the equilibrium pressure gradient ∇p₀ is related to the gravitational scale height, while the perturbed quantities vary on a much shorter spatial scale. Hence, we assume that all the plasma magnitudes (namely the density, the pressure, and the ionization degree) are constant in each zone. Pressure balance along the discontinuity demands that the total pressure for each species must be equal in each side, and since the magnetic field is assumed to be uniform, this means the gas pressures are equal in each side (p₁ = p₂). On the other hand, the temperature, density, and ionization degree on each region are parameters of our model. Since we are interested in studying the RTI, we assume that ρ₂>ρ₁. No ionization-recombination processes are included, so the ionization degree in each region remains constant. Note also that the plasma beta $\mathit{\beta }\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{2}}^{\mathrm{2}}\mathit{/}{\mathit{c}}_{\mathrm{A}\mathrm{2}}^{\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{1}}^{\mathrm{2}}\mathit{/}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma \mu }{\mathit{p}}_{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ is then constant all over the domain.

Another important simplification in this particular configuration is that we can neglect the variations in the coefficients in Eq. (14) during the evolution of the instability. The equilibrium field satisfies ∇ × B₀ = 0, so there are no currents in the reference state and the only contribution of the diffusive terms to the first-order induction equation are those with the coefficients calculated in the reference configuration.

Finally, using this reference configuration, the classical instability criterion (Eq. (1)) can be written in the following form $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{2}\mathrm{(}{\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{2}{\cos }^{\mathrm{2}}\mathit{\theta }\mathrm{-}{{\mathit{c}}_{\mathrm{crit}}^{\mathrm{2}}}^{\mathrm{)}}\hspace{0.17em}{\mathit{k}}^{\mathrm{2}}& & \end{array}$(16)with θ being the angle between the equilibrium magnetic field and the wavevector k and k the wavevector modulus (wavenumber). We define the critical speed as $\begin{array}{ccc}{\mathit{c}}_{\mathrm{crit}}\mathrm{=}{\left(\frac{\mathit{g}}{\mathrm{2}\mathit{k}}\hspace{0.17em}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{2}}\mathrm{-}{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{2}}\mathrm{+}{\mathit{\rho }}_{\mathrm{1}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}& & \end{array}$(17)and the reduced square Alfvén speed as $\begin{array}{ccc}{\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{=}{\left(\frac{{\mathit{B}}_{\mathrm{0}}}{{\mathit{\mu }}_{\mathrm{0}}\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}\mathrm{)}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}& & \end{array}$(18)with the usual definition for the squared Alfvén speed ${\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{=}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathit{\mu }{\mathit{\rho }}_{\mathrm{0}}\mathrm{\right)}$, so ${\overline{\mathit{c}}}_{\mathrm{A}}{}^{-2}\mathrm{=}\mathrm{(}{\mathit{c}}_{\mathrm{A}\mathrm{1}}{\mathrm{)}}^{-2}\mathrm{+}\mathrm{(}{\mathit{c}}_{\mathrm{A}\mathrm{2}}{\mathrm{)}}^{-2}$. It is convenient then to use ${\overline{\mathit{c}}}_{\mathrm{A}}$ as a parameter. Notice that we have ${\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{\approx }{\mathit{c}}_{\mathrm{A}\mathrm{2}}$ in the case of a prominence with ρ₂ ≫ ρ₁, so this averaged Alfvén speed is approximately the Alfvén speed in the prominence.

4. MHD and ambipolar diffusion

It is clear that dealing with the all the terms in Eq. (15) is very difficult. Hence, we first concentrate in the modifications introduced in the ideal MHD theory by the ambipolar term, which has proved to be relevant in solar atmospheric situations (see for example Khodachenko et al. 2004; Arber et al. 2007; Khomenko & Collados 2012, and references therein). We neglect all the magnetic diffusion terms in this section except the ambipolar one, obtaining a simple form of the induction equation, $\begin{array}{ccc}\frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\left({v}\mathrm{\times }{B}\mathrm{+}{\mathit{\eta }}_{\mathrm{A}}\left\{\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\right\}\mathrm{\times }{B}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}\right)\mathit{.}& & \end{array}$(19)The mass and momentum conservation equations are not modified by the presence of ambipolar diffusion. In addition, we assume an adiabatic energy equation plus the contribution from the ambipolar diffusion term (neglecting transport terms such as conduction, radiation of other heating sources). Hence, after deriving the energy term corresponding to the ambipolar diffusion, our system of basic equations is $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathrm{·}\left(\mathit{\rho }{v}\right)\mathit{,}\\ & & \mathit{\rho }\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathit{\rho }{v}\mathrm{·}\mathrm{\nabla }{v}\mathrm{-}\mathrm{\nabla }\mathit{p}\mathrm{+}\frac{\mathrm{1}}{\mathit{\mu }}\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\mathrm{+}\mathit{\rho }{g}\mathit{,}\\ & & \frac{\mathit{\partial }{B}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\mathrm{(}{v}\mathrm{\times }{B}\mathrm{+}{\mathit{\eta }}_{\mathrm{A}}\left\{\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\right\}\mathrm{\times }{B}\mathit{/}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}\\ & & \frac{\mathit{\partial p}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{v}\mathrm{·}\mathrm{\nabla }\mathit{p}\mathrm{-}\mathit{\gamma p}\mathrm{\nabla }\mathrm{·}{v}\\ & & \mathrm{+}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\frac{{\mathit{\eta }}_{\mathrm{A}}}{\mathit{\mu }{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }\left[\left\{\left(\mathrm{\nabla }\mathrm{\times }{B}\right)\mathrm{\times }{B}\right\}\mathrm{\times }{B}\right]\mathit{,}\end{array}$(20)where γ is the adiabatic index and p = p_i + p_e + p_n the total scalar pressure of the fluid.

4.1. Linearized equations

Next, we study linear perturbations from the uniform state. To obtain a general formulation, we label the reference quantities with the subscript 0 and the linear perturbations without subscript (B = B₀x + b). The subscript 0 can be replaced with 1 of 2 when one of the regions in which the physical domain is considered, but otherwise, the deduction is valid for any uniform configuration.

Since no equilibrium flow is present (v₀ = 0), all the advection terms are second order effects. No currents are present in the reference state, so the term in the energy equation coming from the ambipolar diffusion is also a second order effect and can be neglected in the linealized problem. Hence, we are left with a considerably simpler system of differential equations, namely

$\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{=}\mathrm{-}{\mathit{\rho }}_{\mathrm{0}}\mathrm{\nabla }\mathrm{·}{v}\mathit{,}\\ & & {\mathit{\rho }}_{\mathrm{0}}\frac{\mathit{\partial }{v}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathit{p}\mathrm{+}\frac{\mathrm{1}}{\mathit{\mu }}\left(\mathrm{\nabla }\mathrm{\times }{b}\right)\mathrm{\times }{{B}}_{\mathrm{0}}\mathrm{+}\mathit{\rho }{g}\mathit{,}\\ & & \frac{\mathit{\partial }{b}}{\mathit{\partial t}}\mathrm{=}\mathrm{\nabla }\mathrm{\times }\left({v}\mathrm{\times }{{B}}_{\mathrm{0}}\mathrm{+}\frac{{\mathit{\eta }}_{\mathrm{A}}}{{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}}\left\{\left(\mathrm{\nabla }\mathrm{\times }{b}\right)\mathrm{\times }{{B}}_{\mathrm{0}}\right\}\mathrm{\times }{{B}}_{\mathrm{0}}\right)\mathit{,}\\ & & \end{array}$(21)

This system of equations can be reduced to only two by taking the time derivative of the momentum equation and substituting the expressions for the density and pressure from the continuity and energy equations, respectively. We obtain a system of two partial differential equations for the perturbations in velocity and magnetic field, namely $\begin{array}{ccc}& & \frac{{\mathit{\partial }}^{\mathrm{2}}{v}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{=}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\left(\mathrm{\nabla }\mathrm{\times }\frac{\mathit{\partial }{b}}{\mathit{\partial t}}\right)\mathrm{\times }{{B}}_{\mathrm{0}}\mathrm{-}\left(\mathrm{\nabla }\mathrm{·}{v}\right){g}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\nabla }\left(\mathrm{\nabla }\mathrm{·}{v}\right)\mathit{,}\\ & & \end{array}$(22)We have defined the squared sound speed as ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}$. The equilibrium properties of the medium relevant for the RTI are included in the sound speed, the Alfvén speed, and the ambipolar diffusivity coefficient. In our problem, it is not possible to derive a single equation by eliminating the Lorentz force term in the motion equation by using the induction equation as is routinely done in ideal MHD (see for example Roberts 1981; Priest 1982; Goedbloed & Poedts 2004).

We are left with the ambipolar diffusivity η_A (Eq. (14) as the only parameter depending on the ionization fraction. This diffusivity is calculated in the equilibrium state and depends on the ionization degree ξ_n and the neutral friction coefficient (Eq. (6)), $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{n}}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}\mathrm{)}\mathrm{(}{\mathit{\nu }}_{\mathrm{in}}\mathrm{+}{\mathit{\nu }}_{\mathrm{en}}{\mathit{m}}_{\mathrm{e}}\mathit{/}{\mathit{m}}_{\mathrm{i}}\mathrm{)}\mathit{.}& & \end{array}$(23)We can neglect the term with the ratio of the electron and ion mass and use the expression for the ion-neutral collision frequency in a plasma (see, e.g., Braginskii 1965; Soler et al. 2009), $\begin{array}{ccc}{\mathit{\nu }}_{\mathrm{in}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{n}}}{{\mathit{m}}_{\mathrm{n}}}\sqrt{\frac{\mathrm{16}{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\pi }{\mathit{m}}_{\mathrm{n}}}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{in}}\mathit{,}& & \end{array}$(24)where T is the temperature, m_n the neutron mass, k_B is the Boltzmann constant, and σ_in ≈ 5 × 10^-19 m² is the collisional cross section for proton-hydrogen collisions (assuming a hydrogen plasma). Notice that this collision frequency is a theoretical value for hard-sphere collisions between protons and H molecules, while there are hints that actual values can differ from this simple calculation (Mitchner & Kruger 1973; Vranjes & Krstic 2013). A strong thermal coupling is assumed, so the temperature of the different species is the same. Hence, ${\mathit{k}}_{\mathrm{B}}\mathit{T}\mathit{/}{\mathit{m}}_{\mathrm{n}}\mathrm{=}{\mathit{p}}_{\mathrm{n}}\mathit{/}{\mathit{\rho }}_{\mathrm{n}}\mathrm{=}\mathrm{(}\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}{\mathrm{)}}^{-1}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathit{\gamma }$ and the expression for the friction coefficient is $\begin{array}{ccc}{\mathit{\alpha }}_{\mathrm{n}}\mathrm{=}\frac{\mathrm{4}{\mathit{\sigma }}_{\mathrm{in}}}{{\mathit{m}}_{\mathrm{n}}{\left(\mathit{\pi \gamma }\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}\frac{{\mathit{\xi }}_{\mathrm{n}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}\mathrm{)}}{\mathrm{(}\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}& & \end{array}$(25)Finally, we obtain an equation for η_A from Eq. (14) in terms of the equilibrium parameters in each region. $\begin{array}{ccc}{\mathit{\eta }}_{\mathrm{A}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{n}}{\left(\mathit{\pi \gamma }\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{4}{\mathit{\sigma }}_{\mathrm{in}}}\hspace{0.17em}\frac{{\mathit{\xi }}_{\mathrm{n}}\mathrm{(}\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}\hspace{0.17em}\frac{{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{0}}{\mathit{c}}_{\mathrm{s}}}\mathrm{·}& & \end{array}$(26)This expression depends on the medium density. We use ρ₀ = 10¹⁰ kg m^-3, a value representative of typical densities in prominences, for computing this coefficient in the prominence region through the paper, while the value in the corona is adjusted by taking into account the prominence-corona contrast ratio ρ₂/ρ₁ used in each calculation.

4.2. Normal mode analysis

We consider the normal mode decomposition and write the temporal dependence of the perturbation as e^{− iωt}. We use Fourier analysis in the spatial directions, where the medium is uniform, and write the perturbations as e^{ik_xx + ik_yy} with k_x and k_y as the wavenumbers in the x- and y-directions, respectively, and as the wavevector parallel to the surface.

Then, we combine Eqs. (22) and arrive at a system of two coupled equations for v_z, the z-component of the velocity (normal to the surface), and b_x, the x-component of the perturbation of the magnetic field (along the equilibrium magnetic field): $\begin{array}{ccc}& & \mathrm{i}\mathit{\omega }{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\frac{\mathrm{d}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}\mathit{z}}\mathrm{-}\frac{\mathrm{i}\mathit{g}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{\eta }}_{\mathrm{A}}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{)}}{{\mathit{\omega }}^{\mathrm{2}}\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{(}\mathit{\omega }{\mathit{\eta }}_{\mathrm{A}}\mathrm{-}\mathrm{i}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{)}}\hspace{0.17em}{\mathit{b}}_{\mathit{x}}\mathrm{=}\\ & & \mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{)}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{v}}_{\mathit{z}}}{\mathrm{d}{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\mathit{g\omega }\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{)}\frac{\mathrm{d}{\mathit{v}}_{\mathit{z}}}{\mathrm{d}\mathit{z}}\\ & & \\ & & \frac{\mathrm{d}{\mathit{v}}_{\mathit{z}}}{\mathrm{d}\mathit{z}}\mathrm{=}{\mathit{\eta }}_{\mathrm{A}}\mathrm{\{}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathit{\omega }\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}^{}\\ & & \hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}{\mathrm{)}}^{\mathrm{\}}}\mathit{/}\mathrm{\{}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\\ & & \hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{i}\mathit{\omega }{\mathit{\eta }}_{\mathrm{A}}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{\}}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}{\mathit{z}}^{\mathrm{2}}}\mathrm{+}\mathrm{\{}\mathrm{i}\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{)}\mathrm{[}\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\\ & & \hspace{0.17em}\mathrm{-}\hspace{0.17em}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathit{\omega }\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right(}\mathit{\omega }\mathrm{+}{\mathit{k}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{)}{\mathrm{]}}^{\mathrm{\}}}\mathit{/}\mathrm{\{}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{)}}^{}\\ & & \end{array}$We can further operate these equations to obtain a single differential equation, $\begin{array}{ccc}{\mathit{C}}_{\mathrm{4}}\frac{{\mathrm{d}}^{\mathrm{4}}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}{\mathit{z}}^{\mathrm{4}}}\mathrm{+}{\mathit{C}}_{\mathrm{3}}\frac{{\mathrm{d}}^{\mathrm{3}}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}{\mathit{z}}^{\mathrm{3}}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}{\mathit{z}}^{\mathrm{2}}}\mathrm{+}{\mathit{C}}_{\mathrm{1}}\frac{\mathrm{d}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}\mathit{z}}\mathrm{+}{\mathit{C}}_{\mathrm{0}}\hspace{0.17em}{\mathit{b}}_{\mathit{x}}\mathrm{=}\mathrm{0}& & \end{array}$(29)with the following definitions for the coefficients, $\begin{array}{ccc}{\mathit{C}}_{\mathrm{4}}& \mathrm{=}& \mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathit{,}\\ {\mathit{C}}_{\mathrm{3}}& \mathrm{=}& \mathrm{-}\mathit{g\omega }{\mathit{\eta }}_{\mathrm{A}}\mathit{,}\\ {\mathit{C}}_{\mathrm{2}}& \mathrm{=}& \mathrm{i}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathit{\omega }\mathrm{[}{\mathit{\omega }}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}\mathrm{i}\mathit{\omega }\mathrm{-}\mathrm{2}{\mathit{k}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}{\mathrm{)}}^{\mathrm{]}}\mathit{,}\\ {\mathit{C}}_{\mathrm{1}}& \mathrm{=}& \mathrm{-}\mathrm{i}\mathit{g\omega }\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{)}\mathit{,}\\ {\mathit{C}}_{\mathrm{0}}& \mathrm{=}& \mathrm{i}\mathrm{[}\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathit{\omega }\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}{\mathrm{)}}^{\mathrm{]}}\mathit{.}\end{array}$(30)This ordinary differential equation is valid in each zone with the equilibrium quantities ${\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}$, ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, and η_A. The subscripts 1 or 2 are applied to the two regions above and below the contact surface at z = 0, respectively.

Finally, we need the boundary conditions to match the solutions at the boundary z = 0. In ideal MHD, the continuity of the normal component of the velocity perturbation and the continuity of the total pressure are enough, along with the contribution of gravity from the momentum balance at the boundary. However, additional constraints are necessary in this particular problem. We derive these conditions by integrating Eqs. (21) across the surface z = 0 and finding the limit of infinitesimal integration volume. We obtain only four independent jump relations, namely $\begin{array}{ccc}& & \mathrm{[}{\mathit{\rho }}_{\mathrm{0}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}_{\mathit{z}}\mathrm{]}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \mathrm{[}{\mathit{\rho }}_{\mathrm{0}}\mathrm{(}\mathrm{i}\mathit{\omega }{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}{\mathit{b}}_{\mathit{x}}\mathrm{-}\mathrm{i}{\mathit{k}}_{\mathit{x}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}_{\mathit{x}}\mathrm{-}\mathrm{i}{\mathit{k}}_{\mathit{y}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}_{\mathit{y}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{v}}_{\mathit{z}}^{\mathrm{\prime }}\mathrm{-}\mathit{g}{\mathit{v}}_{\mathit{z}}{\mathrm{)}}^{\mathrm{]}}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \left[\mathrm{i}{\mathit{k}}_{\mathit{x}}{\mathit{\eta }}_{\mathrm{A}}{\mathit{b}}_{\mathit{z}}\mathrm{-}{\mathit{v}}_{\mathit{z}}\mathrm{+}{\mathit{\eta }}_{\mathrm{A}}{\mathit{b}}_{\mathit{x}}^{\mathrm{\prime }}\right]\mathrm{=}\mathrm{0}\mathit{,}\\ & & \end{array}$(31)where the prime represents a derivative on the z-direction and [ X] = X₂(0⁺) − X₁(0⁻) stands for the jump of the quantity X across z = 0. Expressing the components of the perturbed velocity in terms of b_x and v_z, we obtain the set of jump relations required for our system (since v_z requires the integral of b_x). The first relation is just the typical boundary condition for the velocity perturbation (since ${\mathit{\rho }}_{\mathrm{0}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}$ is equal in both sides due to the equilibrium pressure balance), and the second is related to the momentum balance, since the first term is related to the magnetic pressure, the next three to the gas pressure and the last one to the gravity force. The remaining two conditions come from the new terms from the induction equation, which do not give any additional information for η_A = 0. In terms of v_x and b_z, our final set of boundary conditions is $\begin{array}{ccc}& & \left[{\mathit{v}}_{\mathit{z}}\right]\mathrm{=}\mathrm{0}\mathit{,}\\ & & [\mathrm{i}\mathit{\omega }{\mathit{\rho }}_{\mathrm{0}}\frac{{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathit{\omega }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{(}\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}\mathrm{)}}{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\hspace{0.17em}{\mathit{b}}_{\mathit{x}}\\ & & \mathrm{+}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{0}}\mathit{g}{\mathit{v}}_{\mathit{z}}\mathrm{+}\frac{{\mathit{\omega }}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}{\mathit{b}}_{\mathit{x}}^{\mathrm{\prime \prime }}}{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}]\mathrm{=}\mathrm{0}\mathit{,}\\ & & \left[\frac{{\mathit{\eta }}_{\mathrm{A}}{\mathit{b}}_{\mathit{x}}^{\mathrm{\prime }}\mathrm{-}{\mathit{v}}_{\mathit{z}}}{\mathit{\omega }\mathrm{+}\mathrm{i}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{A}}}\right]\mathrm{=}\mathrm{0}\mathit{,}\\ & & \end{array}$(32)

We use the standard definition for the linear growth rate of the instability, Im(ω). Hence, Im(ω) > 0 is related to unstable modes, while Im(ω) < 0 marks a damping in the wave. We also define the density contrast ρ₂/ρ₁ and $\mathit{\beta }\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{2}}^{\mathrm{2}}\mathit{/}{\mathit{c}}_{\mathrm{A}\mathrm{2}}^{\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathrm{s}\mathrm{1}}^{\mathrm{2}}\mathit{/}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}\mathrm{=}\mathit{\gamma \mu }{\mathit{p}}_{\mathrm{2}}\mathit{/}{\mathit{B}}_{\mathrm{0}}^{\mathrm{2}}$ as a measure of the magnetic field strength compared with the pressure terms.

4.3. Fully ionized plasma

Before dealing with the full problem, we check the known limit of the ideal MHD. This is achieved by considering a fully ionized plasma and letting η_A → 0. In this case, the equations are highly simplified, and Eq. (29) just becomes $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\frac{\mathit{g}{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{\Omega }}\frac{\mathrm{d}{\mathit{b}}_{\mathit{x}}}{\mathrm{d}\mathit{z}}\mathrm{+}\frac{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}\mathrm{(}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}\mathrm{)}\mathrm{\Omega }}{\mathrm{\Omega }}\hspace{0.17em}{\mathit{b}}_{\mathit{x}}\mathrm{=}\mathrm{0}& & \end{array}$(33)with $\mathrm{\Omega }\mathrm{=}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$. We also obtain the relation $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{v}}_{\mathit{z}}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{i}\mathit{\omega }\frac{\mathit{g}{\mathit{\omega }}^{\mathrm{2}}{\mathit{b}}_{\mathit{x}}\mathrm{-}\mathit{\omega }\mathrm{\Omega }{\mathit{b}}_{\mathit{x}}^{\mathrm{\prime }}}{\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}}^{\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}}\mathrm{·}& & \end{array}$(34)This differential equation describes the propagation of ideal MHD modes. For example, inserting a solution b_x = Ae^ik_zz, we recover the well-known dispersion relation for the MHD fast and slow modes present when gravity is not taken into account. Moreover, the boundary conditions in Eqs. (32) reduce to the continuity of the normal component of the velocity perturbation and the continuity of the total pressure (plus a gravity term) by imposing η_A → 0 with the two extra conditions either identically vanishing or reducing to those; thus, they recover the ideal MHD boundary conditions.

We can study the linear phase regime of the compressible MHD RTI by solving Eq. (33). The solutions of its indicial (characteristic) equation obtained after setting b_x = A_inde^mz (with A_ind an arbitrary constant) are $\begin{array}{ccc}{\mathit{m}}_{\mathrm{\pm }}\mathrm{=}\frac{\mathit{g}{\mathit{\omega }}^{\mathrm{2}}}{\mathrm{2}\mathrm{\Omega }}\mathrm{\pm }{\left({\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}\mathrm{-}\frac{{\mathit{\omega }}^{\mathrm{4}}}{\mathrm{\Omega }}\mathrm{+}\frac{{\mathit{g}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{4}}}{\mathrm{4}{\mathrm{\Omega }}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}& & \end{array}$(35)We need to choose the solution in each region that guarantees the perturbation to vanish far from the discontinuity. Hence, the solution is $\begin{array}{ccc}{\mathit{b}}_{\mathit{x}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(36)with A₁ and A₂ being constants. Applying the remaining boundary conditions in Eq. (32), we obtain the dispersion relation for the system $\begin{array}{ccc}& & {\mathit{\rho }}_{\mathrm{1}}\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\mathit{g}\mathrm{+}\frac{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\mathrm{-}\frac{\mathit{g}{\mathit{\omega }}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}\mathrm{)}}{{\mathit{m}}_{\mathrm{1}}\mathit{g}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{m}}_{\mathrm{1}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{1}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathrm{=}\\ & & \end{array}$(37)There are a couple of interesting limiting cases to this expression. If we set g = 0, we recover the solution for surface MHD waves in an interface (see, e.g., Wentzel 1979; Roberts 1981) with no instabilities, namely $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{1}}\frac{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\mathrm{=}{\mathit{\rho }}_{\mathrm{2}}\frac{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{2}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{2}}}& & \end{array}$(38)with ${\mathit{m}}_{\mathrm{i}}^{\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{y}}^{\mathrm{2}}\mathrm{+}{\mathit{\omega }}^{\mathrm{4}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{i}}$. Another interesting limit is the incompressible case, which is obtained if we set c_si1 → ∞ and c_si2 → ∞. The dispersion relation is then $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{1}}\frac{\mathit{g}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{1}}}\mathrm{=}{\mathit{\rho }}_{\mathrm{2}}\frac{\mathit{g}{\mathit{m}}_{\mathrm{2}}\mathrm{+}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{2}}^{\mathrm{2}}}{{\mathit{m}}_{\mathrm{2}}}& & \end{array}$(39)with m₁ = k and m₂ =  −k from Eq. (35) in this limit. We can obtain an explicit equation for the frequencies of the modes, $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{2}}\mathrm{=}\mathrm{-}\mathit{gk}\frac{{\mathit{\rho }}_{\mathrm{2}}\mathrm{-}{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{1}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}}\mathrm{+}\frac{\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}{\mathit{c}}_{\mathrm{A}\mathrm{2}}^{\mathrm{2}}\mathrm{)}{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}}{\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}\mathrm{)}}\mathit{,}& & \end{array}$(40)which is equivalent to the classical RTI relation in Eq. (1).

Fig. 2 Linear growth rate of the RTI for a fully ionized plasma (ideal MHD) as a function of the Alfvén speed ${\overline{\mathit{c}}}_{\mathrm{A}}$. In the upper panel curves, different values of the propagation angle θ are shown for a fixed value of β = 0.1, while different β are plotted for a fixed value of θ = 40° in the lower panel curves. In all the panels, the values ρ2/ρ1 = 100, g = 270 m s-2, and k = 10-7 m-1 have been used. The dashed curves correspond to the incompressible MHD limit in Eq. (1).

After checking the limiting cases, we proceed to directly solve Eq. (37). The linear growth rate is plotted in Fig. 2, as compared to the predicted rate from the classical formula in Eq. (1) for the incompressible limit. The main conclusions from these results are as follows:

• 1.

  The threshold is not modified by compressibility. This can beeasily demonstrated by recalling that in ideal MHD the frequencyof the modes is either real or pure imaginary (Goedbloed &Poedts 2004; Goedbloedet al. 2010), so the transition from astable to an unstable situation is necessarily at the points in whichω = 0 is satisfied. Inserting this condition in Eq. (37), we immediately recover $\begin{array}{ccc}\mathit{g}\mathrm{=}\frac{{\mathit{k}}_{\mathit{x}}^{\mathrm{2}}\mathrm{(}{\mathit{\rho }}_{\mathrm{1}}{\mathit{c}}_{\mathrm{A}\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{\rho }}_{\mathrm{2}}{{\mathit{c}}_{\mathrm{A}\mathrm{2}}^{\mathrm{2}}}^{\mathrm{)}}}{\mathrm{(}{\mathit{\rho }}_{\mathrm{2}}\mathrm{-}{\mathit{\rho }}_{\mathrm{1}}\mathrm{)}\mathit{k}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{or}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{=}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta ,}& & \end{array}$(41)which matches with the stability criteria from Eq. (1) and Eq. (16). A magnetic field increase has a stabilizing effect, while increasing the angle between the equilibrium field and the wavevector has the opposite effect, as it happens in the incompressible limit.

• 2.

  The incompressible approximation becomes more valid as θ approaches π/ 2. This is caused when the incompressible limit is recovered as ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\to }\mathrm{\infty }$, which implies Ω → ∞, so terms containing the gravity in Eqs. (35) and (37) are negligible. We can see that increasing the longitudinal wavenumber has a similar effect in the definition of m and the dispersion relation, and hence, the incompressible limit is a better approximation as k is increased.

• 3.

  The linear growth rate for the compressible case is always below the incompressible limit prediction. As β is lowered, the linear growth rate is decreased substantially.

The curves in Fig. 2 tend to zero when the magnetic field is very low. This is caused by the choice of the sound speed: since β is fixed in these curves, ${\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{\to }\mathrm{0}$ also implies that c_s → 0. If the sound speed is prevented from tending to zero as B₀ → 0 (implying that β is no longer held constant and tends to zero too), the drop disappears. There is also a real part of the frequency only when the configuration is stable in this ideal MHD regime, but we focus on the imaginary part and the instability. Leaky modes may also be considered (i.e., modes that propagate in the direction across the surface), but these modes do not appear for the range of parameters selected in these plots.

4.4. Partially ionized plasma

Now, we turn to the general problem with the ambipolar diffusion coefficient different from zero. Equation (29) is a fourth-order ordinary differential equation with constant coefficients, whose solutions are a combination of exponentials e^m_kz, with m_k as one of the four solutions to the indicial equation $\begin{array}{ccc}{\mathit{C}}_{\mathrm{4}}{\mathit{m}}^{\mathrm{4}}\mathrm{+}{\mathit{C}}_{\mathrm{3}}{\mathit{m}}^{\mathrm{3}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}{\mathit{m}}^{\mathrm{2}}\mathrm{+}{\mathit{C}}_{\mathrm{1}}\mathit{m}\mathrm{+}{\mathit{C}}_{\mathrm{0}}\mathrm{=}\mathrm{0.}& & \end{array}$(42)Two of these solutions are close to those in Eq. (35), while the other two are typically larger and depend strongly on the exact value of η_A. Equation (42) must be solved in each of the two regions, and then the two solutions that imply evanescence away from the discontinuity are kept. The general expression of the four roots of this fourth order algebraic equation is massive, so we choose to solve it numerically. Our general solution is then $\begin{array}{ccc}{\mathit{b}}_{\mathit{x}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(43)with the A coefficients being arbitrary constants, the subscript of m denoting the order of the real part among the set of m_k, and the superscript of m as the region where it applies.

The boundary conditions must be applied to obtain a dispersion relation, taking into account that v_z must be obtained by integrating Eq. (28) after inserting the solution for b_x in Eq. (43). Using the same notation in Díaz et al. (2012), the four boundary conditions are written in matricial form as $\begin{array}{ccc}\left[\sum _{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{4}}{\mathit{B}}_{\mathit{ij}}{\mathit{b}}_{\mathit{x}}^{\mathrm{(}\mathrm{3}\mathrm{-}\mathit{j}\mathrm{)}}\right]\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(44)where the index i ∈ [ 1,4] stands for each boundary condition in Eqs. (32) and ${\mathit{b}}_{\mathit{x}}^{\mathrm{\left(}\mathit{j}\mathrm{\right)}}$ is the j-th z-derivative of b_x (j = 0 being the function itself without derivatives and j =  −1 the first integral of the function). The coefficients in this matrix are given in the Appendix. Inserting the solutions from Eq. (43) in Eq. (44), we obtain a system of equations for the A-coefficients, $\begin{array}{ccc}\sum _{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{4}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{h}}_{\mathit{k}}}{\mathit{A}}_{\mathit{k}}{\mathit{B}}_{\mathit{ij}}^{\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right)}}{\mathrm{(}{{\mathit{m}}_{\mathit{k}}^{\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right)}}}^{\mathrm{)}}}^{\mathrm{3}\mathrm{-}\mathit{j}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(45)where the m-coefficients are defined in Eq. (42) with the requirement that the exponentials are bounded at z → ± ∞. In this expression, h_k = 1 for k = 1,2, and h_k = 2 for k = 3,4. The dispersion relation of the system is obtained by requiring the determinant of such system to vanish, namely $\begin{array}{ccc}\left|{\mathit{C}}_{\mathit{ik}}\right|\mathrm{=}\mathrm{0}& & \end{array}$(46)with the C-matrix defined as $\begin{array}{ccc}{\mathit{C}}_{\mathit{ik}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{h}}_{\mathit{k}}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}\mathit{,}\mathrm{4}}{\mathit{B}}_{\mathit{ij}}^{\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right)}}{\mathrm{(}{{\mathit{m}}_{\mathit{k}}^{\mathrm{\left(}{\mathit{h}}_{\mathit{k}}\mathrm{\right)}}}^{\mathrm{)}}}^{\mathrm{3}\mathrm{-}\mathit{j}}\mathit{.}& & \end{array}$(47)The imaginary part of the solutions to Eq. (46) are plotted in Fig. 3 with the growth rates predicted by the overplotted incompressible and compressible RTI.

Fig. 3 Linear growth rate of the RTI for a partiallyl ionized plasma (taking into account ambipolar diffusion) as a function of the Alfvén speed ${\overline{\mathit{c}}}_{\mathrm{A}}$. The dashed line corresponds to the incompressible limit given in Eq. (1) and the black solid line to the compressible MHD results from Sect. 4.3. The values ρ2/ρ1 = 100, β = 0.1, θ = 40, k = 10-7 m-1, ccrit = 33 km s-1, and ξn1 = 10-4 have been used.

The following points must be emphasized:

• 1.

  A similar plot to Fig. 3 with higher values of θ would draw the collisional plasma results closer to the incompressible MHD limit and modify the critical speed. Thus, the incompressible limit is a much better approximation as θ increases, as has happened in the collisionless plasma.

• 2.

  The instability threshold is no longer the one predicted by Eq. (1), namely ${\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{=}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }\mathrm{=}\mathrm{47.4}$ km s^-1 for the parameters used in the plot. The configuration is unstable for all the values of ξ_n and magnetic field. Considering the ambipolar term as the only effect of the presence of neutrals is enough to render this configuration (a heavier partially ionized fluid above a lighter fluid) unstable no matter how strong the magnetic field is.

  Fig. 4 Linear growth rate of the RTI for a partiallyl ionized plasma as a function of perturbation wavenumber k. The dashed line corresponds to the incompressible limit given in Eq. (1). The values ρ2/ρ1 = 100, θ = 40, ${\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{=}\mathrm{30}$, and ξn1 = 10-6 have been used. Solid lines are calculated with β = 0.1 and dot-dashed lines with β = 0.5, while red lines have a value for the neutral fraction ξn2 = 0.9, blue lines ξn2 = 0.5, and purple lines ξn2 = 0.1.

• 3.

  For parameters which are classically unstable (${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{<}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$), the linear growth rate is reduced a lot with respect to Eq. (1) because of the compressibility. Increasing the ambipolar coefficient slightly raises the growth rate, but the effect of the ambipolar term is small compared with compressibility in this range.

• 4.

  For parameters which are classically stable (${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{>}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$), the ambipolar diffusion still drives the instability, but the linear growth rate in this regime is an order of magnitude smaller than in the classically unstable range.

• 5.

  Close to the stability threshold (${\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{\approx }{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$), the differences induced by the ambipolar diffusion term are relatively higher.

• 6.

  In any case, the growth rate in all the parameter space is significantly lower than that of an uncoupled neutral gas that is subject to the hydrodynamic RTI (Eq. (1) with B₀ = 0, ρ_n1 and ρ_n2), which would be Im[ ω] = 0.0053 s^-1 for the parameters in the plot. The collisional coupling between neutrals and charged particles prevents the neutrals from fully developing their instability, even for values of ξ_n2 close to 1.

• 7.

  As the ionization fraction is raised (and thus ξ_n2 and ξ_n1 are lowered), the curve resembles the MHD limit: there is a bifurcation very close to the critical value which can not be clearly seen in the scale of these plots and its value tends to c_crit/cos θ as the neutral fraction tends to zero.

With the inclusion of the ambipolar term, the frequencies of the modes are no longer restricted to be either pure real or pure imaginary as in the ideal MHD limit (collisionless plasma). The solutions plotted in Fig. 3 have a real conterpart Re[ ω], which are not shown in the plot. This real part is close to the compressible MHD results when ${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{>}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$ and much smaller than Im[ ω] when ${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{<}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$. We can check that there is no critical value of ${\overline{\mathit{c}}}_{\mathrm{A}}$ for which the system becomes stable: if we require ω → 0, the only real solution is c_crit = 0, confirming the numerical results in Fig. 3 and the absence of a stable region in the parameter space.

Another important parameter is the perturbation wavenumber k. So far, we have fixed a value of k = 10^-7 m^-1, following the typical wavenumbers from the fast transversal MHD modes of a prominence thread that was used in previous studies in prominence seismology and RTI instability in threads (Díaz et al. 2002, 2012; Terradas et al. 2012). We can further explore the effects of the initial perturbation. One important consequence is obtained after scaling the problem: it can be shown that the curves can be rescaled using the variables ω/ (ck) and ${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{/}\mathit{c}$ in the ideal MHD limit (with c a characteristic speed, such as Alfvén speed in the prominence, for example). This scaling involves the adimensional quantity η_Ak/c if the ambipolar diffusion is included. Hence, increasing the wavenumber perturbation has the direct effect of increasing the relevance of the ambipolar term. This is expected, since it is known that the ambipolar diffusion grows as the typical length scale is reduced. On the other hand, we can also plot the linear growth rate as a function of the perturbation wavenumber (Fig. 4). We see the same main effects: there is no stable regime; the compressibility lowers the growth rate for ${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{<}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$, and the ambipolar diffusion slightly raises it as ξ_n2 is increased. It is also interesting to study this dependence near the incompressible limit with values of θ close to 90° (Fig. 5). Since compressibility is no longer dominant, the inclusion of the ambipolar term raises the growth rate over the classical RTI (as reported in Shadmehri et al. 2013 for the same assumptions that Díaz et al. 2012 uses in the context of local bubble of the solar system).

Fig. 5 Same plot as in Fig. 4 for the values of θ = 89 and ${\overline{\mathit{c}}}_{\mathrm{A}}\mathrm{=}\mathrm{100}$ (near the incompressible limit). Here, the values of the beta are β = 0.01 in the solid lines and β = 0.1 in the dashed lines.

Finally, we can plot the growth rate vs. the ambipolar diffusivity (Fig. 6). As mentioned above, the rate is slightly modified if the configuration is unstable in the ideal MHD limit (${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{<}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$, upper panel) unless a very unrealistical high value of the ambipolar diffusivity is assumed. For configurations that are close to the critical value (middle panel), the dependence on the ambipolar diffusivity is more important. In the stable range (${\overline{\mathit{c}}}_{\mathrm{A}}\mathit{>}{\mathit{c}}_{\mathrm{crit}}\mathit{/}\cos \hspace{0.17em}\mathit{\theta }$, lower panel), the linear growth rate is never zero (so strictly speaking the configuration is unstable), but the linear growth rate is at least about an order of magnitude lower than the values in classically unstable regime for typical values of η_A, so the instability would take an excessively long time to develop in practice.

Fig. 6 Linear growth rate of the RTI for a partially ionized plasma as a function of the ambipolar diffusion coefficient ηA2 with some values of the neutral fraction ξn2 for which the values of ηA2 are obtained and also included in the plot as large points. The values ρ2/ρ1 = 100, β = 0.1, k = 10-7 m-1, θ = 40°, and ξn1 = 10-6 have been used (so ccrit = 47 km s-1). The upper panel corresponds to a classically unstable configuration, the middle panel to a marginally stable configuration, and the lower panel to a classically stable configuration. The dotted line in the upper panel corresponds to the MHD limit (in the middle and lower panel the classical limit is zero).

5. Full-induction equation

We have dealt in the previous section with the effect of the induction term and the ambipolar diffusion term alone in the generalized induction equation for partially ionized plasmas (Eq. (15)). Considering only these terms has allowed us to analytically solve the linearized equations, but the effect of the other supposedly smaller terms must be taken into account. In this case, obtaining analytical solutions can be much more demanding, even in the linearized problem.

5.1. Linearized equations

In order to solve this problem, we first need to eliminate the partial pressures and the density of each species in the induction equation. This can be done by using the definition of the partial pressure p_α = n_αk_BT_α and invoking that the temperature of each species is the same (T_i = T_e = T_n). Hence, $\begin{array}{ccc}\mathit{p}\mathrm{=}{\mathit{p}}_{\mathrm{i}}\mathrm{+}{\mathit{p}}_{\mathrm{e}}\mathrm{+}{\mathit{p}}_{\mathrm{n}}\mathrm{=}\mathrm{2}{\mathit{p}}_{\mathrm{e}}\mathrm{+}{\mathit{p}}_{\mathrm{n}}\mathrm{=}{\mathit{p}}_{\mathrm{e}}\mathrm{(}\mathrm{2}\mathrm{+}{\mathit{\xi }}_{\mathrm{n}}\mathit{/}{\mathit{\xi }}_{\mathrm{i}}\mathrm{)}\mathit{,}& & \end{array}$(48)and we obtain $\begin{array}{ccc}{\mathit{p}}_{\mathrm{e}}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}{\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}\hspace{0.17em}\mathit{p}\mathrm{;}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{p}}_{\mathrm{n}}\mathrm{=}\frac{{\mathit{\xi }}_{\mathrm{n}}}{\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}\hspace{0.17em}\mathit{p}\mathrm{;}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{e}}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{e}}}{{\mathit{m}}_{\mathrm{i}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}\mathrm{)}\mathit{\rho ,}& & \end{array}$(49)so we can operate Eq. (9) to obtain an expression for the pressure gradients in the battery term. That is, $\begin{array}{ccc}{G}& \mathrm{=}& \frac{{\mathit{\xi }}_{\mathrm{n}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}\mathrm{)}}{\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}\mathrm{\nabla }\mathit{p,}\\ \frac{\mathit{\epsilon }{G}\mathrm{-}\mathrm{\nabla }{\mathit{p}}_{\mathrm{e}}}{\mathit{e}{\mathit{n}}_{\mathrm{e}}}& \mathrm{=}& \frac{{\mathit{m}}_{\mathrm{e}}}{\mathit{e}{\mathit{\rho }}_{\mathrm{e}}}\frac{\mathrm{1}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}{\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}\mathrm{(}\mathit{\epsilon }{\mathit{\xi }}_{\mathrm{n}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{\nabla }\mathit{p}\mathrm{=}\frac{{\mathit{m}}_{\mathrm{i}}}{\mathit{e}}\hspace{0.17em}\frac{\mathit{\epsilon }{\mathit{\xi }}_{\mathrm{n}}\mathrm{-}\mathrm{1}}{\mathrm{2}\mathrm{-}{\mathit{\xi }}_{\mathrm{n}}}\hspace{0.17em}\frac{\mathrm{\nabla }\mathit{p}}{\mathit{\rho }}\mathrm{·}\end{array}$Here it is explicit that the G combination vanishes for a fully ionized plasma (ξ_n = 0) and the only contribution to the battery term in that case is the well-known form of the electron pressure gradient.