Annihilation of Magnetic Islands at the Top of Solar Flare Loops

Magnetic reconnection, a fundamental energy-releasing process in magnetized plasma, is believed to be the core mechanism driving solar eruptions, including solar flares, coronal mass ejections (CMEs), and other solar activities. According to the standard flare model, the sheared magnetic field is stretched and forms a large-scale thin current sheet (CS) in the wake of the eruption (Kopp & Pneuman 1976; Lin & Forbes 2000; Lin et al. 2015). The reconnection in the CS in turn drives a CME eruption and the formation of flare loops. Theoretically, a squeezed thin CS becomes unstable as the tearing mode instability grows (Priest & Forbes 2000; Pucci & Velli 2013; Loureiro & Uzdensky 2015). In particular, in a high-Lundquist-number corona, the plasmoid instability will dominate the dynamics of the magnetic reconnection (Bhattacharjee et al. 2009; Samtaney et al. 2009; Huang & Bhattacharjee 2010; Huang et al. 2011; Shen et al. 2011; Huang & Bhattacharjee2013, 2012; Loureiro & Uzdensky 2015; Huang et al. 2017; Zhao & Keppens 2020). The classical CS predicted by Sweet–Parker is thus split and the reconnection enters into a fast scheme and generates a number of magnetic islands.

Magnetic islands are believed to be closely related to fine structures in the CS (Shen et al. 2011; Lin et al. 2015). The coalescence of magnetic islands can further produce secondary CSs and islands (Bárta et al. 2011). Various researchers have analyzed the cascading law of magnetic island flux theoretically or numerically, which predicts that the distribution function of the magnetic island flux $f\left(\psi \right)$ in the CS follows a power law between ψ⁻¹ to ψ⁻² (e.g., Uzdensky et al. 2010; Huang & Bhattacharjee 2012; Shen et al. 2013; Lynch et al. 2016; Ye et al. 2019). The fragmented and turbulent characteristics of the CS have also been suggested by recent observations (e.g., Cheng et al. 2018).

The flare loop top, interacting with the reconnection downflows, also shows complex dynamics. The turbulent characteristics of the loop top have been studied by high-resolution simulations (e.g., Cai et al. 2019; Ye et al. 2020; Cai et al. 2021). The generation of flare quasi-periodic pulsations (QPPs) and the supra-arcade downflows above the flare loop top is interpreted to be closely related to the fast reconnection downflows (Takasao & Shibata 2016; Guo et al. 2014). Abundant MHD shocks have been predicted by simulations (Takasao et al. 2015). In particular, as an important prediction of the standard solar flare model, fast mode termination shocks (TSs), formed when the downflows hit the relatively steady high-density structure at the loop top, are believed to be an effective mechanism to accelerate particles (Tsuneta & Naito 1998; Shen et al. 2018; Kong et al. 2019, 2020; Ruan et al. 2020), which are supported by radio imaging observations (Aurass & Mann 2004; Chen et al. 2015). The downflow magnetic islands also affect the dynamics of the flare loops (Cai et al. 2019; Ye et al. 2020), the formation of TSs (Shen et al. 2018), and even the acceleration of energetic electrons (Kong et al. 2020).

Nevertheless, mainly owing to the limitation of observations and simulations, magnetic reconnection in the CS and its relation to the dynamics at the flare loop top are still far from being fully understood. In this paper, we perform high-resolution 2.5-dimensional MHD simulations of the CS reconnection in a high-Lundquist-number and low-β coronal environment and focus on the dynamics of magnetic islands evolving toward flare loops. We find that the downflow magnetic islands quickly annihilate after colliding with the tip of cusp-shaped flare loops. Such a second reconnection process is characterized by horizontal CSs, smaller-scale islands, and a comparable reconnection rate with that in the main CS. In Section 2, we describe our method. Section 3 presents the main results, which are followed by a summary and discussion.

In this work, the MHD equation we solve is as follows:

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}\displaystyle \frac{{\rm{\partial }}\rho }{{\rm{\partial }}t}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{u}}\right) & = & 0,\\ \displaystyle \frac{{\rm{\partial }}\left(\rho {\boldsymbol{u}}\right)}{{\rm{\partial }}t}+{\rm{\nabla }}\cdot \left(\rho {\boldsymbol{u}}{\boldsymbol{u}}-{\boldsymbol{BB}}+{P}^{\ast }\right) & = & 0,\\ \displaystyle \frac{{\rm{\partial }}e}{{\rm{\partial }}t}+{\rm{\nabla }}\cdot \left[\left(e+{P}^{\ast }\right){\boldsymbol{u}}-{\boldsymbol{B}}\left({\boldsymbol{B}}\cdot {\boldsymbol{u}}\right)\right] & = & {\rm{\nabla }}\cdot \left({\kappa }_{\parallel }\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}T\right),\\ \displaystyle \frac{{\rm{\partial }}{\boldsymbol{B}}}{{\rm{\partial }}t}-{\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}\right) & = & -{\rm{\nabla }}\times \left(\eta {\boldsymbol{J}}\right),\\ {\boldsymbol{J}} & = & {\rm{\nabla }}\times {\boldsymbol{B}},\end{array}\end{array}\end{eqnarray} \tag{ 1 }$$

where P ^* = (p + B²/2) I , $e=p/\left(\gamma -1\right)+\rho {u}^{2}/2+{B}^{2}/2$, κ_∥ = κ₀ T^2.5, and standard notations of variables are used. In our model, only the thermal conduction (TC) along the magnetic field is considered, which is much more important than the perpendicular components, and κ₀ is set as 10⁶ erg · s⁻¹ cm⁻¹ K^−3.5 (Yokoyama & Shibata 2001). Here, all variables have been normalized according to constant units. The unit of space is chosen as the scale of simulation region L₀ = 5 × 10⁹ cm, which is comparable with the value for a typical observed flare event (see Cheng et al. 2018). We also assume that the coronal plasma is composed of fully ionized hydrogen with an electron density n_e = 10¹⁰ cm⁻³ and an average particle mass $\bar{m}=0.5{m}_{p}\,=8.36\times {10}^{-25}\,{\rm{g}}$, which gives a unit of mass density ${\rho }_{0}\,=2{n}_{e}\bar{m}=1.67\times {10}^{-14}\,{\rm{g}}/{\mathrm{cm}}^{3}$ (Priest & Forbes 2000). Here, m_p is the proton mass. The unit of magnetic strength is set as a typical coronal value, namely B₀ = 20 Gauss (also see Chen et al. 1999; Ye et al. 2020). Based on L₀, $\bar{m}$, ρ₀, and B₀, normalized units of other variables are deduced as ${u}_{0}={B}_{0}/\sqrt{{\mu }_{0}{\rho }_{0}}=4.36\,\times {10}^{7}\,\mathrm{cm}/{\rm{s}}$, t₀ = L₀/u₀ = 114.61 s, ${p}_{0}={\rho }_{0}{u}_{0}^{2}=3.18\,\mathrm{Pa}$, ${T}_{0}=\bar{m}{u}_{0}^{2}/{k}_{B}=11.52\,\mathrm{MK}$, ${\kappa }_{\parallel 0}={k}_{B}{\rho }_{0}{u}_{0}{L}_{0}/\bar{m}=6.02\,\times {10}^{11}\,\mathrm{erg}/{\rm{s}}\cdot \ \mathrm{cm}\cdot \ {\rm{K}}$, and η₀ = L₀ u₀ = 2.18 × 10¹⁷ cm²/s, where μ₀ is the magnetic permeability in vacuum and k_B denotes the Boltzmann constant. The ratio of specific heat is γ = 5/3. Under this configuration, the unit time, velocity, and temperature are also on the same order as the observational ones given by Cheng et al. (2018).

The initial magnetic field is set according to the CSHKP model (Chen et al. 1999; Ye et al. 2020)

$$\begin{eqnarray}\begin{array}{rcl}{B}_{x} & = & 0,\\ {B}_{y} & = & \left\{\begin{array}{ll}\sin \left(\displaystyle \frac{\pi x}{2\lambda }\right), & \left|x\right|\leqslant \lambda ,\\ 1, & x\gt \lambda ,\\ -1, & x\lt -\lambda ,\end{array}\right.\\ {B}_{z} & = & \sqrt{1-{B}_{x}^{2}-{B}_{y}^{2}},\end{array}\end{eqnarray} \tag{ 2 }$$

where λ denotes the half-width of the CS, which is set as 0.1. To approximate the chromosphere, transition region, and corona, the initial density distribution follows Takasao et al. (2015), namely,

$$\begin{eqnarray}\begin{array}{rcl}\rho \left(y\right) & = & {\rho }_{{chr}}+\displaystyle \frac{{\rho }_{{cor}}-{\rho }_{{chr}}}{2}\\ & & \left[\tanh \left(\displaystyle \frac{y-{h}_{{chr}}}{{l}_{{tr}}}+1\right)\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where ρ_cor = 1, ρ_chr = 10⁵, h_chr = 0.1, and ${l}_{{tr}}=0.02$. Fast reconnection is initially triggered by a localized anomalous resistivity, and thus the resistivity distribution is

$$\begin{eqnarray}\eta =\left\{\begin{array}{ll}{\eta }_{b}+{\eta }_{a}\exp \left[-\displaystyle \frac{{x}^{2}+{\left(y-{h}_{\eta }\right)}^{2}}{{l}_{\eta }^{2}}\right], & t\leqslant {t}_{\eta },\\ {\eta }_{b}, & t\gt {t}_{\eta },\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where η_a = 5 × 10⁻⁴, h_η = 0.5, l_η = 0.03, t_η = 5, and the background resistivity is uniformly set as η_b = 5 × 10⁻⁶ to simulate a high-Lundquist-number environment. Physically, the localized anomalous resistivity can be caused by microscopic instabilities in the CS such as lower-hybrid drift and/or ion acoustic instabilities, which can boost the process of reconnection (Yokoyama & Shibata 2001). Because η_a ≫ η_b , the anomalous resistivity dominates the evolution during t < t_η , which provides similar initial states at t = t_η for the cases listed in Table 1. The background pressure p_b is uniform to obtain a static initial state and is assigned as shown in Table 1. The boundary conditions are arranged as follows. The left (x = −1) and right (x = 1) are free boundaries, the top (y = 4) is the no-inflow boundary, and the bottom (y = 0) is the symmetric boundary. To reduce the influences of numerical boundaries, we only use data in the region defined by $x\in \left[-0.5,0.5\right]$ and $y\in \left[0,2\right]$ for analysis.

The above system is numerically solved with Athena 4.2 (Stone et al. 2008). We use the HLLD Riemann solver (Miyoshi & Kusano 2005), three-order piecewise parabolic flux reconstruction algorithm (Stone et al. 2008), and the Corner Transport Upwind method (Gardiner & Stone 2008) to solve the conservation part of Equation (1). The resistivity and TC are calculated by the explicit operator splitting method. To suppress numerical dissipation and get a uniformly high resolution (Shen et al. 2011; Ye et al. 2020), we set high-precision uniform Cartesian grids, namely 3840 and 7680 grids on the x and y directions, respectively. The corresponding pixel scales in two directions are Δx = Δy = 26 km. The maximum simulation time in our simulation is ${t}_{\max }=15$, which corresponds to 28.65 minutes in physical time.

3.1. Overview

3.2. Reconnection in the Main CS

We first estimate the reconnection rate in the main CS as shown by the area R₂ in Figure 1(b4). It is represented by the reconnection electric field at the principal X-point (see also Yokoyama & Shibata 2001), namely,

$$\begin{eqnarray}&&{ \mathcal R }\left(t\right)=\displaystyle \frac{\max \left({\eta }_{b}\left|{J}_{\mathrm{zxp}}\left(t\right)\right|\right)}{{B}_{\mathrm{in}}\left(t\right){u}_{\mathrm{in}}\left(t\right)},\end{eqnarray} \tag{ 5 }$$

where J_zxp denotes the out-of-plane current density at the X-points, $\max \left(\cdot \right)$ means taking the maximum value in the target CS region, and ${ \mathcal R }$ is normalized by the product of inflow magnetic strength B_in and Alfvén speed u_in, which are averaged in the adjacent region of the principal X-point defined by $x\in \left[-0.05,0.05\right]$ and $y\in \left[{y}_{\mathrm{xp}}-0.05,{y}_{\mathrm{xp}}+0.05\right]$ (see Figure 1(b2)). Here, y_xp denotes the y-coordinate of the principal X-point. The null points are determined by the method developed by Parnell et al. (1996).

For each simulation case, the reconnection rate in the main CS, ${{ \mathcal R }}_{\mathrm{cs}}$, shows a similar trend of evolution (Figure 2). To be specific, the average value of ${{ \mathcal R }}_{\mathrm{cs}}$ first rises quickly and then reaches its peak value, which is on the order of 0.01. We can see significant oscillations in the ${{ \mathcal R }}_{\mathrm{cs}}$ curves, which correspond to the formation of magnetic islands (Yokoyama & Shibata 2001). The overall evolution of the reconnection rate in the main CS is consistent with previous numerical results of 2D reconnection (e.g., Ni et al. 2012; Huang & Bhattacharjee 2016; Zenitani & Miyoshi 2020). Changes to β mainly affect the initial rising stage before t = 7. With the increase of β, the reconnection rate at this stage becomes smaller. However, when the reconnection is fully developed (t > 10), the evolution of ${{ \mathcal R }}_{\mathrm{cs}}$ is less affected by β. After including TC, the reconnection rate before t = 6 gets smaller and the average value also decreases slightly after t = 10 (Figure 2).

Zoom In Zoom Out Reset image size

3.3. Second Reconnection of Magnetic Islands

To study the loop-top annihilation of magnetic islands in detail, we analyze a typical event that appears at $t\in \left[11,11.65\right]$ in Case 4 (third column of Figure 1). As the magnetic island collides with the loop-top flux, a horizontal CS forms (see Figure 1(c3) and Figure 3(a)). Correspondingly, a pair of horizontal jets appear and quickly move along two opposite directions (Figure 3(b)). The peak speed of these jets reaches 0.8, which is of a similar order as the fastest downflow in the main CS. When these horizontal jets encounter the edge of the flare loop, they start to move downward, leading to an increase of downward speed u_y (Figure 3(c)), as observed by Tian et al. (2014). It is thus clearly revealed that the second reconnection enables the magnetic island annihilation at the flare loop top. During the second reconnection of islands, new small-scale magnetic islands are also generated in the loop-top CSs (see Figures 3(a) and (b)), similar to the "fragmenting coalescence" picture proposed by Bárta et al. (2011).

Zoom In Zoom Out Reset image size

To estimate the reconnection rate of the horizontal CS, we trace it using a box with a height of 0.008 and width of 0.05 (see Figure 3(a)). After locating all X-points in this region, we use Equation (5) to estimate ${{ \mathcal R }}_{\mathrm{lt}}$, the reconnection rate in the horizontal CS. It is shown that, as the magnetic island annihilates, the height of the horizontal CS keeps decreasing (Figure 3(d)). While for the reconnection rate within, it first increases, then reaches its peak value, and finally decreases (Figure 3(d)). Obviously, the peak reconnection rate is comparable with that in the main vertical CS. We also analyze a similar event in Case 1 where TC is not considered (see Figures 3(e)–(h)). It is found that the evolution of the horizontal CS, outjets, and the reconnection rate are similar to Case 4.

We further study the energy conversion during the second reconnection at the loop-top region and compare it with that in the main CS. In an area S, the effective changing rates of the magnetic, kinetic, and internal energies are, respectively, defined by (Forbes 1988)

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{P}}_{m} & = & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\displaystyle \int }_{S}\displaystyle \frac{{B}^{2}}{2}{\rm{d}}x{\rm{d}}y\\ & & +{\oint }_{\partial S}\left({\boldsymbol{E}}\times {\boldsymbol{B}}\right)\cdot {\boldsymbol{n}}{\rm{d}}l,\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{P}}_{k} & = & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\displaystyle \int }_{S}\displaystyle \frac{\rho {u}^{2}}{2}{\rm{d}}x{\rm{d}}y\\ & & +{\oint }_{\partial S}\displaystyle \frac{\rho {u}^{2}}{2}{\boldsymbol{u}}\cdot {\boldsymbol{n}}{\rm{d}}l,\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{P}}_{i} & = & \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\displaystyle \int }_{S}\displaystyle \frac{p}{\gamma -1}{\rm{d}}x{\rm{d}}y\\ & & +{\oint }_{\partial S}\displaystyle \frac{\gamma p}{\gamma -1}{\boldsymbol{u}}\cdot {\boldsymbol{n}}{\rm{d}}l,\end{array}\end{eqnarray} \tag{ 8 }$$

where ∂S denotes the boundary, E = η J − u × B is the electric field, n is the unit normal vector on the boundary, and dl denotes the line element. The negative value of ${\tilde{P}}_{m}$ implies the release of magnetic energy. We calculate the evolutions of ${\tilde{P}}_{m}$, ${\tilde{P}}_{k}$, and ${\tilde{P}}_{i}$ at the loop-top region, R₁, and that in the main CS, R₂ (see Figure 1(b4)). It is found that, during most of the time, the magnetic energy release is dominated by the main CS reconnection, especially during the early stage (Figure 4(a)). However, the second reconnection at the loop-top region is able to further release magnetic energy (see negative spikes of solid curves in Figure 4(a)). For example, in Case 4, during $t\in \left[11,11.65\right]$, the effective magnetic energy release rate in R₁ reaches about half of that in the main CS (Figure 4(a)). In R₂, the evolutions and magnitudes of ${\tilde{P}}_{k}$ and ${\tilde{P}}_{i}$ are similar (see dashed curves in Figures 4(b) and (c)), which shows that magnetic energy keeps converting to kinetic and internal energies. However, in R₁ the energy conversion process is relatively complicated. When a magnetic island enters R₁, its velocity is damped by the magnetic tension force of the flare loop. As a result, during the second reconnection, ${\tilde{P}}_{k}$ is mostly negative (see solid curves in Figure 4(b)); the magnetic energy and kinetic energy significantly increase the internal energy (see positive peaks of ${\tilde{P}}_{i}$ in Figure 4(c)). Furthermore, the time-integrated effective changes of magnetic, kinetic, and internal energies during $t\in \left[5,15\right]$ are calculated by ${\tilde{W}}_{m}={\int }_{5}^{15}{\tilde{P}}_{m}{\rm{d}}t$, ${\tilde{W}}_{k}={\int }_{5}^{15}{\tilde{P}}_{k}{\rm{d}}t$, and ${\tilde{W}}_{i}={\int }_{5}^{15}{\tilde{P}}_{i}{\rm{d}}t$, respectively. For Case 4, the ${\tilde{W}}_{m}$ released in R₁ is 13.7% of that in R₂ (see Table 2), which shows that the total release of magnetic energy during solar flares is dominated by the reconnection in the main CS region and further replenished by that at the loop-top region. In R₂, about half of the released magnetic energy is converted to internal energy and half becomes kinetic energy. However, in R₁ the main increment of the internal energy (about 72.5%) is due to magnetic energy release and the rest is from the damping of kinetic energy.

Zoom In Zoom Out Reset image size

Table 2. Time-integrated Increments of Magnetic, Kinetic, and Internal Energies in Regions R₁ and R₂ for Case 4

Region	R1	R2
${\tilde{W}}_{m}$	−0.0153	−0.1115
${\tilde{W}}_{k}$	−0.0091	0.0527
${\tilde{W}}_{i}$	0.0211	0.0582

Note. The normalizing unit of ${\tilde{W}}_{m}$, ${\tilde{W}}_{k}$, and ${\tilde{W}}_{i}$ is W₀ = 7.96 × 10²⁰ erg cm⁻¹. Summations of ${\tilde{W}}_{m}$, ${\tilde{W}}_{k}$, and ${\tilde{W}}_{i}$ in both areas should be zero theoretically but small deviations may be caused by numerical errors.

Download table as:  ASCIITypeset image

The magnetic energy cascading at the loop-top region is closely related to the distribution of the corresponding magnetic islands. The loop-top magnetic islands generated during the second reconnection can be characterized by the evolution of the cumulative distribution function (CDF) of island flux $N\left(\psi ,t\right)$, which measures the number of islands with flux ≥ ψ at moment t (Figure 4(d)). When the downflow magnetic islands initiate the second reconnection, as marked by negative peaks of ${\tilde{P}}_{R1}$, the CDF of the island flux is significantly enhanced (see duration $t\in \left[11,11.65\right]$ in Figure 4(d) as an example). In Figure 5, we analyze the probability distribution functions (PDFs) of magnetic island area and flux at the loop-top region R₁ based on the same method widely used for calculating plasmoid PDFs in CS (e.g., Shen et al. 2013; Lynch et al. 2016; Ye et al. 2019). The PDF of the magnetic island flux, $f\left(\psi \right)$, is determined by $f\left(\psi \right)={\rm{d}}{N}_{\tau }\left(\psi \right)/{\rm{d}}\psi$, where ${N}_{\tau }\left(\psi \right)=\int N\left(\psi ,t\right){\rm{d}}t$ and the integration is taken over duration $t\in \left[5,15\right]$. We also normalize the PDF to satisfy $\int f\left(\psi \right){\rm{d}}\psi =1$. The procedure is similar for calculating $f\left(S\right)$, the PDF of the island area. For most of the sample domain, namely S < 10³ and ψ < 10⁻³, both $f\left(S\right)$ and $f\left(\psi \right)$ follow a power law independent to β and TC. The slopes of the PDFs vary from −1 to −2 as the island scale increases. For the large-scale end, the PDFs slightly vary for different cases, but may be unreliable because the number of large-scale islands we collected is very limited.

Zoom In Zoom Out Reset image size

3.4. Formation of Termination Shocks

As predicted by the standard model, the reconnection downflows may drive a TS (Forbes & Priest 1983; Forbes 1986; Forbes & Malherbe 1986, 1991; Yokoyama & Shibata1997, 1998, 2001; Shen et al. 2018; Kong et al. 2019, 2020). Here, we further study the influences of β and TC on the formation of TSs. Because the TS is highly dynamic, in particular, it could be fragmented as the magnetic islands pass (see Shen et al. 2018); we mainly focus on the principal TS, namely the largest and strongest segment of the TS front (see Figures 3(c) and (g)). Without TC, the scale and strength of the principal TSs both decrease with the increase of β (Figure 6). Further, after TC is considered, although the histogram of M_TS is similar to that without TC (see the blue and purple curves in Figure 6(b)), the distribution of L_TS concentrates toward the lower end. This shows that an increase of β and TC can suppress the scale and strength of TSs.

Zoom In Zoom Out Reset image size

In this paper, we study the loop-top annihilation of magnetic islands formed in the main vertical CS in detail. Besides the reconnection in the main eruption CS, we find that reconnection also occurs at the flare loop top, as illuminated by Bárta et al. (2011). The corresponding reconnection rate is found to be of the same magnitude as that in the main CS. Interestingly, this process is highly similar to the magnetic reconnection taking place in the Earth's magnetotail (Wang et al. 2010, 2016; Lu et al. 2015). We thus suggest that magnetic island annihilation at the flare loop top, as revealed here but neglected in the standard flare model, plays a nonnegligible role in further releasing magnetic energy so as to heat flare plasma and accelerate particles.

The fast magnetic reconnection in the main CS is dominated by the generation of plasmoids, which subsequently move quickly upward and downward along the CS. This is largely in agreement with previous results (e.g., Bhattacharjee et al. 2009; Shen et al. 2011, 2018; Kong et al. 2019, 2020) and implies the existence of fine structures in the CS during solar eruptions. The second reconnection, characterized by horizontal CSs, is initiated when magnetic islands collide with the flux above the flare loops. As the magnetic islands move toward the flare loop, the reconnection rate first grows and then decreases, with the peak being comparable with that in the main CS. The second reconnection also generates smaller-scale secondary magnetic islands, which enhance the turbulence and thus further facilitate the release of magnetic energy at the flare loop top. Moreover, the whole reconnection process seems to be intermittent because of the presence of magnetic islands; it is thus expected that the acceleration of nonthermal particles is intermittent. This might be used to explain the QPPs observed during flares.

The cascading law of magnetic islands at the loop top seems to be less affected by background β or TC. The PDF of the magnetic island flux generated during the second reconnection follows a power law varying from $f\left(\psi \right)\sim {\psi }^{-1}$ to ψ⁻², similar to the PDF of the magnetic island flux in the main CS obtained by Loureiro et al. (2012), Shen et al. (2013), and Lynch et al. (2016). However, for the spatial scale and strength of the TSs formed by the downward outflows or plasmoids, it is found that both decrease with an increase of plasma β or when including TC.

It is worth mentioning that simulations of magnetic reconnection depend on the spatial resolution. The results presented here are obtained under a high-resolution spatial grid (3840 × 7680), which can realize simulations with background S = 10⁶, as discussed by Ye et al. (2020). The thickness scales of CSs in both R₁ and R₂ are on the order of 0.01, which is 0.5% of the x-domain (see Figure 1 and Figure 3). Under the current resolution, the CS thickness can be resolved by about 20 grids. Furthermore, the scale of the internal singular layer in classic tearing mode theory is $\delta \sim {{aS}}_{a}^{-1/4}$ (Pucci et al. 2017), where a denotes the width of a CS, S_a = av_A /η is the local Lundquist, and v_A is the Alfvén speed. The number of grids resolving the inner singular layer can be estimated by

$$\begin{eqnarray}&&\displaystyle \frac{\delta }{{\rm{\Delta }}L}\sim \displaystyle \frac{{{aS}}_{a}^{-1/4}}{{\rm{\Delta }}L}=\displaystyle \frac{{\eta }^{1/4}}{{\rm{\Delta }}L}\displaystyle \frac{{a}^{3/4}}{{v}_{A}^{1/4}},\end{eqnarray} \tag{ 9 }$$

where ΔL denotes the dimensionless size of the mesh grid. Setting ΔL = 2/3840, η = 5 × 10⁻⁶, v_A = 1, and a = 0.01, we have δ/ΔL ∼ 3. It is shown that the small-scale inner singular layer can still be roughly resolved. We have also studied the influence of resolutions on our conclusions by running cases with different grid configurations. It is shown that the conclusions are only valid for relatively high-resolution simulations.

In our numerical model, the fine dynamics of the chromosphere might be inaccurate, because the effects of partially ionized plasma (PIP) are not included. It is true that partial ionization has influences on the properties of shocks, the CS structure, and fractal reconnection (Imada et al. 2011; Hillier et al. 2016). However, in this paper we mainly focus on the reconnection dynamics in the high-temperature corona. Thus, the effects of the partially ionized chromosphere on our results are very limited; the only influence is on the structure and dynamics of the post-flare loops into which the chromospheric plasma is evaporated.

We would like to thank the anonymous referee for valuable suggestions. This research is supported by Natural Science Foundation of China grant Nos. 11722325, 11733003, 11790303, 11790300, and 11805203.