Numerical schemes and dissipative processes in three-dimensional modeling of post-flare loops and their EUV signatures

1 Introduction

Solar flares are explosive events that convert magnetic energy into internal energy and kinetic energy, thereby modifying the plasma and magnetic field conditions in the solar corona. Observational evidence indicates that magnetic reconnection is a key mechanism driving this energy conversion and the subsequent evolution of solar flares. Consistent with this interpretation, the standard reconnection-driven flare model has been developed and subsequently applied in numerous theoretical and numerical studies of solar flare dynamics (e.g., Carmichael, 1964; Sturrock, 1968; Kopp and Pneuman, 1976; Biskamp, 1986; Tajima and Shibata, 2002; Reeves et al., 2010; Priest, 2014; Takasao et al., 2015).

Numerical studies in two-dimensional (2D) and two-and-a-half-dimensional (2.5D) geometries have examined flare reconnection within the framework of non-ideal MHD, including models with anomalous resistivity (Wu et al., 2010; Mei et al., 2012). Xie et al. (2021) performed 2D MHD simulations to investigate oscillations and plasmoid motions within the current sheet (CS) associated with a coronal mass ejection (CME) and its accompanying flare. Similarly, Ruan et al. (2020) combined an MHD model that includes anisotropic thermal conduction, radiative cooling, gravitational stratification, and anomalous resistivity with an analytical treatment of fast electrons to study magnetic reconnection, particle acceleration, and X-ray emission in flares. The role of plasmoids and turbulence within CSs in shaping plasma dynamics and heating has also been analyzed by Ye et al. (2020). In a complementary study, Ruan et al. (2021) reproduced post-flare coronal rain as a consequence of thermal instability in post-flare loops. More recently, González-Servín and González-Avilés (2024) investigated the formation and evolution of solar flares in a localized region of the solar atmosphere using resistive MHD simulations, finding that localized resistivity produces substructure around post-flare loops that may be associated with the Richtmyer–Meshkov instability (RMI). In another 2.5D study, Izquierdo-Guzmán et al. (2023) investigated magnetic reconnection in solar flares, focusing on the effects of magnetic-field strength and thermal conduction at the chromosphere-corona interface. Their results showed that the magnetic-field strength has a strong impact on the reconnection rate, producing variations of up to 100%, whereas thermal conduction plays a more modest role, with fluctuations of about 10%.

State-of-the-art three-dimensional (3D) simulations of solar eruptions can capture much of the complexity of plasma dynamics in the solar corona. For example, Jiang et al. (2021) performed ultra-high-resolution 3D MHD simulations showing that simple shearing of a bipolar arcade can generate a current sheet (CS), trigger reconnection, form a flux rope, and ultimately drive a full eruption, including both a flare and a CME, without requiring pre-existing topological structures. Shi et al. (2021) demonstrated that kink waves can contribute to coronal heating through the Kelvin-Helmholtz instability (KHI). In addition, Terradas and Neukirch (2023) developed a method to construct 3D magnetohydrostatic (MHS) models of solar active regions that include the magnetic field, pressure, and gravity, thereby providing a more realistic representation of the lower atmosphere. Likewise, Ye et al. (2023) investigated the thermodynamic evolution of confined turbulence in a large-scale CME-flare current sheet in 3D, concluding that the combined effects of turbulence, shocks, and anisotropic thermal conduction can account for the high-temperature observational signatures of these structures. More recently, Li et al. (2025) studied electron acceleration and transport during solar flares using 3D MHD simulations that include various acceleration sites such as reconnection CSs, termination shocks (TSs), and supra-arcade downflows (SADs) to explain the observed populations of energetic electrons. Testing numerical models against observations in solar flares requires synthetic diagnostics that can be compared with space-based instruments. Facilities such as AIA/SDO (Lemen et al., 2012; Pesnell et al., 2012), UVCS/SoHO (Kohl et al., 1995), EIS/Hinode (Culhane et al., 2007), and IRIS (De Pontieu et al., 2014) provide key observational constraints on the thermal and dynamic structure of flare plasma. Recent studies have therefore used 3D MHD simulations to generate synthetic observables for direct comparison with these instruments, successfully reproducing features such as multi-wavelength coronal emission, supra-arcade downflows, turbulent reconnection regions, and hot loop-top plasma (Lionello et al., 2008; Xie et al., 2019; Ruan et al., 2023; Shen et al., 2023; Wang et al., 2023; Gou and Reeves, 2025; Xie et al., 2025).

Taken together, these studies demonstrate that numerical MHD models can reproduce a wide range of physical and observational properties of solar flares, from magnetic reconnection and energy transport to post-flare loop formation and synthetic emission. However, the physical interpretation of these simulations depends not only on the physical processes included in the model but also on the numerical methods used to solve the MHD equations. In particular, numerical dissipation introduced by the adopted Riemann solver can influence the conversion and partition of magnetic energy, while radiative losses, coronal heating, and thermal conduction regulate the thermodynamic response of the flare plasma. Therefore, a controlled analysis in which the numerical scheme and the main dissipative processes are varied independently is important for distinguishing numerical effects from physical effects and for assessing their influence on synthetic observables.

Motivated by these considerations, in this paper we investigate the formation and evolution of solar flares, with particular emphasis on post-flare loops, using 3D non-ideal MHD simulations that include magnetic resistivity, optically thin radiative losses, anisotropic thermal conduction along magnetic field lines, and both uniform and spatially dependent coronal heating terms. The objectives of this paper are threefold: (i) to assess the influence of the numerical scheme, particularly the choice of Riemann solver, on energy dissipation during the formation of post-flare loops; (ii) to examine the effects of optically thin radiative losses, coronal heating, and thermal conduction on the thermodynamic evolution and energy balance of post-flare loops; and (iii) to compute synthetic EUV emission maps and their corresponding light curves. By addressing these aspects within the same numerical framework, we aim to distinguish solver-dependent effects from those associated with the physical dissipative processes and to evaluate how both influence the resulting synthetic EUV signatures. Overall, this study provides a controlled comparative analysis of numerical and physical effects in an idealized 3D flare configuration, rather than an event-specific reconstruction or an exhaustive numerical benchmark.

The paper is organized as follows. Section 2 describes the system of MHD equations, the model solar atmosphere, the magnetic field configuration, the dissipative terms, the numerical methods, and the simulation cases. Section 3 presents the numerical results, with emphasis on the effects of the Riemann solvers and dissipative terms. Section 4 describes the generation of synthetic emission maps from the MHD simulations, together with the computation of synthetic light curves obtained by integrating the emission in the SDO/AIA channels. Finally, Section 5 summarizes the main conclusions.

2 Model and numerical methods

2.1 The system of MHD equations

To model the plasma in a localized region of the solar atmosphere, we solve the non-ideal MHD equations, including magnetic resistivity, anisotropic thermal conduction, optically thin radiative losses, and coronal heating. In conservative form, the system is written as follows: where ρ, v, and B denote the mass density, velocity, and magnetic field, respectively. The total pressure is defined as p t = p + B 2 / 2, where p is the thermal pressure and B 2 / 2 is the magnetic pressure. Here, I is the identity matrix. The quantities η, g, J, q, Q rad, and H represent the magnetic resistivity, gravitational acceleration, electric current density, heat flux vector, optically thin radiative losses, and coronal heating, respectively. The gravitational acceleration appearing in Equations 2, 3 is assumed to be constant and uniform and is given by g = [ 0 , − g , 0 ], with g = 2.74 × 1 0 4 cm s − 2.

∂ ρ ∂ t + ∇ ⋅ ρ v = 0 , (1)

∂ ρ v ∂ t + ∇ ⋅ ρ v v − B B + p t I = ρ g , (2)

∂ E ∂ t + ∇ ⋅ E + p t v − B v ⋅ B = ρ v ⋅ g − ∇ ⋅ η J × B + ∇ ⋅ q + Q rad + H , (3)

∂ B ∂ t + ∇ ⋅ v B − B v = − ∇ × η J , (4)

∇ ⋅ B = 0 , (5)

The total energy density, E, is given bywhere the three terms correspond to the internal, kinetic, and magnetic energy densities, respectively, and γ = 5 / 3 is the adiabatic index. Throughout this paper, we use the term internal energy to refer specifically to the energy component p / ( γ − 1 ) defined in Equation 6. The term thermal is retained only in established physical expressions such as thermal conduction, thermal pressure, and thermal structure.

E = p γ − 1 + ρ v 2 2 + B 2 2 , (6)

The system of Equations 1–5 is closed using the ideal gas law, given by Equation 7,where T is the temperature, k B is the Boltzmann constant, and m = μ m H is the mean particle mass. For the solar corona, we assume a fully ionized hydrogen–helium plasma and adopt μ = 0.6, where m H is the mass of the hydrogen atom.

p = k B m ρ T , (7)

2.1.1 Magnetic resistivity

To trigger magnetic reconnection and promote the formation of solar flares and post-flare loops, we adopt a localized resistivity profile given bywhere η 0 is the resistivity amplitude and r = x 2 + ( y − h p ) 2 + z 2. The Gaussian profile is centered at a height h p = 6.0 Mm and has a characteristic width w η = 1.0 Mm.

η x , y , z = η 0 ⁡ exp − r / w η 2 , (8)

2.1.2 Anisotropic thermal conduction

Anisotropic thermal conduction is included through the term ∇ ⋅ q in Equation 3. The heat flux is decomposed into components parallel and perpendicular to the magnetic field (Priest, 2014):where the subscripts ‖ and ⊥ denote conduction along and across the magnetic field, respectively, and b ̂ = B / | B | is the unit vector along the field direction. Following the Spitzer formulation (Spitzer, 1962), the parallel thermal conductivity is given by κ ‖ ( T ) = κ 0 T 5 / 2. For coronal conditions, we adopt κ 0 ≈ 9 × 1 0 − 7 e r g s − 1 c m − 1 K − 7 / 2. Since cross-field conduction is much weaker than field-aligned conduction in the solar corona, we neglect the perpendicular contribution, i.e., κ ⊥ ≪ κ ‖.

q = − κ ‖ T b ̂ b ̂ ⋅ ∇ T − κ ⊥ T ∇ T − b ̂ b ̂ ⋅ ∇ T , (9)

2.1.3 Optically thin radiative losses and coronal heating

The term Q rad in Equation 3 represents optically thin radiative losses and is expressed aswhere n e is the electron number density, n H is the hydrogen number density, and Λ ( T ) is the radiative loss function. For simplicity, we approximate the fully ionized coronal plasma as hydrogen-dominated and take n e ≃ n H. The radiative loss function, Λ ( T ), was computed using CHIANTI version 7 (Landi et al., 2011), whereas the AIA temperature response functions used for the synthetic emission maps were computed with SolarSoft/SSW using CHIANTI v9.0.1.

Q rad = − n e n H Λ T , (10)

In addition, we include an empirical coronal heating term, H, in Equation 3. Two forms of coronal heating are considered. The first is a spatially and temporally uniform heating term applied throughout the computational domain, following Shi et al. (2021):where n e = 1 0 8 cm − 3 and Λ 0 = 2.6 × 1 0 − 22 erg cm 3 s − 1. The second heating term is height-dependent and follows Yang et al. (2013):where F = 3 × 1 0 5 erg cm − 2 s − 1 is the energy flux, λ = 4 × 1 0 8 cm is the damping length scale, and y chr = 3 × 1 0 8 cm denotes the height of the chromosphere.

H 1 = n e 2 Λ 0 , (11)

H 2 = F λ , if y ≤ y chr , F λ exp − y − y chr / λ , if y > y chr , (12)

These two forms of coronal heating were selected to investigate the influence of different spatial distributions of coronal heating on the evolution of the post-flare loops. The uniform heating term, H 1, provides a reference configuration in which the energy input is spatially constant, allowing the effects of a homogeneous background heating to be isolated. In contrast, H 2 adopts a height-dependent profile, remaining constant below the chromosphere and decreasing exponentially in the corona with a characteristic scale length λ. This parametrization is commonly employed in coronal MHD models to maintain a background corona against conductive, radiative, and solar wind losses (Hansteen and Leer, 1995; Aiouaz et al., 2005; Yang et al., 2013; Shi et al., 2021).

2.2 Model of solar atmosphere and magnetic configuration

To represent the density stratification between the chromosphere and the corona, we adopt an idealized atmosphere with an analytical prescription similar to that used by Shibata et al. (2023). The density profile is defined as a smooth transition between chromospheric and coronal values:where h T R = 1 M m and w T R = 0.2 M m define the height and width of the transition region, respectively. Equation 13 provides a continuous stratification that captures the sharp density variation between the chromosphere and the low corona. The corresponding temperature profile is obtained from the ideal gas equation of state, and the gas pressure is initialized as p = 1 dyn cm − 2 throughout the computational domain. Together, these profiles define a numerical domain that includes the chromosphere, transition region, and corona, as illustrated in Figure 1.

ρ y = ρ c h r + ρ c o r − ρ c h r 1 2 tanh y − h T R w T R + 1 , (13)

The initial magnetic-field configuration adopted in this work follows Yokoyama and Shibata (2001); Ruan et al. (2020); Izquierdo-Guzmán et al. (2023); Shibata et al. (2023):

B x = 0 , (14)

B y = B 0 ⁡ tanh x / w , (15)

B z = B 0 / cosh x / w . (16)

The initial magnetic field defined by Equations 14–16 depends only on x and therefore represents a translationally invariant, 2.5D-like force-free current-sheet configuration embedded in the three-dimensional computational domain. The subsequent evolution, however, is not constrained to remain invariant along the z direction. In particular, the localized resistivity prescribed by Equation (Equation 8) depends explicitly on x, y, and z, producing a non-ideal region with a finite extent in all three spatial directions. The localized resistivity, therefore, breaks the initial translational symmetry and allows the plasma and magnetic field to develop a three-dimensional structure during the evolution. Thus, the present model should be regarded as a three-dimensional resistive-MHD evolution initialized from a 2.5D-like force-free magnetic configuration.

Besides, this configuration represents a force-free magnetic field that forms a vertical current sheet (CS). Here, B 0 = 33.63 G is the characteristic magnetic-field strength, and w = 5 × 1 0 7 cm (500 km) defines the width of the initial CS magnetic transition. These values are adopted from Izquierdo-Guzmán et al. (2023) to maintain consistency with our previous flare-reconnection model and are kept fixed in all simulations presented here. The parameter w should be regarded as an idealized model scale for the initial current sheet rather than a direct observational estimate of the width of a flare current sheet.

The adopted magnetic configuration is a force-free Harris-type current sheet. In particular, B y 2 + B z 2 = B 0 2 and ∇ × B = [ sech ( x / w ) / w ] B, such that the electric current is parallel to the magnetic field and J × B = 0. Thus, in the ideal-MHD limit, the magnetic configuration represents a force-free equilibrium. Force-free current sheets of this type can be susceptible to resistive tearing modes (Landi et al., 2015; Del Zanna et al., 2016); however, the present simulations are not designed to investigate the spontaneous linear instability of the initial current sheet. No initial velocity or magnetic perturbation is imposed, and no linear tearing-mode growth rate is measured. Instead, reconnection is initiated in a controlled manner by the localized finite resistivity prescribed in Equation (Equation 8). Because this resistivity is present from the beginning of the simulation, the initial force-free magnetic configuration is no longer an equilibrium of the resistive induction equation. With v = 0 initially, the term − ∇ × ( η J ) immediately produces localized magnetic diffusion, breaks the frozen-in condition within the prescribed non-ideal region, and initiates reconnection. Therefore, the onset of reconnection in the present model is imposed through localized resistivity rather than attributed to the spontaneous development of a tearing instability.

2.3 Numerical methods

We solved equations (Equations 1–5) using the PLUTO code (Mignone et al., 2007). All simulations used a second-order Runge–Kutta time integrator, piecewise-linear reconstruction with the minmod limiter, and a Courant–Friedrichs–Lewy (CFL) number of 0.4. To control numerical errors associated with the solenoidal constraint on the magnetic field, Equation 5, we used the extended generalized Lagrange multiplier (GLM) divergence-cleaning formulation (Dedner et al., 2002). The diffusive source terms in the total energy equation, namely magnetic resistivity and thermal conduction in Equation 3, were integrated using a super-time-stepping scheme, allowing for an efficient and stable treatment of parabolic operators. To assess how the numerical scheme affects energy dissipation during post-flare loop formation, we considered five approximate Riemann solvers: Harten–Lax–van Leer (HLL; Harten et al., 1983), Harten–Lax–van Leer–Einfeldt (HLLE; Einfeldt, 1988), Harten–Lax–van Leer–Contact (HLLC; Li, 2005), Harten–Lax–van Leer–Discontinuities (HLLD; Miyoshi and Kusano, 2005), and Roe (Roe, 1981). These solvers provide a useful framework for evaluating the level of numerical dissipation in the MHD simulations.

All simulations were performed in a three-dimensional Cartesian domain spanning [ − 10,10 ] Mm in the x-direction, [0, 20] Mm in the y-direction, and [0, 20] Mm in the z-direction, using a base grid of 400 × 400 × 50 cells. Outflow boundary conditions were imposed on all variables at x = ± 10 Mm, y = 20 Mm, and z = 0 and 20 Mm. These conditions allow plasma and MHD disturbances to leave the computational domain while minimizing artificial reflections at the boundaries. At the lower boundary, y = 0 Mm, a line-tied condition was applied: the plasma velocity was fixed to zero, the normal magnetic-field component, B y, was treated antisymmetrically, while the density, pressure, and tangential magnetic-field component, B z, were assigned outflow boundary conditions. The governing equations were solved in normalized units using the characteristic length L 0 = 1 0 8 cm and velocity V 0 = 1 0 8 cm s − 1, yielding a characteristic time scale of t 0 = L 0 / V 0 = 1 s. Consequently, the simulation time is directly expressed in seconds. Each simulation was evolved up to a physical time of t = 50 s.

2.4 Simulation cases

We consider nine simulation cases, listed in Table 1. These cases explore different combinations of five Riemann solvers, two coronal heating functions (H 1 and H 2, given by Equations 11, 12, respectively), optically thin radiative losses (Q rad, given by Equation 10), and anisotropic thermal conduction, parameterized through κ 0. In all cases, a spatially dependent magnetic resistivity is included with η 0 = 1 0 16 cm 2 s − 1, while the initial magnetic-field strength is fixed at B 0 = 33.63 G. For the adopted initial conditions, p = 1 dyn cm − 2 and B 0 = 33.63 G, the initial plasma beta, β = 8 π p / B 0 2, is approximately 0.022. This value confirms that the magnetic pressure strongly dominates over the thermal pressure, placing the simulations in the low-plasma-β regime characteristic of the solar corona.

TABLE 1. Summary of the numerical solvers and dissipative terms considered in each simulation.

3 Results of the numerical simulations

We next present representative numerical results. We first examine how different Riemann solvers affect energy dissipation during solar-flare evolution and post-flare loop formation. To verify that the reported differences are not dominated by numerical resolution, we performed a resolution test using three grid resolutions and assessed the robustness of the solutions. We then investigate the effects of optically thin radiative losses, coronal heating, and anisotropic thermal conduction. Finally, we present synthetic EUV intensity maps and light curves derived from the simulations.

3.1 Effects of the Riemann solvers

To assess the influence of the Riemann solver on the evolution of post-flare structures, we consider the five solvers listed in Table 1 and focus on Cases #1-#5. In these cases, thermal conduction, coronal heating, and optically thin radiative losses are neglected. Thus, localized magnetic resistivity, defined in equation (Equation 8), is the only explicit dissipative process included in the model.

To complement the analysis, we performed a quantitative comparison by computing the relative deviations in the temporal evolution of the magnetic, internal, and kinetic energy densities. The relative deviation is defined aswhere E i denotes the energy density obtained with a given numerical solver, and E ref corresponds to the energy density obtained with the Roe solver. Since no analytical solution is available for the resistive MHD problem considered here, the Roe solver is used as a numerical reference rather than as an exact benchmark. Its selection is motivated by its relatively low numerical diffusivity, sharp resolution of discontinuities, and well-established performance in computational MHD. In the present study, the numerical diffusion of the Roe solver is not quantified independently; instead, all comparisons are performed relative to the Roe solution under the same numerical setup. Accordingly, the relative deviations presented in Figure 2 quantify the differences between the HLLD, HLLC, HLL, and HLLE solutions with respect to the Roe reference, rather than providing an absolute measure of numerical diffusion. Therefore, the results presented in this work should be interpreted as relative comparisons against the Roe reference solution rather than as measures of agreement with an exact solution. The relative deviations are computed from the time series of the magnetic, internal, and kinetic energy densities.

E rel = E i − E ref E ref ,

Figure 2 shows the temporal evolution of the relative deviations in the magnetic, internal, and kinetic energy densities. The comparison reveals clear solver-dependent differences among Roe–HLLD, Roe–HLLC, Roe–HLL, and Roe–HLLE. For the magnetic energy density, Roe–HLL and Roe–HLLE show a pronounced increase in relative deviation between t ≈ 11 s and t ≈ 22 s, reaching their maximum values during this interval. This is followed by a slight decrease and a subsequent rise toward secondary peaks. Roe–HLLC exhibits early-time fluctuations and reaches its maximum near t ≈ 50 s, whereas Roe–HLLD follows a similar trend but with smaller variations and decreases after t ≈ 25 s.

For the internal energy, the relative deviations increase monotonically for all solvers. Roe–HLL and Roe–HLLE evolve almost identically; Roe–HLLC rises more steeply and approaches the Roe–HLLE values at late times, while Roe–HLLD remains systematically lower throughout the simulation. For the kinetic energy density, the relative deviations show early-time fluctuations, with particularly sharp peaks for Roe–HLL and Roe–HLLE. These two cases later stabilize and fall below the other schemes.

For the kinetic energy density, the relative deviations show pronounced fluctuations during the early stages of the simulation. The largest initial peaks are observed for Roe–HLLC, Roe–HLL, and Roe–HLLE. Roe–HLL and Roe–HLLE exhibit a second local maximum at approximately t ≈ 16 s, after which their relative deviations decrease substantially. Roe–HLLC undergoes several early-time fluctuations, followed by a gradual increase after approximately t ≈ 28 s. In contrast, Roe–HLLD shows less pronounced early-time fluctuations and reaches its first maximum relative deviation at approximately t ≈ 12 s. After t ≈ 25 s, however, its relative deviation increases steadily, eventually becoming the largest among the four solvers by the end of the simulation.

Overall, HLLD (Case #2) shows the smallest relative deviations from the Roe reference solution and the least diffusive behavior in reproducing the evolution and redistribution of the magnetic, internal, and kinetic energy densities. In contrast, HLL (Case #4) and HLLE (Case #5) exhibit the largest relative deviations, indicating the highest levels of numerical diffusion among the solvers considered. HLLC (Case #3) shows intermediate behavior, with more pronounced early-time fluctuations than HLLD but generally smaller relative deviations than HLL and HLLE throughout most of the simulation.

Table 2 summarizes the computational efficiency of the five Riemann solvers corresponding to Cases #1–#5 listed in Table 1. All measurements were obtained using 16 computational nodes. Although the Roe solver provides the least diffusive solution in the present study, it is also one of the most computationally expensive solvers. Roe and HLLD require an average of approximately 12.2 s per time step, whereas HLLC, HLL, and HLLE require only about 3.5–4.0 s per step, while all solvers exhibit nearly identical memory consumption (∼ 330 MB). Therefore, the choice of a Riemann solver represents a trade-off between numerical accuracy and computational efficiency. For the present configuration, HLLD achieves a level of accuracy comparable to Roe while retaining the robustness characteristic of HLL-type solvers, although it does not provide a significant reduction in computational time. In contrast, HLLC, HLL, and HLLE reduce the computational time by approximately a factor of three, making them attractive for large-scale parameter studies or long-duration simulations where computational efficiency is prioritized over the accurate resolution of fine-scale magnetic structures.

TABLE 2. Computational performance of the five Riemann solvers. The table reports the total execution time, average computational time per time step, and peak memory consumption for each simulation.

We further quantify the partition of the total energy density into magnetic, internal, and kinetic components for each numerical solver. This provides a useful diagnostic of the energy-conversion processes associated with magnetic reconnection. At each time step, the total energy density is defined as the sum of these three components, and their relative contributions are evaluated at the final simulation time.

For the Roe solver (Case #1), the final energy partition is 50.05% magnetic energy, 43.54% internal energy, and 6.40% kinetic energy. HLLD (Case #2) yields a very similar distribution, with 49.05% magnetic energy, 44.93% internal energy, and 6.02% kinetic energy, indicating good preservation of magnetic energy and a balanced redistribution among the three components. HLL (Case #4) and HLLE (Case #5) show similar behavior, although with slightly higher internal energy fractions and lower kinetic contributions. In contrast, HLLC (Case #3) produces the lowest magnetic fraction, 45.94%, and the highest internal energy fraction, 48.69%, indicating a larger redistribution of magnetic energy into internal energy in the final energy partition. Overall, Roe and HLLD provide the most balanced redistribution of magnetic energy into internal energy and kinetic energy components.

To further assess the robustness of the numerical results, we performed simulations at three grid resolutions: 200 × 200 × 25, 400 × 400 × 50, and 800 × 800 × 100, hereafter referred to as Res1, Res2, and Res3, respectively. The lower panel of Figure 3 shows the temporal evolution of the normalized magnetic, internal, and kinetic energy components for these three resolutions. These results correspond to Case #2 in Table 1, which uses the HLLD solver. At early times, the curves are very similar across all resolutions, indicating that the global behavior of the system is already captured at a relatively low resolution. As the system evolves, however, a gradual conversion of magnetic energy into internal and kinetic energy becomes evident, consistent with the onset and development of magnetic reconnection. Although the overall trends remain similar, resolution-dependent differences appear at later times.

In particular, the lowest-resolution case, Res1, retains a larger fraction of magnetic energy, whereas the higher-resolution cases, Res2 and Res3, show a more pronounced decrease in magnetic energy, accompanied by a corresponding increase in internal and kinetic energy. This behavior suggests that the energy-conversion process is less effectively captured at lower resolution. It is also consistent with the larger numerical diffusion associated with the coarser grid, which introduces additional smoothing of the numerical solution and reduces the effective resolution of the sharp gradients associated with the reconnection region. As a result, a smaller fraction of the magnetic energy is converted into internal energy and kinetic energy. In contrast, the higher-resolution simulations reduce numerical diffusion, better preserve localized gradients, and therefore capture the energy-conversion process more accurately, leading to a larger transfer of magnetic energy into internal energy and kinetic energy. Overall, these results highlight the importance of using sufficient spatial resolution to properly describe the energy transfer and dynamical evolution of the system during the post-flare phase.

3.2 Effects of dissipative terms

3.2.1 Optically thin radiative losses and coronal heating

We examine the influence of optically thin radiative losses and coronal heating on energy dissipation in post-flare structures, focusing on Cases #2, #6, and #7 in Table 1. Case #2 is used as the reference case and does not include either coronal heating or radiative losses. Case #6 includes optically thin radiative losses together with the uniform coronal heating model, H 1, whereas Case #7 includes the same radiative losses but uses the height-dependent heating model, H 2. These three cases, therefore, allow the effects of radiative losses and the two heating prescriptions to be compared directly while all other simulation parameters are kept unchanged. The upper-left panel of Figure 3 shows the temporal evolution of the domain-integrated magnetic, internal, and kinetic energies, normalized by the total energy integral, for the three cases. The upper-right panel shows the ratio of internal energy to magnetic energy, E i n t / E m a g, as a diagnostic of the relative balance between the internal and magnetic energy reservoirs. Up to t ≃ 25 s, the curves are nearly indistinguishable, indicating that the system remains magnetically dominated and that the increase in internal energy is primarily driven by the conversion of magnetic energy during magnetic reconnection. Consequently, neither coronal heating nor optically thin radiative losses significantly affect the global energy budget during this stage.

At later times, around t ≃ 40 s, Case #2 begins to diverge from Cases #6 and #7. Case #2 exhibits a stronger decrease in magnetic energy and a larger increase in internal energy, indicating more efficient retention of the released energy as internal energy. In contrast, in Cases #6 and #7, which include coronal heating and optically thin radiative losses, part of the released magnetic energy is transferred to the internal energy reservoir and subsequently removed by radiation. By t = 50 s, the ratio E i n t / E m a g reaches about 13% in Case #2, but only about 9% in Cases #6 and #7, showing that coronal heating and radiative cooling modify the fraction of energy retained as internal energy during the post-flare evolution.

Figure 4 shows the temperature and density profiles along the vertical sampling line passing through the center of the localized reconnection region at x = 0 Mm and z = 0 Mm, at t = 15 and 35 s for Cases #2, #6, and #7. At t = 15 s, the profiles are nearly indistinguishable, indicating that coronal heating and radiative cooling have not yet significantly altered the thermodynamic structure of the plasma. By t = 35 s, clear differences appear between Case #2 and Cases #6 and #7. These differences are most evident across the upper chromosphere, transition region, and low corona, approximately over 2.0 ≲ y ≲ 7.5 Mm. In this region, Case #2 shows smoother temperature and density gradients, whereas Cases #6 and #7 develop steeper profiles. In the cases including coronal heating and radiative cooling, the density decreases more abruptly and the temperature reaches higher values, consistent with enhanced pressure-driven expansion of the coronal plasma. Despite the different heating prescriptions adopted in Cases #6 and #7, their temperature and density profiles remain nearly indistinguishable throughout the simulated time interval, indicating that the specific form of the coronal heating function has only a minor influence on the plasma stratification under the present conditions. Thus, although all three cases exhibit significant thermodynamic evolution by t = 35 s, the inclusion of coronal heating and optically thin radiative losses produces steeper temperature and density gradients than those obtained in the reference Case #2. Overall, these results indicate that coronal heating and optically thin radiative losses modify the thermodynamic structure of the post-flare plasma, whereas the specific spatial distribution of the heating function has only a minor effect over the simulated time interval.

3.2.2 Thermal conduction

We compare Case #6 with Case #8 to investigate the effect of anisotropic thermal, defined by Equation 9, on energy dissipation during post-flare loop formation. Figure 5 shows the results at t = 15 s, with Case #6 in the left column and Case #8 in the right column. The figure displays maps of temperature, in K, mass density, in g cm − 3, and the z-component of the current density, J z, in statA cm − 2, overlaid with magnetic-field lines.

The top-row panels show the temperature distributions for Cases #6 and #8, respectively. In both cases, the magnetic-field lines bend around the heated region and trace the developing post-flare loop system. However, in Case #8, thermal conduction redistributes energy along the field lines and strongly reduces the localized temperature enhancement relative to Case #6. The maximum temperature decreases from 2.4 × 1 0 9 K in Case #6 to 8.1 × 1 0 7 K in Case #8, indicating that thermal conduction reduces the peak temperature by a factor of about 30. This result emphasizes the role of field-aligned heat transport in limiting localized temperature buildup in the present model.

The middle-row panels show the corresponding mass density distributions. In both cases, the density structure follows the magnetic topology. However, in Case #8, the inclusion of anisotropic thermal conduction transports heat along the magnetic-field lines, thereby redistributing internal energy and reducing localized temperature and pressure gradients. As a result, the plasma evolves more smoothly, and the transition between the dense plasma confined within the post-flare arcade and the surrounding coronal plasma becomes more gradual, producing a less localized density distribution than in Case #6. The bottom-row panels show the corresponding distributions of J z. In both cases, enhanced current density is concentrated within the main reconnecting current sheet and around the cusp region of the post-flare arcade. Significant current concentrations also develop along the low-lying post-flare loops. When anisotropic thermal conduction is included, the enhanced current region surrounding the post-flare arcade extends over a larger area, producing a less localized distribution of J z, while the reconnecting current sheet remains clearly defined.

To isolate the effect of anisotropic thermal conduction, we compare Cases #6, #8, and #9 listed in Table 1. These cases share the same Riemann solver, coronal-heating model ( H 1 ), and optically thin radiative losses, differing only in the value of the thermal conductivity coefficient κ 0. Figure 6 shows vertical profiles of temperature, density, and J z as functions of height y at x = 0. The three cases are compared at two representative times, t = 15 and 35 s, to illustrate the evolution of the vertical stratification. At t = 15 s, all three cases show a temperature enhancement near y ≈ 2.5 Mm. However, Case #6 reaches temperatures of about 1 0 9 K, whereas Cases #8 and #9 remain near 1 0 7 K. This indicates that, in the absence of thermal conduction, heat accumulates locally and produces much stronger temperature enhancements. The density profiles show the same trend: the higher temperatures in Case #6 drive stronger plasma expansion and, therefore, lower densities, while Cases #8 and #9 retain higher densities because thermal conduction redistributes the deposited energy more efficiently. At this stage, the J z profiles are broadly similar in all three cases, although Cases #8 and #9 show a localized enhancement near y ≈ 7.5 Mm.

By t = 35 s, the differences become more pronounced. Case #6 develops sharper temperature and density gradients between y ≈ 5 and 15 Mm, reaches higher temperatures than Cases #8 and #9, and shows a stronger J z peak near y ≈ 10 Mm. This behavior is consistent with a more concentrated current sheet in the reconnection region when thermal conduction is absent. In contrast, the cases including thermal conduction exhibit smoother thermal and density structures, lower peak temperatures, and a broader current-density distribution. These results show that anisotropic thermal conduction plays a key role in regulating the thermal response of the post-flare loop system by redistributing energy along magnetic-field lines. They should nevertheless be interpreted in light of the limitations of the adopted conductive prescription and the idealized flare configuration.

4 Synthetic intensity maps

To connect the MHD results with observable EUV emission, we constructed synthetic intensity maps using the temperature response functions of the selected SDO/AIA channels. Since the simulated post-flare plasma reaches temperatures on the order of 1 0 7 K, we focus on the AIA 94 Å and 131 Å channels, which include hot flare-sensitive response components and are commonly used to diagnose high-temperature flare plasma (Shen et al., 2022). The procedure follows Shen et al. (2022); Xie et al. (2019). For a given AIA channel m, the synthetic intensity is estimated aswhere n e i is the electron number density in cm − 3, T e i is the plasma temperature in K, Δ l i is the thickness of the ith cell along the line of sight, and R m ( T ) is the temperature response function of channel m. In practice, the density and temperature were extracted from each cell of the 3D MHD simulations, and the electron number density was computed aswhere ρ sim is the simulated mass density, μ e is the mean molecular weight per electron, and m H is the mass of the hydrogen atom. The instrumental response functions, R m ( T ), were computed using the SolarSoftWare (SSW) package in IDL. ¹ Figure 7 shows the response functions for the AIA 94 and 131 Å channels used in this work.

I m = ∑ i n e i 2 R m T e i Δ l i , (17)

n e = ρ sim μ e m H , (18)

Synthetic emission maps, computed using Equations 17, 18, are presented for Case #8, listed in Table 1, because this case includes the full set of physical processes considered in this study, namely coronal heating, optically thin radiative losses, and anisotropic thermal conduction. It therefore provides the most representative configuration for assessing the observable EUV signatures of the simulated post-flare loops. Figure 8 shows the corresponding synthetic emission maps for the AIA 94 Å and 131 Å channels at t = 15 and 35 s. The morphology of the synthetic emission is qualitatively consistent with the standard flare scenario. Bright emission near the loop-top region indicates strong plasma heating above the arcade. As the flare evolves, downward plasma outflows along newly reconnected magnetic-field lines compress and heat the plasma near the loop apex, producing the cusp-shaped morphology characteristic of post-flare loop systems. The increase in brightness and spatial extent of the emission from t = 15 to 35 s reflects the continued formation of reconnected loops and the accumulation of hot plasma within the growing flare arcade. The cusp and current-sheet structures are more sharply defined in the AIA 131 Å channel because this channel is more sensitive to high-temperature plasma. Consequently, the hottest regions are more clearly identified in the AIA 131 Å channel than in the AIA 94 Å channel.

Figure 9 shows the synthetic light curves and flux ratio for the AIA 94 Å and 131 Å channels. The light curves were obtained by integrating the synthetic emission maps over the full computational domain at each simulation time, yielding the total integrated emissions F 94 and F 131 for the AIA 94 Å and 131 Å channels, respectively. The corresponding flux ratio, F 131 / F 94, is shown in the bottom panel. Throughout the evolution, the emission in the 131 Å channel remains substantially higher than that in the 94 Å channel. The flux ratio ranges from approximately 12–37, indicating that the synthetic emission is dominated by hot flare plasma. This behavior is consistent with the synthetic emission maps, in which the cusp and current-sheet structures are more clearly defined in the AIA 131 Å channel. According to the upper panel of Figure 9, the temporal evolution of the synthetic light curves can be divided into three phases. During the initial phase, both channels exhibit a rapid increase in emission as magnetic reconnection rapidly heats the plasma to temperatures within the response range of the AIA 94 Å and 131 Å channels. This is followed by an overall decrease accompanied by moderate temporal fluctuations, reflecting the continuous evolution of the post-flare arcade. During the final phase, the integrated emission increases rapidly as the bright emitting region expands around the growing post-flare arcade and loop-top region, increasing the volume of plasma contributing to the total emission.

Although the magnetic configuration and thermodynamic conditions in the simulation are not intended to reproduce a specific observed flare, the synthetic emission maps recover the main morphological features commonly observed in post-flare loop systems. At the same time, the integrated emission is systematically higher in the AIA 131 Å channel than in the AIA 94 Å channel, with flux ratios ranging from approximately 12–37. These results indicate that the synthetic emission is dominated by hot plasma, consistent with the temperature distribution expected in flare arcades.

Figure 10 compares the synthetic emission maps obtained with the Roe and HLLC Riemann solvers at t = 15 s for the AIA 94 Å and 131 Å channels. These cases correspond to Case #1 and Case #3 of Table 1, respectively. These two cases are selected to illustrate how differences in the numerical treatment of the MHD equations are reflected in the synthetic EUV emission. In both cases, the synthetic emission reproduces the characteristic morphology of a post-flare arcade, including the cusp-shaped structure, the loop-top region, and the elongated current sheet above the arcade. However, differences are evident in the spatial distribution of the emission. The emission maps obtained with the Roe solver preserve a more compact and sharply defined emitting structure around the cusp and current-sheet regions, whereas those obtained with the HLLC solver exhibit a more diffuse emission pattern around the cusp, and the emission associated with the current-sheet region also appears more diffuse.

The Roe solver preserves a more compact and sharply defined emitting structure around the cusp and current-sheet regions, whereas the HLLC solver produces a more spatially extended emission pattern, particularly around the post-flare arcade and loop-top region.

Although the synthetic emission maps obtained with different Riemann solvers show qualitative morphological differences, the effects of numerical dissipation are more clearly seen in the total emitted flux. To quantify this effect, we computed synthetic light curves for the AIA 94 Å and 131 Å channels using the five Riemann solvers considered in this work, corresponding to Cases #1-#5. Figure 11 shows the resulting light curves. Throughout the simulation, the HLLE solver exhibits markedly different behavior from the other solvers. The light curves obtained with Roe, HLLD, HLLC, and HLL remain relatively close to each other and display a similar temporal evolution up to approximately t = 30 s, whereas the HLLE curve departs significantly from this trend. At later times, the HLLD light curve remains systematically below those obtained with Roe and HLLC, while the latter two show a more pronounced increase in the emitted flux.

Overall, the synthetic emission analysis shows that the simulations reproduce the main observational signatures of post-flare loop systems, including a cusp-shaped arcade, enhanced loop-top emission, and an elongated current sheet. The systematically stronger emission in the AIA 131 Å channel, together with flux ratios ranging from approximately 12–37, indicates that the synthetic EUV emission is dominated by high-temperature plasma generated during magnetic reconnection. Although different Riemann solvers produce broadly similar large-scale morphologies, the corresponding light curves reveal significant differences in the total emitted flux. These results show that numerical dissipation affects not only the thermodynamic evolution of the plasma but also the observable EUV signatures derived from the simulations. Therefore, the choice of numerical solver should be carefully considered when performing quantitative comparisons between MHD models and observations.

5 Conclusion

Our results show that the evolution of post-flare loops in three-dimensional MHD simulations is controlled not only by the underlying magnetic processes but also by the numerical treatment used to resolve them. Within the fixed numerical setup considered here, the choice of Riemann solver directly affects the partition and redistribution of magnetic energy, internal energy, and kinetic energy; therefore, it influences the physical interpretation of the modeled flare evolution. Among the solvers examined, HLLD exhibits the least diffusive overall behavior, with energy partitions closest to those of the Roe reference solution. This result indicates that the numerical method plays an active role in the simulated dynamics and in the resulting formation of post-flare structures.

From a physical perspective, radiative losses, coronal heating, and anisotropic thermal conduction substantially modify the thermodynamic evolution of the modeled flare plasma. Radiative losses reduce the fraction of released magnetic energy retained as internal energy and help prevent unrealistically high coronal temperatures, while the adopted heating prescriptions influence the arcade structure at later stages. Most importantly, anisotropic thermal conduction plays a central role in transporting energy along magnetic-field lines, smoothing the temperature and density distributions within the loops and reducing the maximum temperature reached in the reconnection region by a factor of about 30. These results show that the evolution of modeled post-flare loops is strongly affected by the inclusion of the main dissipative processes.

The synthetic EUV intensity maps provide an additional diagnostic of the modeled post-flare structures. Although the simulations are not intended to reproduce a specific observed flare on a one-to-one basis, they recover the main morphological features commonly seen in AIA observations, including cusp-shaped post-flare arcades, enhanced loop-top emission, and the gradual development of hot loop systems. These qualitative similarities are further supported by the synthetic light curves derived from the emission maps. The integrated emission remains systematically higher in the AIA 131 Å channel than in the AIA 94 Å channel, indicating that the synthetic flare arcade is dominated by high-temperature plasma.

The results of this study provide a useful three-dimensional framework for examining how numerical solvers and dissipative processes influence post-flare loop evolution in an idealized non-ideal MHD configuration. More generally, they demonstrate that the physical interpretation of three-dimensional flare simulations depends not only on the dissipative processes included in the model but also on the numerical scheme used to evolve the system. This work provides a foundation for applying validated three-dimensional numerical flare models to real observations. Future studies will build on this framework by incorporating observation-based magnetic-field configurations, enabling more realistic simulations and a closer connection between theoretical predictions and observed flare behavior.

Finally, we note that the conclusions drawn here are restricted to the fixed numerical framework adopted in this study, including the grid resolution, reconstruction scheme, limiter, divergence-control method, and boundary conditions. Accordingly, the ranking of Riemann solvers reported in this work should be interpreted as setup-dependent rather than universal. Although a broader robustness assessment, including systematic resolution studies and sensitivity tests involving alternative numerical choices, is beyond the scope of the present paper, the emergence of consistent solver-dependent trends across several independent diagnostics supports the internal coherence of the comparison. In this sense, the present results provide a physically motivated reference for the adopted configuration, while a more general assessment of solver performance will require future studies spanning different numerical setups and observation-based magnetic-field geometries.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. The authors verify and take full responsibility for the content of this manuscript. Generative AI was used to assist with language editing and improvement of clarity in the text. All scientific content, analysis, and conclusions were developed and validated by the authors.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Abhishek Kumar Srivastava, Indian Institute of Technology (BHU), India

Reviewed by: Alexey A Kuznetsov, Institute of Solar-Terrestrial Physics (RAS), Russia

Satyam Agarwal, The University of Alabama in Huntscille, United States

Footnotes

https://www.lmsal.com/solarsoft/ssw/