© ESO 2020

1. Introduction

The hydrodynamic Rayleigh–Taylor instability (RTI) was first explored by Rayleigh (1900) and Taylor (1950), for the incompressible case. This concerns the instability of an interface that separates two fluids with different properties, where a dense fluid is supported above a lighter fluid, in the presence of a transversal gravitational force (or equivalently, an acceleration of the fluid system in the direction of the denser fluid). The properties and consequences of this instability are well established. Whether a similar instability might occur in magnetic fluids (plasmas) was studied initially by Kruskal & Schwartzschild (1954), for example. Their investigation showed that the magnetic field provides a stabilising effect to perturbations in the direction of the magnetic field by supplying magnetic tension, which counteracts some of the destabilising effects of buoyancy. Clearly, this statement remains true as long as the equilibrium magnetic field has a horizontal component. Perturbations of the interface in a magnetic field that is perpendicular to the interface, however, behave similarly to the hydrodynamic case, but with an altered magnitude. The particular effect of the inclination of the wave propagation with respect to the magnetic field direction was first investigated by Chandrasekhar (1961). Perturbations of the horizontal equilibrium surface (or interface) give rise to oscillatory surface motion. The displaced fluid (essentially a fluid interchange) experiences a potential energy change that destabilises the interface. Numerical investigations of RTI show the evolution of an ordered system into a stage where the two fluids mix and small-scale structures (or filaments) are developed.

The problem of RTI in magnetic fluids (also called magnetic RTI, or MRT instability) received special attention in several astrophysical applications. For example, the filamentary structure of the Crab nebula could well be due to this instability, as was corroborated by comparisons made by Hester et al. (1996) between the observed structuring and MRT instability simulations made by Jun et al. (1995). Other examples where MRT instabilities may develop are in supernovae (e.g., Fryxwell et al. 1991; Jun & Norman 1996; Chevalier 1982), accretion discs (e.g., Wang & Nepveu 1983; Kulkarni & Romanova 2008), buoyant magnetic bubbles in galaxies (e.g., Robinson et al. 2004; Jones & De Young 2005) for (2D) simulations, and O’Neill et al. (2009) for (3D) and relativistic jets, see also Maatsumoto & Masad (2013).

In parallel with the increase in high-resolution observational capabilities, MRT also received special attention in solar and solar-terrestrial physics. Isobe et al. (2005, 2006) suggested that MRT instabilities are responsible for the filamentary structure in emerging flux regions. Most notably, solar prominences are likely to become Rayleigh–Taylor unstable because they are composed of cool and dense plasma that is suspended above much lighter coronal plasma. The formation of filamentary threads of prominences was modelled analytically by Terradas et al. (2012), for example, who found that the RTI caused filaments in quiescent prominences to have very short lifetimes, but that the magnetic tension in active regions may be sufficient to stabilise these prominences. MRT instabilities in prominences have been invoked by Ryutova et al. (2010) to explain the formation of plumes and even to determine the magnetic field strength from the wavelength and growth rate of instabilities. This was expanded upon by Innes et al. (2012), who used the critical wavelength at which waves become unstable to estimate the Alfvén speed. Carlyle et al. (2014) used the most unstable mode analysis to find the magnetic field strength in fragmenting eruptions of filaments. More recently, several studies have been devoted to the investigation of the combined effects of Kelvin–Helmholtz instability upon so-called Rayleigh–Taylor fingers. This has in particular been applied to observations of plumes or bubbles within prominences (Berger et al. 2010, 2017; Mishra et al. 2018; Mishra & Srivastava 2019). For a more comprehensive review of RTI within prominences, see Hillier (2017).

Over the years, the theory of MRT instabilities has been advanced in several aspects. Compressibility effects have been widely studied in different context by e.g., Vandervoort (1961), Shivamoggi (1982), Bernstein & Book (1983), Ribeyre et al. (2004), Livescu (2004), and these studies showed unequivocally that compressibility is able to stabilise the interface against the MRT instability. More recently, several studies focussed on the effect of a non-parallel magnetic field, that is, magnetic shear at a tangential discontinuity e.g., Ruderman et al. (2014) and Ruderman (2017). The growth time derived by these authors was found to be dependent on the shear angle, and the instability increment obtained by these authors allowed them to determine that the optimal shear angle of the magnetic field is 13 degrees, based on Hinode observations (see Okamoto et al. 2007). Because the temperatures of solar prominences are often not high enough to completely ionise the plasma, the effect of partial ionisation on the generation and evolution of MRT instabilities was also studied by e.g., (Diaz et al. 2014; Ruderman et al. 2017). These investigations showed that the MRT instability only becomes sensitive to the degree of plasma ionisation for plasmas with low values of plasma beta and in a very weakly ionised state. Moreover, perturbations are only unstable for wavenumbers that are below a cut-off value.

While these investigations into RTI are far reaching, almost all theoretical investigations are concerned with a tangential discontinuity, where the magnetic field is parallel to the interface. However, the magnetic field in the solar atmosphere is in many cases not strictly parallel to separating interfaces (the penumbras of sunspots, solar prominences, plumes or interplumes in coronal holes, etc.), meaning that a new approach is needed to analyse how instabilities are generated at these interfaces. The theory of discontinuities includes tangential discontinuities in which the magnetic field is parallel to the separating interface. In this case, the normal velocity and total pressure (kinetic and magnetic) components are conserved quantities across the interface. Density, kinetic pressure, and tangential component of the magnetic field can be discontinuous across the interface. In contrast, when the magnetic field permeates the interface, the discontinuities are contact discontinuities, where the kinetic pressure, the magnetic field, and the velocity are all continuous quantities and only mass density and temperature are allowed to have a jump. In an incompressible plasma, the only allowable discontinuous quantity is the mass density.

The aim of the present work is to expand the MRT instability theory to the case of contact discontinuities. This has applications to arcade-type prominences with their typical “dipped” magnetic field structure that intersects the denser upper plasma, and so contact discontinuities become pertinent. We analytically calculate the growth rate in terms of equilibrium parameters for the case of a contact discontinuity of an incompressible plasma.

The structure of this paper is as follows. Section 2 formulates the problem and sets out the equilibrium configuration we explore Sect. 3 is dedicated to the analytical derivation of the dispersion relation by solving the governing equations on both sides of the interface and then connecting the solutions using specific boundary (jump) conditions. Our analytical results are further investigated numerically in Sect. 4 for a number of parameters, and we determine the variation of the instability growth rate with respect to the wavenumber, inclination angle of the magnetic field, propagation angle, density ratio, and magnetic field strength. In Sect. 5 we apply our theoretical results to diagnose the inclination angle of the magnetic field based on high-resolution observations of prominence threads. Finally, we conclude and summarise our results in Sect. 6.

2. Problem formulation

The properties of the MRT instability were investigated by assuming an interface in the (x, y)-plane situated at z = 0 that separates two distinct plasma regions. The density changes sharply at this interface according to

$$\begin{array}{c}\hfill {\rho }_0(z)=\{\begin{array}{cc}{\rho }_1(z),\hfill & z<0,\hfill \\ {\rho }_2(z),\hfill & z>0.\hfill \end{array}\end{array}$$(1)

The contact discontinuity was attained by considering an equilibrium magnetic field that intersects the density interface at an angle, θ, therefore we took the equilibrium magnetic field to be B₀ = B₀(cos θ, 0, sin θ). Because the magnetic field strength is required to be continuous across the interface, the background magnetic field has the same value in both plasma regions and B₀ is constant. The equilibrium configuration is shown in Fig. 1.

Fig. 1. Sketch of the equilibrium configuration. The equilibrium state consists of a surface separating two regions, each with different density. The magnetic field is uniform throughout both regions and is inclined within the (x, z)-plane at an angle, θ, with respect to the x-direction. The configuration is invariant in the x- and y-directions. Perturbations are described by the wave vector, k, in the (x, y) plane.

The dynamics of the plasma is described using the incompressible, linearised, and ideal set of magnetohydrodynamics (MHD) equations,

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial \mathit{v}}{\partial t}=-\mathrm{\nabla }p+(\mathrm{\nabla }\times \mathit{b})\times \frac{{\mathit{B}}_0}{{\mu }_0}-\rho g\widehat{\mathit{z}},\end{array}$$(2)

$$\begin{array}{c}\hfill \frac{\partial b}{\partial t}=({\mathit{B}}_0·\mathrm{\nabla })\mathit{v},\end{array}$$(3)

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathit{v}=0,\mathrm{\nabla }·\mathit{b}=0,\end{array}$$(4)

where all quantities with subscript 0 denote the equilibrium state, b = (b_x, b_y, b_z) is the perturbation of the magnetic field, v = (v_x, v_y, v_z) is the velocity perturbation, p and ρ are the perturbations of kinetic pressure and density, μ₀ is the permeability of free space, g = 274 m s⁻² is the constant gravitational acceleration, and $\widehat{\mathit{z}}$ is the unit vector along the z-axis. Equation (4) describes the incompressible limit and the solenoidal condition, respectively.

We considered wave-like perturbations directed along the interface situated at z = h(x, y, t) and wrote all perturbed variables as f ∼ f̂ exp[i(k_xx+k_yy−ωt)]. Because gravity renders the plasma inhomogeneous, we took the background pressure and density to depend on height alone, that is, ρ₀ = ρ₀(z), p₀ = p₀(z). Because the equilibrium magnetic field crosses the density interface, the surface h(x, y, t) becomes a contact discontinuity. In order to connect the solutions of both sides across the discontinuity, we imposed boundary conditions that were different from those used in the case of tangential discontinuity. In our case, first of all, we required the continuity across the interface of all three components of the velocity and the magnetic field. In the linear approximation, they are given by

$$\begin{array}{c}\hfill \left[\left[\mathit{v}\right]\right]=0,\left[\left[\mathit{b}\right]\right]=0,\end{array}$$(5)

where ⟦⋅⟧ denotes the jump of a quantity across the discontinuity. To determine the final continuity condition, we recall that the equilibrium pressure balance is valid, that is,

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}_0}{\mathrm{d}z}=-g{\rho }_0.\end{array}$$

Assuming small changes in pressure (linear approximation), we can expand the pressure perturbation at the interface in Taylor series as

$$\begin{array}{c}\hfill p(z=h)=p(0)+h\frac{\mathrm{d}{p}_0}{\mathrm{d}z}{|}_{z=0}=p(0)-hg{\rho }_0(0).\end{array}$$(6)

By differentiating with respect to time the condition for continuity of pressure at z = h(x, y, t), we therefore arrive at

$$\begin{array}{c}\hfill \left[\left[g{\rho }_0{v}_z-\frac{\partial p}{\partial t}\right]\right]=\left[\left[i\omega p+g{\rho }_0{v}_z\right]\right]=0.\end{array}$$(7)

In what follows, we solve the system of MHD equations on both sides of the interface and connect the solutions at the interface using the boundary (jump) conditions. The resulting equation constitutes the dispersion relation that we investigate below.

3. Dispersion relation of waves propagating at the density interface

When all perturbations proportional to the exponential ansatz introduced in Sect. 2 are considered, Eqs. (2)–(4) can be written as

$$\begin{array}{c}\hfill -i\omega {\rho }_0{v}_x=-i{k}_{x}p+\frac{B_0}{\mu }sin\theta (\frac{\partial b_x}{\partial z}-i{k}_x{b}_z),\end{array}$$(8)

$$\begin{array}{cc}\hfill -i\omega {\rho }_0{v}_y& =-i{k}_{y}p+\frac{B_0}{\mu }sin\theta (\frac{\partial b_y}{\partial z}-i{k}_y{b}_z)\hfill \\ \hfill & +\frac{B_0}{\mu }cos\theta (i{k}_x{b}_y-i{k}_y{b}_x),\hfill \end{array}$$(9)

$$\begin{array}{c}\hfill -i\omega {\rho }_0{v}_z=-\frac{\partial p}{\partial z}+\frac{B_0}{\mu }cos\theta (i{k}_x{b}_z-\frac{\partial b_x}{\partial z})-\rho g,\end{array}$$(10)

$$\begin{array}{c}\hfill -i\omega \mathit{b}=B_0(i{k}_{x}cos\theta \mathit{v}+sin\theta \frac{\partial \mathit{v}}{\partial z}).\end{array}$$(11)

These equations can be combined into two relations that connect v_x and v_z as

$$\begin{array}{c}\hfill i\frac{k^2}{k_x}f(v_x)=-f\left(\frac{\partial v_z}{\partial z}\right)\end{array}$$

$$\begin{array}{c}\hfill \frac{1}{k_x}[f\left(\frac{\partial v_x}{\partial z}\right)+{\omega }^2\frac{{\rho }_0\prime }{{\rho }_0}{v}_x]=i[f(v_z)+g\frac{{\rho }_0\prime }{{\rho }_0}{v}_z],\end{array}$$(12)

where the operator function f is defined as

$$\begin{array}{c}\hfill f=v_{\mathrm{A}}^2{sin}^2\theta \frac{\partial {}^2}{\partial z^2}+i{k}_x{v}_{\mathrm{A}}^{2}sin2\theta \frac{\partial }{\partial z}+({\omega }^2-k_x^2{v}_{\mathrm{A}}^2{cos}^2\theta ).\end{array}$$(13)

Using Eq. (8), we can now express the kinetic pressure perturbation as

$$\begin{array}{cc}\hfill p=& [\frac{1}{k_x}\frac{{\rho }_0{v}_{\mathrm{A}}^2}{\omega }{sin}^2\theta \frac{\partial {}^2}{\partial z^2}+i\frac{{\rho }_0{v}_{\mathrm{A}}^2}{2\omega }sin2\theta \frac{\partial }{\partial z}+\frac{\omega }{k_x}{\rho }_0]v_x\hfill \\ \hfill & +{\rho }_0{v}_{\mathrm{A}}^2[-\frac{i}{\omega }{sin}^2\theta \frac{\partial }{\partial z}+\frac{k_x}{2\omega }sin2\theta ]v_z.\hfill \end{array}$$(14)

When the vertical perturbation scales are much smaller than the gravitational scale-height, such as at the onset of instability in the given problem, we can take the long scale-height limit and so treat the plasma either side of the interface as uniform. Given this approximation, we simplified our treatment by considering ${\rho }_0^{\prime }=0$, therefore we combined the equations given by Eq. (12) into the governing equation of the wave propagation in the incompressible plasma,

$$\begin{array}{c}\hfill (k^2-\frac{\partial {}^2}{\partial z^2})f(v_z)=0.\end{array}$$(15)

We assumed that the solution of this equation is of the form v̂_z ∼ e^Γz, where the quantity Γ is a complex quantity. Equation (15) then simplifies to

$$\begin{array}{c}\hfill ({\mathrm{\Gamma }}^2-k^2)(v_{\mathrm{A}}^2{sin}^2\theta {\mathrm{\Gamma }}^2+i{k}_x{v}_{\mathrm{A}}^{2}sin2\theta \mathrm{\Gamma }+{\omega }^2-k_x^2{v}_{\mathrm{A}}^2{cos}^2\theta )=0,\end{array}$$(16)

which admits the roots

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\pm k,m_{\pm },\end{array}$$

where the expression of m_± is given by

$$\begin{array}{c}\hfill m_{\pm }=\frac{-i[v_{\mathrm{A}}{k}_{x}cos\theta \pm \omega ]}{v_{\mathrm{A}}sin\theta }·\end{array}$$(17)

The form of the governing equation derived in Eq. (12) implies that the wave frequency is complex, with the imaginary part of ω describing the temporal evolution of the amplitudes of perturbations. According to the temporal variation of perturbations assumed earlier, it is clear that a positive imaginary part of ω describes an unstable amplification of the amplitude. Because we are interested in waves that are localised at the surface, we assumed that far away from the interface, waves are evanescent. The z-component of the velocity is therefore of the form

$$\begin{array}{c}\hfill {\widehat{v}}_z=\{\begin{array}{c}{A}_1{e}^{\mathit{kz}}+B_1{e}^{m_{+1}z},z<0,\hfill \\ A_2{e}^{-kz}+B_2{e}^{m_{-2}z},z>0,\hfill \end{array}\end{array}$$(18)

where A₁, A₂, B₁, and B₂ are arbitrary constants. Below we are only interested in unstable perturbations. Hence we assume that ℑ(ω) > 0. This also ensures that waves are evanescent (through the roots m_±). While the first terms in Eq. (18) represent an exponential decay with z, the second terms describe an oscillatory decay. It is also obvious that when θ = 0, the second terms of Eq. (18) vanish, and we recover the result of a tangential discontinuity.

Our analysis here, where the two plasma regions are treated as homogeneous, is valid when 1/|ℜ(Γ)| is much less than the gravitational scale-height, H, that is,

$$\begin{array}{c}\hfill H\gg \frac{1}{k},\frac{v_{\mathrm{A}1}sin\theta }{{\omega }_i},\frac{v_{\mathrm{A}2}sin\theta }{{\omega }_i}·\end{array}$$

In other words, stratification effects may be ignored when the perturbation wavelength and the product of the instability time with the Alfvén speeds are sufficiently small.

With the help of Eq. (18), we determined the corresponding equations for the pressure and magnetic field perturbations. The solutions obtained for the two regions are connected at the interface using the boundary conditions specified earlier.

3.1. Boundary conditions

Using the surface wave solutions of Eq. (18), we employed the boundary conditions given by Eqs. (5) and (7) to find the dispersion relation of incompressible waves propagating along the interface. We may now express these conditions (including Eq. (7)) in terms of v̂_x and v̂_z.

The incompressibility condition yields

$$\begin{array}{c}\hfill \frac{\partial v_z}{\partial z}=-i(k_x{v}_x+k_y{v}_y).\end{array}$$

Taking into account the behaviour of v_x and v_y at the interface, we can obtain in addition that the condition

$$\begin{array}{c}\hfill \left[\left[\frac{\partial v_z}{\partial z}\right]\right]=0\end{array}$$(19)

has to be also satisfied. Imposing the continuity of b_x across the interface together with Eq. (11) results in

$$\begin{array}{c}\hfill \left[\left[isin\theta \frac{\partial v_x}{\partial z}-k_{x}cos\theta v_x\right]\right]=0.\end{array}$$

Combining this condition with continuity of v_x, we obtain

$$\begin{array}{c}\hfill \left[\left[\frac{\partial v_x}{\partial z}\right]\right]=0.\end{array}$$(20)

Furthermore, based on the same consideration, the boundary condition written for b_y simplifies to

$$\begin{array}{c}\hfill \left[\left[\frac{\partial {}^2{v}_z}{\partial z^2}\right]\right]=0.\end{array}$$(21)

Finally, we can use the expression we derived for the kinetic pressure (Eq. (14)) together with the jump conditions derived above to transform the dynamic boundary condition Eq. (7) as

$$\begin{array}{c}\hfill \left[\left[i{\rho }_0(v_{\mathrm{A}}^2{sin}^2\theta \frac{\partial {}^2{v}_x}{\partial z^2}+{\omega }^2{v}_x)+k_{x}g{\rho }_0{v}_z\right]\right]=0.\end{array}$$(22)

This boundary conditions was used to connect the solutions from both sides of the interface in order to derive the dispersion relation of waves.

3.2. Derivation of the dispersion relation

In order to use the boundary conditions specified above, we first determined the expression for v_x. Using the expression of v_z given by Eq. (18), we solved the differential equation for v_x (see the first equation of system (12)), and we obtain that

$$\begin{array}{c}\hfill {\widehat{v}}_x=\{\begin{array}{cc}i\frac{k_x}{k}{A}_1{e}^{\mathit{kz}}+C_1{e}^{m_{+1}z}\hfill & z<0,\hfill \\ -i\frac{k_x}{k}{A}_2{e}^{-kz}+C_2{e}^{m_{-2}z}\hfill & z>0,\hfill \end{array}\end{array}$$(23)

where C₁ and C₂ are more arbitrary constants.

Using the exponential forms of v_x and v_z specified by Eqs. (18) and (23), we then applied the boundary conditions (Eqs. (19)– (22) along with ⟦v_x⟧ = 0 and ⟦v_z⟧ = 0). After long but straightforward calculations, we obtain a system of linear equations for the unknown six constants that in matrix form can be given as

$$\begin{array}{c}\hfill \left[\begin{array}{cccccc}1& 1& 0& -1& -1& 0\\ i\frac{k_x}{k}& 0& 1& i\frac{k_x}{k}& 0& -1\\ k& m_{1+}& 0& k& -m_{2-}& 0\\ k^2& m_{1+}^2& 0& -k^2& -m_{2-}^2& 0\\ i{k}_x& 0& m_{1+}& -i{k}_x& 0& -m_{2-}\\ {\alpha }_1& {\beta }_1& {\gamma }_1& -{\alpha }_2& -{\beta }_2& -{\gamma }_2\end{array}\right]\left[\begin{array}{c}{A}_1\\ B_1\\ C_1\\ A_2\\ B_2\\ C_2\end{array}\right]=M\left[\begin{array}{c}{A}_1\\ B_1\\ C_1\\ A_2\\ B_2\\ C_2\end{array}\right]=0,\end{array}$$(24)

where the various entries of the matrix M are defined as

$$\begin{array}{c}\hfill {\alpha }_1=k_{x}g{\rho }_1-{\rho }_1\frac{k_x}{k}(v_{\mathrm{A}1}^2{k}^2{sin}^2\theta +{\omega }^2),\end{array}$$

$$\begin{array}{c}\hfill {\alpha }_2=k_{x}g{\rho }_2+{\rho }_2\frac{k_x}{k}(v_{\mathrm{A}2}^2{k}^2{sin}^2\theta +{\omega }^2),\end{array}$$

$$\begin{array}{c}\hfill {\beta }_1=k_{x}g{\rho }_1,{\beta }_2=k_{x}g{\rho }_2,\end{array}$$

$$\begin{array}{c}\hfill {\gamma }_1=i{\rho }_1(v_{\mathrm{A}1}^2{m}_{1+}^2{sin}^2\theta +{\omega }^2),{\gamma }_2=i{\rho }_2(v_{\mathrm{A}2}^2{m}_{2-}^2{sin}^2\theta +{\omega }^2),\end{array}$$

and $k={(}^{k_x^2}$. The non-trivial solution to the homogeneous system of Eq. (24) only exists when det(M)=0, which provides the dispersion relation of wave perturbations propagating along the interface in the presence of gravity and inclined magnetic field as

$$\begin{array}{c}\hfill (m_{2-}-m_{1+})(k-m_{1+})(k+m_{2-})S=0,\end{array}$$(25)

where

$$\begin{array}{cc}\hfill S& =(d^{1/2}+1)(d-1){\omega }^3-2ik{d}^{1/2}{v}_{\mathrm{A}1}sin\theta {(d^{1/2}+1)}^2{\omega }^2\hfill \\ \hfill & +(d^{1/2}+1)[2d{v}_{\mathrm{A}1}^2(k^2-k_y^2{cos}^2\theta )+gk(d-1)]\omega \hfill \\ \hfill & +2i{k}^2{d}^{1/2}{v}_{\mathrm{A}1}gsin\theta (d-1),\hfill \end{array}$$(26)

and d = ρ₁/ρ₂ <  1 is the density ratio. The dispersion relation Eq. (25) has five non-trivial solutions (the first multiplier in the dispersion relation corresponds to a trivial solution) in terms of ω. The second and third brackets admit the two solutions

$$\begin{array}{c}\hfill {\omega }_1=k_x{v}_{\mathrm{A}1}cos\theta -ik{v}_{\mathrm{A}1}sin\theta \end{array}$$

$$\begin{array}{c}\hfill {\omega }_2=k_x\sqrt{d}{v}_{\mathrm{A}1}cos\theta -ik\sqrt{d}{v}_{\mathrm{A}1}sin\theta .\end{array}$$(27)

In order to simplify the discussion, we introduce the propagation angle within the (x, y)-plane, α, such that

$$\begin{array}{c}\hfill k_x=kcos\alpha ,k_y=ksin\alpha .\end{array}$$

In addition, we use one single Alfvén speed (v_A1), and write the Alfvén speed in region 2 as ${\mathit{v}}_{\mathrm{A}2}=\sqrt{d}{\mathit{v}}_{\mathrm{A}1}$.

Clearly, ω₁ and ω₂ have negative imaginary parts, and these solutions were disregarded because we only considered ℑ(ω)> 0. Hence, we considered the three solutions to S = 0. In order to gain more information about the unstable solutions, this may be reformatted in terms of Ω = −iω, such that ℜ(Ω) = ℑ(ω), and the unstable solutions are given by solutions where the real part of Ω is positive. We find that S = −iσ, where

$$\begin{array}{cc}\hfill \sigma =& (d^{1/2}+1)(1-d){\mathrm{\Omega }}^3+2k{d}^{1/2}{v}_{\mathrm{A}1}sin\theta {(d^{1/2}+1)}^2{\mathrm{\Omega }}^2\hfill \\ \hfill & +(d^{1/2}+1)[2d{v}_{\mathrm{A}1}^2(k^2-k_y^2{cos}^2\theta )-gk(1-d)]\mathrm{\Omega }\hfill \\ \hfill & -2k^2{d}^{1/2}{v}_{\mathrm{A}1}gsin\theta (1-d).\hfill \end{array}$$(28)

Because the coefficients of Ω³ and Ω² are both positive (when d <  1), the sum of the roots must be negative, and because the coefficient of Ω⁰ is negative, the product of the roots must be positive. Assuming that at least one solution is unstable, that is, positive Ω, these two conditions lead to the fact that the other two roots for Ω must have negative real parts and so are unphysical amplifying modes. This leaves one physical solution that we display graphically to explore the effects of varying wavenumber, k, inclination angle, θ, density ratio, d, and propagation direction, α, shown in Figs. 2–5.

Fig. 2. Solutions for the dispersion relation for waves propagating in the (x, y)-plane for two propagation directions (α = 0 solid line, and α = π/4 dashed line) with respect to the wavenumber, k. The density ratio is d = 0.5, the Alfvén speed of the lower plasma is vA = 10 km s−1, and the magnetic field inclination angle is θ = 5°. The upper panel shows the real part of the frequency, while the imaginary part is plotted in the lower panel. The limiting value of ℑ(ω) when k → ∞ is shown by the grey horizontal lines in the lower panel. For illustration we also show the real and imaginary part of the frequency in the case of a tangential discontinuity (θ = 0), plotted here in green.

Fig. 3. Imaginary part of the frequency for waves propagating in the (x, y) plane for several propagation directions with respect to the magnetic field inclination, θ. The density ratio is assumed to be d = 0.5, the reference Alfvén speed is vA1 = 10 km s−1, and the value of the wavenumber is fixed at k = 10−5 m−1.

Fig. 4. Growth time for k = 10−5 m−1, d = 0.5 vA = 10 km s−1 in terms of propagation angle, α, and field inclination, θ.

Fig. 5. Growth time for k = 10−5 m−1, d = 0.5 θ = 5° in terms of propagation angle, α, and Alfvén speed of the lower plasma, vA.

4. Solutions and results

The variation of the only physically acceptable root of S(ω) = 0 with various physical parameters is investigated numerically. In Fig. 2 we plot the variation of the real and imaginary part of the frequency with respect to the wavenumber, k, for a small field inclination angle. In order to compare our results with the well-known results obtained in the case of a tangential discontinuity, we plot the results we obtain for θ = 0 (green lines). We chose two values for the propagation direction: α = 0 (propagation along the x-axis, solid line) and α = π/4 (dashed line). The smallest instability increment (longest amplification time) is obtained for propagation parallel to the x-axis, while when the direction of propagation is increased away from α = 0, we see an increase in the instability rate. The maximum of the instability rate is obtained at a smaller wavelength than in the case of a tangential discontinuity, and in general, the maximum of the rate is higher than the one obtained for tangential discontinuity. This result is easy to interpret because at contact discontinuity, only the horizontal component of the equilibrium magnetic field tends to stabilise the plasma.

In the case of a tangential discontinuity, there is always a critical wavenumber below which the solution is purely imaginary and hence gives rise to instability. However, above this critical value, the solution is real and the wave is propagating. In contrast, for the contact discontinuity, this is no longer the case; solutions are unstable for all values of wavenumber, k. Moreover, for very large wavenumbers, the imaginary part of the frequency tends towards a fixed value given by

$$\begin{array}{c}\hfill \omega \to \frac{igsin\theta (1-d)}{v_{\mathrm{A}1}{d}^{1/2}(1-{sin}^2\alpha {cos}^2\theta )(d^{1/2}+1)}·\end{array}$$(29)

These values are shown in grey in Fig. 2. The figure also shows that the instability rate decreases for higher wavenumber.

Interestingly, the behaviour of the imaginary part of the frequency (and the disappearance of the critical wavenumber) is similar to the results obtained by Diaz et al. (2014), who studied the effect of partial ionisation on MRT instability in a single-fluid approximation. They attributed the change in the imaginary part of the frequency to the ambipolar diffusion in the induction equation, that is, the modification was due to the presence of neutrals that can diffuse in the perpendicular direction to the ambient magnetic field. This suggests that changes occurring in the transversal direction (relative to the interface) notably modify the behaviour of the instability increment. One important difference is that for the partially ionised case considered by Diaz et al. (2014), the instability rate tends to zero for high wavenumbers, but in our case, this quantity never reaches the zero value.

In the other limit, when the propagation direction is along the y-axis (α = π/2), the only wave that can propagate is the gravity surface wave. The magnetic field begins to have an effect on the propagation characteristics of waves because the propagation direction inclines towards the x-axis.

We investigated how the inclination angle of the magnetic field, θ, affects the stability of incompressible waves propagating along the interface. Now, the value of the wavenumber was fixed at k = 10⁻⁵ m⁻¹ (a typical value for oscillations in prominences), and we chose a number of characteristic values for the propagation angle, α. The numerical result of our analysis is shown in Fig. 3. Again, the mode that has the smallest instability rate is the one that propagates strictly along the x-axis, and this rate shows a pronounced increase for lower values of θ, after which this rate saturates and becomes independent of the inclination angle of the magnetic field. With increasing propagation angle, the instability rate increases, meaning that the amplification time is reduced. This result can be understood when we keep in mind that with increasing propagation angle, the magnetic tension has less effect on the stabilisation of the interface. It is clear that regardless of the propagation angle of the waves, the instability increments tend to a steady value of approximately ℑ(ω) = 0.055 s⁻¹.

Figure 3 shows that regardless of the value of the propagation angle of the waves, all modes will tend towards the same instability rate for a large field inclination angle. In order to translate our results into observable quantities, we show the contour plots of the inverse instability rate (the growth time, in minutes) for a particular wavenumber, density ratio, d, and Alfvén speed in terms of the inclination angle, θ, and propagation angle, α (see Fig. 4), that is, we plot the pair of angles that satisfy the given growth time. This plot shows two distinct behaviours. While left-hand side of the plot shows that for shorter growth time we require higher inclination angle, the mode that appears on the right-hand side shows a different behaviour for very high values of inclination angle. However, this mode is not an Alfvén mode; instead, it is the surface gravity mode that appears for a nearly perpendicular propagation. Similarly to the findings shown in Fig. 3, for a large inclination angle of the magnetic field, the instability growth time becomes independent of θ.

It is well known that magnetic field can stabilise the unstable interface. We therefore investigated the variation in growth time with respect to the wave propagation angle and Alfvén speed, keeping the wavenumber, density ratio, and magnetic field inclination constant (see Fig. 5). Because the density ratio is constant, changing Alfvén speed means a change in the intensity of the magnetic field. An increase in the magnetic field intensity means that the growth time increases, that is, the magnetic field has a stabilising effect for any propagation direction, as expected. In addition, the value of the propagation angle becomes important for stronger magnetic fields.

Finally, we investigated the effect of a changing density ratio on the growth rate of unstable modes, and we obtain that the growth rate decreases monotonically towards zero as the density ratio is increased when the d <  1 condition is maintained. In the case of the tangential MRT instability (but also true for the hydrodynamic case), as the density of the upper plasma increases in comparison to the density of the lower plasma (d decreasing), the system becomes less stable. Taking the limit of d → 0, the instability rate tends towards a fixed value. By taking the solution of S(ω) = 0, when d → 0, the limiting value of the growth rate is found to be $\Im (\omega )\to \sqrt{\mathit{gk}}$. This is clearly independent of either inclination angle or magnetic field strength.

5. Applications to solar prominences

Solar prominences are magnetic features that are suspended in the solar corona. They consist of dense, cold plasma surrounded by tenuous and hot coronal plasma. High-resolution observations show that prominences present threads along which plasma can flow and waves propagate. Such an observation was carried out with the Hinode/SOT instrument by Okamoto et al. (2007), who observed an active region prominence (NOAA AR 10921) in a 0.3 nm broadband region centred at 396.8 nm. They found that a multi-thread prominence was suspended above the main sunspot. Their analysis showed in-phase oscillatory motion with periods 130–250 s. These authors concluded that the observed waves were propagating or standing Alfvén waves.

Terradas et al. (2008) used these oscillations to carry out a seismological study and determined Alfén speeds in the prominence and corona in terms of the density ratio. The threads were considered to be thin flux tubes (observations showed that they have a small radius) in the presence of a flow that proved to slightly influence the period of waves. Using a few simplifications (e.g., straight flux tubes, homogeneous plasma, longitudinal magnetic field, non-stratified plasma, equal-length threads, and linear approximation), these authors obtained that the minimum Alfvén speed varies between 120 and 350 km s⁻¹. Although not analysed in the above studies, Ruderman et al. (2014) found that the typical lifetimes of the threads under investigation is 10 minutes. If the MRT is responsible for these short lifetimes, then we may assume that the instability time is approximately equal to the thread lifetime. This gives us the necessary information to find the magnetic field inclination, θ, in terms of the propagation direction α, for a given density ratio.

Based on the observations by Okamoto et al. (2007), we used the analysis presented by Terradas et al. (2008; their Eq. (3)) to determine the wavenumbers and Alfvén speeds for three possible density ratios (d = 0.1,  0.2,  0.5). With these values, we numerically solved Eq. (26), and this determined the pair of values for magnetic field inclination and propagation direction that satisfies the observed variables and the governing Eq. (26). We note that not all of observable data sets have solutions for any α and θ values, suggesting that only certain density ratios may be possible in these circumstances, thus giving even more information about hard-to-observe variables. The solutions in terms of magnetic field inclination and propagation direction are displayed for all six observed threads in Fig. 6, and the relevant data, including the maximum possible inclination angle (when the wave is directed in the same plane as the magnetic field i.e. x-propagating), are shown in the table.

Fig. 6. Solutions of the dispersion relation for a given growth time in terms of magnetic field inclination (θ) and propagation direction (α) for six observed prominence threads and three possible density ratios. Observational data were adapted from Okamoto et al. (2007).

It is clear that for all threads that the 10 min growth time can be satisfied only for particular values of density ratio and a particular combination of field inclination angle and propagation angle. Thread 1 shows only solutions for a density ratio of d = 0.1, with the maximum field inclination only about half a degree that is attained when the wave propagates along the x-axis. Threads 4, 5, and 6 all have solutions for d = 0.2 with a maximum inclination of approximately 1°. In addition, thread 4 has a solution for d = 0.5 with a noticeably higher maximum inclination of 4°. Thread 2 shows a higher maximum θ for the d = 0.2 solution than the other three, of 3°, and a lower maximum θ for the d = 0.5 solution than thread 4 of approximately 1°. Finally, thread 3 has a solution only for d = 0.5, but with much higher maximum inclination than any of the other threads, with θ = 25°. Although all of the magnetic field inclinations we found are relatively low (except for prominence 3), the inclination is only zero if the propagation direction is almost perpendicular to the magnetic field, which would be very unlikely for MHD waves in the solar atmosphere. This also gives a justification for the study carried by Vickers et al. (2018), where small magnetic field inclinations were considered for MHD waves at a contact discontinuity.

In order to expand our analysis, we also considered the case of plumes, which are seen to develop in prominences. In particular, we studied two plumes observed by Mishra & Srivastava (2019), where growth rates of plumes were estimated to be 1.32  ±  0.29  ×  10⁻³ s⁻¹ and 1.48  ±  0.29  ×  10⁻³ s⁻¹. Using these values, along with the approximate wavelength, densities, and Alfvén speed found in Mishra & Srivastava (2019), we performed an inversion using Eq. (26), and the results are are shown in Fig. 7. For this case, we find a much higher maximum inclination of 22−41°. This shows the variation of inclination that is possible in prominences. It also seems logical that hte inclination would be greater after plumes begin to develop, so that our results are consistent with this assumption.

Fig. 7. Solutions of the dispersion relation for a given growth time in terms of magnetic field inclination (θ) and propagation direction (α) for two observed prominence plumes. Observational data were taken from Mishra & Srivastava (2019).

We should note that in the context of solar prominences, compressibility would have a pronounced effect in stabilising the interface. A model including compressibility is therefore required to give more pertinent results. The density ratios in prominences in general also give lower d values than those considered here, which destabilises the system further. We therefore assume that solutions would be of a similar order of magnitude to those considered here.

6. Conclusions

Through its importance in re-organising the plasma and formation of small-scale structures, the magnetic Rayleigh–Taylor instability is one of the most frequently studied instabilities in solar and solar-terrestrial plasmas. We here studied the effect of an oblique magnetic field that crosses an interface in the development of incompressible MHD Rayleigh–Taylor instability. In this configuration, the sharp interface that separates plasma regions with different densities forms a contact discontinuity.

We showed that in contrast to well-known results for tangential discontinuity, for even a small inclination of the magnetic field, the system is unstable for any wavenumber, that is, there is no longer a critical wavenumber above which the system is stable. The instability rate is also higher than for the tangential case, but it is lower than for the purely hydrodynamic case and is reduced with increasing Alfvén speed because the magnetic tension has a greater component that stabilises against gravity. The instability rate was found to be considerably higher for perturbations perpendicular to the plane in which the magnetic field is inclined than for perturbations in the direction of the field, and there is a smooth change between these two extremes because the magnetic tension has a less pronounced effect as a restoring force.

The high-density gradients at the edges of prominences, between the low-density coronal plasma and the higher density plasma of the prominence itself (readily modelled on the small scale by a sharp interface) is often intersected by inclined magnetic fields. Therefore, we applied our results to model MRT instabilities in solar prominences. These results suggest that even with a high magnetic field strength, instabilities are able to develop for perturbations of any wavelength. Thus, inversion techniques relying on the critical wavelength may not always be valid because there is no such critical wavelength for the contact discontinuity. This also gives us information about where in the prominences plumes are most likely to develop because contact discontinuities appear most frequently at the sides of prominences.

The dispersion relation derived in the current study may be used to perform an inversion technique on observations of oscillating prominences, and this can give us information on the angle of inclination of the magnetic field and direction of the wave propagation. A simple and approximate analysis of six prominence threads observed by Hinode/SOT was performed, and the results confirmed that there is a combination of the two angles that can account for the observed instability growth time. In the cases we studied, the typical inclination angles (where waves propagate in the direction of the plane of inclination) were found to be 1−4°. The inversion also provided information about possible density ratios between the prominence plasma and sparse coronal plasma for each prominence thread. When we applied this inversion to two prominence plumes, observed by Mishra & Srivastava (2019), the average maximal field inclination angle (for perturbations directed in the plane of inclination) was found to be approximately 30°.

This study has provided a first step in the understanding of magnetic Rayleigh–Taylor instabilities at a contact discontinuity. This study used a linear technique, which gives relevant information about the onset of instabilities, but as the instability evolves, non-linear effects will need to be taken into account, including secondary Kelvin–Helmholtz instabilities as the plumes develop, for instance. To take these into account, a time-dependent analysis is required, which, due to its complexity, must be performed numerically. We restricted our investigation to a local analysis, where the effect of gravitational stratification has been neglected, but this approximation might not always be valid because at later stages, as plumes begin to develop, for example, gravitational stratification will become more relevant. Compressibility is known to stabilise many plasma configurations, while partial ionisation has been shown to aid stability in certain circumstances and destabilise in others. These two effects should therefore also be taken into account to obtain a more detailed insight into the effect of field inclination on MRT instabilities.

Acknowledgments

The authors acknowledge the help of M.S. Ruderman. E.V. acknowledges the support by Science and Technology Facilities Council (STFC). R.E. is grateful to Science and Technology Research Facilities Council (STFC, grant number ST/M000826/1) UK for the support received.

References

All Tables

All Figures

Fig. 1. Sketch of the equilibrium configuration. The equilibrium state consists of a surface separating two regions, each with different density. The magnetic field is uniform throughout both regions and is inclined within the (x, z)-plane at an angle, θ, with respect to the x-direction. The configuration is invariant in the x- and y-directions. Perturbations are described by the wave vector, k, in the (x, y) plane.

In the text

Fig. 2. Solutions for the dispersion relation for waves propagating in the (x, y)-plane for two propagation directions (α = 0 solid line, and α = π/4 dashed line) with respect to the wavenumber, k. The density ratio is d = 0.5, the Alfvén speed of the lower plasma is vA = 10 km s−1, and the magnetic field inclination angle is θ = 5°. The upper panel shows the real part of the frequency, while the imaginary part is plotted in the lower panel. The limiting value of ℑ(ω) when k → ∞ is shown by the grey horizontal lines in the lower panel. For illustration we also show the real and imaginary part of the frequency in the case of a tangential discontinuity (θ = 0), plotted here in green.

In the text

	Fig. 3. Imaginary part of the frequency for waves propagating in the (x, y) plane for several propagation directions with respect to the magnetic field inclination, θ. The density ratio is assumed to be d = 0.5, the reference Alfvén speed is vA1 = 10 km s−1, and the value of the wavenumber is fixed at k = 10−5 m−1.
In the text

	Fig. 4. Growth time for k = 10−5 m−1, d = 0.5 vA = 10 km s−1 in terms of propagation angle, α, and field inclination, θ.
In the text

	Fig. 5. Growth time for k = 10−5 m−1, d = 0.5 θ = 5° in terms of propagation angle, α, and Alfvén speed of the lower plasma, vA.
In the text

	Fig. 6. Solutions of the dispersion relation for a given growth time in terms of magnetic field inclination (θ) and propagation direction (α) for six observed prominence threads and three possible density ratios. Observational data were adapted from Okamoto et al. (2007).
In the text

	Fig. 7. Solutions of the dispersion relation for a given growth time in terms of magnetic field inclination (θ) and propagation direction (α) for two observed prominence plumes. Observational data were taken from Mishra & Srivastava (2019).
In the text