# Time-dependent Stellar Flare Models of Deep Atmospheric Heating

Adam F. Kowalski National Solar Observatory, University of Colorado Boulder,
3665 Discovery Drive, Boulder, CO 80303, USA  Department of Astrophysical and
Planetary Sciences, University of Colorado, Boulder, 2000 Colorado Ave, CO
80305, USA  Laboratory for Atmospheric and Space Physics, University of
Colorado Boulder, 3665 Discovery Drive, Boulder, CO 80303, USA.  Joel C.
Allred NASA Goddard Space Flight Center, Solar Physics Laboratory, Code 671,
Greenbelt, MD 20771, USA  Mats Carlsson Institute of Theoretical Astrophysics,
University of Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway  Rosseland
Centre for Solar Physics, University of Oslo, P.O. Box 1029 Blindern, 0315
Oslo, Norway

###### Abstract

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\approx 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
$\gamma$ lines in a powerful flare with echelle spectra. The self-consistent
formation of the wings and nearby continuum level provide 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.

##  1 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 $\approx$ 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\approx 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 optical (e.g., Kowalski
et al., 2019b). The optical color temperature attains the hottest values in
the rise and peak phases (Mochnacki & Zirin, 1980) when the Balmer jumps are
smallest (Kowalski et al., 2013). In the fast decay phase, the optical color
temperatures rapidly decrease and settle to a cooler value around $5500$ K
through the gradual decay phase (see also Kowalski et al., 2018; Kowalski,
2023). The evolution of blackbody temperatures in stellar flares has been
investigated with broadband photometry (Hawley et al., 2003; Robinson et al.,
2005; Zhilyaev et al., 2007; Howard et al., 2020), but the interpretation of
the continuum radiation can be difficult without spectra (Allred et al., 2006;
Kowalski et al., 2019b).

Large heating rates in gas-dynamic or radiative-hydrodynamic (RHD) models
produce moderate-to-large optical depths at continuum wavelengths and have
been successful in explaining some of the observed continuum properties
(Livshits et al., 1981; Katsova et al., 1981, 1991, 1997; Kowalski et al.,
2015; Kowalski & Allred, 2018). However, the beam current densities may be too
large to last during the propagation from the corona to the chromosphere.
Further, electron-beam-generated chromospheric condensations (e.g., Livshits
et al., 1981; Fisher et al., 1985; Fisher, 1989; Gan et al., 1992) in some of
these models produce hydrogen Balmer line spectra that are too broad
(Kowalski, 2022); other recent perspectives and arguments have been provided
by Belova & Bychkov (2019) and Morchenko (2020a, b). Stellar flare RHD models
have considered smaller electron beam heating fluxes, but the model Balmer
jumps are larger than in all known spectral observations; the photosphere is
radiatively backheated by this Balmer continuum radiation, which originates at
chromospheric heights (Allred et al., 2006). Other models of stellar
chromospheric flares include static NLTE slab calculations (Jevremovic et al.,
1998; García-Alvarez et al., 2002; Morchenko et al., 2015), static LTE slab
calculations (Kunkel, 1970; Eason et al., 1992; Heinzel & Shibata, 2018;
Simões et al., 2024), semi-empirical atmospheric adjustments with NLTE
radiative transfer (Houdebine et al., 1991; Christian et al., 2003;
Fuhrmeister et al., 2005, 2010; Kowalski et al., 2012; Schmidt et al., 2012),
and static energy equilibrium loop models with NLTE radiative transfer (Hawley
& Fisher, 1992). The semi-empirical model atmospheres of Cram & Woods (1982)
explored six variations of temperature and electron density adjustments within
the deep chromosphere and photosphere, one of which produced very broad Balmer
line wings and a hot optical blackbody color temperature. For a detailed
review of stellar flares and stellar flare modeling, see Kowalski (2024).

In this study, we present a grid of radiative-hydrodynamic stellar flare
models that reproduce the range of observed optical and near-ultraviolet
continuum and hydrogen Balmer emission line properties. The atmospheres are
driven by heating rates from power-law distributions of electrons. On the Sun,
nonthermal electrons can produce hard X-ray sources that are co-temporal and
co-spatial with the white-light radiation at the footpoints of newly-
reconnected magnetic loops (Hudson et al., 1992; Neidig et al., 1993, 1994; Xu
et al., 2006; Fletcher et al., 2007; Kawate et al., 2016). However, the
parameters of the power-law distributions in our grid are varied through a
larger range than considered for solar flare models (e.g., Carlsson et al.,
2023) so that deep heating achieves optical blackbody color temperatures of
$T>10^{4}$ K in the emergent radiative flux spectra. The grid thus essentially
follows on the static, semi-empirical model predictions of Cram & Woods (1982)
by including the evolution of self-consistent atmospheric opacities, flows,
Coulomb heating rates, and non-equilibrium radiative transfer. The detailed
spectra and the hydrodynamic output are made publicly available to facilitate
comparisons to and interpretations of stellar flare data. The grid spectra
have already been used for a wide range of purposes (Kowalski et al., 2017b;
Namekata et al., 2020; Kowalski, 2022, 2023; Brasseur et al., 2023; Howard et
al., 2023; Monson et al., 2024).

The paper is organized as follows. Section 2 contains an extensive description
of the model setup, the grid of flare heating parameters, and the context of
1D flare loop models. The analysis of the grid models (Section 3) is divided
into three parts. Section 3.1 presents an overview of the time-averaged
optical and NUV spectral quantities calculated from each model. In Section
3.2, the quantities from a subset of the grid are compared to global trends
among impulsive-phase observations of dMe flares. Section 3.3 demonstrates
detailed modeling of the hydrogen Balmer $\gamma$ line in a remarkable dMe
event whose broad lines have not yet been discussed in detail elsewhere. We
discuss the applicability and limitations of the model grid in Section 4, and
Section 5 includes some concluding remarks. The Appendix describes how to
access the model output and additional documentation.

##  2 Radiative-Hydrodynamic Flare Models

###  2.1 Model Setup

We calculate a grid of RHD models with the RADYN code (Carlsson & Stein, 1992,
1995, 1997, 2002; Allred et al., 2015; Carlsson et al., 2023). To simulate
flare heating, we model the energy deposition from a power-law distribution
(hereafter, “beam”) of electrons. The one-dimensional equations of mass,
momentum, internal energy, and charge 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 semi-implicit (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
$\mathcal{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 log10
$g_{\star}/[$cm s${}^{-2}]=4.75$, and the calculated effective temperature is
$T_{\rm{eff}}\approx 3600$ K. The loop apex has an electron density of
$n_{e}\approx 3\times 10^{10}$ cm-3 and plasma temperature of $T\approx 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 corona111We apply a constant
energy deposition of $0.23$ erg cm-3 s-1 at $z\geq 460$ km, which causes
thermal conduction to transition from a net cooling term at $z>1.20$ Mm to a
net heating term at $0.41<z<1.20$ Mm. 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 $\tau(\lambda=500\ \mathrm{nm})=\tau_{500}=1$ in
the preflare state, where $\tau$ is the optical depth. The one-dimensional
formulation of the equations of hydrodynamics is justified if $\beta\ll 1$,
where $\beta$ 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 ($\Delta t\leq 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 $\delta=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., 2020, 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}\approx 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}\approx 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
$\delta=3$ electron-beam power-law distributions that contain a large
percentage, up to $\approx 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}\approx 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 $\approx 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 $\Delta 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.

###  2.2 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, $\Delta t\approx
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, $\Delta t\approx 0.1-3$ s, bursts are
observed in H$\alpha$ 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., 2007, 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) assumes adiabatic222Meaning that the pitch angles change as the
particles bounce around in this region, while the total energy of each
particle remains constant according to the adiabatic invariants. The so-called
betatron heating and acceleration (e.g., Veronig et al., 2006; Birn et al.,
2017) are not considered. 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_{\rm{FWHM}}=2.3$ s, the height of the fully contracted looptops,
$h_{L}=0.64\times 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

$E_{\rm{c}}\approx\frac{0.15}{2}m_{i}c_{A,X}^{2}\approx
0.15\frac{B_{X}^{2}}{8\pi n_{e}},$ (1)

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}\approx 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

$F_{fl}\approx\frac{1}{2\pi}B_{X}^{2}c_{A,X}\sin^{4}(\Delta\theta/4)\approx
0.9-5\times 10^{13}\ \mathrm{erg}\ \mathrm{cm}^{-2}\ \mathrm{s}^{-1}$ (2)

This range of flare fluxes corresponds to a range of shear angles,
$\Delta\theta$, 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
$\approx$ kG fields in the low corona.

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 $\approx 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
$g_{z}=g_{\star}\sin(\theta^{\prime})$. 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).

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, 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\approx 8-12$ kG, may cover $\approx
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”333In a different off-limb flare, Liu et al. 2013
discovered “fast downward contractions” and “slow downward shrinkages”, which
are most easily seen after the hard X-ray impulsive phase. 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., 2023). The 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_{\rm{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 $\delta=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 $\Delta t=2.3$ s). In Figure 2, we also show the
evolution of the detailed emergent radiative flux at a representative
wavelength, $\lambda=3615$ Å, that is just shortward of the ideal Balmer
recombination limit ($\lambda=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).

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 $\Delta t=0.2$ s; (_Bottom_) a pulsed step-function injection in the
cF13-85-3 model at $\Delta t=0.05$ s. The evolution of the emergent radiative
flux at a Balmer continuum wavelength is shown in each panel for comparison.

###  2.3 Emergent Radiative Flux Spectra

RADYN provides spectra from $\lambda=6-40,000$ Å in physical (cgs) units,
which facilitates comparisons to flux-calibrated stellar observations.
Emergent radiative surface flux density spectra (erg cm-2 s-1 $\rm{\AA}^{-1}$;
hereafter, just “flux”) are calculated within RADYN using a Gaussian
quadrature sum of the emergent intensity at five $\mu$ values
($\mu=0.05,0.23,0.50,0.77,0.95$, where $\mu=\cos(\theta)$ and $\theta$ is the
angle with respect to the atmospheric normal), and at 95 continuum
wavelengths. We denote the flux at a specific continuum wavelength $\lambda$
as “C$\lambda(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] $\approx-2$ and $-3$ to gas
temperatures of $T\gtrsim 10^{4}$ 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 extreme-ultraviolet (XEUV) spectral
regions are similar in shapes. The color temperatures that are inferred from
the X-rays at $\lambda<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
(TOPBASE444http://cdsarc.u-strasbg.fr). 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 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 H2 (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_{\rm{eff}}\approx 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.

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 $\lambda<1100$ Å ($2p^{2}\ {}^{3}P$), S I
$\lambda<1199$ Å, Si I $\lambda<1525$ Å ($3p^{2}\ {}^{3}P$), Si I
$\lambda<1682$ Å ($3p^{2}\ {}^{1}D$), Al I $\lambda<2076$ Å ($3p^{2}\ P^{0}$),
and Mg I $\lambda<2513$ Å ($3p\ ^{3}P^{0}$). The detailed bound-free edges are
due to He II $\lambda=227$ Å and He I $\lambda<504$ Å in addition to the
typical H I edges.

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$\alpha$, H$\beta$, H$\gamma$) are modeled
using the Doppler-convolved, TB09+HM88 line profile functions (Vidal et al.,
1970, 1971, 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,
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$\alpha$, H$\beta$, and H$\gamma$ lines are calculated on frequency grids
with 51, 31, and 31 points, respectively. The sampling in the far wings
($|\lambda-\lambda_{0}|\gtrsim 2$ Å) is courser than around the rest
wavelength. To mitigate systematic errors in the line-integrated fluxes of the
H$\gamma$ 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}\phi_{\alpha}$ on a grid of $T$ and $n_{e}$, where $\alpha$ is the
typical detuning parameter $\propto|\lambda-\lambda_{0}|$, and $\phi_{\alpha}$
is the normalized line profile function. This is followed by a four-point,
third-order interpolation scheme on a grid of $\log_{10}\alpha$ and
$\log_{10}\phi_{\alpha}$ (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$\gamma$ and low-resolution spectra of the Balmer H$\beta$ and H$\alpha$
lines. The diagnostic advantages of Balmer H$\gamma$ 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$\alpha$ line than have been previously
demonstrated (Stehle, 1994). However, the new line shape calculations of
H$\gamma$ are essentially identical to the TB09+HM88 profiles that we use.
Extending the _Xenomorph_ calculations of H$\alpha$ to lower electron
densities, which are relevant for flare atmospheres, will be the subject of
future work (Kowalski & Gomez, in prep).

###  2.4 XUV Backheating

X-ray and ultraviolet (XUV; $\lambda=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
grid555Notably, however, the preflare emergent intensities of the Balmer lines
are factors of $2-3$ fainter when the starting atmosphere is relaxed with the
method of Allred et al. (2015)..

###  2.5 Optically Thin Radiative Loss Function

The optically thin radiative loss function, $\Lambda_{\rm{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\approx 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 $\approx 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\gtrsim 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\approx
10^{4}$ K, originate primarily666In decreasing order, the six largest
contributions to the bound-free radiative cooling at $T=10^{4}$ K are H I, Si
I, Mg I, S I, C I, and Al I. In decreasing order, the six largest
contributions to the bound-bound radiative cooling at $T=10^{4}$ K are H I, Mg
II, Fe II, Ca II, Si II, and Al II. Note, the bound-bound cooling is generally
much larger than the bound-free. The main contributions to the remaining
radiative cooling in our adopted thin loss function at $T=10^{4}$ K are bound-
bound transitions in N I, S II, and Ni II in decreasing order.  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\approx 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 bound-bound 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 $\lambda 2632$
emission line, the emergent wavelength-integrated emergent intensity is
$\approx 200$x greater than the optically thick calculation. For the faint Fe
II $\lambda 2814.445$ line, the intensity is a factor of $\approx 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.

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).  Table 1: Species excluded from the
CHIANTI optically thin radiative loss function

Type | Species
---|---
free-free (f-f) | H II, 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.

###  2.6 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, 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, $\Delta 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. (2017a, 2022). The auxiliary models include the following: a
recalculation of a constant injection of F13 flux with $E_{c}=37$ keV and
$\delta=3$ for 1.6 s (following Kowalski et al., 2016), a high-flux model with
a very hard power-law index, $\delta=2.5$ (m2F12-37-2.5), and a recalculation
of the soft ($\delta=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 $\Delta
t=2.3$ s, and another for $\Delta 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 $\approx 10^{12}$ and $5\times 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 $\Delta 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 $\approx 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.

##  3 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).

###  3.1 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_{\lambda}$, from each model777In Kowalski (2022), the model surface flux
spectra were abbreviated as $S_{\lambda}$; here, we do not want to confuse
this with the source function, $S_{\nu}$.. 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_{\lambda}$ 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 $\Delta t=10$ s,
beginning when the conductive pulses reach the chromosphere. The flare-only,
time-averaged radiative surface flux spectra are denoted as
$F_{\lambda}^{{}^{\prime}}$, and the flare-only radiative surface flux spectra
at time $t$ are indicated as $F_{\lambda}^{\prime}(t)$. Note that Figure 3
displays $F_{\lambda}(t=1\rm{s})$ 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 3.1 lists several of the
calculations, and others are described in the supplemental document that is
referenced in the Appendix. The columns of Table 3.1 are the following:

  1. 1.

Unique model identifier.

  2. 2.

C4170′: The value of $F_{\lambda}^{{}^{\prime}}$ at the continuum wavelength
of $\lambda=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).

  3. 3.

$F_{\rm{H}\gamma}^{\prime}$: continuum-subtracted $F_{\lambda}^{{}^{\prime}}$
integrated over the wavelengths, $\lambda_{1}$ to $\lambda_{2}$, of the Balmer
H$\gamma$ 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_{\lambda}(t=0)$ has been subtracted.

  4. 4.

$F_{\rm{H}\gamma}^{\prime}$ /C4170′: the ratio of column (3) to (2)

  5. 5.

$T_{\rm{BB}}(1\rm{s})$: blackbody temperature fit at the model continuum
wavelengths over $\lambda=4000-4800$ Å, following Kowalski et al. (2013); this
is calculated from $F_{\lambda}^{\prime}(t=1\rm{s})$.

  6. 6.

$T_{\rm{BB}}(2.2\rm{s})$: blackbody temperature fit at the model continuum
wavelengths over $\lambda=4000-4800$ Å, following Kowalski et al. (2013); this
is calculated from $F_{\lambda}^{\prime}(t=2.2\rm{s})$.

  7. 7.

$T_{\rm{BB}}$: blackbody temperature fit at the model continuum wavelengths
over $\lambda=4000-4800$ Å, following Kowalski et al. (2013); this is
calculated as $F_{\lambda}^{\prime}$ from the time-averaged surface flux
spectrum. The temporal calculations of $T_{\rm{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_{\rm{FcolorR}}$ (not listed in this table), is calculated from
the ratio of $F_{\lambda}^{\prime}$ at $\lambda=4170$ Å and $\lambda=6010$ Å.
Spectra indicate that $T_{\rm{BB}}\neq T_{\rm{FcolorR}}$ (Kowalski et al.,
2016). In Section 3.2, the model values of $T_{\rm{FcolorR}}$ are compared to
comprehensive observational constraints over a broader wavelength range than
fit with $T_{\rm{BB}}$.

  8. 8.

C3615′/C4170′: the Balmer jump ratio, which is the flare-only continuum flux
at $\lambda=3615$ Å divided by column 2; this is calculated from
$F_{\lambda}^{\prime}$.

  9. 9.

H$\gamma$ Eff. Width: the effective width (Kowalski et al., 2022) of the
H$\gamma$ line; this is calculated from $F_{\lambda}$ after a linear fit to
the nearby continuum flux on both sides of the emission line has been
interpolated and subtracted:

$\Delta\lambda_{\rm{eff}}=\int_{\lambda_{1}}^{\lambda_{2}}\frac{F_{\lambda}(\lambda)-F_{\lambda,\rm{cont}}(\lambda)}{F_{\lambda}(\lambda_{\rm{max}})-F_{\lambda,\rm{cont}}(\lambda_{\rm{max}})}d\lambda$
(3)

where $\lambda_{1}=4312$ Å and $\lambda_{2}=4368$ Å.

  10. 10.

$T_{\rm{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_{\rm{max}}$.

  11. 11.

$T_{\rm{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_{\lambda}(t=t_{\rm{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 $\lambda=8$ Å to that at $\lambda=35$ Å. An X-ray color temperature
($T_{\rm{X,max}}$; not listed) is calculated from the ratio of the emergent
flux at $\lambda=0.5$ Å to that at $\lambda=4$ Å.

The values in Table 3.1 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\approx 15$ MK to $T\gg 50$ MK (e.g., Güdel, 2004; Mullan et al.,
2006; Liefke et al., 2010; Osten et al., 2005, 2010, 2016). Second, a few of
the models produce maximum plasma temperatures that are extremely high,
$T_{\rm{plasma,max}}\approx 40-90$ MK. However, the thermal bremsstrahlung
temperatures that are calculated from XEUV continuum spectra over
$\lambda=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_{\rm{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_{\rm{X,max}}-T_{\rm{plasma,max}}\approx 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}\approx 1750$ km s-1 evaporation reaches the looptop at $t=6$
s, and then there is a jump in the apex temperature from $T\approx 50$ MK to
$\approx 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\approx 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\times 10^{11}$ cm-3 and very hot,
$T\approx 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\approx 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 ($\Delta t=5$ or $10$ s) to provide the average quantities in
Table 3.1. 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.

Model ID | C4170′ | $F_{\rm{H}\gamma}^{\prime}$ |  $F_{\rm{H}\gamma}^{\prime}$ /C4170′ | $T_{\rm{BB}}(1\rm{s})$ | $T_{\rm{BB}}(2.2\rm{s})$ | $T_{\rm{BB}}$ | C3615′/C4170′ | H$\gamma$ Eff. Width | $T_{\rm{plasma,max}}$ | $T_{\rm{XEUV,max}}$
---|---|---|---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
| erg cm-2 s-1 Å-1 | erg cm-2 s-1 | Å | K | K | K |  $\cdots$ | Å |  $10^{6}$ K |  $10^{6}$ K
$t=0$ s | 2.5e+05 | 6.24e+06 | 25 |  $\cdots$ |  $\cdots$ | 4490 | 0.7 | 0.49 | 5.0 | 4.2
main group |  |  |  |  |  |  |  |  |  |
mF10-85-3 | 1.54e+04 | 1.07e+07 | 691.2 | 6200 | 5400 | 5600 | 40.62 | 1.06 | 5.0 | 4.2
mF10-150-3 | 8.46e+03 | 6.85e+06 | 810.4 | 4500 | 6000 | 6600 | 56.57 | 0.89 | 5.0 | 4.2
mF10-200-3 | 1.28e+03 | 4.98e+06 | 3885.9 |  † $\cdots$ |  $\cdots$ |  $\cdots$ | 309.23 | 0.78 | 5.0 | 4.2
mF10-350-3 | -1.86e+04 | 2.31e+06 | -124.6 |  † $\cdots$ |  $\cdots$ |  $\cdots$ | -15.25 | 0.62 | 5.0 | 4.2
mF10-500-3 | -3.03e+04 | 1.48e+06 | -48.7 |  † $\cdots$ |  $\cdots$ |  $\cdots$ | -8.17 | 0.57 | 5.0 | 4.2
mF11-85-3 | 4.78e+05 | 8.86e+07 | 185.2 | 5200 | 5200 | 5100 | 12.02 | 3.84 | 5.0 | 4.3
mF11-150-3 | 5.39e+05 | 6.93e+07 | 128.7 | 5300 | 5300 | 5200 | 9.32 | 4.27 | 5.0 | 4.2
mF11-200-3 | 5.61e+05 | 6.00e+07 | 107.0 | 5400 | 5400 | 5200 | 7.99 | 4.11 | 5.0 | 4.2
mF11-350-3 | 1.20e+06 | 9.09e+07 | 75.5 | 5300 | 5400 | 5400 | 6.07 | 5.66 | 5.0 | 4.2
mF11-500-3 | 5.25e+05 | 2.98e+07 | 56.8 | 5000 | 5400 | 5200 | 4.84 | 2.53 | 5.0 | 4.2
mF12-85-3 | 8.33e+06 | 7.53e+08 | 90.4 | 6100 | 6200 | 5900 | 4.97 | 6.55 | 5.8 | 5.1
mF12-150-3 | 1.37e+07 | 6.58e+08 | 48.1 | 6400 | 6600 | 6300 | 3.45 | 11.08 | 5.0 | 4.3
mF12-200-3 | 1.61e+07 | 4.85e+08 | 30.1 | 6600 | 7000 | 6500 | 2.63 | 14.41 | 5.0 | 4.3
mF12-350-3 | 1.89e+07 | 2.35e+08 | 12.4 | 7000 | 7900 | 7100 | 1.71 | 23.97 | 5.0 | 4.2
mF12-500-3 | 2.00e+07 | 1.31e+08 | 6.5 | 7200 | 8500 | 7400 | 1.41 | 26.05 | 5.0 | 4.2
m2F12-85-3 | 2.00e+07 | 1.15e+09 | 57.3 | 6700 | 7100 | 6600 | 3.38 | 8.05 | 7.2 | 5.9
m2F12-150-3 | 3.72e+07 | 9.46e+08 | 25.5 | 7600 | 8800 | 7700 | 2.28 | 11.76 | 5.1 | 4.4
m2F12-200-3 | 4.55e+07 | 5.87e+08 | 12.9 | 8300 | 10000 | 8500 | 1.79 | 12.16 | 5.0 | 4.2
m2F12-350-3 | 7.71e+07 | -2.60e+08 | -3.4 | 10000 | 12500 | 10400 | 1.06 | 6.62 | 5.0 | 4.2
m2F12-500-3 | 5.80e+07 | -4.10e+08 | -7.1 | 10600 | 13600 | 11100 | 0.86 | 14.95 | 5.0 | 4.2
mF13-85-3 | 1.14e+08 | 1.62e+09 | 14.2 | 11100 | 12300 | 10800 | 1.73 | 9.13 | 16.1 | 10.1
mF13-150-3 | 1.54e+08 | 1.16e+09 | 7.5 | 12900 | 13400 | 12100 | 1.54 | 10.33 | 7.0 | 5.7
mF13-200-3 | 1.73e+08 | 6.21e+08 | 3.6 | 13800 | 14200 | 12800 | 1.42 | 10.38 | 5.3 | 4.7
mF13-350-3 | 1.93e+08 | -7.08e+08 | -3.7 | 16200 | 16500 | 14900 | 1.12 | 13.78 | 5.0 | 4.2
mF13-500-3 | 1.87e+08 | -1.37e+09 | -7.3 | 18200 | 18900 | 16500 | 0.95 | 14.96 | 5.0 | 4.2
const group |  |  |  |  |  |  |  |  |  |
cF11-85-3 | 8.40e+05 | 1.48e+08 | 176.4 | 5200 | 5300 | 5200 | 11.56 | 4.92 | 5.0 | 4.3
cF11-150-3 | 9.35e+05 | 1.15e+08 | 123.4 | 5300 | 5400 | 5300 | 9.09 | 5.97 | 5.0 | 4.2
cF11-200-3 | 9.68e+05 | 1.00e+08 | 103.6 | 5400 | 5400 | 5300 | 7.84 | 6.08 | 5.0 | 4.2
cF11-350-3 | 9.69e+05 | 6.94e+07 | 71.5 | 5400 | 5500 | 5400 | 5.75 | 5.23 | 5.0 | 4.2
cF11-500-3 | 8.91e+05 | 4.95e+07 | 55.6 | 5400 | 5500 | 5400 | 4.73 | 4.00 | 5.0 | 4.2
c5F11-85-3 | 6.32e+06 | 7.15e+08 | 113.2 | 5900 | 5900 | 5800 | 6.11 | 6.14 | 5.3 | 4.6
c5F11-150-3 | 9.16e+06 | 5.76e+08 | 62.8 | 5900 | 6100 | 5800 | 4.55 | 10.80 | 5.0 | 4.2
c5F11-200-3 | 9.98e+06 | 4.64e+08 | 46.5 | 6000 | 6200 | 5900 | 3.72 | 17.07 | 5.0 | 4.2
c5F11-350-3 | 1.10e+07 | 3.35e+08 | 30.6 | 6200 | 6500 | 6100 | 2.69 | 20.95 | 5.0 | 4.2
c5F11-500-3 | 1.13e+07 | 2.63e+08 | 23.4 | 6300 | 6700 | 6300 | 2.27 | 21.34 | 5.0 | 4.2
cF12-85-3 | 1.45e+07 | 1.14e+09 | 79.1 | 6200 | 6400 | 6200 | 4.30 | 7.86 | 5.8 | 4.9
cF12-150-3 | 2.61e+07 | 1.12e+09 | 42.8 | 6600 | 7200 | 6600 | 2.99 | 11.80 | 5.0 | 4.3
cF12-200-3 | 3.13e+07 | 7.64e+08 | 24.4 | 6900 | 7800 | 7000 | 2.26 | 14.74 | 5.0 | 4.3
cF12-350-3 | 3.75e+07 | 2.22e+08 | 5.9 | 7600 | 9200 | 7800 | 1.39 | 37.62 | 5.0 | 4.2
cF12-500-3 | 3.90e+07 | 4.21e+07 | 1.1 | 8000 | 10000 | 8200 | 1.16 | 23.82 | 5.0 | 4.2
c2F12-85-3 | 3.50e+07 | 1.58e+09 | 45.1 | 7000 | 7700 | 7000 | 2.81 | 10.13 | 7.2 | 5.8
c2F12-150-3 | 6.54e+07 | 1.39e+09 | 21.3 | 8300 | 9700 | 8400 | 2.03 | 12.41 | 5.1 | 4.4
c2F12-200-3 | 8.05e+07 | 8.87e+08 | 11.0 | 9200 | 10800 | 9300 | 1.68 | 13.02 | 5.0 | 4.2
c2F12-500-3 | 1.01e+08 | -8.44e+08 | -8.3 | 12300 | 15100 | 12000 | 0.81 | 14.14 | 5.0 | 4.2
c5F12-85-3 | 9.78e+07 | 1.65e+09 | 16.9 | 9600 | 11000 | 9600 | 1.81 | 11.42 | 11.5 | 7.6
c5F12-150-3 | 1.46e+08 | 1.20e+09 | 8.2 | 11500 | 12600 | 11300 | 1.55 | 12.17 | 5.6 | 4.8
c5F12-200-3 | 1.69e+08 | 7.65e+08 | 4.5 | 12500 | 13300 | 12100 | 1.44 | 12.18 | 5.1 | 4.4
cF13-85-3 | 1.65e+08 | 1.44e+09 | 8.7 | 11900 | 12800 | 11600 | 1.54 | 11.48 | 17.3 | 8.9
cF13-150-3 | 2.10e+08 | 8.73e+08 | 4.2 | 13500 | 14100 | 13200 | 1.42 | 10.38 | 6.9 | 5.6
cF13-200-3 | 2.27e+08 | 5.07e+08 | 2.2 | 14300 | 14900 | 14000 | 1.37 | 9.38 | 5.3 | 4.6
cF13-350-3 | 2.43e+08 | -5.27e+08 | -2.2 | 16500 | 17000 | 15800 | 1.23 | 9.02 | 5.0 | 4.2
cF13-500-3 | 2.64e+08 | -1.19e+09 | -4.5 | 18600 | 19100 | 17500 | 1.14 | 10.12 | 5.0 | 4.2
Ec37 group |  |  |  |  |  |  |  |  |  |
mF10-37-3 | 1.69e+04 | 1.64e+07 | 967.6 | 5100 | 5100 | 5200 | 41.57 | 1.15 | 5.1 | 4.2
m5F10-37-3 | 1.49e+05 | 7.50e+07 | 502.3 | 5600 | 5400 | 5300 | 17.30 | 1.71 | 5.4 | 4.5
mF11-37-3 | 3.21e+05 | 1.24e+08 | 385.6 | 5700 | 5500 | 5400 | 13.33 | 1.95 | 5.8 | 4.9
m2F11-37-3 | 6.55e+05 | 1.93e+08 | 294.2 | 5700 | 5500 | 5500 | 10.88 | 2.31 | 6.7 | 5.6
m5F11-37-3 | 1.67e+06 | 3.41e+08 | 203.8 | 5800 | 5700 | 5600 | 8.47 | 3.05 | 9.4 | 7.2
mF12-37-3 | 3.51e+06 | 5.45e+08 | 155.2 | 5900 | 5800 | 5700 | 6.87 | 3.85 | 13.6 | 8.9
m2F12-37-3 | 7.93e+06 | 9.33e+08 | 117.8 | 6100 | 6100 | 5900 | 5.33 | 4.99 | 21.1 | 10.3
m5F12-37-3 | 2.92e+07 | 1.99e+09 | 68.1 | 6700 | 7400 | 6800 | 3.22 | 8.41 | 39.3 | 15.9
m7.5F12-37-3 | 5.34e+07 | 2.37e+09 | 44.4 | 7400 | 8700 | 7900 | 2.44 | 10.59 | 56.1 | 18.1
mF13-37-3 | 7.94e+07 | 2.22e+09 | 27.9 | 8200 | 10500 | 9500 | 1.96 | 12.11 | 88.2 | 17.8
m5F10-37-5 | 1.37e+05 | 8.54e+07 | 622.9 | 6000 | 5700 | 5600 | 14.03 | 1.47 | 5.5 | 4.7
mF12-37-5 | 1.10e+06 | 4.14e+08 | 375.9 | 6200 | 6100 | 6000 | 9.31 | 2.50 | 17.0 | 9.3
sol group |  |  |  |  |  |  |  |  |  |
mF10-17-3 | 1.60e+04 | 2.47e+07 | 1542.0 | 5400 | 5300 | 5300 | 39.73 | 0.97 | 5.4 | 4.5
mF11-17-3 | 1.58e+05 | 9.39e+07 | 595.7 | 5800 | 5500 | 5500 | 16.27 | 1.51 | 9.7 | 7.3
m2F11-17-3 | 3.37e+05 | 1.56e+08 | 461.9 | 5800 | 5500 | 5500 | 13.15 | 1.84 | 14.1 | 9.1
m5F10-17-5 | 3.95e+04 | 4.67e+07 | 1182.4 | 5900 | 5600 | 5400 | 19.00 | 1.04 | 9.1 | 6.9
m5F11-25-4 | 6.35e+05 | 2.95e+08 | 464.6 | 6000 | 5900 | 5800 | 10.80 | 2.20 | 18.6 | 10.0
c15s-5F11-25-4 | 3.30e+06 | 1.33e+09 | 402.4 | 6000 | 6000 | 5700 | 9.12 | 5.33 | 60.3 | 17.9
thermal group |  |  |  |  |  |  |  |  |  |
thermal-long | 3.17e+04 | 6.62e+07 | 2086.2 | 4400 | 4400 | 5300 | 18.94 | 1.28 | 31.9 | 10.4
thermal-short | 1.36e+03 | 1.31e+07 | 9652.2 | 5700 | 5700 | 4600 | 175.70 | 0.70 | 21.4 | 9.0
auxiliary group |  |  |  |  |  |  |  |  |  |
m2F12-37-2.5 | 1.72e+07 | 8.65e+08 | 50.2 | 6500 | 7200 | 6500 | 2.99 | 5.80 | 15.8 | 9.5
c1.6s-F13-37-3 | 7.40e+07 | 2.31e+09 | 31.3 | 9300 | 8300 | 9000 | 2.13 | 12.02 | 71.4 | 17.9
c2F12-500-7 | 9.76e+07 | -8.25e+08 | -8.5 | 13600 | 14800 | 12200 | 0.91 | 13.23 | 5.0 | 4.2
N+20 group |  |  |  |  |  |  |  |  |  |
tF12-37-3 | 2.41e+06 | 3.64e+08 | 151.0 | 5800 | 5900 | 5700 | 6.85 | 4.02 | 12.3 | 8.0
tF12-37-5 | 6.50e+05 | 2.32e+08 | 357.3 | 6400 | 6100 | 6000 | 9.55 | 2.46 | 16.5 | 8.8
tF10-37-5 | 3.07e+04 | 2.19e+07 | 713.6 | 4900 | 5000 | 4900 | 24.22 | 1.25 | 5.0 | 4.2
thermal-t5F10 | -5.60e+03 | 6.76e+06 | -1209.0 | 5200 | 5200 | 5400 | -35.76 | 0.64 | 16.0 | 7.9
thermal-tF12 | 1.48e+05 | 1.52e+08 | 1028.5 | 4400 | 5100 | 5500 | 14.53 | 2.01 | 89.9 | 9.7

††footnotetext: The mF10-200-3 model has very faint optical continuum flux
(column 2), and thus reliable blackbody fits (columns 5, 6, and 7) can not be
determined. 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.

###  3.2 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 $\lambda\approx 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 3.1) from the main group
of models to the impulsive phase observations from these studies.

The Balmer jump ratio correlates with the H$\gamma$ 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\approx 10,000$ K. In Figure 5(b), the Balmer jump ratios are shown
against the optical color temperatures, $T_{\rm{FcolorR}}$, that are
calculated from the ratio of the flare-only continuum fluxes at $\lambda=4170$
Å and $\lambda=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 spectrum888The blackbody fit to the continuum
regions at $\lambda\approx 3915-4420$ Å (Kowalski et al., 2013) is
extrapolated to $\lambda=6010$ Å. 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.

In summary, the high-flux ($10^{12}-10^{13}$ erg cm-2 s-1), large $E_{c}$ main
models generate heating rates that reproduce the most extreme dMe flare
continuum spectral properties. Models with fluxes between $10^{10}-10^{12}$
erg cm-2 s-1 and smaller low-energy cutoffs between $17-37$ keV (e.g., the
mF12-37-3 and m5F11-25-4 models; Table 3.1) produce Balmer jumps, optical
color temperatures, and Balmer line-to-continuum ratios that are inconsistent
with the cluster of flux ratios from impulsive-phase spectral observations.

Figure 5:  The predictions from the main group of models with large values of
$E_{c}$ and a range of nonthermal energy fluxes are compared to observational
constraints. (a) Balmer H$\gamma$ 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) 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.

###  3.3 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
$\lambda\approx 4300$ Å is $0.027$ Å pix-1. A series of four consecutive
Stokes I spectra covering $\lambda=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 $\Delta\lambda=-0.35$ Å,
which was required to align the pre-flare Ca II K emission line centroids in
the model and observations.

Spectra around Balmer H$\gamma$ and H$\epsilon$ 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_{\rm{H}\gamma}^{\prime}$ /C4170′ increases from $\approx 10$ to $\approx
100$, and the effective width (Eq. 3) of the H$\gamma$ line decreases from
$\Delta\lambda_{\rm{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).

Figure 6:  Model spectra, $F_{\lambda}$, of the H$\gamma$ emission lines
during a flare on YZ CMi. (a) Observations of the Balmer H$\gamma$ line in
three consecutive spectra. (b) Observations of Ca II K, Ca II H, and Balmer
H$\epsilon$. (c) The mF13-85-3 model spectrum (Kowalski, 2022) well-reproduces
the H$\gamma$ 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$\gamma$ line profile.
(e) The tF12-37-3 model spectrum (Namekata et al., 2020) well-reproduces the
H$\gamma$ 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$\gamma$ 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\lambda=327$ (Section 2.3) are used in
this figure.

In Figure 6(b), the H$\epsilon$ 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$\gamma$ line and the nearby continuum flux
(Figure 6(c)). The model spectrum, $F_{\lambda}$, at $n\lambda=327$ wavelength
points (Section 2.3) is multiplied by a scale factor so that it matches the
approximate continuum count fluxes around $\lambda\approx 4320$ and
$\lambda\approx 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$\beta$
and H$\delta$ model predictions are also within the constraints of the
observed wings of these lines (not shown). The H$\gamma$ 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$\gamma$. 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 $\approx 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$\gamma$
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.

####  3.3.1 Line Formation Properties

The mF13-85-3 and tF12-37-3 models satisfactorily explain the H$\gamma$
broadening at $|\lambda-\lambda_{0}|\gtrsim 0.5$ Å in the continuum peak and
continuum decay phase flare spectra, respectively. In this section, we compare
the formation of the H$\gamma$ 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$\alpha$
lines in a superflare on AD Leo.

Figure 7 shows contribution functions to the emergent radiative intensity
across the hydrogen Balmer $\gamma$ line. The preflare formation of the model
Balmer $\gamma$ 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 $\tau(\lambda)=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 $\Delta
z\approx 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}\approx 1-2.5\times 10^{15}$ cm-3 at heights well above the heights
that correspond to $\tau(\lambda)=1$ in the preflare and in the tF12-37-3
models. The formation of the Balmer $\gamma$ line wings in the mF13-85-3 model
is very close to LTE.

The formation of the hydrogen Paschen $\beta$ ($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 $\gamma$ and thus should experience
approximately the same amount of Stark broadening. This transition is also
more optically thin than the Balmer $\gamma$ 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 $\tau=1$ occurs
at selected wavelengths in the top panel. The continuum wavelengths
($\lambda\approx 12150$ Å) near Paschen $\beta$ are more optically thick at a
given height in the atmosphere than the blue continuum wavelengths
($\lambda\approx 4300$ Å) around Balmer $\gamma$. This pushes the formation of
the Paschen $\beta$ 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
$\beta$ line is closer to the gas temperature, indicating that it is formed
closer to LTE than Balmer $\gamma$, 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 $\gamma$. Here, we use the
spectra calculated at $n\lambda=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\times 10^{15}$ cm-3 and
$n_{e}=5\times 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\times 10^{15}$
cm-3 (Figure 7(c)). This discrepancy is due to the large optical depths
($\tau\approx 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\times 10^{15}$ cm-3 slab calculation, but the actual electron
densities where the line is formed are an order of magnitude smaller999All of
the discrepancies that were discussed remain if the RADYN model spectra are
scaled to 0.75 instead of 1.0. (Figure 7(b)).

It is also valuable to compare the time-averaged flux profile ($F_{\lambda}$)
that fits the Balmer wings in the YZ CMi flare (Figure 6) to the intensity and
flux spectra ($F_{\lambda}(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 core101010This is
expected due to the plane-parallel calculation of the radiative transfer. At
smaller $\mu$ values, the optical depth accrues to larger values at a given
height; the formation of the core then occurs at smaller values of the source
function..

Figure 7:  Contribution functions to the emergent $\mu=0.95$ intensity for the
hydrogen Balmer (Ba) $\gamma$ 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) $\beta$ 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_{\nu}$ $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.  Figure 8:  (top) Emergent
continuum intensity spectrum ($\mu=0.95$) for the mF13-85-3 model at $t=1$ s.
(bottom) Column mass corresponding to $\tau(\lambda)=1$ for a representative
subset of 53 continuum wavelengths from the top panel between $\lambda=700$
and $\lambda=20900$ Å. In the bottom panel, the central wavelengths for
hydrogen Balmer $\gamma$ and Paschen $\beta$ are indicated (see text).  Figure
9:  Model spectra around the hydrogen Balmer $\gamma$ 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 $n_{e}=5\times
10^{13},10^{14},5\times 10^{14},10^{15},2\times 10^{15},$ and $5\times
10^{15}$ cm-3 at $T=10^{4}$ K.

##  4 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\gtrsim 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$\gamma$ line is remarkably
well-reproduced.

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., 2023). Another possibility is that
time-dependent transport effects (Kontar et al., 2012) generate an enhancement
of $E\gtrsim 100$ keV electrons in stellar flares. The enhancement produces a
similar local heating maximum at $m_{c}\approx 0.01$ g cm-2 as in the
$E_{c}=85$ keV models with $\delta=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 3.1, which we refer
to as the v1.0 models, simulate flare heating ($Q_{\rm{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 $\lambda$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) can be modeled in
detail. It is also possible to readily synthesize optically thick UV resonance
lines of C II, Mg II, and possibly Si IV (Kerr et al., 2019b) with a
complementary radiative transfer code (e.g., Uitenbroek, 2001; Kerr et al.,
2019a; Zhu et al., 2019).

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 $\gamma$ line (Section
3.3). The impulsive-phase and gradual-phase blue-optical continuum levels
together with the H$\gamma$ 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 $\approx
170$ keV to $\approx 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., 2010, 2016). In our grid, the short durations
($\Delta t\approx 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., 2010, 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\times
10^{10}$ cm-3, the stopping distance (Emslie, 1978) for $E_{0}=75$ keV beam
electron is 170 Mm, and the stopping depth 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 ($\Delta t\sim 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 $\approx 10^{32}$ solar flare occur within an observational cadence of 19 s,
which is a factor of $\approx 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) in 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).

##  5 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\approx 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 $\gamma$ 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\approx 9000$ K blackbody
interpretation of stellar flare white-light continuum radiation.

##  Appendix 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 $\Delta
t=0.2$ s (in most cases). Average spectral quantities, such as those in Table
3.1, 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
(CDF) 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.

We thank an anonymous referee for helpful comments and critiques. AFK thanks
Dr. Kosuke Namekata for providing models used in his paper, Dr. Rachel Osten
for discussions about grid applications and nonthermal particles in stellar
flares, Dr. Jim Drake and Dr. Eduard Kontar for conversations about beam
instabilities and low-energy cutoffs, Dr. Pier Emmanuel-Tremblay for providing
the hydrogen broadening calculations, Dr. Thomas Gomez for helpful
conversations about line shape theory, and Dr. Dana Longcope for helpful
discussions about magnetic field retraction. AFK thanks Dr. Graham Kerr and
Dr. Chris Osborne for recommending Plotly. CHIANTI is a collaborative project
involving the University of Cambridge (UK), the NASA Goddard Space Flight
Center (USA), the George Mason University (GMU, USA) and the University of
Michigan (USA).

## References

  * Abbett (1998) Abbett, W. P. 1998, PhD thesis, MICHIGAN STATE UNIVERSITY
  * Abbett & Hawley (1999) Abbett, W. P., & Hawley, S. L. 1999, ApJ, 521, 906, doi: 10.1086/307576
  * Abrevaya et al. (2020) Abrevaya, X. C., Leitzinger, M., Oppezzo, O. J., et al. 2020, MNRAS, 494, L69, doi: 10.1093/mnrasl/slaa037
  * Alaoui et al. (2021) Alaoui, M., Holman, G. D., Allred, J. C., & Eufrasio, R. T. 2021, ApJ, 917, 74, doi: 10.3847/1538-4357/ac0820
  * Allred (2005) Allred, J. C. 2005, PhD thesis, University of Washington, Seattle
  * Allred et al. (2020) Allred, J. C., Alaoui, M., Kowalski, A. F., & Kerr, G. S. 2020, ApJ, 902, 16, doi: 10.3847/1538-4357/abb239
  * Allred et al. (2005) Allred, J. C., Hawley, S. L., Abbett, W. P., & Carlsson, M. 2005, ApJ, 630, 573, doi: 10.1086/431751
  * Allred et al. (2006) —. 2006, ApJ, 644, 484, doi: 10.1086/503314
  * Allred et al. (2022) Allred, J. C., Kerr, G. S., & Gordon Emslie, A. 2022, ApJ, 931, 60, doi: 10.3847/1538-4357/ac69e8
  * Allred et al. (2015) Allred, J. C., Kowalski, A. F., & Carlsson, M. 2015, ApJ, 809, 104, doi: 10.1088/0004-637X/809/1/104
  * Alvarado-Gómez et al. (2018) Alvarado-Gómez, J. D., Drake, J. J., Cohen, O., Moschou, S. P., & Garraffo, C. 2018, ApJ, 862, 93, doi: 10.3847/1538-4357/aacb7f
  * Arnaud & Rothenflug (1985) Arnaud, M., & Rothenflug, R. 1985, A&AS, 60, 425
  * Arnold et al. (2021) Arnold, H., Drake, J. F., Swisdak, M., et al. 2021, Phys. Rev. Lett., 126, 135101, doi: 10.1103/PhysRevLett.126.135101
  * Aschwanden (1998) Aschwanden, M. J. 1998, ApJ, 502, 455, doi: 10.1086/305890
  * Aschwanden (2004) —. 2004, ApJ, 608, 554, doi: 10.1086/392494
  * Aschwanden (2012) —. 2012, ApJ, 757, 94, doi: 10.1088/0004-637X/757/1/94
  * Aschwanden & Alexander (2001) Aschwanden, M. J., & Alexander, D. 2001, Sol. Phys., 204, 91, doi: 10.1023/A:1014257826116
  * Aschwanden et al. (1998a) Aschwanden, M. J., Kliem, B., Schwarz, U., et al. 1998a, ApJ, 505, 941, doi: 10.1086/306200
  * Aschwanden et al. (1995) Aschwanden, M. J., Schwartz, R. A., & Alt, D. M. 1995, ApJ, 447, 923, doi: 10.1086/175930
  * Aschwanden et al. (1998b) Aschwanden, M. J., Schwartz, R. A., & Dennis, B. R. 1998b, ApJ, 502, 468, doi: 10.1086/305891
  * Ashfield & Longcope (2023) Ashfield, W., & Longcope, D. 2023, ApJ, 944, 147, doi: 10.3847/1538-4357/acb1b2
  * Belova & Bychkov (2019) Belova, O. M., & Bychkov, K. V. 2019, Astrophysics, 62, 234, doi: 10.1007/s10511-019-09577-4
  * Bian et al. (2018) Bian, N., Emslie, A. G., Horne, D., & Kontar, E. P. 2018, ApJ, 852, 127, doi: 10.3847/1538-4357/aa9f29
  * Birn et al. (2017) Birn, J., Battaglia, M., Fletcher, L., Hesse, M., & Neukirch, T. 2017, ApJ, 848, 116, doi: 10.3847/1538-4357/aa8ad4
  * Bochanski et al. (2007) Bochanski, J. J., West, A. A., Hawley, S. L., & Covey, K. R. 2007, AJ, 133, 531, doi: 10.1086/510240
  * Brasseur et al. (2023) Brasseur, C. E., Osten, R. A., Tristan, I. I., & Kowalski, A. F. 2023, ApJ, 944, 5, doi: 10.3847/1538-4357/acab59
  * Canfield et al. (1991) Canfield, R. C., Wulser, J.-P., Zarro, D. M., & Dennis, B. R. 1991, ApJ, 367, 671, doi: 10.1086/169663
  * Carlsson (1998) Carlsson, M. 1998, in Space Solar Physics: Theoretical and Observational Issues in the Context of the SOHO Mission, ed. J. C. Vial, K. Bocchialini, & P. Boumier, Vol. 507, 163, doi: 10.1007/BFb0106920
  * Carlsson & Leenaarts (2012) Carlsson, M., & Leenaarts, J. 2012, A&A, 539, A39, doi: 10.1051/0004-6361/201118366
  * Carlsson & Stein (1992) Carlsson, M., & Stein, R. F. 1992, ApJ, 397, L59, doi: 10.1086/186544
  * Carlsson & Stein (1995) —. 1995, ApJ, 440, L29, doi: 10.1086/187753
  * Carlsson & Stein (1997) —. 1997, ApJ, 481, 500
  * Carlsson & Stein (2002) —. 2002, ApJ, 572, 626, doi: 10.1086/340293
  * Carlsson et al. (2023) Carlsson, M., Fletcher, L., Allred, J., et al. 2023, A&A, 673, A150, doi: 10.1051/0004-6361/202346087
  * Chavali et al. (2022) Chavali, S., Youngblood, A., Paudel, R. R., et al. 2022, Research Notes of the American Astronomical Society, 6, 201, doi: 10.3847/2515-5172/ac9591
  * Chen et al. (2020) Chen, B., Shen, C., Gary, D. E., et al. 2020, Nature Astronomy, 4, 1140, doi: 10.1038/s41550-020-1147-7
  * Cho et al. (2022) Cho, P. B., Gomez, T. A., Montgomery, M. H., et al. 2022, ApJ, 927, 70, doi: 10.3847/1538-4357/ac4df3
  * Christian et al. (2003) Christian, D. J., Mathioudakis, M., Jevremović, D., et al. 2003, ApJ, 593, L105, doi: 10.1086/378217
  * Colgan et al. (2008) Colgan, J., Abdallah, J., J., Sherrill, M. E., et al. 2008, ApJ, 689, 585, doi: 10.1086/592561
  * Collier et al. (2023) Collier, H., Hayes, L. A., Battaglia, A. F., Harra, L. K., & Krucker, S. 2023, A&A, 671, A79, doi: 10.1051/0004-6361/202245293
  * Cram & Woods (1982) Cram, L. E., & Woods, D. T. 1982, ApJ, 257, 269, doi: 10.1086/159985
  * Crosley & Osten (2018) Crosley, M. K., & Osten, R. A. 2018, ApJ, 856, 39, doi: 10.3847/1538-4357/aaaec2
  * Dappen et al. (1987) Dappen, W., Anderson, L., & Mihalas, D. 1987, ApJ, 319, 195, doi: 10.1086/165446
  * Del Zanna et al. (2015) Del Zanna, G., Dere, K. P., Young, P. R., Landi, E., & Mason, H. E. 2015, A&A, 582, A56, doi: 10.1051/0004-6361/201526827
  * Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149, doi: 10.1051/aas:1997368
  * Dhillon et al. (2007) Dhillon, V. S., Marsh, T. R., Stevenson, M. J., et al. 2007, MNRAS, 378, 825, doi: 10.1111/j.1365-2966.2007.11881.x
  * Donati et al. (1997) Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658, doi: 10.1093/mnras/291.4.658
  * Dorfi (1998) Dorfi, E. A. 1998, in Saas-Fee Advanced Course 27: Computational Methods for Astrophysical Fluid Flow., ed. O. Steiner & A. Gautschy, 263
  * Dorfi & Drury (1987) Dorfi, E. A., & Drury, L. O. 1987, Journal of Computational Physics, 69, 175, doi: 10.1016/0021-9991(87)90161-6
  * Dorfi & Feuchtinger (1995) Dorfi, E. A., & Feuchtinger, M. U. 1995, Computer Physics Communications, 89, 69, doi: 10.1016/0010-4655(94)00187-7
  * Doyle et al. (1988) Doyle, J. G., Butler, C. J., Bryne, P. B., & van den Oord, G. H. J. 1988, A&A, 193, 229
  * Doyle et al. (2012) Doyle, J. G., Giunta, A., Singh, A., et al. 2012, Sol. Phys., 280, 111, doi: 10.1007/s11207-012-0025-6
  * Eason et al. (1992) Eason, E. L. E., Giampapa, M. S., Radick, R. R., Worden, S. P., & Hege, E. K. 1992, AJ, 104, 1161, doi: 10.1086/116305
  * Egedal et al. (2015) Egedal, J., Daughton, W., Le, A., & Borg, A. L. 2015, Physics of Plasmas, 22, 101208, doi: 10.1063/1.4933055
  * Emslie (1978) Emslie, A. G. 1978, ApJ, 224, 241, doi: 10.1086/156371
  * Emslie (1983) —. 1983, ApJ, 271, 367, doi: 10.1086/161202
  * Emslie et al. (1998) Emslie, A. G., Mariska, J. T., Montgomery, M. M., & Newton, E. K. 1998, ApJ, 498, 441, doi: 10.1086/305551
  * Fang et al. (1993) Fang, C., Henoux, J. C., & Gan, W. Q. 1993, A&A, 274, 917
  * Feinstein et al. (2022) Feinstein, A. D., France, K., Youngblood, A., et al. 2022, AJ, 164, 110, doi: 10.3847/1538-3881/ac8107
  * Fisher (1989) Fisher, G. H. 1989, ApJ, 346, 1019, doi: 10.1086/168084
  * Fisher et al. (1985) Fisher, G. H., Canfield, R. C., & McClymont, A. N. 1985, ApJ, 289, 434, doi: 10.1086/162903
  * Fleishman et al. (2022) Fleishman, G. D., Nita, G. M., Chen, B., Yu, S., & Gary, D. E. 2022, Nature, 606, 674, doi: 10.1038/s41586-022-04728-8
  * Fletcher et al. (2007) Fletcher, L., Hannah, I. G., Hudson, H. S., & Metcalf, T. R. 2007, ApJ, 656, 1187, doi: 10.1086/510446
  * Fletcher & Hudson (2008) Fletcher, L., & Hudson, H. S. 2008, ApJ, 675, 1645, doi: 10.1086/527044
  * Forbes & Priest (1995) Forbes, T. G., & Priest, E. R. 1995, ApJ, 446, 377, doi: 10.1086/175797
  * Froning et al. (2019) Froning, C. S., Kowalski, A., France, K., et al. 2019, ApJ, 871, L26, doi: 10.3847/2041-8213/aaffcd
  * Fuhrmeister et al. (2008) Fuhrmeister, B., Liefke, C., Schmitt, J. H. M. M., & Reiners, A. 2008, A&A, 487, 293, doi: 10.1051/0004-6361:200809379
  * Fuhrmeister et al. (2005) Fuhrmeister, B., Schmitt, J. H. M. M., & Hauschildt, P. H. 2005, A&A, 436, 677, doi: 10.1051/0004-6361:20042518
  * Fuhrmeister et al. (2010) —. 2010, A&A, 511, A83, doi: 10.1051/0004-6361/200810224
  * Gan & Fang (1990) Gan, W. Q., & Fang, C. 1990, ApJ, 358, 328, doi: 10.1086/168989
  * Gan et al. (1992) Gan, W. Q., Rieger, E., Zhang, H. Q., & Fang, C. 1992, ApJ, 397, 694, doi: 10.1086/171825
  * García-Alvarez et al. (2002) García-Alvarez, D., Jevremović, D., Doyle, J. G., & Butler, C. J. 2002, A&A, 383, 548, doi: 10.1051/0004-6361:20011743
  * Gomez et al. (2016) Gomez, T. A., Nagayama, T., Kilcrease, D. P., Montgomery, M. H., & Winget, D. E. 2016, Phys. Rev. A, 94, 022501, doi: 10.1103/PhysRevA.94.022501
  * Graham et al. (2020) Graham, D. R., Cauzzi, G., Zangrilli, L., et al. 2020, ApJ, 895, 6, doi: 10.3847/1538-4357/ab88ad
  * Güdel (2004) Güdel, M. 2004, A&A Rev., 12, 71, doi: 10.1007/s00159-004-0023-2
  * Güdel et al. (2004) Güdel, M., Audard, M., Reale, F., Skinner, S. L., & Linsky, J. L. 2004, A&A, 416, 713, doi: 10.1051/0004-6361:20031471
  * Güdel et al. (1996) Güdel, M., Benz, A. O., Schmitt, J. H. M. M., & Skinner, S. L. 1996, ApJ, 471, 1002, doi: 10.1086/178027
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951, doi: 10.1051/0004-6361:200809724
  * Haggerty et al. (2015) Haggerty, C. C., Shay, M. A., Drake, J. F., Phan, T. D., & McHugh, C. T. 2015, Geophys. Res. Lett., 42, 9657, doi: 10.1002/2015GL065961
  * Hawley et al. (2014) Hawley, S. L., Davenport, J. R. A., Kowalski, A. F., et al. 2014, ApJ, 797, 121, doi: 10.1088/0004-637X/797/2/121
  * Hawley & Fisher (1992) Hawley, S. L., & Fisher, G. H. 1992, ApJS, 78, 565, doi: 10.1086/191640
  * Hawley & Fisher (1994) —. 1994, ApJ, 426, 387, doi: 10.1086/174075
  * Hawley & Pettersen (1991) Hawley, S. L., & Pettersen, B. R. 1991, ApJ, 378, 725, doi: 10.1086/170474
  * Hawley et al. (1995) Hawley, S. L., Fisher, G. H., Simon, T., et al. 1995, ApJ, 453, 464, doi: 10.1086/176408
  * Hawley et al. (2003) Hawley, S. L., Allred, J. C., Johns-Krull, C. M., et al. 2003, ApJ, 597, 535, doi: 10.1086/378351
  * Heinzel & Shibata (2018) Heinzel, P., & Shibata, K. 2018, ApJ, 859, 143, doi: 10.3847/1538-4357/aabe78
  * Holman (1985) Holman, G. D. 1985, ApJ, 293, 584, doi: 10.1086/163263
  * Holman et al. (2003) Holman, G. D., Sui, L., Schwartz, R. A., & Emslie, A. G. 2003, ApJ, 595, L97, doi: 10.1086/378488
  * Hori et al. (1997) Hori, K., Yokoyama, T., Kosugi, T., & Shibata, K. 1997, ApJ, 489, 426, doi: 10.1086/304754
  * Houdebine et al. (1991) Houdebine, E. R., Butler, C. J., Panagi, P. M., Rodono, M., & Foing, B. H. 1991, A&AS, 87, 33
  * Howard et al. (2020) Howard, W. S., Corbett, H., Law, N. M., et al. 2020, ApJ, 902, 115, doi: 10.3847/1538-4357/abb5b4
  * Howard et al. (2023) Howard, W. S., Kowalski, A. F., Flagg, L., et al. 2023, arXiv e-prints, arXiv:2310.03792, doi: 10.48550/arXiv.2310.03792
  * Hubeny et al. (1994) Hubeny, I., Hummer, D. G., & Lanz, T. 1994, A&A, 282, 151
  * Hubeny & Mihalas (2014) Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres
  * Hudson et al. (1992) Hudson, H. S., Acton, L. W., Hirayama, T., & Uchida, Y. 1992, PASJ, 44, L77
  * Hummer & Mihalas (1988) Hummer, D. G., & Mihalas, D. 1988, ApJ, 331, 794, doi: 10.1086/166600
  * Inglis et al. (2015) Inglis, A. R., Ireland, J., & Dominique, M. 2015, ApJ, 798, 108, doi: 10.1088/0004-637X/798/2/108
  * Ireland et al. (2013) Ireland, J., Tolbert, A. K., Schwartz, R. A., Holman, G. D., & Dennis, B. R. 2013, ApJ, 769, 89, doi: 10.1088/0004-637X/769/2/89
  * Jess et al. (2010) Jess, D. B., Mathioudakis, M., Christian, D. J., et al. 2010, Sol. Phys., 261, 363, doi: 10.1007/s11207-009-9500-0
  * Jevremovic et al. (1998) Jevremovic, D., Butler, C. J., Drake, S. A., O’Donoghue, D., & van Wyk, F. 1998, A&A, 338, 1057
  * Johns-Krull et al. (1997) Johns-Krull, C. M., Hawley, S. L., Basri, G., & Valenti, J. A. 1997, ApJS, 112, 221, doi: 10.1086/313030
  * Kanodia et al. (2022) Kanodia, S., Ramsey, L. W., Maney, M., et al. 2022, ApJ, 925, 155, doi: 10.3847/1538-4357/ac3e61
  * Katsova et al. (1997) Katsova, M. M., Boiko, A. Y., & Livshits, M. A. 1997, A&A, 321, 549
  * Katsova et al. (1981) Katsova, M. M., Kosovichev, A. G., & Livshits, M. A. 1981, Astrophysics, 17, 156, doi: 10.1007/BF01005196
  * Katsova et al. (1991) Katsova, M. M., Livshits, M. A., Butler, C. J., & Doyle, J. G. 1991, MNRAS, 250, 402, doi: 10.1093/mnras/250.2.402
  * Kašparová et al. (2007) Kašparová, J., Kontar, E. P., & Brown, J. C. 2007, A&A, 466, 705, doi: 10.1051/0004-6361:20066689
  * Kawate et al. (2016) Kawate, T., Ishii, T. T., Nakatani, Y., et al. 2016, ArXiv e-prints.  https://arxiv.org/abs/1610.04328
  * Kazachenko et al. (2022) Kazachenko, M. D., Albelo-Corchado, M. F., Tamburri, C. A., & Welsch, B. T. 2022, Sol. Phys., 297, 59, doi: 10.1007/s11207-022-01987-6
  * Kerr et al. (2019a) Kerr, G. S., Allred, J. C., & Carlsson, M. 2019a, ApJ, 883, 57, doi: 10.3847/1538-4357/ab3c24
  * Kerr et al. (2019b) Kerr, G. S., Carlsson, M., Allred, J. C., Young, P. R., & Daw, A. N. 2019b, ApJ, 871, 23, doi: 10.3847/1538-4357/aaf46e
  * Kiplinger et al. (1984) Kiplinger, A. L., Dennis, B. R., Frost, K. J., & Orwig, L. E. 1984, ApJ, 287, L105, doi: 10.1086/184408
  * Kontar et al. (2014) Kontar, E. P., Bian, N. H., Emslie, A. G., & Vilmer, N. 2014, ApJ, 780, 176, doi: 10.1088/0004-637X/780/2/176
  * Kontar et al. (2017) Kontar, E. P., Perez, J. E., Harra, L. K., et al. 2017, Phys. Rev. Lett., 118, 155101, doi: 10.1103/PhysRevLett.118.155101
  * Kontar et al. (2012) Kontar, E. P., Ratcliffe, H., & Bian, N. H. 2012, A&A, 539, A43, doi: 10.1051/0004-6361/201118216
  * Kowalski (2022) Kowalski, A. F. 2022, Frontiers in Astronomy and Space Sciences, 9, 1034458, doi: 10.3389/fspas.2022.1034458
  * Kowalski (2023) —. 2023, ApJ, 943, L23, doi: 10.3847/2041-8213/acb144
  * Kowalski (2024) —. 2024, arXiv e-prints, arXiv:2402.07885, doi: 10.48550/arXiv.2402.07885
  * Kowalski & Allred (2018) Kowalski, A. F., & Allred, J. C. 2018, ApJ, 852, 61, doi: 10.3847/1538-4357/aa9d91
  * Kowalski et al. (2022) Kowalski, A. F., Allred, J. C., Carlsson, M., et al. 2022, ApJ, 928, 190, doi: 10.3847/1538-4357/ac5174
  * Kowalski et al. (2017a) Kowalski, A. F., Allred, J. C., Daw, A., Cauzzi, G., & Carlsson, M. 2017a, ApJ, 836, 12, doi: 10.3847/1538-4357/836/1/12
  * Kowalski et al. (2019a) Kowalski, A. F., Butler, E., Daw, A. N., et al. 2019a, ApJ, 878, 135, doi: 10.3847/1538-4357/ab1f8b
  * Kowalski et al. (2015) Kowalski, A. F., Hawley, S. L., Carlsson, M., et al. 2015, Sol. Phys., 290, 3487, doi: 10.1007/s11207-015-0708-x
  * Kowalski et al. (2010) Kowalski, A. F., Hawley, S. L., Holtzman, J. A., Wisniewski, J. P., & Hilton, E. J. 2010, ApJ, 714, L98, doi: 10.1088/2041-8205/714/1/L98
  * Kowalski et al. (2012) —. 2012, Sol. Phys., 277, 21, doi: 10.1007/s11207-011-9839-x
  * Kowalski et al. (2013) Kowalski, A. F., Hawley, S. L., Wisniewski, J. P., et al. 2013, ApJS, 207, 15, doi: 10.1088/0067-0049/207/1/15
  * Kowalski et al. (2018) Kowalski, A. F., Mathioudakis, M., & Hawley, S. L. 2018, in 20th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, 42, doi: 10.5281/zenodo.1463140
  * Kowalski et al. (2016) Kowalski, A. F., Mathioudakis, M., Hawley, S. L., et al. 2016, ApJ, 820, 95, doi: 10.3847/0004-637X/820/2/95
  * Kowalski et al. (2017b) Kowalski, A. F., Allred, J. C., Uitenbroek, H., et al. 2017b, ApJ, 837, 125, doi: 10.3847/1538-4357/aa603e
  * Kowalski et al. (2019b) Kowalski, A. F., Wisniewski, J. P., Hawley, S. L., et al. 2019b, ApJ, 871, 167, doi: 10.3847/1538-4357/aaf058
  * Krucker et al. (2011) Krucker, S., Hudson, H. S., Jeffrey, N. L. S., et al. 2011, ApJ, 739, 96, doi: 10.1088/0004-637X/739/2/96
  * Krucker et al. (2020) Krucker, S., Masuda, S., & White, S. M. 2020, ApJ, 894, 158, doi: 10.3847/1538-4357/ab8644
  * Kundu et al. (2001) Kundu, M. R., White, S. M., Shibasaki, K., Sakurai, T., & Grechnev, V. V. 2001, ApJ, 547, 1090, doi: 10.1086/318422
  * Kunkel (1970) Kunkel, W. E. 1970, ApJ, 161, 503, doi: 10.1086/150555
  * Lalitha et al. (2013) Lalitha, S., Fuhrmeister, B., Wolter, U., et al. 2013, A&A, 560, A69, doi: 10.1051/0004-6361/201321419
  * Lee & Gary (2000) Lee, J., & Gary, D. E. 2000, ApJ, 543, 457, doi: 10.1086/317080
  * Li et al. (2012) Li, T. C., Drake, J. F., & Swisdak, M. 2012, ApJ, 757, 20, doi: 10.1088/0004-637X/757/1/20
  * Li et al. (2013) —. 2013, ApJ, 778, 144, doi: 10.1088/0004-637X/778/2/144
  * Liefke et al. (2010) Liefke, C., Fuhrmeister, B., & Schmitt, J. H. M. M. 2010, A&A, 514, A94, doi: 10.1051/0004-6361/201014012
  * Liefke et al. (2008) Liefke, C., Ness, J. U., Schmitt, J. H. M. M., & Maggio, A. 2008, A&A, 491, 859, doi: 10.1051/0004-6361:200810054
  * Lin & Forbes (2000) Lin, J., & Forbes, T. G. 2000, J. Geophys. Res., 105, 2375, doi: 10.1029/1999JA900477
  * Liu et al. (2013) Liu, W., Chen, Q., & Petrosian, V. 2013, ApJ, 767, 168, doi: 10.1088/0004-637X/767/2/168
  * Livshits et al. (1981) Livshits, M. A., Badalian, O. G., Kosovichev, A. G., & Katsova, M. M. 1981, Sol. Phys., 73, 269, doi: 10.1007/BF00151682
  * Longcope et al. (2016) Longcope, D., Qiu, J., & Brewer, J. 2016, ApJ, 833, 211, doi: 10.3847/1538-4357/833/2/211
  * Longcope et al. (2018) Longcope, D., Unverferth, J., Klein, C., McCarthy, M., & Priest, E. 2018, ApJ, 868, 148, doi: 10.3847/1538-4357/aaeac4
  * Loyd et al. (2018a) Loyd, R. O. P., Shkolnik, E. L., Schneider, A. C., et al. 2018a, ApJ, 867, 70, doi: 10.3847/1538-4357/aae2ae
  * Loyd et al. (2018b) Loyd, R. O. P., France, K., Youngblood, A., et al. 2018b, ApJ, 867, 71, doi: 10.3847/1538-4357/aae2bd
  * MacGregor et al. (2020) MacGregor, A. M., Osten, R. A., & Hughes, A. M. 2020, ApJ, 891, 80, doi: 10.3847/1538-4357/ab711d
  * MacGregor et al. (2018) MacGregor, M. A., Weinberger, A. J., Wilner, D. J., Kowalski, A. F., & Cranmer, S. R. 2018, ApJ, 855, L2, doi: 10.3847/2041-8213/aaad6b
  * MacGregor et al. (2021) MacGregor, M. A., Weinberger, A. J., Loyd, R. O. P., et al. 2021, ApJ, 911, L25, doi: 10.3847/2041-8213/abf14c
  * Mann et al. (2013) Mann, A. W., Brewer, J. M., Gaidos, E., Lépine, S., & Hilton, E. J. 2013, AJ, 145, 52, doi: 10.1088/0004-6256/145/2/52
  * Martínez Oliveros et al. (2012) Martínez Oliveros, J.-C., Hudson, H. S., Hurford, G. J., et al. 2012, ApJ, 753, L26, doi: 10.1088/2041-8205/753/2/L26
  * Mathioudakis et al. (2006) Mathioudakis, M., Bloomfield, D. S., Jess, D. B., Dhillon, V. S., & Marsh, T. R. 2006, A&A, 456, 323, doi: 10.1051/0004-6361:20054752
  * Mathioudakis et al. (1999) Mathioudakis, M., McKenny, J., Keenan, F. P., Williams, D. R., & Phillips, K. J. H. 1999, A&A, 351, L23
  * McClymont & Canfield (1986) McClymont, A. N., & Canfield, R. C. 1986, ApJ, 305, 936, doi: 10.1086/164307
  * Melrose & Brown (1976) Melrose, D. B., & Brown, J. C. 1976, MNRAS, 176, 15, doi: 10.1093/mnras/176.1.15
  * Mochnacki & Zirin (1980) Mochnacki, S. W., & Zirin, H. 1980, ApJ, 239, L27, doi: 10.1086/183285
  * Moffett (1974) Moffett, T. J. 1974, ApJS, 29, 1, doi: 10.1086/190330
  * Monson et al. (2024) Monson, A. J., Mathioudakis, M., & Kowalski, A. F. 2024, arXiv e-prints, arXiv:2401.02261, doi: 10.48550/arXiv.2401.02261
  * Morchenko et al. (2015) Morchenko, E., Bychkov, K., & Livshits, M. 2015, Ap&SS, 357, 119, doi: 10.1007/s10509-015-2347-y
  * Morchenko (2020a) Morchenko, E. S. 2020a, Astrophysics, 63, 91, doi: 10.1007/s10511-020-09617-4
  * Morchenko (2020b) —. 2020b, Astrophysics, 63, 440, doi: 10.1007/s10511-020-09648-x
  * Mullan et al. (2006) Mullan, D. J., Mathioudakis, M., Bloomfield, D. S., & Christian, D. J. 2006, ApJS, 164, 173, doi: 10.1086/502629
  * Murphy et al. (2007) Murphy, R. J., Kozlovsky, B., Share, G. H., Hua, X. M., & Lingenfelter, R. E. 2007, ApJS, 168, 167, doi: 10.1086/509637
  * Namekata et al. (2020) Namekata, K., Maehara, H., Sasaki, R., et al. 2020, PASJ, doi: 10.1093/pasj/psaa051
  * Neidig et al. (1994) Neidig, D. F., Grosser, H., & Hrovat, M. 1994, Sol. Phys., 155, 199, doi: 10.1007/BF00670740
  * Neidig et al. (1993) Neidig, D. F., Kiplinger, A. L., Cohl, H. S., & Wiborg, P. H. 1993, ApJ, 406, 306, doi: 10.1086/172442
  * Notsu et al. (2024) Notsu, Y., Kowalski, A. F., Maehara, H., et al. 2024, ApJ, 961, 189, doi: 10.3847/1538-4357/ad062f
  * Osten et al. (2007) Osten, R. A., Drake, S., Tueller, J., et al. 2007, ApJ, 654, 1052, doi: 10.1086/509252
  * Osten et al. (2006) Osten, R. A., Hawley, S. L., Allred, J., et al. 2006, ApJ, 647, 1349, doi: 10.1086/504889
  * Osten et al. (2005) Osten, R. A., Hawley, S. L., Allred, J. C., Johns-Krull, C. M., & Roark, C. 2005, ApJ, 621, 398, doi: 10.1086/427275
  * Osten & Wolk (2015) Osten, R. A., & Wolk, S. J. 2015, ApJ, 809, 79, doi: 10.1088/0004-637X/809/1/79
  * Osten et al. (2010) Osten, R. A., Godet, O., Drake, S., et al. 2010, ApJ, 721, 785, doi: 10.1088/0004-637X/721/1/785
  * Osten et al. (2016) Osten, R. A., Kowalski, A., Drake, S. A., et al. 2016, ApJ, 832, 174, doi: 10.3847/0004-637X/832/2/174
  * Paudel et al. (2021) Paudel, R. R., Barclay, T., Schlieder, J. E., et al. 2021, ApJ, 922, 31, doi: 10.3847/1538-4357/ac1946
  * Paulson et al. (2006) Paulson, D. B., Allred, J. C., Anderson, R. B., et al. 2006, PASP, 118, 227, doi: 10.1086/499497
  * Penn et al. (2016) Penn, M., Krucker, S., Hudson, H., et al. 2016, ApJ, 819, L30, doi: 10.3847/2041-8205/819/2/L30