Coronal Rain in Randomly Heated Arcades

Our simulations are performed by solving MHD equations that include gravity, field-aligned heat conduction, parameterized heating terms, and radiative cooling:

$$\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 (\rho {\boldsymbol{v}}{\boldsymbol{v}}+{p}_{\mathrm{tot}}{\boldsymbol{I}}-\displaystyle \frac{{\boldsymbol{B}}{\boldsymbol{B}}}{{\mu }_{0}})=\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

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

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

where v , ρ, I , and B are the velocity, plasma density, unit tensor, and magnetic field, respectively. Assuming that plasma is completely ionized and the abundance of hydrogen and helium is 10:1, then we can obtain density ρ = 1.4m_p n_H with the number density of hydrogen n_H and the proton mass m_p . The total pressure p_tot = p + B²/2μ₀ consists of thermal pressure p and magnetic pressure, where p is obtained by the ideal gas law p = 2.3n_H k_B T. The gravitational acceleration is usually calculated by ${\boldsymbol{g}}=-{g}_{\odot }{r}_{\odot }^{2}/{({r}_{\odot }+y)}^{2}\hat{y}$ where r_⊙ is solar radius, and g_⊙ = 274 m s⁻² is the solar-surface gravitational acceleration. The E = p/(γ − 1) + ρ v²/2 + B²/2μ₀ is total energy density with the ratio of specific heats γ = 5/3. The anisotropic thermal conduction along the magnetic field lines is quantified by the term that contains κ = κ_∥ΓΓ, where the Spitzer conductivity κ_∥ equals 10⁻⁶ T^5/2 erg cm⁻¹ s⁻¹ K⁻¹, and b = B /B is the unit vector along the magnetic field. The radiative cooling term $R=1.2{n}_{H}^{2}{\rm{\Lambda }}(T)$ is composed of the number density of hydrogen and the radiative loss function for optically thin emission Λ(T), which is interpolated from the data given in Colgan et al. (2008) as demonstrated in previous simulations (Xia et al. 2012, 2014; Keppens & Xia 2014; Fang et al. 2015). Λ(T) is set to be zero below 10,000 K, because the plasma there is no longer fully ionized and becomes optically thick. Many optically thin cooling curves have been implemented in MPI-AMRVAC, and Hermans & Keppens (2021) studied the influence of these cooling curves on condensation formation and found that, for different cooling curves, condensations may have different formation times and morphology. H is the heating term that includes the background heating and localized heating H = H_bgr + H_loc. The background heating can resemble the photospheric source for coronal heating and balance the energy losses to help the corona reach a steady state. Here we adopt a vertical exponentially decaying background heating rate as follows:

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

where c₀ = 10⁻⁴ erg cm⁻³ s⁻¹ and λ₀ = 50 Mm.

For the localized heating H_loc, it is added after the system reaches a quasi-equilibrium. Here we adopted a randomized heating term to imitate the energy input from the lower solar atmosphere that has turbulent convection and moving magnetic features, similar to recent previous work (Zhou et al. 2020). The heating function is randomly distributed both in space and time, which can be expressed as H_loc = c₁ H(x)H(y)H(t), where c₁ equals 10⁻² erg cm⁻³ s⁻¹. The randomized spatial distribution is imposed on the x direction only, while H(y) is set to be

$$\begin{eqnarray}&&H(y)=\left\{\begin{array}{l}\exp (-{\left(y-{y}_{c}\right)}^{2}/{\lambda }_{1})\quad y\gt {y}_{c}\\ 1\quad y\leqslant {y}_{c}\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where λ₁ = 2 Mm² and y_c = 3 Mm, representing the height of the TR in the final relaxed quasi-equilibrium system. The heating term in the x direction is a superposition of a series of waves of different wavelengths from the grid size to the simulation box: $H{(x)={{\rm{\Sigma }}}_{i}{A}_{i}\sin (2\pi x/{\lambda }_{i}+{\phi }_{i}))}^{2}$, where A_i is the amplitude, and ϕ_i is the phase angle; and the wavelength is defined by ${\lambda }_{i}=i({\lambda }_{\max }-{\lambda }_{\min })/{10}^{3}+{\lambda }_{\min }$ where ${\lambda }_{\min }$ is the grid size, ${\lambda }_{\max }$ is twice the width of the domain, and i = 1, 2, 3, ..., 10³. The heating in the solar atmosphere is considered to have a power-law spectrum; as a consequence, we adopt A_i ∝ λ^s , where the spectral index s is chosen to be 5/6 (Matsumoto & Kitai 2010). Zhou et al. (2020) also mentioned that the choice of the spectral index did not have too much impact on the results when 0.2 ≤ s ≤ 2. Observational works provide evidence that heating in the solar corona is highly episodic (Testa et al. 2014, 2020; Reale et al. 2019). Testa et al. (2020) observed that two different peaks in 94 and 131 Å channels appear with 5 minutes intervals, which mark the presence of two distinct heating episodes. Large-scale corona simulations also reveal that chromosphere and hot corona could be maintained self-consistently when the convection zone is included in the simulation domain, and heating shows timescales from 2 to 5 minutes (Hansteen et al. 2015). Here the temporal distribution of the heating profile is an episodic Gaussian profile that has 5 minutes ± 75 s durations, attributed to the typical timescale and correlation time of the lower solar atmospheric motions. The localized heating can then be expressed as the sum of a series of Gaussian functions with respect to time under these assumptions. To normalize the equations, we adopt normalization units of length L₀ = 10 Mm, temperature T₀ = 10⁶ K, and number density n₀ = 10⁹ cm⁻³. The other normalization units can be deduced accordingly.

Following Fang et al. (2013, 2015), we adopt a linear force-free magnetic configuration where all field lines make a constant angle of θ₀ = 30° with the neutral line (x = 0, y = 0) for the initial magnetic field:

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

Here B₀ = 20 G, and L = 60 Mm is the horizontal size of our domain from −30 to 30 Mm.

In the initial state, the thermal structure below the TR height of h_tr = 2 Mm is a vertical distribution of temperature of about T_bo = 8000 K at the bottom with smooth connection to T_top = 1.5 × 10⁶ K at the top using a hyperbolic tangent function:

$$\begin{eqnarray}&&\begin{array}{l}T(y)={T}_{\mathrm{bo}}+0.5({T}_{\mathrm{top}}-{T}_{\mathrm{bo}})\\ \quad \times (\tanh ((y-{h}_{\mathrm{tr}}-a)/{w}_{\mathrm{tr}})+1)\,y\leqslant {h}_{\mathrm{tr}},\end{array}\end{eqnarray} \tag{ 8 }$$

where w_tr = 0.2 Mm and a = 0.027. The temperature profile above the TR height is given by

$$\begin{eqnarray}&&T(y)={[3.5{F}_{c}(y-{h}_{\mathrm{tr}})/{\kappa }_{\parallel }+{T}_{\mathrm{tr}}^{7/2}]}^{2/7}\,y\gt {h}_{\mathrm{tr}},\end{eqnarray} \tag{ 9 }$$

where the vertical thermal conduction flux is F_c = 2 × 10⁵ erg cm⁻² s⁻¹, and T_tr = 1.6 × 10⁵ K (Zhou et al. 2018). The number density at the bottom is set to be 1.151 × 10¹⁵ cm⁻³, and then the initial density is derived under the assumption of hydrostatic equilibrium. The initial velocity is set to be zero.

For boundary conditions, two grid layers exterior to the domain are used as ghost cells. We use a fixed zero velocity, fixed magnetic field, and fixed gravity stratification of density and pressure predetermined in the initial condition for the bottom boundary. At the left and right physical boundaries, density, energy, y component of momentum, and B_y are set symmetrically; meanwhile x and z components of momentum, B_x , and B_z are adopted asymmetrically to ensure zero face values. As for the top conditions, we use a fixed zero velocity and employ a separate pressure-density extrapolation from the gravity stratification as introduced in the initial condition. We determine B_x and B_z using a zero-gradient magnetic field extrapolation and derive B_y with a second-order centered difference evaluation of zero-divergence constraint.

We use MPI-AMRVAC ¹ to solve the governing equations adopting a three-step Runge–Kutta time integration with a third-order slope-limited reconstruction (Čada & Torrilhon 2009) and Harten-Lax-van Leer Riemann solver. To fulfill that the magnetic field divergence is close to zero, we employ the upwind constrained transport (CT) method that uses staggered grids for the magnetic field, which was applied in the companion general relativistic MHD code variant BHAC (Gardiner & Stone 2005; Porth et al. 2017). Radiative cooling and field-aligned thermal conduction are handled with the exact integration method (Townsend 2009) and the Runge–Kutta Legendre super-time-stepping scheme (Meyer et al. 2012, 2014), respectively. As mentioned in Section 1, we adopt the TRAC method to correct the evaporation. The basic idea is defining a cutoff temperature T_c , and then modifying the local evaluations of radiative cooling and thermal conduction terms at the places that have a temperature lower than T_c (Lionello et al. 2009; Mikić et al. 2013; Johnston et al. 2017a, 2017b, 2020; Johnston & Bradshaw 2019). In this paper, we employ the basic 1D TRAC method, and the local cutoff temperature inside of each block is calculated separately when it is used for 2D or 3D simulations. It is worth mentioning that a field line-based TRACL method and a block-based TRACB method (Zhou et al. 2021), as well as the LTRAC method proposed by Iijima & Imada (2021), can now also be used for multidimensional simulations inside MPI-AMRVAC.

To help the system with the above initial setup reach a quasi-equilibrium state, we run the simulation with H = H_bgr for 171.6 minutes and make sure that the maximal residual velocity is below 5 km s⁻¹. Beginning with this equilibrated system, we reset the time to zero and include the localized heating H =H_bgr + H_loc to the simulation. The localized heating evaporates plasma and causes upflows into the corona, and the densities around the loop apexes continue to increase. Meanwhile, the temperature first increases but then starts to reduce slowly (see Figure 1). This reduction happens as radiative cooling starts to dominate over heating and thermal conduction, which leads to the onset of thermal instability. After about 53 minutes of localized heating, near the loop top region, a small condensation segment with temperature of around 0.02 MK and number density of 6.8 × 10¹⁰ cm⁻³ suddenly forms, which is shown in the left panels in Figure 1. Figure 2 exhibits the details of a loop top region denoted by the green rectangle in Figure 1(a1). In this loop top region, temperature and gas pressure decline significantly, causing matter to be drawn in from the surroundings. The induced siphon flows, coming from two sides with speeds up to 57 km s⁻¹, form the condensation segment, which has a length of 1.5 Mm (see Figure 2). Similar condensation processes continuously take place and extend into strands on both sides of the first affected coronal loop, which was ascribed to a common footpoint heating process (Antolin & Rouppe van der Voort 2012; Fang et al. 2013, 2015). The large-scale condensations therefore exhibit a zigzag-shape as shown in panels (a2) and (b2) of Figure 1. As the condensations continue to develop, they lose their balance spontaneously and start to descend slowly along either side of the magnetic field lines and split into several smaller blobs. After this large-scale condensation process, condensations continuously occur at nearly every place on the coronal loops, which seems more disorganized and less vigorous (see panel (a4)).

Zoom In Zoom Out Reset image size

During the first condensation, we observe similar phenomena such as siphon flows and rebound shocks compared to earlier simulations in Xia et al. (2011, 2012, 2017) and Fang et al. (2013, 2015), as presented by the temporal evolution of number density, gas pressure, velocity magnitude, and temperature maps in Figure 3. Panels (a1) and (b1) show the dense plasma blob formed by the siphon flows, and its pressure is much higher than the surroundings. The pressure at the right side of the blob is much lower than the left side. The velocity magnitude map in panel (c1) reveals that the siphon flow at the left of this condensation has a higher speed of 67 km s⁻¹, compared to the right siphon flow at a speed of around 10 km s⁻¹. From 52.955 minutes to 54.386 minutes, the location of this condensation moved to the right for about 1.0 Mm in the presence of pressure gradient and stronger siphon flow at its left, as shown by the blue arrows in panel (a2). The solid arrow indicates the present position of the condensation while the dashed one denotes the position at the previously shown snapshot. Meanwhile, one left rebound shock was generated and propagated against the continuous siphon flows along the magnetic field. The shock front displays a fan-shaped structure, as the processes of condensation and shock formation on adjacent magnetic field lines start at different times, which first happen at the center and then propagate outward from the blob. At this instant, the shock compresses the plasma from 1.8 × 10⁹ to 3 × 10⁹ cm⁻³, decreases the inflow speed from about 67 km s⁻¹ to around 20 km s⁻¹, and heats the plasma from 0.1 to 0.6 MK (see panels (a2)−(d2)). About 1.43 minutes later, the right-directed rebound shock is also generated and begins to propagate and heat the plasma (see panels (a3)−(d3)). As the left rebound shock occurs earlier, its front moves farther than the right one. The plasma near both fronts becomes hotter, as denoted by the black arrows in panel (d3).

Zoom In Zoom Out Reset image size

As the blobs lose their balance and move downwards, they eventually hit the TR and enter the lower chromosphere, and we find that nearly all of them trigger concurrent upflows that follow the same loops as the downflows. Panels (a1)−(a4) in Figure 4 show the number density map between 120.222 minutes and 124.516 minutes when a falling blob enters the TR, while panels (b1)−(b4) and (c1)−(c4) display the corresponding pressure and temperature change. The blue arrows point out the positions of the blob we focus on. As the blob falls downwards, it compresses the plasma on its way, which makes the pressure in front of it increase, and the pressure in its wake is greatly reduced. After the blob reaches the TR, the local density and pressure increase, which leads to an upward force, and plasma is sent back to the corona. These upflows move along the same coronal loops as the coronal rain blobs, and the plasma in the loops is heated to a higher temperature of 1.51 MK, as shown in panel (c4). The similar phenomenon has been observed in observations and simulations (Tripathi et al. 2009; Antolin et al. 2010; Kleint et al. 2014; Fang et al. 2015). Fang et al. (2015) found that before the rebound shock and upflows reach far into the loop, the temperature of the low-density loop increases due to the background heating. Then the rebound shocks and upflows, produced by the impact of the rain, heat the low-density loops efficiently to a high temperature, and the heating follows the movement of the upflows. Further upflows coming from evaporation due to the footpoint heating will also heat the loop afterward, which can be seen from the animation associated with Figure 1.

Zoom In Zoom Out Reset image size

After the first condensation process, the whole system gets into a constant cycle where both small and large blobs coexist, as displayed in Figure 1(a4). The green and blue arrows point to a small condensation and a large condensation, respectively. Large clumps consisting of numerous condensations are usually denoted as coronal rain showers (Antolin & Rouppe van der Voort 2012). To compare the results with the observed images from the Atmospheric Imaging Assembly (AIA) instrument on board the Solar Dynamics Observatory, Figure 5 displays the number density and synthetic EUV images of the coronal loops and coronal rain at 228.995 minutes in our simulation. The technique to obtain synthetic images like this is described in Xia et al. (2014). Here we choose the 171, 193, and 304 Å wave channels, which have main contribution temperatures of around 0.8, 1.5, and 0.08 MK, respectively. In the number density map (see panel (a)), we can clearly see that condensed plasma is connected to others and combined into large clumps. The hot 171 and 193 Å channel images in panels (b) and (c) distinctly show the bright coronal loops and the dark coronal rain clumps embedded in them. In the 304 Å channel, the coronal rain clumps look bright at the boundary and dark inside, implying that these clumps have a core with a temperature around or even lower than 10⁴ K, which is consistent with the observations of coronal rain showers (Antolin et al. 2015).

Zoom In Zoom Out Reset image size

To explore the dynamics of the coronal rain blobs, we count all the blobs that appeared during the simulation and study their properties statistically. We use the former standard (Hirayama 1985; Müller et al. 2005; Antolin et al. 2010) to define a grid cell located above the chromosphere-corona-transition-region as a component of the coronal rain blobs if its number density is larger than 7 × 10⁹ cm⁻³ and has a temperature beneath 1 × 10⁵ K. Regarding all connected cells as one coronal rain blob, we identify the coronal rain blobs and collect their characteristics, such as number, mass, width, length, area, centroid location, and centroid velocity at each saved snapshot. Using this identification method, we find that the total number of blobs is 6060 in all 427 snapshots (the number of snapshots for which we chose to save all data at a cadence of ${\rm{\Delta }}t=1.431\,\min$ from the condensation start until the end of simulation), which does of course count blobs that moved location multiple times. There is a maximum of 24 blobs in one snapshot. The temporal change of the number of coronal rain blobs at any instant is shown in Figure 6(a). Panel (b) displays the variation of the mass of hot, warm, and cool plasma in the corona over time. The "cool" plasma refers to the plasma with a temperature below 0.02 MK, and the "warm" plasma is the plasma between 0.02 MK and 0.5 MK, and the plasma hotter than 0.5 MK is identified as "hot" plasma. As soon as the simulation starts, the mass of hot plasma begins to increase continuously, and it increases by about 4.5 × 10⁴ g cm⁻¹, almost tripled within 53 minutes. At around 53 minutes, as the mass of hot plasma decreases, the mass of warm plasma begins to increase, and then the mass of cool plasma also increases, which illustrates that catastrophic cooling processes cool down hot plasma in situ in the corona, and the first condensation formed. From t ∼ 75 minutes, blobs begin to enter the TR, the mass of cool plasma reaches its first peak and then starts to decrease. After blobs hit the chromosphere, there are upflows entering the relatively empty loops, heating the loop as presented in Figure 4, so the mass of hot plasma has a tendency to increase. Then the number and mass curves both seem to have several cycles, which has been investigated in previous observational works (Antolin & Rouppe van der Voort 2012; Kohutova & Verwichte 2016; Froment et al. 2020) and simulations (Müller et al. 2003; Fang et al. 2015). Here we do not see any time gap between these cycles, i.e., the coronal rain blobs continue to occur all the time, different from the results from Fang et al. (2015). This is likely a direct result of the randomized heating pattern we adopted; as a result, the chromospheric evaporation rate is not constant and has a different spatial distribution. This difference indeed confirms earlier studies on the basis of observations, where one speculated that the temporal behavior of coronal rain cycles may well reflect the heating prescription at work, and thus may provide clues to the coronal heating paradigm (Antolin et al. 2010). Although the continuous but randomized heating imposed on our arcade footpoints has a spectral behavior with a peak of about 5 minutes, the time sequence displayed in Figure 6 rather shows peaks at much longer timescales: Fast Fourier Transform (FFT) results of the cool plasma mass signal (blue curve) and the hot plasma mass signal (red curve) show main periodicities of 83 minutes (see panels (c) and (d)). Longer ones such as the peak at 332 minutes may be meaningless, as the FFT is not accurate enough due to the insufficient simulation time.

Zoom In Zoom Out Reset image size

We track several different field lines during the whole simulation, and the temperature and number density evolution along the magnetic field line passing the footpoint at x = 16 Mm are demonstrated in Figures 7(a) and 7(b), respectively. The blue curves in panel (b) represent the change of localized heating at two footpoints, with the lower curve exhibiting the original value at the left footpoint, and the upper curve is drawn by 0.3 minus the value at the right footpoint. We can see that, as soon as the localized heating starts, the loop is heated to a temperature of several MK. Around tens of minutes later, the loop begins to cool down, and the temperature drops to the TR and chromospheric temperature in approximately 20 minutes; meanwhile the loop density gradually increases. Suddenly the condensations form, move inside the loop, and eventually enter the lower atmosphere. Then the loop continuously undergoes heating and cooling, and condensations occur from time to time. Panel (c) displays the time evolution of the mass of coronal rain blobs on this magnetic field line, from which we can obtain that the coronal rain in a single field line experiences several cycles. The FFT result in panel (d) illustrates that the coronal rain mass on this field line has several period components, similar to the period of the whole system. Comparing the results of different field lines, we find that different field lines have similar periods, which implies that a synchronizing mechanism is acting across the field lines, making coronal rain occur near-simultaneously across neighboring field lines rather than purely driven by the field line's own heating, as suggested by Antolin & Rouppe van der Voort (2012) and Antolin (2020). This synchronizing mechanism could be purely due to the magnetic and gas pressure variations across the field lines, as shown by Fang et al. (2013, 2015). A single loop or adjacent magnetic loops experience typical heating-condensation cycles repeatedly: heating leads to chromospheric evaporation and then radiative cooling accompanied by thermal instability and condensation of the plasma. For the whole system, it undergoes a continuous heating-condensation cycle, including heating, radiative cooling resulting in the subsequent plasma condensations, and long competing phases during which the heating and cooling keep happening, and the condensations occur all the time. Changes in plasma properties are communicated across field lines by fast magnetosonic waves.

Zoom In Zoom Out Reset image size

Recent observational studies reveal that EUV intensity pulsations with periods of several hours are a common phenomenon in solar coronal loops (Auchère et al. 2014, 2016; Froment et al. 2015), and they are interpreted as signatures of evaporation—condensation cycles (Froment et al. 2017, 2018). In order to compare our results with observations, we also calculate the intensity of the synthetic EUV images. The total flux of the AIA 94, 171, 193, and 304 Å emissions in the corona as time curves are presented in Figure 8 panels (a)–(d). The AIA 94 Å channel has strong responses to temperatures at about 6.8 MK, so it reflects the change of the mass of hot plasma. As mentioned before, the 193 and 171 Å channels also have the main contribution from hot and warm plasma. As a result, the AIA 94 and 193 Å flux increase immediately after the localized heating starts and have similar morphological changes with the mass curve of hot plasma in Figure 6 (b). Compared to these hot channels, the start time of the 171 Å emission is delayed. The intensity of AIA 304 Å reflects the emissions of cold coronal rain, which increases drastically after tens of minutes. All these channels display cycles with similar periods compared to the mass of coronal rain (see panels (e)–(h)).

Zoom In Zoom Out Reset image size

Figure 9 displays the distributions of width, length, area, and mean number density of our simulated coronal rain blobs. The width of each blob is calculated from the diameter of its biggest inscribed circle. Most blobs have widths less than 800 km, as shown in panel (a). The length of each blob is calculated by half the perimeter of the blob minus its width, to mimic the curved length along the clump's spine. The lengths of blobs distribute broadly from a few hundred km to 25 Mm, with the bulk peaking at 1.5–2.0 Mm. As there are large condensations and irregularly shaped blobs, the area may describe the scale of the blob more accurately, as shown in panel (c). The areas of most blobs are under 5 Mm², and blobs are dominated by smaller ones with areas under 0.5 Mm². The mean number density of the blobs is ranging from 7 × 10⁹–4.4 × 10¹⁰ cm⁻³ (see panel (d)). The results are in good agreement with previous studies (Fang et al. 2013; Antolin et al. 2015; Antolin 2020). Due to the limit of the simulation resolution, small-scale condensations may be not observable, which is considered as the tip-of-the-iceberg effect and therefore suggests a number increase for blobs with smaller sizes for higher resolution (Antolin & Rouppe van der Voort 2012; Fang et al. 2013; Scullion et al. 2014). We may study even more highly resolved rain dynamics in follow-up work, with the distinct possibility to identify small-scale, secondary fluid instabilities at work.

Zoom In Zoom Out Reset image size

Dynamic properties of coronal rain plasma are analyzed and shown in Figure 10. The distributions of velocity and acceleration are displayed in panels (a) and (b). The velocities of the coronal rain blobs have a wide distribution varying from several km s⁻¹ to above 120 km s⁻¹, with a mean speed of 22.4 km s⁻¹. Most blobs have a speed under 40 km s⁻¹. The acceleration of the majority of the blobs is between −200 and 200 m s⁻², and is clustered around 50 m s⁻². The measured acceleration is typically below the gravitational-freefall acceleration, which has been pointed out and explained by the pressure-mediated effects in Müller et al. (2003), Antolin et al. (2010), and Fang et al. (2015). Panel (c) shows the percentage of downward blobs (V_y < 0) as a function of time. We can see that, although most blobs are falling downward, there are blobs moving upwards at basically every moment. Sometimes the number of rising blobs is even higher. The scatter plot in panel (d) displays the y-velocity components and the mean number densities of the blobs. We can see that, as the number density increases, fewer blobs are moving upward, and their upward speeds also decrease. For the denser blobs, they hardly move upwards. This can be explained since a blob of larger density requires a larger pressure gradient to overcome the acceleration of gravity and move upwards, which has also been mentioned in previous studies (Oliver et al. 2014; Martínez-Gómez et al. 2020).

Zoom In Zoom Out Reset image size

As we have the properties such as positions and velocities of all blobs in the simulation, we can track the movement of individual blobs. Here we use the position and velocity of a blob to predict the position where it may appear at the next moment, and look for the blob appearing near this position. The spatial distribution of all blobs in the simulation is displayed in Figure 11. We can see that the blobs appear nearly everywhere in the magnetic loops. In order to make the figure more concise, only the trajectories of the blobs between 133.103 minutes and 194.646 minutes are indicated by the red curves. Most blobs fall down along the magnetic loops, while there are still many blobs moving upwards as quantified in Figure 10(c), and some of them even pass the top of the loop to the other side. Figure 12(a) is the number density map showing the movement of a blob in detail, denoted by the black arrow in Figure 11 and blue arrows in Figure 12. This blob is formed at 133.103 minutes and located at (−22, 16) Mm. Panel (b) displays the time evolution of the number density along the field line where the blob is located, which clearly shows the movement of this blob. The change of the x-velocity component and the y-velocity component of this blob from 133.103 minutes (panel (a1)) to 186.058 minutes (panel (a6)) are denoted by the blue and red lines, respectively. Panel (c) displays the time evolution of the pressure along this field line, and the blue curve shows the trajectory of this blob. The trajectory and the velocities are obtained from the tracking method mentioned before, which is consistent with the slice plot along the field line. The blob first moves upwards and then downwards, and from the pressure change in panel (c), we can see that the change of the movement direction of this blob is closely related to the change of the pressure gradient on the upper and lower sides of the blob. At first, the pressure below the blob is higher than the pressure above it, and the blob moves upwards lifted by the pressure gradient; then the blob moves downwards as the pressure ahead becomes bigger, but it immediately changes direction as the pressure gradient changes direction again. The change of the pressure gradient may result from a piston effect due to the compression of the moving blob and the upflows caused by the heating from the footpoints (Adrover-González et al. 2021). After 174.609 minutes (indicated by the green dashed line in panel (b)), the blob starts to accelerate downwards and catches up with the condensation next to it, after which they merge and finally fall down to the chromosphere (see panels (a7) and (a8)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 13 shows substructures of this coronal rain blob at 141.691 minutes, and its field of view (FOV) is outlined by the green rectangle in Figure 12(a2). As pointed by the black dashed vertical line in Figure 12(b), the blob is moving downwards at this moment, while at the next snapshot it is moving upwards. To better quantify this, we make two lines crossing the center of the blob shown in panels (a) and (b), and the temperature, gas pressure, and radiative loss along these lines are displayed in panels (c) and (d). The solid line is parallel to the field line, while the dashed line is perpendicular to it. As mentioned in the last paragraph, along the field line, the pressure below the blob is much higher than the pressure above it (see the blue curve in panel (c)), as the descending blob and the evaporation at the loop footpoint may compress the plasma between them. In agreement with, e.g., Oliver et al. (2014), the strong pressure gradient supports the blob and slows down the blob in its descent, and even lifts the blob to move upwards, which explains why the blob changes the direction of movement at the next snapshot. The coronal rain blob is surrounded by an area where the temperature declines from more than 0.5 MK outside the blob to 0.02 MK inside the blob. This area is usually considered as PCTR structure in prominences and has been found in coronal rain as well (Antolin et al. 2015; Fang et al. 2015; Antolin 2020). The radiative loss curves exhibit two peaks near this area, indicating that catastrophic cooling mainly occurs at the boundary of blobs, ensuring the condensation continues to grow.

Zoom In Zoom Out Reset image size

Except for simple O-shaped blobs, there are many condensations that show V-like structures during propagation. Figures 14(a1)–(a6) is the number density map showing the movements of blobs during 94.460 minutes to 125.947 minutes. The blue arrows in panels (a1) and (a2) indicate two blobs we focus on, which are formed at the same time and located at nearly the same magnetic loops. We choose a magnetic field line passing these two blobs, as denoted by the dashed curve in panel (a5), and the number density and pressure changes along this field line are presented in panels (b) and (c), respectively. These blobs move upwards and approach each other, as the gas pressure between them is much lower than the pressure outside (see panel (c)), which is caused by catastrophic cooling and the lower density in between them (see panel (b)). During this period, the condensations extend from the position where the first condensation occurs, forming larger blobs shown in panel (a2). As the condensations at different loops are experiencing different pressure gradients and gravity, the velocity of the condensation segments are different, and the blobs change their shapes accordingly. Both blobs have faster siphon flows at the position closer to their center, consequently they become deformed and exhibit V-shaped patterns shown in panel (a3). At 114.498 minutes, these two blobs meet and merge, then they move toward the right footpoints of the coronal loops together and enter the chromosphere ultimately. After they merge together, the density of the new blob becomes bigger (see panel (b)), and it moves downwards against the pressure gradient, as the increased pressure underneath it cannot balance out the gravity anymore. The radiative cooling and pressure maps at the green square in panel (a3) are displayed in Figures 15(a) and 15(b). The temperature, gas pressure, and radiative loss along the solid and dashed lines are plotted in Figures 15(c) and 15(d), respectively. We can see that the V-shaped blobs also have the PCTR area around them. In this area, the temperature changes rapidly from a coronal temperature to a chromospheric temperature, the radiative cooling is stronger, and the pressure is lower. During the propagation of these blobs, catastrophic cooling further happens in their heads and tails, resulting in the elongation of the blobs. We can see a higher radiative loss at the heads or tails of both O-shaped and V-shaped blobs in Figures 13 and 15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size