CORONAL RAIN IN MAGNETIC ARCADES: REBOUND SHOCKS, LIMIT CYCLES, AND SHEAR FLOWS

Our numerical setup follows our previous 2.5D thermodynamic MHD simulation from Fang et al. (2013), which includes gravity, field-aligned heat conduction, and radiative cooling and parametrized heating terms, on a rectangular plane with horizontal extension −40 Mm $\leqslant \;x\leqslant 40$ Mm and vertical extension $0\leqslant y\leqslant 50$ Mm. The governing equations are as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{vv}}+{p}_{\mathrm{tot}}{\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{BB}}}{{\mu }_{0}}\right)=\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot \left(E{\boldsymbol{v}}+{p}_{\mathrm{tot}}{\boldsymbol{v}}-\displaystyle \frac{{\boldsymbol{v}}\cdot {\boldsymbol{B}}}{{\mu }_{0}}{\boldsymbol{B}}\right)\\ &&\quad =\rho {\boldsymbol{g}}\cdot {\boldsymbol{v}}+{\rm{\nabla }}\cdot ({\boldsymbol{\kappa }}\cdot {\rm{\nabla }}T)-Q+H,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot ({\boldsymbol{vB}}-{\boldsymbol{Bv}})=0,\end{eqnarray} \tag{ 4 }$$

where T, $\rho ,{\boldsymbol{B}},{\boldsymbol{v}}$, and ${\boldsymbol{I}}$ are respectively temperature, density, magnetic field, velocity, and unit tensor, with the ratio of specific heats γ = 5/3, and a total energy density of $E=p/(\gamma -1)+\rho {v}^{2}/2+{B}^{2}/2{\mu }_{0}$; ${p}_{\mathrm{tot}}\equiv p+{B}^{2}/2{\mu }_{0}$ is the total pressure, consisting of magnetic pressure and thermal pressure p; g $=\;{g}_{0}{R}_{\odot }^{2}/{\left({R}_{\odot }+y\right)}^{2}\hat{y}$ is the gravitational acceleration with the solar surface gravitational acceleration ${g}_{0}=-274\;{\rm{m}}\;{{\rm{s}}}^{-2}$ and the solar radius ${R}_{\odot }$; H and Q are respectively the heating and radiative loss rates; and ${\boldsymbol{\kappa }}$ is the thermal conductivity tensor. Assuming a 10:1 abundance of hydrogen and helium of completely ionized plasma, we obtain $\rho =1.4\;{m}_{p}{n}_{{\rm{H}}}$, where ${m}_{p}$ is the proton mass and ${n}_{{\rm{H}}}$ is the number density of hydrogen. We use the ideal gas law $p=2.3{n}_{{\rm{H}}}{k}_{{\rm{B}}}T$, where ${k}_{{\rm{B}}}$ is the Boltzmann constant. We also adopt $Q=1.2{n}_{{\rm{H}}}^{2}{\rm{\Lambda }}(T)$ as the radiative cooling term, where ${\rm{\Lambda }}(T)$ is the radiative loss function for optically thin emission, quantified by Colgan et al. (2008) using a recommended set of quiet-region element abundances, as used in our previous work (Xia et al. 2011, 2012; Fang et al. 2013; Keppens & Xia 2014). In the calculations, Colgan et al. (2008) used a complete and self-consistent atomic data set and an accurate atomic collisional rate over a wide temperature range. Below 10,000 K, we set ${\rm{\Lambda }}(T)$ to vanish because the plasma there is optically thick and no longer fully ionized. We use the exact integration method introduced by Townsend (2009) to evaluate the radiative loss term. The use of explicit, (semi-)implicit, and exact integration methods in grid-adaptive simulations has been compared in van Marle et al. (2011). The term containing ${\boldsymbol{\kappa }}={\kappa }_{| | }\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}$ quantifies the anisotropic thermal conduction along the magnetic field lines, composed of the unit vector $\hat{{\boldsymbol{b}}}$ along the magnetic field and the Spitzer conductivity ${\kappa }_{| | }$ as ${10}^{-6}\;{T}^{5/2}$ erg cm⁻¹ s⁻¹ K${}^{-3.5}$.

We employ a linear force-free magnetic field for the initial magnetic configuration, which is characterized by a constant angle ${\theta }_{0}$ as follows:

$$\begin{eqnarray}{B}_{x} & = & -{B}_{0}\mathrm{cos}\left(\displaystyle \frac{\pi x}{{L}_{0}}\right)\mathrm{sin}{\theta }_{0}\mathrm{exp}\left(-\displaystyle \frac{\pi y\mathrm{sin}{\theta }_{0}}{{L}_{0}}\right),\\ {B}_{y} & = & {B}_{0}\mathrm{sin}\left(\displaystyle \frac{\pi x}{{L}_{0}}\right)\mathrm{exp}\left(-\displaystyle \frac{\pi y\mathrm{sin}{\theta }_{0}}{{L}_{0}}\right),\\ {B}_{z} & = & -{B}_{0}\mathrm{cos}\left(\displaystyle \frac{\pi x}{{L}_{0}}\right)\mathrm{cos}{\theta }_{0}\mathrm{exp}\left(-\displaystyle \frac{\pi y\mathrm{sin}{\theta }_{0}}{{L}_{0}}\right),\end{eqnarray} \tag{ 5 }$$

with ${\theta }_{0}=30^\circ$; the arcade makes a 30° angle with the neutral line ($x=0,y=0$). ${L}_{0}=80$ Mm is the horizontal size of our domain from −40 to 40 Mm, and when adopting ${B}_{0}=12$ G, our magnetic arcade has a total box averaged field strength of 2.9 G.

For the initial thermal structure, we set a uniform temperature of 10,000 K below a height of 2.7 Mm and choose a temperature profile with height ensuring a constant vertical thermal conduction flux (i.e., $\kappa \partial T/\partial y$ = 2 × 10⁵ erg cm⁻² s⁻¹) above this height, as also exploited by other authors (Fontenla et al. 1991; Mok et al. 2005). The initial density is then derived by assuming hydrostatic equilibrium with the number density of 1.2 × 10¹⁵ cm⁻³ at the bottom and the initial velocity field of all plasma is static. Since the corona needs to achieve a self-consistent thermal structure, we employ a background heating rate decaying exponentially with height into the whole system all the time,

$$\begin{eqnarray}&&{H}_{0}={c}_{0}\mathrm{exp}\left(-\displaystyle \frac{y}{{\lambda }_{0}}\right),\end{eqnarray} \tag{ 6 }$$

where ${c}_{0}={10}^{-4}$ erg cm⁻³ s⁻¹ and ${\lambda }_{0}=50$ Mm. This heating is meant to balance the radiative losses and heat conduction related losses of the corona in its steady state. The slight difference in heating scale height between 50 Mm in Equation (6) above and 36 Mm in Equation (2) in Fang et al. (2013) improves numerical stability at the top boundary and prevents it from cooling down during the longer timescale run performed here. With the above initial setup the whole system is now out of thermal equilibrium. We integrate the governing equations in time with heating $H={H}_{0}$ active until the system achieves a quasi-equilibrium state. After 72 minutes, the above configuration reaches a quasi-equilibrium state shown in Figure 1, which represents a 3D impression of the numerical box quantifying the temperature and number density profile and selected magnetic field lines. The t = 0 in Figure 1 means that after reaching the quasi-equilibrium state, we reset the time of the system back to zero for the next stage of simulation. As seen in Figure 1, the numerical relaxation phase leads to some thermodynamic structuring in the final arcade. Some chromospheric plasma is quickly evaporated into coronal loops at the beginning of the relaxation, but this material gradually loses its kinetic energy. As a result, the final relaxed state of the system is identified as the time when the maximal residual velocity in the simulation is less than 5 km s⁻¹. In this end state, Figure 1 shows a relatively thin TR located at heights between 3 and 5 Mm, which connects the chromosphere to corona. This TR is higher above the neutral line, due to less downward thermal flux there because of the strong horizontal magnetic field. The plasma beta is 0.06 at 20 Mm height above the neutral line while the temperature and number density there are around 1.7 MK and $3.5\times {10}^{8}\;{\mathrm{cm}}^{-3}$. The total mass (per unit length in the ignored dimension) of hot plasma in the corona is around $3.2\times {10}^{4}$ g cm⁻¹.

Zoom In Zoom Out Reset image size

Following this equilibrated system, we turn on a relatively strong heating H₁. This extra heating is localized near footpoints in the chromosphere with formula (Fang et al. 2013):

$$\begin{eqnarray}{H}_{1}=\left\{\begin{array}{ll}{c}_{1} & \mathrm{if}\;y\;\lt \;{y}_{c}\\ \mathrm{and}\quad A({x}_{1},0)\;\lt \;A(x,y)\;\lt \;A({x}_{2},0) & \\ {c}_{1}\mathrm{exp}(-{(y-{y}_{c})}^{2}/{\lambda }^{2}) & \mathrm{if}\;y\;\geqslant \;{y}_{c}\\ \mathrm{and}\quad A({x}_{1},0)\;\lt \;A(x,y)\;\lt \;A({x}_{2},0) & \end{array}\right.\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&A(x,y)=\displaystyle \frac{{B}_{0}{L}_{0}}{\pi }\mathrm{cos}\left(\displaystyle \frac{\pi x}{{L}_{0}}\right)\mathrm{exp}\left(-\displaystyle \frac{\pi y\mathrm{sin}{\theta }_{0}}{{L}_{0}}\right),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\lambda }^{2}=\displaystyle \frac{a\left(A(x,y)-A({x}_{2},0)\right)}{A({x}_{2},0)-A({x}_{1},0)}+b\quad ({\mathrm{Mm}}^{2}),\end{eqnarray} \tag{ 9 }$$

where ${c}_{1}={10}^{-2}$ erg cm⁻³ s⁻¹, ${y}_{c}=3$ Mm, ${x}_{1}=26$ Mm, ${x}_{2}=14$ Mm, a = 0.8 Mm², and b = 1.2 Mm². This choice of strong base heating contrast (${c}_{1}/{c}_{0}=100$) between H₁ and H₀ can mimic flare related chromospheric evaporation. However, this H₁ heating decreases with height to very small values and reaches one-tenth of H₀ heating at 10 Mm. As a result, the H₁ heating dominates heating in the chromosphere and the TR, while the H₀ heating plays a more important role in the heating in the corona. The parameter ${y}_{c}=3$ Mm represents the height of the TR in the quasi-equilibrium system. $A(x,y)$ is the magnetic potential depending on the location and decaying exponentially with height into the whole system. Because the magnetic potential along a single magnetic field line is constant, we add extra heating at both feet of all magnetic field lines identified by $A(x,y)$ in the range of ${x}_{1}\lt | x| \lt {x}_{2}$. Since catastrophic cooling is very sensitive to the heating decay scale and the length of magnetic field lines (Xia et al. 2011), the heating decay scale λ is set to larger values for longer field lines by the above formulae.

We use the MPI-parallelized Adaptive Mesh Refinement (AMR) Versatile Advection Code MPI-AMRVAC (Keppens et al. 2012; Keppens & Porth 2014; Porth et al. 2014) to run the simulation. An effective resolution of 4096 × 2560 or an equivalent spatial resolution of 20 km in both directions is obtained through six AMR levels. This represents an effective fourfold improvement in resolution with respect to our earlier model (Fang et al. 2013). Our numerical strategy to advance the governing partial differential equations uses a three-step Runge–Kutta type scheme. For flux computations, a third-order-accurate limited reconstruction (Čada & Torrilhon 2009) is introduced to calculate the variable evaluation from cell center to cell edge. We adopt a suitably mixed prescription between a diffusive total variation diminishing Lax–Friedrichs and contact-resolving Harten–Lax-van Leer with contact restored (HLLC) scheme (Meliani et al. 2008).

For the boundary treatment, we employ two grid layers exterior to the domain as ghost cells to prescribe cell center values. Considering the left and right physical boundary, density, energy, y and z momentum components, B_y and B_z are set symmetrically, while v_x and B_x are adopted antisymmetrically to ensure zero face values. In the bottom boundary ghost cells, we use the primitive variables ($\rho ,{\boldsymbol{v}},p,{\boldsymbol{B}}$) to set all velocity components antisymmetrically to enforce both no-flow-through (vertical) and no-slip (horizontal), while the ${\boldsymbol{B}}$ are fixed to the initial analytic expressions of Equation (5), and the stratification of density is kept at pre-determined values from the initial condition, as well as the pressure. We always resolve the bottom region up to y = 0.5 Mm at the maximum resolution. For the top conditions, we set all velocity components as antisymmetric, and adopt a discrete pressure-density extrapolation from the top layer pressure with a maximal temperature ${T}_{\mathrm{top}}=2\times {10}^{6}$ K. For the magnetic field, we use a two-cell zero-gradient extrapolation to determine ${\boldsymbol{B}}$ in the ghost cells and improve B_y from a second order one-sided centered difference evaluation of ${\rm{\nabla }}\cdot {\boldsymbol{B}}=0$.

The forming process of the first condensation in our 2.5D simulation is shown in Figure 2, which presents the temporal evolution of number density (left columns), temperature (middle columns), and gas pressure (right columns) at t ≈ 101.2, 101.5, and 102.2 minutes. When we compare these results with the corresponding Figures 5 and 7 of 1D hydrodynamic simulations in Xia et al. (2011), we conclude that all three parameters behave similarly in the forming process, as the number density increases rapidly from 10⁸ cm⁻³ to 10¹⁰ cm⁻³, while the temperature decreases down to 0.01 MK. Along each arched field line, this is analogous to the sudden thermal instability onset in 1D runs. This similarity confirms the applicability of restricted 1D model efforts that assume a rigid 1D loop under the prevailing plasma β conditions, which takes on a local value of around 0.06. The middle panel in the right column of Figure 2 also shows significantly increased gas pressure inside the condensation and a layer of low gas pressure surrounding it after its formation. In the bottom panels of Figure 2, we notice that density, temperature, and gas pressure all reveal a front propagating as expanding wings on both sides of the condensation. This phenomenon is because fast siphon inflows are driven into the forming condensation by a strong pressure gradient between the lower gas pressure around the condensation and relatively higher gas pressure away from the condensation, as seen in the middle panel in the right column of Figure 2. These two siphon flows meet up with the blobs and dynamically impact on the blob to generate two rebound shocks. Hence, while thermal instability and runaway cooling triggers a growing condensation, one also forms two rebound shock fronts that propagate away from the blob. The slightly different formation times at different parts of the condensation on adjacent magnetic field lines (Fang et al. 2013), which are due to gradual variations in length and chromospheric footpoint conditions, are the reason that these two expansion shock fronts display a fan-shaped structure, forming earliest in the blob center and spreading away from the blob. This fan-shaped structure of the rebound shocks is also clearly observed in Xia et al. (2012).

Zoom In Zoom Out Reset image size

However, not every condensation realizes this nearly left-right symmetric situation as seen near the loop apex for this first condensation from Figure 2. Due to slightly asymmetric conditions already prevailing after the numerical relaxation process and due to perturbations from existing condensations, most of the following condensations initiate in loop limbs (also shown in the online movies of Fang et al. 2013). The field-projected gravity force on the limbs leads to asymmetric plasma distributions, as seen, e.g., in the number density map in panel (a) of Figure 3 at $t\approx 113.3$ minutes, the moment when local catastrophic cooling begins there (about 10 minutes after the first condensation). The higher central gas pressure indicates the initial forming location of this condensation in panel (b) of Figure 3. Due to its limb-loop location, the number density map points out that the right side of the condensation holds a relatively denser (3 × 10⁹ cm⁻³) and wider plasma distribution than the left ($1.5\times {10}^{9}$ cm⁻³). Still, strong pressure gradients drive siphon flows from both sides toward this condensation. After a short time at $t\approx 113.7$ minutes, the denser and heavier plasma at the right of this condensation has a (left-directed) siphon flow with a slower speed of 23 km s⁻¹, compared to the left siphon flow (which is right-directed) at a speed of 42 km s⁻¹, shown by the velocity magnitude map in panel (c) of Figure 3. As discussed above, the impact of siphon flows on the condensation naturally generates rebound shocks, whose speeds are determined by the original speeds of the siphon flows and the mass contrast between the condensation and the siphon flows. The slower and heavier siphon flow on the right of the blob here leads to a much slower rebound (right-directed) shock seen to separate at 7 km s⁻¹, while the left one (left-directed) travels at 21 km s⁻¹. These two rebound shocks are identifiable in the gas pressure map in panel (e) of Figure 3 at $t\approx 114.8$ minutes. The condensation itself has a velocity of 5 km s⁻¹, meaning that basically the right rebound shock barely can sweep up and heat little siphon flow plasma. Because the central condensation keeps sucking in plasma from nearby and the rebound shock at the right of the blob is too slow to sweep and heat up plasma, the gas pressure there does not rise to a higher value and keeps a strong pressure gradient at the right of the blob, as shown in panel (e) of Figure 3. About 3 minutes later at $t\approx 117.6$ minutes, this persistent pressure gradient at the right of the blob accelerates the left-directed siphon flow to a higher speed of 52 km s⁻¹ (shown in panel (i) of Figure 3); therefore, the corresponding right-directed rebound shock finally speeds up to 28 km s⁻¹ and begins to sweep and shock-heat the plasma on its way. In short, initial asymmetric situations on the condensation can lead to a complicated thermal and dynamical evolution and result in a delay of rebound shocks spreading at one side of the condensation.

Zoom In Zoom Out Reset image size

Additionally, we also find another special case, namely blob A in Figure 4, which has only one rebound shock on its left side. Figure 4 shows the gas pressure map (a) and (b) at $t\;\approx \;134.8$ and 137.6 minutes, with a dotted isocontour of the number density at 7 × 10⁹ cm⁻³ overplotted. This density contour at 7 × 10⁹ cm⁻³ is one of the criteria that identifies whether a cell contains cool plasma belonging to coronal rain, as used later on. Panel a in Figure 4 indicates a similar situation for blob A as in the second row of Figure 3 in which only the left rebound shock spreads out. In contrast to what happens in the third row of Figure 3, for blob A we do not find a right rebound shock in Figure 4 until the collision and merging of blob A with blob B. The reason is that when the thermal instability triggers the condensation labeled there as blob A, another existing condensation labeled as blob B in the same coronal loop has already depleted the plasma between these two blobs. Therefore the small pressure gradient in the emptied loop between the two blobs cannot drive a fast siphon flow to create a strong rebound shock for blob A, even though the gas pressure on the right of blob A is low enough (panel (a) in Figure 4). Afterwards, when blob A catches up and merges with blob B because of the strong pressure gradient outside these two blobs, the rebound shock at the right side of blob A is still not fast enough to show clear separation and propagation.

Zoom In Zoom Out Reset image size

We also observe the details of a gas pressure substructure within these shock-bounded regions around the condensation in the simulations. These reveal the establishment of a PCTR-like structure around all blobs. The gas pressure substructure around the first condensation consists of three components shown in panel (a) of Figure 5 and panel (i) of Figure 2, namely a high gas pressure outside of the condensation, a low gas pressure at the boundary of the condensation, and a higher gas pressure in the center of the condensation. Actually, it is not only this first condensation in Figure 2 that has this kind of gas pressure substructure, but all the blobs which establish a dynamic equilibrium around themselves also have it, e.g., all the blobs in Figures 3 and 4. To better quantify this, we identify a field line crossing the center of the blob shown in panel (a) of Figure 5 and plot gas pressure, temperature and radiative loss along this field line in panel (c) of Figure 5. The temperature declines from a coronal temperature of 0.35 MK to a cool plasma temperature of 0.01 MK in 200 km and density increases from 1 × 10⁸ cm⁻³ to 1 × 10¹⁰ cm⁻³; therefore this 200 km area could basically be considered a PCTR. Within this area, we find that two highly radiative loss peaks exist, introduced by a temperature around 0.02 MK. This corresponds to the two dips of gas pressure at the boundary of the blob. These two strong radiation areas also indicate the location in which catastrophic cooling takes place, ensuring that the condensation keeps growing. Indeed, the two dips in gas pressure always relocate with the boundary of the blobs, coincident with the strong emissive loss. Although the temperature of 0.01 MK inside the blob is lower than in the surrounding coronal plasma, a much higher density at the center of the condensation (5 × 10¹⁰ cm⁻³) leads to a slightly higher gas pressure there. The high gas pressure outside of the condensation reflects the post shock conditions prevailing there after the rebound shocks run against the condensation inflows. Note that our resolution is such that we have about seven grid points along the field line through the PCTR at each side of the blob in Figure 5, clearly resolving the PCTR around the blob in our simulation.

Zoom In Zoom Out Reset image size

The difference in gas pressure between the inside and outside of the condensation is found to persist throughout the lifetime of the blobs and plays a role in the movement of the blobs. Especially when the blobs fall along the field lines toward footpoints, the gas pressure and temperature ahead of the descending blob increase as shown in panel (b) of Figure 5 due to the blob compressing the plasma ahead of it in the loop and the strong evaporation at the loop footpoints. We also identify a field line crossing the center of the blob shown in panel (b) of Figure 5 and plot gas pressure, temperature, and radiation loss along this field line in panel (d) of Figure 5, which also shows an obvious PCTR. Due to the gravity variation and the strong gas pressure gradient between the two sides of the blob, the lower part of this blob has a higher density distribution, which naturally leads to a higher radiative loss. This strong gas pressure gradient slows down the acceleration of the blob in its descent. This was also pointed out in Fang et al. (2013), where we stated that sometimes, it can even lift lighter blobs to cross the loop apex.

Panel (a) of Figure 6 shows the temporal evolution of the total mass of cool (solid) and hot (dashed) plasma in the corona, and panel (b) of Figure 6 shows the number of blobs for the entire time interval of our 2.5D simulation. The criteria to identify whether a cell contains cool plasma belonging to coronal rain are that (i) the number density is higher than $7\times {10}^{9}$ cm⁻³, (ii) the temperature is lower than $2\times {10}^{4}$ K, and (iii) the location is above the chromosphere-corona-transition-region. We dynamically locate the height of the TR at each x-position as ${y}^{\mathrm{tr}}(x,t)$ by searching the vertical position of the (first) maximum gradient value of temperature from the bottom boundary. Each blob is defined as a collection of neighboring cells which hold cool plasma. However, if the number of grid cells in one blob is smaller than 10 at our highest resolution, we remove this blob from the blob list to avoid counting spurious transient features that do not collect into a clearly resolvable blob, and also to mimic the observational resolution. As stated before, we adopt a four times higher resolution than in Fang et al. (2013), but also extend the simulation to a two times longer time of around 370 minutes (previously 190 minutes). By running our 2.5D simulation for these much longer times, we find that the whole coronal rain process shows limit cycles, which has been discussed in earlier 1D simulations (Müller et al. 2003), as well as in observational work (Antolin & Rouppe van der Voort 2012). This is the first time that we can report limit cycles of coronal rain in a multidimensional simulation, which confirms that constant heating conditions that provide enough energy can form secondary (or even more) coronal rain cycles in a single arcade. From panels (a) and (b) of Figure 6, we find the time interval between the first and secondary cycle to be around 175 minutes, when measured between successive maxima in cool mass matter. Panel (a) of Figure 6 shows the temporal evolution of total mass of hot coronal plasma which is the mass in the corona, excluding the cool plasma identified by the above criteria. We find that at $t\approx 140$ minutes, the total mass of cool plasma reaches its peak in the first cycle in panel (a) of Figure 6, while at $t\approx 143$ minutes the catastrophic cooling process has cooled down most of the hot plasma in the corona shown in panel (a) of Figure 6. From about $t\approx 133$ minutes, blobs begin to fall into the TR, then the evaporation of plasma in the chromosphere driven by the extra heating H₁ fills the evacuated loops left by blobs that have already sank into the chromosphere. From this moment, until the onset of the secondary cycle of our coronal rain shower at $t\approx 250$ minutes, it takes about 120 minutes, which is of similar duration to the time for the first cycle to reach its onset (about 100 minutes). So although we infer from the total mass evolution of cool plasma in panel (a) that there is only about 50 minutes between the ending of the first and the beginning of the second cycle, actually the continued heating at the chromosphere has already spontaneously began to fill the empty corona 70 minutes before. We also see that the total mass of hot plasma before the onset of the secondary cycle is higher than in the first cycle (panel (a)), which leads to a longer lasting secondary cycle with more mass in condensations. Panel (a) of Figure 6 indicates that at $t\approx 130$ minutes (before the first blob falls into the chromosphere), there is at least 9 × 10³ g cm⁻¹ of cool plasma in the corona, which originally was hot plasma. Meanwhile panel (a) of Figure 6 also suggests that compared with the corona before the onset of catastrophic cooling at $t\approx 100$ minutes, the decrement in the same time of total mass of hot plasma at $t\approx 140$ minutes is only 5 × 10³ g cm⁻¹. The difference between the increase in cool plasma and the decrease in hot plasma indicates that during these 30 minutes since onset at $t\approx 100$ minutes, the evaporation in the chromosphere evaporates 4 × 10³ g cm⁻¹ into the corona, i.e., at an evaporation rate of 2.2 g cm⁻¹ s⁻¹. We can similarly estimate an evaporation rate of 2.3 g cm⁻¹ s⁻¹ between the onset of the secondary cycle and the moment its first blob falls into the TR. Until the onset of the first cycle at $t\approx 100$ minutes, the increment of total hot plasma from turning on the extra heating H₁ is about $13.2\times {10}^{3}$ g cm⁻¹ in total, further confirming this evaporation rate of 2.2 g cm⁻¹ s⁻¹. Based on these estimates, we infer that anywhere in both simulated cycles the constant extra heating H₁ leads to a nearly constant evaporation rate. We can thus extrapolate to even more cycles expected further on, and interpret these limit cycles as a chronological sequence of mass recycling between chromosphere and corona: heating in the chromosphere brings plasma to the corona by evaporation, where it ultimately triggers catastrophic cooling; the cooling process manages the coronal plasma into coronal rain where the plasma drains back to the chromosphere; and persistent heating causes the chromospheric material to evaporate again toward the corona.

Zoom In Zoom Out Reset image size

Although the duration and peak value of the total mass in both computed cycles are similar, their initial condensation rates (in contrast to the previously discussed evaporation rate) computed from the temporal variation of their total mass curve work out to be 6.7 and 4.5 g cm⁻¹ s⁻¹, respectively, and thus are different. It is known from linear thermal instability theory (Field 1965) and 1D simulation results in Xia et al. (2011) that this initial condensation rate in catastrophic cooling depends on parameters controlling the energy input from heating. One notices that the condensation rate (the local derivative of the solid curve in panel (a) of Figure 6) varies dramatically even within one cycle, despite a constant heating energy input in our multidimensional simulation. We now will interpret the reason for the changes seen in the condensation rate by surveying especially the process of growth for the first condensation, which forms under a relatively simple and almost symmetric condition.

The solid line in panel (a) of Figure 7 shows the temporal evolution of the mass accumulation for this first condensation (the one from Figure 2) from $t\approx 100$ to 110 minutes. Its nearly linear behavior quantifies that the condensation rate remains almost constant in this time interval at a value of about 2.3 g cm⁻¹ s⁻¹. We deliberately do not discuss what happens to the first condensation after $t\approx 110$ minutes, since afterwards it breaks into two smaller blobs. In the same figure, panel (a), the dashed line displays the growth of the total mass of cool plasma as seen on a single field line through the center of the first condensation, i.e., in a 1D fashion. To show this, we identify the group of grid points that are passed by the field line. The total mass of cool plasma determined on the single field line keeps growing in time, but its growth rate is much smaller than that for the whole 2D condensation. Panel (b) of Figure 7 quantifies the temporal evolution of typical lengths for the first condensation, where we quantify both the length parallel to the magnetic field lines and the length perpendicular to the magnetic field lines. This indicates that blob growth in the perpendicular direction is much faster than in the parallel one, which can be seen visually as well in all columns in Figure 2 and Figure 3. As discussed in Fang et al. (2013), the low pressure region surrounding the first condensation onset leads to magnetic restoring forces on adjacent loops. These in turn influence the location where the catastrophic cooling will take place on the adjacent loops, which are all close to the thermal instability onset. The different growth rates found for blob sizes in these two directions then relate to the fairly fast "growth" along the perpendicular direction due to sympathetic runaway cooling onset, versus the slower growth seen in the parallel one, which is the only one found in 1D setups. The average density of each cell of the condensation is quantified in panel (c) of Figure 7, and this density stays basically the same in the forming process, meaning that the total mass of the condensation is just proportional to the increasing number of neighboring grid cells that contain cool plasma. While the number of cells in the condensation increases in both directions, the larger condensation rate of the whole blob in panel (a) of Figure 7 again directly reflects the faster growth in size in the perpendicular direction. We conclude that the growth of the total mass of individual blobs in our simulation is mainly determined by the onset of catastrophic cooling in neighboring loops rather than the growth along the loops in which catastrophic cooling gets triggered. We can indeed verify this 2D growth aspect by further showing a correlation between the total mass of cool plasma and two other measures, which holds up for an even longer time than the first 10 minutes, i.e., when several blobs have started to form. This is shown in panel (d) in Figure 7 where we plot the temporal evolution of the total mass of cool plasma, the size of the onset TR, and the total blob region width. The total blob region width indicates the total width of all magnetic loops where catastrophic cooling takes place on. The size of the onset TR means the corresponding width as found at the TR height, of all the loops undergoing catastrophic cooling. Because of the magnetic arcade configuration adopted, these size measures for the affected loops give higher values for higher locations, i.e., the total blob region width always exceeds the (field aligned remapped) onset TR size. The latter size of the onset TR shows a nice correlation with the total mass of the cool plasma evolution.

Zoom In Zoom Out Reset image size

In the simulation, we observe plenty of blobs hitting the TR and disappearing into the lower chromosphere, which is also known to occur in observations (Antolin et al. 2010; Antolin & Rouppe van der Voort 2012). Tripathi et al. (2009) observed high-speed downflows and concurrent upflows in coronal loops close to the footpoints and argue in favor of upflows in coronal loops at higher temperatures. Antolin et al. (2010) confirmed that the high-speed downflows represent the cool plasma, which correspond to the falling blobs in our simulation (see also the movie in Fang et al. 2013). Meanwhile, our 2.5D simulation also shows the possibility of triggering concurrent upflows as observed by Tripathi et al. (2009) and Kleint et al. (2014). Panel (a) in Figure 8 shows the number density map at $t\approx 143.7$ at a moment when falling blobs sink into the TR and compress the plasma on its way (at about $x\approx -2.1$ Mm). Panel (b) in Figure 8 shows the vertical velocity map at the same location and instant, which clearly displays the concurrent upflows rising at the tails of the declining blobs in the same field line bundles. Hence, this answers the question in Kleint et al. (2014) whether the upflows can flow along the same field lines as the downflows. These upflows in our simulation are actually rebound shocks from the impact of the blobs on the TR. They arise immediately when the blobs impact on the TR and spread from one footpoint to another footpoint in around 5 minutes with a velocity of around 50 km s⁻¹. From Panel (a) we can see the enhanced density left after the passage of these rebound shocks. However, panel (c) in Figure 8, which quantifies temperature, indicates that the temperature in the loop already increases before the rebound shocks have reached far into the loop, since the parametrized background heating H₀ very efficiently heats the low density loops left by falling blobs. Panel d also shows the temperature, but now on a larger domain and at a later time. It shows that afterwards the rebound shocks heat the low density loops to an even higher temperature of 2.0 MK. After the rebound shocks reach the other footpoint, the loops are at high temperature of about 2 MK but with a low number density of 1 × 10⁸ cm⁻³. We distinguish this from further upflows coming from evaporation due to the extra strong heating H₁ located in the chromosphere. This enhances the density to 1 × 10⁹ cm⁻³ again and the temperature to 2.3 MK. However, these upflows from evaporation rise with a much slower velocity of around 15 km s⁻¹.

Zoom In Zoom Out Reset image size

To quantify even further the detailed fate of a blob when it hits and descends into the TR, Figure 9 shows the temporal evolution of the mass, density, velocity, kinetic energy, momentum, and temperature of the first coronal rain blob to hit the TR from the corona and to sink down into the chromosphere. The vertical dashed line in each panel of Figure 9 points at $t\approx 132$ minutes when this blob hits the TR. Because the density and temperature of plasma in the TR is comparable with those of the blobs, we can no longer use only the density and temperature as criteria to distinguish blobs when they are near or partially below the TR. In order to identify plasma belonging to the blob as it hits and descends in the TR after $t\approx 132$ minutes, we change our criteria to require the local velocity to be larger than 3 km s⁻¹ and the locations are below the TR line ${y}^{\mathrm{tr}}(x,t)$ after $t\approx 132$ minutes. Since the velocity of the plasma in the TR is almost zero, this velocity-based criterion captures the location of sinking blobs. In panel (a) of Figure 9, we find that the mass identified as blob material by the above criteria begins to increase at t ≈ 132.7 minutes. This is because the mass detected not only includes the blob itself, but also counts mass compressed and accelerated by the blob impact. At $t\approx 136$ minutes, the total mass affected reaches its peak at six times the original blob mass. After $t\approx 136$ minutes, due to the combined influence of reflection–transmission processes at the TR and the higher gas pressure from the impact, the velocities in much of the blob-impacted area decrease to values smaller than the criterion 3 km s⁻¹. This is then seen as a mass decrease in our panel (a). In panel (b) of Figure 9, the density versus time profile keeps rising while the blob hits the TR. As we know, this blob impact compresses the TR plasma swept up by the blob and transfers momentum from the sinking blob to the impacted plasma, and therefore in panel (c) of Figure 9 we find that the average velocity of the region identified keeps decreasing during the whole process, as well as the kinetic energy shown in panel (d). Panel (e) of Figure 9 shows the total momentum of the mass identified. Due to the gravitational acceleration, the blob momentum keeps increasing until it reaches its maximum value at t ≈ 136 minutes, then it reduces quickly. This is a combination of the mass evolution in panel (a) and the velocity information from panel (c). The momentum and velocity decreasing after the impact relate to momentum transfer to the surrounding TR and upper chromosphere plasma, until the regions selected by the velocity-based criteria vanishes: the local conditions settle to static chromosphere conditions. Panel (f) of Figure 9 shows the average temperature evolution during the blob impact. The temperature increases before hitting the TR due to the compressional heating when the blob descends through the higher gas pressure region just above the TR. After the impact, since also more cooler material gets identified as impacted matter, one settles back to upper chromospheric temperature values.

Zoom In Zoom Out Reset image size

The impact speed of blobs in Figure 9 is around 30 km s⁻¹, and the highest impact speed of all blobs in our simulation is around 60 km s⁻¹ and number densities range from 4 to 6 × 10¹⁰ cm⁻³. Our maximum impact speed is much lower than the falling speeds reported in Kleint et al. (2014) which went up to 200 km s⁻¹. They report that these coronal rain events with high impact speeds are correlated with local brightenings which probably indicate an increase of density and temperature in the TR. Panel (b) of Figure 9 and panel (a) of Figure 8 confirm the expected increase of the number density of impacts in our simulation.

We also find another interesting phenomenon in our numerical simulation, namely the self-consistent establishment of counter-streaming flows. Such anti-parallel flows are very commonly found in solar observations, especially in prominences (Alexander et al. 2013). Panels (a), (b), and (c) in Figure 10 respectively show the signed velocity magnitude map (with the sign taken from the horizontal velocity component), the gas pressure map, and the number density map at $t\approx 132.6$ minutes. Panel (d) shows the signed velocity magnitude map as in panel (a), but at a later time, namely at $t\approx 146.9$ minutes. These four panels in Figure 10 display many cases of counter-streaming flows established on neighboring field line bundles in our simulation and allow the origins of counter-streaming flows to be explained. After thermal instability inducing a runaway catastrophic cooling and initial growth in an almost static state, the condensations lose their delicate force balance and begin to slide toward one footpoint along magnetic field lines. Whether a particular condensation segment slides to the left or right is influenced by its initial location and local total force balance (gravity, gas pressure gradient, and magnetic field force). Once in motion, they are accelerated by the field-projected gravitational force. Meanwhile catastrophic cooling keeps taking place around the condensations. As discussed in Section 3.1, the initial catastrophic cooling process depletes the local plasma and sucks in fast inflows, then the spontaneous rebound shocks heat the plasma and increase the gas pressure. Afterwards, no stronger inflows are driven again due to the increased gas pressure. However, there can be several blob pairs lying in the same or neighboring field line bundles, as shown in panel (c) of Figure 10. Figure 11 shows gas pressure maps with magnetic field lines at $t\approx 109.4$ and 113.0 minutes. The black contour relates to the temperature distribution and is an isocontour at 0.1 MK. Both the gas pressure and temperature in panel (a) in Figure 11 indicate the clear PCTR around the blob as previously discussed in Figure 5. After 3 minutes, panel (b) of Figure 11 shows two low (white) pressure sections after the blob breaks into three segments. These low pressure sections slant through the field lines and they are the elongated PCTR cross sections from the original parts of the whole blob in panel (a) of Figure 11. Because of the strong radiation in the PCTR, the temperatures inside these elongated regions remain low during their deformation. As a result, we could consider these cross sections to undergo isothermal expansion. Because the condensed mass in these narrow regions grows much slower than their areal growth due to elongation, the densities inside these elongated cross sections decrease faster as well as the gas pressure. This leads to blob sequences with low pressure sections in between them. This is also seen in panel (b) of Figure 10 where a sequence of blobs show up with low (white) pressure sections in between them. The depleted areas trigger siphon inflows that refill these regions. Then this pair of siphoned fast inflows establish the counter-streaming flows between the pair of neighboring blobs. Panels (a) and (b) in Figure 10 also show that falling condensations with larger velocities induce larger density depletions and lower gas pressure areas on their way down, which lead to faster inflows than those found for more static condensations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Panel (d) in Figure 10 indicates a different origin of counter-streaming flows at $t\approx 146.9$ minutes. As we discuss in Section 3.3, we observe that after blobs decline into the TR, concurrent upflows rise up toward the loop apex. Upflows labeled A in panel (d) in Figure 10 are the concurrent upflows shown in panel (b) of Figure 8, but about 3 minutes later (concurrent with the later temperature panel (d) of Figure 8). Upflows arising from blob impacts also have the chance to establish a counter-streaming flow if there is an opposite flow pattern in the neighboring loops. The difference between these two different origins for counter-streaming flows is that the one based on depleted sections between a pair of blobs lying on neighboring field line bundles can last through the whole falling process of blobs, or on timescales of about 10 minutes, while the other ones will vanish after the upflows refill the loops, typically in a shorter timescale of about 5 minutes.

The sheared flows that are established by the detailed blob dynamics could also in return influence the further evolution of the condensations. An example is shown in panel (a) of Figure 12, where we show a signed velocity map, with overlaid contours of the density distribution of condensations at levels of 7, 25, and 50 ×10⁹ cm⁻³ at $t\approx 123.7$ minutes. Concentrating on the density feature labeled A, after its initial formation, sheared flows already begin to take shape. About 10 minutes later, this segment A is seen as segment A1 and A2 in panel (b) of Figure 12, and the condensation has broken into two distinct segments with increasing separation between them. Segment A2 is also going to break into two segments a little later. At t ≈ 123.7 minutes in panel (a) of Figure 12, this segment A feature is more like one whole elongated condensation. However, by $t\approx 132.6$ minutes in panel (b) in Figure 12, several condensations behave totally separately from each other. Another example is the one of segment B in panel (a) and panel (b). This breaks up into segment B1 and segment B2 in panel (c) at t ≈ 136.2 minutes. Then segment B1 further breaks into segment B1 and B3 in panel (d) at t ≈ 139.8 minutes. This gradual change from one elongated dense blob or strand breaking into several segments, surrounded by fast sheared flows, hints at the influence of Kelvin–Helmholtz instabilities (KHI). However, there is no clear vorticity pattern emerging in our simulation that would clearly identify KHI development, which may not have enough time to develop. We speculate that other KHI-related substructure may well arise under different parameter settings (field strength, heating scale height), but already establish that sheared flows contribute to the breaking up of elongated condensations into smaller fragments.

Zoom In Zoom Out Reset image size