The Lorentz Force at Work: Multiphase Magnetohydrodynamics throughout a Flare Lifespan

The hour-long, gradual phase of solar flares is well observed across the electromagnetic spectrum, demonstrating many multiphase aspects, where cold condensations form within the heated post-flare system, but a complete 3D model is lacking. Using a state-of-the-art 3D magnetohydrodynamic simulation, we identify the key role played by the Lorentz force through the entire flare lifespan, and show that slow variations in the post-flare magnetic field achieve the bulk of the energy release. Synthetic images in multiple passbands closely match flare observations, and we quantify the role of conductive, radiative, and Lorentz force work contributions from flare onset to decay. This highlights how the non-force-free nature of the magnetic topology is crucial to trigger Rayleigh–Taylor dynamics, observed as waving coronal rays in extreme ultraviolet observations. Our C-class solar flare reproduces multiphase aspects such as post-flare coronal rain. In agreement with observations, we find strands of cooler plasma forming spontaneously by catastrophic cooling, leading to cool plasma draining down the post-flare loops. As there is force balance between magnetic pressure and tension and the plasma pressure in gradual-phase flare loops, this has potential for coronal seismology to decipher the magnetic field strength variation from observations.


INTRODUCTION
Flares occur frequently in solar and stellar atmospheres, and signal energetic events in association with compact objects and black holes of all sizes (e.g.Shibata & Magara 2011;Günther et al. 2020;Aliu et al. 2012).The sudden brightenings seen across X-ray (or even more energetic) wavebands result from rapid release of magnetic energy through magnetic reconnection, where detailed kinetic plasma processes boost photon emissions and generate near-light-speed particles (Aschwanden et al. 2017;Pontin & Priest 2022).The bulk of the plasma radiates in extreme ultraviolet (EUV) and soft X-ray (SXR) channels, and this thermal plasma emission can be enhanced by several orders of magnitude when a solar flare occurs.Solar flares dramatically impact the ionosphere of the Earth and can harm radio communication and navigation systems (Tsurutani et al. 2009).The study of solar flares led to the standard flare scenario based on lightcurve variations and detailed spatio-temporally resolved images.This standard flare model has been modeled frequently in 2D magnetohydrodynamic settings where a vertical current sheet is evolved to a flare loop configuration, and is recently revisited with a number of advanced 3D simulations (Shen et al. 2022;Ruan et al. 2023;Wang et al. 2023), as well as with 2.5D models where self-consistent two-way coupling between electron beams and the multidimensional MHD scenario is incorporated (Ruan et al. 2020;Druett et al. 2023b), or where post-flare coronal rain is obtained in multi-D settings (Ruan et al. 2021).It should be noted that these models deliberately adopt a simple initial magnetic topology, and as yet do not involve the actual full flux rope eruption process, but rather concentrate on the realistic thermodynamic and energetic evolution in and below the current sheet region.This should be contrasted with the simulations (in zero-beta settings, i.e. without thermodynamics) of actual eruptions that have provided insight into the evolving topological evolutions where an initial unstable flux rope is present (Aulanier et al. 2012).Such eruptive models, extended to full thermodynamics and concentrating on the current sheet fine-structure (Ye et al. 2023), are currently feasible thanks to adaptive mesh refinement strategies provided by our open-source code MPI-AMRVAC (Keppens et al. 2023).
A typical solar flare has two phases: the impulsive phase lasting few minutes and the gradual phase, lasting up to several hours (Fletcher et al. 2011;Hudson 2011;Benz 2017).Magnetic reconnection rapidly forms dense (electron number density N e ∼ 10 10 cm −3 ) and hot (temperature T ∼ 10 MK) coronal flare loop systems at the impulsive phase, during which a large amount of high energy electrons are generated and the electrons produce strong hard X-ray emission through the Bremsstrahlung mechanism.The density and temperature slowly return to pre-flare levels (N e ∼ 10 9 cm −3 , T ∼ 1 MK) in the long-duration gradual phase.Although most of the EUV and SXR radiation is released during this extended gradual phase, we still know very little about it, but remote observations across EUV and SXR wavebands present us with many dynamically evolving aspects, revealing clear multi-phase, multi-thermal plasma evolutions.A recent spectroscopic study of coronal rain formation in an X2 solar flare provides new clues on how cool condensations form in the hot flare system and then dynamically evolve (Brooks et al. 2024).Here, we aim to model this intruiging aspect of the long-term flare evolution in a self-consistent manner.Recent radiative MHD models from photosphere to corona, where interacting sunspot regions are forced to pass one-another, also demonstrate flares when the sunspots have close encounters, with indications of rain forming (Rempel et al. 2023), but our current standard flare setup allows for a clear cause-and-effect study that is fully reproducable and more easily analyzed.
During the gradual phase, flare loops are thought to gradually cool and become denser due to an interplay of thermal conduction and radiative losses, and eventually cold coronal rain forms in the cooling, dense postflare loops and this rain drains along the field to the solar surface (e.g.Schmieder et al. 1995;Mason & Kniezewski 2022).Recent fieldaligned flare simulations show that coronal rain cannot form in flare loop systems when the magnetic field configuration does not vary, which indicates that multi-dimensional effects are crucial for gradual-phase flare development (Reep et al. 2020(Reep et al. , 2022)).Since the duration of the gradual phase can be several hours (e.g.Reep & Knizhnik 2019), with a cooling rate lower than theoretically predicted (e.g.Bian et al. 2016b), additional heating mechanisms or factors inhibiting cooling may be present during the gradual phase.On EUV and SXR images during the impulsive and the gradual phase, a fan of bright spike-like rays is frequently seen above flare loops, and these rays demonstrate waving (e.g.McKenzie & Hudson 1999;Khan et al. 2007).The coronal rays are thought to be associated with the so-called supra-arcade downflows (SADs) phenomenon, which represents downward-moving dark blobs that are often seen between the bright rays and are interpreted as sub-Alfvénic downflows (e.g.McKenzie & Savage 2009;Savage & McKenzie 2011;Warren et al. 2011;Reeves et al. 2017;Samanta et al. 2021;Xie et al. 2022;Tan et al. 2023).Waving coronal rays suggest that gradual flare loop systems are quite dynamic, pointing towards three-dimensional magnetohydrodynamic processes at play.
We present a three-dimensional (3D) magnetohydrodynamic simulation of a solar (or stellar) flare event which covers both the impulsive and the gradual phases.Our simulation captures many phenomena seen during the gradual phase, including coronal rain and waving coronal rays.We show that both rain and waving rays form spontaneously under the action of the Lorentz force, in combination with radiative losses allowing thermal instability.The energy stored in the non-force-free magnetic field of the flare loop system is slowly released through work done by the Lorentz force in the gradual phase, a process that has the potential to extend the duration of the flare.Our paper is organised as follows.We discuss our model ingredients in Sect. 2. All our results and analysis are performed in Sect.3. We conclude in Sect. 4.
The MHD simulation is carried out using the open-source MPI-AMRVAC code (Keppens et al. 2023).The evolution of a solar flare is simulated in a square box with a domain of −50 Mm ≤ x ≤ 50 Mm, −50 Mm ≤ y ≤ 50 Mm and 0 ≤ z ≤ 100 Mm.Here z is the vertical direction and the x − y plane is parallel to the solar surface.An equivalent resolution of 512 × 512 × 512 is achieved with a block-adaptive mesh, where a 4-level adaptive refinement starts with a level-1 mesh of resolution 64 × 64 × 64.The cells in the flare loop systems and the reconnection current sheet have a size of 195 km × 195 km × 195 km.The contributions of thermal conduction, radiative loss, and gravity have been taken into account.The governing equations are given by where ρ, v, B, e, p tot , T , J are density, velocity vector, magnetic field vector, energy density (sum of kinetic energy, thermal energy and magnetic energy), total pressure (sum of thermal pressure and magnetic pressure), temperature, and current density vector, respectively.The gravity acceleration is given by g = −274R 2 s /(R s + z) 2 ẑ m s −2 , where R s = 696.1 Mm is the radius of the sun.A thermal conductivity tensor κ = κ ∥ T 5/2 bb is employed, where κ ∥ = 8 × 10 −7 erg cm −1 s −1 K −7/2 and b is a unit-tangent to the magnetic field vector.Saturation of thermal flux is employed (Xia et al. 2018).The optically thin radiative loss is given by where N e is electron number density, N i is ion number density and Λ(T ) is an optically thin cooling curve from Colgan et al. (2008).Where the height is greater than 6 Mm, an extension of the cooling curve at low temperatures is activated (see curve 'Colgan DM' in Hermans & Keppens (2021)), allowing the coronal rain plasma to cool to 10, 000 K or below.In this simulation, we assume that the number ratio of hydrogen to helium is 1 : 0.1, so we have N i = 1.1NH and N e = 1.2NH , where N H is the number density of hydrogen (protons).A background heating H b is adopted to compensate for the radiative loss outside the flaring region.The background heating is much weaker than the radiative loss, thermal conduction, or Lorentz work and therefore does not play an important role in the development of the flare loop systems, as will be clear from our detailed analysis.This background heating is necessary to sustain the corona against the radiative losses, in the region away from the flare site where we settle on a stationary stratified chromosphere to corona.Ohmic dissipation is included to reproduce the magnetic reconnection process, where η is the resistivity.In a finite volume approach, the conservative form of MHD equations (Eqn.( 1) -( 4)) is numerically solved using the classical 'HLL' Riemann solver (initials of authors Harten, Lax and van Leer in Harten et al. (1983)).Flux limiters are employed during the calculation of cell fluxes to avoid unphysical oscillations near high gradient regions.The robust second-order limiter 'vanleer' (from van Leer (1974)) is employed at the flaring regions and in the low atmosphere (the local mesh has a refinement level greater than three there).The high-precision third-order limiter 'cada3' (from first author of Čada & Torrilhon (2009)) is employed in the corona outside the flaring region (where the local mesh refinement level is less than or equal to three).The field-line-based transition region adaptive conduction (TRAC) method is employed for the low atmosphere (below height h ≤ 5 Mm) in our simulation to ensure that coronal density and temperature evolve correctly (Johnston et al. 2020;Zhou et al. 2021).

Initial conditions
Our model covers the chromosphere, the transition region, and the corona.The initial density and temperature profiles are obtained from a 2.5D atmosphere relaxation simulation with a uniform upward magnetic field, in which the C7 temperature profile from Avrett & Loeser (2008), a density profile calculated based on hydrostatic equilibrium, and an adaptive background heating are employed.More details on the relaxation are available in Ruan et al. (2021) and the results are demonstrated in Fig. 1.The background heating shown in Fig. 1 is our source term H b .The initial magnetic field is also quantified in Fig. 1.The initial magnetic field is a force-free field given by the following equations: B y = 0, (7) where B 0 = 30 G and λ = 10 Mm.The magnetic field reverses at the plane y = 0, which allows for magnetic reconnection.

Resistivity prescription
A time-varying resistivity is used in our simulation, which is given by where η 1 = 10 −2 , η 2 = 10 −3 , w ηy = 10 Mm, w ηz = 15 Mm, and h η = 30 Mm.A localized resistivity at the first stage t < 450 s dissipates the wide initial current sheet, resulting in the formation of closed arcades at the low atmosphere below h η and a thin current sheet with a large magnetic field gradient in higher atmospheric regions.Thereafter, sustained magnetic reconnection occurs at the thinned current sheet, where we can adopt a uniform, low value of resistivity η 2 and let the reconnection evolve by a combination of fully resolved and numerical contributions.
Anomalous resistivity, where a critical current threshold activates resistivity much stronger than that from Spitzer's theory, is often used in magnetic reconnection simulations to reproduce fast reconnection (e.g.Yokoyama & Shibata 2001;Ruan et al. 2020).In magnetic reconnection, the physics behind this anomalous resistivity is that microscopic instabilities can lead to this resistivity model (Yokoyama & Shibata 2001).It is still unclear whether this also works outside the reconnecting current sheet, such as in flare loops.In this work, we focus on flare loop systems in their gradual phase, rather than on reconnection in the current sheet.We therefore do not use anomalous resistivity after the initial phase of 450 seconds in this simulation, avoiding strong heating in the flare loop systems caused by such anomalous resistivity.Instead, we employed a relatively weak uniform resistivity in the gradual phase.In the simulation period t < 450 s, our adopted resistivity is only spatially located, and its magnitude η ∼ 10 −2 obviously dominates the reconnection process.Given our effective resolution, our effective numerical resistivity can be estimated to be in the order of magnitude of 10 −3 or lower.We employ an explicit uniform resistivity η = 10 −3 in the period t ≥ 450 s to control the overall numerical resistivity at play (combined explicit as source terms and numerical).

Boundary conditions
Periodic conditions are employed at the x-boundaries.Symmetric conditions are used for density, v x , v z , pressure, B x and B z at the sideways y-boundaries, and anti-symmetric conditions are employed for v y and B y at y-boundaries.
Here, a symmetric condition indicates Γ(s) = Γ(−s), while anti-symmetric implies Γ(s) = −Γ(−s), for scalar Γ and where s = 0 is at the boundary.Fixed conditions are employed at the bottom boundary (at z = 0): the initial temperature at the lower boundary ghost cells is extrapolated linearly from the C7 temperature profile; the density is obtained based on hydrostatic equilibrium; the velocity is set to zero; and the magnetic field is given by Eqn. ( 6) to (8).Zero gradient extrapolation is adopted for density at the upper boundary.Equal gradient extrapolation is employed for thermal pressure in the upper boundary ghost cells, where p/(dp/dh) = 50 Mm as function of height h = z.Zero gradient extrapolation is employed for velocity and magnetic field components at the period t < 540 s, which leads to an open boundary.The velocity and magnetic field in the upper boundary ghost cells are treated differently from t = 540 s onwards, when we switch to an effectively closed top boundary prescription: velocity is set to zero, a symmetric condition is employed for B x and B z , and the anti-symmetric condition is employed for B y .This time-specific modification of the upper boundary is employed to limit the otherwise unrealistic infinite vertical expansion of the flare loop system.If the boundary conditions are not changed, the reconnection magnetic energy release rate will gradually decrease as the magnetic field strength and topology of the current sheet change.However, the expansion of the flare loops would not stop when they reach our upper boundary.The closed upper boundary forces the current sheet to represent a structure that is narrow in the middle and wider at both ends.In this case, due to the bending of the magnetic field, there will be an outward Lorentz force at the reconnection inflow region, thus suppressing the reconnection process and flare loop expansion.In real flares, the existence and associated expulsion of a flux rope will naturally introduce such effects.

Observational reference data and forward modeling
We compare our simulated flare evolution in various passbands, and use selected reference data from (different) flare events to compare with: to that effect we also produce synthetic views on our numerical flare evolution.We will do so for EUV and Hα images.
Observations of solar flare loop systems in EUV passbands (Fig. 2A-C) are obtained by SDO/AIA (Lemen et al. 2012).The AIA/EUV images have a cadence of 12 s and a spatial sampling of 0.6 ′′ /pixel.The aia prep routine in the Solar Software (Freeland & Handy 1998) is applied to align AIA images.The AIA 304 Å image in Fig. 2A was observed at 18:58:08 UT on June 22, 2015.The AIA 211 Å image in Fig. 2B was observed at 16:16:37 UT on September 10, 2017.The AIA 131 Å image in Fig. 2C was observed at 14:30:34 UT on May 22, 2013.We note that these different events also sample different flare classes, but their overall morphological evolution can be considered similar, as also evidenced through standard flare model efforts presented in Druett et al. (2023a).
The possibility to generate EUV emission synthesis has been included in the open-source version of the MPI-AMRVAC 3.0 code (Keppens et al. 2023).Corresponding functions from the CHIANTI atomic database (version 8) are employed in the synthesis (Del Zanna et al. 2015).Optically thin assumptions are employed.In the high-density lower atmosphere (z < 5 Mm) that deviates from the optically thin condition, the EUV emissivity is set to zero.A point spread function has been used to include scattering effects by the instrument, and the synthesized images have the same resolution as the observations.The method given in Pinto et al. (2015) is employed in the calculation of GOES SXR flux.
When synthesizing an Hα image, we assume that at a height of 5 Mm there is Hα radiation of unit intensity propagating along a supposed line of sight (LOS) direction.The Hα light is assumed to be partially absorbed when it meets the cold rain materials, and the absorption rate is calculated with the method in Heinzel et al. (2015).The impact of plasma motion is neglected in this calculation.An observation of flare-driven coronal rain at Hα wavelengths is found in Jing et al. (2016).

Overall evolution from flare onset to decay
We adopt a standard 3D flare model setup that evolves a current-sheet to its impulsive phase, including its realistic thermodynamic stratification in solar atmospheric conditions.After about seven minutes, an induced fast reconnection sets in, and the flare goes into its impulsive phase.By the 12th minute, the thermal SXR emission reaches its maximum value (Fig. 2G).After that, the SXR flux gently drops and the flare moves into its gradual phase.The simulated flare is a C-class flare, according to the classification of solar flares based on the peak GOES SXR flux. Figure 2 and the associated movies demonstrate that we recover multiple typical observational behaviour across EUV wavebands throughout the entire gradual phase.Note how the AIA 304 and 211 views are representative for the late gradual phase, where we have obvious condensations in direct agreement with the selected flare observations.In the AIA 131 view, taken in the early decay phase of the flare, our side-on view clearly shows coronal rays above the arcade system.Our animated views (Movie Fig2) show the temporal evolution of the flare loop system in 304 (from 10 to 29 s of the movie) and in 131 and 193 (0 to 10 s of the movie) wavebands, while also addressing the role played by the LOS, by rotating the scene.
Typical temperature, density, and magnetic field configurations in the gradual phase are displayed in Fig. 3. Magnetic reconnection occurs above the flare loop systems.While the loops at the outer layers are newly formed, those at the inner layers were generated before or during the early impulsive phase.Newly formed loops get filled with hightemperature and high-density evaporative flows, and their temperature gradually decreases due to heat conduction and radiative cooling.The loops become shorter during this cooling process (Movie Fig3).During the gradual phase, the reconnection rate is slow, and the rate of energy injection is lower than the rate of energy loss caused by heat conduction and radiative cooling (this will be shown in Section 3.3), so the temperature of the coronal flare region generally decreases slowly.Simultaneously, the density of this region drops over time, due to the continuous loss of material from the corona to the chromosphere (Bradshaw & Cargill 2005).

Non-force-freeness of the evolving arcades
The newly formed magnetic arcades at the outer layer of the flare region are far from a force-free stage, and the Lorentz force (J × B) is pointing inward.Under the influence of this downward Lorentz force, the arcades are forced to shrink, meanwhile compressing the plasma below, until the Lorentz force and the pressure gradient (−∇p) reach a force balance.This can be quantified by introducing suitably averaged force measures, where we acknowledge that our initial setup is invariant along the x-direction (i.e.identical in all y − z planes): the magnetic field has initially mostly z components, and only within |y| < λ is the field showing out-of-plane x-components B x (y, t = 0).All variations in the x direction are due to spontaneously ensuing turbulent fluctuations, as also our initial anomalous resistivity is y − z-dependent only.After the initial 450 seconds, low-lying loops formed with again essentially purely y − z variations.Therefore, we can meaningfully average out the x-variation at all times in our study, to get the essential force-balance achieved across the arcade.This is true throughout the 3D flare evolution, as our magnetic topology closely follows the cartoonized 2.5D standard flare model.Indeed, cross-sectional views at varying x values at a fixed time look overall similar, while differing in details.Fig. 4 shows the distributions of the mean Lorentz force, gravity, and pressure gradient force on the y − z plane.Since an average in the x direction is involved, we refer to these forces as mean forces.Averaging not only better shows the force distribution, but also smooths out the influence of the occuring turbulence.Here we give the details of the calculations for these mean forces.
The calculation of the mean Lorentz force is based on the mean magnetic field configuration.The mean magnetic field is given by for field components along s = x, y, z, where (x min , x max ) is the domain in the x-direction.The vertical component of the mean Lorentz force shown in Fig. 4A is then given by FB,z = Jx By − Jy Bx = ( where Jx and Jy are current density components calculated from the mean magnetic field, and the x-variation of the mean magnetic field is zero (in accord with our statements above on the similarity across x).The mean gravity shown in Fig. 4B is given by Fg,z (y, z) = −g(z) where ρ is plasma density, g(z) = 274 (z/R sun ) 2 m s −2 , and R sun = 696.1 Mm is the solar radius.The mean pressure gradient force is calculated from the average thermal pressure, hence it is given by p(y, z) The vertical component of the mean pressure gradient force shown in Fig. 4C is then simply Moreover, a line-based analysis is also shown in Fig. 4, where the black line in Fig. 4D is a 'magnetic field line' derived from the average magnetic field ( By , Bz ), and the blue line is a vertical line from (y = 0, z = 3 Mm) to (y = 0, z = 40 Mm).Fig. 4E shows the density distribution (black crosses) along this magnetic 'field line', derived from the temperature curve and gravity stratification (ρg − ∇p = 0).This computed (black crosses) density profile is given by where the symbol × indicates our derived variable, l and l is the integral pathlength (tangent of the 'field line'), s denotes location along the path, where s 0 is the integral starting point located at z = 3 Mm on the 'field line'.Furthermore, g = −ge z , R is the universal gas constant and the average temperature is computed from Eq. ( 15) uses the following derivation of the pressure profile: The perpendicular 'field line' density distribution (blue triangles) derived from the temperature curve and pressure gradient force-Lorentz force balance (J × B − ∇p = 0) is also given in Fig. 4E.This predicted vertical (blue triangle) density profile is given by where the integral path is the blue vertical line in Fig. 4D, and symbol △ indicates our derived variable, where s 0 is the endpoint of the line located at z = 3 Mm.Our movie accompanying Fig. 4 quantifies these mean forces for 7 minutes starting after t = 14 minutes.At all times, vertical Lorentz forces balance pressure gradients, so a non-force-free state is evidently realized throughout.Indeed, our model shows that the pressure gradient and the Lorentz force in gradual-phase flare loops are closely balanced (J × B − ∇p ≈ 0), while gravity is much weaker than these two forces (Fig. 4A-D and its Movie).Such a force balance controls the cross-magnetic-field density distribution (Fig. 4E), with typically higher density in the inner layers (also seen in Fig. 3B).The field-aligned density distribution is governed by a hydrostatic balance instead, namely that of gravity and pressure gradient forces (Fig. 4E).As a result, density along a loop tends to decrease with height.
The cool and dense loops in the lower layers represent the ideal environment to form coronal rain.Radiative loss eventually dominates the cooling in these loops (see section 3.3).Radiative-loss-driven catastrophic cooling initiates in the dense loops after they cooled to about 1 MK (Field 1965;Cargill & Bradshaw 2013).As a result, sudden condensations occur, creating colder (∼ 0.01 MK), denser rain blobs (Fig. 3 and its movie).These rain blobs appear as bright spots in EUV 304 Å images (Fig. 2A,D and its Movie Fig2A) and dark in Hα images.In Fig. 5 and its movie, we show how the development of the rain is clearly seen in isotemperature animations (showing regions with temperature lower than 25, 000 K and proton number density higher than 10 9 cm −3 ).Fig. 6 and its movie translates these simulations to an Hα proxy, as explained in Section 2.5.These views qualitatively agree with the stranded nature of the field-aligned post-flare rain condensation observed in many flare events, e.g. in the work by Jing et al. (2016).Through the formation of coronal rain, the inner loops gradually become empty (see also Movie Fig3), as suggested by current flare models (Priest & Forbes 2002).
We note that there are other ways to calculate the mean Lorentz force, e.g. one could calculate the Lorentz force from the original current and magnetic field and then average the Lorentz force (Appendix A).We verified that for our setting (typified by an invariant direction initially), these two methods qualitatively agree, although they obviously differ in details (especially where turbulence is active).The main force balance conclusions we noted are completely identical between both methods.The advantage of the method used in our article is that it can directly show the relationship between the average magnetic field and the average pressure.In flare observations, we may thus be able to infer the gradual-phase magnetic field based on this relationship.

Energetic analysis
The delicate balance between Lorentz force and pressure gradient evolves dynamically during the gradual phase.Influenced by downward enthalpy fluxes, thermal conduction, and radiative losses, the pressure in the flare region is continuously decreasing (Bradshaw & Cargill 2010).The Lorentz force will thus temporarily be stronger than the pressure gradient, and the magnetic arcades move down, compressing the plasma and reaching a new equilibrium.The Lorentz force does work during this procedure, releasing the energy in the magnetic field and making up for the internal energy lost.Instead of being converted directly into internal energy (as is the case for Ohmic heating in reconnection regions), the postflare magnetic energy is first transformed into kinetic energy and then, by compression, partially into internal energy.The Lorentz work, which is of the same magnitude as the radiative losses and thermal conduction (Fig. 7A,C), thereby helps to slow down the cooling of the flare region.Our C-class flare adopted a coronal magnetic field strength of B 0 = 30 G, but in X-class flares (e.g.Fleishman et al. 2020), magnetic field strengths can be an order of magnitude higher.The energy provided by the work done by the Lorentz force (∼ B 2 0 ) will then be two orders of magnitude higher and may effectively prolong the duration of the gradual phase.This is consistent with the finding that in X-class flares, the magnetic flux from the footpoints of flare loop systems has a positive correlation with flare durations (Reep & Knizhnik 2019).On the other hand, prolonging the duration of lower-class flares might be dominated by other reasons, such as on-going reconnection and long-lasting cooling (e.g.Cargill & Priest 1983;Liu et al. 2013;Li et al. 2014;Toriumi et al. 2017;Reep & Toriumi 2017), since the duration of small flares appears to be independent of the magnetic flux (also see Reep & Knizhnik (2019)).The way we quantify the detailed energy evolution is given below.
The distributions of Lorentz force work, radiative loss, and thermal conductive heating/cooling on the y − z plane are shown in Fig. 7A,B.The method to calculate the dimensionally-reduced distribution of the energy conversion rates is as follows: FRadiative (y, z) FConductive (y, z) Note that the conductive heating or cooling rates given in Fig. 7B are here estimated in post-processing, without accounting for specific numerical effects, like the heat flux limiter or cooling flux adjustments (TRAC aspects), as well as flux fixes active at the borders of AMR grids.Fig. 7C shows the instantaneous 3D-space-integrated energy release rates, and their calculation happens as follows:  Birn et al. (2009) emphasizes the role of the Lorentz force work in converting reconnection energy.We also here briefly analyse the energy conversion in the combined system during the early impulsive phase, where it is known that an MHD approximation does not capture essential non-thermal aspects of the flare evolution (Ruan et al. 2020).In our 3D MHD simulation, the Lorentz force work at the sideways bounding slow shocks in region S2 is primarily responsible for the energy release in magnetic reconnection, consistent with the original Petschek steady, fast reconnection model.In panels (B) and (C) of Fig. 8, respectively, the spatial distributions of the Lorentz force work and Ohmic heating at the impulsive phase are displayed.In panel (D) of Fig. 8, integrals of Lorentz force work and Ohmic heating in region S2 are compared.The release of magnetic energy is much more effectively accomplished by the work done by the Lorentz force than by the Ohmic heating, as seen in Fig. 8, panels (B)-(D).Figure 8 shows that the effect of ohmic heating caused by the uniform resistivity is much weaker than that of Lorentz force work.The effect of this ohmic heating is negligible for the evolution of flare loop systems.
For completeness, the energy evolution of the flare loop system in region S1 is analysed simultaneously.The evolution of the thermal energy in flare loop systems tracks the evolution of the energy released through reconnection very well during our impulsive phase, and the thermal energy is close to half of the released magnetic energy (Fig. 8E), as the other half would go into the (eventually ejected) fluxrope system above.The thermal energy shown is given by and the released magnetic energy is quantified by This ratio is reasonable when considering that roughly half of the released energy is transported upward and leaves the simulation domain.To further understand the intricate relationship between the energy release in the reconnection region and the energy evolution of the flare loop systems, a detailed analysis of energy conversion and transport can be performed (e.g., Birn et al. 2009;Reeves et al. 2010).However, our main effort here concentrates on the gradual phase, and especillay the multi-thermal aspects at play that drive condensation formations, which is why we emphasize the role played by the non-adiabatic effects included in our study: radiation and conduction.
We conclude that magnetic energy can be slowly released in the flare loops through Lorentz work during the gradual phase.However, the role that this released energy plays in extending the duration of the flare remains to be evaluated.More detailed analysis, where we also vary this energy budget across different classes of flares is left for future work.

Instability analysis
Although the Lorentz force and pressure gradient establish a dynamically evolving force balance inside the gradual phase flare loops, this mechanical equilibrium is actually (linearly) unstable in the flare region where it is fully dominated by these two forces.We find that Rayleigh-Taylor instability inevitably arises (Fig. 9E), in an environment where the cross-field density distribution deviates from purely gravitational stratification (Fig. 4D and 4E) and the magnetic field is not potential.Previous studies invoked the high-speed reconnection outflow for inducing Rayleigh-Taylor instability in the impulsive phase and explaining the formation of SADs (Guo et al. 2014;Shen et al. 2022).Here we show instead that the occurrence of Rayleigh-Taylor instability in flare regions is unavoidable throughout the gradual phase, and does not require the participation of high speed reconnection outflows: the detailed 3D magnetic, density and pressure variations induced by the gradually evolving postflare field configuration inherently bring the entire region to Rayleigh-Taylor induced mixing.An uneven distribution of EUV emissivity reflects the instability effects on temperature and density distributions (Fig. 9 and Movie Fig9).Waving spike-like coronal rays appear in EUV images as a result (Fig. 2F and Movie Fig2).In our simulation covering the entire postflare gradual phase, internal energy feeds the Rayleigh-Taylor instability, and a significant portion of the internal energy comes from the chromosphere in the impulsive phase through evaporations (Fig. 7C).Therefore, an energy transfer pathway of "reconnection -fast electron deposition/thermal conduction -evaporations -coronal rays" is at play.In contrast, the main energy transfer path in Guo et al. (2014) and Shen et al. (2022) is stated as "reconnection -fast outflows -SADs".
We make a brief analysis of the Rayleigh-Taylor process in what follows.Flare loop systems generally have high temperatures and thus have a large scale height.Assuming a loop temperature of T = 5 MK, then the corresponding scale height is where k B = 1.38×10 −16 cm 2 g s −2 K −1 is Boltzmann constant, the mean mass of particle m = 8×10 −25 g is supposed to be half that of a proton, and gravity is supposed to be a constant g = 2.74 × 10 4 cm s −2 .As a typical peak proton number density inside flare loops is 10 10 cm −3 and that above flare loops is 10 9 cm −3 , the distance one estimates from a balance between gravity and pressure gradients in an isothermal setting for these density contrasts leads to a length estimate of Such a length is much larger than the size of flaring regions.Hence, the observed density contrasts on the smaller scale of the loop system must be maintained by means of the Lorentz force, where it balances a non-negligible pressure gradient.Although the force balance can be achieved by the Lorentz force, this balance is unstable and can easily lead to Rayleigh-Taylor instability.
Here we take Fig. 9 as an example to briefly analyze the Rayleigh-Taylor instability in the gradual-phase flare loops.Fig. 9A shows the density distribution at slice y = 0, and this distribution is shown more clearly in Fig. 10A.Finger-like structures are continually generated at a height of 20 Mm due to Rayleigh-Taylor instability (also see Movie Fig9).From the velocity distribution in Fig. 10B, it can be seen that there is no high-speed flow of hundreds of km s −1 near the roots of the structures, indicating that the formation of the structures is independent of SADs.The pressure gradient is large at this height (Fig. 10C,D), and the pressure gradient force promotes the upward growth of the structure.When there is a huge density difference on both sides of a interface, the growth rate of the Rayleigh-Taylor instability can be written simply as where k is wavenumber, a is acceleration and a = −|g| in gravity dominated Rayleigh-Taylor instability.Here, the pressure gradient force takes the role of gravity and a = ∇p • e z /ρ, where ρ is density.With a wave vector k perpendicular to the magnetic field, the magnetic field has little inhibitory effect on this instability.Suppose k = 1/(10 Mm), dp/dz = −10 −8 Ba cm −1 (refers to Fig. 10D) and ρ = 10 −14 g cm −3 (N H ≈ 6 × 10 9 cm −3 ), then an acceleration a = −10 6 cm s −2 and a growth rate of |ω| ≈ 1/30 s −1 are obtained.The instability hence grows on a sub-minute timescale.

SUMMARY
In summary, our magnetohydrodynamic model captures many of the mysterious multi-phase aspects witnessed in multi-frequency views of solar flare events, especially those occurring in the extended, hour-long gradual phase.We showed that models which do not account for the work delivered by Lorentz forces (such as those assuming pure fieldaligned processes) entirely miss the most important ingredient to explain the flare sustained and gradual energy release across EUV and SXR channels.We reproduce three-dimensional multi-thermal processes in unprecedented detail.We showed there is sustained force balance between the magnetic pressure and tension and the plasma pressure inside the slow-evolving gradual-phase flare loops.This could be used to figure out the magnetic field strength variation from EUV and SXR observations, if both plasma density and temperature can be deduced.Since the Lorentz force scales quadratically with field strength, stellar flaring events should also demonstrate a clear correlation between their gradual phase duration and their magnetic field strength.
We argue theoretically and demonstrate with our simulation that flare loop systems are in turbulent states due to Rayleigh-Taylor instability, when the loop density remains high.The coronal rays observed in the gradual phase may be the external manifestation of the Rayleigh-Taylor instability.We suggest that sub-Alfvénic downflows (SADs) can lead to the generation of coronal rays, but not all coronal rays are related to such flows.It is foreseeable that transverse waves will be continuously excited in the turbulent loop systems and transport a large amount of energy from the corona along the loops to the underlying atmosphere, accelerating the cooling of the flare loops.On the other hand, the turbulence has the potential to reduce downward thermal conductive heat flux, which in turn slows down the loop cooling (Bian et al. 2016a,b;Emslie & Bian 2018;Xie & Reeves 2023).The impact of the turbulence on the flare loop systems development is still to be studied, as turbulence inherently requires even more extreme resolution simulations.
In the observation reported in Jing et al. (2016), flare-driven coronal rain initially appears in the inner loops, and then the formation locations expands outward over time.Such a process is well captured in our simulation (see the Movie of Fig. 5).Our simulation demonstrates that the high density flare loop systems are gradually emptied from the inside out, due to the formation of coronal rain (see the Movie of Fig. 3).Possibly due to our limited (order 200 km) simulation resolution, the coronal rain in our simulation does not exhibit many fiber-like fine structures as the observation.For the next step, higher resolution will be employed in the numerical study of the gradual phase to capture more details, and a more realistic magnetic field topology will be employed, such as a magnetic field topology inspired from Aulanier et al. (2012) or a topology extrapolated from a magnetogram.Future work could improve on the neglect of non-thermal plasma aspects that are vital during the impulsive (and still present in the gradual) flaring phase, along the lines of self-consistent models unifying fast electron beam physics with multi-dimensional MHD modeling (Ruan et al. 2020;Druett et al. 2023b,a).
In this article we show the importance of the Lorentz force during the flare evolution, especially in the gradual phase and in the lower-lying postflare loops where ultimately postflare rain developed.For understanding how complex magnetic topologies may be liable to flare onsets, non-linear forcefree field extrapolations remain valuable, and nonlinear forcefree field approaches are very useful in modelling coronal magnetic fields (Wiegelmann & Sakurai 2012).However, the Lorentz force and how it acts to redistribute energy should be important in regions where pressure changes by orders of magnitude on a small length scale, like flare loops.Such kind of regions may only occupy a small part of the corona.13).Panel (D) gives profile of the mean pressure at the vertical line in (C), and the corresponding region in (A) and (B) is the region between the horizontal dotted lines.A fitting of the pressure profile is also given in (D), as well as the pressure gradient from the fitting.The corresponding time is the same as that in Fig. 9, that is, t = 21.5 min.
) where S1 indicates the 2D area where high-density flare loop systems are located in the y − z plane (referring to the region inside the solid curves in panels (A)-(C) of Fig. 8), S2 indicates the area around the reconnection current sheet in the y − z plane (referring to the region inside the dashed curves in panels (A)-(C) of Fig. 8), and S1b (a 2D horizontal surface) indicates the lower border of the region S1, Rec stands for magnetic reconnection, Lor stands for Lorentz force work, Rad stands for radiative loss, Con stands for thermal conduction and Eva stands for upward chromospheric evaporations.These time-varying integral regions S1 and S2 are shown in the Movie of Fig. 7.

Figure 1 .
Figure 1.Initial conditions of the simulation.The initial magnetic field lines are displayed in panel (B), where the background is density.Panel (A) shows how the magnetic field components change in the y-direction.The variations of proton number density (black solid line) and plasma temperature (red dashed line) in the vertical z-direction are given by panel (C).The time-invariant background heating is demonstrated in panel (D), where the initial profile of radiative loss is also given.

Figure 2 .
Figure 2. Solar flare loop systems at multiple passbands and multiple viewing angles.(A) An M6.5 flare observed at EUV 304 Å wavelength by the Atmospheric Imaging Assembly (AIA) onboard the Solar Dynamics Observatory (SDO)(Lemen et al. 2012).This wavelength is sensitive to plasma temperature of 0.05 MK.This case has been reported inMason & Kniezewski (2022).(B) An X8.2 flare observed at EUV 211 Å wavelength by AIA.The EUV 211 Å wavelength is sensitive to plasma temperature of 2 MK.This case has been reported in Martínez Oliveros et al. (2022).(C) An M5.0 flare observed at EUV 131 Å wavelength by AIA.This wavelength is sensitive to plasma temperatures of 0.4 and 10 MK.This case has been reported in Tan et al. (2023).(D) Synthetic EUV 304 Å image from our simulation.(E) Synthetic EUV 211 Å image from our simulation.(F) Synthetic EUV 131 Å image from our simulation.(G) Synthetic GOES SXR 1-8 Å flux from our simulation.In the AIA images, the distance between the axis tickmarks is 20 arcsec (14.5 Mm).Flare-driven coronal rains and fan of coronal rays are phenomena frequently observed in C, M, and X-class flares (e.g., Khan et al. 2007; Mason & Kniezewski 2022).An animation of this figure (Movie Fig. 2) is provided.The movie has a real-time duration of 29 s.The period 0 to 10 s of the movie shows the synthetic 131 Å view as well as 193 Å view (not shown in the figure) from simulation time 14 min to 21 min.The direction of the LOS is rotated after the time reaches 21 min.The period 10 to 29 s of the movie shows the synthetic 304 Å view from simulation time 35 min to 57 min.The direction of the LOS is rotated after the time reaches 57 min.

Figure 3 .
Figure3.Multi-phase aspects in the gradual-phase flare loop systems.We show temperature distribution (T in panel (A)), proton number density (NH in panel (B)), and magnetic field configuration (solid lines in both panels).The coloring of magnetic field lines is according to the local temperature (panel (A)) or density (panel (B)).Time is shown at the top of panel (A).Coronal rain blobs, which have proton number densities greater than 10 10 cm −3 and temperatures less than 15, 000 K, are visualized as isosurfaces seen in both panels.The surface coloring of these rain blobs gives the local plasma beta (β), which quantifies the ratio of thermal to magnetic pressure.An animation of this figure (Movie Fig3) is provided.The animation covers ∼19 minutes of physical time starting at t = 28 min (real-time duration 9 s).

Figure 4 .
Figure 4. Force analysis in the flare loop system.(A to C) The vertical (z) components of the (x-averaged) Lorentz force, gravity, and pressure gradient force, respectively.The unit of the forces is 10 −9 dyn cm −3 .(D) The summed vertical components of pressure gradient and Lorentz force.The innermost loops are stressed differently than the outer loops, as the inner loops are low-density loops that formed before the impulsive phase.The red dashed and the black solid line in panel (E) represent the average temperature and density profiles along the black solid magnetic field line in panel (D).Black crosses indicate the expected density distribution derived from the temperature and gravity stratification along that field line.The x-averaged temperature and density profiles in the direction perpendicular to the magnetic field lines (i.e.along the vertical blue solid line in panel (D)) are shown in panel (E) by the purple dotted and the blue dashed-dotted line, respectively.The blue triangles represent the expected density distributions derived from the temperature variation and the pressure gradient-Lorentz force balance.The forces along the vertical line in panel (D) are compared in panel (F).An animation of this figure (Movie Fig4) is provided.The animation covers ∼7 minutes of physical time starting at t = 14 min (real-time duration 5 s).

Figure 5 .
Figure5.Post-flare coronal rain.We show regions in darkgrey isosurfaces with temperatures below 25,000 K and of higher than 10 9 cm −3 proton number density, as well as the lower temperature variation across the transition region and selected post-flare arcade field lines.An animation of this figure (Movie Fig5) is provided, starting at 28.6 minutes, and lasting till 57 minutes.

Figure 6 .
Figure 6.Post-flare coronal rain in an Hα proxy.We show an angled view on the flare loop system.An animation of this figure (Movie Fig6) is provided, starting at 35.8 minutes, and lasting till 57 minutes.

Figure 7 .
Figure 7. Energetic throughout the postflare system.Panel (A) compares spatial distributions of radiative loss and Lorentz force work, while panel (B) compares radiative loss with conductive cooling or heating.For coordinate axis information, refer to Fig. 4. The dimensional heating or cooling rate is given in erg cm −3 s −1 , and uses a multi-color representation that shows the dominant component visually.The time evolution of radiative cooling, conductive cooling, Lorentz force work, energy release rate by reconnection, and the energy flux brought into the flare loop systems by the chromospheric evaporation are shown in panel (C).The cyan vertical dotted line in panel (C) is at the time corresponding to panels (A) and (B).The orange solid curve in panel (C) gives the GOES SXR flux.For detailed computations, we refer to Sect.3.3, while an animation (Movie Fig7) is provided.The animation covers ∼64 minutes of physical time starting at t = 0 (real-time duration 11 s).

Figure 8 .
Figure 8. Energetic analysis of our impulsive flare phase.Proton number density (A), Lorentz force work (B), and Ohmic heating (C) distributions in the impulsive phase.The area is as shown in Fig. 1B.The time shown in panels (A)-(B)-(C) refers to the vertical dotted line in panels (D) and (E).The time-evolving spatially-integrated (over area 2 in (A)) energy release rates due to Ohmic heating (I Ohm ) and Lorentz force work (ILor) are depicted in panel (D), where the unit is erg cm −3 s −1 .Comparison of the thermal energy (E th ) in flare loop systems (over area 1 in (A)) with half of the energy released by magnetic reconnection (half Erc).The GOES SXR flux is given in both panels (D) and (E) by the orange curve.

Figure 9 .
Figure 9. Rayleigh-Taylor instability inside flare loop systems.The cartoon in panel (E) illustrates the mechanism.The proton number density, temperature, thermal pressure, and vertical component of the Lorentz force are shown in panels (A) to (D).Magnetic field lines are represented by the curves in panels (A) to (D), and the surface color of the field lines in panels (A) to (C) shows the local vertical component of the Lorentz force.An animation (Movie Fig9) is provided.The animation covers ∼14 minutes of physical time starting at t = 11 min (real-time duration 5 s).The proton number density, temperature, thermal pressure, and vertical component of the Lorentz force, EUV 131 Å emissivity, and EUV 193 Å emissivity are shown in panels (A) to (F) of the movie.

Figure 10 .
Figure 10.Instability analysis of gradual phase flare loops.Proton number density (A) and vertical velocity (B) at slice y = 0. Panel (C) shows mean pressure distribution on y-z plane given by Eqn.(13).Panel (D) gives profile of the mean pressure at the vertical line in (C), and the corresponding region in (A) and (B) is the region between the horizontal dotted lines.A fitting of the pressure profile is also given in (D), as well as the pressure gradient from the fitting.The corresponding time is the same as that in Fig.9, that is, t = 21.5 min.