Exploring self-consistent 2.5 D flare simulations with MPI-AMRVAC

Context. Multi-dimensional solar flare simulations have not yet included detailed analysis of the lower atmospheric responses such as down-flowing chromospheric compressions and chromospheric evaporation processes. Aims. We present an analysis of multi-dimensional flare simulations, including analysis of chromospheric up-flows and down-flows that provide important groundwork for comparing 1D and multi-dimensional models. Methods. We follow the evolution of an MHD standard solar flare model including electron beams, where localized anomalous resistivity initiates magnetic reconnection. We vary the background magnetic field strength, to produce simulations that cover a large span of observationally reported solar flare strengths. Chromospheric energy fluxes, and energy density maps are used to analyse the transport of energy from the corona to the lower atmosphere, and the resultant evolution of the flare. Quantities traced along 1D field-lines allow for detailed comparison with 1D evaporation models.


Introduction
In this paper we present multi-dimensional simulations of solar flares with a focus on the lower atmospheric dynamics that result from the coronal energy release.
The first solar flare models based on magnetic reconnection Sweet (1958); Petschek (1964) were developed in the mid 1900s (Parker 1963;Carmichael 1964;Sturrock 1966), and led to the so-called standard CSHKP flare model (Sturrock 1992;Shibata 1996).These models explained generalised observed behaviours of solar flares using the release of stored magnetic energy through magnetic reconnection events (see Fig. 1).This cartoon model has been expanded upon by many works, which are summarised well in Priest & Forbes (2002).Subsequent 3D models have stressed the role of current concentrations and quasi-separatrix layers (Aulanier et al. 2012(Aulanier et al. , 2013;;Janvier et al. 2013Janvier et al. , 2014)).
The 2.5D simulations in this work present much simpler magnetic topologies than those in 3D models, but will account for all the relevant thermodynamic processes and the effects of electron beams.Common features of such cross sections through solar flare models, are: 1.A flux rope that runs perpendicular to the 2.5D cross section, and an underlying sheared flare loop arcade with a re-connection site in between (Sturrock 1992;Shibata 1996).In eruptive flares this flux rope is ejected outward (upward).2. A magnetic x-point below the flux rope, where successive field lines are drawn inward, reconnected, and ejected outwards.3. Fast reconnection out-flow jets from the x-point, in the orientation of the current sheet, one of which is directed towards the solar surface (Petschek 1964).These shocks can take a lobster claw form during the initial ejection (Zenitani & Miyoshi 2011).4. Reconnected loops, ejected towards the surface, group together below the x-point and form hot coronal loop arcades (Sturrock 1966;Shibata 1996). 5.The out-flow jets become super-Alfvénic, and establish slow-mode shocks at their edges.A fast-mode termination shock forms at the interface between the core of the outflows and the hot coronal loop arcade below (Shibata 1996;Yokoyama & Shibata 2001;Ruan et al. 2020;Shen et al. 2022;Ruan et al. 2023).6.Energy is liberated from the magnetic field at the reconnection x-point and potentially in a number of other larger volumes throughout the flare loop system, for example, in slow shocks in the reconnection jets (Petschek 1964;Priest & Forbes 2002).The liberated energy is converted into many .Later (at 16:18:07 UT) the flux rope has left the field of view and we see emission that is interpreted to originate from the long bright current sheet (red arrows), and dense coronal flare loops below (green arrows).These images use the short exposure AIA filters, but still also display significant fringe artifacts emanating outward from the bright loop arches.
forms including heat, with models showing a significant fraction can be converted into the acceleration of high energy non-thermal particles (Rowan et al. 2017;Werner et al. 2018;Hoshino 2023).Energy is transported from the neighbourhood of the x-point, and the hot coronal loop-tops, along magnetic field lines and towards chromospheric footpoints near the surface of the Sun.There are many proposed processes for this transport (see review by Zharkova et al. 2011).A key process is energetic particle acceleration (electrons and ions) along the field-lines (Syrovatskii & Shmeleva 1972;Emslie 1978;Holman et al. 2011;Kong et al. 2022) with the acceleration mechanism falling into three broad categories: Direct Current electric field acceleration by an electric field above the Dreicer limit (Dreicer 1959(Dreicer , 1960)), shock acceleration (Ellison & Ramaty 1985), and turbulent/stochastic acceleration (Miller et al. 1996;Cargill et al. 2012;Kontar et al. 2017), including by the Kelvin Helmholtz Instability in turbulent loop tops (Fang et al. 2016a;Ruan et al. 2018).Other energy transport mechanisms that will be present include thermal conduction predominantly parallel to the field-lines (Spitzer 1962;Spitzer & Tomasko 1968;Spitzer & Scott 1969), and Alfvén waves (Fletcher & Hudson 2008).7. Flare ribbons are the chromospheric locations where much of the released magnetic energy is deposited.This energy heats and excites the plasma and producing increased emission.Hard X-ray (HXR) sources are generally most intense in footpoints of the flare-loops.HXR sources are believed to result from bremsstrahlung of non-thermal (high-energy) electrons that lose their energy in collisions with the am-bient thermal plasma.These energetic non-thermal distributions of particles are believed to accelerated near the x-point, and may also be accelerated in the termination shock and the turbulent reconnection that occurs in the tops of the loop arcades or in multi-stage acceleration processes (Holman et al. 2011;Zharkova et al. 2011;Kong et al. 2019Kong et al. , 2020)).8.The energy released in the flare loop foot-points causes hot up-flows (chromospheric evaporation) and cooler downflowing "chromospheric compressions" towards the solar surface.The chromospheric evaporation fills the flare loops with hot dense plasma causing bright emission in the UV lines and soft X-ray (SXR) spectrum (Bruzek 1969;Hirayama 1974;Syrovatskii & Shmeleva 1972;Fisher et al. 1985;Polito et al. 2016;Druett et al. 2017;Polito et al. 2023).
In coronal physics a "condensation" refers to material that is dramatically cooling (with temperatures decreasing by an order of magnitude or so) and, typically, dropping from a higher ionisation state into a lower or partially unionised state.When this occurs the material suddenly becomes visible in chromospheric spectral lines, appearing to "condense"."Chromospheric evaporation" refers to an up-flow of chromospheric material with simultaneous heating (temperature increasing by an order of magnitude or so), typically including a large increase in the ionisation degree of a particular state for an element.Although this "evaporation" is technically an ablation, the material appears to "evaporate" from chromospheric spectral lines.Therefore, we reserve the words "evaporation" and "condensation" for processes that have, at the very least, some semblance of a change in state.We note that the "down-flowing chromospheric compressions" in flares have often been referred to in the literature as "chromospheric condensations".However the conditions for the "condensation" analogy do not hold true for these phenomena.Therefore, we refer to this phenomenon as a "down-flowing chromospheric compression" for the rest of the paper (including quotation marks as a reminder to the reader).
Recently multi-dimensional magnetohydrodynamic (MHD) models have investigated plasma flows in flares.Cheung et al. (2019); Rempel et al. (2023) used sub-photospheric velocity driving to build up and release energy in the corona.Kong et al. (2019Kong et al. ( , 2020Kong et al. ( , 2022)); Shen et al. (2022) have reproduced and interpreted supra-arcade down-flows in the corona immediately above the flare loops, with Kong et al. (2020Kong et al. ( , 2022) ) inspecting the energetic electron acceleration and transport without yet transferring energy out of, and back into the MHD simulation.Kerr et al. (2020) used a set of 1D loop models to build up a 3D flaring volume.However, the main focus of multi-dimensional flare simulations has been on coronal dynamics.The investigation presented in this work includes a significant focus on the lower atmospheric dynamics and provide methods that, for the first time, enable clear comparisons of results from 1D and multidimensional modelling.
There have also been attempts to reproduce the standard model in multi-dimensional MHD models.Yokoyama & Shibata (2001) initialised their 2D model as a weak bipolar magnetic field region within a vertically stratified approximation of solar conditions.They used anomalous resistivity to trigger reconnection above the polarity inversion line, which resulted in a loop arcade forming beneath.The combination of the reconnection out-flows impacting on the lower atmosphere and the thermal conduction front drive chromospheric evaporation from the footpoints of the loops.The resultant dense coronal loop arcade matches the general evolution scheme of a solar flare.Ruan et al. (2020) self-consistently built on the model of Yokoyama & Shibata (2001) in 2.5D using the message passing interface-adaptive mesh refinement versatile advection code (MPI-AMRVAC, Keppens et al. 2012;Porth et al. 2014;Xia et al. 2018;Keppens et al. 2023), and expanded it significantly by using the Ohmic heating term in selected regions of anomalous resistivity as an energy reservoir to accelerate non-thermal electrons along field-lines.This energy is redistributed along these field-lines using analytical solutions for the 1D thick-target modelling (Emslie 1978) with remotely deposited energy subsequently re-interpolated onto the automated, block-adaptive grid.However, the agents causing chromospheric evaporation (hot upflows of plasma with chromospheric densities into the corona) in the original paper are thermal conduction and the impact of the out-flows on the lower atmosphere, transported down the flare loop arcade.
A companion paper to this study presents the first selfconsistent multi-dimensional model of this kind reproducing chromospheric evaporation via energetic particle beams (Druett et al. 2023).Ruan et al. (2023) presents a simulated flare in 3D using this modelling suite, but without including beam electrons, to study the formation of MHD turbulence in the flare loop-tops.
In this paper we explore the 2.5D models including beam electrons described in Ruan et al. (2020), in particular how variations of the coronal field strength affect the resultant coronal and lower atmospheric dynamics.This investigation provides the first solid basis for the comparison between 1D radiation hydrodynamic flare simulations and multi-dimensional flare modelling results.Thereby, we also lay groundwork to assess the need for the inclusion of critical field-aligned 1D physics to be built into truly multi-dimensional flare models.

MHD model description
The setup of the models in these experiment is comprehensively described in Ruan et al. (2020).Here we only recall the equations used, along with an overview of how this simulation reproduces features of the standard solar flare model.The MHD equations are, The equations are written in a dimensionless format.ρ, v, t, B, e are the plasma density, velocity, time, magnetic field, and energy density.g is the gravitational acceleration, which acts vertically downwards.This is calculated via the equation g = −274R 2 s /(R s + y) 2 ŷ m s −2 , where R s is the solar radius.J is the current density defined by J = ∇ × B, and η is the resistivity, with anomalous forms described in Ruan et al. (2020) and section 2.2.
Equation ( 1) is the continuity equation, expressing the conservation of mass.Equation ( 2) is the equation of motion also writable as, The Lorentz force J × B has been brought to the left hand side and decomposed into the magnetic tension and pressure making the total pressure p tot = p gas + (B 2 /2) = p gas + p mag .The form of these equations is discussed at length in section 4.3 of Goedbloed et al. (2019).Equation ( 3) is the energy equation, where the total energy density is e = (ρv 2 )/2 + p gas /(γ − 1) + p mag with γ = 5/3.The first term on the right side of equation (3) represents gravitational potential energy and additional sink or source terms on the right side express heat conduction (with a thermal conductivity tensor κ), resistive effects (The term ∇ • (B × ηJ) results from the inclusion of Ohmic dissipation, but is not Ohmic heating.Instead it shows that Ohm's law specifies the comoving electric field to  d) shows the typical flare loop system that is formed via this process, through a plot of absolute magnetic field strength.Energetic electron transport is switched on after 31.2 s of the simulation.Thus, in all cases presented in this manuscript, the electrons are switched on before the impact of the reconnection jet on the lower atmosphere.On the lower panels, electron acceleration sites are shown in green, with energy deposition locations displayed in yellow.The energy deposition locations are saturated at very low values to help indicate their paths through the experiment.In fact, these beams deposit the majority of their energy in the lower atmosphere at the footpoints of the flare loops.
be ηJ, and that total energy remains conserved in resistive MHD, see section 4.4.2 of Goedbloed et al. (2019)), optically thin radiative losses (Q r ) and an artificial background heating that maintains the quiet sun coronal temperature (see equation 5 of Ruan et al. (2020)).
In the second phase of the resistivity description, we take the Joule heating term |ηJ 2 | out of the local energy equation and use it as a reservoir of energy for the acceleration of energetic electrons (see section 2.3).The term Q e represents this energy that is lost from the sites of energetic electron acceleration.The term H e represents the heating of the plasma by these energetic electrons, which is non-local and transferred along the magnetic field lines.
Equation (4) shows the induction equation, which governs the advection of the magnetic field with the plasma.A source term on the right side describes the effects of resistive field dif-fusion, misaligned currents, and resistivity gradients (see section 4.4.2 and equation 4.132 of Goedbloed et al. (2019)).Coupled with the energy equation described above, this acts to convert magnetic energy into internal energy at sites of resistivity.The system of equations is closed by an ideal gas law as the equation of state.

Anomalous resistivity description
The anomalous resistivity prescription in these simulations is a two-stage model.Both are calculated as described in Ruan et al. (2020) and were based on the earlier model of Yokoyama & Shibata (2001).The first stage is used to trigger the reconnection at the x-point in the corona and takes the form η 0 = 0.05 is the maximum value of the anomalous resistivity in this stage.r η = 2.4 Mm is the radius over which the anomalous resistivity drops to zero, with distance r away from the point (0, 50)Mm.The second phase is activated at times after t η > 31.2 s, which is t η = 0.4 units of experiment time, This generates anomalous resistivity in locations where the calculated electron drift velocity, v d (x, y, t) = J eN e , is greater than a critical value v c = 1000u v , where u v = 128 km s −1 , e is the electron charge, and α = 1 × 10 −4 .The resistivity produced is given a maximal value of 0.1.It is greatest at heights near the original resistivity patch, h η = 50 Mm, and decreases with a scale height of h s = 10 Mm.
Ohmic heating results from the combination of the second stage anomalous resistivity description and the current in the current sheet through the domain about x = 0.This Ohmic heating term, η|J| 2 , is used to model unresolved microscopic instabilities, and in the second stage of our resistivity description is the energy reservoir Q e used to accelerate energetic electrons.
Any number of other resistivity schemes are available.These different schemes will impact the MHD evolution of the system.Investigation of large deviations from the scheme above is outside the scope of this work.It will be investigated in other works, including Druett et al. (2023).

Beam model description
The beam electron modelling used in this work is unchanged from the description of Ruan et al. (2020), who describe the approach as a generalisation of the 1D treatment provided by Emslie (1978).
Firstly, we define the flaring region by tracing magnetic fieldlines forward and backward from points to the left and the right of the reconnection x-point, at (±2.5, 50) Mm.Within this flaring region field lines are traced from their photospheric footpoints.A check is made that neighbouring field-lines do not separate from each other by more than the separation of the grid-cell centres.If they do, additional field lines are added.The points where these field lines pass closest to the line x = 0 are defined as the starting points for the subsequent 1D energy redistribution by beam electrons.
The input energy flux along each field-line is calculated as the total of the Ohmic heating terms of cells that intersect with each field-line, as described in Equation 8.The energy transferred to the field lines is subtracted from the MHD energy equation as Q e .Where multiple field lines intersect a cell that produces energetic particles, the energy from that cell is shared equally between the field-lines.
Our numerical treatment of purely resistive MHD is fully conservative by construction as we use a finite volume approach, but partially open boundaries can modify exact conservation.Moreover, the additional presence of electron beams imply that at any instant, the beam injected energy is stored on (evolving) field lines, and given back to the plasma by interpolation from field lines to grid cells at a remote location.A consideration of how overall conservation is approximately achieved for the beam electrons at all times is given in Appendix A.

Domain and solution methods
The equations in section 2.1 are solved in a spatial domain spanning −50 Mm< x < 50 Mm and 0 Mm< y < 100 Mm, using the open-source MPI-AMRVAC code (Keppens et al. 2012;Porth et al. 2014;Xia et al. 2018;Keppens et al. 2023).The hierarchical, block-adaptive grid used has a block size of 16 by 16 cells, with a minimum of 64 cells (4 by 4 blocks) spanning the domain in each dimension, at refinement level 1.The grid is refined by splitting a block into 4 sub-blocks for each increase of refinement level.The maximum refinement level for the experiment is 6.This means that at lowest refinement the grid cell separation is 100/64 = 1.5625Mm, and at maximum refinement the resolution is 100/64/2 5 = 48.8km.
The refinement level is forced to be maximal below y = 3 Mm, and for blocks within the box containing the dynamically tracked magnetic field-lines (with the regions as described in Appendix B of Ruan et al. 2020).Additional automatic refinement and de-refinement is switched, on to ensure accurate shock capturing in locations away from the user-prescribed refinement areas.This is implemented with a weighting of 1:2:2 (as described in Keppens et al. 2012) between the conserved variables of mass density, vertical magnetic field, and internal energy density respectively.
We employ a three step time-stepping scheme.The flux scheme uses the HLL approximate Riemann solver (Harten et al. 1983) and as in Ruan et al. (2020), a mixture of high-order slope limiters is used: a third-order limiter ( Čada & Torrilhon 2009) is employed in regions of low refinement, i.e. the background corona, and a second-order limiter (van Leer 1974) in the regions of high grid refinement (greater than level 3), namely the lower atmosphere, reconnection region, and flare loop.The various limiters and all solvers available are discussed in Keppens et al. (2023).

Differences to previous studies
There a few differences between the experiments presented here and that produced in Ruan et al. (2020), namely, -The heat saturation model parallel to the field lines has been enacted using the MPI-AMRVAC thermal conduction module Xia et al. (2018), with a monotonized central limiter as per Woodward & Colella (1984).-The side boundaries of the experiment (x-direction) are now treated as open, rather than the previously employed combination symmetric and asymmetric boundary conditions in the ghost zones.This avoids reflection of shocks that emanate from the flare, and these shocks simply exit the domain from the sides of the experiment.In Ruan et al. (2020) shocks reaching the side boundaries were reflected and returned and interact with the flaring region.As a result we now allow (negligible) mass loss in the chromosphere and the corona via the open boundaries over the duration of the experiment.
The upper and lower boundaries are unchanged from Ruan et al. (2020).-The magnetic field vectors are split.There is a constant background component with a distribution as in the model of Yokoyama & Shibata (2001), and we solve for a time varying component that is zero at the start of the experiment.In Ruan et al. (2020) the background part of the field was given a magnitude of B 0 = 35G.In this investigation, 4 different values are used (B 0 = 20G, 35G, 50G, 65G) to explore the impact of different coronal field strength on the flare evolution.
2.6.How this experiment reproduces features of the standard solar flare model Figure 2 shows how these experiments reproduce features of the standard solar flare model.The experiment is initialised with a low current-density vertical current sheet in the centre, where the vertical background magnetic field components transition from positive (left side of the experiment) to negative values (right side).This bipolar field region undergoes magnetic reconnection due to the anomalous resistivity region inserted at a height of 50 Mm.The current sheet thins and grows stronger in the reconnecting locations (Fig. 2a).
The reconnection, and associated expansion of the plasma due to heating, drives out-flows from the x-point.One of these out-flows is directed towards the surface of the Sun (see the blue patch in Fig. 2b), and the other is directed upwards (the red jet in Fig. 2b).Note that there is no overhead flux rope contained in this magnetic field configuration.We replicate only the portion of the standard solar flare model below the overlaying flux rope.In these experiments the magnetic field strength is chosen to reproduce realistic values in the corona, near the reconnection site, rather than at chromospheric and photospheric heights.Indeed, our maximal field values are on the order of 2B 0 in the lower corona.The chromospheric field strength values will be made more realistic in future works.
Electrons are accelerated from the reconnection site, in the out-flow regions, and around magnetic islands/plasmoids (see the green regions in Fig. 2c and d).Such plasmoids are often caused by tearing events in thin current sheets in 2D simulations.The transport time for these electrons is considered to be shorter than the hydrodynamic time-steps of the simulation (Siversky & Zharkova 2009) and thus their energy transport is modelled as instantaneous.The energy deposition sites are shown in the lower panels using yellow colouration.This colouration is saturated at relatively low intensities to highlight the entire paths of the electrons, however, the energy deposition is actually focused in fairly concentrated kernels at the chromospheric foot-points of the flare loops.Particle trapping is possible in our beam model due to mirroring, and depends on the adopted beam pitch angle (Ruan et al. 2020).In such cases the energetic electrons remain on portions of the field-line in the next time-step of the simulation.In practice, the (yellow) beam visualizations of the energetic electrons act as a proxy to highlight the separatrices of the reconnected field-lines from those which are not currently reconnected.The inner regions of the electron energy deposition (i.e. x-locations closer to zero) are due to recently reconnected fieldlines still accelerating electrons in the X-point out-flows or in plasmoids, but also shows loops that have retained some of their energy from earlier times due to trapping of energetic electrons.Note that there are no specific mechanisms in these simulations to replicate particle acceleration in the termination shock or in turbulent flare loop-tops.However, this could be enacted in the future by judiciously generalizing the current heuristic recipes for the beams.
In these models the reconnection progresses rapidly from the start of the experiment, thus the out-flow jets impact somewhat directly on the chromosphere (see Fig. 2 c).In a pure-MHD (no beam electrons), but 3D model, Ruan et al. (2023) first initiated a gentle reconnection phase that led to the formation of a loop arcade before the impulsive phase began.The impulsive out-flow under these circumstances impacts first upon the loop-tops of this arcade before reaching down the field lines to the chromosphere.For ease of comparison to the experiments of Ruan et al. (2020) we replicate their set-up, resulting in a more direct impact of the out-flows on the lower atmosphere, and leave investigation of variations to a separate paper (Druett et al. 2023).The impact of the reconnection out-flow on the lower atmosphere, and the heating of the lower atmosphere due to other processes such as thermal conduction, causes chromospheric evaporation.This is the heating of initially cool chromospheric material up to coronal values, and its associated expansion and up-flow into the coronal flare loops.
Material ejected downwards from the reconnection outflows, as well as upflows from the chromospheric evaporation, increase the densities in the hot flare loops, and turbulence is also seen in the loop-tops below the termination shock (Fig. 2d).Although the chromosphere is only treated in single-fluid, nonradiative MHD here, we will also inspect the downward propagating shocks in these models that are the equivalents of "downflowing chromospheric compressions", and inspect the momentum they supply to the photosphere.We also calculate the SXR and HXR outputs, but present only a few relevant parameters for our analysis.The X-ray periodicity, light-curves, and other synthetic observables will be discussed in detail in a future work.

Free parameters
In the models presented here, there are several free parameters.The electron beams have energy profiles defined via a spectral index, lower cutoff energy, and initial mean pitch angle distribution.All of these are currently set to pre-determined values as in Ruan et al. (2020) (δ = 4, E c = 20 keV, and θ = 18 • respectively) and will be updated to be based on relevant atmospheric parameters in a future work.
The evolution of the model is also controlled by the description of the anomalous resistivity involving a switch between resistivity schemes described in Ruan et al. (2020Ruan et al. ( , 2023)).Manipulation of resistivity parameters are presented in Druett et al. (2023), where we demonstrate how these can lead to electron beam-driven evaporation.
The geometry of the flare is determined by the strength of the background magnetic field strength (B 0 ), the height of the resistivity patch, whether or not we insist on left-right symmetry, and thermodynamic values that initialise the atmosphere and the magnetic field structure.In this work we will focus on the variations of the background magnetic field strength, B 0 .

Results
In section 3.1 we discuss how the variations of the background magnetic field strength impacts the magnetic reconnection and out-flows.In section 3.2 we analyse the impacts on the lower solar atmosphere of the beam energetics (sec 3.2.1)and the reconnection out-flows which has, to date, largely been overlooked in multi-dimensional flare simulations.We analyse down-flows and chromospheric compressions (sec 3.2.2),up-flows and chromospheric evaporation (sec 3.2.3).In section 3.3 we discuss the formation of the flare loop arcade and turbulence in the looptops.Section 3.4 analyzes 1D cuts along field-lines to investigate the evolution of the flare simulation in selected flux tubes.Much of the physics in flares is magnetic-field aligned, and there is a long history of detailed one-dimensional (1D) flare simulations.This section establishes a basis for the comparison of results in 1D with multi-dimensional research.In multi-dimensional experiments there is a diversity of flux tube configurations.Thus, in sections 3.4.1,3.4.2,and 3.4.3we present the results for three such flux tubes with foot-points at x = −10 Mm, x = −12.5 Mm, x = −15 Mm respectively, to provide a more complete picture of the field-aligned physics.

Reconnection and out-flow
The initial atmospheres of each of the four experiments have identical thermodynamic variables.For the subsequent evolution, it is only difference in background magnetic field strengths that affects the release of magnetic energy.
The first hydrodynamic sign of the reconnection in each simulation is in the conversion of magnetic energy into internal energy along the reconnecting field-lines via Joule heating, which occurs before the electron acceleration process is switched on.By summing energy components over the entire domain in each simulation (not shown here for brevity) we see that the conversion of magnetic energy into internal energy continues relatively gently, while the conversion into kinetic energy accelerates with the development of the reconnection out-flow.Much of this kinetic energy escapes through the top boundary of the models or is re-converted into internal energy when the reconnection outflow impacts the lower atmosphere.These down-flows also transiently and locally raise the magnetic energy when the out-flow compresses the magnetic loops down onto the lower atmosphere, before they rebound.
The newly reconnected magnetic configuration, generates a Lorentz force.It is the combination of the heating and the suddenly altered pressure and Lorentz force values which drives the subsequent acceleration of the plasma away from the reconnection X-point.This out-flow (see Fig. 3) forms a "lobster claw" shape for reasons discussed in Zenitani & Miyoshi (2011), namely that in the fast-mode shock the density is highest in the central location and decreases away from the centre.This feature can be seen in the out-flow velocity plots (top row of Fig. 3) with the velocity increasing with background field strength, due to the faster rate of energy release.The heating (second row, showing temperatures) is concentrated in the tails of these out-flows (see Fig. 3 central panels, green arrows), and to a lesser extent in the outer edges of the out-flows (blue arrows).The high density regions (bottom row) are concentrated in the central locations and some distance behind their leading edge.This core of high density material is more compact for simulations with stronger background magnetic field strengths (For context, see also the Alfvén Mach numbers of different sections of the "lobster claw" out-flow jets presented in Fig. 4).
Once the energetic electrons are activated in these models, the Joule heating energy term is removed from the energy equation and instead put into the acceleration of energetic particles in regions where the drift velocities of the electrons exceed a threshold value.In these locations the out-flows are still generated by the Lorentz force and other heating that results from the magnetic realignment and energy release, for example, shock heating and adiabatic compression (equation 3: ∇ • (p gas v) and thermal conduction (equation 3: ∇ • (κ∇T )).
In the bottom panels of Fig 3 the energetic electron acceleration regions (green) and energy deposition regions (yellow) are shown using a logarithmic colour-bars, so that the areas in which they are present is well highlighted.Again, these energy quantities are higher for the experiments with greater background magnetic field strengths, as the Joule heating term increases with the liberated magnetic energy.Note that we chose to visualize the 4 experiments in Fig. 3 at different experiment times, but at similar magnetic morphological times as seen in the selected field lines shown.Thus the panel of this figure showing the B 0 = 65 G experiment does not have energetic beam electrons because they are switched on at t = 31.2s in all of the simulations.
The dense core of the lobster claw shock accelerates towards the local Alfvén speed (see Fig. 4, upper row, red arrows), with the claws travelling at significantly sub-Alfvénic speeds (see Fig. 4, upper row, magenta arrows).Out-flows from the continuing reconnection increase in velocity to become a fast-mode shock.The fast shock forms in this tail of the out-flow (see Fig. 4, green arrows).Independent of the background magnetic field strength, the out-flows reach a similar Alfvén Mach number of 9-10 in each simulation (see the solid lines in Fig. 4, lower panel).
To examine whether this maximal out-flow and Mach number is persistent, we also varied the free parameters that determine the anomalous resistivity values.Examples of the B 0 = 35 G experiment were run for 160 s of solar time with the anomalous resistivity a factor two greater and smaller values, and maximum threshold values as described in Ruan et al. (2020), equation (11) also increased or decreased by the same factor.Results of these experiments are included in Fig. 4, lower panel and confirm that the limiting Alfvén mach number of the fast shock in the out-flows of the flare obtain similar maximum values independent of the resistivity or background magnetic field strength.Spiky behavior is due to the turbulent region overlapping with the diagnosed area, which was a fixed spatial box across all experiments, based on the typical region of the reconnection outflow jet.Variations of the height of the resistivity patch, or asymmetries could also impact the maximum Alfvén mach number of the out-flow, which will be addressed in future work.
The timing of the out-flow jet reaching this Mach number does not coincide with the timing of the maximum out-flow velocities reached in each experiment (compare the peaks in the solid and dashed lines in Figure .4).In the stronger flare models the current sheet thins further and plasmoids form due to tearing instability, including a case of plasmoid coalescence in the experiment with B 0 = 65 G (see e.g., Keppens et al. 2013;Sen & Keppens 2022).
The reconnection rate in the corona can be characterised using the ideal electric field, which is given by −v × B in the reconnection region and which, for the set-up used here, has a magnitude |v x B y |.Ruan et al. (2020) found that the sweeping of the footpoints in their simulations, located using the peak value in the footpoint HXR signal, related to this coronal reconnection flow via the relationship presented in Isobe et al. (2002), 5 shows the values of |v x B y | in the corona (solid lines) and |v x(HXR) B y | shows values for the chromospheric footpoints (dashed lines).The units are converted into those of CGS ideal electric field to aid comparison with reconnection and acceleration studies which often use these units.We automated the calculations, in contrast to the previous hand-made calculations of Ruan et al. (2020) and v x(HXR) values were taken from the apparent horizontal motion.To account for the slow movement in terms of grid cell number in the footpoints, a ten-second moving average was used for the (signed) value of v x(HXR) .Once the flare loop system is stabilised a clear periodicity of the measured reconnection rates appears in the simulation with B 0 = 35 G.This occurs in both the foot-points and the x-point, with the foot-point reconnection measure (about 10 s) varying at half the period of the loop-top measure (about 20 s).
The relationship |v x B y | CORONA = |v x(HXR) B y | FOOT POINT holds much more consistently at later times in the experiment, once the flare loop system is formed, as was the case in Ruan et al. (2020).It is in reasonable agreement across the range of background magnetic field strengths, B 0 .This is consistent with Isobe et al. (2002), which presented findings only after the initial phase of the flare, although we note numerous caveats for this comparison including the very different timescales involved.
The spikes for the B 0 = 65 G experiment at around t = 80 − 90 s are produced by a passing plasmoid which increases the velocities and field strengths in the corona, as well as effecting the HXR footpoint locations in the chromosphere.The photospheric and chromospheric magnetic field strengths in our experiments are lower than reported solar flare field strengths by more than an order of magnitude.Typical solar magnetic field has a strong vertical gradient that is not present in our experiment.Thus, for more realistic flaring atmospheres one would expect much faster reconnection inflows in the corona than we find, if the footpoint sweeping speed was similar.

Electron beam energetics
1D models of flares with energetic electron heated lower atmospheres generally do not use self-consistent energy schemes, instead injecting fluxes of high-energy electrons as functions of time at the tops of the models.The energies of these fluxes can be fixed to particular values or time profiles (Allred et al. 2005, 2015; Druett & Zharkova 2018, 2019), or can be driven by observational constraints (Druett et al. 2017;Polito et al. 2023).Figure 6 shows the chromospheric electron beam heating in our models, which can be compared with values used in 1D models like RADYN (Allred et al. 2005(Allred et al. , 2015)), HYDRO2GEN (Druett & Zharkova 2018, 2019), and FLARIX (Varady et al. 2010;Heinzel et al. 2017).To compare a 1D beam model that uses a time-profile injected input heating with our multi-dimensional models, one should take a slice at a constant position (vertical slice) and read off the variations in footpoint heating flux.From the figure, one can see that our models have characteristic beam  injection duration times of around 5-20 seconds, in line with some "elementary burst" models used in 1D simulations.
It is clear that the model with B 0 = 20 G represents a very weak beam injection, with "F8" energy fluxes, i.e. an input energy flux F 0 on the order of 10 8 erg cm −2 s −1 , peaking at values greater than 10 9 erg cm −2 s −1 at only a few locations within the domain.The B 0 = 35 G experiment is a reasonable analogue of a weak "F9" elementary burst model at most locations, although there is an absolute maximum flux value throughout the domain of 1.4×10 10 erg cm −2 s −1 .The stronger B 0 = 50 G and B 0 = 65 G models represent elementary burst models with duration of 5 to 20 s with moderate electron beam fluxes on the order of "F10", i.e. with F 0 ≈ 10 10 erg cm −2 s −1 .On the basis of 1D modelling results in the literature one would expect the stronger flare models to produce some evaporation (hot up-flows) and cooler "down-flowing chromospheric compressions" signatures as a result of the beam electrons.Figure 7 shows kinetic energy maps of the flare experiments at times after the electrons are switched on and before the impacts of the reconnection jets on the lower atmosphere.The kinetic energy signatures are shown in red with the electron energy deposition locations shown in blue.Before the impact of the reconnection out-flow jets there are indeed (minor) signatures in the red kinetic energy plots of up-flows and down-flows in these experiments (with locations indicated by blue arrows in Fig. 7).However, the reconnection out-flow jets which arrive and impact a bit later completely swamp these beam-driven evaporation signatures.The weaker flares have a much longer time window between the start of the beam heating and the impact of the reconnection jet, making it appear as if they have a greater influence on the lower atmosphere.When we instead look at similar times after the switching on of the beam, the stronger flares have stronger beams that exert a greater rate of influence, in line with what would be expected from their higher beam fluxes.In a separate paper we adapt the models presented here to investigate chromospheric evaporation driven primarily by electron beams (Druett et al. 2023).

Chromospheric down-flows
In our MHD models the lower atmosphere is highly simplified.It is treated as a fully-ionised hydrogen plasma with a simple radiative loss function.The photospheric field strengths are of order 50 G, in rather stark contrast to the typical observationally derived values of a flare's lower atmospheric field strength, which are several kiloGauss.This seems to be a common situation for flare simulations derived to model coronal conditions, for example Ruan et al. (2020Ruan et al. ( , 2023)); Shen et al. (2022).Simulations developed originally from photospheric models that have been extended to accurately reproduce coronal conditions do not have this proviso, for example Cheung et al. (2019)  Chromospheric beam heating (erg cm 2 s 1 ) Fig. 6.Electron beam heating in the chromospheric foot-points.Each panel shows results for an experiment with a different background magnetic field strength.These are shown as functions of time (y-axis), and footpoint location (x-axis).The heating at each footpoint is computed by integrating the source term for the electron beam heating rate.We integrate this quantity over a vertical distance in the spatial domain that spans from the lower boundary of the experiment up to (but not including) the grid cell where the temperature first exceeds 50,000K.The figure then shows the electron beam flux density that is applied down through the top of the "chromospheric" material and deposited at each footpoint.The colourmap saturates to red at the low end.This occurs at a beam strength of F8 (F 0 = 10 × 10 8 erg cm −2 s −1 ), and thus the red colour indicates negligible or zero beam heating.(2023).We shall refer to the low temperature, high density lower atmosphere region as the chromosphere despite its simplicity in our models, and to the region at the very base of our model as the photosphere, although we do not accurately reproduce this region of the Sun.
The instant the energetic electrons are switched on, they reach the lower atmosphere, as per the modelling assumptions.In all but the weakest field strength modelled, B 0 = 20 G, there are electrons that reach our photosphere, i.e. the base of the model, see Figs. 8 and 9, top panels.The beam model used here, when implemented in 1D models, generally results in electrons being stopped at greater heights in the chromosphere (Emslie 1978;Allred et al. 2005Allred et al. , 2015)).The electrons do not impart directed momentum on the plasma in these simulations, acting only through a source term in the energy equation (3).
Before the reconnection jets impact the chromosphere, the electron beams already heat the chromospheric footpoints from initial ∼ 6000 K temperatures up to ∼ 20, 000 K over a range of heights that extends down to around 1.5Mm, in the B 0 = 20 G experiment (Fig. 8, left panels).For the B 0 = 65 G ex- periment there is heating of plasma by the beam electrons to T ∼ 50, 000 K just above 2 Mm, and to T ∼ 20, 000 K at heights as low as 1 Mm.Moreover, for the B 0 = 65 G experiment this occurs before the impact of the reconnection out-flow arrives, despite the only 2 to 3 seconds delay between the switching on of the electron beams and the arrival of the reconnection out-flow jets (see B 0 = 65 G experiment at t = 34 s in video form of the Fig. 9).However, the beam heating does not cause significant upflows in any of the experiments presented here.In the stronger flares there is not significant time for up-flows to form before the impact of the reconnection out-flow jet.In weaker flares the heating and expansion of the plasma has time to cause a gentle chromospheric up-flow (with chromospheric densities), reaching up to around 3 Mm (Fig. 8, top-left panel, indicated by a blue arrow), but this does not qualify as chromospheric evaporation as it does not rise above this height.
The different stages of the down-flows in the lower atmosphere can be seen for each model at t = 80 s in Fig. 8.In the left panels (weak B 0 ) we see a lower atmosphere after the beam electrons have started heating it, but before the reconnection out-flows impact it.When the out-flows from the reconnection do impact the lower atmosphere they transfer downward momentum and kinetic energy, as well as increasing the pressure.There is conduction of thermal energy along field-lines due to the temperature gradient.These processes heat the lower atmosphere and push it downwards.Figures 8 and 9 (top panels) show dense impact fronts at the top of the hot flare chromospheres highlighted with red arrows.Above the flare chromosphere, in simulations with stronger background magnetic field, there is a stronger conversion of chromospheric material to hot plasma that up-flows into the coronal loops.This will be discussed in the next section, and these up-flows can be seen as grey patches of increased number density in the coronas of the top panels.They are also visible as coronal kinetic energy signatures in the lower panels, with both signatures indicated with blue arrows in the figures.
Meanwhile, below, the down-flow starts to slow and cool as it travels to the photosphere (Fig. 8).Some downward travelling material is heated up to around ∼ 20, 000 K (green arrows) and below this there is very significant kinetic energy that travels down to the photosphere (magenta arrows).The left panels (a-b-c) of Fig. 10 shows the downward fluxes of kinetic energy, the maximum kinetic energy density and the maximum downward velocities at different heights through the atmosphere of the simulations with B 0 = 65 G.These heights were chosen at least five grid-points away from the experiment lower boundary, to avoid significant influence from boundary effects.The peak of the downward kinetic energy flux at a height of 300 km above the photosphere across the simulations with different B 0 values ranges from 2 × 10 25 erg s −1 (B 0 = 20 G) to 3 × 10 27 erg s −1 (B 0 = 65 G) with peak flux densities from 5 × 10 6 erg cm −2 s −1 to 4 × 10 8 erg cm −2 s −1 .These fluxes were essentially unchanged across in simulations where we kept B 0 constant and varied the initial mean pitch angles (not presented here), due to the relatively low energy fluxes achieved via the energetic particles.
The "down-flowing chromospheric compressions" travel initially as acoustic shocks (at speeds greater than the sound speed just below them, which has typical values of 8-10 km s −1 ).The "down-flowing chromospheric compression" in the B 0 = 20G flare ceases to be a shock in the mid chromosphere, when its downward velocity drops below 8 km s −1 .This transition occurs deeper into the model for increasing B 0 , but even in the strongest flare, the compression is travelling below the sound speed by the time it reaches a height of 200-300 km above the photosphere.Therefore, the "down-flowing chromospheric compressions" in these simulations would not be expected to cause a sunquake (Kosovichev & Zharkova 2001;Macrae et al. 2018;Zharkova et al. 2020) when they move below the photosphere.

Evaporation
Fig. 9 illustrates that the area of the chromosphere that gets heated, compressed, or evaporated due to the flare is greater with increasing field strength, both in depth (due to the higher energy fluxes) and in lateral area (due to the faster speed of the leading edge of the flare ribbon that results from the faster reconnection rates).The fraction of the mass at heights between 300 km and 5 Mm that is at temperatures less than 20,000K is quantified over time in Fig. 10d, alongside the chromospheric mass flux through a horizontal line at 5 Mm through the experiment (Fig. 10e).These values vary co-temporally, and don't show any in-phase behaviours with times of large photospheric mass fluxes.
Due to the mass density gradient with height that occurs between the photosphere and the chromosphere, the general plasma motions in our simulations have larger mass fluxes through the bottom boundary than through the tops of the chromosphere.These are especially larger when the front from the chromospheric compression reaches the lower boundary.These fluxes were checked and the chromospheric mass fraction in panel d of Fig. 10 showed in-phase variations with the chromospheric mass flux but no in-phase variations with these photospheric fluxes.Thus, we can be confident that the evolution of chromospheric material in panel d is due to chromospheric rather than photospheric changes.
The mass flux via evaporation (panel e), and the decrease in chromospheric material over time (panel d), both increase with background magnetic field strength between B 0 = 35 G and B 0 = 65 G.In the case of B 0 = 20 G we see that the net mass flux is at first consistently downward from the corona to the chromosphere, and the change in the fraction of material at chromospheric heights that has temperatures greater than 20,000 K varies relatively negligibly.Thus, it is shown to be possible  for high density flare loop systems to form that are fed primarily from the coronal reconnection jets and not from the chromospheric footpoints, since dense flare loops form in all the simulations including the B 0 = 20 G experiment that shows net downward mass flux through the upper layer of the chromosphere.
There is some chromospheric evaporation in each simulation, despite the net downwards mass flux through the 5 Mm height plane in the simulation with B 0 = 20 G.The maximum velocity of these evaporations at heights of 5 Mm are shown in Fig. 10f.These velocities scale from around 200 km s −1 to around 500 km s −1 over the simulations in reasonable agreement with typical values derived from ultra-violet line observations of chromospheric evaporations and 1D flare models (Kennedy et al. 2015;Polito et al. 2016;Druett et al. 2017).There are small peaks of high velocity that exceed the typically observed evaporation speeds in ultra-violet spectral lines, which reach up to nearly 800 km s −1 in the experiment with B 0 = 65 G.

Flare loop-tops
Figure 11 presents kinetic energy density maps which highlight the current sheets and flare loops.We use them to illustrate the evolution of the flare loop-tops in the simulations with different magnetic field values.The top row shows the impact of the reconnection out-flow jets onto the chromosphere.Strong lateral flows away from the center of the impact move along the top of the chromosphere, and a shock-wave expands as a dome over the whole coronal domain outside of the flare loops, centered on those loops.This phenomenon is also present in the experiment with B 0 = 20 G but with kinetic energy densities values close to the lower saturation limit that make it hard to see in the figure .The impact and rebound of the reconnection out-flow on the lower atmosphere, as well as heat conduction, create strong evaporation flows up the flare loops from their footpoints.The densities and velocities of these evaporation flows scale with the density and velocity of the impacting reconnection out-flow jet, and thus with the background field strength as described in section 3.2.3.We note the absence of strong rebound up-flows in the B 0 = 20 G model, which is consistent with the lack of variation in chromospheric mass fraction and negative chromospheric mass flux for this model seen in Fig. 10(d-e).Evaporation upflows do also begin for the experiment with B 0 = 20 G at some time after the time shown in Fig. 11a, as can be seen in panel Fig. 11e.
The down-flows from the loop-tops that are concurrent with the onset of chromospheric evaporation are slow shocks caused by the compression of the loop-top region during the impact, and the resulting negative pressure gradient outwards from the central position of the loop-top in the directions along the fieldlines.
The second row of diagrams shows the flare loop arcade after the rebound of the reconnection out-flow on the lower atmosphere.In each experiment the loop-tops display oblique and horizontal fast shocks, and potentially multiple termination shocks, as described in Takasao et al. (2015); Takasao & Shibata (2016).This pattern of shocks is a consistent feature across all the experiments.However, these typical behaviours as described in Takasao et al. (2015); Takasao & Shibata (2016) can be disrupted by plasmoids which are ejected downwards from the current sheet in the experiments with stronger background field strength (see e.g.Fig. 11g, with B 0 = 50 G).
The magnetic tuning fork process (Takasao & Shibata 2016) is a flare loop-top oscillation that is controlled by the back-flow of the reconnection out-flow.It is evident across all the experiments.This is shown in the lower two rows of panels, in which turbulent eddies form on each side of the termination shock, with a dominant extent that alternates between the left and right sides.These eddies are also associated with pulses of shocks that propagate out into the surrounding plasma and can be identified as sets of fringes in the kinetic energy density that move away from the loop-top locations in the lowest two rows of panels in Fig. 11.
The magnetic tuning fork process and the plasmoids are both capable of sending high density flows outward around the turbulent loop-tops.Sometimes these flows are dense enough to inter-cept a significant fraction of the energy in the energetic beams of electrons before they can reach the chromosphere.This phenomenon can be seen in the video version of Fig. 11 in which the beam energy deposition (blue colour) increases in the coronal region of the experiment (see B 0 = 65 G before and after t = 100 s, B 0 = 50 G at around t = 123 to t = 129 s, and B 0 = 35 G after t = 200 s).This interception of beam electrons is also visible in the chromospheric beam heating rates at the footpoints of the models in Fig. 6.It cna be seen through the decreases in the electron beam energy reaching the chromosphere, due to the deposition of a significant portion of the beam fluxes in the corona.
In the B 0 = 35 G experiment we see a periodic brightening of the beam energy deposition rate in the coronal loops after t = 220 s.This appears to be associated with the magnetic tuning fork process which periodically (around every 16 s) emits shocks into the surrounding loops.These shocks periodically increase the densities in field-lines that have recently reconnected, reaching maximum around the time that the tuning fork pulse passes from one side of the loop-top to the other, and occurring on both sides of the experiment.This increased density, in turn, removes energy from the beams of electrons before they reach the footpoints, with a periodicity of around 16 s.Faint traces of this process are visible in Fig. 11n and the video versions of this figure.Thus the magnetic tuning fork can directly affect the loop-top X-ray emission as described in other papers (Takasao & Shibata 2016;McLaughlin et al. 2018;Zimovets et al. 2021).For example, the waves emitted from the loop-top align with the concept of a periodic fast-mode magnetoacoustic wave as per the analysis of Takasao & Shibata (2016), which used similar models to those presented here.It has been argued that such waves propagating towards reconnection x-points may also generate quasi-periodic-pulsations (QPPs, McLaughlin et al. 2009).Here we describe a multidimensional secondary process which can also explain simultaneous QPPs in footpoint and coronal loop HXR sources.Pure MHD models without beam electrons cannot self-consistently quantify these HXR or beam-related effects and their related to QPP variations.However, the lack of particle acceleration modelling for non-thermal electrons in the termination shock and turbulent reconnection loop-top regions means that the origins of QPPs in HXR sources cannot be definitively discerned from our models.In follow-up work we will produce a more rigorous analysis with direct synthetic images in wavebands relevant to these processes.
The experiment with B 0 = 50 G also exhibits periodic pulsing of beam energy deposition in the coronal loops, but in this and the B 0 = 65 G case the main factor in the disruption of the flow of electron beams from the X-point to the chromosphere is plasmoids.The reduction in beam electron footpoint heating in the case of B 0 = 50 G is evident at 120 s, just after a plasmoid strikes the loop-tops (Fig. 6(c)), with the energy deposition around this time shown in Fig. 11(o).For the simulations with B 0 = 65 G plasmoids can be seen approaching the loop-tops at around 95-100 s in Fig. 11(p), and the subsequent reduction in beam heating of the chromospheric footpoints is evident for the rest of the experiment in Fig. 6(d).
Within a post-processing full kinetic model that was evolved over the backdrop of a resistive MHD-simulated atmosphere, Kong et al. (2020) found that the interaction of plasmoids with the termination shock in the loop-tops significantly modulated and softened (increased the negative exponent of the power law energy distribution) the electron beams accelerated.This attenuation and shift towards lower energy electrons was related to the collisions of the plasmoids with the looptops, because these collisions increased the number of grid cells showing compressions.Such modelling is highly valuable and may contribute to future works including recipes for the accelerated electron spectrum in longer, larger scale simulations at affordable computing costs.Our work does not include a detailed kinetic model for the accelerated spectrum of electrons, but the effect of the plasmoids is, rather coincidentally, similar for the thresholding of electrons reaching the chromosphere.In our experiment the collision of the plasmoid with the looptops causes density and magnetic field variations in the corona that result in both increased beam en-ergy deposition near the termination shock (with reductions at the footpoints) and of greater particle trapping.
We generated SXR curves for each flare model (Not shown in figures) in Watts per metre squared, in order to produce GOES classifications.For this we use the method outlined in Ruan et al. (2018) based on the work of Del Zanna et al. (2015);Pinto et al. (2015); Fang et al. (2016b).Their formula expresses fluxes in photon cm −2 s −1 .We adapt this by multiplying the integrand by the photon energy to produce a result in erg cm −2 s −1 .The integral is taken between limits with energies corresponding to those of the GOES 1-8 Å band and then converted to W m −2 .It is the peak of the flux in this GOES channel that defines the standard X-ray classification of a solar flare (Baker 1970).These values need to be multiplied by a representative depth in 2.5D models.For consistency with the work of Ruan et al. (2023) we choose this depth to be 100 Mm.Assuming an arcade of this depth we obtain the data seen in table 1.The flare classifications are spread across a reasonable span of the observed range on the Sun, but do not reach the X-class flare classification.

Flows along a field-line
In section 3.2.1,statistics for the electron beam deposition were presented at each base point of the models.Now that we have presented the multi-dimensional evolution of the lower and upper atmosphere we inspect the variations in 1D of parameters along individual field-lines.Many of the physical processes in flare loop systems are field-aligned, and so there is significant value to inspecting the dynamics along such cuts.Moreover, this analysis provides a much greater basis than currently exists in the literature to enable the comparison of results of flare simulations in multiple dimensions with decades of research results derived from detailed 1D radiation hydrodynamic modelling of flare loops.
For this investigation, we inspect the strongest (M2 class) flare with B 0 = 65 G and a maximum beam electron energy flux over all space and time of 4.7 × 10 10 erg cm −2 s −1 .Field-lines with footpoints at x = −10, −12.5 and −15 Mm are selected as representatives of the variety of locations available within this multi-dimensional morphology.
Plasma number densities, vertical velocities, temperatures, and kinetic energy densities are shown as functions of distance from the loop apex (y-axis) and time (x-axis) in Figs. 12, 13, and 14.The electron fluxes along these field-lines are shown as functions of time, over-plotted as red lines on these images.These views form the direct counterpart of restricted 1D hydrodynamic models (e.g.Fig. 6 of Unverferth & Reep 2023), and can be compared readily.
Before magnetic reconnection the tracked field-lines exit the experiment at the top boundary, and after reconnection they reach to a footpoint on the opposite side of the flare loop system.Thus, there is a sharp disconnect between the top halves of these panels at times before and after the reconnection, as this portion of the field-line tracks completely different plasma.After the field-lines reconnect, they rapidly retract and shorten.This can be seen by the rapidly decreasing total length of the fieldlines between the photospheric footpoints of the field-lines.The footpoints are plotted as green lines at the tops and bottoms of the panels.From figure 12 one can see that this field-line reconnects at t = 67 s.This is 4 s after t = 63 s, when the beam electrons switch on for this field line, accelerated in regions away from the x-point.The onset of the evaporation can be seen from the lower footpoint at t = 65 s.This requires some explanation with respect to the standard flare model.Before reconnection there is significant heating (to around T = 10 MK) on this field-line from energy supplied by neighbouring plasma via methods including heat diffusion and the expansion of the neighbouring flux tubes which generates compressive heating of the flux tube we are inspecting.The heating intrudes at lengths/heights around s = −25 Mm at t = 53 s.This causes upward and downward flows to expand outward from this point.More dramatic heating occurs at t = 60 s just before the reconnection time, both at the similar heights the previous heating and near the reconnection region higher up the open field-line.This heating results in heat conduction and velocity flows.It is this conduction front approaching the chromosphere which drives evaporation up from the footpoints.The evaporated material continues to be heated and expand, driving further acceleration up to around 300 km s −1 by the time the evaporation reaches y = 20 or s = −10 Mm, at temperatures of around 2 MK.
Meanwhile, the top of the reconnected loop collapses downwards at high velocity, and shortens in total length.This process is visible as a dark blue horizontal stripe in vertical velocity (down-flow) immediately after the reconnection event, seen in panel (b).The apex of the line shortens until it is around s = 20 to s = 30 Mm away from the footpoints.The compression of the loop-tops drives a series of hot (50-100 MK) outward (downward) flows from the apex, starting at around t = 80 s.The velocity plot, kinetic energy plot, and temperature structure of the loop-tops show that the region is undergoing turbulence as well as heating events.
The downward flows driven from the loop-tops meet the rising chromospheric evaporation.The evaporation and downward travelling loop-top flows shock when they meet, reducing the velocities of both streams and heating the upward moving material significantly.This can be seen in the changes of gradients of the rising and falling density fronts in figure 12a at around s = ±8 Mm, at t = 80 s.In the beam driven evaporation model of Druett et al. (2023) there was no strong reconnection out-flow to compress the loop tops and drive downward flow that meet the rising evaporation, in that study a significant portion of the evaporated plasma passed over the loop-tops and down to the other side of the arcade.In contrast, in this study the upward evaporations do not directly pass over the loop-top region, but get caught in the turbulence until gentle down-flows form at around t = 140 s.The direct passing of up-flows over the loop apex and down towards the opposite footpoint can also be seen in 1D loop models such as in the central and right panels of Fig 6 in Unverferth & Reep (2023) which also shows no loop-top turbulence, due to the 1D nature of the model.
Along this field-line, pulses of extra density and kinetic energy rise upwards from both footpoints after around t = 85 s which become broader, slower and weaker over time.This occurs after the beam electron processes have ceased along the field-line, possibly indicating some wave-like behaviour.

B
Figure 13 shows that, predictably, the field line further out reconnects later (t = 85 s).Both the beam acceleration and the  1.F 0 and GOES classifications of the simulated flares.The left column shows the background field strength of each simulation.The second gives the F 0 classification, and the third the flux of the electron energy deposition at its maximum value.The time and location of this maximum is shown in the 4th and 5th columns.The sixth column gives the GOES SXR classification as per the scheme of Baker (1970), as described in Pietrow (2022), and the seventh gives the peak of the flux in the GOES 1-8Å channel, assuming a third dimension depth of 100 Mm, which scales the flux linearly.The final column shows the time at which this maximum GOES flux occurs.
evaporation processes from the footpoints start at 80 s.This is also co-temporal with the arrival of the fast downward propagating hot jet due to heating and expansion of neighbouring material that heats this flux tube around s = −25 Mm beginning at around t = 75 s.The driving of the evaporation is broadly similar to that described for the footpoint at x = −10 Mm.After reconnection the loop-top collapses downwards at greater velocities, and halts with a wider span of s values.The simulation ends before the gentle down-flows from the loop-tops can reach the footpoints.
Figure 14 shows that this field-line experiences chromospheric evaporation (at t = 94 s) well before reconnecting (at t = 102 s), and before the beam electrons reach the chromosphere (at t = 104 s).The heat driven expansion from the nearby loops begins at positions near s = −25 Mm.Again, the hot plasma expands in upwards and downwards directions from s = −25 Mm.This flow reaches the chromosphere at t = 93 s, and immediately drives chromospheric evaporation which achieves similar speeds of around 300 km s −1 .The evaporation collides with down-flows at positions of around s = 20 Mm.For this field-line which is situated in the outer region of the flare, by the time it reconnects, significant loop-top turbulence is already present.This can be seen through the high-speed, direction-varying velocities in the loop-tops (panel b).Also, before the time of the reconnection of the field-line there are already some slightly higher-density features in the loop tops, as well as the rising chromospheric evaporation fronts.The beam electrons deposit a significant fraction of their energy in these evaporation and loop-top features.Thus the beam flux reaching the chromosphere is significantly reduced.A dedicated future investigation will be necessary to understand these effect in detail, including in experiments with evaporation driven by beam electrons (Druett et al. 2023).

Summary and discussion
In this paper we presented a study of the 2.5 D MPI-AMRVAC flares including beam electrons, first reported in Ruan et al. (2020).By varying the background magnetic field strength by a factor of 3.25 in these simulations between B 0 = 20 G and 65 G, the GOES classification of the simulation changes by over 2 orders of magnitude, between B1.5 and M2.3 (assuming a 100 Mm arcade length in the third dimension, see table 1).
The flux of energy supplied by energetic electrons at any given chomospheric footpoint has a characteristic duration of between 5 and 20 seconds, usually with a relatively triangular profile of flux against time, peaking earlier in the profile (see Figs. 6,12,13,and 14).The peak beam heating flux at a footpoint in each experiment varies between 1F9 for the case B 0 = 20 G (F 0 = 1 × 10 9 erg cm −2 s −1 ) and 5F10 for the case B 0 = 65 G (F 0 = 5 × 10 10 erg cm −2 s −1 ), over and order of magnitude difference.This is the first paper reporting the details of chromospheric beam fluxes and their evolution in multi-dimensional simulations.
In all simulations, bi-directional reconnection out-flow jets are formed in the corona at heights of 50 Mm where the initial reconnection x-point forms.The out-flows have a classic "lobster claw" shape (Zenitani & Miyoshi 2011).A fast shock exists in the tail of this feature and stabilises some time after the out-flow impacts the lower atmosphere.The maximum speed achieved in these flows scales by a similar amount to the background field, from around 1000 km s −1 to 3200 km s −1 across the experiments with B 0 = 20 G to B 0 = 65 G respectively.As a result, after the loop system settles, the maximum out-flow speed is approximately constant across the simulations when expressed as an Alfvén mach number (see Fig. 4).It is possible that this O(10) maximum value would alter based on the variation of other simulation parameters, such as the vertical position of the resistivity patch, which would provide a longer or shorter reconnection outflow jet if placed higher or lower in the atmosphere respectively.The maximum out-flow Alfvén mach number was insensitive to changes of the maximum anomalous resistivity.
We perform the first detailed investigation of chromospheric response to the impacts of reconnection out-flow jets in multidimensional models, including the changes in these responses across a variety of flare strengths.The impact of the reconnection out-flow jets and the heat conduction front on the lower atmosphere heats the chromospheric material.This generates hot up-flows (T ∼ 2 MK, chromospheric evaporation, see Figs.12 to 14, panel c) and cooler down-flows (T ∼ 20, 000 K, "downflowing chromospheric compressions", see Figs. 8 and 9).Chromospheric material is also heated from around 6000 K to temperatures around 20,000 K within a second of the beam electrons being switched on.Heating of the plasma to temperatures around 20,000 K extends downwards from the tops of the chromosphere to depths of 1.5 Mm for the B 0 = 20 G simulation and to 1.0 Mm for the B 0 = 65 G case.There are noticeable kinetic energy imprints of the beam electrons at the chromospheric footpoints of the flares after the acceleration is switched on, but these are swamped by the reconnection out-flow jet in the present simulation suite (see Fig. 7).This is investigated in more detail in a companion paper (Druett et al. 2023).
At heights of 300 km above the photosphere the downward kinetic energy flux densities reach 5 × 10 6 erg cm −2 s −1 for the experiment with B 0 = 20 G, and up to 4 × 10 8 erg cm −2 s −1 for B 0 = 65 G.This demonstrates a significant transfer of energy and momentum to the photosphere, however even the "downflowing chromospheric compression" for the strongest flare presented drops below the local sound speed at heights between 200 and 300 km above the photosphere.This implies that we would fer whether this process would also produce pulsations in visible and near-UV sources.
Our investigation of the flows along 1D field-lines reveal several multi-dimensional effects that are not accounted for in 1D studies of solar flare loops, even those attempting to recreate multi-dimensional effects such as Kerr et al. (2020); Unverferth & Reep (2023).Firstly, the reconnection and loop top turbulence sources are intimately linked to the multi-dimensional nature of the simulation.Secondly, a shock occurs when chromospheric evaporation meets downward loop-top sources that have been forced along the field-lines due to the high pressure in the loop tops.Thirdly, there are sources of chromospheric evaporation caused by heat conduction via field-aligned transport from the loop-top sources, but also from neighbouring flux tubes via processes including heat diffusion and compression.In strong flares these processes cause the leading edges of the flare ribbons to begin evaporating material by thermal conduction before the associated field lines have reconnected.These regions completely lack beam electrons and HXR sources at these times (Fig. 14).This matches the pattern of up-flows from ultraviolet satellite observations by the Interface Region Imaging Spectrometer (IRIS, De Pontieu et al. 2014), as reported in Polito et al. (2023), who note that the leading edges of the flare ribbons show significantly lower evaporation of chromospheric material than the main bodies of the flare ribbons.Polito et al. (2023) interpret that the beams of electrons in the leading edges of the flare ribbon are significantly weaker (presumably associated with acceleration processes near the X-point, in the current sheet above the flare loops) than those beams inside the body of the flare ribbon (which may be associated with acceleration sources nearer the termination shock, the loop-top turbulence, or the flare looptops themselves).Our models demonstrate that additional and alternative interpretations should be considered to those that can be provided by 1D analysis.However, modelling advances are required in our simulations to include particle acceleration in the termination shock and turbulent loop-tops before we can provide a comprehensive answer.
Finally, we note a number of improvements that can be made to these models and the benefits these would bring.This lineage of simulations has focused on accurately reproducing coronal conditions, and future simulations should improve the lower atmosphere by increasing the magnetic field strength (eventually by several orders of magnitude at the base), and making the density profile accurately represent chromospheric and photospheric values.This would improve the accuracy and credibility of interpretations derived from spectral line synthesis regarding the motions of the "down-flowing chromospheric compressions", ribbon formation, and evaporation processes, thereby better constraining the energetics and fundamental flare processes responsible for visible, UV and SXR emissions.
The beam model should be improved so that it can selfconsistently be the principle agent driving evaporation.This has recently been achieved for the first time in multiple dimensions in a companion paper, Druett et al. (2023).The energy spectrum and mean pitch angles of these beams can be parameterised based on atmospheric quantities and acceleration statistics from detailed studies, including approaches by Bakke et al. (2018);Frogner et al. (2020);Frogner & Gudiksen (2022).Effects such as self-induced electric field and return currents could also be included in the transport model (Zharkova et al. 1995).Energetic protons could also be considered in the particle transport modelling (Zharkova & Zharkov 2015).
Detailed non-Local Thermodynamic Equilibrium radiative transfer (Allred et al. 2005;Carlsson & Leenaarts 2012;Allred et al. 2015;Hong et al. 2022) and non-equilibrium ionisation (Leenaarts et al. 2007) can be included in the lower atmosphere.This could be combined with optically thick spectral line synthesis (Osborne & Milić 2021;Jenkins et al. 2023).The combination of these advances applied to our models would provide multi-dimensional context to our understanding of inhomogeneities, flows, energy delivery, and ribbon substructures observed in a large number of currently debated flare phenomena based on observations in the visible and near ultraviolet emissions of flares (Ichimoto & Kurokawa 1984;Osborne & Fletcher 2022;Pietrow 2022;Polito et al. 2023), including investigating highly broadened flare ribbon spectral line profiles (Zharkov et al. 2020;Druett et al. 2021;Kowalski et al. 2022;Kerr et al. 2024), and the mechanisms responsible for sunquakes (Kosovichev & Zharkova 2001;Quinn et al. 2019).
Further studies using this model would be instructive.An investigation of the effects of varying the height of the resistivity patch and introducing asymmetries into the simulation would examine the robustness of relationships derived using our standard set-up.Data driven simulations and a systematic comparison of our multi-dimensional simulations with 1D simulations that include detailed physics would allow us to determine the most potent admixture of these approaches to use in future studies.Furthermore, an equivalent study for simulations with evaporation driven by beam electrons, and a 3D version of the simulation would allow us to interpret observational signatures of the different evaporation mechanisms, thereby determining which processes principally drive evaporation and other fundamental flare phenomena on the Sun.

Fig. 1 .
Fig. 1.The standard solar flare model.Left: The standard solar flare model in 2D, with features labelled as per the numbered points in section 1, namely (1) a flux rope running normal to the plane of the cross section, (2) reconnecting field-lines at an X-point, (3) reconnection out-flow jets, (4) reconnected hot loops that gather below the X-point, (5) a termination shock, (6) field-aligned transport of energy down to the lower atmosphere (7) flare ribbons, and (8) chromospheric evaporation and "down-flowing chromospheric compressions".Right: Images of a solar flare over the solar limb taken with AIA aboard SDO (Lemen et al. 2012) on Sept 10th 2017, showing the form of the standard solar flare model early in the flare (at 15:51:43 UT) in the 193 Å channel with the loop arches at the base (green arrows) and a flux rope above it (blue arrows).Slightly later (at 15:53:20 UT) we see the flux rope erupting (moving upward).Later (at 16:18:07 UT) the flux rope has left the field of view and we see emission that is interpreted to originate from the long bright current sheet (red arrows), and dense coronal flare loops below (green arrows).These images use the short exposure AIA filters, but still also display significant fringe artifacts emanating outward from the bright loop arches.

Fig. 2 .
Fig. 2. The simulated analogue of the standard solar flare model.Panel (a) shows B y /B 0 .All the simulations presented here start from a bipolar field region with an anomalous resistivity patch at a height of y = 50 Mm, over the polarity inversion line, which causes the field to reconnect and release stored magnetic energy.Panel (b) shows the typical vertical velocity (v y ) structure at a time before the impact of the reconnection out-flow on the lower atmosphere.Positive velocity (red) represents upward motion and negative represents downwards motions (blue).The reconnection out-flow jets form from the X-point where reconnection occurs, one flowing downward and the other upward jet leaves to domain via the top boundary.Panel (c) shows the vertical velocity at the time when the reconnection jet impacts the lower atmosphere.Panel (d)  shows the typical flare loop system that is formed via this process, through a plot of absolute magnetic field strength.Energetic electron transport is switched on after 31.2 s of the simulation.Thus, in all cases presented in this manuscript, the electrons are switched on before the impact of the reconnection jet on the lower atmosphere.On the lower panels, electron acceleration sites are shown in green, with energy deposition locations displayed in yellow.The energy deposition locations are saturated at very low values to help indicate their paths through the experiment.In fact, these beams deposit the majority of their energy in the lower atmosphere at the footpoints of the flare loops.
seen in their Fig.4.The reconnection inflow |v x B y | CORONA was calcu-

Fig. 3 .
Fig. 3.The reconnection out-flow jets of the experiments with different background magnetic field strengths.These are shown at similar morphological stages of the experiment evolution.The columns show results for different background magnetic field strengths from B 0 = 20 G on the left to B 0 = 65 G on the right.The top panels show the vertical velocities of the models, and the central row shows temperatures.In the temperature panels, green arrows indicate the concentration of high temperature in the tails of the reconnection out-flows, the blue arrows indicate the hotter areas at the rear of the lobster claw forms that lead the reconnection out-flows.The bottom row shows plasma number density.Each panel also shows magnetic field lines in black.These are traced from footpoints at x = −2.5, −5.0, and −7.5 Mm in instances where these lines are being processed by the 1D field-line routines.In the number density panels (bottom row), beam electron acceleration sites are shown in green, and the locations where energetic electrons deposit their energy are shown in yellow.

Fig. 4 .
Fig. 4. The Alfvén mach numbers of the out-flow reach similar ranges and values at similar evolution epochs, independent of field strength.These are shown just prior to the impact of the out-flow on the chromosphere (top row), during the compression of the flare loops by the generation of the termination shock (second row), and after the rebound of the impact once the flare loops have settled (third row).In each row a green arrow highlights the fast mode shock in the reconnection out-flow.In the top row a red arrow highlights the high-density core of the lobster claw out-flow formation, and a magenta arrow points to the sub-alfvénic "claws" of this structure.In the lower panel the logarithm of the maximum downward out-flow speeds are shown with dashed lines (right axis), and the maximum of the Alfvén Mach numbers in these out-flows is shown with solid lines.

Fig. 5 .
Fig. 5.The reconnection inflow rate (solid lines) and footpoint sweeping (dashed lines) expressed as |v x B y |, in electric field units E = −v × B. The different colours show the variations in these quantities as functions of time for the experiments with different strengths of background magnetic field B 0 .

Fig. 7 .
Fig. 7. Signatures of the electron beam effects on the lower atmosphere.Kinetic energy maps of the flares are shown in red, with the lobster claw reconnection out-flows approaching the lower atmosphere for the experiments with B 0 values of (a) 20 G at t = 73 s, (b) 35 G at t = 48 s, (c) 50 G at t = 50 s, and (d) 65 G at t = 32 s.Overlays show the electron acceleration densities in green and the energy deposition regions in blue.Note that the lower panels have energy deposition rates masked and scaled to values 100 times greater than the upper panels, in order not to completely cover the footpoint kinetic energy signatures, which are seen as small red blobs (Kinetic energy) next to the blue footpoints blobs (electron energy deposition) and highlighted with blue arrows.

Fig. 8 .
Fig. 8.The heating, down-flows, and evaporation of chromospheric material.These are shown for simulations with increasing background magnetic field strengths in the columns from left to right, each at the same time during the simulation (t = 80 s).The top row shows the plasma density in grey-scale, with electron energy deposition sites overlaid in yellow.Blue arrows indicate the locations of up-flows from the chromosphere, evaporation in the case of all but the B 0 = 20 G experiment.Red arrows indicate the high density impact fronts of down-flowing material at the top of the chromosphere.The central panels show the temperature of the atmosphere, saturating to black at 6000K and to white at 30,000K.In this row green arrows indicate hot chromospheric material.The bottom rows show the kinetic energy densities with chromospheric evaporation signatures highlighted using blue arrows and significant energy transfer to photospheric levels indicated by magenta arrows.

Fig. 9 .
Fig. 9.The heating, down-flows, and evaporation of chromospheric material in the models.The formatting is similar to that of Fig. 8, including those features that are highlighted by arrows.This figure shows the simulations with different background magnetic field strengths, each at similar stages in the evolution of the flare, after the impact and rebound of the reconnection out-flow jet.

Fig. 10 .
Fig.10.Impact of the flare on the dynamics of the lower atmosphere.The left column shows results for the B 0 = 65 G simulation, while the right panels compare all four experiments.Shown are: (a) the total kinetic energy flux (assuming a 3rd dimension depth of 100 Mm) (b) the maximum kinetic energy density directed downward, and (c) the maximum downward velocity, all shown at various heights near the photosphere.In the right column of panels we show: (d) the fraction of the lower atmospheric material that is "chromospheric" as functions of time for the different experiments, (e) the mass fluxes (assuming a 3rd dimension depth of 100 Mm) and (f) the maximum upward velocities of material, each taken at 5 Mm height.

Fig. 11 .
Fig. 11.Kinetic energy density color maps.Kinetic energy density is shown in red and highlights the flare loops, the current sheets, and shocks travelling through the surrounding atmosphere.Overlays show the electron acceleration energy densities in green and the energy deposition density in blue.Energy densities for all panels of this figure have a lower saturation limit of 10 −1 erg cm −2 s −1 , and an upper limit of 10 3 erg cm −2 s −1 .The columns from left to right show results for the experiments with background magnetic field strengths of B 0 = 20 G, 35 G, 50 G, and 65 G respectively.The top row shows the simulations at the time after the impact of the leading edge of the reconnection out-flow jet.This impact and its reflection causes hot up-flows from the chromosphere.The second row shows the time at which the flare loops have rebounded after their compression during the impact.The third and fourth rows show later times when the loop-tops settle, and exhibit turbulent eddies on alternating sides of the central line, i.e. a magnetic tuning fork process.An animated version of the figure is provided in the online materials.

Fig. B. 1 .
Fig. B.1.Kinetic energy fluxes through a height of 5 Mm in the simulations, as functions of time.These are shown for the B 0 = 35 G and B 0 = 50 G experiments in the top and bottom panels respectively assuming an arcade depth (3rd dimension depth) of 100 Mm.Results are shown for the main simulation presented in this paper alongside results for a simulation with identical parameters, except for setting the beam heating term, H e = 0 at all times.