© ESO 2021

1. Introduction

The Rayleigh-Taylor instability (RTI) is present in many astrophysical systems, and it is frequently observed in solar prominences (Berger et al. 2008, 2017). RTI is also present in laboratory plasmas and can play an important role in inertial and magneto-inertial confinement fusion experiments (see review by Zhou 2017a, 2017b, and references therein).

On the Sun, prominences are defined as clouds of chromospheric material sustained in the solar corona by magnetic forces. They are colder and denser than the surrounding corona. As a prominence sits above the lighter corona, a small perturbation at the interface between the prominence and the corona may grow without bound, giving rise to the RTI. Small-scale upflows, also called plumes, have been observed at the top border of prominence bubbles by Berger et al. (2017). The plumes are highly dynamic columns of plasma, which rise from below the prominence and propagate upwards with speeds of ≈15 km s⁻¹. It is believed that they appear due to the RTI (Berger et al. 2010, 2017; Orozco Suárez et al. 2014).

The cold prominences contain a large fraction of neutral particles that get ionized when they enter the hotter corona during the development of the RTI. Some observational studies show hints in support of the decoupling in the velocity of ions and neutrals in prominences (Gilbert et al. 2007; Khomenko et al. 2016; Anan et al. 2017; Wiehr et al. 2019), however, these observations are at the edge of the observational capabilities of the Sun. In Earth’s ionosphere, plasma is also partially ionized, and many phenomena related to the uncoupled behavior of ions and neutrals are similar in the solar atmosphere and in the ionosphere (see e.g., the review in Leake et al. 2014). In this work, we apply a two-fluid model to study the effect of the neutrals on the development of the RTI at the interface between the corona and a solar prominence thread, with a smooth transition between them, in order to resolve the effects of charge-neutral interactions.

One of the important parameters of the RTI, both in partially and in fully ionized plasma, is its growth rate. The growth of the RTI can be separated into two phases: linear and nonlinear. During the initial linear evolution, different modes grow independently without interaction; however, mode-coupling effects become important with growing mode amplitudes and the instability evolves into the nonlinear phase when the bubbles (structures that rise) and spikes (structures that fall) form.

Most of the analytical studies of the linear phase of the RTI have been carried out based on the assumption of a discontinuous density profile. In the absence of the magnetic field (the hydrodynamic case), when a system composed of a heavier fluid with mass density, ρ₂, which is initially on top of a lighter fluid with mass density, ρ₁, is perturbed, an instability develops for all the wavelengths of the initial perturbation (see Chandrasekhar 1961). With homogeneous densities, under the incompressible and inviscid assumption, and when the transition between the two fluids has negligible width (discontinuous density profile), the growth rate in the linear approximation is a monotonically increasing function of the wave number k:

(1)

where,

(2)

is the Atwood number (Chandrasekhar 1961).

Under the same assumptions, a component of the magnetic field in the perturbation direction, considered to be the x direction (see Fig. 1, which illustrates the setup used in our simulations and uses the same conventions), stabilizes small scale perturbations with wavelengths smaller than the critical wavelength (Chandrasekhar 1961). The influence of magnetic field shear at the discontinuous density interface was recently studied by Ruderman et al. (2014, 2018).

Fig. 1. Sketch of the initial configuration. The simulation domain (filled with blue color) is contained in the XOZ plane, where x varies between −L0 and L0 and z varies between −4L0 and 4L0. The magnetic field, B (red lines), is contained in the XOY plane (filled with gray color). At the height z = 0, the magnetic field only has the component, By, and it is slowly rotated above and below z = 0 by an angle which depends on height, keeping its modulus constant. The value of the angle is anti-symmetric with respect to z = 0. The maximum angle obtained at the top and the bottom of the atmosphere is indicated by θ0. When the field is perpendicular, the direction of B is along the y-axis for all the heights, θ0 = 0°. For the sheared magnetic field configurations θ0 = 1°. The gravity, g, points in the negative z direction and is indicated by a green line. The direction of the perturbation, k, is shown by a blue line.

In order to capture the two-fluid effects, we have to resolve scales comparable to the collisional mean free path between ions and neutrals. These very small scales are typically below the density gradient scale length for any solar structure. Thus, the density profile cannot be approximated as discontinuous in two-fluid calculations.

The general case of RTI with arbitrary non-uniform equilibrium density profile with a continuous transition in the vertical direction does not have a known analytical solution. However, in the absence of the magnetic field and under the assumption of incompressibility, for an exponential density profile ∝exp(βz) with uniform β, Chandrasekhar (1961) finds the following expression for the growth rate:

(3)

where it is assumed that the fluid is confined between the heights z = 0 and z = d, and m is the vertical mode number. The stratification is unstable when β >  0, with the maximum growth rate for m = 1. We note that the growth rate obtained in Eq. (3) is a monotonically increasing function of k bounded by −ω² = gβ. Thus, the linear growth rate of the RTI, when the transition region between the heavy and light fluids if of a finite width, exhibits a different k-dependence from that calculated with a discontinuous profile discussed above, when the growth rate increases linearly with the wave number as shown in Eq. (1).

The nonlinear effects which occur from interaction of the modes cannot be accounted for analytically in the general case. There are analytical approaches which give qualitative results about the nonlinear (and linear) phase of the RTI in a 3D geometry (Hillier 2016) and statistical results obtained in the late nonlinear (turbulent) phase using a self-similar ansatz. However, the work of Hillier (2016) does not consider partial ionization effects and the hypothesis used in their analytical derivations are restrictive. In the general case, only numerical solutions allow for a substantial tracking of the evolution of the RTI into the nonlinear regime.

Magnetic field strength in prominence threads can be inferred from observations of density contrast of the RTI bubbles in an approximate way. The shear flow at the bubble boundary gives rise to subsequent plumes as a result of the Kelvin-Helmholtz instability (KHI, Berger et al. 2017). The growth rate of the associated KHI depends on the density profile and flow velocities above and below the interface and the magnetic field. By comparing the measured growth rate of the instability to analytic expressions, the magnetic field can be inferred. Berger et al. (2017) found a magnetic flux density across the bubble boundary of ≈10 G at an angle of ≈70° to the prominence plane. The review in Hillier (2018) shows more examples of the determination of the field strengths from the growth rates in several contexts, along with the uncertainties associated with these methods.

The degree of collisional coupling between charged particles and neutrals can influence the linear growth rate and the nonlinear course of the instability, as well as the appearance of bubbles and spikes. However, these aspects have been investigated only scarcely. Most of the analytical studies reported in the literature neglect the partial ionization effects on the development of the RTI in solar prominences. An analytical calculation of the linear growth rate of the RTI in the two-fluid model has been done for the discontinuous density profile in a 3D geometry by Díaz et al. (2012). In the regime studied in that work, the collisions between neutrals and charge were found to decrease the growth rate at all scales. Viscosity is also known to act in a similar way (Chandrasekhar 1961; Livescu 2004). In a later work, Díaz et al. (2014) introduced interaction between neutrals and charge through the ambipolar term in the generalized induction equation. The authors concluded that the neutrals decouple from the charge at small scales, below the collision scale, eliminating the stabilizing effect of the magnetic field on the neutrals, which removes the cut-off imposed by the magnetic field.

Most of the numerical simulations of RTI have been done using the single-fluid magnetohydrodynamic (MHD) approach (Anuchina et al. 2004). Only a handful of numerical studies have included the partial ionization effects through the ambipolar term (Arber et al. 2007; Khomenko et al. 2014). The latter studies were done in the context of magnetic flux emergence in the chromosphere and for the RTI developed at the prominence-corona transition region (PCTR).

Numerical studies of RTI using the two-fluid approach with application to solar prominences are extremely scarce (Leake et al. 2014). In the work of Leake et al. (2014), only one case was studied and only the initial evolution of the instability was considered, with no in-depth analysis of the results. Recently, Hillier (2019) studied Kelvin-Helmholtz instability numerically in the two-fluid approach, but this idealized study cannot be applied to the prominence conditions and the type of the instability investigated by Hillier (2019) is different. Most of the studies, analytical and numerical, omit the effects of viscosity, thermal conductivity, and ionization-recombination.

The overarching purpose of this study is to explore the influence of and possible observables associated with incomplete collisional coupling between charge and neutrals in the partially ionized prominence plasma during the linear and nonlinear phases of the RTI. To do so, we employ a 2.5D two-fluid model of the plasma and, for the first time, we perform a systematic study of the RTI including elastic collisions, ionization and recombination processes, viscosity, and thermal conductivity in the model. In order to capture the two-fluid effects that occur on small spatial scales, we perform simulations at an adequate resolution and consider the equilibrium with a realistic transition between the prominence and the coronal part of the structure, where the properties of the plasma change in a continuous way. There have been a limited number of analytical and numerical studies examining properties of RTI in smoothly varying equilibria in the presence of magnetic field, and thus the novelty of this work includes incompressible MHD linear analysis of RTI in such equilibria. This single-fluid analysis is necessary to begin to understand the more complex two-fluid RTI evolution in nonlinear solar prominence simulations.

This study starts with the same background model as in Leake et al. (2014). We study how the development of RTI in the two-fluid model depends on the parameters of the background atmosphere, such as the density contrast and the magnetic field strength, because these parameters can be, in principle, extracted from observations. From the simulations, we extract the linear growth rates, study the development of bubbles, spikes and the associated magnetic structures in the nonlinear phase, as well as of the flow decoupling between the neutral and charged fluids.

We describe the configuration and plasma parameters of the simulations in Sect. 2, followed by an overview of the simulation dynamics and evolution of different scales in Sect. 3. We explore the analytical solution limits under the single-fluid incompressible MHD approximation for smoothly varying density and magnetic field profiles in Sect. 4. We then study the linear behavior of RTI using semi-analytical calculations for the assumed background atmosphere and in the two-fluid simulations in Sect. 5. In Sect. 6, we study the decoupling in the velocity fields and in Sect. 7, the magnetic structure formation. We discuss our results and present our conclusions in Sect. 8.

2. Description of the problem

We model the RTI at the border of a thin prominence thread. To approximate this situation, we adopt a configuration similar to that described by Leake et al. (2014), shown in Fig. 1. We use 2.5D geometry where the prominence is contained in the XOZ plane (filled with blue), above the height z = 0, and sheared magnetic field, contained in the XOY plane (red lines). This configuration is described below by Eq. (13). The other configuration used for the simulations presented in this paper is the magnetic field, which is unidirectional along the y-axis. For this configuration, θ₀ = 0°, as in the sketch in Fig. 1.

2.1. Background atmosphere

The equilibrium atmosphere is uniform in the x-direction. The equations that describe the number density of the ions and neutrals, the background temperature, the magnitude of the magnetic field, and the pressures in equilibrium configurations are (Leake et al. 2014) the following:

(4)

(5)

(6)

(7)

(8)

where the ion number density at z = 0 is n₀ = 10¹⁵ m⁻³; the peak neutral number density, reached at z = L₀/2 is n_n00 = 10¹⁶ m⁻³; the background temperature of the corona T_b = 2.02 × 10⁵ K; the neutral number density corresponding to the corona temperature (T_b) is n_nb = 3.5 × 10⁹ m⁻³; B₀₀ = 10⁻³ T is the gravitational scale height of the charged particles of H_c = 2k_BT_b/(m_Hg). The characteristic length scale is L₀ = 1 Mm. The plasma is calculated at z = 0 using B₀₀ as the value of the magnetic field and has the value of β_p ≈ 1.4 × 10⁻². The function f is defined as follows:

(9)

The temperature profile, shown in Eq. (6), is described using the function f. The value of L_t = 20, which appears in Eq. (9), has been chosen so that the temperature is closest to the values observed in the Sun and the ionization fraction ξ_i = ρ_c/(ρ_c + ρ_n) = 0.091 remains small (see Leake et al. 2014). Figure 2 shows the background profiles for the density, pressure, temperature and the module of the magnetic field. The neutrals and charged species have the same temperature profile, and are assumed coupled enough by collisions, so that the prominence, composed essentially by neutrals, is supported by the gradient in the magnetic pressure. The equilibrium for the density of neutrals represents a density enhancement at z = L₀/2, where it attains the maximum, n_n00, as from Eq. (5). The minimum of the temperature is also attained at z = L₀/2, and the maximum temperature is considered to be the coronal temperature (T_b).

Fig. 2. Parameters of the initial model atmosphere as a function of the vertical coordinate, z. Top left: density of neutrals in the models with the original (red) and enhanced (dashed green) density contrasts, and density of charged particles (blue). Top right: same for the pressure of neutrals and charged particles. Bottom left: temperature of both the charged particles and neutrals. Bottom right: modulus of the magnetic field in the models with the original (red) and enhanced (green) density contrasts.

The equilibrium density of the charged particles is gravitationally stratified. Using these density and temperature profiles, the pressure profiles are obtained from the ideal gas law, Eq. (8). The magnetic field profile is obtained by integrating the magneto-hydrostatic equation where the total density and pressure are used. The integration constants permits the scaling of the magnetic field around some chosen value, B₀₀ in Eq. (7) shown above.

2.2. Perturbation

To initiate the instability, we use a perturbation in the density of neutrals, located at height z = 0:

(10)

where δ is the perturbation amplitude and r_x is the form of the a perturbation which depends on the x coordinate. The equilibrium atmosphere is homogeneous in the x direction and the dependence of the perturbation on the x coordinate gives rise to the choice of the ansatz of separable solutions ∝f(z)exp(iωt)exp(−ik_xx). Thus, the perturbation is propagated in the x direction, indicated by the vector k = (k_x, 0, 0) in Fig. 1.

While the initial perturbation should not have any influence on the linear growth rates of individual modes, the evolution in the nonlinear phase may depend on it. In most of the simulations reported in this paper, we use a white noise perturbation, defining r_x in Eq. (10) as a random number in the interval [−1,1] for all the points. In addition, we compare the linear growth rate of modes with the white noise (“WN”) perturbation with a simulation where the perturbation in the x-direction is generated as a multi-mode perturbation (“MM”) specified by:

(11)

where m_x is the number of points in the x direction and ϕ_n is a random phase.

2.3. Two-fluid evolution equations

We solve the two-fluid set of equations, described in Sect. 2.1 of Popescu Braileanu et al. (2019), referred to in this subsection as EqP. To do so we use the two-fluid version of the MANCHA3D code, MANCHA3D-2F. The single-fluid version of this code is described in Felipe et al. (2010), González-Morales et al. (2018), and the two-fluid version preserves several numerical features of MANCHA3D, such as hyper-diffusive algorithms, filtering, and PML layer, see Popescu Braileanu et al. (2019).

The continuity equations of neutrals and charged species (Eq. (1) in EqP) are coupled by the inelastic collisional term (Eq. (2) in EqP) which models ionization and recombination processes. The momentum (Eq. (3) in EqP) and energy (Eq. (9) in EqP) equations are coupled by collisional terms (Eqs. (6) and (10) in EqP, respectively), which include all the elastic and inelastic contributions. The charge-exchange reactions are introduced via the elastic collisional parameter (Eq. (A.8) in EqP). We include the viscosities in the pressure tensor (Eq. (8) in EqP). The thermal conductivity is considered only in the neutral energy equation. The expressions for the viscosity and thermal conductivity coefficients are given in Eq. (A.9) from Appendix A.2 in EqP. We use the ideal Ohm’s law (right hand side of Eq. (14) in EqP is zero).

2.4. Parameters of the simulations

Table 1 describes all the simulations used in this work and includes the abbreviated names of the simulations, the figures where they appear, and the summary of their parameters. There are four input parameters that are changed throughout the simulations:

Magnetic field inclination. The background magnetic field is contained in the XOY plane. We use two types of configurations, namely, the magnetic field perpendicular to the plane defined by the direction of the perturbation and the gravity XOZ, and a slightly sheared magnetic field, that is,

the perpendicular magnetic field

(12)

These simulations are referred in the table with the keyword “P;”

the sheared magnetic field

(13)

We use a fixed value of θ₀ = 1° and L_s = L₀ in Eq. (13) throughout the simulations. These simulations are marked in the table by “L1”.

Density contrast. The value of the peak density of neutrals is changed in some of the simulations by increasing the value of n_n00 in Eq. (5) by a factor of 10. These simulations appear marked with the “NN” keyword in their name. The Atwood number values calculated using the total density at the heights −L₀/2 and L₀/2, that is, ρ₂ = ρ₀(z = L₀/2) and ρ₁ = ρ₀(z = −L₀/2), are A = 0.72 for the original density case and A = 0.85 for the enhanced density case.

Magnetic field magnitude. The value of B₀₀ in Eq. (7), which describes the magnetic field configuration, is increased by a factor of 3 (the “B” keyword).

Amplitude of the perturbation. The amplitude of the perturbation is varied through the parameter δ in Eq. (10) above. We use the value δ = 10⁻² in most of the simulations. But for the study of the linear growth rate, we use a smaller amplitude by a factor of 100, δ = 10⁻⁴, in order to avoid the nonlinear effects. The names of these simulations in the table have an additional suffix “−S”.

Figure 2 shows the background profiles for both the original and enhanced value of the density contrast. We note that when the density peak is increased, larger magnetic pressure gradient is needed to support the prominence. Therefore, the magnetic field profile also changes, as shown in Fig. 2 in the bottom right plot.

We choose the domain size L_x = 2L₀, L_z = 8L₀. We use a uniform grid for the two dimensions and a resolution of 512 × 2048 points; this gives the values of the grid cell sizes: d_x = d_z = 3.9 km. We choose the point (0,0) in the middle of the plane XOZ, so that x varies from −L₀ to L₀, and z varies from −4L₀ to 4L₀.

In all our simulations, we avoid using hyper-diffusion to stabilize the simulations. In order to remove small scale noise, we use the filtering method described in Parchevsky & Kosovichev (2007). We use periodic boundary conditions in the x direction, reflecting boundary conditions for the vertical velocities of neutrals and charged particles (velocities are zero at the boundaries) and open boundary conditions for the rest of the variables (space derivatives in the vertical direction are zero at the boundaries) in the z direction. Similarly to Terradas et al. (2015), the reflecting boundary condition for the velocities at the bottom of the atmosphere plays an important role in supporting the prominence against gravity.

2.5. Collisional dissipation coefficients and scales

Our simulations are set to have sufficient resolution to capture the scales impacted by diffusion due to neutral viscosity or neutral-charge collisions, or both. We evaluated these scales using the parameters of the background model atmospheres by computing the coefficients of the viscosity of charged particles () and neutrals (), and the neutral thermal conductivity (), as well as the neutral-ion (ν_ni) and ion-neutral (ν_in) collision frequencies. The latter two are also compared to the time-scales associated with the fast magneto-acoustic wave (ω_fast), the in-plane Alfvén wave (), and the neutrals sound wave (ω_cn) within the background atmospheres as characteristic of dynamical time-scales in the system. These are calculated as follows,

(14)

where γ is the adiabatic index, and ρ_tot = ρ_n + ρ_c is the total density. The collisional parameter α has been defined in Eq. (A.8) and ξ_α and K_α in Eq. (A.9) in Popescu Braileanu et al. (2019). The density of the neutrals is very low in the corona, and the temperature is high, therefore, the coefficients and are very high. To overcome the very small time-step imposed by such large coefficients of neutral viscosity and thermal conductivity, the values of these coefficients have been limited to 10⁹ m² s⁻¹ in order to speed up the computation. This limit does not significantly impact the regions of the domain where the instability develops and therefore does not make qualitative difference for the outcome of the computations. An improved implementation of the viscosity and thermal conductivity operators for low collisionality systems will be considered in future work.

The left panel of Fig. 3 shows the values of the neutral viscosity and thermal conductivity coefficients, limited as described above, and the viscosity coefficient of the charged particles for the two background atmosphere neutral density profiles. We observe that the neutral viscosity and thermal conductivity coefficients have almost the same profile and the charge viscosity coefficient is an order of magnitude smaller than the neutral viscosity coefficient.

Fig. 3. Left: coefficients of the viscosity of charged particles (, blue line) and neutrals (, orange lines), and the thermal conductivity of neutrals (, green lines), in units of m2 s−1, as from Eq. (14). The neutral viscosity and thermal conductivity coefficients are limited to 109 m2 s−1. Right: ion-neutral (νin, blue lines), neutral-ion (νni, orange line) collision frequencies, and frequencies associated with waves: fast magneto-acoustic (ωfast, green lines), in-plane Alfvén for “L1” shear (, red lines), acoustic in the neutral fluid (ωcn, black line) as functions of height. Values for the model with the original density contrast model are shown with solid lines, and for the enhanced density contrast with dotted lines.

Right panel of Fig. 3 shows the ion-neutral, ν_in, and the neutral-ion, ν_ni, collision frequency profiles as functions of height for the background atmospheres. To compare these to ideal MHD timescales of interest, we also plot the frequencies associated with the largest scale wave modes in the x-direction for the fast magneto-acoustic wave, ω_fast, the in-plane component of the Alfvén wave for “L1” magnetic field shear, , and the neutrals sound wave, ω_cn.

We observe that in the immediate neighborhood of the prominence thread between z ≈ −1 Mm and z ≈ 2 Mm, collision frequencies between the charged and neutral fluids are higher than the large-scale ideal MHD frequencies, and, thus, large-scale flows are expected to be well coupled. However, it is also apparent that the frequency separation is not so large within the thread, and disappears altogether outside the core of the thread, so that smaller scale and higher frequency ideal MHD dynamics is likely to be strongly impacted by the two-fluid effects.

3. Global overview of simulation dynamics

Figure 4 shows time series of snapshots of the four simulations, L1-WN, L1-WN-NN, L1-WN-NN-B, and P-WN-NN. This figure illustrates the development of the instability in cases with different density contrast, magnetic field strength and orientation. The moments where the snapshots are shown to correspond to times indicated in the following Fig. 5 and taken at approximately similar evolution stages of the RTI in each case. The three simulations with a sheared magnetic field (L1-WN, L1-WN-NN, L1-WN-NN-B) have been initialized with exactly the same perturbation, so that it is easy to compare them. The velocities of neutrals and charged particles are indicated by arrows in the panels of the corresponding densities. For the simulations with the sheared field, the bottom panels for each simulation show the in-plane magnetic field lines plotted as iso-contours of the absolute value of the y-component of the magnetic potential obtained by numerical integration. The iso-contour lines span uniformly the range between 0.6 and , where is the maximum value of |A_y| and they are darker for larger values.

Fig. 4. Time series of the snapshots of neutral and charged densities for the simulations L1-WN (top left group of panels), L1-WN-NN (top right group of panels), L1-WN-NN-B (bottom left group of panels), P-WN-NN (bottom right group of panels). The in plane (x and z components) velocities of neutrals and charged particles are represented by orange arrows in the plots which show the densities of neutrals and charged particles, respectively. For the simulations with sheared magnetic field (L1-WN, L1-WN-NN, L1-WN-NN-B), the iso-contours of the absolute value of the y-component of the magnetic field potential, |Ay|, are shown at the bottom panels. The iso-contours are equally spaced between 0.6 and , where is the maximum value of |Ay|.

Fig. 5. Temporal evolution of individual harmonics, computed from the vertical velocity of neutrals, as a function of time. Different color lines show the Fourier amplitudes of the harmonics, grouped in sets of four, from n = 1 to n = 48. The amplitudes are averaged between heights −L0 and L0. Different panels are for different simulations: L1-WN (top left); L1-WN-NN (top right); L1-WM-NN-B (bottom left); P-WN-NN (bottom right). The vertical dashed pink lines mark the times of the snapshots shown in Fig. 4.

By comparing the snapshots of L1-WN to the snapshots of L1-WN-NN, we can observe that small-scale structures are more abundant when the density contrast is increased. Increasing the magnitude of the magnetic field has the opposite effect. As the density of the neutrals increases with height across the interface, and that of charged species decreases, the spikes contain mostly neutral fluid and the bubbles contain mainly charged particles. We observe this effect in the density contrast between charged particles and neutrals in the snapshots for all the simulations in Fig. 4. The collisional coupling is strong and the neutrals drag the charged particles during the development of the instability.

In order to see how the changes in the equilibrium parameters influence the development of different scales in the system, we measure the time evolution of the horizontal modes averaged over the vertical extent of the unstable region. We compute the amplitudes of the modes as the amplitudes of the Fourier coefficients obtained by performing a FFT of a given quantity (for example vertical velocity of neutrals or charged particles), in the x direction, at each time moment. We then average the amplitudes in the vertical direction between z = −L₀ and z = L₀, and over four consecutive modes. Figure 5 shows the natural logarithm of the Fourier amplitudes for the first 48 modes. The average over four modes and the limit of 48 modes are chosen for clarity of presentation.

Noting that the time extent of the horizontal axes in all the panels of Fig. 5 is different, by comparing top-left (L1-WN) and top-right panels (L1-WN-NN) we observe that the increase in the density contrast increases the growth rate of all modes. By comparing top-right (L1-WN-NN) and bottom-right (P-WN-NN) panels, we observe that the magnetic field shear has a stabilizing effect; and by comparing top-right (L1-WN-NN) and bottom-left panel (L1-WN-NN-B), we observe that, in the presence of shear, the increase in the magnetic field magnitude similarly decreases the growth rate of all modes.

In all cases shown in Fig. 5, the modes corresponding to the largest scales (modes 1–8) have a clear linear phase. Small scales (modes 45–48) do not have an initial linear growth for any of the four simulations. As discussed further in Sect. 4.2 below, these are the modes affected by the viscosity and ion-neutral collisions. The largest mode number n that grows in the linear phase in the enhanced density case, L1-WN-NN (top right), is n = 33, while for the original density case, L1-WN (top left), it drops to n = 25. Increasing the magnetic field causes an opposite effect to increasing the density contrast, the largest n growing in the linear phase becomes n = 29. We note that in all cases, at some point in the simulation, the smallest scales begin to grow. This reversal marks the end of the linear phase and the beginning of the nonlinear phase. The small scales involved, which are suppressed during the linear phase by the dissipation effects, are then driven by the energy cascade from low n to high n in the nonlinear phase.

4. Analytical approach

We now explore the expected linear properties of the RTI for the given equilibrium analytically in order to better understand and interpret the simulations, as well as to identify the impacts of the collisional effects on the development of the instability. The following derivations use a single-fluid ideal incompressible MHD approximation of the plasma.

We note that while the derivation below does not take into account the effects of compressibility, it has been shown that compressibility can have a destabilizing effect (Vandervoort 1961; Bhatia 1974; Newcomb 1983; Ribeyre et al. 2005; Xue & Ye 2010). The analytical calculation in the hydrodynamic case for our density profile (not presented in this study) also shows that the compressibility has a destabilizing effect, which is larger for smaller Atwood numbers, similarly to the conclusion of Ribeyre et al. (2005).

Under the assumption of incompressibility, substituting the zero-divergence condition on velocity for the energy equation, the equations become:

(15)

We consider a 2.5D case, where there are no variations in the y-direction, the equilibrium atmosphere is homogeneous in the x direction but not uniform in the z direction, and B_z0 = 0. We use the ansatz of separable solutions with Fourier decomposition of the dependent variables in time and in the x spatial direction, looking for solutions of the form:

(16)

with the general solution in the linear regime being a superposition of such single-mode solutions. Linearizing Eq. (15) and substituting perturbations of the form defined in Eq. (16), we have:

(17)

After manipulating the above we obtain a second order differential equation for v_z:

(18)

where a is defined as:

(19)

where ρ₀ is the total density, calculated from our two-fluid model as:

(20)

where n_i0 and n_n0 have been defined in Eqs. (4) and (5), respectively.

When B_x0 = B_y0 = 0, we recover Eq. (42) from Chandrasekhar (1961). In the hydrodynamic case, neglecting the spatial derivatives in Eq. (18) recovers the Brunt–Väisälä frequency in the incompressible limit.

Equation (18) does not have a general known analytical solution for generic continuous density and magnetic field profiles. In order to obtain an approximate solution, Eq. (18) is numerically solved as a generalized eigenvalue problem with fixed boundary conditions, v_z(−4L₀)=v_z(4L₀)=0. In doing so, we look for the fastest growing mode for a given k_x. We refer to this method as semi-analytical and use it for exploring the linear properties of RTI in our equilibrium in Sect. 5.

4.1. Analytical expressions for the growth rate

The solution for the growth rate above is obtained semi-analytically by solving a general eigenvalue problem. However, this method does not give an analytical expression for the growth rate. We can estimate the growth rate by making further assumptions in some limiting cases.

4.1.1. Discontinuous density profile

In our case, the transition between the prominence and the corona is smooth, and the minimum of the density gradient scale length, is non-zero. However, for modes with k_xL_d/2π ≲ 1, we can approximate our vertical profile as a discontinuous profile, using values of ρ₂ = ρ₀(z = L₀/2)≈1.83 × 10⁻¹¹ kg m⁻³, and ρ₁ = ρ₀(z = −L₀/2)≈2.93 × 10⁻¹² kg m⁻³, with ρ₁ and ρ₂ as in Eq. (2).

It can be shown (see Chandrasekhar 1961) that in a magnetized plasma with discontinuous density profile and a uniform magnetic field component parallel to the direction of perturbation, considered here to be the x direction (see Fig. 1), the growth rate becomes:

(21)

The condition for stability −ω² <  0 gives the expression for the critical wavelength λ_c, so that the scales with wavelength λ = 2π/k_x <  λ_c do not grow:

(22)

where is the characteristic in-plane Alfvén speed at the interface. Thus, Eq. (21) can be used to estimate the expected RTI growth rate for modes with wavelength λ >  max(L_d, λ_c).

Figure 6 shows the cutoff wavelengths λ_c, defined in Eq. (22), as a function of height calculated using the values of ρ₁ and ρ₂ specified above for different magnetic field shear length scales, L_s = [L₀/10, L₀/2, L₀]. These values provide a visual representation of the local magnetic field stabilization effect. The density gradient scale length (red dotted line), attains its minimum of L_d ≈ 409 km at z ≈ 0, and corresponds to the mode number n_d = L_x/(2πL_d)≈1. The density gradient becomes negative in the region z >  L₀/2. We observe that for RTI modes with λ >  L_d that can be treated using the approximation of a discontinuous interface, the horizontal magnetic field may be expected to begin to stabilize such modes for shear scales of L_s ≈ L₀/10 and smaller.

Fig. 6. Cutoff wavelengths (black lines) due to the in-plane component of the magnetic field as a function of height, computed from Eq. (22) for the equilibrium profile with the values of nn00 = 1016 m−3 and B00 = 10−3 T in Eqs. (5) and (7). Cutoff wavelengths for the following four shear lengths: Ls = [0, L0/10, L0/2, L0] are shown with increasing width of the line for increasing Ls. The red dotted line corresponds to the density gradient scale length Ld.

4.1.2. Exponential density profile

Another limiting case is the case of an exponential density profile, ρ₀ ∝ exp(βz), with a constant, β >  0, extended over height, d. In this case, the growth rate can be calculated exactly (Chandrasekhar 1961) with ω(k) given by Eq. (3) and the vertical extent of the mode v_z(z) given by:

(23)

Using the same heavy and light fluid density values as in the sub-section above, ρ₂ = 1.83 × 10⁻¹¹ kg m⁻³ and ρ₁ = 2.93 × 10⁻¹² kg m⁻³; we match them to an exponential density profile ρ₀ ∝ exp(βz), for z between −L₀/2 and L₀/2, with constant β. We obtain β = ln(ρ₂/ρ₁)/L₀ = 1.83 × 10⁻⁶ m⁻¹. Alternatively, we can calculate β from the local density gradient scale length at z = 0 and we obtain a value of β = 2.44 × 10⁻⁶ m⁻¹. Figure 7 shows the dispersion relations obtained by substituting the two values of β and d = 8L₀, corresponding to the vertical extent of the computational domain, into Eq. (3). We note that for long wavelength modes with λβ ≳ 1, the growth rate begins to decrease with decreasing n due to the nature of an exponentially stratified density profile and the finite vertical extent of the modes. However, for modes with wavelength much smaller than both the density gradient scale length and the vertical extent of the domain, the growth rate becomes independent of n.

Fig. 7. Growth rate obtained from the exact solution for the exponential density profile, Eq. (3) with uniform β = 1.83 × 10−6 m−1 (blue curve marked with “o”) and β = 2.44 × 10−6 m−1 (orange curve marked with “x”). The values of n are between n = 1/64 and n = 50.

The growth rate calculations with an exponential density profile give a rather satisfactory result compared to the semi-analytical calculations presented below in Sect. 5. On small scales, for β = 2.44 × 10⁻⁶ m⁻¹, the value of the growth rate in the high n limit is the Brunt–Väisälä frequency ≈2.58 × 10⁻² s⁻¹.

While the analytical approximation of the dispersion relation given by Eq. (3) appears to be a good match to that obtained for the equilibrium configuration considered here, that does not imply that the RTI eigenmodes given by Eq. (23) are also the eigenmodes for an equilibrium with a vertically limited extent of the RTI-unstable density profile. In particular, we note that the functional form of Eq. (23) for the fastest growing m = 1 mode specifies v_z(z) with non-negligible magnitude over the full extent of the domain and shows no dependence on k_x. At the same time, it is easy to show from Eq. (18) that for an equilibrium with B_x0 = 0:

(24)

and the most unstable eigenmodes v_z(z) will be those with v_z ≈ 0 in the regions of stable stratification, (dρ₀/dz)≤0, corresponding to a maximal value of the right-hand-side integral. It follows that for ρ₀ as given by Eq. (20), we can expect the fastest growing RTI modes to have v_z ≈ 0 for |z|≳L₀/2, and for the functional form of v_z(z) to depend on k_x.

4.1.3. Stabilization by sheared magnetic field

As shown in Sect. 4.1.1 above, for a RTI with k ⋅ B ≠ 0, magnetic field acts to reduce the growth rate or entirely eliminate the instability due to the magnetic field-line tension counteracting the gravity force driving the instability. Thus, in magnetized 3D systems with unidirectional B-field, the fastest growing modes will always be those with k ⋅ B = 0.

For more realistic systems where the magnetic field may be sheared, such as in some of the background equilibria considered here, the situation is more complex. To explore the impact of a sheared field on the growth rate for most unstable RTI modes in a 2D system with k = k_xx̂, we note that the vertical extent of an RTI eigenmode in a magnetized plasma may be limited by the magnetic field tension away from the height where B_x0 = 0.

We now consider spatially varying horizontal magnetic field B_x0(z) such that:

(25)

where b₀ is constant and L_s is the characteristic length of magnetic field shear. We assume the density profile to be exponentially stratified, ρ₀ = ρ₀₀exp(βz), as in Sect. 4.1.2 above.

For sufficiently large b₀, the vertical extent of an RTI eigenmode for a given k_x can then be estimated by balancing the gravitational force driving the instability against the magnetic tension in the z-momentum equation (Eq. (17)) and looking for z = z_d where:

(26)

It follows that

(27)

and it is easy to show that positive solutions of Eq. (27) for z_d only exist when

(28)

implying that the RTI mode is not vertically confined by sheared magnetic field in the long wavelength limit such that Eq. (28) is not satisfied. Equation (28) can be approximately evaluated for our equilibrium parameters showing that, for example, for field shear of L_s = L₀ modes with k_xL₀ <  16.3 are not expected to be confined by the sheared magnetic field.

In the opposite limit of (k_x/β)≫1, Eq. (27) has two roots for z_d with that of (βz_d)≪1 as the one of interest. In this limit, it follows from Eq. (27) that:

(29)

implying that the vertical extent of unstable RTI modes in the presence of sheared magnetic field is inversely proportional to k_x for (k_x/β)≫1. Equation (29) can be evaluated for our equilibrium parameters and field shear of L_s = L₀ to give k_xz_d ≈ 4.91, with proportionally smaller values for smaller L_s.

We now substitute the expression for from Eq. (29) in place of d² in Eq. (3), as that would be the estimated extent of the RTI eigenmode as stabilized by the sheared magnetic field, and take the (k_x/β)≫1 limit, resulting in:

(30)

with m = 1 for the fastest growing mode.

The estimate of sheared field stabilization provided by Eqs. (26)–(30) can be further refined by again applying the ∫ v_z * […]dz operator to Eq. (18), which is now in the presence of a sheared magnetic field. In observing that for a mode that is vertically confined by magnetic field tension (or stable density stratification outside of the instability layer), the boundary integral terms can be neglected, we get the following relationship:

(31)

where the first two integrant terms on the right hand side represent the balance between stabilizing magnetic field tension and destabilizing gravitational stratification discussed above, and the third term on the right hand side represents an additional stabilizing effect of compression of in-plane magnetic field within a vertically confined mode. We note that this magnetic field compression effect is present even under the assumption of incompressible flow used to derive Eq. (18), as incompressible vortex flow across the field can close via flow along the in-plane component of magnetic field thus locally compressing the field without compressing the fluid. We also note that this effect is independent of k_x and is purely stabilizing.

An estimate of the importance of the magnetic field compression effect for magnetic field configurations with different shear scales on a mode with wave number n = (k_xL₀)/π can be obtained by assuming mode’s vertical extent to be z_d given by Eq. (29) and approximating |dv_z/dz|≈v_z/z_d while assuming v_z(z) ≈ ṽ_z of some arbitrary constant magnitude for z ∈ (− z_d, z_d), and negligible otherwise. Further approximating (dρ_z/dz)≈βρ₀₀ for z ∈ (− z_d, z_d) and substituting for B_x0 from Eq. (25), under these assumptions, the right hand side of Eq. (31) can be approximated as:

(32)

indicating that such a mode is stable for any field shear length , with the stability criterion independent of k_x as long as (k_x/β)≫1, as assumed above. For the equilibrium configuration used in this work, this stability criterion is L_s <  0.14L₀.

The above derivation shows that presence of sheared magnetic field with finite shear length can stabilize RTI for sufficiently small shear length, but the stabilization effect is largely independent of the wavelength of the mode as long as the magnetic field is strong enough to confine the mode within the RTI-unstable layer. On the other hand, unlike the derivation in Sect. 4.1.1 for uniform B_x0, if the magnetic field shear length is sufficiently large, in the limit of (k_x/β)≫1, the RTI growth rate appears to be independent of k_x, with the asymptotic value reduced from the unmagnetized value of due to the presence of the magnetic field.

4.2. Modes affected by collisional dissipation

Dissipation terms affect the growth of linear modes on scales similar to their characteristic length scales. We can evaluate which modes will be affected by collisions by estimating the mode numbers at which the collisional dissipative terms in the corresponding equations become comparable to the dominant convective terms.

As discussed in Sect. 2.5, neutral viscosity is the dominant viscous dissipation mechanism for the plasma parameters considered in this work. Taking the sound speed of the neutrals, c_n, as the characteristic flow velocity of neutrals and comparing the convective derivative term to the viscous term in the momentum equation for neutrals, the length scale, and mode number associated with the neutral viscosity can be estimated as:

(33)

The length scales associated with the effects of ion-neutral collisions can be similarly estimated in the limit of strong coupling by balancing the magnetic plus thermal pressure gradient force against the inertial terms and the ion-neutral friction term. Doing so leads to:

(34)

where ν_in, ν_ni, and v_fast have been defined in Eq. (14).

Figure 8 shows the mode numbers n corresponding to the ion-neutral decoupling scale and to the neutral viscosity scale. The modes associated with the scales mentioned above are plotted for the original (solid lines), enhanced density contrast (dashed lines), and for the atmosphere where both the density contrast and the magnetic field are increased (dotted lines). We limit the horizontal axis between z = −1.5 Mm and z = 1.5 Mm, these values are considered to be the maximum extent of the structures that form during the instability development in our simulations. The horizontal dashed and dotted lines are n = 48 and n = 256, and represent the maximum mode number which appear in Fig. 5 and the maximum mode number that can exist in the domain (sampled by two points), respectively.

Fig. 8. Modes corresponding to the neutral viscosity (blue lines), ion-neutral collision scale (orange lines), and neutral-ion collision scale (green lines) as a function of height, computed for the original density profile (solid lines), for the enhanced density profile (dashed lines), and for the enhanced density and magnetic field profiles (dotted lines), from Eqs. (33)–(34). The horizontal dashed line marks the maximum mode number shown in Fig. 5, n = 48, and the dotted line represents the largest mode number that can be resolved in the domain, using two points, n = 256.

We observe that both the viscosity and the ion-neutral collisions affect modes that can exist in our domain, that is, the modes below the horizontal dotted line. However, the scales associated with neutral-ion collisions likely provide the strongest damping mechanism. We note that higher neutral density contrast and, thus, higher peak neutral density, is expected to reduce the impact of collisional dissipation. We also note that stronger background magnetic field is expected to lead to more damping due to collisions between ions and neutrals. Both of these effects are particularly consistent with the expectation from damping due to ambipolar resistivity in a single-fluid approximation.

5. Linear phase of the RTI

In this section, we study the linear phase of the RTI. We compare the linear growth rates obtained semi-analytically in the single-fluid incompressible ideal MHD limit to the analytical estimates presented in Sect. 4 and linear growth rates obtained from the simulations.

The left panel of Fig. 9 shows growth rates as functions of the wave number obtained semi-analytically for the original density contrast and magnetic field strength, but different magnetic field shear scales, as well as for a uniform magnetic field rotated by 1° with respect to the y-axis. To compare these to the analytical estimated in both long wavelength and short wavelength limits, we show asymptotic behavior for n ≪ 1, as well as for n ≫ 1. We observe that for long wavelength modes with k_xL₀ = πn <  1, RTI growth rates become insensitive to the presence of an in-plane magnetic field and converge towards ω = 0 for n → 0. This is consistent with lack of magnetic field confinement of the modes for small k_xL₀ as shown in Eq. (28). We can compare the growth for the perpendicular case (the orange line labeled “perp”) to the dispersion relations shown in Fig. 7. We note that the growth rate curves match the analytical estimate well, based on the local density gradient scale β = 2.44 × 10⁻⁶ m⁻¹ giving a particularly good approximation.

Fig. 9. Growth rates as a function of the mode number, n, in logarithmic scale up to n = 50. The first value is n = 1/64 for the left panel and n = 1 for the right panel. Left: growth rate obtained semi-analytically for the perpendicular magnetic field (orange); sheared field with Ls equal to L0 (red), L0/2 (violet), L0/3 (cyan), L0/4 (lime), L0/5 (blue), L0/10 (dark green); and for uniformly rotated magnetic field by 1° (black). Right: growth rate obtained from the simulations P-MM-S (blue), P-WN-S (orange), L1-MM-S (green), L1-WN-S (red).

Fig. 10. Comparison between the growth rates obtained semi-analytically in the incompressible MHD approximation (left) and the growth rates obtained from simulations (right) and as functions of the mode number, n, in logarithmic scale, limited to n = 50. In both panels, the perpendicular magnetic field cases are shown with solid lines and sheared case are shown with dashed lines. Orange lines: P-WN-S; black lines: P-WN-NN-S; red lines: L1-WN-S; green lines: L1-WN-NN-S, and blue lines for L1-WN-NN-B-S.

The cases where the equilibrium magnetic field has shear (and, thus, an in-plane component of the field in the domain) show the peak value of the growth rate decrease with decreasing shear scale, L_s. This is consistent with the analytical calculation of RTI stabilization by sheared magnetic field described in Sect. 4.1.3 and shown in Fig. 6 for a discontinuous density profile. For k_xL₀ >  1, as we consider equilibria with shorter and shorter shear scale, L_s, we observe that the growth rate is decreasing with L_s but is independent of k_x for a range of L_s values, and is fully stabilized for L_s between 0.2L₀ and 0.1L₀, becoming essentially equivalent to the L_s → 0 limit with uniformly rotated magnetic field. This result is again fully consistent with the analytical derivations presented in Sect. 4.1.3 where the stability criterion was estimated to be L_s <  0.14L₀ due to the stabilization of RTI by a combination of magnetic field tension and compression effects.

In order to calculate the growth rate from the simulations, we perform the following procedure. First, we re-run the simulations with smaller perturbation amplitudes to ensure a long linear phase, see Table 1. For each simulation, we identify a time interval of the linear phase, usually observed to be between 100 s and 200 s. We then compute the Fourier amplitude of v_zn(x) along the x direction at every height, for every snapshot of the series. These Fourier amplitudes were averaged in the interval of heights z₁ <  z <  z₂, defined as the region where the normalized Fourier amplitude is larger than 0.3 by the end of the linear phase. We thus obtain average amplitudes of the harmonics as a function of time. The growth rate is obtained as a linear fit to the natural logarithm of the amplitudes as a function of time.

The resulting growth rates of the first 50 modes are plotted in the right panel of Fig. 9 for two of the equilibrium configurations: “P” with magnetic field perpendicular to the plane of the simulation and “L1” with sheared field with L_s = L₀. To verify the numerical accuracy of the simulations, we use two different types of perturbations and show that the spectrum of linear growth rates of the RTI depends on the equilibrium but is independent of the form of the perturbation. The right panel of Fig. 9 compares simulations with exactly the same parameters other than the initial perturbations being the white noise “WN” (as used everywhere else in this paper) or the multi-mode perturbation “MM” described by Eq. (11). The multi-mode and white noise perturbations give growth rate spectra that match almost perfectly for both magnetic field configurations, verifying our numerical procedure.

By comparing the linear growth rates obtained from the simulations to the semi-analytical calculations, we observe that there is a cutoff for the growth rates obtained from the simulations which is not present in the semi-analytical calculation. However, we also note that the growth rate values match very well for the low n mode numbers for both “P” and “L1” magnetic field configurations. The difference between the idealized semi-analytical calculations and the simulations is due to dissipative effects such as viscosity, thermal conductivity, and collisions between neutrals and ions, which are included in the simulations and act on smaller scales. The reduction in the growth rates around n ≈ 10 and cutoff near n ≈ 30−40 is consistent with the estimates provided in Sect. 2.5 and shown in Fig. 8.

Figure 10 compares the growth rates obtained semi-analytically (left) and numerically (right) for the cases with different density contrasts, magnetic field strength and shear. We can observe that for smaller mode numbers, n ≤ 10, there is good agreement between the simulated and the semi-analytical results, both in the magnitude of the growth rates obtained, and in the ordering of the curves for different cases. For higher n we observe that, similarly to Fig. 9, the growth rates decrease and a cut-off is produced for all simulated cases, unlike the semi-analytical cases.

We observe that the increase in density contrast leads to the increase in both the peak growth rate and the cut-off wave number, n_cd. Increasing the magnetic field has an opposite effect, that is, the values of the growth rates and n_cd become smaller. Collisions between neutrals and ions contribute to the creation of the cut-off observed in the simulations in Fig. 10. As shown in Fig. 8, the lowest modes estimated to be impacted by neutral-ion and ion-neutral collisions, n_νni and n_νin, are below those for the viscosity scale, n_viscn. Therefore, it is expected that the observed cut-off is due to collisions between neutrals and ions. In fact, the ordering of the n_cd values for different atmosphere profiles in the simulations with the sheared magnetic field (L1-WN-S, L1-WN-NN-B-S, L1-WN-NN-S, from the smallest to the largest) is the same as the ordering of n_νin in Fig. 8 (yellow curves).

We obtain further information about the RTI instability in a smoothly stratified atmosphere by studying the eigenfunctions of the instability. In particular, the eigenfunctions of the system provide information on the vertical spatial extent of the perturbations for different modes and can therefore be used to estimate which of the modes may be more affected by the collisional effects in the simulations. The left panel of Fig. 11 shows the vertical eigenfunctions, v_z, obtained semi-analytically for the perpendicular field case for several values of n with the original values of the density contrast and the magnetic field. The eigenfunctions are normalized to unity. They correspond to the fastest growing modes for the wave numbers n indicated in the legend, and we can observe that the high n eigenmodes are more spatially localized in the vertical direction. The right panel of Fig. 11 quantifies the vertical extent of the eigenmodes by plotting half width at half maximum, Δ, of the fastest growing modes as a function of n for different configurations of the magnetic field corresponding to the growth rates shown in the left panel of Fig. 9. A linear fit of log(Δ) vs. log(n) is made for high-n modes for each of the magnetic field configurations.

Fig. 11. Left: eigenfunctions corresponding to the largest growth rates for the modes n = 1 (black line), n = 2 (blue line), n = 16 (green line), and n = 128 (red line), normalized to unity, as a function of height for the perpendicular magnetic field configuration. Right: Log–log plot of half width at half maximum, Δ, of the eigenfunctions as a function of the mode number n for the range 1 ≤ n ≤ 50 for different configurations of the magnetic field: perpendicular (orange), sheared with shear length equal to: L0 (red), L0/2 (violet), L0/3 (cyan), L0/4 (lime), L0/5 (blue). The slopes of the curves calculated from a linear fit for the modes 10 ≤ n ≤ 50 and the uncertainties of the fit are indicated in the legend.

We observe from the right panel of Fig. 11 that there is a clear distinction in the Δ(n) functional form for the case of purely perpendicular magnetic field from those with sheared field. The configuration with perpendicular magnetic field has Δ ∝ n^−1/2, while all of the RTI-unstable sheared configurations consistently show Δ ∝ n⁻¹. The wavelength dependence of the vertical extent of the RTI modes in the case with no in-plane field and a vertically limited RTI-unstable region with (dρ₀/dz)> 0 was discussed in Sect. 4.1.2 and the Δ ∝ n^−1/2 dependence is consistent with the conclusions derived from Eq. (24). In turn, the Δ ∝ n⁻¹ dependence for the sheared field configurations is as estimated in Eq. (29), with the eigenfunctions becoming progressively more localized for smaller shear lengths. This implies that high-n modes are more prone to collisional dissipation effects; the effect should be substantially more pronounced for the sheared field configurations, which is consistent with the temporal behavior of different n modes observed in Fig. 4.

6. Decoupling

Here, we consider the decoupling effects as observed in the simulations. To do so, we quantify the differences in the horizontal and vertical velocity components between the neutrals and charged particles. We recall that the driver for the RTI instability is in the gravitational free energy of a prominence thread with a high density of neutrals; that is, the instability acts on the neutral fluid. At the same time, electromagnetic forces act only on charged particles, so it can be expected that incomplete collisional coupling may result in some differences in the velocities of the charged and neutral species. The viscous force is another important agent since in the plasmas under consideration, the viscosity of neutrals is much greater than viscosity of charged particles. The decoupling can occur when the differentiating viscous or electromagnetic forces (or both) become larger than the collisional coupling between neutrals and charged particles.

The normalized decoupling in the vertical velocity (and equivalent for the horizontal velocity) is computed according to the following steps: (i) the Fourier amplitude of the difference between the vertical velocities of charged particles and neutrals, v_zc − v_zn, is computed for all modes n as a function of height; (ii) the Fourier amplitude of the vertical velocity of charged particles, v_zc, is computed for all modes n as a function of height; (iii) the above quantities are averaged separately for the heights between −L₀ and L₀ and over four consecutive modes; (iv) the normalized decoupling is obtained as the natural logarithm of the ratio between the two quantities.

Figure 12 compares the vertical normalized decoupling between flows of neutrals and charged particles for the simulations L1-WN (top left), L1-WN-NN (top right), L1-WN-NN-B (bottom left), and P-WN-NN (bottom right) for the first 48 modes. Similarly to Fig. 5, we choose the averaging over four consecutive modes and the upper limit of n = 48 for clarity of the presentation.

Fig. 12. Normalized vertical decoupling in vertical velocity for different harmonics as a function of time. Different color lines correspond to different harmonics grouped in sets of four, from n = 1 to n = 48. Top left: L1-WN simulation; top right: L1-WN-NN; bottom left: L1-WN-NN-B; bottom right: P-WN-NN.

In all four cases considered, the normalized decoupling for all n has an initially high value that rapidly drops early in the simulations. The initial magnitude of the decoupling is due to the instability being seeded by perturbations in the density of neutrals, without a corresponding perturbation in the density of charged particles. This leads to the generation of neutral flows that are then rapidly followed by corresponding flows in the charged fluid. We can observe that the amount of decoupling strongly varies with n for the sheared magnetic field configurations, while for the perpendicular magnetic field case all modes behave nearly identically during the linear phase. This is explained by the effect of magnetic tension on the development of RTI in the “L1” sheared magnetic field cases, which is absent in the “P” case. Magnetic tension acts only on charged particles, thus seeding decoupling, and more so on the low n modes that are less localized and therefore sample regions of stronger in-plane magnetic field. Thus, the L1-WN-NN simulation (top right) shows larger decoupling than the P-WN-NN one (bottom right). In a similar vein, since the Lorentz force is larger for larger magnetic field, the decoupling in the L1-WN-NN-B case (bottom left) is generally greater than in the L1-WN-NN case (top right).

In the nonlinear phase, in all cases, the decoupling increases with the mode number n. Such behavior is expected because the decoupling should occur for scales below the ion-neutral (or neutral-ion) collisional scales as shown in Fig. 8. The increase of the density contrast produces lower decoupling (L1-WN vs L1-WN-NN case). This is also in agreement with the analysis of the ion-neutral collisional scales from Fig. 8, as n_νin for the increased density case is much larger than for the original density case. Moreover, the viscosity is smaller for the increased density contrast case, therefore neutrals are less affected, and the decoupling can be expected to be less pronounced.

The above conclusions can also be confirmed by a visual inspection of the snapshots of the decoupling shown in Fig. 13. The locations where the decoupling increases, both in horizontal and vertical velocity, locally coincide with locations where the magnetic field gets compressed (corresponding to higher density of the magnetic field potential contours). This greater decoupling is produced as a result of different movements on the part of neutrals and charged particles across the field lines. The decoupling in the vertical velocity is preferentially negative in a larger area right below the instability eddies. A negative value of the decoupling implies that the downward velocities of neutrals are larger than the downward velocities of charged particles. This behavior is consistent with the fact that neutrals drive the instability as they cannot be sustained by the magnetic forces other than through collisions with charged particles. Nevertheless, the absolute value of the normalized decoupling in a horizontal velocity is generally larger (except for the case of the enhanced density L1-WN-NN, top right). The simulation L1-WN shows the largest decoupling among the four cases, reaching values as high as 5−10% of the corresponding velocity. When the density contrast is increased (top right), the values of the vertical decoupling become larger than the values of the horizontal decoupling. On the contrary, when the magnetic field strength is increased (bottom left), the horizontal decoupling relative to the horizontal velocity becomes larger than the vertical decoupling relative to the vertical velocity.

Fig. 13. Pairs of snapshots of the horizontal and vertical decoupling in velocities for the simulations P-WN (top left), P-WN-NN (top right), L1-WN-NN-B (bottom left), and P-WN-NN (bottom right). The snapshots are taken at the time indicated at the right-hand-side image of each pair. The maximum absolute values of the corresponding decoupling and velocity are indicated at the bottom of each image. The black lines represent the iso-contours of the absolute value of the y-component of the magnetic field potential, |Ay|. The iso-contours are equally spaced between 0.6 and , where is the maximum value of |Ay|.

7. Magnetic structures

To get further insights on the formation of magnetic structures, Fig. 14 shows the power spectra of magnetic energy for each of the components of the magnetic field calculated at the last available time for the four simulations L1-WN, L1-WN-NN, L1-WN-NN-B, and P-WN-NN. The power spectra have been constructed as follows: (i) different components of the magnetic energy are computed as functions of x and z; (ii) Fourier transform of these components in x direction is done to obtain their Fourier amplitudes as a function of the vertical direction; (iii) the Fourier amplitudes are averaged at heights between z = −L₀ and z = L₀.

Fig. 14. Power spectra of magnetic energy for each of the components of the magnetic field calculated at the final time for the simulations L1-WN (top left), L1-WN-NN (top right), and L1-WN-NN-B (bottom left), and P-WN-NN (bottom right). We note that the simulation P-WN-NN has no magnetic energy in the x and z directions. Solid lines are individual components of the magnetic energy, (blue), (red), (orange). Dashed lines are for the in-plane component, (green), and the total energy, (black).

Figure 14 shows qualitatively similar behavior of the magnetic energies for the three L1 simulations. In particular, in all three cases and for all modes n, the mode amplitude of total in-plane magnetic energy (, dashed green lines) is greater than in each of the in-plane components of the magnetic energy ( and , blue and orange solid lines). This implies that the in-plane components of magnetic field energy vary approximately when they are in-phase with each other.

On the other hand, the mode amplitudes of the total in-plane magnetic field energy and the energy in the out-of-plane B-field component (, red solid lines) appear to be additive at low n, but not at higher n. At low n, the mode amplitude of the total magnetic energy is approximately a sum of the corresponding amplitudes of its components. The value of n where this behaviour breaks is different for the three L1 cases, with the largest shown for the enhanced density case (upper right panel). For the higher n modes, the mode amplitudes of the in-plane and out-of-plane B-field energy components converge to nearly the same values and both greatly exceed the mode amplitudes of the total magnetic energy; that is, black dashed lines are below the green dashed and the red solid lines, which begin to overlap. This indicates that the high-n modes of the in-plane and the out-of-plane components of magnetic field energy have the same amplitude and are nearly exactly out of phase.

The scale where the total in-plane, , and out-of-plane, , components transition from being additive to canceling each other seems to be related to the largest size of the in-plane plasmoids seen in Fig. 4 for the same time moment. We can observe that the structures formed in the L1-WN case are the largest between three L1 simulations, followed by the sizes in the L1-WN-NN-B case, and then by the L1-WN-NN case. The ordering between these scales is the same as the ordering of the scales where the out-of-plane and the in-plane components of magnetic energy in Fig. 14 become out of phase. Some authors found breaks in the power spectra at a certain scale (Freed et al. 2016; Hillier et al. 2017; Leonardis et al. 2012), which might be related to the formation of the magnetic structures, which will be discussed in a forthcoming paper. The 3D configuration of the magnetic field, with a dominant component in the y-direction being out of phase with the in-plane component, may resemble a flux rope.

8. Discussion and conclusions

The main goal of our paper is to the study of the effects of partial ionization on the development of the RTI at a realistically smooth transition layer between a magnetized prominence thread and the solar corona. We use a two-fluid model, which evolves neutral and charged species, coupled by collisional terms, where we consider both the elastic and inelastic processes. We explore the dependence of the linear and nonlinear evolution of the instability on the magnetic field configuration and mass loading. Our main conclusions can be summarized as follows:

(i) The consideration of a smooth RTI-unstable transition layer as opposed to the typically considered discontinuous interface is both necessary to self-consistently evaluate two-fluid effects on the RTI instability and to demonstrates a different set of RTI characteristics than previously observed in discontinuous interface studies.

(ii) For sheared magnetic field equilibria, the in-plane component of the magnetic field suppresses RTI growth in the linear phase. Both in simulations and in analytical calculations using single-fluid ideal MHD treatment, the growth rate at all scales is smaller when the magnetic field is sheared compared to the perpendicular configuration. The linear growth rate decreases with decreasing shear length, stabilizing the RTI for sufficiently small shear length for wavelength shorter than the density gradient scale; however, sheared equilibria that are unstable to wavelengths of order of the density gradient scale are also unstable with regard to modes of arbitrarily short wavelengths with the same growth rate.

(iii) For a smooth transition layer, low wave-number RTI eigenfunctions have vertical extents that encompass the prominence thread, while higher n modes become localized at the location of highest density gradient. The presence of sheared magnetic field provides further and more rapid vertical localization of the eigenfunctions as a function of n, making high-n modes even more prone to collisional damping effects.

(iv) When mass loading and the density contrast is increased, both semi-analytical ideal MHD calculations and simulations show that the linear growth rate is larger at all scales. Additionally, greater mass loading leads to stronger ion-neutral collisions and lower viscosity. Thus, smaller scales develop in the linear stage of the high density contrast simulations, while the relative decoupling between charged particles and neutrals observed during the nonlinear phase is also smaller.

(v) Ion-neutral collisions and, to a smaller extent, neutral viscosity effects appear to suppress growth of short wavelength modes in the linear phase of the RTI instability. This is established by estimating the scales expected to be impacted by the collisional effects in Fig. 8 and then comparing the linear growth rates calculated from the simulations to the linear growth rates obtained semi-analytically under the single-fluid ideal MHD approximation. These effects will be further explored in follow-up work.

(vi) The increase of the magnetic field magnitude decreases RTI growth rate. At the same time, in the nonlinear phase, the flow decoupling between charged particles and neutrals appears to be strongest near the locations of highest magnetic field stress. We thus infer that much of the decoupling observed in the simulations is likely driven by the magnetic activity due to charged particles, but not neutrals, as the former are directly impacted by magnetic field forces.

(vii) In simulations with the sheared magnetic field equilibria, the spatial spectra of the out-of-plane and in-plane components of the magnetic energy show strong out-of-phase coherence at small scales, but not at the largest scales, during the nonlinear phase of the evolution. The transition scale is different in simulations with different mass loading and magnetic field strength, with the L1-WN simulation showing both the strongest magnetic activity and highest degree of such magnetic self-organization across the scales. These effect will be further explored in follow-on work.

It is known that the linear growth rate of RTI is related to the properties of the background atmosphere, such as the density contrast and the magnetic field. In prominences, the growth rate and the degree of mass loading can be measured from observations and such measurements have been used in the past to derive the magnetic field properties based on observations of the RTI. Here, we obtain that other measurable quantities, such as the decoupling in velocity of ions and neutrals, are also related to the background variables. These latter quantities are much more difficult to measure because they require high resolution and high temporal cadence of the data. Such data could possibly be obtained using the large-aperture 4 m solar telescope facilities (DKIST¹, EST²).

In all the simulations performed, we found that ion-neutral interactions and, to a smaller extent, viscosity damp linear RTI development at smallest scales. The comparison of the numerical linear growth rate to the linear growth rates obtained semi-analytically in the single-fluid ideal MHD approximation, neglecting the effect of the collisions, confirmed that the suppression of the small scales is indeed produced by dissipation effects.

In our simulations, the largest scales do not show effects of dissipation due to ion-neutral collisions and viscosity, with the global RTI growth rate remaining the same as in the single-fluid ideal MHD regime. Some of the previous numerical simulations have similarly shown that the characteristics of RTI flows are identical for the viscous and inviscid models (Doludenko et al. 2019). On the other hand, Mitra et al. (2016) and Gerashchenko & Livescu (2016) found an increased growth rate at all wave numbers in the inviscid case, compared to the case with viscosity. The apparent contradiction between the two results may come from the different scales considered in the simulations since the dissipative effects occur only for the scales below a certain scale.

We find that the linear growth rate increases with increasing the density contrast. This result is in agreement with the classical hydrodynamic model from Chandrasekhar (1961) for the discontinuous density profile, as well as for the density profile exponentially stratified with height (Livescu 2004). It has been shown that compressibility causes a destabilizing effect, which becomes more important for smaller density contrasts (Livescu 2004; Ribeyre et al. 2005). In the nonlinear phase, it has been shown that increasing the density contrast results in a higher growth rate, thinner fingers, and less roll-up at the tip of the fingers due to secondary KHI (Youngs 1984; Gardner et al. 1988; Jun et al. 1995). The magnetic field affects the growth rate at small scales, but even in the magnetic RTI, the growth rate is still larger for a larger density contrast (Dolai et al. 2016).

In our case, the viscosity coefficient is smaller in the model with enhanced density contrast, and the ion-neutral collisional frequency, is higher. This reduction of viscosity and the importance of ion-neutral interaction act in favor of increasing of the growth rate, together with the proper density contrast enhancement.

We find that the viscous force acting dominantly on neutrals along with the electromagnetic force acting only on the charged particles produce the observed decoupling of the respective flows in the nonlinear phase of RTI development. This decoupling decreases with decreasing importance of these forces. For example, when the density contrast is enhanced, the magnitude of the Lorentz force relative to the gravitational force decreases and the decoupling also decreases. Similarly, when the density contrast is enhanced, the neutral viscosity becomes smaller, leading to lower values for the decoupling. Overall, we observe in our simulations that the decoupling can reach up to 5% of the individual flow velocities. There is a correlation between the locations with larger decoupling and the locations of greatest magnetic field stress (i.e., the locations where the Lorentz force is increased).

The magnetic field parallel to the perturbation plane has a stabilizing effect on the RTI. Early analytical studies by Chandrasekhar (1961) have shown that a component of the magnetic field parallel to the direction of the perturbation suppresses the small-scale perturbations with wavelengths smaller than a cutoff wavelength, proportional to the square of the in-plane field strength, as seen from Eq. (22). In the nonlinear phase, Carlyle & Hillier (2017) measured the growth of the magnetic RTI by means of numerical simulations in a 3D setup. These authors found that the growth rate decreases with increasing magnetic field and they attributed it to the effects of higher magnetic tension, requiring greater energy for the frozen-in plasma to move. In our simulations, we found that the linear growth rate is smaller when the field is sheared compared to the perpendicular field case. The ideal MHD semi-analytical calculations for the sheared magnetic field configurations similarly showed that the peak growth rate decreases with decreasing shear length L_s; furthermore, for sufficiently small shear lengths, all modes with a wavelength shorter than the density gradient scale are stabilized due to a combination of magnetic tension and in-plane magnetic pressure effects. However, both analytically estimated and semi-analytical calculations have also shown that for a range of shear lengths, in the absence of collisional effects, the RTI growth rate becomes wavelength-independent in the short wavelength limit.

Using a single-fluid approximation with two-fluid effects introduced through the ambipolar term, Díaz et al. (2014) found that collisions may remove the stabilizing effect of the magnetic field on the small scales. Thus, the linear growth rate at small scales would become larger when the collisions are taken into account. In our case, the small scales are suppressed when the collisions are taken into account through neutral-charge interaction and viscosity. This is due, primarily, to our consideration of a smooth rather than a discontinuous prominence-corona interface with a gradually rotating magnetic field, such that the scales associated with magnetic field stabilization are much smaller than those stabilized by collisional dissipation, as shown in Fig. 9.

At small scales, suppressed by dissipation effects in the linear phase, are introduced in the early nonlinear stage by energy cascade. In the nonlinear phase, we observe a continued production of magnetic energy at intermediate and small scales due to the twisting of the field lines. The in-plane plasmoids observed in the snapshots, together with the dominant perpendicular component of the magnetic field, create 3D structures resembling flux ropes. This is confirmed by the Fourier analysis of the three components of the magnetic energy. The strong out-of-phase correlation between the in-plane and the out-of-plane magnetic energy for intermediate and high n modes hints at magnetic self-organization processes taking place during nonlinear RTI evolution. The scales where the two components of the magnetic energy become coherent seem to be related to the largest scales of the in-plane magnetic structures. Further study of this magnetic field dynamics will be presented in future work.

Acknowledgments

This work was supported by the Spanish Ministry of Science through the project PGC2018-095832-B-I00 and the US National Science Foudation. It contributes to the deliverable identified in FP7 European Research Council grant agreement ERC-2017-CoG771310-PI2FA for the project “Partial Ionization: Two-fluid Approach”. The authors wish to acknowledge the contribution of Teide High-Performance Computing facilities and LaPalma Supercomputer to the results of this research. TeideHPC facilities are provided by the Instituto Tecnológico y de Energías Renovables (ITER, SA). URL: http://teidehpc.iter.es

References

All Tables

All Figures

Fig. 1. Sketch of the initial configuration. The simulation domain (filled with blue color) is contained in the XOZ plane, where x varies between −L0 and L0 and z varies between −4L0 and 4L0. The magnetic field, B (red lines), is contained in the XOY plane (filled with gray color). At the height z = 0, the magnetic field only has the component, By, and it is slowly rotated above and below z = 0 by an angle which depends on height, keeping its modulus constant. The value of the angle is anti-symmetric with respect to z = 0. The maximum angle obtained at the top and the bottom of the atmosphere is indicated by θ0. When the field is perpendicular, the direction of B is along the y-axis for all the heights, θ0 = 0°. For the sheared magnetic field configurations θ0 = 1°. The gravity, g, points in the negative z direction and is indicated by a green line. The direction of the perturbation, k, is shown by a blue line.

In the text

Fig. 2. Parameters of the initial model atmosphere as a function of the vertical coordinate, z. Top left: density of neutrals in the models with the original (red) and enhanced (dashed green) density contrasts, and density of charged particles (blue). Top right: same for the pressure of neutrals and charged particles. Bottom left: temperature of both the charged particles and neutrals. Bottom right: modulus of the magnetic field in the models with the original (red) and enhanced (green) density contrasts.

In the text

Fig. 3. Left: coefficients of the viscosity of charged particles (, blue line) and neutrals (, orange lines), and the thermal conductivity of neutrals (, green lines), in units of m2 s−1, as from Eq. (14). The neutral viscosity and thermal conductivity coefficients are limited to 109 m2 s−1. Right: ion-neutral (νin, blue lines), neutral-ion (νni, orange line) collision frequencies, and frequencies associated with waves: fast magneto-acoustic (ωfast, green lines), in-plane Alfvén for “L1” shear (, red lines), acoustic in the neutral fluid (ωcn, black line) as functions of height. Values for the model with the original density contrast model are shown with solid lines, and for the enhanced density contrast with dotted lines.

In the text

Fig. 4. Time series of the snapshots of neutral and charged densities for the simulations L1-WN (top left group of panels), L1-WN-NN (top right group of panels), L1-WN-NN-B (bottom left group of panels), P-WN-NN (bottom right group of panels). The in plane (x and z components) velocities of neutrals and charged particles are represented by orange arrows in the plots which show the densities of neutrals and charged particles, respectively. For the simulations with sheared magnetic field (L1-WN, L1-WN-NN, L1-WN-NN-B), the iso-contours of the absolute value of the y-component of the magnetic field potential, |Ay|, are shown at the bottom panels. The iso-contours are equally spaced between 0.6 and , where is the maximum value of |Ay|.

In the text

Fig. 5. Temporal evolution of individual harmonics, computed from the vertical velocity of neutrals, as a function of time. Different color lines show the Fourier amplitudes of the harmonics, grouped in sets of four, from n = 1 to n = 48. The amplitudes are averaged between heights −L0 and L0. Different panels are for different simulations: L1-WN (top left); L1-WN-NN (top right); L1-WM-NN-B (bottom left); P-WN-NN (bottom right). The vertical dashed pink lines mark the times of the snapshots shown in Fig. 4.

In the text

Fig. 6. Cutoff wavelengths (black lines) due to the in-plane component of the magnetic field as a function of height, computed from Eq. (22) for the equilibrium profile with the values of nn00 = 1016 m−3 and B00 = 10−3 T in Eqs. (5) and (7). Cutoff wavelengths for the following four shear lengths: Ls = [0, L0/10, L0/2, L0] are shown with increasing width of the line for increasing Ls. The red dotted line corresponds to the density gradient scale length Ld.

In the text

	Fig. 7. Growth rate obtained from the exact solution for the exponential density profile, Eq. (3) with uniform β = 1.83 × 10−6 m−1 (blue curve marked with “o”) and β = 2.44 × 10−6 m−1 (orange curve marked with “x”). The values of n are between n = 1/64 and n = 50.
In the text

Fig. 8. Modes corresponding to the neutral viscosity (blue lines), ion-neutral collision scale (orange lines), and neutral-ion collision scale (green lines) as a function of height, computed for the original density profile (solid lines), for the enhanced density profile (dashed lines), and for the enhanced density and magnetic field profiles (dotted lines), from Eqs. (33)–(34). The horizontal dashed line marks the maximum mode number shown in Fig. 5, n = 48, and the dotted line represents the largest mode number that can be resolved in the domain, using two points, n = 256.

In the text

Fig. 9. Growth rates as a function of the mode number, n, in logarithmic scale up to n = 50. The first value is n = 1/64 for the left panel and n = 1 for the right panel. Left: growth rate obtained semi-analytically for the perpendicular magnetic field (orange); sheared field with Ls equal to L0 (red), L0/2 (violet), L0/3 (cyan), L0/4 (lime), L0/5 (blue), L0/10 (dark green); and for uniformly rotated magnetic field by 1° (black). Right: growth rate obtained from the simulations P-MM-S (blue), P-WN-S (orange), L1-MM-S (green), L1-WN-S (red).

In the text

Fig. 10. Comparison between the growth rates obtained semi-analytically in the incompressible MHD approximation (left) and the growth rates obtained from simulations (right) and as functions of the mode number, n, in logarithmic scale, limited to n = 50. In both panels, the perpendicular magnetic field cases are shown with solid lines and sheared case are shown with dashed lines. Orange lines: P-WN-S; black lines: P-WN-NN-S; red lines: L1-WN-S; green lines: L1-WN-NN-S, and blue lines for L1-WN-NN-B-S.

In the text

Fig. 11. Left: eigenfunctions corresponding to the largest growth rates for the modes n = 1 (black line), n = 2 (blue line), n = 16 (green line), and n = 128 (red line), normalized to unity, as a function of height for the perpendicular magnetic field configuration. Right: Log–log plot of half width at half maximum, Δ, of the eigenfunctions as a function of the mode number n for the range 1 ≤ n ≤ 50 for different configurations of the magnetic field: perpendicular (orange), sheared with shear length equal to: L0 (red), L0/2 (violet), L0/3 (cyan), L0/4 (lime), L0/5 (blue). The slopes of the curves calculated from a linear fit for the modes 10 ≤ n ≤ 50 and the uncertainties of the fit are indicated in the legend.

In the text

	Fig. 12. Normalized vertical decoupling in vertical velocity for different harmonics as a function of time. Different color lines correspond to different harmonics grouped in sets of four, from n = 1 to n = 48. Top left: L1-WN simulation; top right: L1-WN-NN; bottom left: L1-WN-NN-B; bottom right: P-WN-NN.
In the text

Fig. 13. Pairs of snapshots of the horizontal and vertical decoupling in velocities for the simulations P-WN (top left), P-WN-NN (top right), L1-WN-NN-B (bottom left), and P-WN-NN (bottom right). The snapshots are taken at the time indicated at the right-hand-side image of each pair. The maximum absolute values of the corresponding decoupling and velocity are indicated at the bottom of each image. The black lines represent the iso-contours of the absolute value of the y-component of the magnetic field potential, |Ay|. The iso-contours are equally spaced between 0.6 and , where is the maximum value of |Ay|.

In the text

Fig. 14. Power spectra of magnetic energy for each of the components of the magnetic field calculated at the final time for the simulations L1-WN (top left), L1-WN-NN (top right), and L1-WN-NN-B (bottom left), and P-WN-NN (bottom right). We note that the simulation P-WN-NN has no magnetic energy in the x and z directions. Solid lines are individual components of the magnetic energy, (blue), (red), (orange). Dashed lines are for the in-plane component, (green), and the total energy, (black).

In the text