When hot meets cold: post-flare coronal rain

Most solar flares demonstrate a prolonged, hourlong post-flare (or gradual) phase, characterized by arcade-like, post-flare loops (PFLs) visible in many extreme ultraviolet (EUV) passbands. These coronal loops are filled with hot -- $\sim 30 \,\mathrm{MK}$ -- and dense plasma, evaporated from the chromosphere during the impulsive phase of the flare, and they very gradually recover to normal coronal density and temperature conditions. During this gradual cooling down to $\sim 1 \,\mathrm{MK}$ regimes, much cooler -- $\sim 0.01 \,\mathrm{MK}$ -- and denser coronal rain is frequently observed inside PFLs. Understanding PFL dynamics in this long-duration, gradual phase is crucial to the entire corona-chromosphere mass and energy cycle. Here we report a simulation in which a solar flare evolves from pre-flare, over impulsive phase all the way into its gradual phase, which successfully reproduces post-flare coronal rain. This rain results from catastrophic cooling caused by thermal instability, and we analyse the entire mass and energy budget evolution driving this sudden condensation phenomenon. We find that the runaway cooling and rain formation also induces the appearance of dark post-flare loop systems, as observed in EUV channels. We confirm and augment earlier observational findings, suggesting that thermal conduction and radiative losses alternately dominate the cooling of PFLs.


INTRODUCTION
Solar flares represent explosive phenomena in the solar atmosphere, where 10 28 −10 32 ergs of energy originally stored in the solar magnetic field can suddenly be released via magnetic reconnection (Sweet 1958;Shibata & Magara 2011). The time development of a solar flare event can be divided into three phases: a preflare phase, a sudden, impulsive phase and a gradual or post-flare phase (Kane 1974). The magnetic energy is released rapidly in the impulsive phase within a typical timescale of (tens of) minutes. A large fraction of this released energy is transported from the tenuous and hot corona downwards to the denser and colder solar chromosphere via thermal conduction and by energetic electrons (Antiochos & Sturrock 1978;Ruan et al. 2020). This deposition of energy in the chromosphere leads to a sudden heating of the local plasma, and causes upward evaporation of the plasma to form super hot (∼10 MK) and dense (∼ 10 10 cm −3 ) arcade-like loop systems at coronal heights. The loops return to their usual coronal conditions (∼1 MK, ∼ 10 8 − 10 9 cm −3 ) in the following, gradual phase, where field-guided thermal conduction and radiative losses generally contribute to the cooling process (Cargill et al. 1995;Aschwanden & Alexander 2001). These loops, visible in extreme ultraviolet (EUV), but also in Hα images, are usually called post-flare loops (PFLs) (Bruzek 1964).
Thanks to dramatically increased spatio-temporal resolutions in observations, this gradual phase of solar flares is now known to show PFLs which spontaneously develop fine-scale coronal rain (Scullion et al. 2014;Martínez Oliveros et al. 2014;Jing et al. 2016;Scullion et al. 2016). In the multi-thermal coronal rain events, cool and dense rain blobs form in-situ in the hot corona, to fall to the chromosphere with speeds up to 100 km s −1 . Coronal rain in flaring coronal loops has been reproduced in 3D simulations by Cheung et al. (2019) and Chen et al. (2021), although the dynamics of this rain has not been analyzed in any detail. Coronal rain is also observed in non-flaring coronal loops, is frequently found in loops of active regions (Leroy 1972;Levine & Withbroe 1977;Schrijver 2001;O'Shea et al. 2007; Antolin & Rouppe van der Voort 2012; Ahn et al. 2014;Antolin et al. 2015), and this type of coronal rain has been studied previously using magnetohydrodynamic (MHD) simulations (Fang et al. 2013(Fang et al. , 2015aMoschou et al. 2015;Xia et al. 2017;Kohutova et al. 2020). Kohutova et al. (2020) successfully reproduces coronal rain in a self-consistent 3D radiative MHD simulation, in which the radiative loss in chromosphere and corona is offset by Ohmic and viscous heating. It has been generally accepted that these rain blobs are generated in a catastrophic cooling process, essentially caused by thermal instability (Parker 1953;Field 1965;Claes & Keppens 2019;Claes et al. 2020). In a catastrophic cooling event, local temperatures drop from 1 MK to below 0.1 MK within one minute, while local densities can increase by orders of magnitude (Scullion et al. 2016). Observations of flare-driven coronal rain demonstrate that this catastrophic cooling can also happen in PFLs, but this has thus far never been modeled in detail. Another phenomenon thought to have a close relationship with these sudden condensations are the so-called dark post-flare loops (DPFLs), where some PFL loops suddenly vanish from specific EUV passbands, e.g. at 17.1 nm, 30.4 nm and 21.1 nm (Song et al. 2016;Jejčič et al. 2018;Heinzel et al. 2020). In these DPFLs, an EUV loop which was bright for a while suddenly darkens for several minutes, so effectively disappears between the adjacent EUV loops seen at the same height. It has been suggested that the formation of cool and dense coronal rain may contribute to EUV emission and absorption, inducing DPFL formation (Jejčič et al. 2018;Heinzel et al. 2020). However, this suggestion must still be confirmed in an ab-initio model.
Here we perform an MHD simulation of a flare event from its pre-flare phase all the way into the gradual phase. This simulation finally allows us to understand the complex thermodynamic evolutions of PFLs. We (1) reproduce postflare coronal rain; (2) quantify the chromosphere-corona mass and energy cycles during PFLs; and (3) demonstrate the intricate relationship between condensations and the disappearing EUV loops, or DPFLs.

METHODS
We perform the simulation with the open-source MPI-AMRVAC code (Xia et al. 2018;Keppens et al. 2021). The simulation is 2.5D, where the domain is 2D but all vector quantities have three components. The simulation domain is given by -75 Mm ≤ x ≤ 75 Mm and 0 ≤ y ≤ 100 Mm. This simulation box has an initial resolution of 96 × 64, but an equivalent high resolution of 6144 × 4096 is achieved with our block-adaptive mesh. The governing equations are the magnetohydrodynamic (MHD) equations with effects of gravity, thermal conduction, radiative loss and magnetic field dissipation due to resistivity included, also shown in Ruan et al. (2020) (the source terms related to fast electrons have not been activated here). The new multi-dimensional field-line-based transition region adaptive conduction (TRAC-L) method is adopted to properly handle the chromosphere-corona interaction at affordable resolution (Zhou et al. 2021).
A relaxation has been done to obtain a static pre-flare atmosphere before we perform the flare simulation. A background heating is required to offset the energy losses due to radiative cooling and thermal conduction. Inspired by Ye et al. (2020), this background heating is a function of initial and spatio-temporally evolving values, given by where h tra = 3 Mm, h a = 0.1 Mm, N e indicates electron number density, subscript 0 indicates the t = 0 initial value and Λ(T ) is the radiative cooling curve adopted in our simulations. The cooling curve from Colgan et al. (2008) is used. Many other optically thin cooling curves, such as from Schure et al. (2009), are also available in MPI-AMRVAC. Hermans & Keppens (2021) showed that changing the cooling curve does not truly change whether condensations happen, but only impacts the growthrate of the thermal instability and the morphology of the formed condensations. The initial vertical temperature T 0 (y) profile in Avrett & Loeser (2008) (model C7) is employed. The number density at y = 40 Mm is set to 2 × 10 9 cm −3 and the initial density profile is calculated based on hydrostatic equilibrium. In the relaxation stage, a uniform vertical magnetic field is adopted. A numerically static atmosphere is obtained after a relaxation time corresponding to 3.5 hours. The local number density at y = 40 Mm decreases to about 10 9 cm −3 after this relaxation. The final instantaneous background heating rate of this relaxed stage is saved and then used in the subsequent flare simulation. Such a background heating (∼ 10 −4 to ∼ 10 −3 erg cm −3 s −1 ) is crucial for keeping the temperature profile outside the flare region and the PFLs, but it is not impacting strongly the condensation process. We verified that condensations can still happen without this added background heating, with an additional simulation, which can be found on the website of ERC PROMINENT project (https://erc-prominent.github.io//media/Ruan/).
After the relaxation, the magnetic field configuration is changed. The magnetic configuration from Ye et al. (2020) is employed then, given by where B 0 = 30 G is the new initial magnetic field strength and λ = 10 Mm. Such a configuration allows magnetic reconnection. A three-stage resistivity strategy is then activated. A spatially localized resistivity inside the initial current sheet triggers magnetic reconnection in the first stage. This localized resistivity is given by Mm and t η1 = 31 s. In the second stage, we use an anomalous resistivity given by and v c = 128, 000 km s −1 . The resistivity is set to zero in the third stage t > t η2 to force the flare to enter the gradual phase, where only numerical dissipation happens. This means that the explicit influence of Joule heating on the cooling of PFLs and the triggering of the condensation is neglected. The resistivity strategy used in our first and second stage is similar to that in Yokoyama & Shibata (2001). We employ symmetric boundary conditions for number density, pressure and magnetic field components, while anti-symmetric conditions are employed for velocity components at the left and right boundaries. At the upper and bottom boundaries, density and pressure are fixed to their initial values. The magnetic field components at the bottom boundary are also fixed to the initial values. An anti-symmetric condition is applied for the x-component of the magnetic field at our upper boundary, while the other two components of the magnetic field employ symmetric conditions there. The velocity at the upper boundary is set to zero, as we are not interested in an erupting fluxrope structure here, but focus on the PFLs. Anti-symmetric conditions are employed for the velocity components at the bottom boundary.
The SXR emission is calculated with the method reported in Pinto et al. (2015). The EUV emissions are calculated with the contribution function provided by the CHIANTI database and the optically thin assumption (Del Zanna et al. 2015). Fig. 1 demonstrates the full evolution in magnetic topology, from pre-flare current sheet to the formation of PFLs at the impulsive phase (Fig. 1b), extended to the entire evolution of the PFLs through the gradual phase (Fig. 1c,d).

POST-FLARE LOOP FORMATION AND EVOLUTION
There is a vertical current sheet separating regions of opposite field directions at the beginning of our simulation (Fig. 1a). Magnetic reconnection inside this current sheet produces closed magnetic arcades below and a flux rope above a reconnection site in the impulsive phase (t 8 min). The magnetic energy is converted into heat at the reconnection site and then conducted into the chromosphere. Deposition of the thermal energy there leads to a sudden increase of local pressure and this produces upward evaporation flows. They fill the generated coronal loops with hot plasma (∼10 MK) and increase the electron number density by one order of magnitude in the PFLs (Fig. 1b,f). Thereafter, the flare enters the gradual phase when we have a rapidly decreasing magnetic reconnection speed (t 8 min). The PFL temperature decreases slowly due to thermal conduction and radiative losses in this gradual phase, but suddenly triggers thermal instability near t 35 min. The loop density also decreases in this period, but is still much higher than the external coronal density (Fig. 1c). Energetic electron beam deposition can also lead to chromospheric evaporations, but this mechanism is less efficient in weak flares where the beam energy flux is small (Fisher et al. 1985). Electron beam heating has not been included into our simulation for this reason, though this heating has been successfully reproduced in a 2.5D flare simulation in our previous work (Ruan et al. 2020). The peak flux in our simulation is about 4 × 10 −7 W m 2 when assuming that the loop depth in the third direction is 100 Mm (see Fig. 1i), therefore the simulated flare is a B level flare.
In the impulsive phase, the energy gain of the flare loop is determined by the reconnection process whereas the energy loss is mainly handled by thermal conduction, which overwhelms the radiative cooling by more than one order. Therefore, the flare temperature can be predicted from the coronal magnetic field strength, density and flare loop length or simply from coronal plasma β with the scaling laws provided in Yokoyama & Shibata (1998, 2001 and Takasao et al. (2015). This flaring temperature estimated from equation (1)   Catastrophic cooling driven by thermal instability condenses local plasma in-situ and leads to high density (close to 10 11 cm −3 ) and cold (close to 0.01 MK) structures in the coronal PFLs (Fig. 1d,h). Fig. 1i shows how the SXR flux reaches its peak value at the impulsive phase, to then gradually decrease as the loop temperature drops. The temporal evolution of the total coronal rain material (i.e. T e < 0.1 MK, N e > 10 10 cm −3 and y > 5 Mm) is also illustrated in Fig. 1e. We note that sudden condensations happen in two successive events during the entire simulated period. These are located in different loop systems, with the second rain event appearing at a higher altitude.
Once condensations happen within PFLs, the formed cold and dense plasma structures will likely fall down from coronal heights due to gravity and hence appear as observed coronal rain blobs. This is fully reproduced in our simulation as demonstrated in Fig. 2. Cold plasma is formed at a PFL looptop at the beginning of the runaway condensation (Fig. 2a). Thereafter, the cold structure extends to lower and higher loops, meanwhile sliding down to one side along the magnetic post-flare arcade (Fig. 2b,c). This falling cold plasma gets accelerated to a speed of ∼100 km s −1 by gravity before it enters the chromosphere (Fig. 2d-f). Such a speed is close to that found in coronal rain observations (Martínez Oliveros et al. 2014). Considering an acceleration timescale of 10 minutes, the average acceleration rate is lower than the acceleration of gravity. A detailed analysis of rain blob acceleration process for non-flaring (or quiescent) coronal rain was given in Fang et al. (2013). Different from Fang et al. (2013), the guide field B z of the loops where condensation happens is nearly zero. But since coronal rain is formed through thermal instability which depends only slightly on the magnetic field, whether or not there is a guide field should have little influence on our results.

MASS AND ENERGY CYCLES DURING THE GRADUAL PHASE
Here we investigate the mass and energy cycles in our entire 100-minute simulation of the gradual phase and the role of condensations in it. To do so, we track mass and energy budgets into the coronal part (y > 5 Mm) of a loop section in which the first round of condensation happens. This loop section is always bounded by (a) the evolving magnetic field line with a fixed footpoint at x = −25 Mm at our lower y = 0 boundary; and by (b) a similarly evolving field line with footpoint at x = −15 Mm (Fig. 3a). As seen in Fig. 2 and Fig. 3a, this region gets emptied during the first round of condensations, as matter collects along the field lines into localized rain blobs. This entire loop system continuously moves downward as a result of the above reconnection dynamics and the corresponding area of the selected region continuously decreases during the simulation as quantified by the solid line in Fig. 3b. The area-integrated total energy (kinetic, thermal and magnetic combined, dashed line in Fig. 3b), area-integrated mass (solid line in Fig. 3c) and the area-averaged temperature (solid line in Fig. 3d) of this region also continuously decrease before the first condensation (t 35 min). Plasma which evaporated upwards into the corona in the previous impulsive phase now leaks back to the chromosphere in this period, seen in the downward mass flux (dashed line) in Fig. 3c. The decrease of temperature leads to a decrease of the atmospheric scale height. However, this downwards leakage of coronal plasma is severely reduced when the condensation happens after t 35 min, since the gas pressure in the coronal part of the loop then decreases rapidly. Therefore, we see a drop of the downward mass flux in Fig. 3c near t ≈ 38 min, when condensations fully formed. Later on, the downward mass flux experiences a sudden increase, exactly when the cool coronal rain material goes through the lower y = 5 Mm boundary of the studied region. The rain material is dense and heavy, so the material can still drop even if there is an upward pressure gradient below it. The area-integrated mass reaches its minimum value when all cool material leaves and enters the chromosphere. Thereafter, plasma from the chromosphere is injected to the loop again, due to the low pressure inside the loop, leading to an upwards mass flux. The total mass then gradually returns to its value before condensation. The changing coronal energy budget shows a similar tendency with the changing coronal mass cycle. The energy also experiences a decrease in the pre-condensation and during the condensation phase, to then experience an increase after the condensations merged into the chromosphere due to renewed plasma injection (Fig. 3b). It has been suggested that thermal conduction determines the PFLs energy loss at the beginning of the gradual phase of a flare and that subsequently radiative losses will become dominant (Cargill & Klimchuk 2004). Our simulation shows that this suggestion is correct before and also during the occurring condensations. The efficiencies of radiative cooling and of thermal conduction to the energy loss in the selected region are compared in Fig. 3d. The contribution of thermal conduction is greater than radiative losses for t 25 min, but conductive losses drop gradually owing to the decreasing temperature gradient. At t ≈ 25 min, the average temperature is about 3 MK, and the efficiency of radiative losses becomes most prominent. However, conductive losses become stronger than radiative losses again when the condensations vanished from the loop system, as collisions between re-injected flows from both footpoints make the loop hot again and the radiative losses drop for a while due to the decrease of loop density.

CATASTROPHIC COOLING AND RAIN-INDUCED QPP
The first round of condensation happens near t ≈ 35 min. The temporal evolutions of instantaneous maximum/minimum temperature/number density in the condensation region (the same as marked in Fig. 3a) are illustrated in Fig. 4. Triggering of thermal instability switches the radiative cooling process from linear to nonlinear, and then leads to a catastrophic cooling of local plasma. We get an average temperature decreasing rate of -9000 K s −1 in the catastrophic cooling phase. In contrast, the cooling rate before catastrophic cooling is -3000 K s −1 . As a result of catastrophic cooling, the local temperature decreases from 0.2 MK to 0.02 MK within half a minute (Fig. 4a), while the local number density increases by one order, from 10 10 cm −3 to 10 11 cm −3 (Fig. 4b). Noting that the cooling rate is density dependent, the cooling rate we get in the catastrophic cooling phase is in the same order as that observed in Scullion et al. (2016) (−22700 K s −1 ). A quasi-periodic pulsation (QPP) with a period of ∼ 3 minutes appears in the maximum density curve, just after the rain condensation disappeared from the PFL system. This QPP is caused by the injected flows mentioned in the previous section, refilling and reheating the PFL. The density variation due to these flows along a field line is shown in Fig. 4c. Injected flows propagating from one footpoint to the other produce reflected slow mode waves. Such a process has previously been studied in Fang et al. (2015b) for isolated loop systems. The sharp density changes in Fig. 4c are shocks ahead of the injection flows and the wave fronts of the slow mode waves. Such QPPs hence reflect density variations in the low corona due to combined flows and wave propagation. The period of our QPP is close to the time for the slow mode wave to propagate from one footpoint to the other, as the wave speed is about 300 km s −1 .

CORONAL RAIN AND DARK POST-FLARE LOOPS
In synthesized EUV images of our simulation, coronal loop(s) appear in the gradual phase. An example is shown in Fig. 5 at the 17.1 nm waveband. Interestingly, this loop disappears for about 10 minutes during the evolution, as demonstrated in panels (a) to (e) of Fig. 5. The sudden darkening of the bright EUV loop in our simulation resembles the dark post-flare loop (DPFL) phenomenon, previously observed at the same 17.1 nm passband and with similar timescale of 10 minutes (Song et al. 2016). Observed DPFLs and the disappearing EUV coronal loop in our simulation also share the same time evolution of integral EUV flux: the EUV flux reaches its minimum value when the darkening happens (compare Fig. 3 in ref. Song et al. (2016) and our Fig. 5f). The formation of a darkened coronal EUV loop needs to satisfy one or both of the following conditions: (1) an emission drop in an existing bright EUV loop; (2) an absorption of the background EUV emission (Anzer & Heinzel 2005). Here we explain how these conditions can be satisfied based on our simulation results. The role of cool and dense coronal plasma in EUV emission and absorption leading to DPFLs has been emphasized previously (Jejčič et al. 2018;Heinzel et al. 2020), with coronal rain observations (Scullion et al. 2014;Martínez Oliveros et al. 2014;Jing et al. 2016;Scullion et al. 2016) and our simulation results showing how this cool and dense plasma can be generated in PFLs. The drop in the loop emission is understood from our simulation: loop temperature changes relate to nearby rain condensations. Indeed, the temperature of bright loops in this 17.1 passband (about 10 5.8 K) is not far from the critical temperature for the onset of catastrophic cooling (this is density and temperature dependent, but generally happens below 2 MK according to Cargill & Bradshaw (2013)). Condensations can be triggered by thermal instability near bright loops and these suddenly formed structures grow fastest across magnetic field lines (counterintuitive due to the field-aligned thermal conduction, but see Fang et al. (2013); Claes & Keppens (2019); Claes et al. (2020); Zhou et al. (2020)). Rain that forms near ( Fig. 5a-b), and ultimately inside (Fig. 5c) the bright coronal loops, is thus causing the darkening as illustrated in Fig. 5c-d. Once condensation happens, a lot of plasma will collect into a small region, so the plasma density elsewhere in the loop decreases. To maintain pressure balance, these evacuated loop regions will increase in temperature, so EUV brightness decreases due to these combined temperature and density changes. Ultimately, the loop refills and brightens once more (Fig. 5e).

CONCLUSION AND DISCUSSION
To fully understand coronal rain in PFLs and the mass and energy budget in the gradual phase of solar flares, we performed a flare simulation from onset all the way into the long duration post-flare phase. Post-flare coronal rain successfully and, at least in our simulation, repeatedly forms. The flare-induced rain is a result of catastrophic cooling by thermal instabilities, and our simulation shows successive rain formation at increasing heights in the PFL configuration. Falling rain blobs into the chromosphere lead to sudden mass drops in PFLs, but their mass increases again due to spontaneously forming injection flows. Therefore, the coronal rain events do not accelerate the PFL mass loss in the longer term. Such longer term mass loss is more determined by the change of the gravity scale height due to the cooling of PFLs.
Both thermal conduction and radiative losses contribute to the energy budget in PFLs. Thermal conduction dominates the PFL energy loss at the beginning of the gradual phase. Thereafter, it becomes less efficient than radiative losses, owing to decreases in loop temperature and in temperature gradient. However, thermal conduction can efficiently recover again after a condensation falls to the chromosphere, as the loop reaches again a high temperature. In this phase, an emptied loop refills and can show a slow-wave related QPP.
We showed that the formation of DPFLs can result from post-flare rain condensations. Condensations change loop temperatures and can make existing bright EUV loops temporarily disappear for several minutes. This timescale of EUV loop darkening is identical to observed DPFLs.
A recent paper by Reep et al. (2020) suggests that additional energy supplement (heating) is necessary in the formation of flare-driven coronal rain based on a 1D parameter study. In their simulations, a 1D flare loop is filled with hot/dense plasma by non-thermal electron-beam-driven chromospheric evaporations at the beginning and then gradually cools down due to radiative loss. The loop will directly cool down to chromospheric temperature without condensations happening if there is no extra heating provided during the cooling, in which the coronal density also drops by orders of magnitude during the cooling. It is interesting to note that condensations can still happen inside our 2D post-flare loop, even if the background heating is switched off (see the additional simulation mentioned in section 2. Compression caused by thermal pressure of the surrounding environment or by the Lorentz force may serve as additional energy supplement, as it has been shown in Fig. 3b that the loop area decreases gradually during the cooling phase. When this interaction with the surrounding environment is taken into account, the coronal temperature and density are not likely to drop by orders of magnitude at the same time as shown in Reep et al. (2020). This difference between 1D simulations and 2D simulations of post-flare loop development will be studied in future work.