© ESO, 2012

1. Introduction

Thermal or condensational modes have been extensively investigated in magnetized plasmas (e.g., Parker 1953; Field 1965; Heyvaerts 1974). As explained by Parker (1953), thermal instability can occur in a diffuse medium due to imbalance between temperature-independent energy gains, i.e., heating and temperature-dependent radiative losses. Parker (1953) arrived at the qualitative criterion that instability can be present when radiative losses decrease as the temperature increases. This mechanism is important in the context of prominences since it allows the formation of cool plasma condensations in a medium of high temperature. Therefore, unstable thermal modes may play an important role in the formation of solar prominences and in the evolution of the prominence plasma. Field (1965) and Heyvaerts (1974) investigated the phenomenon in more detail than Parker (1953) and derived more accurate instability criteria. Subsequent papers that investigate both linear and nonlinear thermal instabilities, mainly in the field of prominences, include Hildner (1974), Oran et al. (1982), Dahlburg & Mariska (1988), Karpen et al. (1989), Cargill & Hood (1989), Van der Linden & Goossens (1991), Carbonell et al. (2004), and Soler et al. (2011), among others.

As shown in the instability criterion derived by Field (1965), an accurate description of the radiative loss function is crucial to ascertaining the stability of thermal modes. However, determination of the radiative loss function in prominence plasmas depending on the values of temperature and density is a difficult task that requires complicated numerical solutions of the radiative transfer equations for nonlocal thermodynamic equilibrium. Alternatively, several semi-empirical parametrizations of the radiative loss function for prominence and coronal conditions are available in the literature. These parametrizations enable us to incorporate radiative losses in a consistent way in the theoretical models of prominence plasmas without the need of solving the radiative transfer problem.

One of the first parametrizations of this kind can be found in Hildner (1974), who performed a piecewise fit as a function of the temperature of the computations of radiative losses available by then. Subsequent authors have proposed different parametrizations that update Hildner’s fit (e.g., Rosner et al. 1978; Milne et al. 1979), although Hildner’s function is still used in some works nowadays. Another function frequently used in the literature is the so-called Klimchuk-Raymond function (see, e.g., Klimchuk & Cargill 2001) that is a better representation of the radiative losses for prominence-corona transition region and coronal temperatures. The shape of the loss function depends on the completeness of the atomic model used for the calculation, on the atomic processes included, on the ionization equilibrium, and element abundance assumed. More recent loss functions that incorporate accurate atomic data information are, e.g., the loss function used by Parenti et al. (2006) and Parenti & Vial (2007), computed from the CHIANTI database (Dere et al. 1997; Landi et al. 2012), and the loss function computed by Schure et al. (2009) using the SPEX package (Kaanstra et al. 1996).

The purpose of this paper is to investigate the stability properties of thermal modes in cool prominence plasmas. We compare the results using up-to-date computations of radiative losses derived from the CHIANTI v7 atomic database with results that assumed two of the most used loss functions in the literature. These radiative losses are obtained by assuming an optically thin plasma, while the core of prominence cannot completely satisfy this condition. We will discuss our assumption later on. This paper is organized as follows. Section 2 contains a description of the equilibrium and the basic equations. The instability criterion for thermal modes is derived in Sect. 3. Parametric studies of the linear growth rate are done in Sect. 4. Finally, Sect. 5 contains the summary and discussion of the results.

2. Basic equations

2.1. Equilibrium

Our equilibrium configuration is a uniform plasma of infinite extent. We assume a partially ionized hydrogen plasma composed of ions, electrons, and neutrals. We use Cartesian coordinates, and all quantities are expressed in MKS units thorough this paper. The magnetic field, B, is uniform and oriented along the z-direction, namely B = Bê_z, with B constant. We denote the equilibrium mass density, temperature, and gas pressure by ρ, T, and p, respectively. The set of basic nonadiabatic MHD equations governing the plasma dynamics in the single-fluid approximation are (see, e.g., Braginskii 1965) where is the material derivative for time variations following the plasma motion, v the plasma velocity, μ the magnetic permittivity, γ the adiabatic index, κ the thermal conductivity tensor, L(ρ,T) the heat-loss function, η and η_C are the coefficients of Ohm’s and Cowling’s diffusion, respectively, R is the ideal gas constant, and is the mean atomic weight, with ξ_i the relative fraction of ions. This parameter ranges from ξ_i = 1 in a fully ionized plasma and ξ_i = 0 in a neutral gas. In Eqs. (1)–(5), the effects of gravity and viscosity have been omitted. Equation (3) is the induction equation. In Eq. (3) we have neglected some minor terms that are several orders of magnitude smaller than Ohm’s and Cowling’s terms in partially ionized prominence plasmas (see an expression for the complete single-fluid induction equation in, e.g., Forteza et al. 2007; Zaqarashvili et al. 2011). The neglected terms are Hall’s term, the diamagnetic current term, and Biermann’s battery term. Soler et al. (2009b) showed that the effect of Hall’s term on the waves is negligible in prominence conditions. The diamagnetic current term and Biermann’s battery term are only relevant when large pressure gradients are present, a situation that is more representative of stellar interiors. The omission of these terms in the present work is therefore justified. In a partially ionized plasma, Cowling’s diffusion represents an enhanced magnetic diffusion due to ion-neutral collisions. The expressions for η and η_C are where m_e is the electron mass, e the electron charge, lnΛ_C Coulomb’s logarithm (see, e.g., Priest 1982), and α_n is the neutral friction coefficient, given by (8)with m_p the proton mass and k_B Boltzmann’s constant. In the fully ionized case, η_C = η and the third term on the right-hand side of Eq. (3) is absent.

We denote by κ_∥ and κ_⊥ the parallel and perpendicular scalar components of the thermal conductivity tensor with respect to the magnetic field direction, which can be expressed as (9)with κ_e, κ_i, and κ_n the conductivities by electrons, ions, and neutrals, respectively. In a fully ionized medium, κ_∥ is governed by electrons, whereas κ_⊥ is caused by ions. In a partially ionized plasma, the contribution of neutrals, κ_n, has to be added to both scalar conductivities because thermal conduction by neutrals is isotropic. Since κ_i ≪ κ_n, the conductivity by ions can be neglected in κ_⊥. We use (see, e.g., Soler et al. 2010) (10)Regarding the heat-loss function, L(ρ,T), we consider the following expression, (11)where χ^∗ and α are functions of the temperature, and h is an arbitrary heating function. In the equilibrium we assume that radiative losses balance heating, so that h is defined to satisfy L(ρ,T) = 0 in the equilibrium state. Several parametrizations of χ^∗ and α for prominence conditions are available in the literature. In this work, we use two of the most familiar loss functions, namely the well-known parametrizations by Hildner (Hildner 1974) and Klimchuk-Raymond (Klimchuk & Cargill 2001). In addition, we use a parametrization of the radiative loss function computed from the CHIANTI v7 atomic database (Landi et al. 2012) assuming coronal abundances (Mazzotta et al. 1998), ionization equilibrium, and a constant pressure of 6.64 mPa. The parameters χ^∗ and α corresponding to these three different fits are given in Table 1. Figure 1 displays the loss rate per unit volume as a function of the temperature for the three fits. The main differences in the three curves are found at low (T ≲ 3  ×  10⁴ K) and high (T ≳ 10⁶ K) temperatures. For the region of our interest, which is the low temperature range, the peak at about 1.5  ×  10⁴ K in the CHIANTI calculation is mainly due to the H and He emissions. Besides this difference, it should be noted that the CHIANTI database is still incomplete at these temperatures, so that further increase in the loss function is expected for more complete calculations.

Fig. 1 Fit of the energy loss rate per unit volume as function of the temperature computed from the CHIANTI atomic database assuming coronal abundances, ionization equilibrium, and a constant pressure of 6.64 mPa (solid line) compared to Hildner’s (dotted line) and Klimchuk-Raymond’s fits (dashed line). The parameters χ∗ and α of the three fits are given in Table 1.

2.2. Dispersion relation for linear perturbations

We take the plasma initially at rest and superimpose linear perturbations on the equilibrium state. Equations (1)–(5) are linearized. We write perturbations proportional to exp(ik·r − iωt), where r = (x,y,z) is the position vector, the wavenumber vector, and ω the frequency. For simplicity, we choose the reference frame so that we can set k_y = 0 and consider wave propagation in the xz-plane only. We focus our study on nonadiabatic magnetoacoustic and thermal modes. Alfvén modes are not discussed in the present investigation. We combine the linearized Eqs. (1)–(5) and, after some algebraic manipulations, we arrive at the dispersion relation for nonadiabatic magnetoacoustic and thermal modes (see details in, e.g., Carbonell et al. 2004; Soler et al. 2010; Barceló et al. 2011), which can be written in a compact form as (12)where is the square of the wavevector modulus, cosθ = k_z/k is the cosine of the angle between k and B, and is the square of the nonadiabatic sound speed (see Soler et al. 2008), with the effective nonadiabatic index defined as (13)with (14)(15)In addition, Γ² is the square of the modified Alfvén speed (see Soler et al. 2009a) defined as (16)with the square of the ideal Alfvén speed. Both Λ² and Γ² are functions of ω and k. In the absence of nonadiabatic effects, and , so that Λ² becomes the square of the adiabatic sound speed, . In the absence of magnetic diffusion, η_C = η = 0 and Γ² becomes the square of the ideal Alfvén speed, . In such a case, Eq. (12) reverts to the well-known dispersion relation for ideal, adiabatic magnetoacoustic waves in a plasma of infinite extent (see, e.g., Lighthill 1960).

For positive and real k and θ, Eq. (12) has five solutions of ω. In general, ω is complex, namely ω = ω_R + iω_I, with ω_R and ω_I the real and imaginary parts of ω, respectively. The solutions of Eq. (12) are discussed in detail by Carbonell et al. (2004) and Barceló et al. (2011). Of the five solutions of ω, two complex conjugate solutions correspond to damped slow modes, and the two other complex conjugate solutions correspond to damped fast modes. The remaining solution is purely imaginary, i.e., ω_R = 0, and corresponds to the thermal mode. The thermal mode is the subject of our investigation.

3. Approximate instability criterion

We perform a first-order expansion of Eq. (12) for a low-β plasma, where β is the ratio of the gas pressure to the magnetic pressure, and obtain two different approximate dispersion relations, namely (17)for fast modes, and (18)for slow and thermal modes. Fast modes are weakly affected by nonadiabatic mechanisms and are not investigated further in the present work. We focus on Eq. (18) and explore the approximation for thermal modes in more detail. In the low-β regime the propagation of slow and thermal modes is almost parallel to the magnetic field. By assuming Γ² ≫ Λ², i.e., the Alfvén speed is much greater than the sound speed as is typical of prominence conditions, Eq. (18) simplifies to (19)We then use the definitions of Λ² and Γ² and expand Eq. (19) as a third-order polynomial in ω. Since the thermal mode is a purely imaginary solution, we write ω = ω_R + is, where s is the thermal mode growth rate, and set ω_R = 0. Because of the temporal dependence exp( − iωt), thermal mode perturbations are proportional to exp(st), meaning that perturbations grow in time for s > 0. This behavior corresponds to instability. In contrast, for s < 0 thermal disturbances are damped. To obtain an approximation for s, we neglect terms with in the polynomial expansion of Eq. (19). The neglected terms are related to the slow mode and produce minor corrections to the thermal mode growth rate. After some algebraic manipulations, we obtain the approximate growth rate, namely (20)In the absence of magnetic diffusion, η_C = 0 and Eq. (20) becomes (21)Equation (21) is similar to the expressions found by Van der Linden & Goossens (1991) and Soler et al. (2011) for thermal continuum modes, and by Carbonell et al. (2009) and Soler (2010) for the imaginary part of the frequency of propagating thermal waves in a flowing medium. In the case without thermal conduction, i.e., κ_∥ = κ_⊥ = 0, the approximate growth rate is (22)Equation (22) is independent of k, θ, and the magnetic field strength and orientation.

Fig. 2 Approximate thermal mode growth rate vs. temperature in the absence of thermal conduction and for Hildner’s (dotted), Klimchuk-Raymond’s (dashed), and CHIANTI-based (solid) loss functions. The shaded area denotes the region of instability at low temperatures obtained with the CHIANTI-based loss function that is not present for the other parametrizations.

Fig. 3 a) Numerically computed growth rate (solid) and approximate value obtained from Eq. (21) (dashed) vs. wavelength for the CHIANTI-based loss function and different values of the temperature in the case of a fully ionized plasma, i.e., ξi = 1. b) Same as panel a) but for T = 16   000 K and different values of ξi.

Equation (20) provides us with the instability criterion. By taking the definitions of , , ω_T, and ω_ρ (Eqs. (14) and (15)) into account, we find that the combination of parameters to have s > 0 in Eq. (20) must satisfy the condition (23)To the best of our knowledge, Eq. (23) is the first instability criterion for thermal modes that accounts for the effect of Cowling’s diffusion. We deduce from Eq. (23) that Cowling’s diffusion plays a stabilizing role since Cowling’s diffusion incorporates a lower bound in the instability criterion. To compare Eq. (23) with previous instability criteria derived in the literature, we set η_C = 0, and Eq. (23) simplifies to (24)Equation (24) agrees with the instability criterion provided by Field (1965) in his Eq. (25a) if the term accounting for perpendicular thermal conduction is added to Field’s expression. For simplicity, we take the case η_C = 0 and use the expression of the heat-loss function L (Eq. (11)) to rewrite Eq. (24) in terms of parameters χ ^∗  and α. Then, instability is present for values of α satisfying (25)In the absence of thermal conduction (κ _∥  = κ _⊥  = 0) or for long wavelengths (k → 0), the second term on the right-hand side of Eq. (25) vanishes and the instability criterion reduces to α < 1. Taking the values of the parameter α given in Table 1 into account, the lowest thermally unstable temperatures according to this criterion are T ≈ 8 × 10⁴ K in Hildner’s fit, T ≈ 9.33 × 10⁴ K in Klimchuk-Raymond’s fit, and T ≈ 1.58 × 10⁴ K in the CHIANTI-based fit. Important is that we find that the threshold temperature for the thermal instability is substantially reduced using the CHIANTI-based radiative losses in comparison to Hildner’s and Klimchuk-Raymond’s functions.

We do a simple computation without thermal conduction (κ_∥ = κ _⊥  = 0) and without Cowling’s diffusion (η_C = 0). In this case, the approximate growth rate is given by Eq. (22). Figure 2 displays the growth rate computed from Eq. (22) as a function of the temperature. A constant gas pressure of 6.64 mPa is assumed and the density is computed accordingly. As predicted by the instability criterion, there is a region of instability at low temperatures obtained with the CHIANTI-based loss function that is not present for the other parametrizations (see Fig. 2). This region of instability is present for temperatures between 1.58  ×  10⁴ K and 3.16  ×  10⁴ K. At these temperatures the prominence plasma is only partially ionized (see, e.g., Gouttebroze & Labrosse 2009) and the roles of thermal conduction by neutrals and Cowling’s diffusion may be relevant. These effects are investigated numerically in Sect. 4.

Additionally, in Fig. 2 we notice the abrupt jumps in the growth rate at the boundaries of the regions where different values of χ^∗ and α are used. The reason for these jumps is that, although the cooling function L(ρ,T) is continuous, the derivatives with respect to density and temperature are discontinuous where the values of parameters χ^∗ and α change. Additional comments on this issue are given in Van der Linden et al. (1991) and Soler et al. (2011).

4. Numerical results

Here we compute the thermal mode growth rate by solving the general dispersion relation (Eq. (12)) by standard numerical methods. We focus on the region of instability at low temperatures obtained with the CHIANTI-based loss function (Fig. 2). In the following computations we use B = 10 G, p = 6.64 mPa, and θ = π/4. These values of the equilibrium magnetic field strength and pressure are typical prominence parameters (see, e.g., Labrosse et al. 2010). We compute the growth rate, s, as a function of the perturbation wavelength, λ = 2π/k. Wavelengths typically observed in prominences are roughly between 10³ km and 10⁵ km (see Oliver & Ballester 2002). This range of wavelengths corresponds to disturbances usually detected from Doppler time series. We do not claim that all the observed disturbances are thermal modes. It is very difficult to determine the nature of the waves in the absence of additional information, such as the velocity polarization or the magnetic and density perturbations. In particular, a clear distinction between slow and thermal modes may be very complicated (see a discussion on this issue in Carbonell et al. 2009). In this paper we use the observed wavelengths given by Oliver & Ballester (2002) as the most probable range of wavelengths for thermal modes.

4.1. Fully ionized plasma

First, we investigate the fully ionized case and set ξ_i = 1. Therefore, the effect of thermal conduction by neutrals is absent, and Cowling’s diffusion becomes classical Ohm’s diffusion, i.e., η_C = η. Figure 3a displays s versus λ. Since we focus on the unstable behavior of the thermal mode, only positive values of s are displayed in Fig. 3a. We consider three different temperatures within the region of instability in Fig. 2. We compare the numerically computed growth rate with the approximation of Eq. (21) in the absence of diffusion, i.e., η = 0. The approximate result is in good agreement with the actual growth rate. Different test computations with and without Ohm’s magnetic diffusion (not displayed here for simplicity) indicate that magnetic diffusion has almost no impact on the value of the growth rate. This means that we can safely neglect the effect of diffusion and use the approximate growth rate given in Eq. (21) for the case η = 0.

Figure 3a also shows that the thermal mode is stabilized for short wavelengths. The stabilization is due to thermal conduction by electrons. Thermal conduction becomes important as the wavelength decreases, so that the growth rate is reduced, and the instability is suppressed for short enough wavelengths. For long wavelengths the growth rate saturates and becomes independent of λ.

Regarding the temperature, we find that the growth rate decreases as the temperature increases within the region of instability. For T = 16   000 K the maximum growth rate is s ≈ 0.6 s^-1, while for T = 30   000 K the maximum growth rate decreases to s ≈ 0.05 s^-1. These values of the linear growth rate indicate that the thermal instability operates on short timescales and suggest that the effect of the instability in prominences may be observable on the dynamics and evolution of cool plasma condensations. Nonlinear studies beyond the present linear analysis are needed in order to assess the actual impact of the instability on the condensation dynamics.

4.2. Partially ionized plasma

Here we incorporate the effects of partial ionization, namely thermal conduction by neutrals and Cowling’s diffusion. In these computations we fix the temperature to T = 16   000 K. Figure 3b displays s versus λ for different values of ξ_i. First of all, we find that, when ξ_i decreases, the critical wavelength for stabilization increases due to thermal conduction by neutrals. At cool prominence temperatures, thermal conduction by neutrals is more efficient than conduction by electrons. Thus, the critical wavelength for stabilization is now determined by the conductivity of neutrals. We also see in Fig. 3b that the growth rate increases when ξ_i is reduced. The growth rate for ξ_i = 0.1 is about four times higher, approximately, than for ξ_i = 1. We can explain these results by taking into account that a constant gas pressure is assumed in the computations, so when ξ_i decreases, the effective plasma density grows because of the increase in the amount of neutrals. The increase in the growth rate is a consequence of the increase in the effective density.

Also, we compare in Fig. 3b the numerical results with the approximate growth rates given by Eq. (21). Equation (21) misses the effect of Cowling’s diffusion. Nevertheless, a reasonably good agreement between both results is obtained. The differences get larger when ξ_i becomes small. As happens for Ohm’s diffusion, we find that Cowling’s diffusion has little influence on the thermal mode growth rate. Although Cowling’s diffusion imposes a lower bound for the instability criterion (see Eq. (23)), its influence on the growth rate is almost irrelevant for realistic values of η_C.

5. Discussion

In this paper we have investigated the stability of thermal modes for partially ionized prominence plasmas in the single-fluid approximation. We restricted ourselves to the linear phase when deriving an instability criterion that takes the effects of thermal conduction by electrons and neutrals and Cowling’s diffusion into account. We applied the instability criterion using Hildner’s (Hildner 1974) and Klimchuk-Raymond’s (Klimchuk & Cargill 2001) loss functions, which are frequently used in the literature, along with a new loss function derived from the CHIANTI atomic database (Parenti et al. 2006; Parenti & Vial 2007). Results using both Hildner’s and Klimchuk-Raymond’s loss functions predict that there is thermal instability for temperatures higher than 10⁵ K, approximately. However, the threshold temperature is significantly reduced for about an order of magnitude when the more up-to-date function based on the CHIANTI database is used. In particular, instability can occur at temperatures as low as about 15 000 K. Effects of partial ionization, especially thermal conduction by neutrals, are relevant at these low temperatures.

Focusing on the region of instability at low temperatures obtained with the new CHIANTI-based loss function, we performed a parametric study of the linear growth rate by numerically solving the general dispersion relation. For constant gas pressure, we find that the growth rate decreases as the temperature increases. In addition, the growth rate increases as the amount of neutrals grows. We also find that thermal conduction reduces the growth rate for short wavelengths, with conduction by neutrals being more efficient than conduction by electrons in partially ionized plasmas. However, neither Ohm’s nor Cowling’s diffusion strongly influences the growth rate.

A few remarks should be made about the applicability of our results. Assuming optically thin plasma might not be valid at low prominence temperatures. Owing to finite optical thickness, the actual radiative losses of the plasma would be reduced in comparison to the optically thin case, so that the thermal mode growth rate would decrease consistently (see Carbonell et al. 2006). Some attempts to incorporate the effect of finite optical thickness in Hildner’s parametrization can be found in, e.g., Rosner et al. (1978) and Milne et al. (1979). To our knowledge, the effect of finite optical thickness has not been incorporated in more up-to-date loss functions. We recall that the CHIANTI database is still incomplete at temperatures of T ~ 10⁴ K. Further increase in the loss function at low temperatures is expected for more complete calculations. Thus, these two effects, namely the decrease in the radiative losses due to finite optical thickness and the increase due to the incorporation of additional line emissions, would determine the actual shape of the loss function at low temperatures, hence the actual value of the growth rate of the unstable thermal modes. Also, the growth rates obtained here may be compared to the characteristic time scale for atomic processes. This means that the assumptions of ionization equilibrium of the plasma and electron thermal distribution might not be strictly valid on these short time scales. We plan to investigate this aspect further in the near future.

We conclude that thermal instabilities may take place in prominences at lower temperatures than those predicted with previously existing loss functions. This may be important for the

dynamics and energy balance of the prominence plasma. The obtained growth rates suggest that this low-temperature instability may have an observable effect on prominences. For example, this low-temperature instability may help form density enhancements in regions where the plasma is already cool, such as the prominence threads observed in Hα images. However, nonlinear studies beyond the present normal mode approach are needed to assess the actual impact of the instability on the evolution of the prominence medium. This is relegated to future works.

Acknowledgments

The authors acknowledge discussion within the ISSI Team on Solar Prominence Formation and Equilibrium: New data, new models. R.S. acknowledges support from a Marie Curie Intra-European Fellowship within the European Commission 7th Framework Program (PIEF-GA-2010-274716). R.S. and J.L.B. acknowledge the financial support received from the Spanish MICINN under project AYA2011-22846. R.S. and J.L.B. also acknowledge the financial support from CAIB through the “Grups Competitius” scheme. S.P. acknowledges the support from the Belgian Federal Science Policy Office through the ESA-PRODEX program. CHIANTI is a collaborative project involving George Mason University, the University of Michigan (USA), and the University of Cambridge (UK).

References

All Tables

All Figures

	Fig. 1 Fit of the energy loss rate per unit volume as function of the temperature computed from the CHIANTI atomic database assuming coronal abundances, ionization equilibrium, and a constant pressure of 6.64 mPa (solid line) compared to Hildner’s (dotted line) and Klimchuk-Raymond’s fits (dashed line). The parameters χ∗ and α of the three fits are given in Table 1.
In the text

	Fig. 2 Approximate thermal mode growth rate vs. temperature in the absence of thermal conduction and for Hildner’s (dotted), Klimchuk-Raymond’s (dashed), and CHIANTI-based (solid) loss functions. The shaded area denotes the region of instability at low temperatures obtained with the CHIANTI-based loss function that is not present for the other parametrizations.
In the text

	Fig. 3 a) Numerically computed growth rate (solid) and approximate value obtained from Eq. (21) (dashed) vs. wavelength for the CHIANTI-based loss function and different values of the temperature in the case of a fully ionized plasma, i.e., ξi = 1. b) Same as panel a) but for T = 16   000 K and different values of ξi.
In the text