Issue 
A&A
Volume 646, February 2021



Article Number  A93  
Number of page(s)  19  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202039053  
Published online  12 February 2021 
Twofluid simulations of RayleighTaylor instability in a magnetized solar prominence thread
I. Effects of prominence magnetization and mass loading
^{1}
Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain
email: beatriceannemone.popescubraileanu@kuleuven.be
^{2}
Departamento de Astrofísica, Universidad de La Laguna, 38205 La Laguna, Tenerife, Spain
^{3}
National Science Foundation, Alexandria, VA 22306, USA
Received:
28
July
2020
Accepted:
8
December
2020
Solar prominences are formed by partially ionized plasma with interparticle collision frequencies, which generally warrant magnetohydrodynamic treatment. In this work, we explore the dynamical impacts and observable signatures of twofluid effects in the parameter regimes when ionneutral collisions do not fully couple the neutral and charged fluids. We performed 2.5D twofluid (charge – neutrals) simulations of the RayleighTaylor instability (RTI) at a smoothly changing interface between a solar prominence thread and the corona. The purpose of this study is to deepen our understanding of the RTI and the effects of partial ionization on the development of the RTI using nonlinear twofluid numerical simulations. Our twofluid model takes into account viscosity, thermal conductivity, and collisional interaction between neutrals and charge: ionization or recombination, energy and momentum transfer, and frictional heating. In this paper, we explore the sensitivity of the RTI dynamics to the prominence equilibrium configuration, including the impact of the magnetic field strength and shear supporting the prominence thread, and the amount of prominence massloading. We show that at small scales, a realistically smooth prominencecorona interface leads to qualitatively different linear RTI evolution than that which is expected for a discontinuous interface, while magnetic field shear has the stabilizing effect of reducing the growth rate or eliminating the instability. In the nonlinear phase, we observe that in the presence of field shear the development of the instability leads to formation of coherent and interacting 2.5D magnetic structures, which, in turn, can lead to substantial plasma flow across magnetic field lines and associated decoupling of the fluid velocities of charged particles and neutrals.
Key words: Sun: atmosphere / instabilities
© ESO 2021
1. Introduction
The RayleighTaylor 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 magnetoinertial 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. Smallscale 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^{−1}. 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 twofluid 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 chargeneutral 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, modecoupling 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, ρ_{2}, which is initially on top of a lighter fluid with mass density, ρ_{1}, 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:
where,
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 −L_{0} and L_{0} and z varies between −4L_{0} and 4L_{0}. 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, B_{y}, 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 antisymmetric 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 yaxis 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 twofluid 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 twofluid calculations.
The general case of RTI with arbitrary nonuniform 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:
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 −ω^{2} = 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 kdependence 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 selfsimilar 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 KelvinHelmholtz 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 twofluid 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 cutoff imposed by the magnetic field.
Most of the numerical simulations of RTI have been done using the singlefluid 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 prominencecorona transition region (PCTR).
Numerical studies of RTI using the twofluid 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 indepth analysis of the results. Recently, Hillier (2019) studied KelvinHelmholtz instability numerically in the twofluid 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 ionizationrecombination.
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 twofluid 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 twofluid 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 singlefluid analysis is necessary to begin to understand the more complex twofluid 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 twofluid 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 singlefluid incompressible MHD approximation for smoothly varying density and magnetic field profiles in Sect. 4. We then study the linear behavior of RTI using semianalytical calculations for the assumed background atmosphere and in the twofluid 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 yaxis. For this configuration, θ_{0} = 0°, as in the sketch in Fig. 1.
2.1. Background atmosphere
The equilibrium atmosphere is uniform in the xdirection. 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:
where the ion number density at z = 0 is n_{0} = 10^{15} m^{−3}; the peak neutral number density, reached at z = L_{0}/2 is n_{n00} = 10^{16} m^{−3}; the background temperature of the corona T_{b} = 2.02 × 10^{5} K; the neutral number density corresponding to the corona temperature (T_{b}) is n_{nb} = 3.5 × 10^{9} m^{−3}; B_{00} = 10^{−3} T is the gravitational scale height of the charged particles of H_{c} = 2k_{B}T_{b}/(m_{H}g). The characteristic length scale is L_{0} = 1 Mm. The plasma is calculated at z = 0 using B_{00} as the value of the magnetic field and has the value of β_{p} ≈ 1.4 × 10^{−2}. The function f is defined as follows:
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_{0}/2, where it attains the maximum, n_{n00}, as from Eq. (5). The minimum of the temperature is also attained at z = L_{0}/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 magnetohydrostatic equation where the total density and pressure are used. The integration constants permits the scaling of the magnetic field around some chosen value, B_{00} 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:
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_{x}x). 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 xdirection is generated as a multimode perturbation (“MM”) specified by:
where m_{x} is the number of points in the x direction and ϕ_{n} is a random phase.
2.3. Twofluid evolution equations
We solve the twofluid 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 twofluid version of the MANCHA3D code, MANCHA3D2F. The singlefluid version of this code is described in Felipe et al. (2010), GonzálezMorales et al. (2018), and the twofluid version preserves several numerical features of MANCHA3D, such as hyperdiffusive 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 chargeexchange 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:
Simulations used in the paper.
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
These simulations are referred in the table with the keyword “P;”
the sheared magnetic field
We use a fixed value of θ_{0} = 1° and L_{s} = L_{0} 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_{0}/2 and L_{0}/2, that is, ρ_{2} = ρ_{0}(z = L_{0}/2) and ρ_{1} = ρ_{0}(z = −L_{0}/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_{00} 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^{−2} 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^{−4}, 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_{0}, L_{z} = 8L_{0}. 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_{0} to L_{0}, and z varies from −4L_{0} to 4L_{0}.
In all our simulations, we avoid using hyperdiffusion 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 neutralcharge 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 neutralion (ν_{ni}) and ionneutral (ν_{in}) collision frequencies. The latter two are also compared to the timescales associated with the fast magnetoacoustic wave (ω_{fast}), the inplane Alfvén wave (), and the neutrals sound wave (ω_{cn}) within the background atmospheres as characteristic of dynamical timescales in the system. These are calculated as follows,
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 timestep imposed by such large coefficients of neutral viscosity and thermal conductivity, the values of these coefficients have been limited to 10^{9} m^{2} s^{−1} 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 m^{2} s^{−1}, as from Eq. (14). The neutral viscosity and thermal conductivity coefficients are limited to 10^{9} m^{2} s^{−1}. Right: ionneutral (ν_{in}, blue lines), neutralion (ν_{ni}, orange line) collision frequencies, and frequencies associated with waves: fast magnetoacoustic (ω_{fast}, green lines), inplane 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 ionneutral, ν_{in}, and the neutralion, ν_{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 xdirection for the fast magnetoacoustic wave, ω_{fast}, the inplane 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 largescale ideal MHD frequencies, and, thus, largescale 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 twofluid effects.
3. Global overview of simulation dynamics
Figure 4 shows time series of snapshots of the four simulations, L1WN, L1WNNN, L1WNNNB, and PWNNN. 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 (L1WN, L1WNNN, L1WNNNB) 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 inplane magnetic field lines plotted as isocontours of the absolute value of the ycomponent of the magnetic potential obtained by numerical integration. The isocontour 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 L1WN (top left group of panels), L1WNNN (top right group of panels), L1WNNNB (bottom left group of panels), PWNNN (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 (L1WN, L1WNNN, L1WNNNB), the isocontours of the absolute value of the ycomponent of the magnetic field potential, A_{y}, are shown at the bottom panels. The isocontours are equally spaced between 0.6 and , where is the maximum value of A_{y}. 
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 −L_{0} and L_{0}. Different panels are for different simulations: L1WN (top left); L1WNNN (top right); L1WMNNB (bottom left); PWNNN (bottom right). The vertical dashed pink lines mark the times of the snapshots shown in Fig. 4. 
By comparing the snapshots of L1WN to the snapshots of L1WNNN, we can observe that smallscale 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_{0} and z = L_{0}, 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 topleft (L1WN) and topright panels (L1WNNN) we observe that the increase in the density contrast increases the growth rate of all modes. By comparing topright (L1WNNN) and bottomright (PWNNN) panels, we observe that the magnetic field shear has a stabilizing effect; and by comparing topright (L1WNNN) and bottomleft panel (L1WNNNB), 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 ionneutral collisions. The largest mode number n that grows in the linear phase in the enhanced density case, L1WNNN (top right), is n = 33, while for the original density case, L1WN (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 singlefluid 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 zerodivergence condition on velocity for the energy equation, the equations become:
We consider a 2.5D case, where there are no variations in the ydirection, 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:
with the general solution in the linear regime being a superposition of such singlemode solutions. Linearizing Eq. (15) and substituting perturbations of the form defined in Eq. (16), we have:
After manipulating the above we obtain a second order differential equation for v_{z}:
where a is defined as:
where ρ_{0} is the total density, calculated from our twofluid model as:
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_{0})=v_{z}(4L_{0})=0. In doing so, we look for the fastest growing mode for a given k_{x}. We refer to this method as semianalytical 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 semianalytically 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 nonzero. However, for modes with k_{x}L_{d}/2π ≲ 1, we can approximate our vertical profile as a discontinuous profile, using values of ρ_{2} = ρ_{0}(z = L_{0}/2)≈1.83 × 10^{−11} kg m^{−3}, and ρ_{1} = ρ_{0}(z = −L_{0}/2)≈2.93 × 10^{−12} kg m^{−3}, with ρ_{1} and ρ_{2} 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:
The condition for stability −ω^{2} < 0 gives the expression for the critical wavelength λ_{c}, so that the scales with wavelength λ = 2π/k_{x} < λ_{c} do not grow:
where is the characteristic inplane 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 ρ_{1} and ρ_{2} specified above for different magnetic field shear length scales, L_{s} = [L_{0}/10, L_{0}/2, L_{0}]. 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_{0}/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_{0}/10 and smaller.
Fig. 6. Cutoff wavelengths (black lines) due to the inplane component of the magnetic field as a function of height, computed from Eq. (22) for the equilibrium profile with the values of n_{n00} = 10^{16} m^{−3} and B_{00} = 10^{−3} T in Eqs. (5) and (7). Cutoff wavelengths for the following four shear lengths: L_{s} = [0, L_{0}/10, L_{0}/2, L_{0}] are shown with increasing width of the line for increasing L_{s}. The red dotted line corresponds to the density gradient scale length L_{d}. 
4.1.2. Exponential density profile
Another limiting case is the case of an exponential density profile, ρ_{0} ∝ 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:
Using the same heavy and light fluid density values as in the subsection above, ρ_{2} = 1.83 × 10^{−11} kg m^{−3} and ρ_{1} = 2.93 × 10^{−12} kg m^{−3}; we match them to an exponential density profile ρ_{0} ∝ exp(βz), for z between −L_{0}/2 and L_{0}/2, with constant β. We obtain β = ln(ρ_{2}/ρ_{1})/L_{0} = 1.83 × 10^{−6} m^{−1}. Alternatively, we can calculate β from the local density gradient scale length at z = 0 and we obtain a value of β = 2.44 × 10^{−6} m^{−1}. Figure 7 shows the dispersion relations obtained by substituting the two values of β and d = 8L_{0}, 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 semianalytical calculations presented below in Sect. 5. On small scales, for β = 2.44 × 10^{−6} m^{−1}, the value of the growth rate in the high n limit is the Brunt–Väisälä frequency ≈2.58 × 10^{−2} s^{−1}.
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 RTIunstable 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 nonnegligible 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:
and the most unstable eigenmodes v_{z}(z) will be those with v_{z} ≈ 0 in the regions of stable stratification, (dρ_{0}/dz)≤0, corresponding to a maximal value of the righthandside integral. It follows that for ρ_{0} as given by Eq. (20), we can expect the fastest growing RTI modes to have v_{z} ≈ 0 for z≳L_{0}/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 fieldline tension counteracting the gravity force driving the instability. Thus, in magnetized 3D systems with unidirectional Bfield, 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_{x}x̂, 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:
where b_{0} is constant and L_{s} is the characteristic length of magnetic field shear. We assume the density profile to be exponentially stratified, ρ_{0} = ρ_{00}exp(βz), as in Sect. 4.1.2 above.
For sufficiently large b_{0}, 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 zmomentum equation (Eq. (17)) and looking for z = z_{d} where:
It follows that
and it is easy to show that positive solutions of Eq. (27) for z_{d} only exist when
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_{0} modes with k_{x}L_{0} < 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:
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_{0} to give k_{x}z_{d} ≈ 4.91, with proportionally smaller values for smaller L_{s}.
We now substitute the expression for from Eq. (29) in place of d^{2} 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:
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:
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 inplane 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 inplane 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_{x}L_{0})/π 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)≈βρ_{00} 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:
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_{0}.
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 RTIunstable 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:
The length scales associated with the effects of ionneutral 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 ionneutral friction term. Doing so leads to:
where ν_{in}, ν_{ni}, and v_{fast} have been defined in Eq. (14).
Figure 8 shows the mode numbers n corresponding to the ionneutral 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), ionneutral collision scale (orange lines), and neutralion 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 ionneutral collisions affect modes that can exist in our domain, that is, the modes below the horizontal dotted line. However, the scales associated with neutralion 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 singlefluid 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 semianalytically in the singlefluid 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 semianalytically 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 yaxis. 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_{x}L_{0} = πn < 1, RTI growth rates become insensitive to the presence of an inplane magnetic field and converge towards ω = 0 for n → 0. This is consistent with lack of magnetic field confinement of the modes for small k_{x}L_{0} 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^{−6} m^{−1} 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 semianalytically for the perpendicular magnetic field (orange); sheared field with L_{s} equal to L_{0} (red), L_{0}/2 (violet), L_{0}/3 (cyan), L_{0}/4 (lime), L_{0}/5 (blue), L_{0}/10 (dark green); and for uniformly rotated magnetic field by 1° (black). Right: growth rate obtained from the simulations PMMS (blue), PWNS (orange), L1MMS (green), L1WNS (red). 
Fig. 10. Comparison between the growth rates obtained semianalytically 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: PWNS; black lines: PWNNNS; red lines: L1WNS; green lines: L1WNNNS, and blue lines for L1WNNNBS. 
The cases where the equilibrium magnetic field has shear (and, thus, an inplane 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_{x}L_{0} > 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_{0} and 0.1L_{0}, 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_{0} 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 rerun 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_{1} < z < z_{2}, 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_{0}. 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 multimode perturbation “MM” described by Eq. (11). The multimode 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 semianalytical calculations, we observe that there is a cutoff for the growth rates obtained from the simulations which is not present in the semianalytical 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 semianalytical 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 semianalytically (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 semianalytical 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 cutoff is produced for all simulated cases, unlike the semianalytical cases.
We observe that the increase in density contrast leads to the increase in both the peak growth rate and the cutoff 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 cutoff observed in the simulations in Fig. 10. As shown in Fig. 8, the lowest modes estimated to be impacted by neutralion and ionneutral collisions, n_{νni} and n_{νin}, are below those for the viscosity scale, n_{viscn}. Therefore, it is expected that the observed cutoff 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 (L1WNS, L1WNNNBS, L1WNNNS, 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 semianalytically 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 highn 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: L_{0} (red), L_{0}/2 (violet), L_{0}/3 (cyan), L_{0}/4 (lime), L_{0}/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 RTIunstable sheared configurations consistently show Δ ∝ n^{−1}. The wavelength dependence of the vertical extent of the RTI modes in the case with no inplane field and a vertically limited RTIunstable region with (dρ_{0}/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^{−1} 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 highn 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_{0} and L_{0} 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 L1WN (top left), L1WNNN (top right), L1WNNNB (bottom left), and PWNNN (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: L1WN simulation; top right: L1WNNN; bottom left: L1WNNNB; bottom right: PWNNN. 
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 inplane magnetic field. Thus, the L1WNNN simulation (top right) shows larger decoupling than the PWNNN one (bottom right). In a similar vein, since the Lorentz force is larger for larger magnetic field, the decoupling in the L1WNNNB case (bottom left) is generally greater than in the L1WNNN 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 ionneutral (or neutralion) collisional scales as shown in Fig. 8. The increase of the density contrast produces lower decoupling (L1WN vs L1WNNN case). This is also in agreement with the analysis of the ionneutral 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 L1WNNN, top right). The simulation L1WN 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 PWN (top left), PWNNN (top right), L1WNNNB (bottom left), and PWNNN (bottom right). The snapshots are taken at the time indicated at the righthandside 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 isocontours of the absolute value of the ycomponent of the magnetic field potential, A_{y}. The isocontours are equally spaced between 0.6 and , where is the maximum value of A_{y}. 
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 L1WN, L1WNNN, L1WNNNB, and PWNNN. 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_{0} and z = L_{0}.
Fig. 14. Power spectra of magnetic energy for each of the components of the magnetic field calculated at the final time for the simulations L1WN (top left), L1WNNN (top right), and L1WNNNB (bottom left), and PWNNN (bottom right). We note that the simulation PWNNN 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 inplane 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 inplane magnetic energy (, dashed green lines) is greater than in each of the inplane components of the magnetic energy ( and , blue and orange solid lines). This implies that the inplane components of magnetic field energy vary approximately when they are inphase with each other.
On the other hand, the mode amplitudes of the total inplane magnetic field energy and the energy in the outofplane Bfield 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 inplane and outofplane Bfield 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 highn modes of the inplane and the outofplane components of magnetic field energy have the same amplitude and are nearly exactly out of phase.
The scale where the total inplane, , and outofplane, , components transition from being additive to canceling each other seems to be related to the largest size of the inplane plasmoids seen in Fig. 4 for the same time moment. We can observe that the structures formed in the L1WN case are the largest between three L1 simulations, followed by the sizes in the L1WNNNB case, and then by the L1WNNN case. The ordering between these scales is the same as the ordering of the scales where the outofplane and the inplane 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 ydirection being out of phase with the inplane 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 twofluid 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 RTIunstable transition layer as opposed to the typically considered discontinuous interface is both necessary to selfconsistently evaluate twofluid 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 inplane component of the magnetic field suppresses RTI growth in the linear phase. Both in simulations and in analytical calculations using singlefluid 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 wavenumber 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 highn modes even more prone to collisional damping effects.
(iv) When mass loading and the density contrast is increased, both semianalytical ideal MHD calculations and simulations show that the linear growth rate is larger at all scales. Additionally, greater mass loading leads to stronger ionneutral 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) Ionneutral 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 semianalytically under the singlefluid ideal MHD approximation. These effects will be further explored in followup 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 outofplane and inplane components of the magnetic energy show strong outofphase 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 L1WN simulation showing both the strongest magnetic activity and highest degree of such magnetic selforganization across the scales. These effect will be further explored in followon 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 largeaperture 4 m solar telescope facilities (DKIST^{1}, EST^{2}).
In all the simulations performed, we found that ionneutral 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 semianalytically in the singlefluid 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 ionneutral collisions and viscosity, with the global RTI growth rate remaining the same as in the singlefluid 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 rollup 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 ionneutral collisional frequency, is higher. This reduction of viscosity and the importance of ionneutral 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 smallscale perturbations with wavelengths smaller than a cutoff wavelength, proportional to the square of the inplane 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 frozenin 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 semianalytical 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 inplane magnetic pressure effects. However, both analytically estimated and semianalytical calculations have also shown that for a range of shear lengths, in the absence of collisional effects, the RTI growth rate becomes wavelengthindependent in the short wavelength limit.
Using a singlefluid approximation with twofluid 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 neutralcharge interaction and viscosity. This is due, primarily, to our consideration of a smooth rather than a discontinuous prominencecorona 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 inplane 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 outofphase correlation between the inplane and the outofplane magnetic energy for intermediate and high n modes hints at magnetic selforganization 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 inplane 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 PGC2018095832BI00 and the US National Science Foudation. It contributes to the deliverable identified in FP7 European Research Council grant agreement ERC2017CoG771310PI2FA for the project “Partial Ionization: Twofluid Approach”. The authors wish to acknowledge the contribution of Teide HighPerformance 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
 Anan, T., Ichimoto, K., & Hillier, A. 2017, A&A, 601, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anuchina, N., Volkov, V., Gordeychuk, V., et al. 2004, J. Comput. Appl. Math., 168, 11 [Google Scholar]
 Arber, T. D., Haynes, M., & Leake, J. E. 2007, ApJ, 666, 541 [Google Scholar]
 Berger, T. E., Shine, R. A., Slater, G. L., et al. 2008, ApJ, 676, L89 [Google Scholar]
 Berger, T. E., Slater, G., Hurlburt, N., et al. 2010, ApJ, 716, 1288 [Google Scholar]
 Berger, T., Hillier, A., & Liu, W. 2017, ApJ, 850, 60 [Google Scholar]
 Bhatia, P. K. 1974, Ap&SS, 26, 319 [Google Scholar]
 Carlyle, J., & Hillier, A. 2017, A&A, 605, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon Press) [Google Scholar]
 Díaz, A. J., Soler, R., & Ballester, J. L. 2012, ApJ, 754, 41 [Google Scholar]
 Díaz, A. J., Khomenko, E., & Collados, M. 2014, A&A, 564, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dolai, B., Prajapati, R. P., & Chhajlani, R. 2016, Phys. Plasmas, 23, 113704 [Google Scholar]
 Doludenko, A. N., Fortova, S. V., Shepelev, V. V., & Son, E. E. 2019, Phys. Scr., 94, 094003 [Google Scholar]
 Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [Google Scholar]
 Freed, M. S., McKenzie, D. E., Longcope, D. W., & Wilburn, M. 2016, ApJ, 818, 57 [Google Scholar]
 Gardner, C. L., Glimm, J., McBryan, O., et al. 1988, Phys. Fluids, 31, 447 [Google Scholar]
 Gerashchenko, S., & Livescu, D. 2016, Phys. Plasmas, 23, 072121 [Google Scholar]
 Gilbert, H., Kilper, G., & Alexander, D. 2007, ApJ, 671, 978 [Google Scholar]
 GonzálezMorales, P. A., Khomenko, E., Downes, T., & de Vicente, A. 2018, A&A, 615, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, A. 2016, MNRAS, 462, 2256 [Google Scholar]
 Hillier, A. 2018, Rev. Mod. Plasma Phys., 2, 1 [Google Scholar]
 Hillier, A. 2019, Phys. Plasmas, 26, 082902 [Google Scholar]
 Hillier, A., Matsumoto, T., & Ichimoto, K. 2017, A&A, 597, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jun, B.I., Norman, M. L., & Stone, J. M. 1995, ApJ, 453, 332 [Google Scholar]
 Khomenko, E., Díaz, A., de Vicente, A., Collados, M., & Luna, M. 2014, A&A, 565, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Khomenko, E., Collados, M., & Díaz, A. J. 2016, ApJ, 823, 132 [Google Scholar]
 Leake, J. E., DeVore, C. R., Thayer, J. P., et al. 2014, Space Sci. Rev., 184, 107 [Google Scholar]
 Leonardis, E., Chapman, S. C., & Foullon, C. 2012, ApJ, 745, 185 [Google Scholar]
 Livescu, D. 2004, Phys. Fluids, 16, 118 [Google Scholar]
 Mitra, A., Roychoudhury, R., & Khan, M. 2016, Phys. Plasmas, 23, 024503 [Google Scholar]
 Newcomb, W. A. 1983, Phys. Fluids, 26, 3246 [Google Scholar]
 Orozco Suárez, D., Díaz, A. J., Asensio Ramos, A., & Trujillo Bueno, J. 2014, ApJ, 785, L10 [Google Scholar]
 Parchevsky, K. V., & Kosovichev, A. G. 2007, ApJ, 666, L53 [Google Scholar]
 Popescu Braileanu, B., Lukin, V. S., Khomenko, E., & de Vicente, Á. 2019, A&A, 627, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ribeyre, X., Tikhonchuk, V. T., & Bouquet, S. 2005, Ap&SS, 298, 75 [Google Scholar]
 Ruderman, M. S., Terradas, J., & Ballester, J. L. 2014, ApJ, 785, 110 [Google Scholar]
 Ruderman, M. S., Ballai, I., Khomenko, E., & Collados, M. 2018, A&A, 609, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Terradas, J., Soler, R., Luna, M., Oliver, R., & Ballester, J. L. 2015, ApJ, 799, 94 [Google Scholar]
 Vandervoort, P. O. 1961, ApJ, 134, 699 [Google Scholar]
 Wiehr, E., Stellmacher, G., & Bianda, M. 2019, ApJ, 873, 125 [Google Scholar]
 Xue, C., & Ye, W. 2010, Phys. Plasmas, 17, 042705 [Google Scholar]
 Youngs, D. L. 1984, Phys. D: Nonlinear Phenom., 12, 32 [Google Scholar]
 Zhou, Y. 2017a, Phys. Rep., 720, 1 [Google Scholar]
 Zhou, Y. 2017b, Phys. Rep., 723, 1 [Google Scholar]
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 −L_{0} and L_{0} and z varies between −4L_{0} and 4L_{0}. 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, B_{y}, 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 antisymmetric 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 yaxis 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 m^{2} s^{−1}, as from Eq. (14). The neutral viscosity and thermal conductivity coefficients are limited to 10^{9} m^{2} s^{−1}. Right: ionneutral (ν_{in}, blue lines), neutralion (ν_{ni}, orange line) collision frequencies, and frequencies associated with waves: fast magnetoacoustic (ω_{fast}, green lines), inplane 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 L1WN (top left group of panels), L1WNNN (top right group of panels), L1WNNNB (bottom left group of panels), PWNNN (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 (L1WN, L1WNNN, L1WNNNB), the isocontours of the absolute value of the ycomponent of the magnetic field potential, A_{y}, are shown at the bottom panels. The isocontours are equally spaced between 0.6 and , where is the maximum value of A_{y}. 

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 −L_{0} and L_{0}. Different panels are for different simulations: L1WN (top left); L1WNNN (top right); L1WMNNB (bottom left); PWNNN (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 inplane component of the magnetic field as a function of height, computed from Eq. (22) for the equilibrium profile with the values of n_{n00} = 10^{16} m^{−3} and B_{00} = 10^{−3} T in Eqs. (5) and (7). Cutoff wavelengths for the following four shear lengths: L_{s} = [0, L_{0}/10, L_{0}/2, L_{0}] are shown with increasing width of the line for increasing L_{s}. The red dotted line corresponds to the density gradient scale length L_{d}. 

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), ionneutral collision scale (orange lines), and neutralion 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 semianalytically for the perpendicular magnetic field (orange); sheared field with L_{s} equal to L_{0} (red), L_{0}/2 (violet), L_{0}/3 (cyan), L_{0}/4 (lime), L_{0}/5 (blue), L_{0}/10 (dark green); and for uniformly rotated magnetic field by 1° (black). Right: growth rate obtained from the simulations PMMS (blue), PWNS (orange), L1MMS (green), L1WNS (red). 

In the text 
Fig. 10. Comparison between the growth rates obtained semianalytically 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: PWNS; black lines: PWNNNS; red lines: L1WNS; green lines: L1WNNNS, and blue lines for L1WNNNBS. 

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: L_{0} (red), L_{0}/2 (violet), L_{0}/3 (cyan), L_{0}/4 (lime), L_{0}/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: L1WN simulation; top right: L1WNNN; bottom left: L1WNNNB; bottom right: PWNNN. 

In the text 
Fig. 13. Pairs of snapshots of the horizontal and vertical decoupling in velocities for the simulations PWN (top left), PWNNN (top right), L1WNNNB (bottom left), and PWNNN (bottom right). The snapshots are taken at the time indicated at the righthandside 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 isocontours of the absolute value of the ycomponent of the magnetic field potential, A_{y}. The isocontours are equally spaced between 0.6 and , where is the maximum value of A_{y}. 

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 L1WN (top left), L1WNNN (top right), and L1WNNNB (bottom left), and PWNNN (bottom right). We note that the simulation PWNNN 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 inplane component, (green), and the total energy, (black). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.