EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 565 (May 2014) A&A, 565 (2014) A45 Full HTML
Free Access
Issue
A&A
Volume 565, May 2014
Article Number A45
Number of page(s) 15
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201322918
Published online 05 May 2014

© ESO, 2014

1. Introduction

Solar prominences are composed of cool, dense, and partially ionized plasma and their large-scale magnetic structures remain stable for days, or even weeks, in the solar corona. Prominence material is believed to be supported by the magnetic field (see reviews by Tandberg-Hanssen 1995; Mackay et al. 2010). There are several large-scale models that address the problem of the global stability of prominences, and of the origin of their mass that may explain observational properties (Pneuman 1983; van Ballegooijen & Martens 1989; Priest et al. 1989; Antiochos & Klimchuk 1991; Rust & Kumar 1994; Antiochos et al. 1994; Aulanier & Demoulin 1998; DeVore & Antiochos 2000; Gibson & Fan 2006; Aulanier et al. 2006).

Against the background of global stability, prominences are extremely dynamical at small scales, showing a variety of shapes and moving with vertical and horizontal threads (Berger et al. 2010; Ryutova et al. 2010). The dynamical appearance of prominences depends on whether they are located above active or quiet regions, and on the relative orientation with respect to the observer. Quiescent prominences usually occur in quiet regions at high latitudes. They usually reach larger heights than active region prominences, and often show very characteristic vertical threads of less than 1 Mm in size, and downflow velocities of about 10−20 kms-1 along them1. The origin of these downflowing drops has been addressed in many studies (Gilbert et al. 2002, 2007; van Ballegooijen & Cranmer 2010; Haerendel & Berger 2011), aiming to explain the apparent material motion across the magnetic field lines, the speed of the drops, and the amount of the mass loss and gain.

Berger et al. (2010) find large-scale 20−50 Mm arches, expanding from the underlying corona into the prominences. At the top of these arches, at the prominence-corona transition region (PCTR), there are observed dark turbulent upflowing channels of 4−6 Mm maximum width with a profile typical of the Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instabilities. The upflows rise up to 15−50 Mm, with an average speed of 13−17 kms-1, decreasing at the end. Lifetimes of the plumes are about 300−1000 s. These numbers fit well into the theoretical predictions from the classical theory of plasma instabilities (Isobe et al. 2005; Berger et al. 2008, 2010, 2011; Heinzel et al. 2008; Ryutova et al. 2010) and can therefore be used to derive plasma parameters of the prominences (Hillier et al. 2012b). Plumes and spikes are seen at any time and any possible orientation of the limb portion of a prominence, but they are especially evident in quiescent hedgerow prominences. An alternative explanation for the upflowing plumes was recently proposed by Dudík et al. (2012) as due to the presence of separatrix layers and reconnection, arguing that RTI cannot happen for the magnetic field orientation and plasma parameters expected for prominences.

From the theoretical point of view, the existence of the instabilities at the PCTR is easily explained since the two media clearly have different densities, temperatures, and relative velocities. The analytical linear magnetohydrodynamic (MHD) theory of these instabilities is well developed (e.g., Chandrasekhar 1981; Priest 1982). Numerical simulations in the non-linear regime have been performed for different astrophysical contexts in two and three dimensions (Jun & Norman 1995; Jun et al. 1995; Arber et al. 2007; Stone & Gardiner 2007a,b; Isobe et al. 2006). Recent numerical MHD simulations of the RTI in prominences by Hillier et al. (2012a,c), including a rising buoyant tube in a Kippenhahn-Schlüter prominence model, show good agreement with observations.

Prominences are relatively cool and dense objects, with values of temperature and density in the chromospheric range (Tandberg-Hanssen 1995). Therefore, a prominence material is expected to be only partially ionized. The presence of a large number of neutrals must affect the overall dynamics of the plasma, since neutrals are not affected by the magnetic Lorentz force directly, but only through the collisional coupling to ions. The aim of our work is to model the dynamics of the Rayleigh-Taylor instability in the partially ionized prominence plasma in the non-linear regime. The linear theory of the Rayleigh-Taylor and Kelvin-Helmholtz instabilities in a partially ionized plasma has been recently developed by Soler et al. (2012) and Díaz et al. (2012, 2014). Different approaches have been used, including a two-fluid and a single-fluid modeling. It is known from the ideal MHD that the magnetic field parallel to the perturbation interface stabilizes the system, up to a threshold wavelength λc, (1)where B0cosθ is the value of the magnetic field in the plane of the perturbation, and ρ2 − ρ1 is the density contrast between the two media. For a given value of the magnetic field, perturbations with a wavelength shorter than λc are stable, while large wavelength perturbations remain unstable. The results of the linear theory show that, in a partially ionized plasma, there is no critical wavelength, and perturbations in the whole wavelength range are always unstable (Soler et al. 2012; Díaz et al. 2012, 2014). The growth rate of a given perturbation in a partially ionized plasma depends on (among other parameters) the ionization fraction and is rather small for the ionization fractions expected for the prominence plasma.

The linear instability theory has the advantage of being analytical, but then a number of simplifications are necessary, limiting the range of its validity. In order to consider the fully non-linear evolution of the Rayleigh-Taylor instability, and to include more complex physics, numerical simulations have to be performed. Here we report on 2.5D numerical simulations of the RTI at the boundary between the hot corona and a cool prominence taking into account plasma partial ionization by means of the generalized Ohm’s law.

thumbnail Fig. 1

Initial configuration. The instability develops in the XZ plane containing the perturbation vector k; the initial magnetic field B0 forms an angle θ with the XZ plane. B is initially lying in the XY plane, parallel to the interface separating the prominence and coronal plasma, with Bz = 0, Bx = B0cos(θ), and By = B0sin(θ).

2. Initial configuration

Observations of quiescent prominences reveal low field strengths, in the range of 3−30 G (Leroy 1989). The field direction is deduced from many observations to be mainly horizontal (making an acute angle with the axis of the prominence), although there is no unique opinion about this (Mackay et al. 2010).

In the following, we assume the simplest plasma and magnetic configuration, summarized in Fig. 1. We choose an initially homogeneous horizontal magnetic field with a strength of B0 = 10 G, in prominence and corona. We consider two magnetic field orientations, B0 at θ = 90° to the X axis (i.e., normal to the XZ plane), and B0 at θ = 89° to the X axis (slightly skewed from the normal to XZ plane). The plasma is most unstable for these field orientations. For the perturbations forming a larger angle θ with the magnetic field the instability would take a longer time to develop, otherwise leading to similar simulation results (see Khomenko et al. 2013).

Since our aim is to study the effects of non-ideal plasma terms due to neutrals, we choose a high spatial resolution of 1 km in both the X and Z directions. We simulate only a small portion of the interface between the prominence and the corona, with a size of 250 × 1000 km. The equilibrium magnetic field is homogeneous and does not influence the force balance. Pressures and densities are obtained from the hydrostatic equilibrium condition, given the temperature stratification. The temperatures of the prominence and corona are initially constant with a smooth transition. This equilibrium is different from what is usually assumed for prominences, where the magnetic field is expected to exert a force on the plasma and prevents it from falling. Despite this, since we only simulate a small portion of the interface, we expect that the effects of the curvature of the equilibrium magnetic field are not significant. The results presented here represent an initial exploratory study and this configuration allows us to make a direct comparison with the linear theory where a similar configuration is also assumed (Díaz et al. 2014).

Despite the advantage of its simplicity, the homogeneous current-free field configuration is unfavorable for the manifestation of non-ideal plasma effects. The importance of these effects depends on currents (terms proportional to ηAJ⊥ in the Eqs. (2) below). As the instability is being developed, the currents (J⊥) appear only because of perturbations in the magnetic field, at the interface between the two media, and their influence is limited to a small region (PCTR), unlike in a non-potential magnetic field configuration. This lack of currents is one of the reasons why partial ionization effects are relatively local and have not been considered so far in theoretical studies. On the other hand, the ionization degree of the prominence plasma can be very low because of its low temperature, and the large fraction of neutrals may partially compensate for the weakness of currents and increase to significant values the diffusion coefficient due to neutrals, i.e., ηA, so the terms proportional to ηAJ⊥ become significant. Table 1 summarizes the parameters of our equilibrium configuration. For these parameters, Eq. (1) gives the critical wavelength λc = 38 km for θ = 89°. There is no critical wavelength in the θ = 90° case.

3. Numerical solution

We solve numerically the non-ideal MHD equations of conservation of mass, momentum, internal energy, and the induction equation (Khomenko & Collados 2012),

(2)where the gravity is equal to g = 2.74 × 102 m s-2 and the adiabatic index γ = 5/3. The variables ρ and p are summed over the plasma components (electrons (e), ions (i), and neutrals (n)), the center of mass velocity u is an average over the velocity of individual species, (3)and J⊥ is the component of the current perpendicular to the magnetic field, (4)In these equations we neglect the non-diagonal components of the pressure tensor and assume low drift velocities wα = u − uα (α = e,i,n), i.e., | wα | ≪ | u | , | uα | (see Sect. 4.3 for the discussion of the value of the ion-neutral drift momentum in our simulations). The latter condition allows us to remove the terms containing . In the Ohm’s law we neglected the time variation in the relative ion-neutral drift velocity ui − un. This can be safely done since the expected instability growth rate is significantly below the ion-neutral collisional frequency νin ~ 104 s-1 (Díaz et al. 2014). We only keep Ohmic and ambipolar2 diffusion terms. The linear analysis reveals that for the typical conditions of prominence plasma the ambipolar diffusion term dominates by orders of magnitude over the neglected Hall and battery terms (Díaz et al. 2014). The ambipolar diffusion coefficient is equal to (Braginskii 1965), (5)The Ohmic diffusion coefficient is given by (6)and theoretical values of the collisional frequencies are (Spitzer 1968; Braginskii 1965) (7)where e is the electron charge, ne is the electron number density, min = mimn/ (mi + mn) and men = memn/ (me + mn). The respective cross-sections are σin = 5 × 10-19 m2 and σen = 10-19 m2. The Coulomb logarithm Λ for T< 50 eV (T< 600.000 K) is equal to: (8)Equations (2) were solved by means of our code mancha (Khomenko et al. 2008; Felipe et al. 2010; Khomenko & Collados 2012) with the inclusion of the physical ambipolar diffusion term in the equation of energy conservation and in the induction equation. We evolve the electron number density ne in time by means of the Saha equation. The simulations are done in 2.5D, meaning that the derivatives are done only in the X and Z directions, but the vector quantities (velocity and magnetic field) are allowed to have three components. We use periodic horizontal boundary conditions and closed vertical boundary conditions.

Table 1

Parameters of the equilibrium configuration.

The code mancha evolves (non-linear) perturbations to an equilibrium state. In order to initiate the instability in the current experiment we proceeded in the following way. We constructed two models: model 1 is a horizontally homogeneous coronal model in hydrostatic equilibrium at coronal temperature Tcor from Table 1 and model 2 is a model of corona with prominence on the top. The vertical hydrostatic equilibrium in model 2 is imposed in each column with temperatures varying from Tcor (at the lower part of the domain) to Tprom (at the upper part) and the corona−prominence interface corrugated according to Eq. (9). The difference between model 2 and model 1 was used as a perturbation to initiate the instability. The position z where the temperature changes from coronal to prominence values in model 2 varies with horizontal location x, otherwise the perturbation is in equilibrium and does not evolve. We used the multi-mode perturbation to the corona–prominence interface from Jun et al. (1995)(9)where L = 250 km and am,bm,ϕm are randomly generated amplitudes and phases, and zi = 800 km is the location of the interface. The amplitudes am and bmvary in the range 0−10 km. We used M = 25 modes from 10 km to 250 km in wavelength. In the θ = 89° simulation, under ideal MHD conditions, three-quarters of these modes with λ<λc are expected to be damped due to magnetic field stabilization effects.

Our code uses hyperdiffusivity for stabilizing the numerical solution. One has to be extremely careful with the choice of the hyperdiffusivity amplitude and model since hyperdiffusive terms are numerical analogs of the physical dissipative terms, such as viscosity, conductivity, and Ohmic diffusion (Stein & Nordlund 1998; Caunt & Korpi 2001; Vögler et al. 2005; Felipe et al. 2010) and their action may resemble the action of the physical terms. In particular, the numerical diffusive terms are able to amplify the growth rate of small-scale modes of the RTI, the same as the physical diffusion would do. To avoid this problem in our simulations, we kept the amplitude of artificial hyperdiffusive terms 2−3 orders of magnitude smaller than the physical ambipolar diffusion. This way we have assured that the characteristic timescales of action of the hyperdiffusive terms are orders of magnitude larger than that of the physical ambipolar diffusion term and that their influence is not significant during the temporal interval covered by the simulations.

In the following we will refer to the simulations with ambipolar diffusion term “on” as AD simulations; and the ones with the ambipolar term “off” as MHD simulations.

thumbnail Fig. 2

Time evolution of density (top) and pressure (bottom) in the AD simulation with θ = 90°. The size of each snapshot is 250 × 1000 km, the time elapsed is given at the bottom of each panel. The velocity field is indicated by arrows. We note the asymmetry between the large-scale rising bubbles and small-scale downflowing fingers in the density images.

thumbnail Fig. 3

Same as Fig. 2 but for the MHD simulation.

4. Results

4.1. Simulations with B normal to the perturbation plane

Figures 2 and 3 give the time evolution of density and pressure in the AD and MHD simulations with θ = 90°. Prior to discussing this simulation we verified that the behavior of the instability in the non-linear regime follows the well known pattern from previous simulations by other authors. Figure 4 shows the dominant spatial scale in the simulations with θ = 90° (both AD and MHD cases are shown) as derived via a Fourier transform over the whole domain and calculating the wavelength (mode) of maximum power at each time moment. In this calculation we performed a 2D Fourier transform and then averaged the power corresponding to the vertical direction. Figure 4 can be directly compared with, e.g., Jun et al. (1995) and shows that large scales tend to dominate in the non-linear stage of the instability since smaller bubbles merge into the larger structures, in agreement with the classical behavior seen in earlier studies (Youngs 1984; Gardner et al. 1988). This tendency is similar in the AD and MHD cases. The situation with B normal to the perturbation plane is equivalent to a purely hydrodynamical case as there is no cutoff wavelength. In the linear regime, the growth rate of small-scale modes is faster than of the large-scale modes (see, e.g., Chandrasekhar 1981; Priest 1982). However, in the non-linear regime large-scale modes have higher terminal velocities and with time overtake the small-scale modes and dominate the evolution of the system (Jun et al. 1995; Wang & Robertson 1985; Isobe et al. 2006; Stone & Gardiner 2007a,b; Hillier et al. 2012a).

The time evolution of pressure and density in Figs. 2 and 3 shows the development of very small scales. We notice the presence of inverse mushroom shapes (especially evident in the pressure images from 26 s to 59 s), typical for the secondary Kelvin-Helmholtz instability (Youngs 1984; Gardner et al. 1988; Jun et al. 1995). The comparison of the AD and the MHD simulations reveals a different particular form of the turbulent flows (but statistically equivalent, as we show below), since the ambipolar diffusion, acting on small scales, forces their different evolution. Particularly interesting is the difference in scales between the density and pressure variations. Pressure (bottom panel of Fig. 2) is initially maintained constant across the border between the prominence and the corona, forced by the hydrostatic equilibrium condition. However, the pressure scale height in the prominence (upper) part of the domain is small, causing the visible decrease in the pressure toward the upper boundary (blue). The turbulent flow of the instability mixes up the fluids and tends to smooth the transition. The density variations (upper panel) have a pronounced asymmetry between the upward rising bubbles and downflowing fingers. The material flows down at small scales. This asymmetry is caused by mass conservation. The hot and less dense coronal material rises up and easily expands forming a large bubble. To move the same amount of material, the dense prominence plasma has to occupy a much narrower volume on its downward motion. The density images can be taken as a proxy for the Hα observational images. We can speculate that in this type of imaging observations, the upward rising bubbles should be significantly more visible (as dark structures) because of their larger scales, producing a visible asymmetry between the upward and downward moving material.

thumbnail Fig. 4

Mode number (m, see Eq. (9)) with maximum Fourier amplitude of the pressure perturbation in the AD (filled circles) and MHD (open circles) simulations with θ = 90°. The wavelength of the modes varies from 10 km (mode 25) to 250 km (mode 1). The amplitude of the modes is indicated by color from blue (smaller) to red (larger). Dotted line is an exponential function ~exp( − t/ 16) shown for illustration purposes.

thumbnail Fig. 5

Velocity histograms over the whole simulation time for the θ = 90° MHD and AD runs. We note the asymmetry of the histogram and the slight tendency of the AD simulation to develop higher velocities.

thumbnail Fig. 6

Growth rate of the RTI from the linear theory according to Díaz et al. (2014) as a function of the wave number of the perturbation along the discontinuity. The parameters of the atmosphere assumed for the linear calculation are those from Table 1; the magnetic field is inclined with respect to k by θ = 89° (dashed line) and by θ = 90° (dotted line). Solid line shows the case of θ = 89° with the ambipolar diffusion term taken into account.

Figure 5 provides a histogram of the vertical velocities (positive indicates upward motion) around the PCTR3 in the AD simulation with θ = 90° (red line), compared to the MHD simulation (black line). Both histograms are asymmetric in the sense that upflows are fewer (occupying only 1/3 of all points) but stronger, and the average velocity is zero at any time. We note that the prevalence of downflowing points does not contradict the statement that the downflowing regions appear as small-scale fingers on the density images in Figs. 2 and 3, since the less dense coronal material around the fingers is also involved in the downflowing motion, but is not apparent on the density images. There is a slight, but indicative, difference between the AD and the MHD histograms, particularly apparent in their wings. The red line goes systematically above the black line, meaning that there are more cases of extreme velocities in the AD case. To quantify this difference, we measured the percentage of points with the absolute value of velocity above 5 kms-1, obtaining 12% for the MHD simulation and 16% for the AD simulation. This result indicates that the RTI in the non-linear regime develops slightly higher velocities when taking into account the influence of the neutral component, compared to the purely MHD case.

As a next step we analyzed the scales developed in the non-linear regime in the simulation with θ = 90°. As mentioned above, the case of the magnetic field normal to the plane of the instability is analogous to a purely hydrodynamical case since there is no cutoff wavelength. In the linear regime, the small-wavelength perturbation should grow faster, but in the non-linear regime, the large-scale perturbations become more important due to non-linear interaction of the turbulent flows. This behavior is reproduced both in the MHD and AD simulations, as can be seen in Fig. 4. In the non-linear regime, even for the θ = 90° case, the plasma motion generates currents in the plane of the perturbation. These currents dissipate through the corresponding term in the energy equation and cause diffusion effects, especially affecting small scales in the simulation. The dissipative action of the ambipolar term cannot be taken into account in the linear theory since the dissipative term in the energy equation is non-linear and is removed from the linear analysis (Díaz et al. 2014).

The instability growth rate from the linear theory for θ = 90° is given by the dotted line in Fig 6, according to the calculations presented in Díaz et al. (2014). This figure indicates that, depending on the wave number, the timescale of the development of the instability is on the order of tens of seconds. For wavelengths around λ = 2π/k = 100 km the timescale is t = 2π/ Im(ω) ≈ 50 s, and decreases to t ≈ 25 s for λ ≈ 30 km. Figures 2 and 3 show that the instability is already outside of the linear regime at about 10 s of the simulation. Thus, our numerical results are not directly comparable to the linear theory, but we can take the linear theory as an indication.

Figure 7 shows the Fourier analysis of scales of the AD simulation, relative to the MHD simulation. In view of the above considerations we have selected two zones of the domain, one of them around the discontinuity (PCTR), and another one above it. The ambipolar diffusion coefficient is only significant in the prominence part of the domain above the discontinuity, and therefore we do not consider the coronal part of the domain for the Fourier analysis. At each time moment we calculated the power as a function of horizontal wave number k along the discontinuity by 2D Fourier transform in space of the corresponding zone of the snapshot of pressure variations (as the variable that develops more fine structures), and by averaging the power for the vertical wave number to decrease the noise. At each time moment we re-normalized the total power to maintain it constant in time. This way we obtain a power as a function of horizontal wavelength and time at two domains, at and above the discontinuity. This procedure was applied separately to the AD and MHD simulation. We then divided the AD power map by the corresponding MHD map in each domain. The result is given in the left panels of Fig. 7.

Above the discontinuity (bottom-left panel in Fig. 7), the ambipolar diffusion term acts as a pure diffusion quickly removing all small scales and amplifying large scales relative to them. At the PCTR, the behavior of Fourier harmonics is strikingly different. There we do not observe any significant change in power between the AD and MHD cases. The relative power fluctuates in time, but the average remains around one for all harmonics. This result indicates that the ambipolar term can act in a very different way in situations where the instability is being developed at PCTR and outside this region.

4.2. Simulations with B at 89° to the perturbation plane

The time evolution of density and pressure in the AD simulation with B skewed from the normal direction is given in Fig. 84. By just rotating the field by 1°, the scales developed in the simulation completely changed: the small scales disappear and only one big drop is developed after about 100 s of the simulation. The small scales cannot develop the instability because of the cutoff induced by the magnetic field, λc ≈ 38 km for the parameters of our atmospheric model in equilibrium. While the growth rate of small-scale modes (λ<λc) is not zero as in the purely MHD case, it is still smaller compared to large scales (see Fig. 5 in Díaz et al. 2014). One can appreciate from the figure that at t ~ 10 s the dominant wavelength is around λ = L/ 6 ≈ 40 km, while at t ~ 47 s it becomes λ = L/ 2, and finally after t ~ 100 s the dominant wavelength is equal to the whole size of the box in the X direction, L. Similar to the θ = 90° case, this behavior is a consequence of the non-linear interaction of flows, powering energy from small to large scales (see, e.g., Jun et al. 1995; Wang & Robertson 1985).

Several interesting phenomena happen during the evolution of the RTI in Fig. 8. After being formed at about 60−70 s of the simulation, the big drop starts falling at a speed of about 3−5 kms-1. Similar ranges of velocities are also found in the simulations of Hillier et al. (2012a,c). On its way down, the drop deforms the magnetic field lines, compressing them. One can get a false impression from the figure that the magnetic field lines are significantly twisted and even become vertical. This is not the case because the magnetic field lines are actually outside of the plane XZ, directed towards us. Figure 9 provides the 3D rendering of the magnetic field lines in a single snapshot, showing that, as the drop evolves, the magnetic field lines are piled up together and slightly inclined and twisted, otherwise remaining essentially horizontal.

thumbnail Fig. 7

Fourier analysis of scales developed in the θ = 90° simulations. Panels on the left show the relative power in the AD simulation divided by the power in MHD simulation as a function of horizontal wave number along the discontinuity (k) and time. Red means that the AD simulation has more power. Panels on the right give the time average of the power ratio from the panels on the left.

thumbnail Fig. 8

Time evolution of density (top) and pressure (bottom) in the AD simulation with θ = 89°. The size of each snapshot is 250 × 1000 km; the time elapsed is given at the bottom of each panel. Green lines are projections of the magnetic field lines into the XZ plane; the velocity field is indicated by arrows. We note the difference in scales compared to the θ = 90° case.

The time evolution of the velocity averaged over the drop is shown in the upper-right panel of Fig. 10. We observe that after t ~ 100 s the downward motion of the drop is slowed down and stops at t ~ 160 s. Then the drop starts rising again, extending in the horizontal direction (Fig. 8). It continues to oscillate around the equilibrium position and the amplitude of this oscillation decreases in time, until the motion finally ceases. The forces responsible for this behavior are given in the left panels of Fig. 10. At the beginning of the evolution, the magnetic pressure gradient force almost balances the gas pressure gradient force, the magnetic tension force is small, and the gravity force dominates and causes the downward acceleration of the drop. The magnetic pressure gradient force remains the leading one to compensate the gravity over the whole evolution. The increase in the magnetic pressure is produced as the drop compresses the magnetic field lines on its fall. The variations in the magnetic pressure essentially cause the oscillation of the total acceleration (bottom left panel) and of the velocity (upper right panel). At the end of the evolution both gas pressure and magnetic pressure gradient forces act together to balance the gravity of the drop and stop its motion. We note that the force equilibrium produced in this simulation is different from the one found in Hillier et al. (2012c) where the main force to compensate the gravity was the magnetic tension force. This difference is caused by the different initial equilibrium of the system. In our case, the magnetic field is initially homogeneous and does not exert any force on the plasma, while in the Kippenhahn-Schlüter model of the prominence adopted by Hillier et al. (2012c), the plasma is supported against gravity by the magnetic tension.

Besides the force balance, an increase in the horizontal field component below the drop causes an increase in the effective cutoff wavelength of the instability. The bottom right panel of Fig. 10 shows the λc estimated from Eq. (1) taking the values of the horizontal field component, Bcos(θ), immediately below the drop. One can see that λc quickly increases and becomes larger than the size of the domain, L, already at t ~ 80 s. Although λc given by Eq. (1) is only a qualitative estimate, it indicates that perturbations with wavelength below L cannot grow further.

The histograms of the velocities in the simulation of θ = 89° are similar to that shown in Fig. 5 except that the range of the velocity variations is narrower, ± 10 kms-1. We do not show the corresponding figure. Instead, we investigate if there are systematic differences between the velocity field in the AD and MHD simulations. The structures that develop in both simulations are similar (unlike the turbulent flows in the θ = 90° simulation); therefore, we can directly compare the snapshots of the vertical velocity in both cases. The left panel of Fig. 11 gives the difference between the standard deviations of the vertical velocity measured around the discontinuity in the AD and the MHD runs. The difference is measured in percentage change from the value in the MHD case. From the very beginning of the instability the standard deviation of the velocity in the AD run is 2−3% larger. The difference is not large, but it is systematic, and the velocities in the AD case are slightly higher. A similar conclusion is obtained from the right panel of Fig. 11 where we show the histogram of the differences in the absolute velocity values around the discontinuity between the AD and the MHD runs. The histogram shows the prevalence of positive values (62%) so that the absolute value of the velocity is larger in the AD simulation, i.e., the instability develops slightly faster in this case.

To evaluate the responsible scales and to demonstrate the higher growth rate in the AD run, we perform the same Fourier analysis as in Fig. 7, but for the θ = 89° simulation. The results are given in Fig. 12. as in Fig. 7, we divide the simulation domain into two zones, at the discontinuity and above it. The power ratio in the prominence part of the domain above the discontinuity behaves in the same way as in the θ = 90° run, the ambipolar diffusion term acting as pure diffusion and quickly removing small scales. At the discontinuity, the power ratio is significantly different from the θ = 90° case. Now there is a clear increase in the growth rate of the small-scale harmonics in the AD simulation compared to the MHD case. The growth rate of the large-scale harmonics is the same in both cases. The change in the behavior happens at about λ ~ 30 km, as is expected, because this value is close to the cut-off wavelength λc. This result confirms and extends the conclusion from the linear theory (Díaz et al. 2014); see Fig. 6. The solid line in Fig. 6 gives the growth rate from the linear theory for the parameters of our model and field orientation of θ = 89°, taking the ambipolar diffusion term into account. Compared to the MHD case (dashed line in Fig. 6), all RTI modes become unstable when the presence of neutral atoms is accounted for in the analysis. We obtain up to a 50% increase in the small-scale harmonics growth rate, in good agreement with the linear analysis.

Figure 13 shows the temperature difference between the snapshots of the AD and MHD simulations. It follows that appreciable temperature differences exist at all stages of evolution of the RTI. Initially, there are larger temperatures at the discontinuity and in the prominence part of the domain in the AD case due to current dissipation and heating (see Khomenko & Collados 2012). The interior of the drop is maintained about 5−10% hotter in the AD at times up to 160 s. After the initial stage, the neutral and ionized material is mixed, the shape of the drops starts to be progressively more different in both runs and this causes a large temperature difference at the transition between the drop and the corona. All in all, the temperature differences are not small, reaching at maximum ~60 kK, or about 30% of the actual value.

4.3. Ion-neutral drift

The relative velocity between the neutral and the ionized plasma components (drift velocity) can be used as a parameter to evaluate the importance of non-ideal plasma effects due to the presence of neutrals. An approximate expression for this velocity can be derived by combining the momentum equation for electrons, ions, and neutrals and neglecting the electron inertia and the time derivative of the ion-neutral drift velocity: (10)For the derivation of Eq. (10), see, e.g., Díaz et al. (2014). The meaning of the variables is the following: (11)is the partial pressure gradient term; ξn = ρn/ρ and ξi = ρi/ρ are neutral and ion fractions, respectively; pe and pn are the electron and neutral pressure, respectively; and αn is the neutral collisional frequency, calculated from the expression (12)where ni is the ion number density, and the collisional frequencies νin and νen are given by Eq. (7). The small parameter ϵ is defined as (13)Keeping in mind all approximations used to derive Eq. (10), once the drift velocity is obtained, we can recover the individual velocities of ions and neutrals from the linear system of Eqs. (3) and (10):

(14)The last two terms in Eq. (10) are small, and we keep only the two leading ones for the analysis. The drift velocity defined this way is produced by cross-field currents (first term) and by gradients of the partial pressures (second term). A high drift velocity means that the species are not entirely coupled together by collisions.

thumbnail Fig. 9

Three-dimensional representation of the magnetic field lines in a snapshot of the AD simulation with θ = 89° at t = 160 s. The vertical Z axis is 1000 km long. The density of the field lines is proportional to the magnetic field strength and their origins in XZ plane is random. Projections of the field lines into XZ and YZ planes are indicated in green color.

The individual velocities of species do not carry much information themselves, since they do not reflect the amount of material moving with a given velocity. This information is contained in their momenta: (15)For this reason, following Pandey & Wardle (2008), we define the ion-neutral drift momentum as (16)To preserve the sign of the drift momentum we do not use the square of this quantity, unlike, e.g., Pandey & Wardle (2008) and Martínez-Sykora et al. (2012). We have calculated the neutral, ion, and drift momenta a posteriori from the parameters of the AD θ = 89° simulation and the result is displayed in Fig. 14.

At the beginning of the evolution (t = 10 s), the upper prominence part of the domain has a very low ionization fraction, and the lower coronal part is completely ionized, but the absolute value of density in the corona is about 100 times lower than in the prominence. Significant velocities appear only at the interface and, as a consequence, both ion and neutral momenta are only significant at the PCTR. Later on (t = 45 − 120 s) significant negative neutral momentum appears in the prominence part of the domain, pushing the drop downwards. The sharp transition in the ionization fraction between the prominence and the corona disappears, producing a smooth increase in the ionization fraction at the interface and negative ion momentum just at the border of the drop. At times t = 160 − 195 s, the negative neutral momentum ceases, and there is a positive ion momentum surrounding the drop and pushing it upwards. Later, at t = 235 − 275 s, the ionized and the neutral fluids inside the drop mix together, and the neutral momentum becomes small due to the increase in the ionization fraction inside the drop.

thumbnail Fig. 10

Top left: acceleration as a function of time due to the gravity force (dash-dotted red line), gas pressure gradient force (triple-dot-dashed green line), magnetic pressure gradient force (solid black line), and magnetic tension force (dashed black line). The values are integrated over the big drop from Fig. 8 (AD run) after it forms at about 70 s of the simulation time. Bottom left: acceleration due to the sum of all above forces. Top right: vertical velocity of the big drop. Vertical dotted lines indicate time instants of zero velocity. Bottom right: effective cutoff wavelength calculated from Eq. (1) for the horizontal field component Bcos(θ) at the bottom of the drop. Horizontal dotted lines indicate the size of the computational domain L = 250 km.

thumbnail Fig. 11

Left: difference between the standard deviations of the vertical velocity near the discontinuity in the AD and MHD θ = 89° simulations during the first 120 s. Right: histogram of the difference between the absolute values of the vertical velocity near the discontinuity in the AD and MHD θ = 89° simulations during the first 120 s.

thumbnail Fig. 12

Fourier analysis of scales developed in the θ = 89° simulations. The format of the figure is the same as Fig. 7. The vertical dotted lines in the right panels mark the cutoff wavelength λc = 38 km.

thumbnail Fig. 13

Difference between the temperature variations in the AD and MHD simulations with θ = 89°, (TA − TMHD) /TMHD. The size of each snapshot is 250 × 1000 km; the time elapsed is given at the bottom of each panel. Red colors mean that the AD model is hotter.

thumbnail Fig. 14

Top: time evolution of the vertical component of the neutral momentum from Eq. (14) in the θ = 89° AD run. Negative values mean downward motion. The contours of vertical drift momentum and −10-9 kg m-2 s-1 (see Eq. (16)) are overplotted in green and gray, respectively. Bottom: the same for the ion momentum. The size of each snapshot is 250 × 1000 km; the time elapsed is given at the bottom of each panel, the color scheme is constant in time, and is the same in both panels.

The contours in Fig. 14 show the time evolution of the drift momentum calculated according to Eq. (16). This quantity is large at the PCTR. The absolute values of the drift momentum reach 1 − 2 × 10-8 kg m-2 s-1 which means about 10−15% of the values of the individual ion and neutral momenta. At the first instants of the instability, pD is positive. Combined with the negative ion and neutral momenta, positive values of pD mean that neutrals fall faster than the ions. The drift momentum remains positive at the bottom part of the drop at all instants. At t = 160 s, when the drop starts rising, both ion and neutral momenta are positive, and pD> 0 below the drop mean faster upward motion of ions. This behavior can be expected since the magnetic field does not act on neutrals directly, but only via the collisions with charged particles. In the limiting case of the absence of collisions, the instability in the neutral gas would develop entirely hydrodynamically, i.e., without cut-off wavelength. The collisional interaction with ions prevents this from happening and forces the collective behavior of the plasma. Nevertheless, the coupling becomes less strong at the PCTR, where strong gradients in all parameters, including the collisional frequency, are present. The non-ideal partial ionization effects on the RTI essentially originate from this narrow layer.

5. Discussion and conclusions

We have studied how the presence of neutral atoms in a prominence plasma influences the development of the Rayleigh-Taylor instability at the coronal-prominence interface by means of numerical 2.5D simulations. Our approach consisted in solving the single-fluid quasi-MHD equations including the physical diffusion term (ambipolar diffusion) in the induction and energy equations. This approach is justified in the regime of strong collisional coupling of the plasma. Our main goal was to investigate the RTI in partially ionized plasmas in the non-linear regime to verify and extend the conclusions from the linear theory.

Our main finding is that the configuration containing neutral atoms is always unstable on small scales. While in the completely ionized plasma the growth rate of the small-scale harmonics is lowered (or becomes zero) due to the action of magnetic forces, this does not occur if neutral atoms are present. We obtain an increase in the growth rate of the small-scale harmonics of up to 50% when partial ionization effects are taken into account (Fig. 12). This result is in good agreement with the linear theory (Soler et al. 2012; Díaz et al. 2012, 2014, see our Fig. 6). We show that by relaxing the approximations of the linear analysis, the growth rate is still large in the non-linear regime of the RTI. A statistical comparison reveals that this fast growth rate at small scales produces, on average, 2−3 per cent higher flow velocities in the model with the ambipolar diffusion term “on” compared to the purely MHD model (Fig. 11).

Another action of the ambipolar diffusion is the dissipation of currents in the direction perpendicular to the magnetic field. This dissipation allows the transformation of magnetic energy into internal energy (Khomenko & Collados 2012; Martínez-Sykora et al. 2012), and this results in the heating of the plasma in the prominence part of the domain. The temperature increase due to this effect is about 5−10% inside the prominence, and about 10−20% at the PCTR, compared to the model without ambipolar diffusion (Fig. 13). While the heating is produced, we observe that the amplitude of perpendicular currents, J⊥, becomes progressively lower in the model with ambipolar diffusion, reaching a 20−30% difference with the pure MHD case. The ambipolar diffusion introduces an anisotropy in the system since, with time, perpendicular currents dissipate and the longitudinal ones are unaltered. Such anisotropic dissipation tends to create force-free structures, as was shown in the simulations by Leake & Arber (2006) and Arber et al. (2009).

As the instability evolves, the initially sharp interface gets smoother and is filled with a mixture of coronal and prominence material. In this transition layer, the density, ionization fraction, collisional frequency, and other parameters vary strongly from prominence to coronal values. Because of the decrease in density, the collisional coupling becomes less strong, and as a consequence, a non-negligible drift momentum appears at this layer (Eq. (16), Fig. 14) with a value of 10−15% of the individual ion and neutral momenta. The sign of the drift momentum indicates that the neutral atoms at the bottom part of the downflowing drop move with slightly higher downward velocities. The neutrals feel the presence of the magnetic Lorentz force only through the collisions with the ionized particles and, once the collisional coupling weakens, this relative motion between the components becomes possible. Perpendicular currents also reach maximum values inside the transition layer between the prominence and coronal material. Therefore, the action of partial ionization effects on the RTI is localized in space to a small zone. This explains why the inclusion of these effects in our current modeling only alters the global parameters of the flow (thermodynamic parameters, velocities, magnetic field) by no more than a few tens of percent. In a different, non current-free equilibrium configuration, the action of partial ionization effects could be significantly amplified.

Our initial equilibrium configuration is the simplest possible, purely hydrodynamical with a homogeneous magnetic field. This configuration is different from the equilibrium usually thought to exist in prominences producing their large-scale stability (Tandberg-Hanssen 1995; Mackay et al. 2010). The support of the prominence material is thought to be provided by the magnetic tension. Contrarily, we neglect the effects of the magnetic field curvature in our analysis. Our choice was motivated by two factors. On the one hand, this equilibrium is a natural choice for the 2.5D modeling that does not allow us to include the effects of magnetic field curvature perpendicular to the perturbation plane (like the Kippenhahn-Schlüter model, see Hillier et al. 2012a,c). We have presented here an exploratory study to investigate the importance of partial ionization effects, and 2.5D simulations instead of full 3D allow for faster calculations. On the other hand, we have also pursued with this work an adequate comparison with the linear theory, where the same equilibrium is adopted (Díaz et al. 2014).

Our equilibrium model causes several limitations. Since the temperature in the prominence part of the domain is rather low, the pressure scale height is only a few hundred km, and we are unable to extend the model in height to have a large slab of prominence material, since the pressure and density quickly drop to coronal values. This limits the size of the computational box, making impossible the direct comparison to observations of the whole prominence structure. The current-free equilibrium is another drawback since, as already mentioned above, the action of the partial ionization effects is only limited to a narrow transition zone, decreasing its potential influence.

As a consequence of our equilibrium model, the force balance developed in the simulations is different from the one found in Hillier et al. (2012a,c). The main force to balance gravity in the downflowing drop in Fig. 8 is found to be the magnetic pressure force (Fig. 10), unlike the magnetic tension force in the aforementioned papers. Since we have not introduced any buoyant rising material, as in Hillier et al. (2012a,c), the distribution of upflows and downflows is significantly more symmetric in our case (Fig. 5). However, because of mass conservation, the upflowing rising bubbles of the coronal material have significantly larger sizes than the downflowing fingers of dense prominence material (Figs. 2 and 8). Therefore, assuming the density to be a proxy for the intensity in Hα imaging observations, observational detection of upflows is easier and could cause the observed asymmetry (Berger et al. 2010; Ryutova et al. 2010).

The simulations described above are only done for two orientations of the magnetic field, one normal to the perturbation plane, θ = 90°, and another one skewed by just one degree from the normal, θ = 89°. It might seem surprising that the structures developed in both cases are so different. Such behavior is nevertheless expected from the properties of the non-linear flow cascade and are observed in many other simulations of the RTI in different astrophysical contexts in 2D and 3D (Jun et al. 1995; Wang & Robertson 1985; Isobe et al. 2006; Stone & Gardiner 2007a,b; Hillier et al. 2012a). The increase of the dominant scale with time (Fig. 4) is in good agreement with other works (Jun et al. 1995). Since the maximum perturbation wavelength was of the size of the computational domain, λ = 250 km, it is natural that with time the perturbation of this scale dominates. The results of our initial calculations in a larger box (Khomenko et al. 2013) show that the same perturbation develops several drops on the same timescale. There, a simulation with θ = 88° was also analyzed and essentially lead to the same conclusions as for the enhanced growth rate of small-scale modes with the ambipolar diffusion term “on”. The disappearance of small scales in the θ = 89° simulation is caused by the cutoff wavelength introduced by the component of the magnetic field in the perturbation plane and is also in agreement with other numerical works (Stone & Gardiner 2007a). One also has to keep in mind that the size of our computational box is relatively small compared to the cutoff wavelength introduced by a magnetic field skewed by just 1° away from normal, λc = 38 km.

The main issue of our modeling is the use of the Saha equation to update the electron density and neutral fraction at each time step of the simulations. At prominence temperatures of the plasma, the deviations from the instantaneous ionization equilibrium can be significant, and the use of the Saha equation may lead to an underestimation of the electron density. A more appropriate approach would be to consider a time-dependent ionization balance. Leenaarts & Wedemeyer-Böhm (2006); Leenaarts et al. (2007) and Wedemeyer-Böhm & Carlsson (2011) have shown that, while the ionization process happens rapidly, the recombination is slow, so the ionization fraction is maintained almost constant in time and space even when significant temperature fluctuations are present. The calculations of the impact of time-dependent ionization balance on the development of the RTI in partially ionized plasma needs a thorough study. However, such calculations require significant computational resources and are beyond the scope of the present explorative work.

Besides the Rayleigh-Taylor instability in itself, the posterior time evolution of the downflowing drop from Fig. 8 deserves a separate discussion. After the downward motion of the drop stops at t ≈ 160 s, it continues to oscillate around the equilibrium position, extending in the horizontal direction. The evolution of the drop brings the magnetic field lines close together in the coronal part of the domain and reconnection happens at about 200 s of the simulation. The plasmoid formed after this reconnection is visible at time 237 s in Fig. 8. The layer of plasma linking the drop to the main part of the prominence becomes thinner and finally another reconnection happens in the chromospheric part of the domain. Another magnetic island is formed at the location of the drop, and the drop becomes almost completely disconnected from the rest of the prominence, extending even more in the horizontal direction and forming a thread (not shown in Fig. 8). The reconnections and the formation of the horizontal thread will need further analysis in a separate paper. The process of the drop falling across the horizontal magnetic field lines under the action of gravity, its disconnection from the magnetic field of the rest of the prominence, and its forming a finite plasma island may resemble the mechanism proposed by Haerendel & Berger (2011).

Summarizing all the above, we conclude that partial ionization effects on the Rayleigh-Taylor instability in prominences are non-negligible and have to be taken into account in models of prominence dynamics. In the future, we will consider larger simulation boxes in three dimensions to perform the comparison with observations, and will investigate the development of the instability for different ionization fractions and initial equilibrium configurations.

1

Because of their peculiar appearance, quiescent prominences with vertical threads are also known as hedgerow prominences.

2

See Appendix A for the discussion of the common nomenclature for this term used in the literature.

3

Here and below we take as PCTR the region from 780 km to 940 km, with the largest ambipolar diffusion term.

4

We do not show the corresponding MHD simulation in the same format since the developed structures are similar and the differences are much better appreciated after the Fourier analysis or just by subtracting the individual snapshots, see Figs. 11−13.

Acknowledgments

This work is partially supported by the Spanish Ministry of Science through projects AYA2010-18029 and AYA2011-24808. This work contributes to the deliverables identified in FP7 European Research Council grant agreement 277829, “Magnetic connectivity through the Solar Partially Ionized Atmosphere”.

References

  1. Antiochos, S. K., & Klimchuk, J. A. 1991, ApJ, 378, 372 [NASA ADS] [CrossRef] [Google Scholar]
  2. Antiochos, S. K., Dahlburg, R. B., & Klimchuk, J. A. 1994, ApJ, 420, L41 [Google Scholar]
  3. Arber, T. D., Haynes, M., & Leake, J. E. 2007, ApJ, 666, 541 [NASA ADS] [CrossRef] [Google Scholar]
  4. Arber, T. D., Botha, G. J. J., & Brady, C. S. 2009, ApJ, 705, 1183 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aulanier, G., & Demoulin, P. 1998, A&A, 329, 1125 [NASA ADS] [Google Scholar]
  6. Aulanier, G., DeVore, C. R., & Antiochos, S. K. 2006, ApJ, 646, 1349 [NASA ADS] [CrossRef] [Google Scholar]
  7. Babel, J., & Michaud, G. 1991, A&A, 248, 155 [NASA ADS] [Google Scholar]
  8. Bai, X., & Stone, J. M. 2011, ApJ, 736, 144 [NASA ADS] [CrossRef] [Google Scholar]
  9. Berger, T., Testa, P., Hillier, A., et al. 2011, Nature, 472, 197 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  10. Berger, T. E., Shine, R. A., Slater, G. L., et al. 2008, ApJ, 676, L89 [Google Scholar]
  11. Berger, T. E., Slater, G., Hurlburt, N., et al. 2010, ApJ, 716, 1288 [NASA ADS] [CrossRef] [Google Scholar]
  12. Black, D. C., & Scott, E. H. 1982, ApJ, 263, 696 [NASA ADS] [CrossRef] [Google Scholar]
  13. Boyd, T. J. M., & Sanderson, J. J. 2003, The Physics of Plasmas (Cambridge: Cambridge University Press) [Google Scholar]
  14. Braginskii, S. I. 1965, Reviews in plasma physics, 1, 205 [Google Scholar]
  15. Brandenburg, A., & Zweibel, E. H. 1994, ApJ, 427, L91 [NASA ADS] [CrossRef] [Google Scholar]
  16. Brandenburg, A., & Zweibel, E. H. 1995, ApJ, 448, 734 [NASA ADS] [CrossRef] [Google Scholar]
  17. Caunt, S. E., & Korpi, M. J. 2001, A&A, 369, 706 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Chandrasekhar, S. 1981, Hydrodynamic and Hydromagnetic stability (Oxford: Clarendon Press) [Google Scholar]
  19. Cheung, M. C. M., & Cameron, R. H. 2012, ApJ, 750, 6 [NASA ADS] [CrossRef] [Google Scholar]
  20. Chitre, S. M., & Krishan, V. 2001, MNRAS, 323, L23 [NASA ADS] [CrossRef] [Google Scholar]
  21. Ciolek, G. E., & Mouschovias, T. C. 1993, ApJ, 418, 774 [NASA ADS] [CrossRef] [Google Scholar]
  22. Cowling, T. G. 1956, MNRAS, 116, 114 [NASA ADS] [Google Scholar]
  23. DeVore, C. R., & Antiochos, S. K. 2000, ApJ, 539, 954 [NASA ADS] [CrossRef] [Google Scholar]
  24. Díaz, A. J., Soler, R., & Ballester, J. L. 2012, ApJ, 754, 41 [NASA ADS] [CrossRef] [Google Scholar]
  25. Díaz, A. J., Khomenko, E., & Collados, M. 2014, A&A, 564, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Dudík, J., Aulanier, G., Schmieder, B., Zapiór, M., & Heinzel, P. 2012, ApJ, 761, 9 [NASA ADS] [CrossRef] [Google Scholar]
  27. Dudorov, A. E. 1991, Soviet Astron. Lett., 17, 222 [NASA ADS] [Google Scholar]
  28. Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659 [NASA ADS] [CrossRef] [Google Scholar]
  29. Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [Google Scholar]
  30. Fiedler, R. A., & Mouschovias, T. C. 1992, ApJ, 391, 199 [NASA ADS] [CrossRef] [Google Scholar]
  31. Flower, D. R., & Pineau des Forêts, G. 2003, MNRAS, 341, 1272 [NASA ADS] [CrossRef] [Google Scholar]
  32. Fontenla, J. M., Avrett, E. H., & Loeser, R. 1990, ApJ, 355, 700 [NASA ADS] [CrossRef] [Google Scholar]
  33. Forteza, P., Oliver, R., & Ballester, J. L. 2007, A&A, 461, 731 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Gardner, C. L., Glimm, J., McBryan, O., et al. 1988, Phys. Fluids, 31, 447 [NASA ADS] [CrossRef] [Google Scholar]
  35. Gibson, S. E., & Fan, Y. 2006, J. Geophys. Res. 111, 12103 [Google Scholar]
  36. Gilbert, H., Kilper, G., & Alexander, D. 2007, ApJ, 671, 978 [NASA ADS] [CrossRef] [Google Scholar]
  37. Gilbert, H. R., Hansteen, V. H., & Holzer, T. E. 2002, ApJ, 577, 464 [NASA ADS] [CrossRef] [Google Scholar]
  38. Greaves, J. S., & Holland, W. S. 1999, MNRAS, 302, L45 [NASA ADS] [CrossRef] [Google Scholar]
  39. Haerendel, G., & Berger, T. 2011, ApJ, 731, 82 [NASA ADS] [CrossRef] [Google Scholar]
  40. Heinzel, P., Schmieder, B., Fárník, F., et al. 2008, ApJ, 686, 1383 [NASA ADS] [CrossRef] [Google Scholar]
  41. Heitsch, F., Zweibel, E. G., Slyz, A. D., & Devriendt, J. E. G. 2004, ApJ, 603, 165 [NASA ADS] [CrossRef] [Google Scholar]
  42. Hillier, A., Berger, T., Isobe, H., & Shibata, K. 2012a, ApJ, 746, 120 [NASA ADS] [CrossRef] [Google Scholar]
  43. Hillier, A., Hillier, R., & Tripathi, D. 2012b, ApJ, 761, 106 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hillier, A., Isobe, H., Shibata, K., & Berger, T. 2012c, ApJ, 756, 110 [NASA ADS] [CrossRef] [Google Scholar]
  45. Hiraki, Y., Krishan, V., & Masuda, S. 2010, ApJ, 720, 1311 [NASA ADS] [CrossRef] [Google Scholar]
  46. Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2005, Nature, 434, 478 [NASA ADS] [CrossRef] [Google Scholar]
  47. Isobe, H., Miyagoshi, T., Shibata, K., & Yokoyama, T. 2006, PASJ, 58, 423 [Google Scholar]
  48. Jones, A. C., & Downes, T. P. 2011, MNRAS, 418, 390 [NASA ADS] [CrossRef] [Google Scholar]
  49. Jun, B.-I., & Norman, M. L. 1995, Ap&SS, 233, 267 [NASA ADS] [CrossRef] [Google Scholar]
  50. Jun, B.-I., Norman, M. L., & Stone, J. M. 1995, ApJ, 453, 332 [NASA ADS] [CrossRef] [Google Scholar]
  51. Khodachenko, M. L., Arber, T. D., Rucker, H. O., & Hanslmeier, A. 2004, A&A, 422, 1073 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Khodachenko, M., Rucker, H., Oliver, R., Arber, T., & Hanslmeier, A. 2006, Adv. Space Res., 37, 447 [NASA ADS] [CrossRef] [Google Scholar]
  53. Khomenko, E., & Collados, M. 2012, ApJ, 747, 87 [NASA ADS] [CrossRef] [Google Scholar]
  54. Khomenko, E., Collados, M., & Felipe, T. 2008, Sol. Phys., 251, 589 [NASA ADS] [CrossRef] [Google Scholar]
  55. Khomenko, E., Díaz, A. D., de Vicente, A., Collados, M., & Luna, M. 2013, in IAUS300: Nature of solar prominences and their role in Space Weather, eds. J. M. Malherbe, & B. Schmieder (Cambridge University Press), in press [Google Scholar]
  56. Krall, N. A., & Trivelpiece, A. W. 1973, Principles of plasma physics, International Student Edition − International Series in Pure and Applied Physics (Tokyo: McGraw-Hill Kogakusha) [Google Scholar]
  57. Krasnoselskikh, V., Vekstein, G., Hudson, H. S., Bale, S. D., & Abbett, W. P. 2010, ApJ, 724, 1542 [NASA ADS] [CrossRef] [Google Scholar]
  58. Leake, J. E., & Arber, T. D. 2006, A&A, 450, 805 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Leake, J. E., & Linton, M. G. 2013, ApJ, 764, 54 [NASA ADS] [CrossRef] [Google Scholar]
  60. Leake, J. E., DeVore, C. R., Thayer, J. P., et al. 2013 [arXiv:1310.0405] [Google Scholar]
  61. Leenaarts, J., & Wedemeyer-Böhm, S. 2006, A&A, 460, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Leroy, J. L. 1989, in Dynamics and Structure of Quiescent Solar Prominences, Astrophys. Space Sci. Lib., ed. E. R. Priest, 150, 77 [Google Scholar]
  64. Li, P. S., McKee, C. F., & Klein, R. I. 2006, ApJ, 653, 1280 [Google Scholar]
  65. Mac Low, M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726 [NASA ADS] [CrossRef] [Google Scholar]
  66. Mackay, D. H., Karpen, J. T., Ballester, J. L., Schmieder, B., & Aulanier, G. 2010, Space Sci. Rev., 151, 333 [NASA ADS] [CrossRef] [Google Scholar]
  67. Martínez-Sykora, J., De Pontieu, B., & Hansteen, V. 2012, ApJ, 753, 161 [NASA ADS] [CrossRef] [Google Scholar]
  68. Mestel, L., & Spitzer, Jr., L. 1956, MNRAS, 116, 503 [Google Scholar]
  69. Mouschovias, T. C. 1977, ApJ, 211, 147 [NASA ADS] [CrossRef] [Google Scholar]
  70. Mouschovias, T. C., & Paleologou, E. V. 1981, ApJ, 246, 48 [NASA ADS] [CrossRef] [Google Scholar]
  71. Myers, P. C., & Goodman, A. A. 1988, ApJ, 329, 392 [NASA ADS] [CrossRef] [Google Scholar]
  72. Nakano, T. 1977, PASJ, 29, 197 [NASA ADS] [Google Scholar]
  73. Oishi, J. S., & Mac Low, M.-M. 2006, ApJ, 638, 281 [NASA ADS] [CrossRef] [Google Scholar]
  74. Osterbrock, D. E. 1961, ApJ, 134, 270 [Google Scholar]
  75. Paleologou, E. V., & Mouschovias, T. C. 1983, ApJ, 275, 838 [NASA ADS] [CrossRef] [Google Scholar]
  76. Pandey, B. P., & Wardle, M. 2008, MNRAS, 385, 2269 [NASA ADS] [CrossRef] [Google Scholar]
  77. Pandey, B. P., Vranjes, J., & Krishan, V. 2008, MNRAS, 386, 1635 [NASA ADS] [CrossRef] [Google Scholar]
  78. Pneuman, G. W. 1983, Sol. Phys., 88, 219 [NASA ADS] [CrossRef] [Google Scholar]
  79. Priest, E. R. 1982, Solar magneto-hydrodynamics (Dordrecht: D. Reidel) [Google Scholar]
  80. Priest, E. R., Hood, A. W., & Anzer, U. 1989, ApJ, 344, 1010 [NASA ADS] [CrossRef] [Google Scholar]
  81. Roberge, W. G., Hanany, S., & Messinger, D. W. 1995, ApJ, 453, 238 [NASA ADS] [CrossRef] [Google Scholar]
  82. Rust, D. M., & Kumar, A. 1994, Sol. Phys., 155, 69 [NASA ADS] [CrossRef] [Google Scholar]
  83. Ryutova, M., Berger, T., Frank, Z., Tarbell, T., & Title, A. 2010, Sol. Phys., 267, 75 [Google Scholar]
  84. Sakai, J. I., & Smith, P. D. 2009, ApJ, 691, L45 [NASA ADS] [CrossRef] [Google Scholar]
  85. Sakai, J. I., Tsuchimoto, K., & Sokolov, I. V. 2006, ApJ, 642, 1236 [NASA ADS] [CrossRef] [Google Scholar]
  86. Scott, E. H. 1983, ApJ, 275, 836 [NASA ADS] [CrossRef] [Google Scholar]
  87. Singh, K. A. P., Shibata, K., Nishizuka, N., & Isobe, H. 2011, Phys. Plasmas, 18, 1210 [Google Scholar]
  88. Smith, M. D., & Mac Low, M.-M. 1998, Ap&SS, 261, 191 [NASA ADS] [CrossRef] [Google Scholar]
  89. Soler, R., Oliver, R., & Ballester, J. L. 2009, ApJ, 699, 1553 [NASA ADS] [CrossRef] [Google Scholar]
  90. Soler, R., Oliver, R., & Ballester, J. L. 2010, A&A, 512, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Soler, R., Díaz, A. J., Ballester, J. L., & Goossens, M. 2012, ApJ, 749, 163 [NASA ADS] [CrossRef] [Google Scholar]
  92. Spitzer, L. J. 1968, Diffuse matter in space (New York: Interscience) [Google Scholar]
  93. Spitzer, L. J. 1978, Physical processes in the interstellar medium (New York: Interscience) [Google Scholar]
  94. Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
  95. Stone, J. M., & Gardiner, T. 2007a, Physics of Fluids, 19, 094104 [Google Scholar]
  96. Stone, J. M., & Gardiner, T. 2007b, ApJ, 671, 1726 [NASA ADS] [CrossRef] [Google Scholar]
  97. Stone, J. M., & Zweibel, E. G. 2010, ApJ, 724, 131 [NASA ADS] [CrossRef] [Google Scholar]
  98. Tandberg-Hanssen, E. 1995, The nature of solar prominences, in Astrophys. Space Sci. Lib. 199 [Google Scholar]
  99. Tsap, Y. T., Stepanov, A. V., & Kopylova, Y. G. 2012, Astron. Rep., 56, 138 [NASA ADS] [CrossRef] [Google Scholar]
  100. van Ballegooijen, A. A., & Cranmer, S. R. 2010, ApJ, 711, 164 [NASA ADS] [CrossRef] [Google Scholar]
  101. van Ballegooijen, A. A., & Martens, P. C. H. 1989, ApJ, 343, 971 [NASA ADS] [CrossRef] [Google Scholar]
  102. Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Wang, Y.-M., & Robertson, J. A. 1985, ApJ, 299, 85 [NASA ADS] [CrossRef] [Google Scholar]
  104. Wedemeyer-Böhm, S., & Carlsson, M. 2011, A&A, 528, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Youngs, D. L. 1984, Phys. D Nonlinear Phenomena, 12, 32 [NASA ADS] [CrossRef] [Google Scholar]
  106. Zweibel, E. G. 1988, ApJ, 329, 384 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Notes on the nomenclature

The term proportional to ηAJ⊥ in Eqs. (2) describes the coupling of neutral particles to the charges. It causes the magnetic field to diffuse through neutral gas because of the collisions between neutrals and charged particles. In an astrophysical context, especially in the fields of interstellar medium, star formation, and stellar atmosphere modeling, this process is usually referred to as “ambipolar diffusion”. This notation was first used by Spitzer in his books “Diffuse matter in space” (1968) and “Physical processes in the interstellar medium” (1978). Ambipolar diffusion makes star formation possible from gravitational collapse in magnetized protostellar clouds as was suggested in the pioneering work by Mestel & Spitzer (1956). In the subsequent years this terminology was adopted as standard in this field. An incomplete list of works where the decoupling of neutral particles from plasma is called “ambipolar diffusion” includes Osterbrock (1961); Nakano (1977); Mouschovias (1977); Mouschovias & Paleologou (1981); Black & Scott (1982); Scott (1983); Paleologou & Mouschovias (1983); Zweibel (1988); Myers & Goodman (1988); Fontenla et al. (1990); Dudorov (1991); Babel & Michaud (1991); Fiedler & Mouschovias (1992); Ciolek & Mouschovias (1993); Brandenburg & Zweibel (1994, 1995); Mac Low et al. (1995); Roberge et al. (1995); Smith & Mac Low (1998); Greaves & Holland (1999); Flower & Pineau des Forêts (2003); Heitsch et al. (2004); Oishi & Mac Low (2006); Li et al. (2006); Duffin & Pudritz (2008); Stone & Zweibel (2010); Bai & Stone (2011), and Jones & Downes (2011).

On the contrary, in the field of plasma physics, the term “ambipolar diffusion” refers to the diffusion of positive and negative particles at the same rate due to their Coulomb interaction, which maintains the charge neutrality at scales larger than the Debye length (see the standard plasma physics books by Boyd & Sanderson 2003; Krall & Trivelpiece 1973). The use of the same terminology applied to these two different phenomena may lead to confusion.

In the context of solar physics, the additional diffusion due to neutrals is frequently referred to as “Cowling resistivity”, following the notation by Braginskii (1965; see the original work by Cowling 1956). Strictly speaking, the Cowling resistivity includes both Ohmic and neutral diffusive terms, i.e., ηC = ηA + η; however, the Ohmic term is orders of magnitude lower in the solar atmosphere and is usually neglected. Recent citations where the Cowling resistivity terminology is used include, e.g., Khodachenko et al. (2004, 2006); Leake & Arber (2006); Forteza et al. (2007); Arber et al. (2007, 2009), and Sakai & Smith (2009). Alternatively, this term is also called “Pedersen resistivity” in works by Sakai et al. (2006); Krasnoselskikh et al. (2010); Leake & Linton (2013); Leake et al. (2013); Cheung & Cameron (2012), and others. Nevertheless, many other works in solar physics use the terminology introduced by Spitzer, i.e., “ambipolar diffusion” (see Chitre & Krishan 2001; Pandey et al. 2008; Soler et al. 2009, 2010; Hiraki et al. 2010; Singh et al. 2011; Cheung & Cameron 2012; Martínez-Sykora et al. 2012; Soler et al. 2012; Tsap et al. 2012; Díaz et al. 2012; Khomenko & Collados 2012).

In this paper we use ambipolar diffusion in the astrophysical sense, as a widely spread nomenclature in this research field.

All Tables

Table 1

Parameters of the equilibrium configuration.

In the text

All Figures

thumbnail Fig. 1

Initial configuration. The instability develops in the XZ plane containing the perturbation vector k; the initial magnetic field B0 forms an angle θ with the XZ plane. B is initially lying in the XY plane, parallel to the interface separating the prominence and coronal plasma, with Bz = 0, Bx = B0cos(θ), and By = B0sin(θ).
In the text
thumbnail Fig. 2

Time evolution of density (top) and pressure (bottom) in the AD simulation with θ = 90°. The size of each snapshot is 250 × 1000 km, the time elapsed is given at the bottom of each panel. The velocity field is indicated by arrows. We note the asymmetry between the large-scale rising bubbles and small-scale downflowing fingers in the density images.
In the text
thumbnail Fig. 3

Same as Fig. 2 but for the MHD simulation.
In the text
thumbnail Fig. 4

Mode number (m, see Eq. (9)) with maximum Fourier amplitude of the pressure perturbation in the AD (filled circles) and MHD (open circles) simulations with θ = 90°. The wavelength of the modes varies from 10 km (mode 25) to 250 km (mode 1). The amplitude of the modes is indicated by color from blue (smaller) to red (larger). Dotted line is an exponential function ~exp( − t/ 16) shown for illustration purposes.
In the text
thumbnail Fig. 5

Velocity histograms over the whole simulation time for the θ = 90° MHD and AD runs. We note the asymmetry of the histogram and the slight tendency of the AD simulation to develop higher velocities.
In the text
thumbnail Fig. 6

Growth rate of the RTI from the linear theory according to Díaz et al. (2014) as a function of the wave number of the perturbation along the discontinuity. The parameters of the atmosphere assumed for the linear calculation are those from Table 1; the magnetic field is inclined with respect to k by θ = 89° (dashed line) and by θ = 90° (dotted line). Solid line shows the case of θ = 89° with the ambipolar diffusion term taken into account.
In the text
thumbnail Fig. 7

Fourier analysis of scales developed in the θ = 90° simulations. Panels on the left show the relative power in the AD simulation divided by the power in MHD simulation as a function of horizontal wave number along the discontinuity (k) and time. Red means that the AD simulation has more power. Panels on the right give the time average of the power ratio from the panels on the left.
In the text
thumbnail Fig. 8

Time evolution of density (top) and pressure (bottom) in the AD simulation with θ = 89°. The size of each snapshot is 250 × 1000 km; the time elapsed is given at the bottom of each panel. Green lines are projections of the magnetic field lines into the XZ plane; the velocity field is indicated by arrows. We note the difference in scales compared to the θ = 90° case.
In the text
thumbnail Fig. 9

Three-dimensional representation of the magnetic field lines in a snapshot of the AD simulation with θ = 89° at t = 160 s. The vertical Z axis is 1000 km long. The density of the field lines is proportional to the magnetic field strength and their origins in XZ plane is random. Projections of the field lines into XZ and YZ planes are indicated in green color.
In the text
thumbnail Fig. 10

Top left: acceleration as a function of time due to the gravity force (dash-dotted red line), gas pressure gradient force (triple-dot-dashed green line), magnetic pressure gradient force (solid black line), and magnetic tension force (dashed black line). The values are integrated over the big drop from Fig. 8 (AD run) after it forms at about 70 s of the simulation time. Bottom left: acceleration due to the sum of all above forces. Top right: vertical velocity of the big drop. Vertical dotted lines indicate time instants of zero velocity. Bottom right: effective cutoff wavelength calculated from Eq. (1) for the horizontal field component Bcos(θ) at the bottom of the drop. Horizontal dotted lines indicate the size of the computational domain L = 250 km.
In the text
thumbnail Fig. 11

Left: difference between the standard deviations of the vertical velocity near the discontinuity in the AD and MHD θ = 89° simulations during the first 120 s. Right: histogram of the difference between the absolute values of the vertical velocity near the discontinuity in the AD and MHD θ = 89° simulations during the first 120 s.
In the text
thumbnail Fig. 12

Fourier analysis of scales developed in the θ = 89° simulations. The format of the figure is the same as Fig. 7. The vertical dotted lines in the right panels mark the cutoff wavelength λc = 38 km.
In the text
thumbnail Fig. 13

Difference between the temperature variations in the AD and MHD simulations with θ = 89°, (TA − TMHD) /TMHD. The size of each snapshot is 250 × 1000 km; the time elapsed is given at the bottom of each panel. Red colors mean that the AD model is hotter.
In the text
thumbnail Fig. 14

Top: time evolution of the vertical component of the neutral momentum from Eq. (14) in the θ = 89° AD run. Negative values mean downward motion. The contours of vertical drift momentum and −10-9 kg m-2 s-1 (see Eq. (16)) are overplotted in green and gray, respectively. Bottom: the same for the ion momentum. The size of each snapshot is 250 × 1000 km; the time elapsed is given at the bottom of each panel, the color scheme is constant in time, and is the same in both panels.
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.