Numerical MHD simulations of solar flares and their associated small-scale structures

Using numerical simulations, we study the formation and dynamics of solar flares in a local region of the solar atmosphere. The magnetohydrodynamics (MHD) equations describe the dynamic evolution of flares, including space-dependent and anomalous magnetic resistivity and highly anisotropic thermal conduction on a 2.5 D slice. We adopt an initial solar atmospheric model in magnetohydrostatic equilibrium, with a magnetic configuration consisting of a vertical current sheet, which helps trigger the magnetic reconnection process. Specifically, we study three scenarios, two with only resistivity and the third with resistivity plus thermal conduction. The main results of the numerical simulations show differences in the global morphology of the flares, including the post-flare loops and the current sheet in three cases. In particular, localized resistivity produces more substructure around the post-flare loops that could be related to the Ritchmyer-Meshkov Instability (RMI). Furthermore, in the scenario of anomalous resistivity, we identify the formation of a plasmoid and a jet at coronal heights. On the other hand, in the scenario with resistivity plus thermal conduction, the post-flare loops are smooth, and no apparent substructures develop. Besides, in the 𝑧 − component of the current density for the Res+TC case, we observe the development of multiple magnetic islands generated due to the Tearing instability in the non-linear regime.


INTRODUCTION
Solar flares are the most energetic explosions in the heliosphere, and they power space weather events that can affect the magnetic environment of our Earth.Solar flares efficiently convert magnetic energy into internal and kinetic energy, producing emissions in different wavelengths.One of the standard models of the solar flares has been discussed in Tajima & Shibata 2002;Forbes 2010;Priest 2014a, where the authors extended the cartoon view of magnetic reconnection using observational facts and, therefore, identify processes such as fast evaporation upflows, X-ray loops, and hard X-ray (HRX) source regions (Masuda et al. 1994;Su et al. 2013).One of the most relevant 2D MHD simulations of a solar flare based on a magnetic reconnection model with anisotropic thermal conduction was performed by Yokoyama & Shibata (2001), who studied the thermal evolution in the post-flare loops to determine the flare temperature.Similar MHD simulations of magnetic reconnection have been performed to study coronal mass ejections and associated giant arcades (Shiota et al. 2005), and chromospheric evaporation jets (Miyagoshi & Yokoyama 2004).More recently, Ruan et al. (2020) presented a self-consistent solar MHD flare model that demonstrated the observationally suggested relationship between flux swept out by the hard X-ray footpoint regions and the actual reconnection rate at the Xpoint, which is a significant unknown in flaring scenarios.In solar flares, various plasma flows above the post-reconnection flare arcade ★ E-mail: jgonzaleza@enesmorelia.unam.mxhave been reported by Liu et al. (2013); Lin et al. (2005); Warren et al. (2018).In particular, flare current sheets have rapidly descended, dark finger-like features commonly referred to as supra-arcade downflows (SADs), often observed in extreme ultraviolet (EUV) soft X-ray images.The slow speeds of SADs have been discussed in numerical models containing the bursty jet (Cécere et al. 2015) or Rayleigh-Taylor-type (RT) instabilities in reconnection downflow regions (Guo et al. 2014).Other features observed in post-flare loops are the formation of magnetic islands at the top of the loops and along the current sheet (CS).These magnetic islands are believed to be closely related to fine structures CS in the wake of the eruption (Shen et al. 2011;Lin et al. 2005).Moreover, the magnetic islands' coalescence can further produce secondary CSs and islands (Bárta et al. 2011).More recently, Wang et al. (2021) performed 2.5D MHD simulations to study the dynamics of magnetic reconnection in the solar CS, and they suggest that annihilation of magnetic islands at the flare loop top, which is not included in the standard flare model, plays a nonnegligible role in releasing magnetic energy to heat flare plasma and accelerate particles.
In this paper, we study the formation and evolution of a solar flare, including the post-flare loops, the CS, and the development of associated substructures such as shocks, MHD instabilities, and small-scale magnetic islands, including the effects of localized and anomalous resistivity, and highly anisotropic thermal conduction on the formation and dynamics of a solar flare following the main ideas of the classical model developed by Yokoyama & Shibata (2001).In particular, we base our simulations on three cases: i) Res case, where we only consider the effects of a localized resistivity; ii) anomalous resistivity case; and iii) Res+TC, where we include the effect of the highly anisotropic thermal conduction additional to the localized resistivity.We briefly describe the results of case ii) corresponding to the anomalous resistivity profile.
We organize the paper as follows.First, in section 2, we describe the system of MHD equations, the solar atmospheric model, the magnetic field configuration, the resistivity and thermal conduction models, and the numerical methods used to solve the MHD equations.Next, section 3 contains the results of the two study cases, including the analysis of MHD instabilities, shock development, and magnetic island formation.Finally, we draw our concluding remarks in section 4.
We organize the paper as follows.First, in section 2, we describe the system of MHD equations, the solar atmospheric model, the magnetic field configuration, the resistivity and thermal conduction models, and the numerical methods used to solve the MHD equations.Next, section 3 contains the results of the two study cases, including the analysis of MHD instabilities, shock development, and magnetic island formation.Finally, we draw our concluding remarks in section 4.

The system of MHD equations
We consider a gravitationally stratified solar atmosphere described by a plasma obeying the resistive MHD equations with anisotropic thermal conduction to model the plasma that constitutes the solar flare.In particular, we write the equations in a conservative dimensionless form that helps apply the numerical methods easily: (1) where  is the mass density, v represents the velocity field, B is the magnetic field,  is the temperature,  is the magnetic resistivity, and F c is the thermal conduction flux.
Besides,   =  + B 2 /2 is the total (thermal + magnetic) pressure, I is the unit matrix, and  is the total energy density, that is, the sum of the internal, kinetic, and magnetic energy densities, For the fluid, we consider the adiabatic index  = 5/3.The system of equations ( 1)-( 5) is closed with the ideal gas law where  is the temperature of the plasma, m =   is the particle mass specified by a mean molecular weight value  =0.6, for a fully ionized gas, and   is hydrogen's mass, and   is Boltzmann's constant.The gravitational source term on the right-hand side of equations ( 2) and ( 3) is given by g = [0, −] with magnitude  = 2.74 × 10 4 cm s −2 , which represents an average over the solar surface.In this paper, we use a 2.5D model, meaning all the state variables depend on  and , where  is a horizontal coordinate, and  represents height.Besides, a nontrivial magnetic field −component   varies with  is allowed.We solve the 2.5D version of the resistive MHD equations with thermal conduction on the  plane with resolution Δ and Δ.

Model of the solar atmosphere
At the initial time of the simulations, we assume the solar atmosphere is in hydrostatic equilibrium.In particular, we choose the C7 semiempirical model to describe the temperature field of the chromosphere-transition region (Avrett & Loeser 2008), smoothly extended to the solar corona as in González-Avilés et al. (2017).In Fig. 1, we display the equilibrium temperature and mass density profiles as functions of height , where the steep gradient at the transition region  ∼ 2.1 Mm is discernible.

The magnetic field
Initially, we consider a force-free magnetic field ((∇ × B) × B = 0) similar to that in Yokoyama & Shibata (2001); Ruan et al. (2020), which essentially represents a vertical CS defined as follows: =  0 tanh(/), ( 9) here  0 = 22.42 G is the magnetic field strength and  = 5 × 10 7 cm is the width of the CS.We show the 2D maps of the   and   components in Fig. 2, where it is notable that both components only vary in  and are constant along the vertical direction .However, initially, there is a non-negligible − component of the current density, which produces a CS that helps trigger the magnetic reconnection at late times in the simulation.Additionally, we plot the plasma , which is a parameter that estimates the ratio between gas pressure and magnetic pressure, i.e.,

𝛽(𝑥
Here, the pressure () is given by the hydrostatic equilibrium, and  2 = ( 2  +  2  ).We display the plasma  at the initial time in the right of Fig. 2, where it is evident that  is about 10 in the bottom of the photosphere  = 0 Mm, while in the solar corona ( > 2.1 Mm) is approximately 10 −3 .This parameter will help us explain the dominance of pressures in regions where the magnetic islands form, as described in Section 3.Moreover, the solar atmospheric model, which is in hydrostatic equilibrium, and the magnetic field configuration, which is force-free, satisfying the magnetohydrostatic equation, J × B − ∇ −  = 0. Therefore, the whole system achieves the magnetohydrostatic equilibrium at the initial time of the simulations.

Spatially dependent and anomalous resistivity profiles
To perturb the atmosphere in magnetohydrostatic equilibrium and trigger the magnetic reconnection, we set a spatial dependent function for the resistivity, (, ).In particular, we use a similar localized where  0 = 10 16 cm 2 s −1 is the strength of the resistivity,  = √︃  2 − ( − ℎ  ) 2 , with ℎ  = 6 × 10 8 cm (6 Mm) is the vertical location of the Gaussian profile, and   = 1 × 10 8 cm (1 Mm) is its width.This profile is time-independent, but it can induce a fast and quasi-steady magnetic reconnection with a single X-point (see, e.g., Ugai 1992).
We also include an anomalous resistivity function.Microscopic instabilities physically cause this resistivity.Specifically, we adopt the following profile: where,   = 1×10 −4 ,   =0.1, and   = 10 s.In the expression above,   (, ) = J/(  ) represents an instantaneous distribution that quantifies the relative ion-electron drift velocity, where J is the current density magnitude,  is the electron charge and   is the electron number density.Profile of equation 13 is similar to the used in Ruan et al. (2020); however, we do not couple with the localized profile of equation 12 and consider a small threshold drift velocity   .

Anisotropic thermal conduction
In the equation (3), the thermal conduction flux F c defines a fluxlimited expression that smoothly varies between the classical and saturated thermal conduction regimes   and   , respectively, that is: The above equation represents highly anisotropic thermal conduction suppressed mainly in the direction transverse to the magnetic field lines.For example, in this paper, we denote b = B/|B| as the unit vector in the direction of the magnetic field; therefore, the classical thermal conduction flux is as follows: where the subscripts ∥ and ⊥ denote, respectively, the parallel and normal components to the magnetic field,  is the temperature,  ∥ and  ⊥ are the thermal conduction coefficients along and across the field.In most astrophysical scenarios, such as in the Sun,  ⊥ / ∥ = 2 × 10 −31  2  /( 3 |B| 2 ) is negligible (Priest 2014b), and the thermal conduction is mainly along the field lines.In this paper, we chose  ∥ = 8 × 10 −7 erg s −1 cm −1 K −1 , which is a typical value in the solar corona (Ruan et al. 2020).Saturated effects are accounted for, making the flux independent of ∇ for substantial temperature gradients, such as the transition region to solar corona in the solar atmosphere, which is a novel ingredient of this paper compared to the similar works where include the thermal conduction.In this limit, the flux magnitude approaches   = 5 3  , where   is the isothermal speed of sound and  < 1 is a free parameter that helps to account for the uncertain numerical coefficient in   (Giuliani 1984), which we take equal to 0.3 in this paper.The equations 1-5 were normalized by the scale factors in CGS units as shown in table 1, which are typical scales of the solar atmosphere.The CGS units are used throughout the manuscript to be consistent with the normalization factors.

Numerical methods
We numerically solve the equations ( 1)-( 5) using the PLUTO code (Mignone et al. 2007).In all simulations, we set the Courant-Friedrichs-Levy (CFL) number equal to 0.4 and choose the secondorder total variation diminishing (TVD) Runge-Kutta time integrator.Additionally, we use the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano 2005) in combination with a linear reconstructor and the minmod limiter.
The numerical evolution of the MHD equations can lead to the violation of the divergence-free constraint equation given by ( 5), developing unphysical results like the presence of a net magnetic charge.We use the extended generalized Lagrange multiplier method (Dedner et al. 2002) to control the constraint violation's growth.This method is robust in highly magnetized regions or low plasma beta scenarios, as implied in the solar corona on the left panel of Fig. 2. PLUTO evolves the magnetic resistivity and thermal conduction separately from advection through operator splitting, using the super-time-stepping technique (Alexiades et al. 1996).This method is robust enough to handle the highly parabolic (diffusion) behavior of space-dependent resistivity and anisotropic thermal conduction in the MHD equations adopted in this paper.Our simulations are carried out in a domain  ∈ [−10, 10],  ∈ [0, 20], in units of Mm (10 8 cm), covered by 800×1000 grid cells.Here,  = 0 Mm represents the bottom of the photosphere.We used the following boundary conditions: outflow boundary conditions at   = −10 Mm,   = 10 Mm, and   = 20 Mm.In contrast, at   = 0 Mm, we employ asymmetric boundary conditions in the −component of the magnetic field   , outflow in the −component of the magnetic field   , mass density and pressure, and zero boundary conditions in the component of the velocity field.These boundary conditions are similar to those used in Yokoyama & Shibata (2001).

RESULTS OF NUMERICAL SIMULATIONS
In this section, we show the most representative results of the numerical simulations performed for the following three scenarios: (i) resistivity case (Res), ii) anomalous resistivity case, and (iii) resistivity plus thermal conduction case (Res+TC).We only include the most representative results for the case of the anomalous resistivity profile.

Resistivity case
In Fig. 3, we show the temperature in Kelvin and mass density in gr cm −3 with the magnetic field lines for the Res case.For example, in panel (a), we show the temperature at  = 25 s.There, we note the following plasma structures that are consistent with the standard flare model (Yokoyama & Shibata 2001): i) formation of a hot flare loop (∼ 10 8 K) due to the magnetic reconnection triggered by the localized resistivity, ii) the reconnection outflow collides with the flare loop, and iii) production of evaporation flows of about 10 6 K in the flare loop.In panel (c), the temperature at  = 50 s shows the formation of a small hot region of about 10 6 K at the top of the loops, and it is also evident that inside the loop region, the chromospheric evaporation develops.In panel (e), we show the mass density at  = 25 s, consistently proving the behavior demonstrated on temperature maps.In particular, chromospheric evaporation appears in dense regions of about 10 −14 gr cm −3 .Finally, in panel (g), the mass density at  = 50 s shows how dense material (∼ 10 −12 ) moves upwards while it forms dense region (∼ 10 −17 gr cm −3 ) around it, and it is also visible the small region (but hot) at the top of it, as shown by the temperature maps.
In Fig. 4, we show snapshots of the temperature in kelvin and mass density in gr cm −3 at times  = 25, 50 s for the anomalous resistivity case.In this scenario, we only turn on the resistivity given by the equation ( 13) and let the system evolve.In panels (a) and (b), we display the temperature, where it is discernible that any postflare structure is formed, especially since we do not identify the post-flare loop structures and the development of a thin CS along the vertical direction.Instead, we can see the formation of a plasma structure that becomes wider and moves upwards while the magnetic file lines reconnect.This kind of plasma feature could be related to a plasmoid.In panels (c) and (d), we note that in the mass density, there are no clues of the formation of post-flare loops, but interestingly, we identify the development of a small jet (3 Mm in length) moving upwards from the transition region.This feature is interesting since it could be related to the magnetic reconnection process at microscopic scales, which is triggered by anomalous resistivity at coronal heights, which indicates a relevant result despite the limitations of the MHD to capture the microscopic mechanism associated with the magnetic reconnection process.Finally, it is discernible that these plots are different from the plots of panels (a), (c), (e), and (g) shown in Fig. 3 for the Res case.
In Fig. 5, we show the −component of current density   in statA cm −2 and plasma pressure in dyn cm −2 with the magnetic field lines overlap for the two simulation cases.Specifically, in panels (a), (c), (e), and (g), we display the results for the Res case at times  = 25 s and  = 50 s.For instance, in panel (a) at  = 25 s, we identify substructures in   .In particular, there are regions with high current density inside the post-flare loops, and a thin region of high   is visible that is part of the CS.While, in panel (c), at  = 50 s, we observe substructures at about  ≈ 12 Mm and  ≈ 17 Mm.In panel (e), at  = 25 s, we identify regions of high plasma pressure (∼ 18 dyn cm −2 ) inside the post-flare loops, and  = 50 s (g), in the plasma pressure is also evident the structure shown in   of panel (c), in this particular region, the plasma pressure is around 1.5 dyn cm −2 .

Resistivity with thermal conduction case
In the right panels of Fig. 3, we display the results for the Res+TC case.For example, in panel (b), we identify the formation of a typical flare structure and its associated post-flare loops; however, in this case, the thermal conduction redistributes the plasma inside the flare loops, which makes the evaporation flow less visible than in the Res case.Furthermore, the maximum temperature reaches ∼ 10 7 K, smaller than the maximum temperature in the Res case, as shown in panel (a).In panel (d), we do not observe any structured region on the top of loops, nor the increase of evaporation of the plasma inside it, as shown in panel (c).This feature could be related to the effect of thermal conduction along the magnetic field lines, which   In panels (f) and (h), the plasma pressure shows the same smooth behavior as   at both times,  = 25 and  = 50 s.Notably, the regions with high plasma pressure (∼ 18 dyn cm −2 ) are mainly around the post-flare loops.
To show how the thermal conduction behaves on the flare and post-flare loop structures, we estimate the magnitude of the heat flux q.In Fig. 6  shown in Fig. 3.In particular, we observe a temperature growth at about  ≈ 3 Mm, in both cases at the two times.Nevertheless, in the Res case, the plasma temperature reaches higher values (∼ 10 9 K).In contrast, in the Res+TC case, the temperature of the plasma maintains an order of magnitude of about 10 7 K.We can interpret these results in terms of the effects that thermal conduction has over the plasma temperature, i.e., in the Res case, the heat continues to affect the fluid due to high dynamics.In contrast, in the Res+TC case, the heat is dissipated and redistributed, so the temperature does not rise considerably.

Instabilities
In general, in a fluid the hydrodynamic instabilities such as Rayleigh-Taylor (RT) and Ritchmyer-Meshkov (RM) frequently appear when turbulence takes place.In particular, the RT instability occurs at an interface between different fluids when the lighter fluid is accelerated into the heavy.While the RM instability may be considered a particular case of RT instability, it develops when the acceleration provided is impulsive, resulting from a shock wave (Zhou et al. 2021).Furthermore, in Fig. 8, we show a zoom region of the temperature and vertical velocity overlap with the velocity vector field at time  = 50 s, for the Res case (panels (a) and (c)) and for the Res+TC case (panels (b) and (d).In panels (a) and (c) of the figure, we identify plasma flowing downward and upward, shown in the velocity vector field; therefore, these flows could generate the 'spike' structures inside the post-flare loop regions.In contrast, in panels (b) and (d), we see that plasma moves predominantly upwards, and the post-flare loop looks homogeneous; this is consistent with the behavior of the heat flow that is high inside the flare loops, as schematized by the maps of |q| in Fig. 6.In the four figures is an orange vertical dashed line, which represents the distance where the pressure, density, and velocity gradients are estimated for the instabilities analysis shown below.
To investigate the possible development of RTI/RMI, we estimate the variations of the pressure gradient, density gradient, vertical velocity gradient, and plasma  along the orange vertical dashed line drawn in the panels of Fig. 8, that is similar to the analysis performed by Shen et al. (2022).Notably, we estimate how the above-mentioned variable changes when crossing the interface around the transition region ( =∼2.1 Mm).In Fig. 9, we show 1D cuts of pressure gradient in dyn cm −3 , density gradient in gr cm −4 , vertical velocity (  ) gradient in s −1 and plasma , corresponding to the Res case in panels (a), (b), (c), and (d), and to the Res+TC case in panels (e), (f), (g) and (h) for times  = 25 s and  = 50 s.
The RMI generally occurs when the sign of the density and pressure gradients are opposite, i.e., ∇ • ∇ < 0, if a shock propagates from the light to the dense gas.Nevertheless, if a shock propagates from the heavy to the light gas, which is the scenario of our simulations, the RMI should occur when ∇ • ∇ > 0. However, in the Res case, we see that pressure and density gradients have the same signs just below the interface ( =∼2.1 Mm) at  = 25 s; therefore, the condition ∇ • ∇ > 0 it seems to be satisfied globally.In most of the points, we get ∇ • ∇ > 0. However, we do not see instabilities in all the maps for the Res+TC case; consequently, this analysis does not apply to this scenario.In panel (g), the vertical velocity gradient oscillates between positive and negative values below and over the interface up to about  = 2.5 Mm.However, at a higher distance, the gradient is positive, which indicates a predominantly upward flow, consistent with the behavior shown in the 2D map of the vertical velocity in panel (f) of Fig. 8. Finally, in panel (h), we note that plasma  is always smaller than 1, which means that magnetic pressure dominates over fluid pressure, and therefore, the plasma is flowing in this for the two times  = 25, 50 s.This behavior is consistent with any observed RMI since there is no transition between magnetically-dominated and fluid-dominated regimes.
According to the previous analysis, we can conclude that the internal structure of the loop region is related to an RMI in the Res case.In the Res+TC case, the thermal conductivity determines the plasma dynamics in the loop region, so the plasma moves through the loops and dissipates the heat, and therefore, any RMI develops.

Shock development
To investigate if there are the physical conditions so that a shock could develop in the post-flare loops structures of any of the two simulation cases, we estimate the sonic Mach number,   =   (  number, and the Alfvénic Mach number,   =   /  , where,   is the local vertical velocity of the plasma,   = √︁  /, is the local speed of sound, and   = |B|/ √︁ 4, is the local Alfvén speed.For example, in panel (a) of Fig. 10, we display the sonic Mach number and the Alfvénic Mach number for the Res case at  = 50 s.There, we see that the sonic Mach number is greater than one, i.e., the plasma is supersonic in the bottom regions of the loops, where the plasma propagates upwards.We also note in panel (b) that the Alfvénic Mach number is greater than one in regions where the plasma propagates upwards, which means that plasma flow is superAlfvénic, and as a consequence of the two later conditions, a shock develops in the Res case.On the other hand, in panels (c) and (d), we show the sonic and Alfvénic Mach numbers for the Res+TC case.In this scenario, it is discernible that plasma is neither supersonic nor superAlfvénic; therefore, the environment is unsuitable for the shock to develop.The results shown for the Res case in Fig. 10 are consistent with those obtained in Chen et al. (2015) in the scenario of particle acceleration by a solar flare termination shock.

Magnetic islands formation
In the Res+TC case, we can zoom in, for example, the   and observe the CS closely to identify multiple magnetic islands, as shown in Fig. 11.In panel (a), we display a zoom region that covers  ∈ [−2, 2] Mm and  ∈ [9, 13] Mm, while in panel (b), we show a zoom region that covers  ∈ [−1, 1] Mm and  ∈ [15, 20] Mm.For instance, in panel (a), we note the development of three magnetic islands along the CS, i.e., in regions where   is high.These structures are also observed along the CS at higher altitudes as schematized panels (b).In both panels, we can identify separatrices regions that indicate the presence of the magnetic islands.These regions represent the most recent reconnected field line that separates regions with magnetic fields of different topology; these magnetic fields may be anti-parallel, while the reconnected fields are found where the oppositely directed fields have been interconnected.Generally, the separatrices meet at the X-line, surrounded by a diffusion region, i.e., a high current density region, as shown in the small plots indicated by the blue arrows in panel (a).These magnetic islands can be formed due to a Tearing-instability in the non-linear regime.The tearing-instability is a linear stability of a CS of infinite (Furth et al. 1963), which is characterized by making a system in equilibrium unstable due to the tearing mode, leading to the formation of X-points and plasmoids during the reconnection process.This instability is part of a basic resistive instability that is linked with a long-wave tearing mode, corresponding to the breakup of the layer along current-flow lines (see, e.g., Furth et al. 1963).It can occur on timescales , such that   <  <   .The latter times are the time last the instability to traverse the CS at the diffusion speed   and the Alfvén speed   , respectively, and it has the effect of creating many small-scale magnetic islands in the CS (Shen et al. 2011).Besides, these magnetic islands are likely to develop in solar flares, (see, e.g., Kliem et al. 2000;Wang et al. 2021).
To complement Fig. 11 and go into more detail in the formation of the magnetic islands due to the Tearing-instability, in Fig. 12, we show a zoom region that covers  ∈ [−2, 2] Mm and  ∈ [9, 20] Mm of the magnetic field lines (a) and the plasma  (b) at time  = 51 s.For example, in panel (a), we note the formation of various small-scale (∼ 0.5 Mm of width) magnetic islands, which can be related to regions of high (∼ 10 4 ) plasma  values.We can also see that magnetic islands are souring by low (∼ 10 −1 ) plasma  values, which means that overall magnetic islands are localized in a mixed region.The latter means that magnetic islands are more likely to form in regions where gas pressure dominates over magnetic pressure, i.e., in regions of high plasma .This result makes sense since the magnetic islands form in X-points or O-points, where, by definition, the magnetic field is zero.Besides, it is interesting that magnetic islands are surrounded by regions where the magnetic pressure dominates over gas pressure ( ∼ 10 −1 ), which implies that the magnetic field confines the material that forms the islands.Specifically, in the Res+TC case, the thermal conduction helps distribute the plasma along the magnetic field lines, leading to the generation of a thinner CS compared to the Res case.
Moreover, it is possible to indicate the evolution of the CS's width, which could help identify the formation of multiple X-and O-points in the CS, (see, e.g., Shen et al. 2011).For instance, in Fig. 13 (neutral) point at  = 10.4Mm, which in turn triggers the magnetic island (see, e.g., Shen et al. 2011).In this scenario, plasma pressure dominates over magnetic pressures inside the magnetic islands, as shown in the plasma  map of Fig. 12.Since magnetic effects are small, a sudden drop in current density is observed in the last seconds of the simulation.The appearance of magnetic islands in the last seconds of the simulation is also consistent with  minimum values (∼ 0.1 Mm).Besides, to identify the formation of six O-points along the CS, in Fig. 13(b), we display the plasma velocity along the y-direction for two times  = 43 s and  = 51 s.According to our simulation results, these times cover a period before and after the formation of magnetic islands.Each arrow in Fig. 13(b) indicates the location of each magnetic island and its respective vertical velocity, which is normalized with the velocity scaled factor  0 = 10 8 cm s −1 .Before forming the magnetic islands, the plasma vertical velocity increases with time and does not show significant variations.For example, at  = 51 s, which is the time when we observe the formation of magnetic islands, and from 10 Mm onwards, the plasma vertical velocity fluctuates.These fluctuations indicate the presence of the magnetic islands, known as plasma blobs (see, e.g., Shen et al. 2011).In Fig. 13(b), we also specify the location and vertical velocities of the magnetic islands themselves.We see an island that moves downstream and is close to  = 10.4Mm with a vertical velocity of about   = −0.17.The other five magnetic islands move upstream at  = 11.5 Mm,  = 14.8 Mm,  = 16.3Mm,  = 17.8 Mm, and  = 19.6Mm.They correspond to vertical velocities of approximately   = 0.08,   = 0.79,   = 0.7,   = 0.71, and   = 0.76, respectively.The plasma vertical velocity around these islands is high, although for the magnetic island of approximately  =16.3 Mm, the vertical velocity of the plasma around and the island are similar.
For other magnetic islands, the vertical velocity of the surrounding plasma is considerably higher, which is consistent with the found in, for example, Shen et al. (2011).
Additionally, to indicate how the variations of the vertical velocity shown in Fig. 13(b) are related to the magnetic island formation, in Fig. 14, we show the temporal evolution of the − component of the current density   at the locations of the first and second magnetic island as schematized by the arrows of Fig. 13(b).Essentially, these plots show that the behavior of   in the two magnetic islands is similar, i.e., all reach their maximum value of around  ≈ 43 s and then suddenly decrease to the final of the simulation.The increase in   is due to the appearance of magnetic islands.In the bottom of Fig. 14, we show   (colored in blue) calculated at  =13.2 Mm, where it is a point without magnetic island.Here, we note that the behavior of   is different from the plots where magnetic islands form since, in the final times of the simulation,   grows instead of suddenly decreasing.

CONCLUDING REMARKS
In this paper, we study the formation and dynamics of solar flares, including the post-flare loops in a local region of the solar atmosphere, for two scenarios, according to the classical 2D model.In the Res case, localized resistivity produces more substructure around the post-flare loops, related to the plasma moving upwards and downwards.In addition, we identify the development of a complex structure, which could be associated with a magnetic island at the top of post-flare loops.On the other hand, in the Res+TC case, the whole solar flare structure looks smooth, including the post-flare loops and no apparent substructures are developed, which means that thermal conductivity redistributes the plasma homogeneously.However, we observe the formation of multiple small-scale magnetic islands from the top of the post-flare loops and along the CS in the vertical direction.Additionally, we include the case of an anomalous resistivity without thermal conduction.In that scenario, any post-flare loops and thin CS are formed.Instead, we see a hot plasma structure forming that becomes wider and moves upwards while magnetic field lines reconnect at the coronal level.This kind of plasma feature could be related to a plasmoid.Furthermore, in the mass density maps, we identify the development of a small jet (∼ 3 Mm of height) moving upwards from the transition region.This feature is interesting since it could be related to the magnetic reconnection process at microscopic scales, which is triggered by anomalous resistivity, given a relevant result despite the limitations of the MHD to capture the microscopic mechanism associated with the magnetic reconnection process.However, it is clear that in the case with anomalous resistivity, we do not observe the formation of a typical solar flare and associated post-flare loops as indicated by the classical 2D model introduced by Yokoyama & Shibata (2001).
The development of the multiple magnetic islands in the Res+TC case is related to the Tearing instability in the nonlinear regime.This instability is characterized by making a system in equilibrium unstable due to the tearing mode, leading to the formation of X-points and plasmoids during the reconnection process.This instability is part of a fundamental resistive instability linked with a long-wave tearing mode, corresponding to the breakup of the layer along current-flow lines.Besides, in the Res+TC case, the magnetic islands form along the CS in the vertical direction, consistent with the results obtained by Wang et al. (2021).Moreover, we observe that magnetic islands tend to form in regions with high plasma  (∼ 10 4 ) but surrounded by low plasma  regions (∼ 10 −2 ), which suggest that on those regions plasma pressure and magnetic pressure must contribute.That is, the magnetic pressure confines the plasma, and the plasma pressure makes the structure of the islands.In the magnetic islands, the current density tends to increase, and the CS's half-width becomes thinner when forming.Finally, when magnetic islands form, we identify the generation of plasma bobs in the vertical velocity, which appears due to the flow generated by the magnetic reconnection process.
In the case of the analysis for the instabilities performed to identify the development of RTI and RMI, we conclude that in the Res case, the internal spike-shaped structures are likely related to an RMI since the general physical conditions that state this instability are satisfied.In the Res+TC case, we show that any condition for RMI or RTI is satisfied since the post-flare loops look smooth.Although our simulation setup for the Res case is more straightforward than the one used in Shen et al. (2022), we identify desirable conditions for developing the RMI/RTI instabilities since we use only a localized resistivity profile.This result is crucial since it could indicate that a simple model with magnetic resistivity can solely produce such instabilities observed in underdense plasma downflows associated with magnetic reconnection in solar flares.
A difference of the present study regarding, for example, Yokoyama & Shibata (2001);Takasao et al. (2015); Ruan et al. (2020); Shen et al. (2022), is the inclusion of a highly anisotropic thermal conductivity that smoothly varies between the classical and saturated thermal conduction regimes, which according to our knowledge, it had not been used in numerical simulations of solar flares, so far.Additionally, we adopt a model for the solar atmosphere that includes the minimum temperature of the chromosphere and the sharp gradient of the transition region, which is more realistic than the one adopted by the abovementioned works since they typically consider a simple tangent function to describe a smooth variation of temperature from the photosphere to the solar corona.Furthermore, in comparison to Mondal et al. (2024), where the authors studied the formation and evolution of plasmoids in coronal CS, subjected to an external velocity perturbation under the condition of uniform resistivity, in our paper, we employed a localized resistivity that triggers the magnetic reconnection and leads the formation of the CS.Besides, Mondal et al. (2024) indicates that CS initially goes in thinning and Petschek-type magnetic reconnection followed by tearing instability and plasmoid formation.In contrast, one of the main objectives of our paper focuses on studying the associated small-scale structures of solar flares without entering into detail about the reconnection mechanism.However, both analyses could complement each other and engage in the context of the highly dynamic behavior observed along the CS in flare-like structures.
Finally, the results of this paper could help to continue persuading the dynamics of the formation of solar flares based on the standard model.In particular, it opens the gate to investigate if the anomalous resistivity itself could generate plasmoids and jets, which complete the picture of the magnetic reconnection process using the MHD approximations.Additionally, the current results could help better determine that solar flares can develop substructures even in simple scenarios where, for example, turbulence and particle acceleration are not considered.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.Equilibrium temperature in Kelvin (left) and mass density in gr cm −3 (right) as functions of height .

Figure 2 .
Figure 2. Magnetic field components   (left) and   (middle) in Gauss, and plasma  (right) at the initial time of the simulations.

Figure 3 .
Figure 3. Temperature maps in kelvin overlap with the magnetic field lines at  = 25 s (panels (a) and (b)) and  = 50 s (panels (c) and (d)).In panels (e), (f), we show the mass density in gr cm −3 for the time  = 25 s, while in panels (g) and (h), we display mass density at  = 50 s.

Figure 4 .
Figure 4. Temperature maps in kelvin and mass density in gr cm −3 overlay with the magnetic field lines at  = 25 s (panels (a) and (c)) and at  = 50 s (panels (b) and (d)) for the anomalous resistivity case.

Figure 5 .
Figure 5.The  component of current density in statA cm −2 overlap with the magnetic field lines for the time  = 25 s (panels (a) and (b)) and  = 50 s (panels (c) and (d)).In panels (e) (f), we show plasma pressure in dyn cm −2 for the time  = 25 s, while in panels (g) and (h), we display mass density in gr cm −3 for  = 50 s, for the Res case (right panels) and the Res+TC (left panels).

Figure 7 .
Figure 7. 1D cuts along the vertical direction  of the −component of the current density   in statA cm −2 (top) and temperature in Kelvin (bottom) at times  = 25 s (left panels) and  = 50 s (right panels) the Res case (blue curves) and the Res+TC case (orange curves).

Figure 8 .
Figure 8. Zoom of the temperature in kelvin (top panels), and vertical velocity   in cm s −1 (bottom panels) overlap with vector velocity field at  = 50 s, corresponding to the Res case (panels (a) and (c)), and the Res+TC case (panels (b) and (d)).

Figure 9 .Figure 10 .
Figure 9. 1D cuts of pressure gradient in dyn cm −3 , density gradient in gr cm −4 , vertical velocity (  ) gradient in s −1 and plasma  along the vertical distance  at  = 0 Mm, corresponding to the Res case in panels (a), (b), (c), and (d), and to the Res+TC case in panels (e), (f), (g) and (h) for times  = 25 s and  = 50 s.The black vertical dashed line represents the interface between the dense and tenuous plasma around the transition region at  ∼2.1 Mm.

Figure 11 .Figure 12 .
Figure 11.Zoom of the  component of the current density in statA cm −2 at  = 50 s that schematize the magnetic islands corresponding to the Res+TC case.

Figure 13 .
Figure 13.Variations of  (magenta curve), and the − component of the current density   (black curve) at  =10.4 Mm as functions of time (a) and the vertical velocity along the vertical distance  at two times  = 43 s and  = 51 s (b) for the Res+TC case.

Figure 14 .
Figure 14.Temporal evolution of the −component of the current density   estimated at the locations of the first and second magnetic islands of Fig. 14 (black curves) and at  =13.2 Mm (blue curve).