© S. Sen and R. Keppens 2022

1. Introduction

Magnetic reconnection is ubiquitous in both laboratory and astrophysical plasmas, where the change in magnetic field topology leads to the conversion of the magnetic energy into thermal and kinetic energies (Biskamp 2000). Magnetic reconnection plays a key role in the fast energy release in solar flares (Giovanelli 1939, 1947, 1948; Priest & Forbes 2000; Hesse & Cassak 2020), in coronal mass ejections into the ambient solar wind medium (Gosling et al. 1995; Schmidt & Cargill 2003; Karpen et al. 2012), and mediates the loss of plasma confinement in tokamak experiments (Günter et al. 2015). The rearrangement of the magnetic field topology occurs in a localized plasma region where non-ideal magnetohydrodynamics (MHD) effects dominate, breaking the frozen-in condition. The Sweet–Parker model (Parker 1957; Sweet 1958) predicts the reconnection rate to scale with the Lundquist number (S_L = Lv_A/η) as $S_{\mathrm{L}}^{-1/2}$, where L is the characteristic length, v_A is the Alfvén velocity, and η is the resistivity of the medium. However, this prediction is too slow to agree with reconnection observations for the solar atmosphere. On the other hand, the Petschek model (Petschek 1964) estimates the reconnection rate to scale as (log S_L)⁻¹. Biskamp (2000) and Huang & Bhattacharjee (2010) state that the Petschek reconnection is achievable only if the local resistivity of the current sheet is enhanced, while Baty et al. (2009) report the occurrence of Petscheck like reconnection for a low uniform resistivity (10⁻³).

The simplest configuration susceptible to magnetic reconnection is a single current layer model formed by a polarity reversal of the magnetic field. This reconnection may be triggered due to the growth of a classical linear resistive instability, known as tearing mode instability (Furth et al. 1963). A current sheet of aspect ratio L/δ ≳ 2π (where L and δ are the characteristic length and thickness of the current layer respectively) can develop magnetic islands due to the growth of linearly unstable perturbations. Pucci & Velli (2014), Landi et al. (2015), and Tenerani et al. (2016) have reported the development of ideal tearing modes in current sheets for a large aspect ratio of $\sim S_{\mathrm{L}}^{1/3}$. On the other hand, double current layer models are also seen to give rise to resistive instabilities known as double tearing modes (DTMs). Double current layer models with two widely seperated current layers can develop a single standard tearing mode on each layer that will influence each other later in the nonlinear evolution stage (Keppens et al. 2013; Paul & Vaidya 2021). In contrast, DTMs are tearing modes that are intrinsically coupled and that co-develop on nearby resonant surfaces. The evolution of DTMs in the nonlinear regime has been seen to lead to an explosive reconnection and a weak dependence on the resistivity (Zhang & Ma 2011; Akramov & Baty 2017 and references therein). DTMs have been studied under various important conditions like external shear flows (Ofman 1992; Bierwage et al. 2007; Wang et al. 2008), bootstrap current (Yu 1997), anomalous electron viscosity (Dong et al. 2003), collisionless plasma (Bierwage & Yu 2007), and Hall effects (Zhang et al. 2008; Zhang & Ma 2009). Pucci et al. (2018) estimated scaling relations between the maximum growth rate, Lundquist number and the aspect ratio, for linear tearing modes in a double current sheet setup.

Thermal instability is an essential mechanism to form condensations in the solar atmosphere. The theory is laid out in classical treatments by Parker (1953) and Field (1965). These works explain how a runaway process of radiative cooling leads to thermal instability (TI) in plasma. The solar corona can be considered to be in a delicate thermal equilibrium balancing the optically thin radiative loss and background heating in combination with thermal conduction. If this balance is perturbed, and the increment of radiative loss cools down the plasma, isobaric, isentropic or isochoric evolutions may self-amplify the radiative losses. This drives the enhancement of the local plasma density, which further increases the energy loss by radiation (because radiative energy loss in an optically thin medium varies with the density squared), which in turn decreases the temperature even more. Hence, a catastrophic runaway process results in a rapid drop in temperature and an increase in plasma density. Field (1965) reported a detailed analysis of the thermal instability in an infinite homogeneous medium that triggers the catastrophic radiative cooling. Later, the analysis was extended to non-uniform slab geometry (van der Linden & Goossens 1991b; van der Linden et al. 1992), and cylindrical flux tubes with solar coronal conditions (van der Linden & Goossens 1991a; Ireland et al. 1995; Soler et al. 2011). The thermal instability theory can be invoked to explain various fascinating features of the solar atmosphere. For example, Smith & Priest (1977) discussed the formation of a solar prominence in a current sheet, and Priest & Smith (1979) extended the study for prominence formation to solar coronal arcades. Forbes & Malherbe (1991) used 2D MHD simulation in a line-tied current sheet with the effect of radiative energy losses to explain the formation of post-flare loops. The ab initio formation of a solar prominence due to chromospheric evaporation and thermal instability is shown by Xia et al. (2012) in a 2.5D simulation, and the dynamical evolution of similar prominence setups is further explored in Keppens & Xia (2014). Later, Xia & Keppens (2016a) developed a 3D model of prominence formation due to the plasma cycle between corona and chromosphere, while the complex 3D dynamics in a twin-layer prominence is reported by Xia & Keppens (2016b). More recently, linear and nonlinear stability analysis of thermal instability due to the interaction of the entropy and slow MHD wave modes is discussed in Claes & Keppens (2019) and Claes et al. (2020), and the effect on the thermal instabilities arising from different radiative cooling models is reported by Hermans & Keppens (2021). The formation of fine structures in the prominence may well relate to the linear magnetothermal modes affected by anisotropic thermal conduction with finite (albeit small) conduction across the magnetic field lines (van der Linden & Goossens 1991a). Recent solar applications include the works by Li et al. (2022), who studied the formation of coronal rain due to the thermal instability of randomly heated arcades; by Jenkins & Keppens (2021), who studied the formation of prominences due to levitation-condensation; and by Zhao & Keppens (2022), who developed a novel plasmoid-fed prominence formation scenario during flux rope eruption. In the latter work, chromospheric plasma collects into a current sheet, which ultimately shows chaotic plasmoid formation, where the cool chromospheric plasma gets trapped and lifted into an erupting prominence structure. This study motivated our current work, where we will investigate more rigorously how thermal and tearing effects can reinforce each other.

In this work we study the tearing and thermal instabilities of a single current sheet model with the non-adiabatic effects of radiative energy loss and background heating in a resistive 2D MHD simulation. Recent theoretical studies by Ledentsov (2021a,b,c) show that the instability growth rate in the linear regime of a tearing mode is modified when the non-adiabatic effects, radiative energy loss, electrical and thermal conductivities are incorporated. This motivated us to explore the growth rate of a thermally influenced tearing mode in the linear and nonlinear domains by means of MHD simulation. The focus of our idealized study is mainly devoted to estimating the characteristic growth rate timescales in the linear and nonlinear regimes and to finding scaling relations with the resistivity in the different phases of the nonlinear evolution. We also quantify the number of generated plasmoids as influenced by different physical parameters, and its scaling relation with the Lundquist number. In various earlier works, Zhang & Ma (2011), Keppens et al. (2013), Akramov & Baty (2017), Paul & Vaidya (2021) have studied the tearing mode effect in explosive reconnection events for adiabatic conditions. By including the radiative energy loss in an optically thin medium and background heating to investigate the effect of thermal instability in a current sheet model, we determine whether explosive reconnection may be triggered in a greater number of Lundquist regimes than in pure resistive MHD alone.

The paper is organized as follows. In Sect. 2 we describe the model setup and the numerical framework along with the initial and boundary conditions. In Sect. 3 the main results of the study and its analysis are reported. Section 4 discusses the significance of the work for a typical coronal medium, summarizes the key findings, and finally concludes how our findings may be useful for future studies.

2. Numerical setup of the model

To understand the evolution of the thermal instability in association with the tearing mode, we use a resistive 2D MHD simulation using the MPI-parallelised Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC)¹ (Keppens et al. 2012, 2021; Porth et al. 2014; Xia et al. 2018). The spatial domain of the simulation setup is −12.8 × 10⁴ km to 12.8 × 10⁴ km along the x- and y-direction with maximum five levels of adaptive mesh refinement (AMR) between −2 × 10⁴ km ≤ y ≤ 2 × 10⁴ km and −12.8 × 10⁴ km ≤ x ≤ 12.8 × 10⁴ km, effectively achieving a maximum spatial resolution of 2048 × 2048, which makes the smallest cell 125 km in size. The triggering of the (de)refinement is based on the errors estimated by the density (gradient) of an instantaneous time step. To explore the non-adiabatic effects in the evolution of the tearing instability in the resistive MHD regime, we solve the following normalised MHD equations numerically:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\mathbf{v}\rho )& =0,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial (\rho \mathbf{v})}{\partial t}+\mathrm{\nabla }·(\rho \mathbf{v}\mathbf{v}+p_{\mathrm{tot}}\mathbf{I}-\mathbf{B}\mathbf{B})& =0,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial \mathcal{E}}{\partial t}+\mathrm{\nabla }·(\mathcal{E}\mathbf{v}+p_{\mathrm{tot}}\mathbf{v}-\mathbf{B}\mathbf{B}·\mathbf{v})& =\eta {\mathbf{J}}^2-\mathbf{B}·\mathrm{\nabla }\times (\eta \mathbf{J})\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & -{\rho }^2\mathrm{\Lambda }(T)+H_{\mathrm{bgr}},\hfill \\ \hfill \frac{\partial \mathbf{B}}{\partial t}+\mathrm{\nabla }·(\mathbf{v}\mathbf{B}-\mathbf{B}\mathbf{v})+\mathrm{\nabla }\times (\eta \mathbf{J})& =0,\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \mathrm{\nabla }·\mathbf{B}& =0,\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill \mathbf{J}& =\mathrm{\nabla }\times \mathbf{B}.\hfill \end{array}$$(6)

Here I is the unit tensor, and the quantities ρ, T, B, v, and η have their usual meanings of mass density, temperature, magnetic field vector, velocity, and resistivity respectively. The total pressure, p_tot, adds the plasma and the magnetic pressure

$$\begin{array}{c}\hfill p_{\mathrm{tot}}=p+\frac{B^2}{2},\end{array}$$(7)

where, the thermodynamic quantities are linked through the ideal gas law: p = ρk_BT/μm_H, where p is the plasma pressure, k_B is the Boltzman constant, μ = 0.6 is the coronal abundance due to the fully ionized plasma of H and He atoms with the abundance ratio of 10 : 1, and m_H is the proton mass. The total energy density is given by

$$\begin{array}{c}\hfill \mathcal{E}=\frac{p}{{\gamma }_{\mathrm{gc}}-1}+\frac{\rho v^2}{2}+\frac{B^2}{2},\end{array}$$(8)

where, γ_gc = 5/3 is the ratio of specific heats for the monoatomic gas. The solenoidal condition of magnetic field, and the current density, J are given by the Eqs. (5) and (6) respectively. The non-adiabatic effect due to the radiative cooling of the optically thin medium, which is relevant for the solar corona, is incorporated by the third term on the right-hand side of Eq. (3). The optically thin cooling due to the radiation depends on the local density, and the temperature-sensitive cooling model Λ(T). In this work, we use the combined cooling model developed by Colgan et al. (2008) and Dalgarno & McCray (1972), which we call the Colgan_DM model, shown in Fig. 1. The details of different radiative cooling curves and their effects on the formation of condensations are reported in Hermans & Keppens (2021). In order to maintain a thermal equilibrium in the initial state, we use the background heating, H_bgr (last term on the right-hand side of Eq. (3)) in such a way that it compensates the radiative loss at the initial state. Hence, we take the background heating as

$$\begin{array}{c}\hfill H_{\mathrm{bgr}}={\rho }_i^2\mathrm{\Lambda }(T_i),\end{array}$$(9)

Fig. 1. Radiative cooling curve due to the Colgan_DM model.

where ρ_i and T_i are the equilibrium density and temperature, respectively, and therefore the background heating is constant with time (but varies in space as explained in the following section). We note that for simplicity here we study an idealized current sheet setup, where we ignore gravity and the role of (anisotropic) thermal conduction.

To study the long-term behavior of the evolution of a current layer configuration subjected to the resistive MHD and non-adiabatic effects of radiative cooling, we set a 2D square domain in the x − y plane that spans from −12.8 to 12.8 (in dimensionless unit) along the x- and y-direction. The unit density, temperature, and length scale that serve to normalize the simulation are set as $\overline{\rho }=2.34\times {10}^{-15}$ g cm⁻³, $\overline{T}={10}^6$ K, and $\overline{L}={10}^9$, cm respectively, which are typical values for the solar corona. The magnetic field in the initial setup is taken as a non-force-free planar field given by

$$\begin{array}{cc}& B_x=B_0\tanh (y/l_s),\hfill \end{array}$$(10)

$$\begin{array}{cc}& B_y=0.\hfill \end{array}$$(11)

This implies that the B_x component realizes a polarity inversion around y = 0, and hence a current sheet is formed at y = 0 that is in the range $x=[-12.8,12.8]\times \overline{L}$. We set the background field amplitude, B₀ = 1 (which corresponds to 2 G in physical unit), and l_s = 0.5, which sets the total width of the current sheet to 2l_s. The planar magnetic field B_x(y) is further perturbed to trigger some tearing-type evolution, as

$$\begin{array}{cc}& \delta B_x=-\frac{2\pi {\psi }_0}{l_y}\cos (\frac{2\pi x}{l_x})\sin (\frac{2\pi y}{l_y}),\hfill \end{array}$$(12)

$$\begin{array}{cc}& \delta B_y=+\frac{2\pi {\psi }_0}{l_x}\sin (\frac{2\pi x}{l_x})\cos (\frac{2\pi y}{l_y}),\hfill \end{array}$$(13)

where the geometric parameters $l_x=l_{\mathit{y}}=25.6\times \overline{L}$ match the domain sizes of the simulation, and the perturbation amplitude ψ₀ = 0.1 is 10% of the magnetic field amplitude B₀. It should be noted that Eqs. (12) and (13) satisfy the condition ∇ ⋅ δB = 0. The initial density profile is taken as

$$\begin{array}{c}\hfill {\rho }_i={\rho }_0+{\cosh }^{-2}(y/l_s),\end{array}$$(14)

where ρ₀ = 0.2 (which corresponds to 4.68 × 10⁻¹⁶ g cm⁻³ in physical units) is the density outside the current sheet. The initial density profile has the peak value of 2.808 × 10⁻¹⁵ g cm⁻³ at y = 0 and gradually converges to ρ₀ for |y|> 0. The initial equilibrium temperature T_i is taken as a constant value of 0.5 MK throughout the simulation domain in order to fix the plasma-β = 0.2, which is below unity, as is appropriate for the solar corona. The initial variation in density, together with the uniform temperature, realizes a pressure variation p(y) that exactly balances the Lorentz force associated with the field B_x(y). This implies that in a simulation where no magnetic field perturbation is applied, and where the resistivity is set to zero (ideal MHD), we actually have a force-balanced and thermally balanced environment. In our simulations below, resistive effects modify the temperature and hence drive the system away from the thermal equilibrium balance between losses and heating H_bgr. We note that the finite resistivity is crucial to allow tearing. If we simulate without non-adiabatic effects included, and just evolve the system in resistive MHD at the used constant and uniform resistivity values, we get an evolution toward a standard reconnection experiment with a central Sweet–Parker-type current sheet within a growing island structure (at both periodic sides).

After this initial setup the system evolves as governed by Eqs. (1)–(6). This set of equations is solved using a three-step Runge-Kutta time integration with a third-order slope limited reconstruction method (Čada & Torrilhon 2009), and a Harten-Lax-van Leer (HLL) flux scheme (Harten et al. 1983). As we encounter fairly extreme density and temperature contrasts in the evolution, we need to enforce an automated recovery procedure to ensure positivity throughout, and we do so by fixing the minimum pressure and density values equal to 10⁻¹⁴ and 10⁻¹² code units, respectively. We follow the current sheet evolution for up to 2550 s (42 min), and usually save data at a cadence of 42.5 s, which gives 60 simulation snapshots. We use periodic boundary conditions in the x-direction, and open boundary conditions in the y-direction. The typical wall clock time of a single run is ≈20 h for parallel computation with 15 CPUs.

3. Results and analysis

3.1. Global evolution

The spatial distribution of the absolute current density, |J_z| is shown in Fig. 2 for four different time stages, with plasma-β  =  0.2 and resistivity η  =  0.001 (equivalent to 1.2  ×  10¹⁴ cm² s⁻¹ in physical units). Figure 2a represents the initial configuration of the current sheet, which is located around y = 0 and extends all along the x-direction. Figure 2b and its animated view clearly shows that the current sheet narrows as a result of the thermodynamic evolution driven by radiative losses and quickly thereafter fragments, forming a pronounced chain of many small-scale plasmoids due to the combination of thermal and tearing instabilities. We calculate the force (per unit volume) along the y-direction, ${\displaystyle F_{\mathit{y}}=-\frac{\partial p}{\partial \mathit{y}}+B_x(\frac{\partial B_{\mathit{y}}}{\partial x}-\frac{\partial B_x}{\partial \mathit{y}})}$. At t = 0, the system maintains the equilibrium condition and hence F_y = 0, but when the system evolves, the equilibrium is violated due the thermally influenced tearing mode instability, where we see F_y has a dominating positive force (upward direction) for y <  0, and a dominating negative force (downward direction) for y >  0 in the vicinity of the current sheet at x = 0 between y = ±0.75 × 10⁴ km at t = 300 and 429 s before the fragmentation stage of the current sheet (see Fig. 3). Therefore, the current sheet is squeezed along the y-direction. As time progresses, the small islands coalesce and merge with each other at later times (see Figs. 2c and d). The evolution of the plasma density, ρ, at the same stages is shown in Fig. 4. The initial configuration of the density distribution shown in Fig. 4a has an enhanced density region present in the vicinity of y = 0 according to Eq. (14). In the later stage, at t = 900 s (Fig. 4b), the current sheet fragments into smaller plasma blobs, and these merge with each other at the later stages shown in Figs. 4c and d, respectively. The same instants but seen in the temperature evolution are shown in Fig. 5. Comparing Figs. 4b, c, and d with Figs. 5b, c, and d shows that the temperature depletion regions are formed at the regions where the plasma materials are condensed. We calculate the radiative loss at x = 0 between y = ±1.25 × 10⁴ km, which covers the entire vertical domain of the current sheet. At t = 0, the radiative loss is equal to the background heating, H_bgr, which is constant with time. However, when the system evolves due to the thermally influenced tearing mode, the radiative loss term dominates over the background heating within the selected region. This is shown in Fig. 6a, where the radiative loss dominates over H_bgr at t = 429 s, which leads to a temperature drop to 21 000 K within that region as shown in Fig. 6b. The spatial distribution of the density on a horizontal cut all along the current sheet (y = 0) is shown in Fig. 7, for two different evolution stages (t = 900 and 2447 s). Due to the relative motions of the plasmoids along the ±x-direction, they merge to form denser plasmoids, separated by density depletion regions, which is reflected in Fig. 7. The temporal evolution of the instantaneous minimum and maximum density and temperature, as computed over the entire domain are shown in Fig. 8. These extrema are virtually always encountered along the central current sheet. It shows that the minimum density (ρ_min) falls gradually with the evolution, whereas the maximum density (ρ_max) increases with time, which means that the overall density contrast of the medium increases with time. It is estimated that the ratio of ρ_max to ρ_min changes by a factor of ≈10¹ to ≈10⁶ between the initial and the final times. The equilibrium temperature of the initial setup is constant (0.5 MK) for the entire simulation domain, and the minimum temperature during the evolution is not allowed to drop below 1000 K, which is well above the minimal temperature of the exploited cooling curve. Hence, the maximum (T_max) to minimum (T_min) temperature ratio of the medium is unity at the initial stage, and increases to a factor of ≈10⁴ at the final time. The local thermodynamic analysis for some typical, selected plasmoids are shown in Figs. 9–11. The zoomed-in version of a local plasmoid during the evolution stage for two different times, t = 900 and 2447 s, are shown in Figs. 9b and d, respectively, where the variation in the mass densities along the horizontal cuts are indicated, and shown in Figs. 10c and f, respectively. The maximum density enhancement along the line for t = 900 s is ∼10⁻¹³ g cm⁻³ (Fig. 10c), whereas for t = 2447 s, the value increases to ∼10⁻¹² g cm⁻³ (Fig. 10f). This implies that the condensations formed by thermal instability within smaller-scale plasmoids are collected into larger condensed regions, as the plasmoids merge by coalescence. Similarly, the temperature variation in the same zoomed-in region is shown in Fig. 10. It shows that the minimum temperature within this domain is ∼10⁴ K, and the temperature depletion regions correspond to the overdense regions of Figs. 9a and c, respectively. The radiative loss for an optically thin medium, which depends on the local density and temperature (ρ²Λ(T)) is shown in Fig. 11. The variation in the radiative loss along the horizontal lines of Figs. 11b and e are shown in Figs. 11c and f, respectively. It should be noted in Figs. 10b and 11b (or in Figs. 10e and 11e) that the radiative loss of the temperature depleted region is greater than in its surroundings. This is because of the density enhancement of that region compared to its surroundings (see Fig. 9b or d), and the ρ² term dominates over the Λ(T) resulting in the increase in radiative loss, ρ²Λ(T). The condensations that are entrapped within the coalescing plasmoids thereby also show rapid variations in the radiative losses across their edges, much like the prominence corona transition region (PCTR), also at play for individual coronal rain blobs. We note that the precise variation in temperature, density, and radiative losses here do not incorporate effects of anisotropic thermal conduction, and therefore we have extremely sharp transitions, also explained in Hermans & Keppens (2021).

Fig. 2. Distribution of the absolute current density, |Jz|, normalized with the instantaneous absolute peak current density, |Jmax|, for different evolution stages. The y-domain is only shown in the range −1.45 × 104 km ≤ y ≤ +1.45 × 104 km, which contains the region of the current sheet. An animation of the figures is available online.

Fig. 3. Variation in ${\displaystyle F_{\mathit{y}}=-\frac{\partial p}{\partial \mathit{y}}+B_x(\frac{\partial B_{\mathit{y}}}{\partial x}-\frac{\partial B_x}{\partial \mathit{y}})}$ along the y-direction in the range y = ±0.75 × 104 km for three different times (t = 0, 300, and 429 s) before the fragmentation stage of the current sheet.

Fig. 4. Same as Fig. 2, but for plasma density, ρ, for different evolution stages. The overplotted white lines represent the magnetic field lines. An animation of the figures is available online.

Fig. 5. Distribution of the temperature, T, within the same domain as Fig. 2, for different evolution stages. An animation of the figures is available online.

Fig. 6. Distribution of (a) radiative cooling and background heating, Hbgr, and (b) temperature, T, along the vertical cut at x = 0 in the range y = ±1.25 × 104 km for t = 0 and 429 s. The minimum temperature obtained within this regime at t = 429 s is Tmin = 21 000 K.

Fig. 7. Spatial distribution of plasma density along the y = 0 cut for two different evolution stages.

Fig. 8. Temporal evolution of the instantaneous maximum and minimum densities (left) and temperatures (right) within the entire simulation domain.

Fig. 9. Analysis of density distribution of local plasmoids. Top panel: (a) Distribution of the density for t = 900 s, within the same domain as Fig. 2; (b) zoomed-in version of the selected region. Bottom panel: Same as top panel, but for t = 2447 s.

Fig. 10. Analysis of temperatutre distribution of local plasmoids. Top panel: (a) Distribution of the temperature for t = 900 s, within the same domain as Fig. 2; (b) zoomed-in version of the selected region; (c) Density and temperature distributions along the horizontal cuts indicated by the dashed lines in Figs. 9b and 10b, respectively. Bottom panel: same as top panel, but for t = 2447 s.

Fig. 11. Same as Fig. 10 but for radiative loss in an optically thin medium.

3.2. Growth rate and scaling relation

The velocity distribution of the plasma motion is shown in Fig. 12, where the magnitude is scaled with respect to the Alfvén velocity, v_A, which is measured based on the magnetic field strength, B₀ = 2 G, and the mass density of the equilibrium current sheet, ρ_c = 2.81 × 10⁻¹⁵ g cm⁻³. The Alfvén timescale is measured by $t_A=\overline{L}/{\mathit{v}}_A$, where, $\overline{L}$ is the unit length of 10⁹ cm. It is evident from Figs. 12a, b and c that the velocity, v_x remains localized in the vicinity of the current sheet (y = 0). We note that in line with our initial single-island magnetic field perturbation, we see a pronounced rightward motion in the right half of the domain, and a leftward motion in the left half. We later see typical Petschek-like signatures in the flow fields between islands, especially about the middle x = 0, with super-Alfvénic outflow speeds bounded by slow shocks. Figure 13a represents the evolutionary nature of a current sheet in adiabatic and non-adiabatic conditions. It is clear from the figure that the instantaneous maximum velocity growth for the non-adiabatic case is more rapid than for the adiabatic conditions. The evolutionary behavior of the current sheet configuration due to thermal and tearing instabilities is shown by the black curve in Fig. 13a for plasma-β = 0.2, and a given resistivity value, η = 0.001, while the evolution for different η values are shown in Fig. 13b. As a diagnostic measurement of the instability, we determine the evolution of the instantaneous maximum absolute velocity, |v_x|_max. Figure 13a shows that this evolution exhibits three distinct phases: (i) the early phase (between t = 0 and 250 s), where the velocity growth occurs exponentially (linearly on the logarithmic-linear scale), which is called the linear growth regime; (ii) the middle phase between t = 250 and 665 s, where the growth rate is slower compared to the linear phase, which is called the Rutherford regime (Rutherford 1903), and (iii) the final phase, which starts at t = 665 s, where the instability suddenly develops in an explosive way, and finally saturates at a later time, which we call the post-Rutherford regime. To infer the evolution rates quantitatively for all the different phases, we calculate the growth rates by scaling them with respect to the Alfvén timescale, t_A. We define the growth rate as γ = d(ln(|v_x|_max))/dt. To estimate the linear growth rate, γ_lin, we calculate the growth rate in the linear regime by taking the mean value of the slope, which gives ${\gamma }_{\mathrm{l}\mathrm{i}\mathrm{n}}=3.76\times {10}^{-1}{t}_A^{-1}$. This value is larger by an order of magnitude compared to the studies of the double current sheet problem (Otto & Birk 1992; Zhang & Ma 2011; Akramov & Baty 2017; Paul & Vaidya 2021), where the radiative cooling effect (or other non-adiabatic effects, e.g., thermal conduction) is not incorporated. This implies that the higher linear growth rate can be ascribed to the non-adiabatic effects of the radiative cooling and background heating. This is also in agreement with our own study for a single current layer model reflected in Fig. 13a, where the average growth rate for the adiabatic medium is lower than the non-adiabatic case. Similarly, we estimate the average growth rates for the Rutherford regime (γ_Ruth) and the post-Rutherford regime (γ_PR) for different resistivity values within the range of η = 0.0001 to 0.005. The velocity evolution for some selected resistivity values are shown in Fig. 13b. This shows that the explosive phase of the evolution starts at later times for higher resistivity values, and converges at the final stage. For a Sweet–Parker type current sheet (where the inverse aspect ratio of the current sheet follows the scaling relation $l_s/L\sim S_{\mathrm{L}}^{-1/2}$), the thickness of the current sheet increases with the resistivity (Loureiro et al. 2007), which reduces the growth rate of the tearing mode when it is normalized with respect to the Alfvén crossing time along the length of the current sheet (in the x-direction in our case). Hence, the explosive phase of the evolution in our simulation starts at later times for higher resistivity values. We estimated the absolute current density, |J_z| (normalized to unity) before the fragmentation stage of the current sheet (t = 214 s) by taking a vertical cut along the y-direction at x = 0 for two different resistivities, η = 0.0001 and 0.001, to confirm that the thickness of the current sheet increases with resistivity (see Fig. 14). The resistivity dependence for the different evolution phases is shown in Fig. 15. Figure 15a shows that γ_Ruth follows a power-law dependence with the resistivity, γ_Ruth ≈ η^−0.1, with a correlation coefficient (CC) of −64.1%. The resistivity scaling relation for the post-Rutherford and the entire nonlinear regimes are shown in Figs. 15b and c respectively. We estimate the growth rate scaling relations for the post-Rutherford regime, γ_PR ≈ η^0.03 (with CC = 59.9%), and the entire non-linear regime, γ_avg ≈ η^0.017 (with CC = 66.7%). Previous studies by Zhang & Ma (2011), Akramov & Baty (2017), and Guo et al. (2017 and references therein) have reported the resistivity scaling relation of the non-linear growth rates for the DTM setup in the adiabatic environment, which have larger power-law indices compared to our estimation. Hence, our study infers that the resistivity dependence on the nonlinear growth rates is weaker when the thermal instability reinforces the tearing mode.

Fig. 12. Distribution of the velocity, vx, within the same domain as Fig. 2, for different evolution stages. The magnitude of the velocity is scaled with respect to the Alfvén velocity, vA. An animation of the figures is available online.

Fig. 13. Maximum absolute value of vx as a function of time. The left panel represents the velocity evolution of the current layer system for two different cases: adiabatic medium, and non-adiabatic medium when radiative energy loss and background heating are incorporated for resistivity, η = 0.001 and plasma-β = 0.2. The right panel is for non-adiabatic evolution for different η values keeping all the other parameters the same. The vertical dashed lines represent the different phases of the evolution.

Fig. 14. Variation in the absolute current density |Jz| (normalized to unity) at t = 214.68 s, at x = 0 along the y-direction between ±3.75 × 103 km, for two different resistivities, η = 0.0001 and 0.001.

Fig. 15. Variation in the average growth rates, γ with the resistivities for (a) Rutherford, (b) Post-Rutherford, and (c) entire nonlinear regions. The growth rates are scaled with respect to the Alfvén timescale. The solid lines represent the linear fit of the growth rates vs. η in the log-log scale. The values of the slope and the correlation coefficients are noted in the corresponding figures.

3.3. Plasmoid distribution and density filling factor

Due to the combined thermal and tearing instability, the current sheet becomes unstable and the magnetic islands are formed as shown in Figs. 2c and d. These islands show the coalescence tendency to merge with the neighborhood companions to form larger plasmoids (see Figs. 4c and d). The mass density of the plasmoids is higher than in the local background medium. Since we wish to quantify some statistical properties on the evolving and coalescing plasmoids and their internal thermodynamics, we need a criterion to identify and count them. We define the plasmoids by a density threshold condition: if the density of a region is higher than a density threshold, ρ_th, then we call it a plasmoid. Here the threshold density ρ_th is defined as 0.03% of the peak density (ρ_max) for an instantaneous time. This means the ρ_th varies with time, according to the temporal variation in ρ_max. This is equivalent to capturing the density regions up to 3σ level (99.97%) of the peak density of a Gaussian distribution and separating out the density enhanced regions from their local background medium. We also calculate the fraction of volume occupied by these density enhanced structures (or plasmoids) with respect to the entire volume of the simulation domain (which is 256 × 256 Mm²), which we call the density filling factor. The number distribution of the plasmoids

and the density filling factor with time for different resistivities are shown in Fig. 16. At the initial time (t = 0), the ratio of ρ_min to ρ_max is 16.67% (which is more than 0.03%). Hence, the threshold density captures the entire volume of the simulation domain and therefore the plasmoid numbers obtained in this method remain at unity unless the density contrast satisfies the density threshold condition as mentioned above (see Fig. 16a). After this phase when the current sheet breaks up, the maximum number of small-scale plasmoids form, which is represented by the highest peak of the distribution curve of Fig. 16a. We also see in Fig. 16a that the current sheet breaks up into more plasmoids for lower resistivity values, and the fragmentation phase occurs at a later time for higher resistivity. We obtain an inverse relation of the maximum plasmoid numbers with resistivity that follows a scaling relation, N ∼ η^0.223 (see Fig. 17). This is similar to the case N ∼ η^3/8 obtained by Loureiro et al. (2007) for an adiabatic medium, where it is reported that the current sheet thickness depends on the resistivity. Thus, our simulation results also imply that the current sheet thickness is a function of the resistivity similar to the Sweet–Parker-type current sheet, which is also reflected in Fig. 14. Similarly, the density filling factor for the plasmoids at the initial phase is unity, as shown in Fig. 16b, and it decreases with time, so the volume fraction occupied by the plasmoids diminishes at the later stage of the evolution. The distribution of the maximum plasmoid numbers with the Lundquist number, S_L, is shown in Fig. 17, which shows that the number scales with $S_{\mathrm{L}}^{0.223}$ (with a correlation coefficient of the linear fit equal to 0.96). The stability analysis by Loureiro et al. (2007) reports that this number distribution scales as $S_{\mathrm{L}}^{0.375}$ in an adiabatic medium. This suggests that the dependence of the Lundquist number in the number distribution of secondary islands is weaker for a non-adiabatic medium. In the Sweet–Parker-type current sheet the inverse aspect ratio of the current sheet, l_s/L, follows the scaling as $l_s/L\sim S_{\mathrm{L}}^{-\alpha }$, where α = 0.5. In our study, for l_s = 0.5, L = 12.8, and S_L in the range from 4.67 × 10³ to 2.34 × 10⁵, we estimate the value of α that varies between 0.26 and 0.38, which is clearly lower than 1/2. This implies that the current sheet is thicker than the Sweet–Parker-type within our explored domain of the Lundquist number. We also estimate the temporal variation for the distribution of plasmoid numbers and the density filling factor for different plasma-β values, which are shown in Figs. 19a and b respectively. To vary the plasma-β = 0.2, 0.1, and 0.02, we tune the initial temperature to T_i = 0.5, 0.25, and 0.05 MK, respectively, keeping the magnetic field strength constant at B₀ = 2 G. From Fig. 19a we see that the maximum fragmentation phase of the current sheet occurs at a later time for higher plasma-β values, but there is no specific trend for the peak values of the plasmoids with plasma-β. We perform the same analysis by upgrading the maximum numerical resolution by 4096 × 4096 with η = 0.001 and plasma-β = 0.2, keeping all the other parameters constant (see Fig. 19). The numerical resistivity of the medium decreases for higher resolution values, and hence the maximum plasmoid number in the fragmentation phase increases, as shown in Fig. 19a, though the results do not change significantly for different resolutions. The result is also consistent with the study of the physical resistivity cases, as shown in Fig. 16a. We also find that the variation in the density filling factor distribution with time is not significantly different for different plasma-β (see Fig. 19b) and numerical resolutions (see Fig. 19b).

Fig. 16. Temporal variation in (a) plasmoid numbers and (b) density filling factor for different η with plasma-β = 0.2.

Fig. 17. Maximum number of the plasmoids vs. Lundquist number, SL, obtained from different simulation runs (for plasma-β = 0.2) as shown by the blue triangles. The solid line represents the linear fit in the log-log scale of the plasmoid numbers vs. SL. The slope and Pearson’s correlation coefficient for the linear fit are 0.223 and 0.96, respectively.

Fig. 18. Temporal variation in (a) plasmoid numbers and (b) density filling factor for different plasma-β with η = 0.001.

Fig. 19. Temporal variation in (a) plasmoid numbers, and (b) density filling factor for different numerical resolution for η = 0.001 and plasma-β = 0.2.

Finally, we estimate the kinetic energy density

$$\begin{array}{c}\hfill \mathrm{KE}=\frac{1}{V}{\displaystyle \int \int \frac{\rho v^2}{2}}\mathrm{d}x\mathrm{d}y,\end{array}$$(15)

for each time step, where we integrate over of the full simulation domain, V = l_xl_y, and also quantify the maximum Ohmic heating, $H_{\mathrm{o}\mathrm{h}\mathrm{m}}=\eta |J_z{|}_{\mathrm{m}\mathrm{a}\mathrm{x}}^2$. The evolution of the kinetic energy density and the maximum Ohmic heating rate for different plasma-β are shown in Figs. 20a and b respectively, which shows that the maximum energy dissipation per unit time by the Ohmic heating is approximately two orders of magnitude less than the kinetic energy. There is more rapid temporal variation in the later merging stages of the plasmoids as seen in the Ohmic heating extremal evolution, and only some modest dependency of the overall energetics on the plasma beta parameter.

Fig. 20. Temporal variation in (a) kinetic energy density and (b) Ohmic heating rate for η = 0.001 and different plasma-β.

4. Discussion and summary

The instability problem addressed in this work can be related to the pre-flare current layer model. A cartoon geometry of the configuration for a current sheet associated with a pre-flare event is shown in Fig. 5 of Ledentsov (2021a), where the theoretical linear stability analysis is carried out with the inclusion of viscosity, electrical and thermal conductivity, and radiative cooling. In the follow-up works by Ledentsov (2021b,c), the effects of the guiding magnetic field and the oblique fragmentation of the current sheet are investigated in linear MHD with an analysis of the growth rate and spatial periodicity scales of the instability. In the present work, instead, we extend the analysis for growth rate for both linear and nonlinear domains with the incorporation of radiative losses, and constant background heating by a series of resistive nonlinear high resolution MHD simulations. The instability that occurs in this work can be seen as a thermal instability enhanced pathway to rapid small-scale tearing behavior, where thereafter coalescing islands evolve and collect small condensations into larger entrapped cool plasma sites within magnetic islands. At the initial time (t = 0) the thermal balance is maintained due to the equal and counteracting effects of the radiative loss and

background heating. However, due to the magnetic field perturbation and the finite resistivity, some of the sections within the current sheet begin to increase in density, and hence lose more heat due to radiation as compared to the constant background heat. Hence, a thermal imbalance occurs, which leads to instability. As a result, the current sheet starts to disintegrate into the form of plasmoids, which move along the current sheet by merging with the neighboring plasmoids. The regions of main solar flare energy release occur from current sheet regions, which can lead to the outburst of a coronal mass ejection and modern observations may detect fine-scale multi-thermal structure in the reconnecting sheet, by the brightening in the ultraviolet (UV) regime (Jing et al. 2016; Wyper & Pontin 2021). As shown in Fig. 4, the secondary islands are present near the magnetic X-points in the form of plasmoids that appear during the explosive phase. The local Lundquist number defined by S_L = l_xv_A/η (where l_x is the length of the current layer), in the series of our simulations is in the range of 4.67 × 10³–2.34 × 10⁵. We see that the formation of the plasmoids occurs over the entire explored Lundquist number range in our simulation. We note that this extends the chaotic reconnection process to a much lower Lundquist range than previously found from purely resistive (but otherwise adiabatic) evolutions. In the previous studies by Loureiro et al. (2007), Samtaney et al. (2009), Bhattacharjee et al. (2009), and Akramov & Baty (2017), it is reported that the plasmoid formation occurs only beyond a minimum threshold value of the Lundquist number, S_L ∼ 10⁴, where those studies were limited for adiabatic regimes. This threshold is also determined by the inhomogeneous inflow and outflow in the reconnecting sheet, which suppresses the growth of tearing in an adiabatic medium (Shi et al. 2018). This means that the formation of plasmoids may occur for lower Lundquist numbers (≲10⁴), due to the non-adiabatic effects of thermal instability, and this lower Lundquist range is easily resolvable numerically. Even so, extreme resolution is warranted because the fine-scale effects this time are due to the thermally unstable nature of the solar coronal plasma. We also investigate an experiment by switching off the physical resistivity (i.e., η = 0), but keeping the adiabatic effects on for a sufficiently high numerical resolution of 2048 × 2048, so that this experiment only evolves due to numerical (unresolved) resistivity. This experiment shows that even for a low (numerical) resistivity value, the formation of plasmoids occurs due to the fragmentation of the current layer, and we get explosive behavior.

We highlight the following novel features and the key results of this work.

1. We set up a numerical experiment of a 2D current sheet model by incorporating non-adiabatic effects of radiative loss and constant background heat in a resistive MHD simulation using MPI-AMRVAC. Due to the magnetic field perturbations (Eqs. (12) and (13)), the equilibrium of the system breaks down, and the instability begins, in the form of thermal and tearing modes. The current layer starts to disintegrate to form secondary islands, which move along the current sheet by merging with the neighboring plasmoids.

2. The thermodynamical behavior of the current sheet region, and for a local plasmoid in particular is analyzed. We see a density enhancement, due to the accumulation of the neighboring plasmoids, or a density drop, due to the disintegration of the current layer. Accordingly, the temperature of the medium drops for the density enhanced regions and increases for the density depleted regions. We also estimate the global behavior of the energy loss in the medium due to the optically thin radiation.

3. We compute the growth rates for the linear and nonlinear phases of the evolution and estimate the scaling relations with the resistivities for different nonlinear phases. We find that the growth rate obtained from our model is faster by an order of magnitude, and weaker with resistivity in comparison with earlier works by Otto & Birk (1992), Zhang & Ma (2011), Akramov & Baty (2017), and Paul & Vaidya (2021), who assume the medium is adiabatic. We find the occurrence of the explosive nature of the evolution within the resistivity domain of η = 5 × 10⁻³ to 10⁻⁴ in our work. This is a regime where thermal instability enhanced fragmentation triggers small-scale tearing effects.

4. The temporal variation in the generated plasmoid numbers and the associated density filling factors are estimated for different η and plasma-β that are relevant for the solar corona. We calculate the scaling relation of the maximum plasmoid numbers with the Lundquist number, S_L, and find it to vary as $S_{\mathrm{L}}^{0.223}$, which is smaller than the value estimated by Loureiro et al. (2007) for an adiabatic medium. This implies a thicker tearing-unstable current sheet than the usual Sweet–Parker-type, and indicates that the thermal instability facilitates the triggering of tearing modes. We also investigate the analysis for higher numerical resolution (4096 × 4096) and see that it does not alter the results significantly.

5. The time evolution of the kinetic energy density and Ohmic dissipation rate are calculated. The comparison between these two energies shows that the energy dissipation per unit time due to Ohmic heating is around two orders of magnitude less than the kinetic energy. The later stages show clear Petschek-like super-Alfvénic outflow regions between the merged larger islands.

We performed a detailed magnetohydrodynamic study of a current sheet model liable to both thermal and tearing instabilities. We did not incorporate the effects due to thermal conduction, a guide magnetic field, or gravity in this model. A more realistic 3D model with the incorporation of these effects can be explored in near future. However, the current idealized study sheds new light on the formation mechanisms of plasmoids, and explosive reconnection behavior of a pre-flare current layer model, which is one of the important aspects of solar coronal heating. Our findings suggest that multi-thermal plasma aspects must be common in flaring regions, where both hot islands with their entrapped cooler condensations must show up cospatially in different wavebands.

Movies

Movie 1 associated with Fig. 2 (fig2_anim) Access here

Movie 2 associated with Fig. 4 (fig4_anim) Access here

Movie 3 associated with Fig. 5 (fig5_anim) Access here

Movie 4 associated with Fig. 12 (fig12_anim) Access here

Acknowledgments

Data visualization and analysis are performed using ParaView (https://www.paraview.org) and python (https://www.python.org/). S.S. and R.K. acknowledge support by the C1 project TRACESpace funded by KU Leuven. R.K. acknowledges the support by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 833251 PROMINENT ERC-ADG 2018) and a FWO project G0B4521N. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government - department EWI. We thank N. Yadav, X. Li, and J. Hermans for useful discussions during the course of this work.

References

All Figures

	Fig. 1. Radiative cooling curve due to the Colgan_DM model.
In the text

	Fig. 2. Distribution of the absolute current density, |Jz|, normalized with the instantaneous absolute peak current density, |Jmax|, for different evolution stages. The y-domain is only shown in the range −1.45 × 104 km ≤ y ≤ +1.45 × 104 km, which contains the region of the current sheet. An animation of the figures is available online.
In the text

	Fig. 3. Variation in ${\displaystyle F_{\mathit{y}}=-\frac{\partial p}{\partial \mathit{y}}+B_x(\frac{\partial B_{\mathit{y}}}{\partial x}-\frac{\partial B_x}{\partial \mathit{y}})}$ along the y-direction in the range y = ±0.75 × 104 km for three different times (t = 0, 300, and 429 s) before the fragmentation stage of the current sheet.
In the text

	Fig. 4. Same as Fig. 2, but for plasma density, ρ, for different evolution stages. The overplotted white lines represent the magnetic field lines. An animation of the figures is available online.
In the text

	Fig. 5. Distribution of the temperature, T, within the same domain as Fig. 2, for different evolution stages. An animation of the figures is available online.
In the text

	Fig. 6. Distribution of (a) radiative cooling and background heating, Hbgr, and (b) temperature, T, along the vertical cut at x = 0 in the range y = ±1.25 × 104 km for t = 0 and 429 s. The minimum temperature obtained within this regime at t = 429 s is Tmin = 21 000 K.
In the text

	Fig. 7. Spatial distribution of plasma density along the y = 0 cut for two different evolution stages.
In the text

	Fig. 8. Temporal evolution of the instantaneous maximum and minimum densities (left) and temperatures (right) within the entire simulation domain.
In the text

	Fig. 9. Analysis of density distribution of local plasmoids. Top panel: (a) Distribution of the density for t = 900 s, within the same domain as Fig. 2; (b) zoomed-in version of the selected region. Bottom panel: Same as top panel, but for t = 2447 s.
In the text

	Fig. 10. Analysis of temperatutre distribution of local plasmoids. Top panel: (a) Distribution of the temperature for t = 900 s, within the same domain as Fig. 2; (b) zoomed-in version of the selected region; (c) Density and temperature distributions along the horizontal cuts indicated by the dashed lines in Figs. 9b and 10b, respectively. Bottom panel: same as top panel, but for t = 2447 s.
In the text

	Fig. 11. Same as Fig. 10 but for radiative loss in an optically thin medium.
In the text

	Fig. 12. Distribution of the velocity, vx, within the same domain as Fig. 2, for different evolution stages. The magnitude of the velocity is scaled with respect to the Alfvén velocity, vA. An animation of the figures is available online.
In the text

Fig. 13. Maximum absolute value of vx as a function of time. The left panel represents the velocity evolution of the current layer system for two different cases: adiabatic medium, and non-adiabatic medium when radiative energy loss and background heating are incorporated for resistivity, η = 0.001 and plasma-β = 0.2. The right panel is for non-adiabatic evolution for different η values keeping all the other parameters the same. The vertical dashed lines represent the different phases of the evolution.

In the text

	Fig. 14. Variation in the absolute current density |Jz| (normalized to unity) at t = 214.68 s, at x = 0 along the y-direction between ±3.75 × 103 km, for two different resistivities, η = 0.0001 and 0.001.
In the text

	Fig. 15. Variation in the average growth rates, γ with the resistivities for (a) Rutherford, (b) Post-Rutherford, and (c) entire nonlinear regions. The growth rates are scaled with respect to the Alfvén timescale. The solid lines represent the linear fit of the growth rates vs. η in the log-log scale. The values of the slope and the correlation coefficients are noted in the corresponding figures.
In the text

	Fig. 16. Temporal variation in (a) plasmoid numbers and (b) density filling factor for different η with plasma-β = 0.2.
In the text

	Fig. 17. Maximum number of the plasmoids vs. Lundquist number, SL, obtained from different simulation runs (for plasma-β = 0.2) as shown by the blue triangles. The solid line represents the linear fit in the log-log scale of the plasmoid numbers vs. SL. The slope and Pearson’s correlation coefficient for the linear fit are 0.223 and 0.96, respectively.
In the text

	Fig. 18. Temporal variation in (a) plasmoid numbers and (b) density filling factor for different plasma-β with η = 0.001.
In the text

	Fig. 19. Temporal variation in (a) plasmoid numbers, and (b) density filling factor for different numerical resolution for η = 0.001 and plasma-β = 0.2.
In the text

	Fig. 20. Temporal variation in (a) kinetic energy density and (b) Ohmic heating rate for η = 0.001 and different plasma-β.
In the text