Time-dependent Stellar Flare Models of Deep Atmospheric Heating

Optical flares have been observed from magnetically active stars for many decades; unsurprisingly, the spectra and temporal evolution are complicated. For example, the shortcomings of optically thin, static slab models have long been recognized when confronted with the observations. A less incorrect—but equally simple—phenomenological T ≈ 9000 K blackbody model has instead been widely adopted in the absence of realistic (i.e., observationally tested) time-dependent, atmospheric models that are readily available. We use the RADYN code to calculate a grid of 1D radiative-hydrodynamic stellar flare models that are driven by short pulses of electron-beam heating. The flare heating rates in the low atmosphere vary over many orders of magnitude in the grid, and we show that the models with high-energy electron beams compare well to the global trends in flux ratios from impulsive-phase stellar flare, optical spectra. The models also match detailed spectral line-shape properties. We find that the pressure broadening and optical depths account for the broad components of the hydrogen Balmer γ lines in a powerful flare with echelle spectra. The self-consistent formation of the wings and nearby continuum level provides insight into how high-energy electron-beam heating evolves from the impulsive to the gradual decay phase in white-light stellar flares. The grid is publicly available, and we discuss possible applications.


INTRODUCTION
Stellar flares are powerful bursts of radiation that are observed from radio frequencies through the X-rays.They are thought to be caused by similar processes as in solar flares: a catastrophic release of magnetic energy into the corona accelerates particles, which either escape to space or return the energy to the atmosphere as they thermalize.The response of the stellar atmosphere initially occurs at lower altitudes within the chromosphere, which cools by radiating away this energy.Early in a flare, radiation can predominantly occur through optical and near-ultraviolet continua (Hawley & Pettersen 1991;Woods et al. 2004;Hawley et al. 1995;Osten & Wolk 2015;Schmitt et al. 2019), thus giving rise to a broadband visible-wavelength response that is sometimes referred to as the white light.
White-light flares dominate the short-duration variability in light curves of magnetically active stars in Kepler, K2, and TESS data.For example, the M dwarf GJ 1243 flares at a detectable level in Kepler for ≈ 30% of the time (Hawley et al. 2014).The short-wavelength extensions of the optical and near-infrared continuum in these flares are important for assessments of exoplanet atmospheric photochemistry (Loyd et al. 2018b;Tilley et al. 2019;Segura 2018) and surface biology (Howard et al. 2020;Abrevaya et al. 2020).However, there is still not a self-consistent physical theory, from magnetic energy release to atmospheric heating, that explains the spectral properties of stellar flare white-light radiation.
A phenomenological T ≈ 8500 − 14, 000 K blackbody model is consistent with some stellar flare spectra in the impulsive phase (Hawley & Fisher 1992;Fuhrmeister et al. 2008;Kowalski et al. 2013;Lalitha et al. 2013;Kowalski et al. 2016), but other flare spectra show clear signatures of Balmer jumps and cooler blackbody temperatures in the conservation are solved on an adaptive grid (Dorfi & Drury 1987) with the equations of plane-parallel radiative transfer and level populations for a six-level hydrogen atom, a nine-level He I/II/III atom, and a six-level Ca II/III ion (see Allred et al. 2015;Carlsson et al. 2023, for details).The numerical method of solution consists of a semiimplicit (multi-dimensional Newton-Raphson) linearization and iteration of the time-centered spatial derivatives with an adjustable time step (Dorfi & Feuchtinger 1995;Carlsson 1998;Dorfi 1998;Abbett 1998).An advantage of this method is that the time step is not limited to the CFL condition, which is very short in the flare transition region.Time steps in the flare simulations in our grid often exceed the Courant step by factors of O(10 3 ), as expected (Dorfi 1998), but they can also be very small and unwieldy.
The model loop geometry is a quarter circle with an arc length of l = 10 Mm and a uniform cross-sectional area.The surface gravity is log 10 g ⋆ /[cm s −2 ] = 4.75, and the calculated effective temperature is T eff ≈ 3600 K.The loop apex has an electron density of n e ≈ 3 × 10 10 cm −3 and plasma temperature of T ≈ 5 MK at t = 0 s, which are in line with those determined from quiescent X-ray spectra of dMe stars (Osten et al. 2006;Raassen et al. 2007;Liefke et al. 2008).Non-radiative heating is applied to the corona1 and to the photosphere to emulate convective energy flux (e.g., Gustafsson et al. 2008), which prevents catastrophic cooling before flare heating begins.The atmosphere begins in a self-consistent state of hydrostatic equilibrium, which is achieved through the steps that are described in the documentation associated with Carlsson et al. (2023).The "height" variable, z, is the arc distance along the loop, and z = 0 km corresponds to τ (λ = 500 nm) = τ 500 = 1 in the preflare state, where τ is the optical depth.The one-dimensional formulation of the equations of hydrodynamics is justified if β ≪ 1, where β is the ratio of gas to magnetic pressure.Our models are thus meant to represent the centers of very strong flux tubes in plage-like environments, where flare ribbons are often observed on the Sun.The model loop is relatively small and represents a compact, low-lying flare loop that appears early within a growing arcade of bright loops in EUV filtergrams (e.g., Veronig et al. 2006;Liu et al. 2013;Chen et al. 2020).A reflecting upper boundary emulates waves originating from an identical conjugate footpoint of the loop.Thus, the flare loops are symmetric with respect to the apex.We also impose a v z = 0 lower boundary condition (as in the F-CHROMA solar flare models; Carlsson et al. 2023), which does not affect the chromospheric dynamics over the short timescales (∆t ≤ 10 s) covered by these simulations.An interesting effect due to the reflected waves at the upper boundary is discussed in Section 3.1.
An electron beam is injected at the loop apex to simulate flare heating.The electron beam pitch-angles are forward-Gaussian distributed with respect to the loop axis.The volumetric flare heating rate is determined by the steady-state solution of the distribution function of energy, pitch angle, and position along the loop.The diffusion and drag of beam particles due to Coulomb collisions with the ionized and neutral components of the background (thermal) plasma are included through the Fokker-Planck formulation, which is solved in a module that was further developed into the standalone FP code that solves the distribution function with the return current force for any beam particle mass and charge (Allred et al. 2020).The heating rate is the gradient of the beam energy flux (see Allred et al. 2015).The heating rate is calculated at every time step in the radiative-hydrodynamic simulation because gradients of the atmospheric density evolve rapidly after t = 0 s in response to the flare heating.Collisional excitation and ionization rates of hydrogen (Fang et al. 1993) and helium (Arnaud & Rothenflug 1985) with the electron beam are included in the detailed rate equation following Allred et al. (2015), but these are rapidly overwhelmed by the thermal collision rates in our simulations.In this first generation of models, the collective force term (i.e., the third term in Eq. 1 of Allred et al. 2020) that accounts for return current electric fields and magnetic mirroring has not yet been included.These physical processes warrant detailed study that is outside the scope of this work, but some of these simplifications are broadly discussed in Section 4. M dwarf flares produce copious radio and millimeter emission, which is usually constrained to be gyrosynchrotron or synchrotron radiation from a population of mildly relativistic electrons (Güdel et al. 1996;Smith et al. 2005;Osten et al. 2005;MacGregor et al. 2018).All of the models in the grid are calculated for injected number-flux distributions with small ("hard") power-law indices in the range of δ = 2.5 − 4, which is consistent with the interpretation of available stellar flare constraints (e.g., Smith et al. 2005;Osten et al. 2007;MacGregor et al. 2020MacGregor et al. , 2021)).The peak injected beam energy flux densities (hereafter, we refer to the flux density as the "flux") above a low-energy cutoff, E c , span four orders of magnitude: 10 10 (F10), 10 11 (F11), 10 12 (F12), and 10 13 (F13) erg cm −2 s −1 .Our grid of M dwarf flare models includes a large range of low-energy cutoff values: E c = 17, 25,37,85,150,200,350, and 500 keV.
The large beam flux densities and high low-energy cutoffs in the grid are necessary to produce deep chromospheric heating in M-dwarf flares.The energy fluxes of 10 12 − 10 13 erg cm −2 s −1 are consistent with the largest values above E c ≈ 20 keV that have been inferred through the collisional thick-target model (CTTM) interpretation of some solar flare hard X-ray events (McClymont & Canfield 1986;Canfield et al. 1991;Wulser et al. 1992;Neidig et al. 1993;White et al. 2003;Krucker et al. 2011).Large, low-energy cutoffs, E c ≈ 70 − 120 keV, are also occasionally inferred in CTTM modeling (Holman et al. 2003;Warmuth et al. 2009).Furthermore, recent numerical advances (Arnold et al. 2021) generate δ = 3 electron-beam power-law distributions that contain a large percentage, up to ≈ 20 − 30%, of the ambient coronal density.Note that coronal electron densities in M dwarf atmospheres outside of flares (e.g., Osten et al. 2006;Liefke et al. 2008) are even larger than fiducial solar coronal values (n e ≈ 10 9 cm −3 ).Section 2.2 further discusses the energetic requirements of high-energy electron beams in stellar flares.
The time-evolution of the electron beam energy injection into a loop is poorly constrained in M dwarf flares, due in part to the lack of direct spatial resolution.However, the injection of high-levels of nonthermal energy into any given flare loop is probably very short in order to be consistent with shortest rise times of ≈ 2 − 10 s in optical, U -band, and UV observations (Moffett 1974;Robinson et al. 1995;Mathioudakis et al. 2006;Kowalski et al. 2016;MacGregor et al. 2021).We choose ∆t = 2.3 s as a representative timescale for injection into the RADYN model loop; this is sufficiently long compared to time-of-flight differences between the high-and low-energy electrons in the beams.The limitations of the applicability of these models due to the choice of short beam injection timescales are discussed further in Section 4. The grid consists of two types of time profiles for the injected beam flux: a constant injection rate and a ramp up/down injection rate.The constant injection models are helpful for isolating effects due solely to the atmospheric evolution, and the ramp up/down evolution is physically motivated by the pulsed-injection model of solar flare hard X-ray bursts (Aschwanden 2004).We now discuss our application of pulsed-beam injection to stellar flares.

Pulsed Beam Injection
There are a large number of assumptions that are needed to simulate the RHD of a stellar flare.Although neither the magnetic field dynamics nor the particle acceleration process are included in our RHD simulations, we relate some of the assumptions in our grid to a larger geometrical context.This provides some insight into the applicability and limitations of 1D RHD flare loop modeling.We draw on state-of-the-art advancements from models of magnetic field geometries in solar flares in order to further motivate the extreme ends of electron beam properties within the grid.
There are many observed timescales of nonthermal hard X-ray variations in solar flares, from tens of milliseconds to tens of seconds (Kiplinger et al. 1984).However, it is not yet clear how each of the respective beam populations contribute to the total thermal response of the chromosphere.Indeed, phenomenological power-law, power-spectra of frequencies adequately explain much of the temporal variations in solar and stellar flares (Inglis et al. 2015).The slow and fast components of solar flare light curves have also been attributed to a "trap-plus-precipitation" model of hard X-rays (Melrose & Brown 1976;Aschwanden 1998).At high temporal resolution, ∆t ≈ 4 − 20 s bursts (Aschwanden 1998;Xu et al. 2006;Penn et al. 2016;Rubio da Costa et al. 2016;Kawate et al. 2016;Collier et al. 2023) and a long-duration, smoothly varying component are often present in hard X-ray and optical solar flares.Even shorter, ∆t ≈ 0.1 − 3 s, bursts are observed in Hα and hard X-ray solar flare kernels (Kiplinger et al. 1984;Aschwanden et al. 1995;Aschwanden 1998;Aschwanden et al. 1998a;Aschwanden 2004;Qiu et al. 2012;Radziszewski et al. 2007Radziszewski et al. , 2011)).A large fraction of the hard X-ray emission is often contained in the smoothly varying component (Aschwanden et al. 1998b), which may originate from a population of electrons that are trapped in loops due to magnetic field convergence at the chromospheric footpoints.Trapping may also occur through some type of plasma-confining effect higher up in the coronal loop volume (e.g., Li et al. 2013;Egedal et al. 2015;Fleishman et al. 2022).In some solar flares, however, there is little evidence for large populations of trapped particles, and both the hard X-rays and gyrosynchrotron radiation can thus exhibit similarly pulsed temporal variations (White et al. 2003, see also Kundu et al. (2001) and Krucker et al. 2020).A wide variety of pulsed injection schemes have been used to drive models of solar flares (Emslie 1983;Lee & Gary 2000;Kašparová et al. 2007;Doyle et al. 2012;Rubio da Costa et al. 2016;Simões et al. 2017;Carlsson et al. 2023).
A step function is the simplest type of pulsed injection.However, we also desire a scheme that exhibits a gradual decrease in the flux injected at the top of a model RADYN loop.The dynamic pulsed injection model of Aschwanden (2004) v assumes adiabatic2 particle motion as the magnetic field lines retract after their reconnection.It is a parametrized model that accounts for the increase in magnetic field from the reconnection X-point to the previously-reconnected looptops that are downstream of the reconnection outflow.As the ratio of the magnetic field in the retracting loop to that in the underlying looptops decreases, the angular size of the loss cone increases, which allows particles with larger pitch angles to escape out of the trap.These particles precipitate into the lower regions of the loops (footpoints) and deposit the energy that drives the flare in the chromosphere.Footpoint convergence of the magnetic field causes further trapping of electrons before they reach the chromosphere, but this effect is not considered in this grid of models.We assume parameter values for the pulse full-width-at-half maximum (FWHM), t w = t FWHM = 2.3 s, the height of the fully contracted looptops, h L = 0.64 × 10 9 cm, and the initial height of the contracting fields, h Y (which is h X in Aschwanden 2004), where we assume that h Y = 2h L .The parameter q a is equal to t w , which is the width of the flux injected at the RADYN looptop, divided by a timescale for the injection of particles from the acceleration region.For the models here, we assume that q a = 1.28, which gives a cumulative number of electrons injected into a retracting loop from the acceleration region that is linearly proportional to t.
One can further obtain an estimate of the low-energy cutoff, E c , to the nonthermal distribution from recent numerical simulations (Haggerty et al. 2015) according to where B X is the magnetic field at the reconnection X-point.For B X = 1 kG and using the RADYN model loop apex electron density, Eq. 1 gives E c ≈ 125 keV.Here, it is assumed that the retraction occurs at the Alfven speed, c A,X , upstream of the reconnection outflow.Detailed calculations show deviations from an Alfvenic speed retraction process (e.g.Longcope et al. 2018), and energy transfer with various turbulent and shock phenomena may also occur farther downstream (Somov & Kosugi 1997;Li et al. 2012;Kontar et al. 2012;Egedal et al. 2015;Kontar et al. 2017;Ruan et al. 2023).Aside from our chosen value of t w = 2.3 s (see Section 4 for further discussion of the chosen timescale and loop sizescales), the values of the parameters from the pulsed injection model are not explored here.There are several ways one can also estimate the energy flux that powers a flare loop.We use a formula for the energy that is released through field line shortening (Longcope et al. 2016), which is This range of flare fluxes corresponds to a range of shear angles, ∆θ, which are proportional to the magnetic field components that are parallel to the active region polarity inversion line.These first-principle assessments, which draw from state-of-the-art numerical modeling, together demonstrate that high beam fluxes and large low-energy cutoffs are not within the realm of fantasy in stellar flares that occur in a region with ≈ kG fields in the low corona.
Figure 1 is an illustration of a possible magnetic field geometry in our application of the pulsed injection model.The reconnecting magnetic field geometry is similar to the analytic (Forbes & Priest 1995;Lin & Forbes 2000;Reeves & Forbes 2005, see also Reeves et al. 2007) and numerical (Stone et al. 2008(Stone et al. , 2020) ) models analyzed in Chen et al. (2020).Large magnetic fields are assumed in the reconnecting current sheet region, but the field strengths are not that much larger in the figure than those inferred in an off-limb solar flare (Chen et al. 2020).The implied photospheric fields are also very strong, but large fields of B ≈ 8 − 12 kG, may cover ≈ 10 − 20% of the surfaces of magnetically active M dwarfs (Shulyak et al. 2019;Reiners et al. 2022).Within this geometry, Chen et al. (2020) indicate regions that correspond to the "reconnecting current sheet", "contracting loops", and "flare arcade"3 In our envisioned scenario, the pulsed injection originates in the contracting loop regime rather than in the retracting loops within the vicinity of the reconnecting current sheet region (where the bulk of particle acceleration plausibly occurs; Arnold et al. 2021).One should keep in mind that there is actually very little agreement on the details in the dynamic processes and energy flow within this type of scenario (e.g., Fletcher & Hudson 2008), but many of the important features are consistent with state-of-the-art observations and 2.5/3D MHD modeling of eruptive events (Sui & Holman 2003;Sun et al. 2012;Liu et al. 2013;Longcope et al. 2018;Krucker et al. 2020;Chen et al. 2020;Rempel et al. 2023;Volpara et al. 2024).The Figure 1.Magnetic field vector lines from the 2D analytic loss-of-equilibrium model described in Reeves & Forbes (2005).The reconnection point ("X") is assumed to occur halfway up the current sheet where the magnetic field strength is ≈ 1 kG in this particular configuration.Three regions to the analogous diagram in Chen et al. (2020) are indicated.The loss cone opening angle increases as the fields contract toward the flare arcade below, allowing electrons with pitch angles within the cone limits to precipitate from the contracting region.The RADYN modeling simplifies the geometry by assuming a quarter-circle loop of constant area and a gravitational acceleration along the loop axis that varies according to gz = g⋆ sin(θ ′ ).In reality, the three regions are not as distinct as drawn here, and the dynamics cannot be shown accurately with a single snapshot (see Chen et al. 2020, for the evolution of fields in a similar geometrical configuration).The streamlines show the magnetic field direction, while the background color scale shows the magnetic flux density.
extensions to stellar coronae are even less established, including whether such models of mass ejections are appropriate (Crosley & Osten 2018;Alvarado-Gómez et al. 2018), possibly due to strong overarching fields as in confined solar flares (Thalmann et al. 2015).Nonetheless, Figure 1 is an adequate starting framework for linking some properties of the RADYN model spectra to particle acceleration and the action of strong, low-lying coronal magnetic fields that are expected in young M-dwarf flare stars.
An example of pulsed, ramp up/down beam energy injection into a RADYN model loop is shown in Figure 2 using the parameter t FWHM = 2.3 s of the Aschwanden (2004) model as described above.This pulse is normalized to a maximum injected flux density of 10 13 erg cm −2 s −1 , the low-energy cutoff is E c = 85 keV, and the power-law index is δ = 3.We introduce the nomenclature for unique model identification as mF13-85-3, where m means that the electron beam flux that is indicated by the F -number is a maximum over the pulsed, ramp up/down injection.Note that time-integrated energy in the pulse in the top panel is a factor of 2.5 larger than the peak flux.Alternatively, c indicates a constant injected flux pulse at t = 0 − 2.3 s and zero flux at t > 2.3 s (unless otherwise indicated in the model identifier, the injection timescale is ∆t = 2.3 s).In Figure 2, we also show the evolution of the detailed emergent radiative flux at a representative wavelength, λ = 3615 Å, that is just shortward of the ideal Balmer recombination limit (λ = 3646 Å).The evolution of the continuum flux indicates that the thermal radiation from the chromosphere responds over slightly longer timescales than the injected nonthermal energy (Allred 2005;Simões et al. 2017).

Emergent Radiative Flux Spectra
RADYN provides spectra from λ = 6 − 40, 000 Å in physical (cgs) units, which facilitates comparisons to fluxcalibrated stellar observations.Emergent radiative surface flux density spectra (erg cm −2 s −1 Å−1 ; hereafter, just "flux") are calculated within RADYN using a Gaussian quadrature sum of the emergent intensity at five µ values (µ = 0.05, 0.23, 0.50, 0.77, 0.95, where µ = cos(θ) and θ is the angle with respect to the atmospheric normal), and at 95  continuum wavelengths.We denote the flux at a specific continuum wavelength λ as "Cλ(t)".Throughout, a prime symbol indicates that the preflare flux is subtracted from the flux at time t.
A representative range of detailed continuum spectra from the grid is shown in Figure 3 for two models, mF13-85-3 and m5F11-25-4.Bound-free discontinuities are present in the spectra because dissolved-level bound-free opacities (Dappen et al. 1987;Hubeny et al. 1994;Hubeny & Mihalas 2014) are not included in RADYN calculations.Away from the bound-free edges, the optical color temperatures, the Balmer jump strengths, and the Balmer continuum slopes vary markedly between the mF13-85-3 and m5F11-25-4 models, which respectively heat column masses of log 10 m c /[g cm −2 ] ≈ −2 and −3 to gas temperatures of T ≳ 104 K.It is not unexpected that the major discriminatory power between these models occurs in the NUV and optical, while the far-ultraviolet (FUV) and X-ray and extremeultraviolet (XEUV) spectral regions are similar in shapes.The color temperatures that are inferred from the X-rays at λ < 60 Å, for example, are comparable between these two models (Section 3.1).Figure 3 also shows how these models improve upon phenomenological, T = 9000 K blackbody modeling of white-light stellar flares.
Background continuum opacity from elements other than those treated in detail (Section 2.1) have been included assuming LTE, with photoionization cross-section data compiled from The Opacity Project dataBASE (TOPBASE 4 ).The background opacities produce many of the faint bound-free edges in the spectra of Figure 3.The H − ion is not included in detail, but its population densities are self-consistent with the non-equilibrium electron densities in The bound-free continuum opacities that generate the ultraviolet edges (Vernazza et al. 1976) in these spectra are due to the following: C I λ < 1100 Å (2p 2 3 P ), S I λ < 1199 Å, Si I λ < 1525 Å (3p 2 3 P ), Si I λ < 1682 Å (3p 2 1 D), Al I λ < 2076 Å (3p 2 P 0 ), and Mg I λ < 2513 Å (3p 3 P 0 ).The detailed bound-free edges are due to He II λ < 227 Å and He I λ < 504 Å in addition to the typical H I edges.
order to give accurate bound-free and free-free opacities.The Ca I populations are also larger than in the lower solar atmosphere.However, only a model Ca II ion file (with a Ca III ionization stage) is available in the detailed calculations.Thus, we consider charge conservation using electrons from the LTE ionization balance of Ca I, II, and III, as for the other background elements, instead of with hydrogen and helium.The bound-free opacities of hydrogen from n = 6 − 8 are included assuming LTE populations relative to the continuum.
In the calculations, the population densities of H I are reduced by including LTE chemical dissociation equilibrium of H 2 (Vardya 1965;Tsuji 1973).In the pre-flare photosphere, the reduction of H I is about 50%.Although this affects the H opacity, the models do not have molecular opacities directly included in detail or in the background.Of course, M-dwarfs with T eff ≈ 3600 K show prominent molecular absorption bands of TiO, MgH, and CaH in their non-flare optical spectra (cf.Fig 4 of Bochanski et al. 2007).It is reasonable to assume that the heating rates in the low chromosphere in our flare models are too large, and the atmospheres heat up too quickly, for the energies of molecular dissociation to be a limiting term in the internal energy equation.
Higher density wavelength grids are used for the calculations across the emission lines of hydrogen, He I, He II, and Ca II that are treated in detail.The hydrogen Balmer line spectra (Hα, Hβ, Hγ) are modeled using the Dopplerconvolved, TB09+HM88 line profile functions (Vidal et al. 1970(Vidal et al. , 1971(Vidal et al. , 1973;;Hummer & Mihalas 1988;Tremblay & Bergeron 2009), which accurately capture the pressure broadening from ambient, thermal electrons and ions in the density regimes of flare chromospheres (Kowalski et al. 2017b(Kowalski et al. , 2022)).The pressure broadening of other Paschen and Lyman hydrogen lines is included as described in Kowalski et al. (2022).Within RADYN, the Balmer Hα, Hβ, and Hγ lines are calculated on frequency grids with 51, 31, and 31 points, respectively.The sampling in the far wings (|λ − λ 0 | ≳ 2 Å) is courser than around the rest wavelength.To mitigate systematic errors in the line-integrated fluxes of the Hγ line, we follow Kowalski (2022) and use a standard Feautrier solver to recalculate the emergent surface flux spectra on a 327 point wavelength grid with the frequency-independent, non-LTE source function from RADYN (see the discussion in Kowalski et al. 2022).We use a bilinear interpolation of log 10 ϕ α on a grid of T and n e , where α is the typical detuning parameter ∝ |λ − λ 0 |, and ϕ α is the normalized line profile function.This is followed by a four-point, third-order interpolation scheme on a grid of log 10 α and log 10 ϕ α (Vidal et al. 1973).
Several of the low-order Balmer, Paschen, and Brackett lines are available in the model output (see Appendix).We provide high-and low-resolution spectra of Hγ and low-resolution spectra of the Balmer Hβ and Hα lines.The diagnostic advantages of Balmer Hγ have been discussed elsewhere (Kowalski et al. 2022) (see also Section 3.2).Additionally, new line-shape calculations (Cho et al. 2022) that use the Xenomorph code (Gomez et al. 2016) have shown that (at high density) the static approximation to the perturbing ions in a plasma and the time-ordering of collisions are more important for the Hα line than have been previously demonstrated (Stehle 1994).However, the new line shape calculations of Hγ are essentially identical to the TB09+HM88 profiles that we use.Extending the Xenomorph calculations of Hα to lower electron densities, which are relevant for flare atmospheres, will be the subject of future work (Kowalski & Gomez, in prep).

XUV Backheating
X-ray and ultraviolet (XUV; λ = 1 − 2500 Å) backwarming is included following earlier RADYN models (Abbett & Hawley 1999;Allred et al. 2005;Kowalski et al. 2015) that use the methods of Gan & Fang (1990), Hawley & Fisher (1992), and Hawley & Fisher (1994).An incident radiation field is calculated from the Astrophysical Plasma Emission Code (APEC), which is part of AtomDB (Smith et al. 2001).Photoionization rates from CHIANTI are used for the calculations of volumetric heating rates.The capability to include the incident radiation field in the equation of radiative transfer (Allred et al. 2015) was not yet incorporated into the M-dwarf version of RADYN when this grid was created.However, the results of Allred et al. (2006) lead us to believe that differences among the treatments of XUV backheating are negligible in comparison to the large electron beam heating rates in this grid5 .

Optically Thin Radiative Loss Function
The optically thin radiative loss function, Λ cool , is an important term in the conservation equation for the internal energy in flare modeling.We use an optically thin loss function from CHIANTI version 8.0.1 (Dere et al. 1997;Del Zanna et al. 2015) to account for radiative losses at T = 10 4 − 10 8 K from species that are not treated in detail.We also exclude species for which the optically thin approximation is likely inaccurate in stellar flare atmospheres at low temperatures, T ≈ 10 4 K. Table 1 lists all of the atoms and ions that are excluded from the adopted thin loss function, which is shown in Figure 4(a).A standard cooling curve that includes these species predicts ≈ 10 3 greater radiative losses at the lowest temperature.A similar cooling curve at low temperatures was used in the RADYN flare simulations described in Kowalski et al. (2015).
The differences at high temperature, T ≳ 10 7 K, in Figure 4 are partially due to excluding the free-free cooling from H II, He II, and He III in the thin losses.In our experience with the grid of models, the thermal conduction, which is of the classical Spitzer form with a flux-limited saturation value (Smith & Auer 1980) in RADYN, is much more important than optically thin radiative losses in the energy equation within the model flare coronae.The treatment of free-free loss terms is even less important in the evolution of the lower, cooler atmospheric regions (the focus of our study) in flare models.We also note that the cooling curve is dependent on the assumed abundances (cf.Section 7.19 of the CHIANTI manual), which are discussed below.The differences among the loss functions at low temperature, T ≈ 10 4 K, originate primarily6 from including or excluding bound-bound transitions in Mg II, Fe II, and Si II ions, which we now further describe and justify.
The lower-atmospheric temperatures (Figure 4(b)) in flare simulations are dramatically different between two versions of the thin loss function with the same injected beam (mF13-500-3).Figures 4(c)-(d) show the contributions to the  internal energy equation in the calculations of panel (b).The thin cooling prevents a temperature increase above T ≈ 10 4 K at large column masses and low altitudes in the atmosphere.The emergent continuum flux spectra (not shown) are also quite different.Thus, our choice of excluding many transitions altogether (Table 1) from the loss curve at T = 10, 000 − 20, 000 K warrants further justification.We compare emergent LTE intensity spectra for a bright and a faint Fe II emission line at t = 1 s in the mF13-85-3 model.One calculation for each line includes the optical depth from the bound-bound transition ("optically thick calculation"), while another calculation excludes the boundbound contribution toward the total optical depth over the wavelengths of the line ("optically thin calculation").The bound-bound emissivity is considered in LTE in both calculations, as are bound-free emissivities and opacities.For the bright Fe II λ2632 emission line, the emergent wavelength-integrated emergent intensity is ≈ 200x greater than the optically thick calculation.For the faint Fe II λ2814.445line, the intensity is a factor of ≈ 10 − 15 larger than in the optically thick calculation.These differences are rather striking and suggest that excluding low-temperature emission lines altogether from the thin loss table, as we have done, is warranted until a more rigorous approach (e.g., as calculated in Carlsson & Leenaarts 2012, for quiet Sun conditions) is possible.Additional assumptions are that the abundances constantly remain at solar coronal levels (Schmelz et al. 2012) and the emissivities in the loss function are calculated from ionization equilibrium (in CHIANTI).These assumptions probably do not hold in flares (e.g., Colgan et al. 2008;Doyle et al. 2012;Warren 2014).We argue that the lack of uniform patterns in abundance changes among M-dwarf flares (e.g., Osten et al. 2005;Paudel et al. 2021), combined with the difficulties associated with metallicity diagnostics of M dwarfs (e.g., Mann et al. 2013), does not warrant varying the abundances away from the solar values at this time.Non-equilibrium rates of the elements not treated in detail may be important, but they must also be considered in the context of large densities that are produced in our models, which may in turn result in non-negligible optical depths in transition region lines (e.g., Mathioudakis et al. 1999;Kerr et al. 2019b).It would also not be surprising if non-LTE photo-ionization effects (e.g., Rathore & Carlsson 2015) for ions in the CHIANTI package are important in flares.All of these limitations are being actively investigated with the intent of improving upon them in future generations of stellar flare model grids.

Summary of the Grid
The model grid consists of 80 simulations that span a large range of heating rates and dynamical evolution.We organize the models into several groups.The main grouping consists of large low-energy cutoffs with a ramp up/down flux injection; the const group consists of large low-energy cutoff models with a constant flux injection; the Ec37 group consists of models with E c = 37 keV and hard power-law indices (Allred et al. 2005(Allred et al. , 2006;;Kowalski et al. 2015), as inferred through CTTM modeling of hard X-rays of the 2002 July 23 X-class flare (Holman et al. 2003;Ireland et al. 2013); the sol group consists of models with smaller low-energy cutoffs and softer (larger) power-law indices.The sol group also includes a model with a longer, ∆t = 15 s, injection of constant electron beam flux; this model is similar to the high-flux c15s-5F11-25-4 flare model in solar gravity that has been analyzed in Kowalski et al. (2017aKowalski et al. ( , 2022)).The auxiliary models include the following: a recalculation of a constant injection of F13 flux with E c = 37 keV and δ = 3 for 1.6 s (following Kowalski et al. 2016), a high-flux model with a very hard power-law index, δ = 2.5 (m2F12-37-2.5),and a recalculation of the soft (δ = 7), fully relativistic (E c = 500 keV) electron beam model (m2F12-500-7) that was used in Kowalski et al. (2017b) as a phenomenological model for a "Vega-like" photospheric spectrum that was reported in the decay phase of a white-light megaflare (Kowalski et al. 2013).Most models in the grid are calculated until t = 10 s, which allows the relaxation of each impulsively energized chromospheric region to be followed.It is also customary to test electron beam models against flare models that include energy transport through thermal conduction, given an ad hoc energy deposition in the corona (Fisher 1989;Reep et al. 2016;Kowalski et al. 2017a;Namekata et al. 2020).We do not specify the energy source, but it could plausibly represent heating across Petschek shocks in the reconnection process (Longcope et al. 2016).Note, however, that gas compression and other physics of magnetic field retraction are obviously not included in the RADYN simulations, which begin with a static, semi-circular loop.For the thermal models, flare heating is simulated as an ad hoc energy deposition of 125 erg cm −3 s −1 at z > 7.5 Mm.One model applies this coronal heating for ∆t = 2.3 s, and another for ∆t = 10 s.
We include five models from the Namekata et al. (2020) grid of models in a separate N+20 grouping; these models were calculated with the same version of the RADYN code.Three triangular (t) pulses have a maximum beam flux injection at t = 2 s.These are the tF12-37-3, tF12-37-5, and tF10-37-5 models.Namekata et al. (2020) also studied the response to triangular, thermal pulses with peak heating rates at t = 8 s at the coronal apex; the model identifications are thermal-tF12 and thermal-t5F10 with spatial integrations of the coronal heating rates that correspond to ≈ 10 12 and 5 × 10 10 erg cm −2 s −1 , respectively.We refer the reader to Namekata et al. (2020) for further details.Several N+20 models are somewhat redundant with models in other groupings.For example, the tF12-37-3 and the mF12-37-3 models are nearly the same except for a ∆t = 1 s between the times of the peak injected fluxes.These two models highlight uncertainties that are linked to small variations in the assumptions about the injected beam flux time profiles.The thermal-tF12 has a maximum heating rate of ≈ 1200 erg cm −3 s −1 , which is one order of magnitude larger than the three other thermal heating models.This model produces the hottest corona achieved among all models (Section 3.1).Thus, there is value in including these five representative models from the Namekata et al. ( 2020) grid.

MODEL GRID ANALYSIS
We divide the model analysis into three parts.First, we summarize the calculated quantities from the optical and near-ultraviolet grid spectra in Section 3.1.From these calculated quantities and the detailed emission line spectra, we demonstrate that the deep heating predictions of the model grid explain many of the observational constraints in the literature (Sections 3.2 -3.3).

Calculated Spectral Quantities
We follow Osten et al. (2016), Kowalski et al. (2017b), and Kowalski (2022) and calculate the temporally averaged radiative surface flux spectrum, F λ , from each model7 .Temporal averaging is a simple method to emulate many heating and cooling loops at different stages in their impulsive evolution on the star in a manner analogous to multi-thread modeling of solar flares (Hori et al. 1997;Reeves & Warren 2002;Warren 2006).The interpretation of F λ is that every model snapshot has an equal filling factor (area) in the flaring source over an exposure time, and all loops in the source are subject to the same heating pulse with uniformly distributed lags (cf. the superposition principle; Aschwanden & Alexander 2001).These are sufficient assumptions to first order, and we expand upon this interpretation in Section 3.3.
The main group of models is averaged over t = 0 − 10 s and the const group is averaged over t = 0 − 5 s.The thermal models are averaged over ∆t = 10 s, beginning when the conductive pulses reach the chromosphere.The flare-only, time-averaged radiative surface flux spectra are denoted as F ′ λ , and the flare-only radiative surface flux spectra at time t are indicated as F ′ λ (t).Note that Figure 3 displays F λ (t = 1s) for the mF13-85-3 and m5F11-25-4 beam models.
From the model spectra, we calculate several quantities to directly compare to stellar flare observations (Section 3.2).Table 2 lists several of the calculations, and others are described in the supplemental document that is referenced in the Appendix.The columns of Table 2 are the following: 1. Unique model identifier.
2. C4170 ′ : The value of F ′ λ at the continuum wavelength of λ = 4170 Å.Note that some small values are negative.The negative continuum fluxes can result from weak heating rates, which elevate the collisional rates in the chromosphere enough to diminish some of the background photospheric radiation (e.g., Allred et al. 2005;Kowalski et al. 2017a).

F ′
Hγ : continuum-subtracted F ′ λ integrated over the wavelengths, λ 1 to λ 2 , of the Balmer Hγ line; this is calculated from the temporally averaged surface flux spectrum (e.g., from t = 0 − 10 s), and the pre-flare value (in the first row) calculated from F λ (t = 0) has been subtracted.
7. T BB : blackbody temperature fit at the model continuum wavelengths over λ = 4000 − 4800 Å, following Kowalski et al. (2013); this is calculated as F ′ λ from the time-averaged surface flux spectrum.The temporal calculations of T BB (cols 5 and 6) give a sense by which the variation of the optical color temperatures differ from the color temperature of the average spectrum (col 7).A similar blackbody color temperature, T FcolorR (not listed in this table), is calculated from the ratio of F ′ λ at λ = 4170 Å and λ = 6010 Å. Spectra indicate that T BB ̸ = T FcolorR (Kowalski et al. 2016).In Section 3.2, the model values of T FcolorR are compared to comprehensive observational constraints over a broader wavelength range than fit with T BB .8. C3615 ′ /C4170 ′ : the Balmer jump ratio, which is the flare-only continuum flux at λ = 3615 Å divided by column 2; this is calculated from F ′ λ .
9. Hγ Eff.Width: the effective width (Kowalski et al. 2022) of the Hγ line; this is calculated from F λ after a linear fit to the nearby continuum flux on both sides of the emission line has been interpolated and subtracted: where λ 1 = 4312 Å and λ 2 = 4368 Å.
10. T plasma,max : the maximum value of the plasma temperature, considering all times in the model.The time at which the maximum plasma temperature occurs is t max .
11. T XEUV,max : the XEUV continuum thermal bremsstrahlung color temperature calculated from the atmospheric snapshot at the same time as in column 10.This is calculated from F λ (t = t max ).We integrate the emissivity from hydrogen and helium at plasma temperatures T > 1 MK.The color temperature is calculated from the ratio of the emergent optically thin flux at λ = 8 Å to that at λ = 35 Å.An X-ray color temperature (T X,max ; not listed) is calculated from the ratio of the emergent flux at λ = 0.5 Å to that at λ = 4 Å.
The values in Table 2 serve as a comprehensive overview of the observables from the models.Before turning to comparisons to observations of optical flare properties, a few comments about the coronal predictions (last two columns) are warranted.First, many of the models with large low-energy cutoffs do not generate enough Coulomb heating in the corona to produce a response at plasma temperatures T > 1 MK.This is certainly inconsistent with many stellar flare observations, which suggest temperatures in the range around T ≈ 15 MK to T ≫ 50 MK (e.g., Güdel 2004;Mullan et al. 2006;Liefke et al. 2010;Osten et al. 2005Osten et al. , 2010Osten et al. , 2016)).Second, a few of the models produce maximum plasma temperatures that are extremely high, T plasma,max ≈ 40 − 90 MK.However, the thermal bremsstrahlung temperatures that are calculated from XEUV continuum spectra over λ = 8 − 35 Å (last column) are much cooler.This systematic discrepancy occurs because the emergent radiation at these wavelengths originates mostly from the cooler, denser material at the base of the evaporation flows.The X-ray color temperatures (T X,max ; not listed) from 0.5 − 4 Å are much closer to the values of the maximum plasma temperatures, though in the hottest model flare coronae, there are still differences of T X,max − T plasma,max ≈ 5 − 20 MK.
The largest coronal plasma temperatures in our models are attained when the evaporation shock fronts reach the looptop.For example, in the mF13-37-3 model, a v z ≈ 1750 km s −1 evaporation reaches the looptop at t = 6 s, and then there is a jump in the apex temperature from T ≈ 50 MK to ≈ 90 MK by t = 7 s.The temperature explosions in the high corona in our models are facilitated by a reflecting upper boundary condition (which we chose to emulate a symmetrical flare loop).We analyze the terms that contribute to changes of the internal energy per gram, and we find that the rapid temperature increase is due to viscosity and compression in the large velocity gradient that occurs as flows from both sides of the loop collide.Notably, these effects have been previously discussed in detail in the high-flux hydrodynamic model of the coronal explosion of CN Leo (Schmitt et al. 2008).We also are able to follow the mF13-37-3 model to t = 20 s.We find that the evolution of the coronal densities and velocities is qualitatively very similar to Fig. 3 of Reeves et al. (2007) on much longer timescales over much larger loop lengths.The main difference is that a secondary chromospheric evaporation wave in our model begins at t ≈ 10 s due to the looptop heat conduction pulse traveling down the loop ahead of the flows from the other side.The secondary evaporation and the flow from the other side collide in the middle corona by about t = 14 s.By t = 20 s, most of the corona is still very dense, n e = 2.5 − 5 × 10 11 cm −3 and very hot, T ≈ 40 MK.It is reassuring that three entirely different numerical codes predict gross commonalities in the coronal evolution driven by large heating rates.
The impulsively-heated RADYN models in our grid are, however, not intended to include the late evolution of the corona long after the end of electron-beam energy injection, which we assume is very short within each loop in stellar flares (see Section 4 for further discussion).The optical and NUV radiation decrease to very faint levels by t ≈ 10 s (Figure 2), at which point the integrated coronal emission is still increasing in some models.Only a few of the models in our grid are extended to t > 10 s, but all are averaged over a short duration (∆t = 5 or 10 s) to provide the average quantities in Table 2.The cooling over the late evolution of the coronae may be further affected by transport physics (Bian et al. 2018;Zhu et al. 2018;Allred et al. 2022;Ashfield & Longcope 2023) that are not considered in this generation of models.
Table 2.The heating in the mF10-350-3 and mF10-500-3 models result in an attenuation of the pre-flare photospheric flux, thus causing negative continuum fluxes and unreliable blackbody fits as well.

Comparison to Global Trends
Kowalski et al. ( 2013) presented global trends in optical spectroscopic observations of dMe flares.The continuum constraints are generally interpreted as evidence for white-light radiating sources that cannot be explained by optically thin thermal recombination radiation or by a superhot free-free source (see also Hawley & Fisher 1992).The spectra include the Balmer jump and extend to λ ≈ 9000 Å in the largest events.The trends have been supplemented with spectra of smaller dMe events in Kowalski et al. (2016) and Kowalski et al. (2019b), and with high-time resolution ULTRACAM (Dhillon et al. 2007) photometry in Kowalski et al. (2016).Figures 5(a-c) compare several calculated quantities (Table 2) from the main group of models to the impulsive phase observations from these studies.
The Balmer jump ratio correlates with the Hγ line-to-continuum ratio in both the models and the observations (Figure 5(a)), but the models predict steeper slopes.The mF12 and m2F12 models satisfactorily reproduce some of the smallest Balmer jumps, which tend to occur in the observations at the peaks of the most impulsive flares.However, these beam fluxes do not account for the relationship between the small Balmer jumps and the optical color temperatures of T ≈ 10, 000 K. In Figure 5(b), the Balmer jump ratios are shown against the optical color temperatures, T FcolorR , that are calculated from the ratio of the flare-only continuum fluxes at λ = 4170 Å and λ = 6010 Å (Section 3.1).The main mF13 models satisfactorily reproduce the range of the largest optical continuum flux ratios that extend to the flare spectra with the Balmer continuum in absorption (open star symbols; Kowalski et al. 2013).The large open circle at (x, y) = (2.3,1.8) is calculated from the peak-phase spectrum8 of the Great Flare of AD Leo (Hawley & Pettersen 1991), which falls close to the prediction of the mF13-150-3 model (Kowalski 2022).
Figure 5(c) contains additional constraints on Balmer jump ratios and optical blackbody color temperatures at the peak phases of a sample of flares with ULTRACAM data (Kowalski et al. 2016).The two largest events in this sample correspond to the IF1 and IF3 events, whose peak colors straddle the prediction of the mF13-85-3 model (Kowalski 2023).Parameterized models of dense chromospheric condensations (purple dashed line; Kowalski & Allred 2018) reproduce the general trend in the peak-flare data.However, the models at the bottom right of the track generate extremely broad Balmer emission line profiles that are inconsistent with observations (Kowalski et al. 2017b;Kowalski 2022).The chromospheric condensation models also do not extend to the lower right end of the observed continuum colors, which are accounted for by the main F13 models in the RADYN grid.

Hydrogen Broadening from Flare Models with Deep Heating
The superposition of a lower-flux beam model with one of the large-E c , high-flux beam models brings many of the discrepancies (e.g., line-to-continuum ratios) further in line with the observations in Figure 5.A linear, least-squares regression analysis with our grid of models was presented in Kowalski (2022), who fit the broadband FUV-to-optical continuum shapes and the broad hydrogen emission lines in the Great Flare of AD Leo (Hawley & Pettersen 1991).In this section, we model the optical spectra of a remarkable white-light stellar flare observed at much higher resolving power, which facilitates rigorous tests of the predictions of the large optical depths and electron densities produced in the large-E c , high-flux models.
An energetic flare from the dM4.5e star YZ CMi was observed on 2008-Jan-21 with the 2.03 m Telescope Bernard Lyot at the Pic du Midi Observatory.The cross-dispersed NARVAL spectropolarimeter (Silvester et al. 2012) was employed with a 2. ′′ 8 slit and 1200 s exposure times.The linear dispersion at λ ≈ 4300 Å is 0.027 Å pix −1 .A series of four consecutive Stokes I spectra covering λ = 3700 Å to 10, 000 Å were obtained from 22:28 to 23:32 UT at an airmass around 1.3.The spectra are not normalized by the continuum, but they are divided by the large-scale sensitivity of a flat field.The data were first presented in the context of other stellar flares with blueshifts in the Balmer lines (Vida et al. 2019).We retrieved these spectra from the PolarBase archive (Donati et al. 1997;Petit et al. 2014), and we uniformly applied a wavelength shift of ∆λ = −0.35Å, which was required to align the pre-flare Ca II K emission line centroids in the model and observations.The ULTRACAM peak-phase colors are annotated for the large events that are referred to as "impulsive flare 1" (IF1) and "impulsive flare 3" (IF3) in Kowalski et al. (2016).The HST-1 and HST-2 data are added from Kowalski et al. (2019b), the F1/IF11 and F2/IF4 flares are added from Kowalski et al. (2013), and the rest of the spectral data in (a) and (b) are obtained from Kowalski et al. (2013) and Hawley & Pettersen (1991).The open star symbols are calculated from "newly-formed" flare flux spectra during secondary events, and the large open circles are obtained from spectra during the impulsive phase of the Great Flare of AD Leo (Hawley & Pettersen 1991).Blackbody color temperatures are indicated in (b) and (c), and the parameterized chromospheric condensation (CC) models from Kowalski & Allred (2018) are included in (c).Error bars on the data indicate marginal 68% uncertainties.In panel (c), the systematic uncertainties are included in the gray error bars, while the black error bars are the statistical uncertainties only (note that the statistical uncertainties are very small for the IF1 and IF3 events).In (b) and (c), flares with only Balmer jump calculations are shown in the right panels in the margins.
Spectra around Balmer Hγ and Hϵ are displayed in Figure 6(a) and Figure 6(b), respectively.The sequential exposures include a pre-flare spectrum, a spectrum with a powerful response throughout the optical continuum and highly broadened hydrogen Balmer lines, and a spectrum with larger line-to-continuum ratios and less broadened hydrogen lines.We refer to the two flare spectra as the "continuum peak" and the "continuum decay" spectra.From the continuum peak to the continuum decay spectra, the value of F ′ Hγ /C4170 ′ increases from ≈ 10 to ≈ 100, and the effective width (Eq. 3) of the Hγ line decreases from ∆λ eff = 4.5 Å to 2.5 Å.Although broadband photometry data are not available to provide context for the spectra, it is conceivable that the long exposure times integrated over flare phases similar to the echelle observations of the giant flare that was analyzed in Fuhrmeister et al. (2008).It is also possible that the Balmer line flux peaks after the continuum, an effect that occurs in some events (Hawley & Pettersen 1991;García-Alvarez et al. 2002;Namekata et al. 2020).
In Figure 6(b), the Hϵ line broadens far more than the Ca II H and K lines.Such differences are well understood to be the result of enhanced charged-particle pressure (Stark) broadening of the hydrogen lines in high-density flare chromospheres (Hawley & Pettersen 1991;Allred et al. 2006;Paulson et al. 2006;Kowalski et al. 2022).We compare to the predictions of the mF13-85-3 model around the Hγ line and the nearby continuum flux (Figure 6(c)).The model spectrum, F λ , at nλ = 327 wavelength points (Section 2.3) is multiplied by a scale factor so that it matches the approximate continuum count fluxes around λ ≈ 4320 and λ ≈ 4360 Å.The comparison demonstrates that the shape of the wings is reproduced, which is remarkable because the model spectrum is a temporal average that consists of the evolution of the model atmospheric densities, temperatures, and heating rates (Section 3.1).The Balmer Hβ and Hδ model predictions are also within the constraints of the observed wings of these lines (not shown).The Hγ profile from the tF12-37-3 model in comparison demonstrates that lower energy beam heating does not reproduce the relative fluxes in the wings and nearby continuum.This model is a much better match to the continuum decay phase spectrum (Figure 6(e)).
The goal of modeling flares with high-flux electron beams and large values of E c is to reproduce the optical continuum properties in M dwarf flares.Then, we assess the extent to which such models agree with emission line spectra.Clearly, the electron beam heating models in Figure 6 do not account for the narrow, core emission in Hγ.The discrepancy is especially evident for the comparison of the continuum peak phase observation to the mF13-85-3 model (Figure 6(c)), which has central dip at line center.The reason for the central dip will be clarified in Section 3.3.1.A least-squares fit of a Gaussian to the extra amount of line core emission in the observations gives best-fit FWHM values of 0.7 − 0.8 Å, and the decay phase spectrum requires a Gaussian with a factor of ≈ 2.5 greater flux at Earth.The total models, which each consist of a primary beam model spectrum and a Gaussian, are shown in Figures 6(d,f).In Kowalski (2022), a similar discrepancy in the Balmer lines was rectified by superposing a lower beam flux model with one of the mF13-85-3, mF13-150-3, or mF13-500-3 models.The interpretation is that lower and higher-energy electron beams impact different areas (see also Cram & Woods 1982;Kowalski et al. 2010), which is in line with detailed analyses of optical (Neidig et al. 1993) and near-ultraviolet (Kowalski et al. 2017b) solar flare kernels at high spatial resolution.However, only a few of our models (mF10-150-3 and mF10-17-3) produce Hγ line profiles with FWHM values that are as small as the Gaussian fits suggest.To account for the narrow line core flux at Earth, the area of the region emitting the mF10-17-3 spectrum would have to be a factor of 6.5 larger than the source area with the higher beam flux heating (mF13-85-3) that accounts for the continuum and the far wing radiation in the impulsive phase.

Line Formation Properties
The mF13-85-3 and tF12-37-3 models satisfactorily explain the Hγ broadening at |λ − λ 0 | ≳ 0.5 Å in the continuum peak and continuum decay phase flare spectra, respectively.In this section, we compare the formation of the Hγ lines in these models.We analyze the tF12-37-3 model, instead of the very similar mF12-37-3 model, because Namekata et al. (2020) used this simulation to understand the formation of the broad Hα lines in a superflare on AD Leo.
Figure 7 shows contribution functions to the emergent radiative intensity across the hydrogen Balmer γ line.The preflare formation of the model Balmer γ intensity (Figure 7(a)) is displayed for the same wavelength and color scales as for the flare simulations (Figure 7(b-c)) in order to emphasize the amounts of the pressure broadening increase in the high-flux beam models.The thin dashed (blue) lines indicate the heights corresponding to 5% and 95% of the emergent intensity, and the thick dashed (blue) lines show the τ (λ) = 1 heights.The margin plots in each panel show the electron density, gas temperature, and the source-function equivalent temperature.Comparisons of these two measures of temperature indicate where the line and nearby continuum formation is close to local thermodynamic equilibrium (LTE) conditions.In both models, the line core formation deviates from LTE.In the mF13-85-3 model, the line core has a central "absorption" because the source function equivalent temperature decreases over a ∆z ≈ 100 km height range following its departure from the gas temperature curve.The line core is formed over much lower electron densities than the wings (Kowalski et al. 2017b;Namekata et al. 2020).In the wings and nearby continuum, the mF13-85-3 model becomes optically thick at electron densities of n e ≈ 1 − 2.5 × 10 15 cm −3 at heights well above the heights that correspond to τ (λ) = 1 in the preflare and in the tF12-37-3 models.The formation of the Balmer γ line wings in the mF13-85-3 model is very close to LTE.The formation of the hydrogen Paschen β (n = 3 to n = 5) transition is shown in Figure 7(d) for the mF13-85-3 model at t = 1 s.This emission line shares the same upper level as Balmer γ and thus should experience approximately the same amount of Stark broadening.This transition is also more optically thin than the Balmer γ transition, which ostensibly allows one to probe deeper into the atmosphere (Kowalski et al. 2022), where there are larger electron densities in the mF13-85-3 model.Figure 8(top) shows the emergent continuum intensity at the wavelengths calculated in detail.The contribution functions and optical depths at continuum wavelengths are obtained from the model output following the method of Kowalski et al. (2017a).Figure 8(bottom) indicates the column mass at which τ = 1 occurs at selected wavelengths in the top panel.The continuum wavelengths (λ ≈ 12150 Å) near Paschen β are more optically thick at a given height in the atmosphere than the blue continuum wavelengths (λ ≈ 4300 Å) around Balmer γ.This pushes the formation of the Paschen β line to smaller column masses in the atmosphere with higher temperatures and lower electron densities (at this particular time in this model).The source function equivalent temperature for the Paschen β line is closer to the gas temperature, indicating that it is formed closer to LTE than Balmer γ, as expected.
To accurately interpret electron densities from the pressure broadening of spectral lines in flares it is crucial to account for the amount of broadening due to optical depths.Optical depth effects in flare models can greatly enhance the broadening beyond the expectation from the largest electron density in the flare atmosphere (e.g., Johns-Krull et al. 1997;Kowalski et al. 2022).How well does the amount of broadening in hydrogen Balmer lines relate to the actual electron densities achieved (cf. the margin plot in Figure 7(c)) in the mF13-85-3 flare atmosphere? Figure 9 compares continuum-subtracted, peak-normalized spectra of Balmer γ.Here, we use the spectra calculated at nλ = 327 wavelength points (Section 2.3).The mF13-85-3 model is compared to several spectra that are calculated from optically thin slabs with a uniform electron density.The two broadest slab model spectra are calculated with n e = 2 × 10 15 cm −3 and n e = 5 × 10 15 cm −3 , respectively.The emergent intensity spectrum at t = 1 s from the mF13-85-3 model is much broader than these, even though the electron densities in the atmosphere are far below n e = 5 × 10 15 cm −3 (Figure 7(c)).This discrepancy is due to the large optical depths (τ ≈ 1) over the heights where the line is formed.The flux spectrum from the tF12-37-3 model at maximum beam injection (t = 2 s) is as broad as the n e = 2 × 10 15 cm −3 slab calculation, but the actual electron densities where the line is formed are an order of magnitude smaller9 (Figure 7(b)).
It is also valuable to compare the time-averaged flux profile (F λ ) that fits the Balmer wings in the YZ CMi flare (Figure 6) to the intensity and flux spectra (F λ (t = 1 s)) at the time of maximum beam heating.The time-averaged flux spectrum from the mF13-85-3 model is narrower than the t = 1 s flux spectrum because the chromosphere experiences the deepest heating and largest ionization around the time of the maximum injected beam flux at t = 1 s (Figure 2).The t = 1 s flux and intensity spectra exhibit minor differences in the wing shapes and larger differences in the core10 .

DISCUSSION
We have calculated a grid of 80 RHD stellar flare models with the RADYN code.The RADYN models extend to large electron beam fluxes in order to reproduce the empirical optical and near-ultraviolet continuum properties in M-dwarf flares.Small electron beam fluxes are included too; these do not reproduce any property to which we compare unless the low-energy cutoff is very large.For large low-energy cutoffs and large beam fluxes, scaling to the optical (or NUV continuum, if available) reproduces the T ≳ 10 4 K blackbody color temperatures and small Balmer jumps in the impulsive phases of several well-studied M dwarf flares.We have shown that if a continuum prediction from one such model is scaled to the observed optical continuum flux in a large dMe flare, the wing shape of the highly broadened Hγ line is remarkably well-reproduced.) velocities that correspond to the wavelength detuning and to the gas velocity (yellow line).Margin plots show the source-function equivalent temperatures (Sν T ) and the gas temperatures (bottom axes) compared to the ambient (thermal) electron densities (top axes).Note the different scales among the electron density axes.As noted in Kowalski (2022), the main mF13 models produce non-negligible chromospheric upflows, but these are not apparent on the velocity scale shown here.
The large low-energy cutoffs in our RHD modeling approach imply that large heating rates are needed over low chromospheric heights in stellar flares.One possible source of deep heating that is not often considered arises from MeV proton beams (Rieger et al. 1996;Emslie et al. 1998;Procházka et al. 2018;Allred et al. 2020;Sadykov et al. 2024).Another possibility is that time-dependent transport effects (Kontar et al. 2012) generate an enhancement of E ≳ 100 keV electrons in stellar flares.The enhancement produces a similar local heating maximum at m c ≈ 0.01 g cm −2 as in the E c = 85 keV models with δ = 3, and the temperatures around this column mass exceed T = 10 4 K if the energy flux is large enough (Kowalski 2023).We note that in some solar flares, there is growing evidence that deeper heating is needed beyond that from single-power law electron beams that are inferred from the standard interpretation of nonthermal hard X-rays (Martínez Oliveros et al. 2012;Warren 2014;Procházka et al. 2017;Kowalski et al. 2019a).The grid of electron beam models that are listed in Table 2, which we refer to as the v1.0 models, simulate flare heating (Q flare ) from binary, beam-plasma Coulomb collisions only.Collective forces on the beam particles have not yet been included in our M dwarf flare models.The effects of return current energy loss and magnetic mirroring introduce a large number of additional free parameters (e.g., turbulent scattering effects and the location of magnetic field convergence; Murphy et al. 2007;Kontar et al. 2014).In the state-of-the-art computational framework (Allred et al. 2020), the return current and magnetic forces are limited to the steady-state solution of the distribution function.This assumption may be valid, but more work is needed to determine how time-dependent transport effects (e.g., Kontar et al. 2012) alter the steady-state distribution.Additionally, nonthermal runaways may be generated in the background plasma due to large return current electric fields (Holman 1985;Alaoui et al. 2021).Nonetheless, it is expected that the summative effects of the return current electric field would prevent much of the beam energy in our models from reaching the low chromosphere, thus precluding any electron beam model from achieving the deep heating that explains the stellar observations.The v1.0 model grid has a large number of valuable applications.These models have already been leveraged in a variety of ways: studies of the merging of the Balmer series (Kowalski et al. 2017b) for NUV exoplanet habitability assessments (Kowalski 2022), the energy partition in NUV and optical photometry (Brasseur et al. 2023), infrared flaring in JWST spectra and the response of the He I λ10830 line (Howard et al. 2023), solar-stellar comparisons (Monson et al. 2024), and the relations among gas, color, and radiation temperatures in flare models (Kowalski 2023).Many further applications of the grid models are possible.For example, the Balmer line profiles in the deep heating models can be comprehensively tested against a recently published large sample of M dwarf flares with high time-resolution and broadband photometry (Notsu et al. 2024).Other uses include the following: FUV and Balmer continuum slopes can be compared to observations (Loyd et al. 2018a;Froning et al. 2019;Kowalski et al. 2019b;Feinstein et al. 2022;Chavali et al. 2022), optically thin coronal and transition region emission lines can be calculated (Allred et al. 2006) from the early stages of chromospheric evaporation, the relative light-curve timescales between the FUV continuum and the U -band (Hawley et al. 2003) can be investigated, more gradual-type white-light events with larger Balmer jumps (Figure 5) can be modeled in their impulsive phases (Kowalski et al. 2019b), and the broadening of the Paschen lines (Fuhrmeister et al. 2008;Schmidt et al. 2012;Kanodia et al. 2022;Howard et al. 2023)  Flux: t = 1.0 s (mF13-85) I λ (µ = 0.95): t = 1.0 s (mF13-85) Flux: t = 2.0 s (tF12-37) Figure 9. Model spectra around the hydrogen Balmer γ line in the two RHD models (mF13-85-3 and tF12-37-3) from Figure 6.The y-axis is the integrand of Eq. 3: the nearby continuum has been subtracted, and the resultant spectra are normalized to the maximum over the line.Doppler-broadened TB09+HM88 line profile functions, which are equivalent to optically thin static slab calculations when peak-normalized, are shown as progressively broader dashed yellow spectra for electron densities of ne = 5 × 10 13 , 10 14 , 5 × 10 14 , 10 15 , 2 × 10 15 , and 5 × 10 15 cm −3 at T = 10 4 K.
The expansive range of input heating rates in the grid facilitates semi-empirical modeling of stellar flares with two beam components (e.g., a linear superposition of model spectra through maximum likelihood parameter estimation; Kowalski 2022).We applied the concept to a large dMe flare with archival echelle observations of the hydrogen Balmer γ line (Section 3.3).The impulsive-phase and gradual-phase blue-optical continuum levels together with the Hγ wing shapes are well-reproduced by high-energy beam models from the grid.In solar flare multithread models that use thermal conduction pulses (Warren 2006), the post-impulsive phase is explained using a superposition of equal-length loops, and the energy flux deposited per loop decreases as the flare progresses.Recently, Rempel et al. 2023 showed the evolution of thermally conducted energy fluxes into the footpoints in one of their 3D solar flare simulations.Although a comparison to 1D multithread models (e.g., Warren 2006) was not discussed, it appears that a similar progression occurs with higher average fluxes into the footpoints during the earlier evolution.The differences in the energy fluxes of high-energy electron beams between the mF13-85-3 model and the tF12-37-3 represents an analogous change from peak to decay in the YZ CMi flare studied here.These two primary beam models are remarkably different, which is justified in accounting for the striking differences in the spectral properties (line-to-continuum ratios, line wing shapes, effective widths) in the observations.The modeling of the continuum and wing in each of the observations further suggests an evolution in the average energy per beam particle from ≈ 170 keV to ≈ 75 keV.As the flare reconnection progresses higher into the corona (which occurs in loss-of-equilibrium models, such as in Figure 1), the flux and average energy per particle may sensibly decrease with the decreasing strength of the magnetic field, which must overcome the increasing length of the current sheet (Lin & Forbes 2000;Reeves & Forbes 2005).All of this supports the idea that the evolution of spatially-integrated stellar flare light curves is not solely driven by changing ribbon areas or cooling times of monolithic loops.
Nonthermal electrons leave signatures in the microwaves well into the gradual decay phase of the white-light in stellar flares (Osten et al. 2005) and after the peak of the hard X-rays in solar flares (Krucker et al. 2020).So, one could imagine that there is a relationship between gyrosynchrotron-emitting electrons and chromospheric heating after the impulsive phase.Multi-wavelength analyses of the empirical optical-radio connection in stellar flares (Osten et al. 2016;Tristan et al. 2023) would be helpful because the high-flux model (tF12-37-3) component in the decay phase of the continuum radiation implies that accelerated electrons are important in atmospheric heating even at late times (see also the RADYN multithread modeling of the decay phase of a megaflare in Kowalski et al. 2017b).However, previous spectral observations of the decay phase in stellar flares show that the Balmer jump ratios are even smaller than the prediction of the tF12-37-3 model (cf.Fig 21 of Kowalski et al. 2013).
We leave the origin of the narrower and fainter line-core model components in each of the impulsive and gradual phase spectra in Figure 6 to future work when other lines can be incorporated into a comprehensive statistical analysis with appropriate priors on the Balmer jump ratios (which are not available in the echelle spectra from NARVAL).The origin of the late-phase peak in Ca II K flux in stellar flares (Hawley & Pettersen 1991;García-Alvarez et al. 2002;Kowalski et al. 2013) is yet unexplained (Kowalski 2022), and detailed comparisons to the non-hydrogen emission lines in the observations (e.g., Figure 6(b)) would thus be valuable.However, there may be unavoidable limitations in the 1D grid of models that we present here.In the Reale et al. (1997) paradigm for modeling the soft X-ray emission of stellar flares, the inferred loop lengths are much larger than 10 Mm (Güdel et al. 2004;Reale et al. 2004;Osten et al. 2010Osten et al. , 2016)).In our grid, the short durations (∆t ≈ 10 s) over which the coronal plasma in a small loop are followed may be partly responsible for the discrepancies around the emission line cores.For example, coronal irradiation from a large volume of late-phase flare loops might be able to sustain Ca II at high levels into the gradual decay phase of the optical and NUV continuum radiation (Hawley & Fisher 1992).Alternatively, we speculate that a temperature increase due to betatron heating / acceleration (Veronig et al. 2006;Birn et al. 2017) during a gradual phase of shrinkage of each loop might contribute to heating that is required to raise the source function around formation layers of the Balmer and Ca II line cores in high-energy, electron-beam models.
The models in our grid are intended to reproduce the bulk of the optical and near-ultraviolet radiation in the impulsive-phase of a variety of dMe flare energies, peak amplitudes, and durations (Figure 5); they are particularly well-suited for the rise and peak phases of flares that exhibit evidence for hot blackbody-like optical continuum radiation and relatively small Balmer jump ratios (e.g., Kowalski et al. 2013;Kowalski 2023).In the early hard X-ray and white-light impulsive phase of solar flares, brightenings elongate the size of ribbons, which tend to closely straddle the magnetic polarity inversion lines (e.g., Qiu et al. 2010Qiu et al. , 2017;;Kazachenko et al. 2022;Tamburri et al. 2024).The projected areas of loop arcades can grow substantially over time (e.g., Aschwanden 2012) as the ribbons spread apart.Thus, chromospheric ribbons in the early phase are connected by the shortest loops within a flare (in the absence of large shear initially).The length chosen for a model loop in the early phase of a flare does not have a large effect on the Coulomb energy losses from representative beam electrons in transit to the chromosphere: for an electron density of n e = 5 × 10 10 cm −3 , the stopping distance (Emslie 1978) for a beam electron with initial kinetic energy E 0 = 75 keV is 170 Mm, and the stopping distance for a E 0 = 170 keV beam electron is 760 Mm.The temporal averaging of radiative flux spectra in models with short beam-injection timescales (∆t ∼ 2 s) approximates the superposition of many sequentially heated chromospheric kernels in stellar flares, which presumably also have an early phase of ribbon elongation and relatively small loop lengths connecting conjugate footpoints.The increasing areas that are inferred through optical blackbody fitting in the rise phase of stellar flares generally supports this analogy (e.g., Kowalski et al. 2013).
We can further justify the short timescales of beam heating by comparing to a recent high-resolution analysis of a solar flare.Graham et al. (2020) calculated that significant amounts of newly brightened chromospheric areas in a ≈ 10 32 erg solar flare occur within an observational cadence of 19 s, which is a factor of ≈ 10 − 20 shorter duration than the total hard X-ray impulsive phase duration.Crudely extending this relation to the U -band impulsive phase timescales (which we take as the t 1/2 values in Table 6 of Kowalski et al. (2013) and Table 6 of Kowalski et al. 2016) of the stellar flares in Figure 5 suggests that it is reasonable to assume that beam heating lasts for only several seconds in each impulsive-phase loop.Better assessments that leverage unique capabilities of the ROSA instrument (Jess et al. 2010) are forthcoming (Kowalski et al 2024, in preparation).

CONCLUSIONS
There has been an explosion of stellar flare observations in the last few decades, but modeling has generally been limited to static, semi-empirical atmospheric calculations or to slabs.When confronted with white-light constraints, time-dependent models (including 2D and 3D magnetohydrodynamic simulations) of X-ray flares do not make predictions at all.The RADYN M-dwarf flare model grid is a modern extension of the static, semi-empirical atmospheric flare models of Cram & Woods (1982), which (to our knowledge) were the first models to produce blackbody temperatures of T ≈ 14, 000 K and broad hydrogen Balmer emission lines.The RADYN models include vertical atmospheric heterogeneity and optical depths that are self-consistent with the equations of RHD in response to Coulomb heating from a power-law distribution of electrons.We analyze observations using superpositions of time-dependent models, which emulate lateral heterogeneity over the spatially unresolved flare regions on other stars.The models reproduce the full range of Balmer jump properties in the observations of the impulsive phase, and we discuss how they are also useful for understanding the role of nonthermal particle heating in the gradual decay phases of light curves of the optical and NUV continuum radiation.Most notably, the hydrogen Balmer γ wing shape that results from very high rates of deep atmospheric heating agrees with symmetric broadening in echelle flare observations (Figure 6(c)), thereby providing a plausible physical interpretation of impulsive-phase, Balmer line "broad components" that are often phenomenologically modeled with a wide Gaussian (e.g., Doyle et al. 1988;Fuhrmeister et al. 2008).In the impulsive phase models, the Balmer wings are formed over large atmospheric electron densities and optical depths, from which the observed blue continuum properties also follow.We conclude that models with deep atmospheric heating provide a robust paradigm that improves upon the ubiquitous T ≈ 9000 K blackbody interpretation of stellar flare white-light continuum radiation.

A. MODEL GRID ACCESS
In this appendix, we describe how to access the model grid output.FITS binary tables contain the detailed spectra and atmospheric variables at ∆t = 0.2 s (in most cases).Average spectral quantities, such as those in Table 2, are contained in modelvals.tave.fits.We distribute a supplementary Jupyter notebook that demonstrates how to load the outputs of the models.The Python routines and this notebook can be obtained through PyPI by issuing the following command: pip install radyn_xtools.We recommend consulting the documentation on the PyPI webpage for radyn_xtools about setting up a standalone conda environment.The user is also advised to consult the RADYN documentation and IDL routines that are linked to the F-CHROMA solar flare model grid paper (Carlsson et al. 2023).The original Common Data Format files are large in size, but they are available upon request to the corresponding author.
The .fits files along with the documentation that describes the installation of radyn xtools and the contents of the files are available at the Zenodo repository, https://doi.org/10.5281/zenodo.10929515.

Figure 2 .
Figure 2. Evolution of electron beam flux injected at the RADYN model loop top for two types of pulses in the model grid: (Top) pulsed ramp up/down precipitation flux from the formulae in Aschwanden (2004) in the mF13-85-3 model at ∆t = 0.2 s; (Bottom) a pulsed step-function injection in the cF13-85-3 model at ∆t = 0.05 s.The evolution of the emergent radiative flux at a Balmer continuum wavelength is shown in each panel for comparison.

Figure 3 .
Figure 3. Radiative surface flux spectra from a representative range of heating rates covered by the flare models in the grid.The pre-flare surface flux spectrum from RADYN and a T = 9000 K blackbody function (scaled by a factor of 2.6) are shown for comparison.The bound-free continuum opacities that generate the ultraviolet edges(Vernazza et al. 1976) in these spectra are due to the following: C I λ < 1100 Å (2p 2 3 P ), S I λ < 1199 Å, Si I λ < 1525 Å (3p 2 3 P ), Si I λ < 1682 Å (3p 2 1 D), Al I λ < 2076 Å (3p 2 P 0 ), and Mg I λ < 2513 Å (3p 3 P 0 ).The detailed bound-free edges are due to He II λ < 227 Å and He I λ < 504 Å in addition to the typical H I edges.
He II, He III, Ca II, Ca III, Mg II, Mg III, Fe II, Fe III, Si II, Si III, C II, Al II free-bound (f-b) H I, He I, He II, Mg I, Mg II, Si I, Si II, C I, C II, Al I bound-bound (b-b) H I, He I, He II, Ca II,Mg II, Fe II, Si II, C I, C II, Al II, O I, S I Note-The f-b species listed refer to the respective bound-states.The f-b transitions in Ca I, Ca II, Fe I, Fe II, and Al II are not in the CHIANTI v8.0.1 database and are thus also excluded.The b-b transitions in Ca I, Mg I, Fe I, Si I, and Cr II are not in the CHIANTI v8.0.1 database and are thus also excluded.We note that in CHIANTI v10.0.1 (the most recent version at the time of writing), b-b transitions in Cr II contribute a significant amount to the radiative losses at T = 10 4 K.

Figure 4 .
Figure 4. (a) Comparison of the optically thin loss function that is used in this work, the one used in Kowalski et al. (2015), and one of the standard versions from the RADYN code.(b) Calculations of the mF13-500-3 model gas temperatures using the standard thin loss function compared to this work.The temperature snapshots are obtained at t = 1 s and correspond to the height ranges in the bottom panels.The preflare temperature is shown as the dotted line.Panels (c) and (d) show net contributions to the changes of internal energy per gram.Individual flux divergence terms (divided by mass density) in the internal energy equation are labeled.The thin loss function in the grid is used in the snapshot in (d), and a standard optically thin loss function is used in (c).Note the large differences between the detailed radiative cooling (which is dominated by wavelengths covering the Balmer continuum range) and the thin radiative cooling in the low atmosphere.The model snapshots in (c) and (d) are the same as in panel (b).

Figure 5 .
Figure 5.The predictions from the main group of models with large values of Ec and a range of nonthermal energy fluxes are compared to observational constraints.(a) Balmer Hγ line-to-continuum ratios over the time corresponding to the optical continuum peak for a sample of flares with spectra.(b) Color-color diagram for these flares.(c) Color-color diagram for the peak phases of the flares with ULTRACAM data.The ULTRACAM peak-phase colors are annotated for the large events that are referred to as "impulsive flare 1" (IF1) and "impulsive flare 3" (IF3) inKowalski et al. (2016).The HST-1 and HST-2 data are added fromKowalski et al. (2019b), the F1/IF11 and F2/IF4 flares are added fromKowalski et al. (2013), and the rest of the spectral data in (a) and (b) are obtained fromKowalski et al. (2013) andHawley & Pettersen (1991).The open star symbols are calculated from "newly-formed" flare flux spectra during secondary events, and the large open circles are obtained from spectra during the impulsive phase of the Great Flare of AD Leo(Hawley & Pettersen 1991).Blackbody color temperatures are indicated in (b) and (c), and the parameterized chromospheric condensation (CC) models fromKowalski & Allred (2018) are included in (c).Error bars on the data indicate marginal 68% uncertainties.In panel (c), the systematic uncertainties are included in the gray error bars, while the black error bars are the statistical uncertainties only (note that the statistical uncertainties are very small for the IF1 and IF3 events).In (b) and (c), flares with only Balmer jump calculations are shown in the right panels in the margins.

Figure 6 .
Figure 6.Model spectra, F λ , of the Hγ emission lines during a flare on YZ CMi.(a) Observations of the Balmer Hγ line in three consecutive spectra.(b) Observations of Ca II K, Ca II H, and Balmer Hϵ.(c) The mF13-85-3 model spectrum (Kowalski 2022) well-reproduces the Hγ wing shape and the relative flux of the nearby continuum radiation in the continuum peak phase observation (dark blue line).(d) An additional narrow Gaussian well-reproduces the entire Hγ line profile.(e) The tF12-37-3 model spectrum (Namekata et al. 2020) well-reproduces the Hγ wing shape and the relative flux of the nearby continuum radiation in the continuum decay phase observations (light blue line).(f ) An additional narrow Gaussian well-reproduces the entire Hγ line profile.In panels (d) and (f), the pre-flare observation has been subtracted from the flare data.The RHD model surface fluxes have been scaled to the count fluxes (in units of counts s −1 ) of the nearby continuum levels in the spectral observations.The model spectra with nλ = 327 (Section 2.3) are used in this figure.

Figure 7 .
Figure 7. Contribution functions to the emergent µ = 0.95 intensity for the hydrogen Balmer (Ba) γ line at t = 0 s (a), at t = 2 s during the tF12-37-3 model (b), and at t = 1 s during the mF13-85-3 model (c).Panel (d) shows the contribution function over the hydrogen Paschen (Pa) β line at t = 1 s during the mF13-85-3 model.The gray scale is logarithmic from 10 −2 (white) to 10 1.5 (black) erg cm −2 s −1 sr −1 Å−1 cm −1 in all panels.The thin dashed (blue) lines indicate 5% (upper) and 95% (lower) of the cumulative of the contribution function, and the top axes indicate the line-of-sight (L.o.S.) velocities that correspond to the wavelength detuning and to the gas velocity (yellow line).Margin plots show the source-function equivalent temperatures (Sν T ) and the gas temperatures (bottom axes) compared to the ambient (thermal) electron densities (top axes).Note the different scales among the electron density axes.As noted inKowalski (2022), the main mF13 models produce non-negligible chromospheric upflows, but these are not apparent on the velocity scale shown here.
can be modeled in detail.It is also possible to readily synthesize optically thick UV resonance lines of C II, Mg II, and

Table 1 .
Species excluded from the CHIANTI optically thin radiative loss function