(middle left) shows an energetic, gradual-type flare event that exhibits
several faster impulsive events superimposed on a more gradual rise phase; the
impulsive phase of this GF event ends around $t=2.3$ hrs and transitions to a
more gradually decaying phase.

###  6.2 Quasi-Periodic Pulsations

High-time resolution data of stellar flares facilitate searches for
periodicity over the light curve evolution. Some events have shown
fluctuations in the broadband flux having periods with a modulated amplitude
or an evolving period over a cycle or two; these are called quasi-periodic
pulsations (QPPs). This topic is of significant interest within the solar and
stellar communities [749], but an exhaustive discussion and reference list are
outside the scope of this review. I summarize a few selected results. For
detailed discussions and recent reviews, see [91] and [801].

QPPs in stellar flares are reported on a variety of timescales, from 10 s to
30 minutes, and they occur in the impulsive and gradual phases. [487] report a
10 s period during the $u$-band flare in the peak phase using a wavelet
analysis [719] and discuss different interpretations, including one in which
the reconnection and particle acceleration are modulated through MHD
oscillations in a nearby large loop [525]. [16] detect a damped 32 minute
oscillation during the decay phase of a $U$-band megaflare [412] on YZ CMi
using Lomb-Scargle and autocorrelation analyses of the detrended light curve.
Doyle et al [184] discuss shorter, $1-4$ minute, QPPs in two giant ($E\approx
5-10\times 10^{33}$ erg) flares observed on YZ CMi with $0.25-0.6$ s sampling.
The AFINO method [353, 354] is another powerful technique for evaluating the
statistical significance of QPPs. [353] apply the AFINO method to the stellar
event in [16] and Kowalski et al [412] using a Fourier power-spectrum analysis
that includes the component of the flare signal that is detrended in the
wavelet technique. This approach accounts for the fact that some types of
power-spectra may naturally give rise to types of bursty light curve behavior
in the temporal domain.

Stellar flare QPPs are reported in data from the optical to the X-ray regimes
on a wide variety of stars. Cho et al [129] investigate a tight correlation
between the damping times and periods of QPPs in a sample of X-ray stellar
flares. White-light QPPs in _Kepler_ data are presented in [50]. [486] report
remarkable oscillations with a 240 s period during the peak phase of a
$U$-band flare on the RS CVn II Peg, and [504] discuss several physical
origins in the damped 750 s period in the X-ray flare on the lower-mass system
AT Mic. For an analysis of multi-band data during the decay phase of an X-ray
flare on the young G star EK Dra, see [92]. QPPs of 320 s and 660 s periods
were reported over the peak phase of a remarkably energetic flare on a PMS M3
star [360].

##  7 Multi-Wavelength Spectral Observations

Multi-wavelength observations of stellar flares are grouped into categories
according to the regions of the atmosphere producing the electromagnetic
response. A predominantly thermal response in the optical and near-ultraviolet
originates from the cool, $T\approx 10^{4}$ K, dense lower atmosphere (Sect.
7.1 – 7.2), while the far-ultraviolet emission lines probe the rapid response
of transition region temperatures around $T\approx 10^{5}$ K (Section 7.3).
Nonthermal radiation originates from accelerated particles gyrating in coronal
magnetic fields (Sect. 7.4 – 7.5), and a thermal response results from the
hot, $T\gtrsim 10^{7}$ K, tenuous upper atmosphere (Sect. 7.6). The temporal
correlations (Sect. 7.7) among these radiative responses justify a solar-like
modeling paradigm with nonthermal electron beams and chromospheric evaporation
and condensation processes (Sects. 3 and 8.2).

###  7.1 NUV and Optical: The Thermal Radiative Response of the Footpoints to
Impulsive Nonthermal Heating from Above

The optical and NUV response is one of the most enigmatic and energetic
aspects of stellar flares. This phenomenon is known for its dramatic impulsive
phase, but it also exhibits a long-duration gradual decay phase that can
extend for many hours after a bright peak flux phase. Long-duration continuum
flare radiation may persist after relatively small peaks as well [see the GF1
event in Kowalski et al 414 reproduced in the middle left panel of Fig. 8 and
the event studied in 298]. The optical and NUV radiation is thought to be
thermal radiation that originates from photospheric and/or chromospheric
heights. For reference, Table 7 summarizes several of the largest broadband
optical flux enhancements in well-studied MVe flares with multi-band
photometry data that include the $U$\- or SDSS $u$-band.

Table 7: Flux enhancements at peak phase for several especially large MVe
events with simultaneous multi-band photometry. The flux enhancement values
($I_{f}+1$; Sect. 2) are reported relative to the pre-flare or quiescence in
each band. The flux enhancements are related to magnitude changes through
Equation 1.

(row) Star | SpT |  $E_{U/u}$ (erg) | $U$ | $u$ | $B$ | $g$ | $V$ | $r$ | $R$ | $I$ | $i$ | 3500† | 4170† | 6010†
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
(1a) YZ CMi | M4.5Ve | $\approx 4\times 10^{33}$ |  |  |  |  |  |  |  |  |  | 30.9 | 7.6 | 1.7
(1b) YZ CMi | M4.5Ve | $\approx 4\times 10^{33}$ |  |  |  |  |  |  |  |  |  | 106 | 28.7 | 3.8
(2) YZ CMi | M4.5Ve | $1.85\times 10^{33}$ |  | 78.2 |  | 6.3 |  | 2.6 |  |  | 1.35 |  |  |
(3) AU Mic | M1(V)e | $1.6\times 10^{33}$ | 3.4 |  |  |  | 1.13 |  |  |  |  |  |  |
(4) AD Leo | M3Ve | $6.5\times 10^{33}$ | 70 |  | 11 |  | 3.4 |  | 1.6 |  |  |  |  |
(5) EV Lac | M3.5Ve | $3.9\times 10^{32}$ | 40 | 40 |  | 4 |  |  |  |  |  |  |  |
(6) EV Lac | M3.5Ve | $6.2\times 10^{32}$ | 8.0 | 7.4 |  | 1.4 |  | 1.18 |  |  |  |  |  |
(7) SDSS J001309 | M6Ve | … |  | 158 |  | 16.9 |  | 8.6 |  |  | 2.0 |  |  |

2929footnotetext: (1a) and (1b): The maxima over the IF1 and IF3 flares,
respectively, from [416]; the combined, estimated $U$-band energy is given
because they are part of the same “Ultraflare” event. (2): The IF3 event from
[414]; (3): Flare # 23 from [722]; (4): The Great Flare from [297]; (5): The
biggest flare from [654] and the most impulsive event (although labeled as
IF10) in the sample of [414]; (6): The GF1 event from [414] that is shown in
the middle left panel of Figure 8. †The 3500, 4170, and 6010 columns refer to
the central wavelengths of the ULTRACAM narrow-band filters (NBF). (7) SDSS
J001309.33-002552 from the SDSS Stripe 82 flare star sample [411] and in Table
1.

####  7.1.1 Overview of Emission Line Properties

An overview of the temporal response of various emission lines and continuum
fluxes in the $U$ band and optical regimes throughout a high-energy flare
($E_{u}>10^{33}$ erg; Table 7) is shown in Fig. 9. The temporal sequence in
the figure is representative of many stellar flares, although quantitative
differences among events occur [329, 327, 414, 420]. The line fluxes respond
rapidly in the impulsive phase, except that Ca II K (and H) peak well into the
gradual decay phase of the optical flare continuum flux. The late peak of Ca
II K is a remarkable phenomenon in stellar flares [297, 330, 327, 251, 414].
The blue continuum flux is the fastest to rise to its maximum and then begin
its decay, followed by the highest order Balmer lines [297]. The Balmer
H$\alpha$ line is the slowest to decay relative to its peak [see also 527]. In
some flares, the Balmer line peaks occur significantly after the peak of the
continuum, which may be empirically related to secondary events [e.g., 297,
251, and see also the IF4/F2 event in Kowalski et al 416]. The neutral helium
lines exhibit a very wide range of decay timescales: The blue He I lines decay
quickly [297], but He I 10830 Å (studied in a different MVe flare) exhibits
one of the most gradually evolving light curves [654].

Quantitative assessments of the time-scales of emission line and continuum
fluxes have been investigated with a time-decrement [414, 420]. The time-
decrement is defined as the empirical trend of the full-width-at-half-maximum
($t_{1/2}$) of each light curve, plotted as a function of the wavelength of
the transition. Are the relative decay timescales in Fig. 9 consistent with a
simple model consisting of a slowly decreasing average temperature over the
flare footpoints on the star [289, 329]? This explanation would be analogous
to the cooling thermal loop sequence in solar flares [28] that is attributed
to the temporal offsets in the maximum of optically thin, XEUV light curves.
The average temperature properties of the flare region are undoubtedly
important in understanding stellar flare light curve evolution in Fig. 9.
However, the ratios of the Balmer line fluxes (which is called the Balmer flux
decrement, or just “Balmer decrement”) and line-to-continuum flux ratios in
the observations suggest that chromospheric optical depths [428, 193],
vertical and transverse spatial inhomogeneities [415, 408], and non-LTE
radiative transfer [295, 9] are probably also critical to interpret the
sustained chromospheric emission line fluxes (and the relative evolution of
the optical continuum fluxes) in stellar flares. In summary, there are not yet
any comprehensive physical models that explain the absolute and relative
timescales of the various optical spectral variations in Fig. 9.

Figure 9:  Representative time-evolution of optical spectral quantities over a
giant flare on the M4.5Ve star, YZ CMi. The peak optical broadband magnitude
enhancements are listed in row (2) of Table 7. The emission line fluxes are
integrated over wavelength, and all light curves are normalized to their
respective peak fluxes. The C4170′ light curve is a proxy for the flare-only
blue continuum flux averaged over $\lambda=4155-4185$ Å and is the fastest to
decay after its peak flux is attained. The H$\alpha$ and Ca II K emission
lines are the slowest. Figure reproduced from [414] with permission.

###  7.2 Spectral Properties of the Optical and NUV Continuum Radiation

The continuum fluxes through the NUV and optical are clearly much more
impulsive than the emission lines in Fig. 9 (see also [297, 251, 420]). The
observed flare continuum radiation from the NUV through the optical is
collectively referred to as the white-light303030The $\lambda=2000-4000$ Å
ultraviolet spectral region can be roughly divided into the NUV
($\lambda=2000-3200$ Å), and the $U$\- (or SDSS $u$-) band
($\lambda=3200-4000$ Å). The UV-C covers $\lambda=2000-2800$ Å, while the FUV
typically refers to wavelengths between $\lambda=1000-2000$ Å. The blue-
optical refers to $\lambda=4000-4800$ Å, the green-optical to
$\lambda=4800-5800$ Å, and the red-optical to $\lambda=5800-7000$ Å. Sometimes
the $U$-band is considered part of the optical because it can be observed from
the ground. The $\lambda=1800-2400$ Å range is sometimes referred to as the
mid-UV. In reference to stellar flare radiation, white-light has traditionally
meant “detected in broadband optical filters”; before the immaculate precision
provided by Kepler, this unambiguously referred to flare continuum radiation.
In solar physics, a common anecdotal use of “white-light flare” refers to an
enhancement of photospheric continuum radiation and thus the excitation of the
solar atmosphere in very deep layers. Adding to the ambiguity, there may be
two types of white-light solar flares, Type I and Type II, which are
classified according to their optical and Balmer jump spectral properties
[see, e.g., 598]. The introduction of [419] provides further discussion about
the ambiguities in terminology surrounding the use of “white-light” in solar
physics and astrophysics. Some researchers colloquially refer to the white-
light (“flare” implied) as an increase in the optical continuum radiation only
(excluding shorter wavelength $U$-band and NUV continuum radiation). Note that
bona-fide “white light” is defined according to the CIE D65 white-light
standard illuminant, which is noon sunlight. See [141] for human color vision
synthesis of various stellar colors. Note, we have synthesized the “true”
colors in Fig. 31 from representative spectra of each source that is
illustrated.. In this review, we adopt the empirical definition (without
regard to its origin in the stellar atmosphere) that a white-light flare is a
broad-wavelength increase in the observed NUV, $U$-band, and optical continuum
stellar flare radiation that sometimes extends into the FUV and NIR. Thus, the
white light would contribute a majority of the integrated flux in optical
broadband filters (e.g., $UBVR$, Kepler, TESS).

In this section, we review observational properties of optical and NUV flare
continuum radiation, as revealed by spectroscopic measurements. We separately
discuss peak/impulsive phase spectra (Sect. 7.2.1) and gradual decay phase
spectra (Sect. 7.2.2). In Sect. 7.2.4 and Sect. 8.2, we review some closely-
related results from analyses of colors in broadband photometry
(“colorimetry”). The M-dwarf flare spectral observations, combined with
detailed modeling, suggest that there are unexpected amounts of heating over a
large (deep) column mass density (hereafter, just “column mass” or $m_{c}$; g
cm-2). Whether the $\lambda=1300-9000$ Å continuum radiation in stellar flares
is caused by heating the photosphere to bona-fide incandescence (i.e.,
isotropic blackbody radiation or a blackbody-like spectral intensity) is an
open question that is discussed in the next section and further in Sect. 8.2.
The data discussed in this review come mostly from dMe flares, which produce
large contrasts against their non-flaring photospheric radiative fluxes at
blue and optical wavelengths. This property facilitates isolating “flare-only”
spectra, which can be readily compared to models without the ambiguities of
subtraction artifacts.

####  7.2.1 Rise and Peak Phase

A commonly reported empirical property of the optical, $\lambda\geq 4000$ Å,
flare-only continuum radiation in the spectra of the impulsive phase of
stellar flares is a color temperature matching a $T\approx 8,000-14,000$ K
blackbody [505, 297, 380, 579, 414, 265, 434, 416]. The blackbody color
temperatures might be purely phenomenological, in which case the color
temperature values are parametrizations of the spectral shapes rather than
analogs to in situ “thermometer readings” (Appendix C.1). A temperature value
is simply easier to relate to than continuum flux ratios, which are unique to
specific spectral windows for a given color temperature. The continuum
radiation that is consistent with $T\approx 10^{4}$ K optical blackbody color
temperatures is thus referred to as “hot, blackbody-like”. As we discuss
further in Sect. 8.2, there are possible physical explanations from optically
thick hydrogen recombination radiation. A peak-phase spectrum is shown in Fig.
10, which was taken at the peak of the event whose evolution is showcased in
Fig. 9. Single-component blackbody function models in the range of $T\approx
9000-12,000$ K can be fit to this flare spectrum at $\lambda\geq 4000$ Å
outside of the major emission lines [414]. Multi-component blackbody curves
with higher and lower temperatures are fit to the entire $\lambda\gtrsim 4150$
Å wavelength range and are compared to a single blackbody fit to the blue-
optical, $\lambda=4000-4800$ Å, range in Fig. 10.

Figure 10:  Flare spectrum at the peak of the large, impulsive-type (IF) event
on YZ CMi in Fig. 9. Several blackbody fits to optical ($\lambda\geq 4000$ Å)
continuum regions are shown. The flare spectrum (black) has the pre-flare
(dotted) spectrum subtracted. Figure reproduced from [414] with permission.

The broadband photometry observations from Hawley and Pettersen [297]
established the 9500 K blackbody hypothesis for the continuum in an energetic
$E\sim 10^{34}$ erg flare from the M3Ve star AD Leo (see Sect. 8.2 for
additional discussion), which was verified in [414] by fitting a blackbody
temperature of $T\approx 11,600$ K to the optical range in the spectra of this
flare. Further spectral investigation at resolving powers of $R\sim 500-1000$
in other events has revealed that a $T_{\rm{BB}}\approx 8500-14,000$ K
blackbody color temperature in the blue-optical, $\lambda=4000-4800$ Å, is not
a unique property of extremely energetic, large amplitude flares: for example,
similar color temperatures were calculated in the impulsive phase flare
spectrum (Appendix A) of the IF event in Fig. 8(bottom), which has an energy
of $E_{U}\approx 1.6\times 10^{31}$ erg [416]. Hot blackbody-like continuum
radiation in the blue-optical regime tends to be more prominent, or more often
reported, in the rise and peak phases of impulsive-type (IF; Sect. 6.1)
events. The flare event in Fig. 10 is a highly impulsive event. We refer the
reader also to the discussions in Kowalski et al [416] and Kowalski [409]
pertaining to the variation in the optical blackbody color temperatures on
$\Delta t=3$ s cadences in the large, IF-type events listed in rows (1a) and
(1b) in Table 7.

In echelle flare spectra with high-resolving power, Schmitt et al [661] and
Fuhrmeister et al [241] have shown that the continuum radiation within the $U$
band at $\lambda=3250-3860$ Å also exhibits hot blackbody color temperatures
of $\approx 11,300$ K in the impulsive phase (Fig. 11a). At lower resolving
power, similarly blue spectral trends at wavelengths shortward of the Balmer
limit have been noted in Kowalski et al [414] for a large sample of flares.
Though spectral observations within the $U$-band benefit from increased flare
contrast, especially afforded by M dwarfs, accurate measurements are generally
very difficult due to the systematic uncertainties of the flux calibration
near the terrestrial atmospheric limit.

Spectral observations with high signal-to-noise (S/N) at $\lambda>3600$ Å
reveal a very important clue about the origin of the white-light continuum
radiation and its hot blackbody-like property. An extrapolation of a blackbody
function that is fit to $\lambda=4000-4800$ Å (thus giving $T_{\rm{BB}}$), or
an extrapolation of a blackbody that is fit to the continuum flux ratio,
$f^{\prime}_{4170\rm{\AA}}/f^{\prime}_{6010\rm{\AA}}=$C4170${}^{\prime}/$C6010′
(thus giving $T_{\rm{FcolorR}}$), fails to account for excess continuum flux
at wavelengths shorter than $\lambda\lesssim 3700$ Å. In other words, there is
a Balmer jump between the blue-optical continuum flux and continuum flux at
shorter wavelengths within the $U$ band. An isothermal blackbody function thus
does not comprehensively explain the panchromatic continuum flux properties. A
Balmer jump was tentatively noted in Kunkel [428] and Moffett [506], and it
was clearly evident and extensively characterized in impulsive-phase flare
spectra using modern instrumentation and high-time resolution in Kowalski et
al [414, 416]. An impulsive phase spectrum during a $E_{U}\approx 10^{31}$
erg, gradual/hybrid-type flare on GJ 1243 is shown in Fig. 11(b) over a
broader spectral range than in panel (a). This is one of two flares analyzed
in [420] that produced similar decreases in the continuum flux into the NUV
range at $\lambda<3200$ Å, as constrained by Hubble Space Telescope data (not
shown here). The Balmer jump in this event is the largest that has been
detected spectroscopically over the impulsive phase of an M dwarf flare. In
flares with very prominent hot blackbody-like continuum radiation at optical
wavelengths, such as the Great Flare of AD Leo and the flare in Fig. 10, the
Balmer jumps are evident as smaller flux excesses above the blackbody
extrapolations to $\lambda\lesssim 3700$ Å.

Figure 11:  (a) VLT/UVES echelle ($R\sim 40,000$) spectral observations at
$\lambda<3860$ Å that nearly extend to the terrestrial atmospheric cutoff.
These spectra show the evolution of a giant white-light flare on the M6Ve star
CN Leo, reproduced from [241] with permission. A $T\sim 11,300$ K blackbody is
fit in the impulsive phase (2nd from top), and the gradual phase spectrum is
fit with a $T\sim 9100$ K blackbody in the third panel. A catalog of the
emission lines in the spectrum in the middle panel, extending to
$\lambda=3060$ Å, is available through VizieR. (b) Other MVe flares in the
impulsive phase show a large Balmer jump and a flatter $U$-band continuum
spectrum that decreases into the NUV at shorter wavelengths than shown here
[420] are not satisfactorily explained by a $T\approx 9000$ K blackbody fit to
the blue-optical continuum wavelengths. Here, the two blue-colored spectra
correspond to the rise/peak phase and early gradual decay phase, respectively,
which are averaged into the black spectrum. The red-colored spectrum is an
optically thin hydrogen recombination continuum model from a RHD simulation
that invokes a solar-type electron beam heating function; the model is scaled
to the observed spectrum at $\lambda\approx 3615$ Å. The spectra have a low-
resolving power, $R\approx 600$, which is clearly sufficient for robust
wavelength-integrated emission line fluxes. Figure reproduced with permission
from Kowalski et al [420].

Impulsive-phase continuum flux ratios, which are “colors” if they are
converted to magnitude differences, are calculated from flare spectra in
[414], [416], and [420] and are shown in Fig. 12. Flux ratios have been
calculated from isothermal static, slab models for $\tau(\lambda)\ggg 1$
(blackbody) and $\tau=0$ (continuous hydrogen LTE emissivity; Appendix C) and
are shown in Fig. 12 for comparison. The flare peak “color-color” diagrams
from the spectra demonstrate that the Balmer jump strengths
($f_{3615}^{\prime}$/$f_{4170}^{\prime}=$C3615${}^{\prime}/$C4170′) imply that
moderate-to-large continuum optical depths develop in stellar flare
atmospheres. Further quantitative understanding about the observed offsets
from the blackbody line in Fig. 12 requires detailed RHD model spectra, which
are represented in the figure by the thick dashed line. These RHD model
calculations include wavelength-dependent opacities and optical depths over
column mass/depth-dependent temperatures and electron densities (Sect. 8.2),
which are not included in blackbody or other simple slab models.

In Appendix B, we discuss comparable color-color diagrams using narrow-band
($\Delta\lambda=50-100$ Å) filters around $\lambda=3500$ Å, $4170$ Å, and
$6010$ Å, which were designed specifically to avoid emission lines in flares
and to measure the Balmer jump with the ULTRACAM instrument [169]. These
narrow-band filters provide high-cadence ($\Delta t\approx 0.3-3$ s),
simultaneous constraints of continuum variations while avoiding degeneracies
in model fits to broadband filters (Sect. 7.2.4). ULTRACAM count-flux ratios
from a large sample of MVe flares were analyzed and converted into flare-only
continuum flux ratios in [416]. The flare color from relative photometry is
the ratio of two values of

$I_{f}(t)\langle f_{q,\lambda}\rangle_{T}$ (12)

in the notation of Equation 7. Several flares were observed simultaneously
with spectra over the impulsive to gradual decay phases, which facilitated
unambiguous interpretation of the filter ratios (_cf_ Figure 3 in Kowalski et
al 416). A relationship between color temperatures fit to blue-optical spectra
and to blue-to-red optical flux ratios is characterized by an offset of
$\Delta T\approx 1900$ K (see Table 5 in Kowalski et al 416), in line with the
various blackbody fits to the peak flare spectrum in Fig. 10.

Figure 12: Compilation of impulsive-phase continuum flux ratios (color-color
diagrams) from spectra [414, 416, 420]. The y-axis is a measure of the Balmer
jump ratio, and the flux ratios on the x-axis correspond to optical color
temperatures. Note that the $x$ and $y$ uncertainties are correlated, but only
marginal error bars are shown. The prime symbols indicate that a pre-flare
flux has been subtracted. The optical continuum colors from the Great Flare of
AD Leo [297] are estimated from extrapolations of blackbody curve fits over
the $\lambda=3800-4440$ Å range [414]. The wavelength range of C4170′ can be
seen in the right panel of Fig. 11. The flux ratio values at the peaks of
flares indicate optical depths between blackbody radiation and optically thin
hydrogen recombination radiation. The general trend is well-reproduced by
approximations [KA18 CC models; 410] to the dynamic atmospheres in RHD
simulations (Sect. 8.2). In the right margin plot, flares with Balmer jump
measurements only are shown. The flares are color-coded to their $U$-band
impulsiveness (IF/HF/GF) classification. Other flare-peak flux ratios from
narrowband ULTRACAM photometry [416] are summarized in Appendix B. The color
temperatures corresponding to the blue-to-red flux ratios in this figure are
typically $\gtrsim 1500$ K less than that obtained from spectral fitting over
the $\lambda=4000-4800$ Å range. All data in this figure are provided in
supplementary online material hosted on Zenodo.

####  7.2.2 The Decay Phase and The Evolution of Optical Spectral Properties

The $T\approx 10,000$ K blackbody-like optical color temperatures do not last
for long during MVe flares [e.g., 414, 409]. After the peak of a broadband
(e.g., $U$-band) light curve, a gradual decay phase almost always follows a
shorter period of faster decay. A characteristic blue-optical, blackbody color
temperature value at the start of the gradual decay phase is
$T_{\rm{BB}}\approx 8000$ K [414]. For equal $\lambda\approx 3615$ Å flare-
only ﬂux levels during the fast rise and fast decay phases, the fast decay
phase is more than $\Delta T\approx 1000$ K cooler. [505] observed similar
trends for optical color temperatures using continuum data during several MVe
ﬂares at $\lambda=4200-6900$ Å (e.g., see their Figure 3). At red-optical
wavelengths, the continuum becomes very flat in the gradual decay phase, with
even cooler temperatures. This was referred to as a red “conundruum” in [414];
the cooler blackbody fits were linked to the increasing fraction of energy
that was previously reported in the redder bands of photometry in the decay
phase of the Great Flare of AD Leo [see Figure 10 and accompanying discussion
in 297] and retrospectively to the deviations from blackbody fits to the
photometry of that flare (bottom right panel of Figure 11 in Hawley and Fisher
295). See Figure 31 of [414] for examples of the stellar flare gradual phase
continuum spectra that exhibit cool $T\lesssim 5000$ K blackbody color
temperatures in large dMe events. In smaller flares, it is difficult to
accurately subtract the preflare spectrum at $\lambda\gtrsim 7500$ Å [420] and
thus determine the NIR continuum properties.

It is not clear how to explain such persistent flat optical continua with
relatively moderate-sized Balmer jump ratios; optically thin (bound-free)
Paschen continua and photospheric backwarming predict cooler color
temperatures but also much larger Balmer jumps, which are not observed – hence
the conundrum. A related phenomenon during the impulsive phase may be the
cause of the lower color temperatures that are calculated over a very wide
optical wavelength range (such as those on the x-axis in Fig. 12) that spans
the blue-optical and red-optical regimes. On the other hand, [241] report
hotter blackbody temperatures over the red-wavelengths (not shown) in the
flare in Fig. 11a. Clearly, further characterization of the red-optical and
near-IR continuum properties is warranted. There may be important
contributions to _Kepler_ flare energetics in the gradual decay phase [304].

At shorter wavelengths within the $U$-band, the Balmer jump ratio increases in
the gradual decay phase of a flare. The fraction of the wavelength-integrated
energy in emission lines (primarily the hydrogen Balmer series in the optical
and within the $U$-band wavelength range) also becomes larger. The energy
partition evolution is shown in Fig. 13 for a large sample of dMe flares. It
was found that the impulsive flare (IF) events were more continuum-dominated
than the gradual flare (GF)-type events in both the impulsive phases and
gradual phases of each event. Thus, an empirical connection among the
broadband impulsiveness, fraction of energy in the Balmer component, Balmer
jump strength (right panel of Fig. 13), and Balmer line-to-continuum ratios
(namely, the H$\gamma$/C4170′ ratios – see Appendix B) was established from
simultaneous spectral and broadband photometry observations [see also 416,
420].

Stellar flares in the NUV generally exhibit a rapid time evolution [87].
However, the spectral properties in the NUV ($\lambda=2000-3200$ Å) are not
nearly as well-observed as in the optical, in either impulsive- or gradual-
type flares. The NUV spectral observations and model predictions are compiled
and reviewed in [88]. [420] investigated the NUV properties in two HF-type
events (that show some optical spectral properties that are more in line with
other GF-type events) with shorter-wavelength spectra ($\lambda=2440-2840$ Å)
from the Hubble Space Telescope / Cosmic Origins Spectrograph. They concluded
that a Balmer continuum enhancement upon a blackbody-like optical continuum
component is necessary to account for the continuum flux down to at least
$\lambda\approx 2500$ Å. International Ultraviolet Explorer (IUE) FUV spectra
covering $\lambda=1900-3100$ Å have been reported in the decaying phases of
two large events on MVe stars [297, 630]. For a review of the timing of the
IUE/NUV and IUE/FUV spectral integrations within the AD Leo Great Flare, see
[408]. Of the few (relatively small) flares with NUV spectra that have been
analyzed through their peak and decay phases, the fractions of energy in the
continuum were calculated in the large range of $\approx 50-90$% [303, 420],
which is much larger than in quiescence (35%). The evolution of the
$\lambda\approx 3300-3800$ Å continuum-to-line energy fractionation is also
noticeable after the peak phase in the M5.5/6Ve event in Fig. 11(a).

Figure 13:  Energy of the hydrogen Balmer component vs. the optical continuum
energy evolution in a large sample of MVe flares (reproduced from Kowalski et
al 414 with permission). The Balmer jump ratios ($\chi_{\rm{flare}}$) increase
into the gradual decay phase and are also generally correlated with the
IF/HF/GF classification scheme. The hydrogen Balmer component includes an
estimate for the Balmer continuum flux above an optical continuum model
extrapolation to $\lambda<3650$ Å. Several more flares with comparable
spectroscopy [416, 420] since this result have been consistent with these
trends.

####  7.2.3 Empirically-determined, Broadband Radiated Energy Budgets Relative
to the $U$ band

Empirical radiative energy fractionation from the X-ray to the radio provides
model-independent313131In some flare models of radiative backwarming [e.g.,
537], the energy in one band is reprocessed energy from other bands.
characterization of stellar flares. The energy in the optical and NUV
continuum provides a large contribution to the fractionation. To date, [562]
is the most comprehensive analysis of the multi-wavelength energy budget in
archival observations of stellar flares across a wide variety of spectral
types and quiescent activity levels. At $\lambda=1200-8000$ Å, 96% of the peak
flare luminosity can originate from the continuum in the impulsive phase
[297]. In the gradual decay phase, the continuum accounts for as much as 83%
when the emission lines remain highly elevated for longer durations (e.g.,
Fig. 9). The energy radiated in soft X-rays (0.04 - 2 keV) is estimated to be
$\sim 11$x the H$\gamma$ energy [104] and to be comparable to that in the
$U$-band, which is roughly 1/6 of the integrated $1200-8000$ Å white-light
energy [297, 298], and a factor of 1/10 of the bolometric flare energy [562].
[297] constrain $E_{U}=0.64E_{2000-3260\mathrm{\AA}}$, however the bandpass-
averaged energy [erg $\rm{\AA}^{-1}$] in the $U$-band is reported as a factor
of 1.2 greater. Several estimates of the $U$-band to Kepler-band (Kp) energy
ratios have been empirically determined to be $E_{U}=(0.4-0.65)E_{\rm{Kp}}$
[304]. [433] determine $E_{U}=1.2E_{B}=1.8E_{V}$. Line-integrated flare
energies in H$\alpha$ are typically factors of $\gtrsim 0.04-0.08$ smaller
than the $U$-band energies, and the H$\gamma$ energies are factors of $\approx
0.08$ of the $U$-band energies over several orders of magnitude in total flare
energy [see the large compilation in Figure 4.2 in 406].

The relative energy in soft X-rays ($E\approx 0.1-5$ keV) is difficult to
compare among flares. This difficulty arises due to different temporal [e.g.,
298] and wavelength integration limits (e.g., by including a significant
portion of the EUV range) among datasets. The X-ray and optical thermal
radiation occur on much different time-scales (Sect. 7.7), and a limiting
assumption [e.g. 38, 39] is sometimes that the broad-band flare energy is
accurately accounted for with a scaled 2-temperature quiescent coronal model
[281]. Nonetheless, many interesting empirical constraints on the X-ray to
optical flare energies are reported in the literature. The ratio of the soft
X-ray (and/or EUV) to the $U$-band energy can take on a large range of values,
from essentially no detectable response in the soft X-rays within a $\pm
10-20$ minute window [177, 567]), down to $1/3$ [298], to as large as a factor
of $\gtrsim 10$ [283, 571]. [527] study the X-ray properties (Sect. 11) of
several X-ray flares with NICER spectra. The $E=0.5-10$ keV soft X-ray
energies fall in the range of $0.3-1\times 10^{32}$ erg, but the event with
complete multi-wavelength coverage in their study showed an H$\alpha$ response
without a detectable $g^{\prime}$-band (white-light) response. In general,
using the energy fractionation relation of
$E_{U}=1/3E_{0.01-10\rm{keV}}=1/3E_{\rm{XEUV}}$ in [562], one obtains
consistent cumulative flare rates in the XEUV and in the $U$-band [e.g., $\sim
1$ per day greater than $10^{32}$ erg in $U$ 590] corresponding to $1$ flare
per day with $E_{\rm{XEUV}}\gtrsim 3\times 10^{32}$ erg [39] for non-
simultaneous flares on AD Leo. [722] calculate the energy ratios of six
simultaneous $U$-band and soft X-ray ($E=0.2-12$ keV) flares to be
$E_{\rm{SXR}}/E_{U}\approx 1.5$. [662] extend the comprehensive, multi-
wavelength analysis of a $E\approx 10^{34}$ erg flare on the K0Ve star AB Dor
from [434]; they also report comparable energies released in soft X-rays and
in an NUV bandpass. This is consistent with a paradigm that a majority [562,
431] of the total bolometric response in many (but not all) white-light
stellar flares occurs from energy deposition in the dense, chromospheric
footpoints (Fig. 2(right)).

The peak radio (3.6 cm) specific luminosity (erg s-1 Hz-1) is $\sim 1/4$ of
the $U$-band, but the bandpass-integrated energy is a factor of
$\mathcal{O}(10^{4})$ smaller [567]. New broadband observations have suggested
much larger energy releases than expected323232TESS results are not included
in this iteration of the review, but we should note briefly the discrepancies
with models are also discussed in [334] and [363]. in the NUV regime relative
to the optical [88], while others are rather well-matched by current RHD
models that invoke large heating rates [571, 420]. Studies of high-cadence,
simultaneous Galex/FUV and Galex/NUV photometry have reported much larger peak
luminosities in the broadband FUV than in the broadband NUV [632] (but see the
discussion of bright Galex sources in Sect. 7.3). In the NIR, Davenport et al
[156] present model predictions for the magnitude enhancements and energies
relative to the $U$ band. [714] used high-cadence, high-precision monitoring
of $\Delta U=-1.5$ mag flares on mid-type M dwarfs to constrain the broadband
NIR response to $|\Delta_{\rm{NIR}}|\lesssim 5-12$ milli-magnitudes.

####  7.2.4 Optical Broadband Flare Colorimetry

Optical flare colorimetry is the study of the temporal tracks of broadband
flare colors or absolute fluxes. For example, the $U-B$ vs. $B-V$ colors of
flares on the dM3.5e star EV Lac were compared to regimes of optically thick
and optically thin hydrogen spectra in [797]. They conclude that the flare
continuum becomes optically thick (and more blackbody-like) in the peak phase
of flares, but in the decay phase, the optical depths decrease. [302] find a
best-fit blackbody temperature of $T\sim$8500 K by fitting to the fluxes in
$UBVR$ filters for a sample of eight moderate-sized flares on AD Leo.
Colorimetry with the traditional Johnson/Bessell bandpasses is a powerful and
convenient method to test a wide variety of flare continuum models [428, 518,
522, 140, 179, 295, 468, 528, 88], and similar methods using broadband filters
in the NUV and FUV have inferred a wide range of even hotter blackbody
temperatures, $T\approx 20,000-50,000$ K [632, 259].

However, there are several important assumptions in broadband colorimetry that
warrant spectroscopic verification or calibration [e.g., 416]. We discuss a
few here. First, a correction to the contribution from emission lines in the
optical bandpasses must be considered [297, 295]. The $V$ band is relatively
free of major flare emission lines, but the contributions from the Balmer
lines vary from flare to flare in the $B$ band (e.g., Fig. 13, Fig. 33 in
Appendix A). The $U$ band integrates over the Balmer jump, a pseudo-continuum
of blended Balmer lines, and to a lesser extent Ca II K and H emission. [9]
demonstrated that a RHD model spectrum with a large Balmer jump and a
discontinuity at $\lambda=3646$ Å exhibits colors that are consistent with a
$T\sim 9000$ K blackbody when convolved with the $UBVR$ filters – but a
blackbody function has no Balmer discontinuity or jump! [420] convolved an
observed flare spectrum with broadband filters and showed that very large
blackbody temperatures $\approx 15,000-22,000$ K could be inferred if the
Balmer jump and high-order Balmer lines are not included in a model for the
$U$-band and/or a constraint on the NUV continuum flux is not available (see
their Figure 16). They showed that the full optical continuum regime may be
flatter than a hot blackbody fit to the blue-optical wavelengths in such
flares. The analysis of [302] included a bona-fide continuum flux constraint
from HST spectra in the FUV, which may have otherwise resulted in much hotter
peak blackbody temperatures, as in [797].

The results from broadband flare colorimetry have largely motivated
spectroscopic investigations [261, 297, 414] of flares and detailed modeling
[295], which will be discussed in Sect. 8.2.

###  7.3 FUV Observations: The Rapidly Evolving Flare Transition Region (TR)

The FUV continuum and emission line response may be a complex mixture of
thermal, nonthermal, and nonlocal radiative processes in the flare transition
region ($T\gtrsim 50,000$ K) and the deeper atmosphere at $T\lesssim 10,000$ K
[297, 591, 302, 43, 237, 408, 182, 183]. Many questions remain about the
origin of the spectral and temporal properties in stellar flares in this
wavelength regime. A few observational properties are summarized in this
section.

Hawley et al [302] presented a comprehensive, high-resolving power
($R=70,000$) study of eight moderate-sized M dwarf flares in the FUV with
contemporaneous optical spectroscopy and broadband photometry. The energy
budgets and timing properties in the brightest spectral lines and continua
were compared in detail. In the biggest flare, the peak line fluxes in C
II$\lambda$1335, Si IV$\lambda$1394, C IV$\lambda$1548, and C III$\lambda$1176
were calculated in the ratio of $\approx$ 25:35:50:100, respectively (note
that the lineshifts in this study are discussed in Sect. 10). The study also
presented a novel timing analysis of the light curves of Si IV, C IV, N V, and
the FUV continuum region (averaged over $\lambda\approx
1266-1295,1420-1452,1675-1710$ Å) compared to the $U$-band by fitting $b$ and
$m$ in the relation $\log_{10}$($F_{X}$) $=b+m$ $\log_{10}$($F_{U}$). For C IV
1551 and Si IV 1403, the slopes are rather consistent with a value of $\approx
1$, but the value of $m_{\rm{FUVcont}}\approx 1.7$ indicates a significantly
more rapid evolution of the FUV continuum compared to the $U$-band, an effect
that has also been noted in a much larger event on AD Leo [297]. The power-law
relationship between the FUV continuum and the $U$ band is not yet explained
by physical models. A linear scaling among neutral silicon FUV continua, C IV,
and C II luminosities among a sample of dMe and RS CVn flares was interpreted
in [591] to be consistent with radiative backheating of the temperature
minimum region.

Remarkable M-star FUV flare spectra were analyzed in [297], [610], [465],
[237], [471], and [216]. In nearly all cases, the continuum spectra are either
very flat [in units of erg cm-2 s-1 Å-1; see also the spectra compiled in 103,
90, 107, 591] or an isothermal blackbody continuum fit to the regions between
emission lines results in $T\gtrsim 15,000$ K, up to $T\approx 30,000-40,000$
K. Of these, [297] and [471] had simultaneous multi-wavelength observations
that were analyzed in detail, and [237] and [465] discuss redshifted line
asymmetries in the FUV. [43] discuss interesting FUV continuum bursts during
the gradual decay phase of a Si IV light curve in an energetic event from the
young solar analog EK Dra (Fig. 14); similar responses during the FUV bursts
do not show up in higher temperature lines. These continuum-only bursts were
interpreted as possible stellar analogs of “Type II” white-light solar flares
[210].

Figure 14:  Figure 7 from [43], reproduced with permission, showing various
FUV spectral quantities from HST/COS spectra during the decay phase of a large
flare on the young solar analog, EK Dra. The bottom $x$-axis is elapsed time
(hr) and top $x$-axis is rotational phase. Note the bursts that only occur in
the FUV continuum ($\lambda=1435$ $\pm$ 25 Å, 1610 $\pm$ 25 Å). Later in the
decay phase, shorter FUV wavelengths were observed by HST/COS, providing
constraints on the continuum at $\lambda=1344.5\pm 6$ Å and $\lambda=1381.0\pm
8.0$ Å, the forbidden coronal line [Fe XXI]$\lambda$1354
($T_{\rm{form}}\approx 10^{7}$ K), and C II$\lambda$1335\.

Broadband comparisons of FUV and NUV responses during flares have been
possible due to the time-tagged photon capability of the Galex mission. [632]
analyzed a giant flare in both Galex bands, and [767] showed similar trends in
other flares: very large ratios of the FUV to NUV fluxes that are far larger
than expected from any known continuum model constrained by optical
spectroscopy (Sect. 7.2). It has been pointed out that non-linearity
corrections are important considerations at these large count rates in Galex
[516, 228]. More recently, the gPhoton tool [502] has been employed to study
events in the NUV without such concerns at large peak count rates, revealing
remarkably fast time variations on a wide variety of stars [87]. Simultaneous
constraints from Kepler and Galex observations favor very large flux ratios
that are, again, far larger than RHD models of optically thick continuum
radiation from high-energy electron-beam heating models [88].

Potential signatures of proton beams are found in both the FUV and in the
gamma rays. [459] and [691] searched for gamma-ray stellar flares, which
result from the higher energy ($E\gtrsim 30$ MeV / nucleon) ions and protons
interacting with the stellar chromosphere. Lower-energy protons ($E\ll 1$ MeV)
in the beams are expected to undergo charge-exchange with the neutral hydrogen
in the chromosphere before it explodes in temperature. The nonthermal hydrogen
beams then emit red-shifted Lyman $\alpha$ photons as a captured electron
falls to the ground state. This is called the Orrall-Zirker effect [558, 388].
[777] reported a broad satellite component that was highly redshifted from
Ly$\alpha$ and was attributed to low-energy proton beams in the early phase of
a stellar flare in AU Mic. Follow-up observations with high-time resolution
spectroscopy of the Ly $\alpha$ red and blue wings have not resulted in
similar detections [629, 631, 216].

###  7.4 Centimeter (Microwave/Radio) Observations: Nonthermal Gyrosynchrotron
Radiation from Mildly Relativistic Electrons

The radio/microwave regime at cm wavelengths (GHz frequencies) is a probe of
nonthermal emission from accelerated electrons in flare coronae. Lower
frequencies, $\lesssim 2$ GHz, tend to be highly circularly polarized, shorter
in duration, and narrowband, indicating plasma or electron-cyclotron maser
emission, while higher frequencies are largely unpolarized, broadband, and
longer in duration. The unpolarized higher frequency radiation is consistent
with gyrosynchrotron radiation from a power-law distribution of mildly
relativistic electrons trapped in the magnetic fields of the flare region. The
peak frequency of the gyrosynchrotron spectrum determines the demarcation at
which the spectrum transitions from optically thin (at frequencies higher than
the peak frequency) to optically thick (at frequencies lower than the peak
frequency). For a homogeneously emitting source, the power-law index of
accelerated electrons is directly related to the power-law index of the
optically thin gyrosynchrotron flux spectrum at Earth; see [195] and Sect. 11
for more discussion and references.

Radio stellar flares have typically been observed in two bands, such as 6 cm
(4.9 GHz) and 3.6 cm (8.4 GHz), which are close to or at lower frequencies
than typical peak frequencies ($\approx 10$ GHz) [note that some solar flares
have peak frequencies constrained to at least $\gtrsim 30$ GHz; 771]. A large
multi-wavelength flare from the dM3.5e star EV Lac was studied in [567] and a
light curve is shown Fig. 15. The peak radio luminosity at 8.4 GHz is
$1.9\times 10^{15}$ erg s-1 Hz-1, but the $U$-band (not reproduced here;
$L_{U,\rm{peak}}\approx 10^{30}$ erg s-1, $E_{U}\approx 5\times 10^{31}$ erg)
peaks about 54 s earlier and decays much more rapidly. The spectral index
($\alpha;S_{\nu}\propto\nu^{\alpha}$) evolution is discussed in [567] and
indicates rapidly varying optical depths at these frequencies over the flare;
in the decay phase, the higher frequency likely becomes optically thin, while
in the peak phase, a reasonable range of expected power-law indices for the
nonthermal electrons cannot explain the values of $\alpha$ for a homogeneous
optically thick source. The lengthening decay timescales were interpreted in
terms of the “trap+precipitation” model (Section 3). Source parameters were
derived during the long decay, indicating very large loops on comparable
scales to the stellar radius with small trapped nonthermal electron densities
($10^{4}-10^{6}$ cm-3), and large magnetic fields (100 G) increasing over time
and decreasing in source size as the higher frequency becomes more optically
thin (R. Osten, priv. communication 2011).

Fig. 16 shows three radio flares from the comprehensive multi-wavelength study
of the RS CVn system HR 1099 [566]. The data show low polarization during the
flares and a similar spectral index evolution patterns among these long-
duration flares. [280] and [687] present radio analyses of stellar flares that
will be discussed below (Sect. 7.7, Sect. 11) in multi-wavelength contexts.

Figure 15:  A remarkable radio flare on the dM3.5e EV Lac from [567] with the
spectral indices calculated from observations at the two observed frequencies.
Reproduced with permission.  Figure 16:  Flux (top row), polarization (middle
row), and spectral indices (bottom row) during three gyrosynchrotron radio
flares from the RS CVn HR 1099, reproduced from [566] with permission.

The emission at lower frequencies results from a plethora of possible
phenomena. Decimetric flares333333[136] review extreme solar decimetric events
and physical processes. can be very luminous and highly circularly polarized
indicating coherent emission from beam-plasma instabilities (“plasma
emission”) or electron cyclotron maser [437, 56, 560, 750, see references and
discussion in Bastian 55, Güdel et al 280, Osten et al 567]. Coherent bursts
can occur at higher frequencies too, and a remarkable example is the 100%
left-hand circularly polarized 30 mJy peak-flux, 3-minute burst at 5 GHz
discussed in [687] from the dM3e AD Leo. Assuming plasma emission, they derive
a coronal electron density of $\approx 3\times 10^{11}$ cm-3; assuming
electron cyclotron, they derive a magnetic field of $B\approx 1.8$ kG. Between
1–2 GHz, stellar analogies to Type III decimetric bursts343434The “plasma
emission” in this case is produced by excitation of Langmuir waves and
conversion to electromagnetic radiation as electron beams travel over the
density gradient of the corona. have been reported in very high time-
resolution data from Arecibo [559]. A surprising lack of Type II radio bursts
has been constrained in the context of extensive ground-based monitoring of
optical flaring [144]. Many recent numerical, theoretical, and observational
efforts have sought to explain the lack of Type II radio bursts from other
stars [14, 521, 15, 776, see also Cliver et al 135].

###  7.5 Millimeter Observations from ALMA: Nonthermal Radiation with Linear
Polarization

Bright nonthermal flare radiation occurs at much higher frequencies around
$\nu\approx 230$ GHz ($\lambda\approx 1.3$ mm) too. Millimeter flares have
been reported in ALMA data from AU Mic and Proxima Centauri [470, 469].
Particularly remarkable is the extremely impulsive flare from Proxima Centauri
with multi-wavelength observations [471]. The dMe mm flares are short in
duration ($2-30$ s), have rather symmetric (Gaussian-like) impulsive-phase
time profiles, and exhibit only a very weak gradually decaying tail. They note
that these properties are rather similar to the FUV continuum response in the
same flares (see Sect. 7.3). The spectral indices, $\alpha$
($S_{\nu}\propto\nu^{\alpha}$), are negative with a significant amount of
linear polarization that “flip-flops” between negative and positive values
during the peak phase.

The millimeter flare emission is interpreted as either an optically thin
extension to the gyrosynchrotron radiation component that peaks around
$\nu=10$ GHz (see the previous section), or as optically thin synchrotron
radiation that undergoes some type of depolarization to the observed levels of
$\pm 20$% [see 469]. The radiation sources are relativistic or ultra-
relativistic electrons, but the magnetic field strengths are largely different
between these hypotheses. Being in the optically thin regime of either
interpretation, these observations have facilitated important inferences of
the power-law indices of the accelerated electrons in the range of
$\delta_{r}\approx 2.8-5.2$, indicating rather hard distributions of
accelerated particles. At much higher frequencies, [66] reported short (FWHM
durations of $\lesssim 1$ s) synchrotron bursts in the optical $U$-band
superimposed on a giant, unpolarized dMe event from UV Ceti. Only lower limits
on the linear polarization were possible, but the analyses of the energetics
of the bursts favored very hard, $\delta<3.4$, power-laws of ultra-
relativistic electrons.

Flares at $\nu\approx 100$ GHz ($\lambda\approx 3$ mm) have been reported from
PMS stars, RS CVn’s, and eccentric binaries with colliding magnetospheres
[592, 82, 647, 648, 100, 4]. In the case of the remarkable periastron event in
the binary T Tauri V773 A [480], the flare was interpreted as synchrotron
radiation. There is a notable observation of an RS CVn event at $\nu\leq 100$
GHz showing a break and a rise toward higher frequencies [57]: i.e., an
opposite spectral index to the dMe events reported in ALMA, and one that is
more similar to some reports in solar flares [383]. Colliding magnetospheres
at periastron are also attributed to some flare-like optical and NUV
variability in the PMS DQ Tau system [648, 259], while the optical variability
properties at periastron at other times follow the expectations of accretion
alone [482, 715].

###  7.6 Soft X-rays and EUV Observations: Loops Filled with Evaporated,
Previously-Chromospheric/TR, $T>10^{7}$ K Thermal Plasma

Stellar flares have been studied in detail in soft X-rays ($E\approx 0.1-10$
keV) for many decades through XMM-Newton, Chandra, Swift, and NICER
spectroscopy and photometry. X-ray flare spectra are a composite of many
optically thin, highly-ionized emission lines and free-free bremsstrahlung
continua that originate from the flaring coronal volume. The bulk of the
radiation is thought to originate in the chromospheric gas that has been
heated to millions of K and ablated into the magnetic loops, which confine the
flows (Sect. 3). One of the highest signal-to-noise flare spectra during a
stellar flare, which occurred in an RS CVn system, is showcased in Fig. 17.

An extensive review of stellar X-ray flares and analysis methods (prior to and
including 2004) are provided by [276] and [277]. A general review of X-ray
spectroscopy of stars is contained in [279]. Since 2004, there have been many
important X-ray flare studies [567, 633, 697, 240, 545, 546, 573, 346, 448,
243, 574, 434, 599, 527, 578, 375]. In this section, we focus on a few key
results from the last several decades pertaining to abundance changes and the
stellar flare Neupert effect. Several other results are included throughout
this review in other contexts where they are especially complementary (e.g.,
multi-wavelength analyses).

Figure 17:  X-ray spectra from Chandra during quiescence (top) and flaring
(bottom) in the active binary $\sigma^{2}$ Coronae Borealis, reproduced from
[565] with permission. The authors note the enhanced continuum flux increasing
towards shorter wavelengths in the flare spectrum.

###  7.7 The Neupert Effect in Stellar Flares

The backbone of the solar-stellar flare connection is the Neupert effect
[540]. This was originally established in three solar flares as a correlation
between the cumulative time-integral of the nonthermal microwave emission
around 2.7 GHz and the thermal soft X-ray flux at 1.87 Å, attributed to the Fe
XXV line emission (1s2–1s2p, which is sensitive to $T\approx 30$ MK). It has
since been the subject of many follow-up studies in solar flares using either
the hard X-rays or the derivative of the GOES 1–8 Å soft X-rays as a proxy for
the nonthermal particle heating in the chromosphere [e.g., 165, 495, 739, 740,
761]. Alternatively, the empirical Neupert effect is expressed as a
relationship between integral of the hard X-ray flux and the luminosity of the
GOES soft X-rays. The relationship between the impulsive, early nonthermal
radiation and the gradual, later-peaking thermal flux is thought to be linked
by the physics of impulsive nonthermal particle propagation and bombardment of
the chromosphere, rapid evaporation/ablation of $T\approx 10-50$ MK
chromospheric mass into the coronal loops, and the eventual cooling to pre-
flare coronal temperatures (Sect. 3).

The Neupert effect is also reported in multi-wavelength observations of
stellar flares. Because stellar flare hard X-ray ($E\gg 10$ keV) emission is
almost always too faint to detect, and radio observations are relatively
difficult to plan and secure, the (thermal) response in the $U$-band is
utilized as a proxy. This approach is justified by the close temporal and
spatial correspondence between hard X-ray and white-light emission on the Sun
[e.g., 370, 534, 345, 535, 231, 425, 426, 392]. A remarkable example of the
stellar Neupert effect is shown for a flare from Proxima Centauri [283] in
Fig. 18.

A stellar flare Neupert effect was first reported in [298] using the
correlation between the extreme-ultraviolet (65–190 Å) flux and the integrated
$U$-band energy during a long-duration, 3.5-hour rise phase of a flare from
the dM3e AD Leo. This was followed shortly thereafter in [280] using a
correlation between the $E=0.5-2$ keV soft X-rays and 3.6 cm radio emission
during a flare on the dM5.5e UV Ceti A. The Neupert effect has subsequently
been reported in stellar flares from RS CVn systems [284, 566], dMe stars
[283, 285, 503, 758, 243, 108, 722], dK-stars [434], T Tauri stars [40]. There
are also notable exceptions to the expected multi-wavelength relations from
the Neupert effect in stellar flares [177, 567]; see these papers for
discussions of possible explanations, which include deep heating from ultra-
relativistic electrons or MeV protons. On the Sun, a clear violation of the
Neupert effect in a “late impulsive [hard X-ray] burst” [adopting the
terminology from 786] has been modeled with a large, low-energy cutoff
$E\gtrsim 100$ keV in the nonthermal electron beam [761], which would
reasonably generate small amounts of collisional heating within and
chromospheric evaporation into the corona [however, the explanation is still
being investigated; 6].

The Neupert effect in Flare D in [280] was investigated in detail and was
compared to a similar multi-wavelength relationship in a solar flare from
[165]. Intriguingly, the microwave response of flare D returned to pre-flare
levels well before the maximum of the soft X-ray light curve, which is not in
strict agreement with the expected Neupert effect. It was shown that similar
deviations in the solar flare could be explained by a temperature sensitivity
on the Neupert effect, such that radiation from hotter plasma more strongly
followed the expectations [see also 495]. [280] thus separately derived a
_generalized Neupert effect for the light curve_ and a _generalized Neupert
effect_ , following energy conservation over the flaring coronal volume, $V$:

$\frac{d}{dt}\Big{(}3n_{e}k_{B}TV\Big{)}=\alpha
f_{\rm{radio}}(t)-n_{e}^{2}V\psi(T)$ (13)

where $n_{e}$ is the ambient/thermal electron density of the heated plasma at
a temperature $T$, $k_{B}$ is Boltzmann’s constant, $f_{\rm{radio}}$ is the
flare radio flux as a function of time $t$, $\alpha$ is a proportionality
constant, and $\psi(T)$ is the optically thin radiative loss function in units
of erg cm3 s-1. The left hand side of Equation 13 is the change of the thermal
energy of the corona due to chromospheric mass evaporation, and the right hand
side is the balance between radio/ microwave flux at Earth (assumed to be
proportional to the kinetic energy of nonthermal electrons causing the
chromospheric evaporation) and the optically thin radiative cooling of the
coronal volume [assuming conductive cooling is overall small; see also 221].
Thus, accounting for the total plasma energetics and the bolometric luminosity
evolution more naturally explains the apparent deviations from the Neupert
effect assessed in a specific X-ray bandpass. An additional finding from the
solar-stellar flare comparison was that the M dwarf flare was “radio-
overluminous”, possibly suggesting a greater efficiency of acceleration of
microwave-emitting ($E\gtrsim 100$ keV) electrons (Sect. 7.4) and/or a lower
efficiency of chromospheric evaporation, among several other possibilities
[see 280].

Figure 18:  A remarkable example of the Neupert effect in a stellar flare in
Proxima Centauri; figure reproduced with permission, from [283]. The $U$-band
light curve is a (thermal) proxy for the nonthermal, impulsive energy
deposition into the lower atmosphere. This flare was subsequently studied and
modeled in [285] and [608].

The long rise times of Ca II K and H in dMe flares are suggestive of a
Neupert-like effect, which could provide insight into their origin [e.g.,
through XEUV radiative backwarming; 295]. [414] investigated the relationship
between the time-integral of the blue continuum flux around $\lambda=4170$ Å
and the Ca II K line flux for a large sample of dMe flares. They found a
Neupert-like relationship for a wide variety of flare types, both impulsive
and gradual-type events, but none of these observations were complemented with
X-ray data.

###  7.8 Abundance Changes

Impulsive heating in the chromosphere evaporates gas into the coronal loops,
but all chemical elements are not evaporated equally. This results in relative
abundance changes in the hot, coronal flare plasma such that elements with low
first ionization potential (FIP) of $<10$ eV, like Fe, Si, and Mg, become
preferentially enhanced above their photospheric values. The fractionation may
be explained by ponderomotive forces [435] that selectively act on singly-
ionized species in the upper pre-flare chromosphere.

Optically thin, coronal equilibrium modeling of X-ray and EUV spectra has been
used to search for the so-called FIP effect (or its inverse, the IFIP effect)
in stellar flares. Detections of changes in low-FIP abundances have been
suggestive but not highly statistically significant [282], as have relative
changes in high-FIP abundances such as Ne [567]. An example of an abundance
analysis of an RS CVn flare in Fig. 17 is shown in Fig. 19 where an abundance
change in Fe was seen, but other expected FIP dependencies were not. Recently,
[578] searched for the (I)FIP effect in individual flares from the dM3.5e star
EV Lac; they discuss the importance of field geometry in the context of the
Laming theory. A homogeneous Sun-as-a-star (i.e., stellar) EUV spectral
analysis of a large sample of solar flare Fe lines revealed coronal abundances
that were remarkably close to photospheric, thus suggesting that chromospheric
heating and evaporation occur deeper in the atmosphere than the standard
models of electron beam heating and elemental fractionation [763], but [436]
briefly contests the interpretation. Nonetheless, these solar results provide
tantalizing connections to the deep heating rates (Sect. 8.2) that are
inferred and hydrogen line blueshifts (Sect. 10) that are reported from
optical observations of M dwarf flares.

Figure 19:  Abundances relative to hydrogen _vs._ the first ionization
potential (FIP) in quiescent (black) and flaring (red) intervals in Chandra
spectra of the RS CVn $\sigma^{2}$ Coronae Borealis. Figure reproduced from
[565] with permission. They note that there is a slight enhancement of Mg, Si,
O, and Ne during the flare, and Fe increases by nearly twice.

##  8 Atmosphere Modeling of Stellar Flares: Overview of Results

The physics of stellar flares is a complicated (and therefore interesting!)
intersection of particle and nuclear physics, radiative transfer, plasma
processes, gas dynamics and shock phenomena. All of the physical processes are
inter-dependent and occur in a highly magnetized, coupled atmosphere, which is
probably also highly turbulent. Flares are transient non-equilibrium
departures of the stellar luminosity, and all layers of the atmosphere are
thought to produce a response to some degree. In this section, we review
different stellar flare modeling techniques. We focus on efforts toward
understanding chromospheric flare processes in the context of the response of
the entire stellar atmosphere. In contrast to purely magneto-hydrodynamic
models, radiative models include detailed, time-dependent spectral
predictions, which are critical for connecting to observations.

Stellar flare modeling can be sub-divided into several types of techniques.
Isothermal, static slab modeling and atmospheric modeling are the two most
general categories. Atmospheric modeling can be further divided into a static
semi-empirical approach [e.g., 140, wherein the atmospheric structure is
adjusted by hand to match observations] and radiative-hydrodynamic (RHD)
forward modeling, whereby the input heating function is data-constrained and
is not fine-tuned by the modeler to match observations [e.g., 9]. A hybrid,
semi-empirical RHD approach is also possible, whereby the injected flare
heating function is varied by hand over a large parameter space to find a
time-dependent spectral response that matches the observations; thus, the
atmospheric response and depth-dependent heating rates are calculated in a
fully self-consistent manner [e.g., 415]. Though atmospheric modeling accounts
for the vertical heterogeneity at a location in a flare source, additional
techniques are required to model the lateral heterogeneity of the flare source
and the source geometries (to be discussed further in Sect. 12).

###  8.1 Slab Models

Uniform, static slab modeling techniques include the following: fitting a
Planck function to spectra or multi-band photometry [298, 302, 797, 414],
calculating LTE continuum spectra for a range of optical depths
$\tau(\lambda)=[0,\infty)$ [428, 197], and modeling the NLTE Balmer line
decrements, optical depths, and continuum flux ratios [193, 364, 251, 512].
The pioneering study of Kunkel [428] found evidence for a prominent Balmer
jump in their spectral observations of M dwarf flares. Their models of
optically thin hydrogen (emissivity) continuum calculations, however, were
largely inconsistent with the colors of M dwarf flares. They concluded that an
additional moderately heated photospheric flare component could better account
for observations (it was not until the radiative-hydrodynamic models of Allred
et al [9] that essentially produced this scenario through a chromospheric
hydrogen recombination spectrum and a radiatively backheated photospheric
spectrum; see below). In Appendix C, we update the LTE hydrogen continuum
emissivity calculations from Kunkel [428] with non-ideal occupational
probabilities around the Balmer limit and extend them to higher temperatures
and shorter wavelengths. These are useful for verifying for oneself how poorly
first-principle simplifications hold up against continuum observations of
actual stellar flares353535Note, the Balmer jump ratios and optical continuum
flux ratios from these slab models are plotted in Figure 12..

The higher quality spectral observations of the Great Flare of AD Leo [297]
were found to (apparently) exhibit little to no Balmer jump in the flux, but
there was unprecedented power in the optical continuum radiation and highly
broadened hydrogen line wings. The broadband colors from the FUV through the
optical ($R$ band) were tested against all of the models available at the time
[295], and a $T\approx 8500-9500$ K blackbody was clearly favored over a
hydrogen recombination model and a hot $T=10^{7}$ K free-free slab model
interpretation; the former model has a slope that is too red, and the latter
interpretation is inconsistent with X-ray to optical energy fractionations
[see Section 4.3.2 of Hawley and Fisher 295, 562, and Sect. 7.2.3 here] and
representative volume emission measures (Sect. 11). The various
interpretations (blackbody _vs._ optically thin Balmer continuum) and reported
characteristics of the properties at the Balmer limit eventually motivated
additional spectral observations [414]. More immediately, they paved the way
toward advances in the realism of models of stellar flares. It became
essential to consider the comprehensive atmospheric response that includes the
complexities of depth-dependent heterogeneities of opacities, temperatures,
and densities. Then, self-consistent non-equilibrium calculations of velocity
fields, mass advection, radiative cooling, and flare heating became possible.
We summarize these efforts next.

###  8.2 Atmospheric Models

A traditional approach to modeling stellar flare spectra is to start with a
hydrostatic equilibrium model atmosphere and semi-empirically adjust the
temperature gradients of the chromosphere, photosphere, and the location of a
transition region. [140] explored the results of this technique using six
atmospheric variations, including two very extreme adjustments where the
chromosphere and transition region are placed at large column mass depths.
They successfully reproduced a $T\approx 14,000$ K blackbody-like continuum
spectral property in one of these models (their model #5), which also had very
broad hydrogen lines and deep central reversals. Similar approaches have been
widely adopted to model EUV, optical, and infrared stellar flare spectra [326,
489, 365, 131, 242, 654]. An example of this modeling technique from [242] is
shown in Fig. 20, where the chromospheric temperature adjustments result in
satisfactory fits to the emission lines in the giant CN Leo flare (Fig. 11(a))
from [241]. Further, a depth-dependent filling factor and missing photospheric
heating were inferred through this technique.

Figure 20:  Demonstration of semi-empirical temperature adjustments used to
model the red- and blue-wavelength optical spectra of a giant flare from CN
Leo. Figure from [242] reproduced with permission from the author. See also
[243] and [239] for similar approaches for modeling flares from Proxima
Centauri and the dM6 star AZ Cnc, respectively.

How to create these deep transition region locations and chromospheres self-
consistently with a flare heating function, $Q_{\rm{flare}}$? One early
hypothesis [522] was that intense irradiation of soft, thermal X-rays from the
flare corona penetrated to the low atmosphere, causing the optical emission
lines and continuum in the gradual decay phase (whereas, the impulsive phase
optical emission was assumed to be caused by particle beams, as was known to
be consistent with hard X-ray and white-light timing in solar flares). [295]
made a major advancement in stellar flare X-ray backheating models by
calculating the full atmosphere, including the coronal structure. The apex
temperatures were set and all else were determined self-consistently through
energy balance (X-ray heating and photoionization _v._ radiative cooling) and
hydrostatic equilibrium. Although the models produced hydrogen lines that were
too narrow, and the optical and NUV continuum radiation was too weak compared
to the observations [297], a large amount of broad, Ca II K line flux was
predicted (perhaps offering clues to understanding its very slow temporal
response). Updates to XEUV irradiation modeling (e.g., the atomic physics and
geometrical assumptions therein) were employed in evolving atmospheres [296,
1, 8, 9, 10] that facilitated comparisons to the role of nonthermal electron
heating. The role of reprocessing much larger X-ray emissivities into lower
atmosphere optical continuum radiation in superflares has been recently
reconsidered in static atmosphere calculations [543].

Significant advances were made in modeling the temporal evolution of the
atmospheric response in gas-dynamic simulations, starting in the 1980s [458,
491, 224, 223, 222, 203]. An accurate treatment of energy deposition in a
partially ionized, nonuniform chromosphere due to an injected nonthermal
electron distribution was developed [296] based on analytic formulae for
Coulomb collisions [199, 200]. With the inclusion of radiative-transfer, this
effort resulted in a series of papers on radiative363636The gas-dynamic
models, of course, include radiative cooling in the energy equation, but it is
treated using optically thin, coronal equilibrium approximations.-hydrodynamic
modeling of solar flares [2, 8] using a version of the RADYN code [113, 114,
115] that was configured for a larger surface gravity, a cooler photosphere,
and a hotter and denser corona of an M dwarf in Allred et al [9]. The RHD
simulations are currently limited to one spatial dimension, which facilitates
resolving short-timescale shock phenomena, non-equilibrium ionization rates of
hydrogen and helium, and detailed radiative transport.

The Allred et al. stellar flare models studied the atmospheric response to
electron beam heating inputs from the first widely studied solar flare
(SOL2002-07-23T00:30 GOES class X4.8) of the RHESSI [451] hard X-ray
spectroscopic imaging era [167]. A double-power law distribution with a low-
energy cutoff of $E_{c}=37$ keV was obtained from the hard X-ray modeling of
[319], and several injected fluxes were investigated. These models reproduced
the broad Balmer lines for the first time and bright optical and NUV continuum
radiation. However, the Balmer jump was even more conspicuous than in the
energy equilibrium models with X-ray flare irradiation [295] and semi-
empirical models [242]. The dominant mode of radiative backheating in the RHD
models takes place through irradiation of the photosphere by the Paschen and
Balmer continuum fluxes of optically thin thermal radiation in response to the
nonthermal electron beam energy deposition in the mid-to-upper chromosphere
(the beam electrons do not penetrate to the photosphere; see Sect. 9). The
detailed RHD spectrum provided an alternative explanation to the $T\approx
9000$ K blackbody interpretation of broadband photometry reported in earlier
observations [e.g., 302], and the general shape of the Balmer continuum in the
$U$-band was later found to be a satisfactory model in the decay phase spectra
of a megaflare from YZ CMi [412]. The lower flux model was run for several
minutes in [9], and it well-reproduced several transition region flare lines
from [302]. Thus, the Allred et al. models provided a comprehensive solar-
flare modeling framework for interpreting many of the observed properties of M
dwarf flares.

Yet, re-analyses of archival spectra and new flare spectra and narrowband
ULTRACAM photometry still did not match the shape and relative strength of the
optical continuum radiation predicted by the highest-beam flux model in the
[9] simulations. Some important details of the line broadening and blending at
the Balmer limit remained unexplained as well (Sect. 10.2.1). The M dwarf gas-
dynamic models of Livshits et al [458] produced a chromospheric condensation
using a large nonthermal electron flux in an initially isothermal atmosphere
and suggested that a large optical depth in the continuum could be produced
well above the photosphere [see also 381]. [249] investigated chromospheric
condensations as a source of white-light in solar flares. Impulsive-phase
chromospheric condensations have been thought to be important in generating
emission line spectral features in solar flares, such as the red-wing
asymmetries in H$\alpha$ [349, 421, 528] and in more optically thin lines such
as Fe II in IRIS NUV spectra [418, 274]. However, stellar radiative-
hydrodynamic models did not produce chromospheric condensations that were
dense enough to explain the observed continuum radiation until recently when
the effects of very large electron beam fluxes were explored in RADYN [415]
and the RH code [726]. (A realistic M dwarf atmosphere has a dense
chromosphere and very efficient radiative cooling that can make the onset of
complete ionization and explosive hydrodynamics a rather energetic threshold
to attain, especially for a beam with a moderately large low-energy cutoff of
$E_{c}=37$ keV, as traditionally employed in the RADYN dMe flare models.) The
large electron beam fluxes were actually motivated by scaling the radiative
backwarming fluxes from the Allred et al [9] models, but when the large beam
fluxes were simulated, the lower chromosphere become too optically thick for
the NUV and optical continuum radiation to penetrate to the pre-flare
photosphere. Instead, the $\tau(\lambda)$ surfaces shift to a dense
chromospheric condensation and the beam-heated layers just below the
condensation (Fig. 2(right)).

A snapshot of the lower flare atmosphere in a RADYN model from [415] is shown
in Fig. 21. The electron beam flux that is injected at the model loop apex is
very large, $10^{13}$ erg cm-2 s-1 (which is hereafter referred to as a beam
flux of “F13”, whereas “F12” refers to an injected beam flux of $10^{12}$ erg
cm-2 s-1, and so on). By the time shown ($t=2.2$ s), the chromospheric
condensation has cooled and accrued mass into a narrow $\Delta z\approx 15$ km
region with a maximum electron density of $n_{e}\approx 5\times 10^{15}$ cm-3.
This was an exciting development because the model produced a small Balmer
jump ratio and an optical color temperature that was in line with the
impulsive phase spectral observations of dMe flares. In these models, the
emergent optical radiative flux spectrum is a result of thermal response to
the nonthermal electron beam energy deposition. The hydrogen recombination
rates are determined by the Maxwellian velocity distribution of recombining
ambient electrons, which are rapidly equilibrated as the beam kinetic energy
is transferred to them (Sect. 9). Non-thermal ionization rates of hydrogen by
the beam electrons (and subsequent nonthermal recombination) [340, 626, 3,
211] are much less than the thermal rates in these atmospheres at all times
except for the first small fraction of a second. Optical and NUV spectral
calculations with dominant non-thermal ionization rates [795] do not seem to
predict the observed blue continuum properties and Balmer jump strengths
(e.g., Fig. 10).

The formation of the optical and NUV continuum radiation in a F13 electron
beam model is illustrated in Fig. 21. The cumulative contribution functions
are shown for $\lambda=3550$ Å, $4300$ Å, and $6690$ Å. The detailed analyses
of emissivity sources and optical depths explain the emergent spectrum as
Balmer and Paschen recombination radiation from optically thick layers [see
also 416, 407, for further analyses of these models]. The $\lambda=3550$ Å and
$6690$ Å radiative fluxes are very optically thick and escape from the top of
the condensation; the blue ($\sim 4300$ Å) continuum flux is less optically
thick and emerges from the condensation and deeper layers that have a lower
electron density of $n_{e}\approx 10^{15}$ cm-3 and a lower gas temperature.
The excitation by the highest energy electrons in the beam, and to a lesser
extent, Balmer continuum backwarming, heat these deep layers and lead to an
emergent spectrum that departs significantly from an optically thin hydrogen
recombination spectrum from a homogeneous slab. [410] reported self-similar
patterns in chromospheric condensation models (using a basic HI373737Human
Intelligence pattern recognition algorithm) and parameterized the RHD
atmospheres at especially interesting times in their evolution. The two main
parameters are a reference column mass ($m_{\rm{ref}}$) and a reference
temperature ($T_{\rm{ref}}$), which are indicated in the cartoon in Fig.
2(right). These parameterized models facilitate LTE estimates of the Balmer
jump and optical continuum ratios over a large range of condensation column
masses and chromospheric temperature gradients. A prediction from [410] is
shown in Fig. 12 for $T_{\rm{ref}}=11,000$ K and a range of $m_{\rm{ref}}$.
The chromospheric condensation models can account for the trend over most of
the impulsive-phase, flare-only Balmer jump ratios and optical continuum flux
ratios in spectral observations.

Figure 21:  RADYN model of the low flare atmosphere after 2.2 s of heating by
a high-flux, F13 ($10^{13}$ erg cm-2 s-1), double power-law electron beam
[415]. (top panel) The cumulative of the contribution function to the emergent
intensity, $I(\lambda=3550\rm{\AA})$, shows that the chromospheric
condensation becomes optically thick in the Balmer continuum. The emergent
blue-optical ($\lambda=4300$ Å) continuum intensity has larger contribution
from the layers below the condensation, which is the large increase in gas
density in the bottom panel. (Bottom panel) Several trajectories of nonthermal
electrons are shown on top of the gas dynamic quantities (the initial kinetic
energies are indicated, and the kinetic energy variations, $E(z)$, are
normalized from 0 to 1 on a linear scale). The initial temperatures and gas
densities are shown as black dashed lines. A movie is available online.

There has been follow-up work by Kowalski et al. on the modeling of
chromospheric condensations to further evaluate their role as the origin of
cooler line and continuum radiation in flares. Here, we summarize some of the
challenges and successes with the hypothesis. First, these model atmospheres
evolve very quickly, and the emergent flux spectra around the Balmer jump
rapidly change on short timescales. Thus, an exposure-time averaged model
spectrum does not necessarily exhibit properties that are consistent with the
observational constraints. An average condensation model that produces a
continuum spectrum that resembles observations requires a very hard,
$\delta\approx 3$, electron beam distribution [e.g., 416, 418]. In other
cases, an average model spectrum has satisfactorily produced the broadband
continuum shape and Balmer jump strength in two moderate-sized dMe flares with
HST/COS observations [420], and success has been achieved in an even broader-
wavelength comparisons to a superflare event [571]. There are also important
theoretical issues with the huge current density and charge displacement due
to an electron beam with a flux of $10^{13}$ erg cm-2 s-1 that have not yet
been addressed.

A further severe challenge has been encountered in dense chromospheric
condensation modeling of dMe flares, as detailed in [418] and [408]. Updated
hydrogen emission line broadening treatments (Sect. 10) reveal that the
optical depths and the ambient (thermal) charge densities in the model
condensations are far too large to be consistent with the symmetric broadening
in archival dMe flare observations. [527] modeled the broad H$\alpha$ emission
lines in a $E>10^{33}$ erg, $\Delta g\approx-1.35$ mag superflare from AD Leo.
They used updated hydrogen broadening profiles in RADYN simulations and
concluded that lower electron beam fluxes with hard power-law distributions
produce optical depths and electron densities in the low atmosphere that
explain the broadening. The models with energy transport through only thermal
conduction were largely discrepant (Fig. 22).

Figure 22:  (Top) Model spectra of H$\alpha$ that were calculated with the
RADYN code and analyzed in [527]. The larger electron beam heating rates
reproduce the broad, symmetric Balmer lines in an AD Leo flare. Figure
reproduced with permission from [527]. (Bottom) Echelle resolving power,
flare-only spectra of the H$\alpha$ line during the decay of a large flare on
YZ CMi. This figure demonstrates the phenomenological Voigt profile $+$
Gaussian modeling that is used to isolate and characterize the asymmetries in
the wings of the Balmer lines. In this spectrum, a Gaussian model is fit to
the broad blue asymmetry; reproduced from Fig. 42(a) in Notsu et al [552],
with permission.

Similarly to [527], the approaches in [418], [408], and Brasseur et al [88]
employ semi-empirical forward modeling with the RADYN and RH codes. These
studies use large low-energy cutoffs ($E_{c}\geq 85$ keV), hard power-law
indices ($\delta=3$), and high energy fluxes ($10^{12}-10^{13}$ erg cm-2 s-1)
to heat the low chromosphere at $\log_{10}m_{c}/[\rm{g\ cm}^{-2}]\gtrsim-2$ to
gas temperatures hotter than $10,000$ K [see also 409]. This approach results
in optically thick continuum radiation at NUV and optical wavelengths and
broad Balmer wings in the early impulsive and peak phase. However, an
additional model component that is attributed to lateral spatial heterogeneity
in the flare heating is required to simultaneously account for the narrow and
broad components of the hydrogen Balmer emission line fluxes. The slow-rising
Ca II K emission line flux is not yet explained with this approach. Large,
low-energy electron beam cutoffs or some other energy transport mechanism
[e.g., 401] that similarly results in a large fraction of the beam heating
over $\log_{10}m_{c}=-2.2$ to $-1.2$ [409] rather than in the upper
chromosphere ($\log_{10}m_{c}\lesssim-3$) seems to be a promising avenue for
explaining some of the more challenging M dwarf flare observations. For
example, a model with $2\times 10^{12}$ erg cm-2 s-1 and $E_{c}=500$ keV
(i.e., a high-flux, fully relativistic electron beam) adequately reproduces
the spectra of the secondary flare events in the decay phase of the YZ CMi
megaflare [412, 414]. These events produced an “A star on an M star”, such
that the flare-only spectra resembled the main sequence A0 V star, Vega. The
Balmer lines and Balmer continuum are “in absorption”. The colors from two
secondary events are shown in Fig. 12 as star symbols. Alternative
explanations are summarized in [16]. [414] claimed that some flares show a
narrowing of the hydrogen Balmer lines in their impulsive-phase peaks that is
consistent with a spatially unresolved component on the star that has the
hydrogen Balmer lines in absorption.

Relatively little attention has been given to the coronal and late phase
predictions in the recent RADYN models of stellar flares. The RADYN models are
generally viewed as simulations of the short, pulsed beam-heating events that
sequentially light up along a two ribbon flare arcade, as in typical solar
flare geometries (see Sect. 12). Ostensibly, the very hot $T\gtrsim 50$ MK
coronal temperatures produced within a short time are expected to eventually
cool down and shine in emission lines that probe cooler temperatures [295,
28]. However, it is imprudent to speculate further on these issues without
X-ray spectral synthesis because of the vast range of electron densities over
the hot temperatures in the model corona and transition region.

Figure 23:  (_Left_) Hydrodynamic modeling of a coronal explosion on CN Leo.
The panel shows temperature, number density, and the emission measure weighted
cooling function. The solid lines show the evolution at $t=0,2,4,8,10$ s when
a flare heating function was applied, while the dashed lines correspond to
$t=13,16$, and $19$ s. (_Right_) Modeled X-ray light curve over the first 20 s
of the model compared to XMM-Newton data of the early phase of the giant X-ray
flare. Figures reproduced from [661] with permission.

[661] used the Palermo-Harvard adaptive regridding code [582, 69] to model the
hydrodynamic response of the giant X-ray CN Leo flare in Fig. 11(a) over the
first 20 s. They assumed a Gaussian heating function centered in the low
corona and extended into the chromosphere. The atmospheric evolution is shown
in Fig. 23. The flare transition region moves downward, and there is an
increase in density that resembles a chromospheric condensation, but low-
temperature spectral synthesis was not the goal of these modeling efforts [see
242, and Fig. 20 here]. They successfully reproduced the X-ray luminosity
evolution by multiplying the model by an area of $10^{19}$ cm2. The model
X-ray luminosity is shown in the right panel of Fig. 23. They explain the
early phase X-ray light curve as a result of heating within the bottom parts
of a flaring loop during the first 10 s, followed by a $v\sim 1000$ km s-1
evaporation upflow that dominates the coronal emission after that.
Compressional heating at the apex occurs later on as well (see their dashed
temperature profiles) and persists after the applied flare heating function is
turned off at 10 s. The Palermo-Harvard code has modeled longer heating
durations to simulate X-ray flares from the M dwarf Proxima Centauri [608] and
the active giant HR 9024 as well [709, 20]. [608] model the second X-ray flare
in the decay of the event in Fig. 18 here as the first in a series of loops.
They develop a two-phase heating model whereby the total energy input consists
of energy deposition near the loop apex and near the low-altitude coronal
footpoints.

##  9 Atmosphere Modeling of Stellar Flares: Overview of Physics and Methods

###  9.1 An Introduction to Radiative-Hydrodynamic Flare Modeling

Radiative-hydrodynamic flare modeling consists of solving the coupled, non-
linear conservation equations for mass, momentum, energy, charge, radiative
transfer, level populations, and grid motion. The atmosphere begins in a
steady-state, and flare heating is introduced through a term in the energy
equation, $Q_{\rm{flare}}$. In this section, I follow the general framework
and modeling assumptions employed in the RADYN flare code, which has unique
capabilities for stellar flare applications (e.g., radiative transfer and non-
solar gravity).

A conservation/continuity equation for a volumetric quantity $w$ ($w$hatever
cm-3), which is a function of space and time, can be written (in laboratory
frame coordinates) as

$\frac{\partial w}{\partial t}+\nabla\cdot(w\mathbf{v})=0$ (14)

in three spatial dimensions (if the time-derivative is 0, then the equation
leads to the steady-state solution), where $\mathbf{v}$ is the macroscopic gas
velocity vector, $w\mathbf{v}$ is the flux383838More precisely, one means the
flux per unit area, or flux density (current density). of the quantity in $w$,
and we have excluded a diffusion term, $-D\nabla w$ [765], in the parentheses.
The second term is the advection term which is a gradient of the flux (while
in more than one spatial dimension, it is the divergence of the flux). From
here on, we discuss the 1D form of Eq. 14 and consider the spatial coordinate
$z$ (“height”) along a hypothetical magnetic loop extending from below the
photosphere to the loop apex in the corona. At the loop apex, the $z$
direction is parallel to the stellar surface, and the surface gravity in the
momentum equation decreases accordingly. In 1D flare modeling, one assumes
that flows and energy transfer are aligned with the magnetic field, meaning
that the magnetic pressure is assumed to be much larger than the gas pressure
everywhere at all times (which may be appropriate for the center of a strong
magnetic flux tube). Note also that the typical 1D form of the continuity
equation in Eq. 14 usually assumes a constant flux tube cross-sectional area;
for a non-variable magnetic field along the flow direction, the adjustments
that must be made to account for funneling effects are written in [203] and
[612].

The conservation equations are solved simultaneously with an adaptive grid
equation [175] using a semi-implicit numerical scheme that allows for
timesteps that are much longer than the Courant step, which is very small in
the transition region during flares. Implicit radiation-hydrodynamics
essentially consists of linearization and iteration (i.e., multi-dimensional
Newton-Raphson) of the conservation equations with a five-point stencil
monotonic upstream scheme [735, 736] to ensure numerical stability of advected
quantities around steep gradients, such as shocks. The spatial derivatives are
written as finite differences and are centered in time between the current and
next time (semi-implicitly) during the linearization of the Taylor-expansion
about the iterative guesses. We do not discuss the details here and instead
direct the interested reader to [394], [491], [112], [1] for general
overviews393939The advantages of Newton-Raphson methods of solution are also
leveraged in stellar evolution calculations, such as with the MESA code
[580].. We especially recommend the lucid and comprehensive lecture on the
topic in [174].

The Eulerian (fixed space) forms of the conservation equations that are solved
in RADYN flare modeling are written in [10]. Here, I highlight the
conservation of internal energy density, $e$, which is the first law of
thermodynamics (i.e., $dU+PdV=Q$) connecting several of the critical
ingredients in flare RHD:

$\frac{\partial e}{\partial t}+\frac{\partial(ev_{z})}{\partial
z}+(P+q)\frac{\partial v_{z}}{\partial
z}=-\frac{\partial\left(F_{c}+F_{r}\right)}{\partial
z}+Q_{\rm{flare}}+\sum_{n}Q_{n}$ (15)

where $P$ is the thermal gas pressure from all species $s$
($P=k_{B}T\sum_{s}n_{s}$) and enters into the $PdV$ work term, $v_{z}$ is the
gas velocity along the loop axis ($z$ direction), $q$ is a viscous stress term
added to aid in numerical stability, $F_{c}$ is the thermal conduction flux
(temperature diffusion), $F_{r}$ is the radiative flux, and $Q_{n}$ can be any
number of external heating terms, which can be flare-related (e.g., beam-
generated electric fields) or non-radiative sources to keep the photosphere
and corona hot in the pre-flare equilibrium state. The form of thermal
conduction is taken from [685] and [224]. The internal energy density is the
sum of thermal and excitation energy densities:

$e=\frac{3}{2}k_{B}T\sum_{s}n_{s}+\sum_{s,i}n_{s,i}\chi_{s,i}$ (16)

where $n_{s,i}$ is the population level density of species $s$ and level $i$
and $\chi$ represents atomic excitation energies ($\chi=0$ for the ground
state). Equation 15 is derived by starting from the continuity equation
(setting Eq. 14 equal to all $Q$ terms, instead of 0) for the energy density,
$\varepsilon$, which is the usual [e.g., 712]:

$\varepsilon=\rho\left(\frac{1}{2}v_{z}^{2}+e^{\prime}+\Phi\right)$ (17)

where $e^{\prime}$ is the internal energy per unit mass ($e=\rho e^{\prime}$),
$\Phi$ is the gravitational potential energy, $gz$. Magnetic energy density is
neglected. The energy flux is

$\mathfrak{F}=\rho\left(\frac{1}{2}v_{z}^{2}+h+\Phi\right)v_{z}$ (18)

where $h$ is the enthalpy per unit mass (“specific” enthalpy;
$h=e^{\prime}+\frac{P}{\rho}$). After substituting the continuity equations
for mass density and linear momentum density, some algebra and reduction, the
form of Equation 15 is readily obtained404040J. Allred, priv. communication
2017..

The equation of radiative transfer is solved simultaneously with the other
equations through complete linearization and a low-order Feautrier method that
is stable in the presence of steep gradients and shocks (M. Carlsson, priv.
communication 2022). The monochromatic specific radiative intensity,
$I_{\nu}(z,\nu,\mu)$, at frequency $\nu$, at an angle $\theta=\cos^{-1}\mu$
from the atmospheric normal, at a height $z$, is assumed to be in steady-state
($\frac{\partial I_{\nu}(z,\nu,\mu)/c}{\partial t}=0$ where $c$ is the speed
of light; _cf._ Section 11.2 of Hubeny and Mihalas 337) within each
hydrodynamic time-step. The solution to the intensity at each depth point is
integrated over angle using a Gaussian quadrature numerical weighting [123],
and then it is integrated over frequency to give the net radiative flux. The
gradient in the energy equation (Eq. 15) thus gives the radiative heating or
cooling rate at each time-step in the flare simulation. Due to this coupling,
the time-history of the atomic level populations, ionization, rates, and
radiative transfer are critical in the self-consistent energy balance,
pressure gradients, and hydrodynamics in a flare simulation.

###  9.2 The Flare Heating Term

The flare heating term, $Q_{\rm{flare}}$, in the energy equation (Eq. 15) may
be a sum from a large number of possible sources: nonthermal electrons,
nonthermal protons, Joule heating, shocks in the reconnection process, and
Alfven waves. Here, we assume the major source of heating in stellar flares is
through nonthermal electrons, as is widely accepted in and borrowed from the
solar flare analogy. The radio gyrosynchrotron radiation from mildly
relativistic electrons in stellar flares also supports this approach (Sect.
11). The following is a synopsis of several relevant points from the
collisional energy loss theory of [199] and a unified model calculation
presented in [10]. Additional formalism about the theory of energy loss on
ionized species and neutrals is contained in [723], [122], [200], [211],
[296], [517], [688]. We consider particle velocities with magnitude, $v$, and
three directional components (whereas in the previous section, one-dimensional
bulk velocities of the ambient/thermalized gas were denoted as $v_{z}$).

A nonthermal (beam) electron penetrating a plasma deposits heat thus affecting
the internal energy of the plasma by increasing the temperature, excitation,
and ionization of the plasma constituents. A beam particle is assumed to begin
with a kinetic energy that is much greater than the average kinetic energy of
the free electrons in the background gas or plasma (thus, the background
plasma is a “cold target”); the beam particle is then followed until it
decreases to an energy that is close ($\sim 2.5k_{B}T$) to the energy of the
thermal distribution, at which point it joins the background/thermal pool of
particles. Energy loss from the beam electron occurs through Coulomb
collisions, and the interactions with all particles in the plasma require an
integral that, when evaluated (with somewhat artificially imposed distance
limits), is known as a Coulomb logarithm [e.g., 61]. The Coulomb logarithms
give a sense of the relative importance of distant to nearby interactions with
the beam particle. A beam particle loses energy to ambient (thermal) electrons
and protons, with an associated Coulomb logarithm414141Most sources refer to
the Coulomb logarithm as $\lambda=\ln\Lambda$, but we follow the abbreviated
notation in [199]. $\Lambda$. The integral that describes collisions with
ambient neutral atoms also has an associated Coulomb logarithm,
$\Lambda^{\prime}$, for each neutral species. For a typical beam electron,
values for $\Lambda$ and $\Lambda^{\prime}$ are $\approx 20$ and 8.5,
respectively. Including the ionization fraction of a hydrogen gas, $x$, the
fraction of energy loss that goes into the neutrals is [626]
$\frac{1-x}{x\Lambda+(1-x)\Lambda^{\prime}}$. For a 50% ionized gas, the
energy loss is about a factor of three more efficient on the ionized component
(i.e., the ambient thermal electrons).

The pitch angle of the beam particle changes through scattering, which can be
elastic (for collisions with ambient protons, electrons, and neutrals) or
inelastic (for collisions with neutrals). The momentum changes in the
direction parallel to the trajectory of the beam particle must include an
additional Coulomb logarithm, $\Lambda^{\prime\prime}$, for elastic scattering
off of neutrals (e.g., collisions with bound electrons that do not undergo
atomic transitions). Distribution-averaged approximations to the Coulomb
logarithms are given in [199] and [296]; [10, 11] gives the energy-dependent,
relativistic formulae. Then, the total differential kinetic energy change per
unit path length [keV km-1] of an electron beam particle interacting with
hydrogen and ambient electrons in a partially ionized gas is, following from
Eq. 24a of [199],

$\displaystyle\Big{(}\frac{dE}{dz}\Big{)}_{\rm{hyd}}$
$\displaystyle=\Big{(}\frac{dE}{dz}\Big{)}_{\rm{hyd,I}}+\Big{(}\frac{dE}{dz}\Big{)}_{\rm{hyd,N}}$
(19) $\displaystyle=-10^{5}\frac{2\pi e^{4}}{\mu E\times
a^{2}}\Big{(}x\Lambda_{ee}+(1-x)\Lambda^{\prime}\Big{)}\gamma n_{\rm{hyd}}$
(20)

which is the sum of the ionized (I) and neutral (N) components of energy loss
in a hydrogen gas. We ignore the $\sim 10^{3}$ less energy loss on ambient
protons, $E$ is the kinetic energy of the electron beam particle in keV,
$a=1.602\times 10^{-9}$ keV erg-1, $\mu$ specifies its pitch angle ($\mu=1$
for beamed along hypothetical magnetic field), $x$ is the ionization fraction
of the hydrogen target, $n_{\rm{hyd}}$ is the number density of neutral
hydrogen atoms plus ambient protons in the target, $\gamma$ is the
relativistic gamma factor of the beam electron, and all of the rest of the
units are in _cgs_. As an example, if $x=0.75$, $n_{\rm{hyd}}=1.7\times
10^{13}$ cm-3, $T\approx 10^{4}$ K, which are representative conditions at the
top of a model M dwarf chromosphere, then $\frac{dE}{dx}=-0.1$ keV km-1 for
$E=40$ keV. Additional energy losses on helium and other gas constituents can
be included.

Figure 24:  A pre-flare and flare model (from the RADYN code) of the
temperature distribution in the low atmosphere compared to Coulomb stopping
depths of beam particles. For reference the (pre-flare) photosphere and upper
photosphere corresponds to $\log_{10}m_{c}>1$ on this plot ($m_{c}$ is the
column mass density in units of g cm-2). The collisional stopping depths of
mono-energetic beam particles are shown on the right using Equation 19. The
kinetic energy decrease of an electron with injected energy of 50 keV is shown
for each state of the atmosphere; these two curves are on a linear scale and
normalized to 1.0 at the top of the loop (where they are injected) outside the
range of the plot on the right. An electron gets stopped at higher altitudes
when the chromosphere has experienced a large amount of hydrogen ionization;
as the atmosphere further evolves, the lower energy electrons in the beam can
be stopped in the evaporation flows (_cf._ Fig. 21), and all the rest of the
beam particles except for the highest energy electrons are stopped in the
chromospheric condensation (Fig. 2(right)). Note, these “cold-target” stopping
depth approximations for the protons are not accurate significantly below
$\lesssim 1$ MeV at which the warm-target transport treatment becomes
essential [11].

Assuming a fully ionized, uniform hydrogen plasma, Eq. 19 can be integrated
from the electron’s initial energy $E=E_{0}$ to $E=0$ to find the stopping
length of an electron, $L_{\rm{stop}}\propto E_{0}^{2}n_{e}^{-1}$, which
defines an effective absorption cross section, $\sigma$, according to $\sigma
n_{e}L_{\rm{stop}}=1$. Higher energy particles penetrate deeper because they
have more energy to lose and because they spend less time losing energy per
collision. The energy loss rate per electron is,
$\frac{dE}{dz}\frac{dz}{dt}=dE/dt=-2\pi e^{4}n_{e}v\Lambda/E$ which can be
integrated over a uniform plasma to give a timescale for energy loss,
$\tau\propto E^{1.5}/n_{e}$ at non-relativistic energies. Note, this is the
same proportionality in the estimate for scattering timescales out of a
magnetic trap (Sect. 3).

The heating from a power-law distribution of electrons and protons is
analytically solved in [199] following [97] and [450] in an integral from,
$Q_{\rm{beam}}(m_{c})$, which is a function of column mass over an atmosphere
with arbitrary, constant ionization fraction. [296] extended the formulae to
atmospheres with non-uniform ionization fractions. The analytic formulae for
$Q_{\rm{beam}}$ were used in RADYN simulations of stellar flares up through
the models in [415, 416].

An alternative method is to numerically solve the distribution function,
$f(z,\mu,E,t)$, of the nonthermal electrons [440, 494], and evaluate the
spatial gradient/divergence. This is the Fokker-Planck solution [637], which
tracks the phase space evolution due to systematic decrease (drag) and
diffusion of the kinetic energies, $E$, and of the cosines of the beam pitch
angles, $\mu$. Generally speaking, the Fokker-Planck solution is required when
some or all of the following are important: time dependence of the beam
propagation, diffusion of pitch angle and energy, external forces (e.g.,
magnetic mirror, synchrotron losses), and beam feedback effects such as return
currents and plasma waves [293, 490, 794, 10]. We refer the reader to the
comprehensive description of the collisional Boltzmann transport equation in
[723] and the application of the Fokker-Planck solution to flare modeling in
[11]. Other helpful references are [712] and [83]. In current RHD simulations,
numerical limitations necessitate a steady-state solution ($\frac{\partial
f}{\partial t}=0$) at each hydrodynamic time-step, thus giving the flare
heating rate as:

$Q_{\rm{flare}}(z)=Q_{\rm{beam}}(z)=\frac{d}{dz}\int_{\mu}\int_{E}\mu
vEf(z,\mu,E)dEd\mu$ (21)

which is the spatial gradient of the beam flux at height $z$ in the
atmosphere.

##  10 Models of Chromospheric & Transition Region Flare Emission Line
Broadening (& Asymmetries)

The emission line properties (strength and shapes) during flares encode much
of the evolution of the atmospheric physics, which may be quite dissimilar to
the density, optical depth, and velocity regimes in the atmosphere before a
flare424242 _cf._ the vast differences in terrestrial atmospheric states and
radiative transfer in a blue-emitting lightning channel and the clear blue
sky.. The interpretation of flare emission lines, especially in the
chromosphere, in general requires forward modeling with non-equilibrium
physics. Thus, there is still much we do not know about the origin of certain
well-observed changes in the spectral line profile shapes during flares. The
efforts to make headway in this area have been stymied because, further, not
all flares show similar behavior in the same emission line (_cf._ in solar
flares, H$\alpha$ shows a rather wide range of empirical properties; Canfield
et al 110). Detailed line shape measurements are further obfuscated by the
heterogeneity of the flare source and long exposure times, which are generally
required of echelle observations.

This section is divided into two parts. First, we categorize the line
broadening mechanisms that are thought to be important in the interpretation
of transition-region and chromospheric line widths and asymmetries during
stellar flares (Sect. 10.1). Then, in Sect. 10.2, we summarize the symmetric
pressure broadening of hydrogen lines in more detail. Observations and
modeling are reviewed together in this section.

###  10.1 An Overview of Broadening Sources

There are many types of possible sources of chromospheric and transition
region emission line broadening in stellar flare spectral observations.
Microscopic broadening occurs over scales that are much smaller than the
optical depth, and each process is generally assumed to be statistically
independent. Thus, the total microscopic broadening enters the line opacity
after a convolution of all the respective probability density functions.
Microscopic broadening is largely a source of symmetric broadening (to first
order). Macroscopic broadening sources can affect the line opacity though the
mean intensity, $J_{\nu}$, but this term typically refers to the broadening
from spatially unresolved, inhomogeneous sources in the plane transverse to
the observer.

Here, we summarize the most important sources of micro- and macroscopic
broadening that are relevant to stellar flare chromospheric, transition
region, and photospheric lines. We exclude rotational and Zeeman broadening,
which affect the broadening of lines in active stars in their quiescent
states. Note that microscopic broadening sources are modeled with an ambient,
thermal distribution of perturbers even when they are considered separately
from thermal Doppler broadening. In Fig. 25 we show several profiles for the
hydrogen Balmer $\gamma$ and $\alpha$ lines, which complement the descriptions
throughout this section. For each type of broadening, we state the probability
distribution that most closely describes the perturbation magnitudes.

Figure 25: Optically thin spectral line profiles from the TB09$+$HM88 pressure
broadening theory for Balmer H$\gamma$ (top two panels) and H$\alpha$ (bottom
two panels) at $n_{e}\approx 10^{15}$ cm-3. These are identical to the Vidal
et al [746] (VCS73) calculations at this density, and they do not include
thermal Doppler convolution in this figure. The 2nd and 4th panels show the
slopes, $\eta$, calculated for each of the profiles in 1st and 3rd panels.
Lorentzian profiles with scaled damping widths are overplotted for comparison.
The slope of the H$\alpha$ pressure broadening profile is between the
Holstmark and Lorentzian because of the higher importance of collisional
broadening from electrons for this transition. The higher order lines,
including H$\gamma$, are more obviously dominated by the quasistatic
(Holtsmark-distributed) perturbations [_cf._ Fig 3 of 745], are less optically
thick, are closer to LTE because of less scattering in this line, and are less
affected by thermal broadening (the FWHM of H$\alpha$ is affected by thermal
Doppler broadening even at such high densities; $1/e$ Doppler widths are
indicated by horizontal red lines). The “self-broadening” of hydrogen is
relatively less in H$\gamma$ and in the higher order lines, but it is
important when the atmosphere is largely neutral [53]. The Balmer H$\gamma$
line is also less affected by thermal, $T=10^{4}$ K ion-dynamic effects [692]:
we show an area-normalized (rather than peak-normalized) profile from [693],
which includes additional Lorentzian damping near line-center from non-
quasistatic components of the perturbing ionic field. Also shown are emergent
spectral intensities (dashed-dotted black), with Doppler broadening included,
from slabs with a physical depth of $dl=500$ m to demonstrate the large
differences that result from curve-of-growth-like optical depth effects.

  * •

H I$+$p+/e-: “Stark broadening”, or electric pressure broadening of hydrogen,
hydrogen-like ions, or Rydberg atoms by ambient charged particles (electrons,
protons, ions) in the flare chromosphere, which is partially ionized by
definition. The perturbations are caused by both electron collisional
“damping” with a Lorentzian probability distribution and quasi-static
splitting of energy levels from ambient positive charged protons (and ions)
with a Holtsmark-like434343More accurately, this is a Hooper distribution,
which accounts for plasma correlations and Debye screening [323, 532]. A
Hooper distribution approaches a Holtsmark distribution at large field
strengths, which translate to the far wings of a spectral line [322].
probability distribution. The splitting of energy levels is due to the first-
order/linear Stark effect which is discussed in more detail in Appendix D.
_Optically thin emission line profiles from a homogeneous slab exhibit a
Holtsmark power-law slope (
$\log_{10}I\propto\log_{10}\mid\Delta\lambda\mid^{-5/2}$) in the far wings_,
if non-ideal effects are not important.

  * •

LEPT: [556] discuss the role of broadening of hydrogen lines due to
chromospheric plasma turbulence. Power-law electron beams drive a cospatial
drifting return current, whose speed may exceed the threshold for the ion-
acoustic / current instability and therefore lead to the development of ion-
acoustic turbulence in the chromosphere (E. Oks, priv. communication 2022).
The low-frequency electrostatic plasma turbulence (LEPT) is broadband, and it
is quasi-static compared to radiative timescales. LEPT is treated as a
statistically independent broadening process to the ionic microfield
broadening and is therefore convolved with the Holtsmark or Hooper
distribution of field strengths that enters into the nominal hydrogen line
broadening (above). [556] analyzed several archival dMe flare observations of
the Balmer series. Without detailed wing shapes, the values of
FWHM/$\lambda_{0}$ from the data are degenerate with either electron densities
of $>10^{15}$ cm-3 or _rms_ LEPT fields of $10-30$ kV cm-1. Lower electron
densities combined with LEPT development in the chromosphere are favored by
considering the Balmer flux decrements or the last resolved Balmer line.
_Profiles dominated by LEPT have an exponential dependence in the far wings,
$I\propto exp(-k\mid\Delta\lambda\mid^{\gamma}),\gamma\approx 1$._

  * •

non-H$+$e-, non-H$+$H: Collisional broadening of non-hydrogenic lines due to
collisions with ambient electrons (the quadratic Stark effect) or with neutral
hydrogen atoms [van der Waals broadening or the more complete ABO theory; 18,
52]. The quadratic (second order) electron impact damping of resonance ion
lines, such as Ca II H and K, causes a small wavelength shift and a broadening
that is typically $10^{3}$ smaller than the linear effect in hydrogen. The
passage of H atoms is relatively fast and so this is a damping. non-H$+$non-H:
Collisions among like atoms or ions results in resonance broadening. This is
an exchange and thus described as a lifetime broadening effect (Bohm 1951).
_Optically thin emission line profiles have a Lorentzian shape (or Voigt if a
Doppler broadening is included)._

  * •

HI$+$HI: Collisional damping of hydrogen lines from transient perturbations
from other hydrogen atoms. van der Waals and resonance forces, collectively
known as self-broadening, drop off rapidly with distance and so this
interaction has a short duration. The neutral fraction of hydrogen is too
small, even at $\sim 50$%, in the flare chromosphere for self-broadening to be
significant in comparison to broadening by charged particles, unlike in the
quiet Sun photosphere [53]. _Optically thin emission line profiles have a
Lorentzian shape (Voigt if thermal Doppler included)._

  * •

Thermal Doppler broadening is important for line opacity near $\lambda_{o}$
but decreases precipitously according to the Gaussian distribution. Note that
one Doppler width is $\approx$ 0.2 Å for H$\gamma$ at $T=10,000$ K. Thus
thermal Doppler broadening cannot explain broad wing emission unless there are
very high temperatures $T\ggg 10,000$ K. _Optically thin emission line
profiles exhibit a Gaussian profile._

  * •

Microturbulence444444A succinct overview of microturbulence and
macroturbulence in reference to quiet-Sun MHD models can be found in the
discussion in Carlsson et al [116]., or nonthermal Gaussian Doppler
broadening. An _ad hoc_ microturbulence parameter is often used to account for
missing broadening sources when all other known broadening sources have been
exhausted. This parameter enters into the opacity (extinction coefficient)
through a convolution of a Gaussian with the thermal Doppler core. (Certain
values of this parameter may propagate through decades of stellar models with
little critical reconsideration of their physical meaning or necessity.) For
recent use in solar flare modeling of Fe II and Mg II lines, see [418] and
[799], respectively. Non-Gaussian microturbulence is discussed often in the
interpretation of solar flare transition region lines [194], but the densities
are rather large and the optical depths are non-negligible even in transition
region lines of quiescent state of active stars [e.g., 485]. [567] report
nonthermal broadening velocities of several tens of km s-1 in transition
region flare lines with a temperature dependence, such that N V ($T>10^{5}$ K)
does not exhibit excess broadening.

Does the microturbulence parameter represent bona-fide gas-dynamic turbulence
(eddies, vortices, random gas motions) in the footpoints of a flare loop
[e.g., 79], and if so, through which fluid instabilities (e.g., Kelvin-
Helmoltz, Rayleigh-Taylor) does the flare turbulence occur? The turbulence
hypothesis can be investigated by considering the energetics inferred through
spectral lines and the emission line flaring area on the star (see Sect. 12).
[303] reported very broad, symmetric wings of Mg II h and k emission that had
a small, bulk blueshift during a NUV flare observed by the Hubble Space
Telescope/STIS. The wings would require electron densities as high as
$10^{17}$ cm-3 if due to collisional broadening. A Gaussian fit to the broad
wings gives the line-of-sight velocity component root-mean-square (rms),
$\sigma_{\rm{los}}$, and the kinetic energy of the chromospheric turbulence is

$KE=\frac{3}{2}\sigma_{\rm{los}}^{2}\rho_{\rm{chrom}}\Delta
z_{\rm{chrom}}A_{\rm{flare}}$ (22)

where $A_{\rm{flare}}$ is the projected area on the star (Section 12) that
produces the emission line flux at Earth. In this flare, the turbulence would
be highly supersonic (200 km s-1) and its kinetic energy would be several
orders of magnitude larger than the radiated NUV energy.

  * •

Opacity broadening, optical depth effects: It has long been known that
simultaneous modeling of electron density and optical depth is required to
interpret the hydrogen lines properly in stellar and solar flares [e.g., 193].
Optical depth effects cause a type of macroscopic broadening in chromospheric
flare lines in the emergent spectral lines from stellar atmospheres. Sometimes
this is loosely referred to as opacity broadening, which is described in
[775], [605], and [294]. Essentially, opacity broadening follows from the
Eddington-Barbier relation for a depth-dependent source function over
optically thick wavelengths of a line454545In laboratory plasma physics,
opacity broadening is also called “self-absorption”, and it is usually
described in the context of a depth-independent source function [e.g., 787,
429]. Stellar atmospheres, and their flaring states, are not homogeneous
slabs, however. . If the source function decreases in the outer layers a
central reversal can form due to an absorbing layer for very optically thick
lines; this is also sometimes called “self-absorption”. An optical depth
broadening effect that is closely related to opacity broadening occurs in the
enhancement of the emergent intensity at wavelengths that correspond to the
transition from a formation over very large optical depths (where there is
opacity broadening, proper) to smaller optical depths in the far wings. This
is similar to the broadening in textbook illustrations of the curve-of-growth
of a spectral line, whereby the Doppler core saturates at the local source
function, and the Lorentzian or Holtsmark wings are very prominently in
emission but are not close to saturation. The core saturation and wing
enhancements of H$\gamma$ and H$\alpha$ are demonstrated in Fig. 25 for a slab
thickness of $dl=500$ m and electron density $n_{e}=10^{15}$ cm-3. An analysis
of these effects in a heterogeneous, non-equilibrium model chromospheric
condensation is presented in Kowalski et al [421]. [527] analyze in detail the
self-absorption and pressure broadening in stellar flare H$\alpha$ emission
lines.

  * •

Red-wing and blue-wing asymmetries cause an increase in the breadth of a
spectral line, especially if the line is spectrally and/or spatially
unresolved. In stellar flares, red- and blue-shifted components of H$\alpha$
have maximum line-of-sight velocity extents of $\approx 100-300$ km s-1 from
the rest wavelength [743]. Reports of line asymmetries in the literature
attribute the asymmetries to a variety of phenomena that involve chromospheric
mass motions: beam-generated impulsive phase chromospheric condensations,
coronal rain (a “coronal condensation”) from post-flare loop plasma draining
and cooling to transition region and chromospheric temperatures, cool filament
eruptions of plasma, absorbing downflows, and chromospheric evaporation of
cool material into magnetic loops [see 784, for an overview]. Many stellar
flare observations, however, with echelle resolving power have rather long
exposure times ($t_{\rm{exp}}\approx 15$ minutes), and some are serendipitous
single-exposure observations from exoplanet radial velocity monitoring.
Additionally, optical depth effects and other microscopic broadening sources
are important for interpretation. Optically thin emission lines can exhibit a
velocity “smearing” in the emergent spectra, simply due to an integration of
the emissivity over the velocity gradient along the line-of-sight, and
optically thick lines are complicated by opacity broadening [527, 784].

There are many reports of line asymmetries in hydrogen, helium I, and metallic
lines throughout the UV, optical, and NIR that have been interpreted with
these phenomena. The spectral observations consist of a bright emission line
component around the rest wavelength and a fainter asymmetric wing component,
or more than one wing component. Most modeling of these lines uses linear
superpositions of Gaussians to characterize the shifts and widths of each
component. [784], [529], and [531] use Voigt/Lorentzian profiles to isolate a
Gaussian fit to a red-shifted components in large stellar flares. The bottom
panel of Fig. 22 demonstrates this method, which effectively isolates a
$\approx-100$ km s-1 blueshifted component in the H$\alpha$ line flux during
the decay phase of a white-light flare [552]. [197] and [241] report broad
H$\alpha$ lines that are blueshifted during the rise and peak spectra during
flares on UV Ceti and CN Leo, respectively. In the event on CN Leo, the
blueshifts were -10 to -20 km s-1 and the Gaussian FWHM values were 4.8 Å; the
event on UV Ceti produced a similar width but a larger blueshift of -70 km/s.
[243] presents detailed line profiles at high-time resolution during flares on
Proxima Centauri with VLT/UVES. [741] and [743] discuss blueshifts observed
across the Balmer series with some temporal resolution. [552] present a large
survey of echelle spectral observations of the H$\alpha$ and H$\beta$ lines in
M dwarf flares; the spectra had relatively short exposure times and were
complemented with simultaneous optical broadband photometry, allowing the
changes in the profile shapes to be characterized in the context of the
impulsive and gradual phases of the white-light response (e.g., Fig.
22(bottom)). [321] report on high-time cadence monitoring of EV Lac and
present persistent blue-wing asymmetries through the entire duration of a
flare [see also 474]. [530] discuss a transient blueshifted H$\alpha$
absorption feature, while [529] analyze broad H$\alpha$ emission lines with
red-shifted components, in flares on the young, solar-type flare star EK Dra,
and Lalitha et al [434] analyze the redshift evolution of a broad H$\alpha$
line component in a flare on the K0Ve star AB Dor. [288] discuss blueshifted
emission components in AT Mic. Of course, the location on the stellar disk is
generally unconstrained in stellar flares, and the inferred velocities are not
corrected for line-of-sight effects for field-confined flows, contrary to what
is possible for solar flare observations [93].

Flux-weighted line-center redshifts in Si IV and C IV during and after the
impulsive phase of several dMe flares were reported in [302]. Chromospheric-
line, red-wing emission in H$\alpha$, Na I D, and He I 5876 Å was described in
Fuhrmeister et al [239]. These observations were tested against scaling
relations (e.g., Eq. 3) from analytic models of impulsive-phase chromospheric
condensations [220]. Late-phase redshifts, on the other hand, are usually
interpreted as condensations that are analogous to the draining of cool
material from the post-flare loops on the Sun. [784] discuss evidence from the
time-evolution of the broadening of the H$\alpha$ line that could favor beam-
generated condensations in the late-phase instead (which we note would be more
in line with the multi-thread RHD modeling of the continuum radiation in the
decay phase of large flares; Sect. 12, Fig. 30). [241] and [371] report on red
wing enhancements in He I 7065 Å and 10830 Å emission lines, respectively, in
two flares. [245] show remarkably broad He I 10830 Å wings during a flare on
the dM3e star EV Lac; correspondingly broad H$\alpha$ profiles were analyzed
in [244] and consist of a spectrally unresolved red wing satellite component.
[465], [237], and [235] report redshifts in other transition region lines,
such as C II, from low-mass stellar flares. [453] and [43] analyze the profile
variations in FUV lines during flares on AU Mic and EK Dra, respectively. The
optical Balmer lines have been decomposed into “separate” broad and narrow
components, which exhibit clearly different Doppler shifts in some flares
[326, 331, 241, 408]. Extreme blueshifted [328] and redshifted [777, 77]
emission has been reported out to several 1000 km s-1 in chromospheric and
transition region lines. [398] discuss blueshifts and redshifts in a large
sample of galactic field stars at larger distances and over a wide range of
spectral types.

  * •

Macroturbulence is a term that refers to a superposition of several Doppler-
shifted profiles that originate from sequentially-ignited, neighboring flare
loop footpoints, as in a solar arcade. The superposition is intended to
emulate spatially unresolved flows, which are termed ‘directed mass motions’
in the literature. The unresolved flows may be preferentially blueshifted,
preferentially redshifted [641], or equally blueshifted and redshifted [177].
To generate symmetric wings, the emergent radiation from downflows (e.g., a
condensation) would have to precisely balance the emergent radiation, and thus
the optical depths, from the upflows (e.g., evaporation) in neighboring loops.
This precarious balance would have to be maintained over long durations, as
for typical exposure times of $60-300$s. _The lineshapes may take on a variety
of symmetric or asymmetric forms that depart from Holtsmark and Lorentzian
wing profiles._

###  10.2 Symmetric Wing Broadening of H Lines

The symmetric (or very nearly symmetric) wing broadening of hydrogen Balmer
lines is a remarkable property of stellar flares. Full widths at 10% maximum
in emission line spectra are as large as $\approx 15-20$ Å [297, 527]. The
source of the symmetric broadening has been debated as a superposition of
directed mass motions (macroturbulence), gas turbulence (as represented by a
microturbulence parameter), and pressure broadening due to the Stark effect
(Sect. 10.1). [177] and [197] argue against the Stark effect in the large dMe
flares that they analyze by comparing to the optically thin wing predictions
of Holtsmark profiles and the [746] profiles for electron densities of
$n_{e}\approx 10^{15}$ cm-3 [see also 330, 331, 287]. The better fits with
multiple Gaussian components in these studies support the hypothesis of some
type464646The solar flare H$\alpha$ profiles in [349] are invoked for
justification, but these spectra are highly asymmetric to the red and are
broadened out to $\Delta\lambda\sim 10$ Å. of supersonic laminar- or
turbulent-emitting flows. The latter corresponds to a microturbulence
parameter of $50$ km s-1 around line center and $150-600$ km s-1 in the wings.
[508] present similar two-component Gaussian fits to broad and narrow
components in H$\alpha$ and H$\beta$ in flares from the K dwarf LQ Hya; the
broad components show relative shifts with respect to the narrow components
and are interpreted in terms of various types of unresolved mass motions.

The alternative hypothesis for the source of the broad wings of hydrogen lines
is pressure (linear Stark) broadening from ambient charged particles, which
imprints unique signatures in the observations. Namely, the pressure
broadening of hydrogen produces larger widths for larger $n_{j}$ within a
series (e.g., within the Balmer, Paschen, or Lyman series) in the absence of
optical depth effects and heterogeneities (here, $j$ refers to the upper level
of a transition). However, the smoking gun signature is generally thought to
be broader hydrogen emission line wings compared to nearby resonance lines of
metallic ions, such as Mg II and Ca II, which experience a much smaller amount
of broadening due to collisions with charged particles [170, 171]. [297] and
[251] present much broader Balmer series lines in comparison to the Ca II K
emission line [see Fig. 4 of 251] in the impulsive phases of two large dMe
flares [see also 434].

We show two remarkable examples of the hydrogen broadening in stellar flare
spectra in Fig. 26. The spectral broadening and emission line shape
differences among Ca II K, H, and Balmer H$\epsilon$ ($n_{j}=7\rightarrow
n_{i}=2$) are clearly evident in these echelle spectral observations. The top
panel shows a serendipitous flare spectrum ($t_{\rm{exp}}=30$ min) from the
inactive dM3.5-4 star Gl 6l 699 (Barnard’s star) that was studied in detail
and modeled with the RADYN and MULTI codes in [579]. The bottom panel displays
a sequence of spectra during a flare from the dM4.5e YZ CMi at the same
times474747We downloaded these data from the PolarBase public spectral archive
[172, 584]. as shown for the H$\gamma$ spectra in Fig. A.1, top panel, of
[743]. The inactive M-star flare shows some qualitative similarities to the
line broadening in the spectrum of the decay phase of the continuum in the YZ
CMi flare. These similarities are remarkable given that the stars have very
different ages ($>7$ Gyr vs. $<500$ Myr, respectively) and quiescent magnetic
activity levels. The impressive FUV flaring activity of Barnard’s star has
been reported recently in [235]. A detailed analysis of the H9 Balmer line was