FORMATION OF SOLAR FILAMENTS BY STEADY AND NONSTEADY CHROMOSPHERIC HEATING

As mentioned above, the plasma beta of the filament environment is believed to be small; therefore, it is generally assumed that the coronal flux tubes, which can support the filaments against gravity, are rigid and the mass flow is channeled along the magnetic field line in the corona.⁴ With such an assumption, the plasma dynamics of the filament threads is simply described by the 1D radiative hydrodynamic equations as follows:

$$\begin{equation} \frac{\partial \rho }{\partial t} + \frac{\partial }{\partial s} (\rho v) = 0, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial }{\partial t} (\rho v) + \frac{\partial }{\partial s} (\rho v^2 + p) = \rho g_\parallel (s), \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial \varepsilon }{\partial t} + \frac{\partial }{\partial s} (\varepsilon v+ p v)=\rho g_\parallel v+H(s)-n_{\rm H} n_{\rm e} \Lambda (T)+\frac{\partial }{\partial s}\left(\kappa \frac{\partial T}{\partial s}\right)\!{,} \end{equation} \tag{ 3 }$$

where ρ is the mass density, T is the temperature, s is the distance along the loop, v is the velocity of plasma, p is the gas pressure, ε = ρv²/2 + p/(γ − 1) is the total energy density, n_H is the number density of hydrogen, n_e is the number density of electrons, and g_∥(s) is the component of gravity at a distance s along the magnetic loop. Furthermore, γ = 5/3 is the ratio of the specific heats, Λ(T) is the radiative loss coefficient for the optically thin emission, H(s) is the volumetric heating rate, and κ = 10⁻⁶ T^5/2 erg cm⁻¹ s⁻¹ K⁻¹ is the Spitzer heat conductivity. As done in the previous works mentioned in Section 1, we assume a fully ionized plasma and adopt the one-fluid model. Considering the helium abundance (n_He/n_H = 0.1), we take ρ = 1.4m_pn_H and p = 2.3n_Hk_BT, where m_p is the proton mass and k_B is the Boltzmann constant. The radiative hydrodynamic Equations (1)–(3) are numerically solved by the Adaptive Mesh Refinement Versatile Advection Code (AMRVAC; Keppens et al. 2003, 2011), where the heat conduction term is solved with an implicit scheme separately from the other terms. To calculate the radiative energy loss, we use second-order polynomial interpolation to compile a high-resolution table based on the radiative loss calculations recently done by Colgan et al. (2008). They calculated the radiative losses for the solar coronal plasma using a recommended set of quiet-region element abundances. In their calculations they used a complete and self-consistent atomic data set and an accurate atomic collisional rate over a wide temperature range. As shown in Figure 1, Λ(T) in our cooling table (solid line), interpolated from Colgan et al. (2008) (square), is generally ∼2 times larger than the Klimchuk–Raymond radiative loss function (dashed line) used in previous works (Karpen et al. 2005, 2006; Karpen & Antiochos 2008; Klimchuk et al. 2010). The figure also demonstrates that our cooling curve better represents the detailed temperature dependence of the radiative loss.

Zoom In Zoom Out Reset image size

Using our cooling table, we then exploit the exact integration scheme (Townsend 2009), rather than traditional implicit or explicit time stepping methods. This method is much faster than an explicit scheme, as it can avoid the numerical limit of the radiative timescale on the simulation timestep. It is also more stable than the implicit schemes based on Newton–Raphson iteration. Below 20,000 K, we set Λ(T) to vanish since the plasma then becomes optically thick and is no longer fully ionized. The use of explicit, (semi-)implicit, and exact integration methods in grid-adaptive simulations has been analyzed recently by van Marle & Keppens (2011).

When using the AMRVAC code, the total Variation Diminishing Lax–Friedrichs scheme using linear reconstruction employing a monotonized central limiter (Tóth & Odstrčil 1996) is chosen for the spatial differentiation, combined with a predictor-corrector two-step explicit scheme for the time progressing. Six levels of adaptive mesh refinement (AMR) in a block-based AMR approach are applied, which leads to a minimum grid spacing of 6.77 km, comparable to the 5–6 km seen in previous works such as Klimchuk et al. (2010). The refinement/coarsening criteria are based on numerical errors estimated using density and its gradient following Löhner's prescription (Löhner 1987). If any error exceeds 0.1, the block is refined. If all errors in the block are less than 0.0125, the block is then coarsened. To include the heat conduction source in the energy equation, we separately solve the heat conduction term in each AMR grid block using the implicit scheme where the central difference is taken for the space derivative of the temperature. This leads to a local tri-diagonal linear system per grid block, where the temperature on the block boundaries is taken from neighboring blocks at the previous timestep. To simulate 10 hr physical time, our implementation uses ∼1.5 hr on four processors.

We adopt a loop geometry with a magnetic dip, which is symmetric about the midpoint. On each side, the loop has a vertical leg 5 Mm in length above the footpoint and a quarter-circular arc 15.7 Mm in length connecting the vertical leg and the dip, which is 218.6 Mm in length, as shown in Figure 2. Note that the geometry of the loop determines the distribution of g_∥(s), which is symmetric about the midpoint and whose left half is described as follows:

$$\begin{equation} g_\parallel (s)=\!\left\{\begin{array}{@{}ll} \dsty {-}g_\odot, & s \leqslant s_1; \cr \dsty {-}g_\odot \cos \left(\frac{\displaystyle \pi }{\displaystyle 2}\frac{\displaystyle s-s_1}{\displaystyle s_2-s_1}\right), & s_1 < s \leqslant s_2; \cr \dsty g_\odot \frac{\displaystyle \pi D}{\displaystyle 2(L/2-s_2)} \sin \left(\pi \frac{\displaystyle s-s_2}{\displaystyle L/2-s_2}\right), & s_2 < s \leqslant L/2, \end{array}\right. \end{equation} \tag{ 4 }$$

where g_☉ = 2.7 × 10⁴ cm s⁻² is the solar gravity, s₁ = 5 Mm, s₂ = s₁ + 15.7 Mm, L = 260 Mm is the total loop length, and D = 0.5 Mm is the dip depth. The value of the total loop length L is suggested by the observations of Okamoto et al. (2007). The dip is very shallow, so the coronal part of the loop is nearly flat. The midpoint of the loop, which is the center of the dip, has a height of 14.5 Mm above the bottom boundary.

Zoom In Zoom Out Reset image size

The initial equilibrium state is obtained by numerically solving Equations (1)–(3). We start with a temperature (T) versus height (h) distribution of T = tanh (h − h₀) with T = 10⁶ K in the corona and 6000 K in the photosphere, which is close to the quiet Sun atmospheric model (Vernazza et al. 1981). The density is determined by balancing the pressure gradient with the gravity, with n_H = ρ/(1.4m_p) = 10⁹ cm⁻³ at the loop center. At this stage, only the background heating H₀(s) is included in the energy equation, namely, H(s) = H₀(s) in Equation (3). H₀(s) is a steady term in order to maintain the hot corona whose physics is still under debate. Considering that the photospheric motions are the source of the energy that is transported upward to heat the chromosphere and the corona, somehow, it is generally conjectured that the heating rate decays with height (Serio et al. 1981; Mok et al. 1990; Aschwanden & Schrijver 2002). Similar to previous works, we assume that H₀(s) decreases exponentially with the distance away from the nearest footpoint along the loop and remains constant in time,

$$\begin{equation} H_0(s)= \left\{ \begin{array}{@{}ll} E_0 \exp (-s/H_m), & s < L/2; \\ \ms E_0 \exp [-(L-s)/H_m], & L/2 \leqslant s < L,\end{array}\right. \end{equation} \tag{ 5 }$$

where the amplitude E₀ = 3 × 10⁻⁴ erg cm⁻³ s⁻¹ (Withbroe & Noyes 1977) and the scale length H_m = L/2 (Withbroe 1988). The prescribed distributions are in force equilibrium, but not in thermal equilibrium, and will evolve to reach a new hydrostatic state. Such a state, whose density and temperature distributions are displayed in Figure 3, serves as the initial conditions for our further simulations. In the initial state, the temperature is highest at the midpoint of the loop, with T = 2.6 × 10⁶ K and n_H = 3.2 × 10⁸ cm⁻³. The thin transition layer between the lower atmosphere and the corona is roughly at a height of s_tr = 6 Mm. In the lower atmosphere, T ranges from 13,000 K to 18,000 K, which is closer to the quiet Sun atmospheric model (Vernazza et al. 1981) than previous works (Karpen et al. 2001, 2006). The number density at the endpoints is about 2 × 10¹⁴ cm⁻³. The simulated chromosphere and photosphere are about twice as thick as those of the real Sun and serve as a mass reservoir for the chromospheric evaporation.

Zoom In Zoom Out Reset image size

For the boundary conditions, we fix the density, velocity, and temperature at the two endpoints of the loop. Because the density in the photosphere is more than four orders of magnitude higher than that in the filaments and the corona, the coronal dynamics has little effect on the photosphere, justifying these fixed boundary conditions.

Similar to Antiochos et al. (1999), in order to simulate the chromospheric evaporation, an extra localized heating H_l(s), which might be due to chromospheric reconnection, is added to the energy equation in addition to the background heating H₀(s), namely, H(s) = H₀(s) + H_l(s) in Equation (3). As described as follows, H_l(s) is uniform in the photosphere and chromosphere and decays exponentially with the distance away from the nearest chromosphere along the loop with a scale length λ:

$$\begin{equation} H_l(s)= \!\left\{\begin{array}{@{}ll} \dsty E_1, & s \leqslant s_{tr}; \\ \ms \dsty E_1 \exp [-(s-s_{tr})/\lambda ], & s_{tr} < s \leqslant L/2; \\ \ms \dsty f E_1 \exp [-(L-s_{tr}-s)/\lambda ], & L/2 < s \leqslant L-s_{tr}; \\ \ms \dsty f E_1, & s > L-s_{tr}, \end{array}\right. \end{equation} \tag{ 6 }$$

where the amplitude E₁ = 10⁻² erg cm⁻³ s⁻¹ (cf. Withbroe & Noyes 1977; Aschwanden 2001), s_tr = 6 Mm is the height of the transition region, and the factor f is the ratio of the localized heating rate near the right footpoint to that near the left. The localized heating H_l(s) is ramped up linearly over 1000 s and maintained thereafter. In this paper, we numerically investigate two situations, with symmetric and asymmetric heating, respectively. The parameters in several typical cases are listed in Table 1. In the symmetric case, a parameter survey is performed, including the effects of λ and E₁.

To better identify in which way our simulations augment the knowledge of prominence formation gained over the last decade, we list the most important parameters in similar works on radiative condensation due to localized heating in Tables 2 and 3. These tables show that our work differs in a variety of aspects connected to the overall loop geometry, to the spatio-temporal prescription of the heating applied, and also notably in the cooling table used to quantify radiative losses. Motivated by observations of active region prominences by Okamoto et al. (2007), our model loop shape represents a low-lying, shallowly dipped loop with a more realistic scale for its vertical legs and chromospheric height region. In this shallow dip configuration, our parametric survey explores a wide range of heating parameters.

As the symmetric localized heating with λ = 10 Mm and f = 1 is introduced in case S1, the chromospheric plasma is heated and evaporated into the corona. As illustrated by Figure 4, both the density ρ and temperature T in the coronal portion increase accordingly. Near the midpoint of the loop, T reaches the maximum value, 3.49 × 10⁶ K, at t = 2664 s, then starts to decline slowly, whereas ρ keeps increasing slowly from the beginning. At t = 10013 s, the temperature and pressure near the midpoint begin to collapse simultaneously, drastically decreasing by nearly one and a half orders of magnitude within 1 minute, creating a low-pressure cold region, which expands to a maximum length of 28.4 Mm. At this stage, the density is still low, increasing gently as seen in the left and middle columns of Figure 5. Since only the density was chosen to automatically refine the mesh and at this stage the density is still smooth, the grid resolution is 108 km per cell. Still, this cold region contains 263 grid cells, which is sufficient to resolve the region. Due to the large pressure gradient that forms at the edge of this cold region, the coronal plasma outside the cold region is driven to move rapidly toward this central cold region. After time t = 10082 s, the density inside this cold well begins to increase rapidly as the inflows from both sides converge toward the midpoint and compress the cold region. The converging velocity reaches 185 km s⁻¹. These inflows are supersonic, with a local Mach number up to 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a result, the inflows collide at the midpoint of the loop where a high-pressure peak appears, exciting two rebound shock waves launched from the midpoint toward the two sides. A small cold condensation region (∼1 Mm in length) is left behind the shocks near the midpoint at t ∼ 10382 s (see the right column of Figure 5). To see the contributions of the various terms in the energy equation, in Figure 6 we plot the absolute value distribution of the energy source terms, including radiative cooling, heat conduction, heating, and gravitational potential, across the magnetic dip at t = 10382 s, when the condensation happens. It is found that the radiative cooling dominates in the coronal parts and the boundaries of the condensation segment, but nearly vanishes inside the condensation. The heat conduction is less important than the heating in most regions except at the boundaries of the condensation. The gravitational potential is always negligible. The pressure in the cold region recovers due to the compression of the inflows from outside. Swept by the outward-propagating rebound shock waves, the depressed pressure outside the cold region also recovers, as illustrated by Figure 7, which shows similar quantities at times later than those of Figure 5. The shock waves are bounced back and forth ∼3 times between the loop footpoint and the loop center, as revealed by the sinusoidal pattern in the right panel of Figure 4 between t = 3 hr and t = 4 hr. During their passage, they dissipate their energy to compress and heat the local plasma. The damping rate is enhanced by thermal conduction and radiation. The plasma condensation remains near the midpoint, with a temperature of 1.8 × 10⁴ K and a density of 1.2 × 10¹¹ cm⁻³. In contrast, the corresponding values in the neighboring corona are 2 × 10⁶ K and 1.03 × 10⁹ cm⁻³, respectively. Note that the condensation temperature is just below 20,000 K, where the radiative loss is set to vanish smoothly. Further tests indicate that if the radiative losses vanished below a lower temperature, the condensation would be cooler accordingly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 depicts the growth of the condensation segment (or the filament thread). It is seen that the onset time of the condensation is at t = 2.8 hr after the localized heating is introduced in case S1. The growing process of the condensation can be described as follows. Its length increases rapidly for ∼20 minutes as the onset of condensation drives fast evaporation flows from the chromosphere, followed by a slight shrinkage of the condensation for ∼10 minutes. As the evaporated plasma flow becomes steady, the condensation length increases linearly with time, with a growth rate of 1511 km hr⁻¹. With such a speed, it would take ∼6.6 hr to form a filament thread with a typical length of 10 Mm. Observations show that active region filaments form within a day (Wang & Muglach 2007). For comparison, in Antiochos et al. (1999), it takes 8.3 hr for a condensation to grow to 10 Mm long. One reason is that they used a deeply dipped magnetic loop, where the gravity scale height was shorter, so that the condensation was strongly squeezed.

Zoom In Zoom Out Reset image size

To investigate the effect of the heating scale length λ, we change its value and perform a series of simulations in 17 runs, with other parameters being the same as in case S1. As seen in Figure 9, the onset time of the condensation roughly increases with λ, with a minimum of 2 hr. However, the growth rate decreases with increasing λ, except for a drop down near λ = 4 Mm. It is noted that if λ is larger than 9 Mm (i.e., 1/28 of the total loop length L), the evolution is similar to case S1 (where λ = 10 Mm) as described above, where only a single condensation forms near the midpoint of the loop. When 3 Mm <λ < 9 Mm (i.e., 1/86 < λ/L < 1/28), two condensation segments first form on the two shoulders of the magnetic dip symmetrically about the midpoint. The two segments move convergently toward the midpoint, during which movement both T and p in the region between the two segments drop down. Under this pressure gradient, the two condensation segments are accelerated from ∼12 km s⁻¹ to 75 km s⁻¹, to finally coalesce near the midpoint. Similar high-speed motion is discussed by Karpen et al. (2006). As λ decreases, the two segments form further away from each other and from the midpoint of the loop. When 2.5 Mm <λ < 3 Mm, the two condensation segments form in the loop legs and then drain down rapidly to the nearby footpoints. When λ < 2.5 Mm (i.e., λ/L < 1/100), no condensation forms, and the loop relaxes to a hydrostatic state in the end. Similar situations happen when λ > 25 Mm (i.e., λ/L > 1/10), which is the same result as mentioned by Klimchuk et al. (2010).

Zoom In Zoom Out Reset image size

Müller et al. (2004) simulated the formation of condensations in a semicircular coronal loop with a length of 100 Mm and a loop-top temperature of 6.8 × 10⁵ K. Their loops are shorter and cooler than our dipped loops. They found that, when λ/L = 1/20, only one condensation forms at the midpoint of the loop and that two condensation segments form at the shoulders of the loop when λ/L = 1/33 or λ/L = 1/50. The transition between one condensation and two condensation segments is somewhere between λ/L = 1/20 and 1/33, which is consistent with our result, that is, 1/28. The transition can be understood as follows. At the beginning, the plasma at the two shoulders of the loop is cooler and denser than that at the midpoint, which means the radiative loss is stronger at the shoulders. Meanwhile, if λ is long enough, the heating at the shoulders is strong enough to slow down the cooling, making the midpoint the fastest cooling place, and only one condensation forms. When λ decreases, the heating at the shoulders is reduced and becomes too weak to obstruct the fast cooling there. Hence, two cold segments are formed at the two shoulders before merging near the midpoint of the loop.

Keeping λ = 10 Mm, we perform another series of simulations with different amplitudes of the localized heating, E₁, namely, 0.005–0.2 erg cm⁻³ s⁻¹. The localized heating still dominates compared to the weak background heating. The evolution is similar to case S1, in that only one condensation forms near the midpoint. As indicated in Figure 10, the onset time of the condensation decreases as E₁ increases, whereas the growth rate of the condensation is maximal at E₁ ∼ 0.01 erg cm⁻³ s⁻¹. It is easy to understand that the condensation forms earlier with a larger E₁ since stronger chromospheric heating leads to stronger evaporation. The growth rate might be determined by the compromise between the evaporation rate and the deposited energy. Stronger heating drives stronger chromospheric evaporation, on the one hand, and restrains the evaporated plasma from cooling, on the other. Another important factor is that the plasma density in the condensation segment increases due to higher compression as E₁ increases. It becomes more difficult for a denser condensation to grow quickly. As a combined result, the growth rate of the condensation peaks at E₁ ∼ 0.01 erg cm⁻³ s⁻¹ and decreases at larger E₁.

Zoom In Zoom Out Reset image size

Comparing the onset time and growth rates deduced from observations may help to further pin down the properties of the localized heating employed. However, to determine the onset time of the thread formation requires combined spectral and imaging observations with high resolution, which are not available yet. The growth speed of the filament thread can be compared with future observations.

A more general case is that in which the localized heating in the chromosphere is not symmetric between the two footpoints of a magnetic loop. In this subsection, we describe three simulations with different f, the ratio of the heating rate at the right footpoint to that of the left.

In case A1, we set f = 0.75, and other parameters are the same as in case S1. The details of the formation process are similar to case S1, as illustrated by Figure 11, which depicts the time evolution of the density (left panel) and temperature (right panel) distributions. At t = 10690 s (i.e., 2.97 hr), a condensation with low temperature and high density is formed in the right part (i.e., the less heated part), of the magnetic dip, 22.5 Mm away from the midpoint. The distributions of various quantities, namely, the temperature (T), the density (n), the pressure (p), and the in situ Mach number (M), across the condensation segment at three times are shown in Figure 12. It is seen that, as the convergent inflows coalesce in the condensation segment, two rebound shock waves are launched, propagating toward the left and right footpoints, respectively. The shock waves are reflected between each footpoint and the condensation several times before fading away, as indicated by the sinusoidal pattern in the right panel of Figure 11 near t = 3.5 hr. A significant difference from case S1, however, is that upon formation, the condensation has a velocity of ∼5 km s⁻¹, moving to the right, or the less heated side, as illustrated by Figure 11. The rightward-moving condensation is accelerated to 15 km s⁻¹ due to the pressure gradient on its two sides, but it is soon made to decelerate by an inversion of the pressure gradient and even falls back by 2 Mm. The pressure gradient inversion is caused by a pressure increase on the right side of the condensation due to interaction with a shock wave that is reflected at the right footpoint, when the left-side shock has not yet reached the left footpoint. After the left-side shock wave is reflected from the left footpoint and catches up with the condensation, the pressure on the left side of the condensation increases and exceeds the pressure on the right side. The condensation is then pushed again by the pressure gradient to move to the right with a velocity of ∼8 km s⁻¹ until it drains down to the right footpoint of the loop. During the travel, the length of the condensation grows from 0.9 Mm to 7.9 Mm. When the condensation impacts the chromosphere, the collision generates a rebound shock wave, which is bounced back and forth between the two footpoints of the loop, as indicated by the sinusoidal pattern in the right panel of Figure 11 near t = 8 hr. The total lifetime of the condensation is ∼4.3 hr. As the simulation goes on, the formation and the drainage of condensation repeat, with a period of 5.5 hr, which is significantly shorter than the 22.8 hr simulated by Karpen et al. (2006). Since the radiative loss coefficient we used is generally ∼2 times larger than theirs, the cooling is stronger and the condensations form more rapidly, which shortens the period of the formation-drainage cycle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned in Section 3.1, two condensations can be formed when the heating scale length λ is small. Therefore, in case A2, we take f = 0.4 and λ = 5 Mm. As illustrated by Figure 13, a condensation forms at t = 2.2 hr near the left shoulder of the magnetic dip (s = 42.1 Mm) with an initial speed of 10 km s⁻¹. It is accelerated during its travel toward the right part of the loop, with its length growing to 5.8 Mm. At t = 2.84 hr (38 minutes later), a second condensation, which is smaller than the first one, is formed near the right shoulder of the magnetic dip (s = 205.5 Mm) with an initial speed of 24 km s⁻¹ toward the left. The left and right condensation segments are accelerated for 24 minutes to 60 km s⁻¹ and 50 km s⁻¹, respectively, and then collide at s = 158 Mm (in the right part of the loop). Two shock waves are generated by the collision, which are then reflected back and forth between the condensation and each footpoint, as revealed by the sinusoidal pattern in the right panel of Figure 13 near t = 3.5 hr. After the collision, the two condensation segments merge into one, with a length of 2.3 Mm. The coalesced condensation moves to the right with an initial velocity of 32 km s⁻¹. It is decelerated to 16 km s⁻¹ at s = 216 Mm and then accelerated to 24 km s⁻¹ with a length of 6 Mm before it drains down to the right footpoint of the loop. The deceleration and acceleration of the condensation are mainly due to the inversion of the pressure gradient caused by the same effect as discussed in case A1, while the gravity effect is small in this shallow dip configuration. The falling down of the condensation excites a shock wave, which is trapped to propagate back and forth in the whole loop as indicated by the right panel of Figure 13 near t = 5 hr. As time goes on, such formation, coalescence, and drainage of condensations repeat with a period of 3.6 hr, which is shorter than that in case A1.

Zoom In Zoom Out Reset image size

As an extreme case, we perform a simulation with f = 0, that is, the localized heating is introduced at the left footpoint only. It is found that no condensation forms in the loop. Instead, we get steady flows along the coronal loop, consistent with previous works (e.g., Patsourakos et al. 2004).

Parker (1953) proposed that some solar features, such as filaments, can be formed by thermal instability. He derived a criterion for thermal instability on the basis of an analysis of the energy equation alone. Field (1965) made a detailed study of the thermal instability for an infinite, uniform, static plasma in initial thermal equilibrium. He pointed out that the criterion given by Parker is based on the isochoric assumption, that is, the density is constant in the whole region, which is not compatible with the force equation since the cooling would lead to a pressure deficit, which would destroy the initial force balance. He derived an isobaric criterion for thermal instability, which is consistent with the force equation. Thermal instability was further studied by many other colleagues (e.g., van der Linden & Goossens 1991; Meerson 1996 and references therein). It was pointed out that these modes are the marginal entropy modes that are advected with the local flow velocity (Goedbloed et al. 2010), driven unstable by non-adiabatic processes. The different criteria may be applicable for different astrophysical environments.

According to our simulations, as mentioned in Section 3.1, during the catastrophic cooling stage, the temperature and the pressure drop rapidly, while the density increases only a little. Significant density enhancement occurs ∼3 minutes after the catastrophic cooling. Therefore, in the simulated coronal loop, the thermal catastrophe is more isochoric than isobaric. So, we use the isochoric thermal instability criterion derived by Parker (1953) as follows:

$$\begin{equation} C\equiv k^2- \frac{1}{\kappa } \left(\frac{\partial H(s)}{\partial T} -\frac{\partial R}{\partial T}\right) < 0, \end{equation} \tag{ 7 }$$

where k is the wave number of the perturbations, κ is the heat conduction coefficient, and R = n_Hn_eΛ(T) is the radiative loss. The heat conduction introduces a stabilizing effect. We numerically calculate ∂R/∂T, using the central difference scheme. ∂H(s)/∂T is zero since the heating depends only on distance in our simulations. Perturbations with any resolvable wavelength exist in the simulations. According to Figure 5, the cool region has a width of ∼10 Mm, therefore, we take the wavelength of the temperature perturbation as 20 Mm in order to quantify k. For small k, the occurrence of the thermal instability is mainly determined by the sign of ∂R/∂T.

Taking case S1 as an example, we plot the temporal evolution of the temperature, the density, the pressure, and the isochoric criterion C at the loop midpoint in Figure 14. Since the initial temperature is 2.63 MK, which corresponds to a negative ∂R/∂T, C = −1.4 × 10⁻¹⁵ cm⁻² is negative, but very close to 0. Besides, the cooling timescale at this stage is ∼10⁴ s. Therefore, the early evolution is dominated by the localized heating and chromospheric evaporation. As more mass fills the corona, radiation is enhanced gradually, which overwhelms the heating after t = 2664 s. The temperature then keeps decreasing slowly. From t = 9850 s to t = 9994 s, C becomes positive for a short interval since T falls in the range where ∂R/∂T is positive. After t = 10013 s, C drops down drastically to −1.2 × 10⁻⁹ cm⁻². Simultaneously, the temperature T, along with the gas pressure, begins to decrease catastrophically, as indicated by the time derivative of T (the dashed line in the top panel of Figure 14). That is to say, thermal instability occurs. The temperature drops from 3.4 × 10⁵ K to 20,000 K in 60 s. However, the density increases by only 20% during this time. Note that C becomes positive out of the plotting range after the catastrophic cooling, corresponding to a thermally stable state. Three minutes later (i.e., at t = 10282 s), the density increases drastically, and a condensation is then formed. We conclude that the isochoric thermal instability may explain the catastrophic cooling. Such a delay is probably due to the difference between the kinematic timescale and the radiative timescale. It takes an extra 3 minutes for the plasma, driven by the pressure gradient, to accumulate in the cooling region.

Zoom In Zoom Out Reset image size

It might be interesting to check the criterion of isobaric thermal instability. According to Equation (25) of Field (1965, see also van der Linden & Goossens 1991), the criterion for thermal instability in the isobaric case is expressed as

$$\begin{equation} C_{\rm isobaric}\equiv \rho \left(\frac{\partial \mathcal {L}}{\partial T}\right)_{\rho } - \frac{\rho ^2}{T}\left(\frac{\partial \mathcal {L}}{\partial \rho }\right)_T+k^2 \kappa < 0, \end{equation} \tag{ 8 }$$

where $\mathcal {L}=(n_{\rm H} n_{\rm e} \Lambda (T)-H(s))/\rho$ is the generalized heat-loss function. For perturbations with a wavelength of 20 Mm, we calculate C_isobaric at the midpoint of the loop and plot its temporal evolution in the bottom panel of Figure 14. The time at which it turns from positive to negative is well before the onset time of the catastrophic cooling. We tried many other perturbation wavelengths and found that the isobaric criterion is crucially dependent on the perturbation wavelength while the isochoric criterion is not. Therefore, we conclude that isobaric thermal instability is not appropriate to explain the catastrophic cooling during condensation formation in the solar corona.

As mentioned in Section 1, it has been demonstrated that the extra heating localized in the lower atmosphere would drive chromospheric evaporation flows, leading to the plasma condensation in the corona due to thermal instability or loss of thermal equilibrium. In previous studies, the localized heating is either continuous (Antiochos et al. 1999) or intermittent (Karpen & Antiochos 2008). In the steady heating case, the condensation can form and grow rapidly, as also demonstrated in this paper. In the successive impulsive heating case, it was found that a condensation can also form steadily when the average interval between heating pulses is less than the coronal radiative cooling time (∼2000 s; Karpen & Antiochos 2008). From the theoretical point of view, the localized strong heating may be due to lower atmospheric activities, such as chromospheric reconnection. It is quite possible that such a heating event has a finite lifetime and not show up again at the footpoint of a given flux tube. Therefore, it is interesting to see the response of the coronal loop to a single heating event. To do that, we make a numerical experiment and stop the localized heating 8 minutes after the condensation is formed at t = 2.87 hr in the symmetric case S1, which means there is only background heating H₀(s) in the heating term H(s) in Equation 3, which changes slightly with the distance along the loop. Note that in order to make the evolution smoother, the heating is turned off linearly over 1000 s. We find that, after the heating ceases, there is still mass upflow from the footpoint to the coronal portion of the loop. For comparison, Figure 15 plots the evolution of the mass flux from the two shoulders of the magnetic dip to the condensation segment in the steady case (top panel) and the finite-time heating case (bottom panel). It is seen that after the condensation is formed at t = 2.87 hr, the mass flux oscillates heavily and then maintains a level of 1.5 × 10⁻⁸ g cm⁻² s⁻¹ in the steady heating case. However, in the finite-time heating case, the mass flux drops abruptly but then still maintains a level of 5.6 × 10⁻¹⁰ g cm⁻² s⁻¹, even after the localized heating is permanently removed. That is to say, the mass flux decreases by 27 times, but does not vanish. Such a mass flux corresponds to a growth rate of the condensation length of 230 km hr⁻¹, a decrease of 6.6 times compared to the steady heating case. The reason why the growth rate does not decrease proportionally with the mass flux is that the plasma density of the condensation is reduced after the heating is switched off.

Zoom In Zoom Out Reset image size

Therefore, it seems that evaporation-driven condensation can serve as the trigger of the formation of filament threads, and there exists a condensation instability. Once the condensation is formed in a small segment near the dip of a coronal loop by the evaporation flow, the condensation will grow, and continual mass supply is siphoned from the chromosphere, although the growth rate is ∼6.6 times smaller than in the steady heating case.

In order to understand the mechanism of the spontaneous siphon flow after the localized heating is halted, we plot in Figure 16 the time evolution of the gas pressure at the midpoint of the loop, where the condensation is located. It is revealed that, after the localized heating is turned off gradually, the pressure at the loop midpoint drops down after t = 3.03 hr. After several hours of small-amplitude oscillations, the pressure remains at 0.2 dyn cm⁻², which is about half of its initial value 0.37 dyn cm⁻² at t = 0. Therefore, our simulation result indicates that, after thermal instability and plasma condensation, the gas pressure of the cold plasma is reduced, compared to that of the hot plasma at the same site in the initial hydrostatic state, which leads to a pressure gradient along the loop. It is such a pressure gradient that drives the spontaneous siphon flow, which makes the condensation continue to grow even after the localized heating is switched off.

Zoom In Zoom Out Reset image size

While the previous section demonstrated that prominence growth will continue even after finite-time localized heating, another aspect worth studying is the fate of prominence condensations subjected to wave buffeting. Since 5 minute solar p-mode oscillations are ubiquitous in the photosphere, it is relevant to investigate whether prominences will be influenced by wave driving, and how they channel (linear) wave modes. Since there is great interest in prominence seismology (see, e.g., the review by Mackay et al. 2010), we now investigate how filaments, once formed, behave under p-mode wave driving, which is denoted as case D1. Case D1 is based on the simulation results of the symmetric case S1, with an initial state taken from 7.1 hr in the evolution shown in Figure 4. We then introduce a sinusoidal velocity perturbation with an amplitude of 1 km s⁻¹ and a period of 5 minutes. We add this perturbation only at the left footpoint of the loop. The p-mode waves propagate upward through the transition region into the corona and steepen into shocks. A snapshot of the velocity distribution after 509 minutes is shown in Figure 17, where two shock fronts can be identified to the left of the central condensation. The shocks damp in the corona while propagating with a speed of 213 km s⁻¹, which is close to the local sound wave speed. These shocks hit the condensation and penetrate into it with little reflection, becoming ordinary linear sound waves while propagating through the filament. When these sound waves hit the other boundary of the condensation, they are mainly reflected, with a weak leakage out to the right part of the loop. This is best visualized using a Schlieren plot of the pressure, as shown in Figure 18. This Schlieren plot of the pressure zooms in around the condensation, and actually quantifies the local value of exp (− 0.01[|∇p| − 500]). The bottom part of Figure 18 shows the first series of shocks hitting the condensation. One notes that the overall thermodynamical changes induce a leftward drift of the central prominence. The symmetry is thus broken due to our asymmetric driving. During the leftward drifting, the filament thread becomes longer due to the chromospheric evaporation as in case S1. The top panel of Figure 18 shows that ultimately the filament settles down at a quasi-permanent location determined by the overall pressure balance, as altered by the periodic driving. The wave mode reflections and transmissions at the left and right edges of the filament thread can be clearly detected. The different slopes of the wave fronts within the filament is due to the lower sound speed there. All these waves, whether inside or outside the filament thread, have a 5 minute period. The impinging waves are clearly nonlinear, while the internal wave modes are primarily linear, and the transmitted waves are further attenuated. During the entire period simulated, the condensation remains thermally stable under the perturbation and energy damping from p-mode waves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our model is appropriate for low-lying, shallowly dipped loops, which have not been studied before. Such shallow dip configurations facilitate higher speeds for displacing the condensations. As advocated by studies of hydrostatic coronal loops by Aschwanden & Schrijver (2002), our background heating uses an exponential height dependence, differing from the uniform prescription used in previous works. In case S1, the fine details of the condensation process are shown for the first time. Our parameter survey finds new complicated cycles of paired filament formation, with high-speed converging condensations in case A2, with strongly asymmetrie, short-scale localized heating. With the improved cooling table, we allow for stronger radiative cooling and therefore find shorter formation timescales compared to those of previous works. Besides continuous localized heating, the case with finite duration, which is more realistic in the solar atmosphere, is investigated for the first time, indicating that additional strong heating is not needed to maintain the growth of condensations after their onset. We studied the dependence of the formation process for a large parameter range of the heating scale length λ and the heating amplitude E₁. It is found that shorter λ or stronger E₁ can make condensations form earlier, say, ∼2 hr after the introduction of the localized heating. Shorter λ also leads to faster growth of the condensation. Moreover, the effect of p-mode waves is studied for the first time in this context.