Stellar flares

Abstract

Magnetic storms on stars manifest as remarkable, randomly occurring changes of the luminosity over durations that are tiny in comparison to the normal evolution of stars. These stellar flares are bursts of electromagnetic radiation from X-ray to radio wavelengths, and they occur on most stars with outer convection zones. They are analogous to the events on the Sun known as solar flares, which impact our everyday life and modern technological society. Stellar flares, however, can attain much greater energies than those on the Sun. Despite this, we think that these phenomena are rather similar in origin to solar flares, which result from a catastrophic conversion of latent magnetic field energy into atmospheric heating within a region that is relatively small in comparison to normal stellar sizes. We review the last several decades of stellar flare research. We summarize multi-wavelength observational results and the associated thermal and nonthermal processes in flaring stellar atmospheres. Static and hydrodynamic models are reviewed with an emphasis on recent progress in radiation-hydrodynamics and the physical diagnostics in flare spectra. Thanks to their effects on the space weather of exoplanetary systems (and thus in our search for life elsewhere in the universe) and their preponderance in Kepler mission data, white-light stellar flares have re-emerged in the last decade as a widely-impactful area of study within astrophysics. Yet, there is still much we do not understand, both empirically and theoretically, about the spectrum of flare radiation, its origin, and its time evolution. We conclude with several big-picture questions that are fundamental in our pursuit toward a greater understanding of these enigmatic stellar phenomena and, by extension, those on the Sun.

1 Introduction

Stars produce energetic bursts of electromagnetic radiation that follows a sudden magnetic energy release into their atmospheres. These electromagnetic bursts are called flares, which occur on a very wide range of timescales, from seconds to days. Stellar flares are observed in all spectral windows, from the high-energy X-rays through the long wavelength radio waves. However, all wavelength regimes do not respond equally in energy and simultaneously in time. Some wavelengths (e.g., microwave) are dominated by nonthermal radiation, while others (optical, low-energy X-rays) are thought to result from primarily thermal radiative processes. The flare energy—both thermal and nonthermal—ultimately originates in the energy that is released from magnetic fields in stellar coronae. These magnetic fields, in turn, originate from turbulent convective energy transport and shear in rotating stellar envelopes, beneath the visible photosphere. The Sun is the best studied flare star due to its proximity to Earth and a fleet of multi-messenger instruments that continuously provide spatially resolved observations.

From Earth, solar and stellar flares are not detectable by the unaided human eye. Stars like the Sun, and other active stars like Algol that are visible to the naked eye, have too much background glare for our eyes to sense the flare brightness increases. In all-sky surveys that are sensitive to the apparent magnitudes of low-mass stars (fainter than a Johnson V-band magnitude of $\sim 9$), however, stellar flares are the most dramatic source of variability. By most dramatic, we mean that flares produce the largest changes in their apparent magnitudes per unit time when compared against nearly all other astrophysical phenomena.¹

In data from the upcoming Vera C. Rubin Observatory’s ten-year Legacy Survey of Space and Time (LSST), stellar flares will constitute a major source of variability (Hawley et al. 2016). Further, many of the largest events that are observed serendipitously (and those from very faint but populous low-mass stars) will be bona-fide transients, whereby the changes in brightness occur from sources that have not been detected in quiescence. Stellar flares are much shorter in duration and much less luminous than extra-galactic transient phenomena, such as supernovae, tidal disruption events, and optical counterparts to neutron star mergers. Most stellar flares can in principle be readily distinguished from these events given enough empirical information, despite some striking spectral similarities (Chang et al. 2020) in the optical spectra of kilonovae (Shappee et al. 2017). Due to their proximity to the solar system, the largest stellar flare events can trigger gamma ray observatories; these events are known as “superflares” in X-rays. Stellar superflares are factors of $\approx 10^2{-}10^4$ more energetic and luminous in comparison to the largest solar flares, which have bolometric energies of $E_{{\rm{bol}}} \approx 3{-}6 \times 10^{32}$ erg (Woods et al. 2004, 2006; Cliver et al. 2022b; Hayakawa et al. 2023a). Such superflares were likely common when the Sun was very young and rotating much more rapidly than today (Maehara et al. 2012). Understanding the physical processes in stellar flares thus provides insight into the heliospheric conditions in the early history of our solar system (Ribas et al. 2005).

It is timely for a review of recent observations and models of stellar flares. The next sunspot cycle maximum is approaching in late 2024 (Upton and Hathaway 2023),² and with it a deluge of flares and eruptions from the Sun. New and upcoming solar observatories, such as the Daniel K. Inouye Solar Observatory (Rimmele et al. 2020), the Expanded Owens Solar Valley Array (Gary et al. 2018), and the Interface Region Imaging Spectrograph (De Pontieu et al. 2014), will provide new avenues for solar-stellar comparisons. The immaculate precision of Kepler (Koch et al. 2010), K2 (Howell et al. 2014), and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) observations has recently facilitated a resurgence in the study of optical broadband stellar flares. These missions provide access to stellar magnetic activity over long time baselines previously not feasible to observe from ground-based campaigns, and over a much larger variety of stars. Models of energy transport, atmospheric response, and emission line broadening have increased in accuracy and sophistication over a large range of heating parameters. As the community looks to next-generation modeling paths, analysis methods, and observational capabilities (such as new space missions and ground-based instruments), synthesizing recent findings and outstanding problems could help to steer these into productive directions. Additionally, stellar flare radiation is now considered an important factor in assessing exoplanet habitability and photochemistry (Shields et al. 2016; Segura 2018), and general interest within the solar and astrophysics communities has grown over the last decade.

1.1 Overview of this review

There have been several reviews on stellar flares prior to ca. 1990 (Pettersen 1989; Byrne 1989; Haisch et al. 1991). The current review supplements these with results from the past three decades. This review features many results from flare studies of low-mass, M-dwarf (dM/MV) stars. Because of their low background glare from non-flaring regions, large contributions to Galactic populations, and inherently high flare rates, the M dwarfs tend to be the most commonly studied flare stars. Flares from active binary systems have also been very well-observed across the electromagnetic spectrum. Results from studies of post main-sequence single stars and young solar-type stars are covered as well. Surprisingly, certain wavelength regimes in M-dwarf flares have been more thoroughly observed than solar flares. Solar flares are not reviewed in any detail here, but a general overview is provided, since the interpretation of and modeling approach to stellar flares are dependent on knowledge of the particle acceleration, magnetic fields, and spatial resolution from the study of the Sun.

The primary purpose of this review is to serve as a compendium of references and to facilitate research in multi-wavelength observations and models of stellar flares. The target audience consists of beginning PhD students or interested scientists in other areas of astrophysics or solar physics. Every current topic in the study of stellar flares is not included here: for example, I leave all results from the TESS mission to the next iteration of this Living Review. Unfortunately, many important results and references cannot be discussed extensively in order to keep this review a reasonable length, while also allowing growth with future addenda. For the same reasons, every important result from the references cited herein regrettably cannot be included at this time. There are several comprehensive introductions to sub-topics in the stellar flare literature; these are indicated throughout instead of being repeated. Broader applications of flare research beyond stellar astrophysics (e.g., exoplanet habitability and exoplanet atmospheric photochemistry) are outside the scope of this review and are not covered. The study of stellar coronal mass ejections and their associated energetic particles is a vast and newly emerging field; these topics are only very briefly alluded to. This review focuses on the electromagnetic response of flaring stellar atmospheres and detailed modeling of the associated physical processes.

This review is organized as follows. We begin with a brief tour of the flare star next door, Proxima Centauri, which is about the same age as the Sun but is otherwise a much different star (Sect. 2). Then we take a detour for a brief introduction to solar flare terminology and phenomenology (Sect. 3). A general overview of all known types of stars that tend to flare is summarized in Sect. 4. The bulk of the review contains a synthesis of stellar flare observations across the electromagnetic spectrum with a focus on the near-ultraviolet (NUV) and optical response, which has not been the subject of most other reviews. Observational studies can be separated into two general topics: flare rates, which are primarily diagnosed through single-band photometry (Sect. 5), and multi-wavelength analyses, which are primarily accomplished through spectroscopy (Sect. 7). The final third of the review summarizes stellar flare modeling approaches (slab, semi-empirical, and radiative-hydrodynamic), beginning in Sect. 8. We go into detail in several areas within stellar flare modeling, focusing on key elements of radiation-hydrodynamics (Sect. 9) and the interpretation of chromospheric spectral line broadening in flares (Sect. 10). A comprehensive analysis of an “ideal” multi-wavelength data set is discussed in Sect. 11, and inferences of the stellar flare geometries are reviewed in Sect. 12. We conclude in Sect. 13 with six questions that we think are currently central in stellar flare research (excluding questions directly related to exoplanets). Throughout, we use the cgs (Gaussian) units system except for wavelengths in Å or in cases that are otherwise noted. In the appendices, we include a visual guide to common filters used in optical studies of stellar flares (Appendix A), some additional observational results from the literature are merged (Appendix B), several common slab modeling practices and assumptions are reviewed (Appendix C), and the basic microphysical processes that symmetrically broaden hydrogen lines in flares is summarized (Appendix D).

2 Our flare star neighbor: Proxima Centauri

The nearest star outside the solar system is Proxima Centauri, which is $\approx$1.3 parsecs, or 268,000 astronomical units (au), from the Earth, and is a low-mass ($\approx 0.12$ M$_{{\rm{Sun}}}$) M dwarf star of spectral type M5.5Ve (dM5.5e) (Hawley et al. 1996). The hydrogen Balmer $\alpha$ line (6562.81 Å) is in emission (e)³ in quiescence outside of detectable flaring events, which occur very frequently. This star is a relatively old ($\approx 5-6$ Gyr) main sequence star, but it is still very magnetically active, consistent with the long magnetic activity lifetimes of low mass stars of its spectral type (West et al. 2008). We describe some of the flare properties of this star from a night of monitoring during an observing run at La Silla Observatory in 2010, while also introducing basic point-source photometry measurements and terminology that are common to almost all analyses of flare star data.

A NUV light curve of Prox Cen (Kowalski et al. 2016) is shown in Fig. 1. The photometry was measured through a custom filter with a full-width-at-half-maximum (FWHM) of 100 Å centered at $\lambda _{{\rm{cen}}} = 3500$ Å. This is a wavelength range over which the star has a very faint luminosity during non-flaring (“quiescent”, denoted by q) intervals. In less than seven hours of continuous 3 s exposures (and 1.5 s exposures at the end of the night), the star produced at least 16 flares at >3$\sigma$ (1$\sigma \approx 6$%) above the mean, which is indicated by a flux enhancement equal to 1.0 in the light curve. In quiescence, the star is very “red” since it has a very cool photospheric temperature ($< 3000$ K), but the flare light causes the star’s total flux at Earth to become much “bluer”. In addition to the time-integrated flare energy (fluence), the peak flux enhancement is sometimes used to characterize the size of a flare. The peak magnitude change of a flare is related to the peak flux enhancement by

$$\begin{aligned} \Delta \text{mag}_{\rm{peak}} =-2.5\ \text{log}_{10}\ \left( C_{\rm{rel}} (t_{\rm{peak}})/C_{\rm{rel},q} \right) = -2.5\ {\text{log}}_{10} \big (I_{\rm{f}}(t_{{\rm{peak}}})+1 \big ) \end{aligned}$$

(1)

where $C_{{\rm{rel}}}$ represents the measured counts per exposure through a bandpass and is calculated relative to the counts in that same exposure from a nearby comparison (non-variable) star or set of stars. The count flux enhancement, $C_{{\rm{rel}}} (t)/C_{{\rm{rel}},q}$, is normalized to a time interval over which the flare star is not clearly varying, usually taken right before a flare. The count flux enhancement has traditionally been denoted as $I_f + 1$ (Gershberg 1972). Because the star is spatially unresolved, the values of $I_f(t) + 1$ during a flare consist of the flux from the non-flaring regions on a flare star in addition to the flaring source; if the flare source has an optical thickness $\tau \gtrsim 1$ at wavelengths in this bandpass, then some or all of the pre-flare flux from the flare area may be diminished at time t (Hawley et al. 1995).

Fig. 1

A light curve of Proxima Centauri in a narrowband (filter FWHM of 100 Å) centered in the near-ultraviolet, U-band wavelength regime. These data are from Kowalski et al. (2016) and can be obtained from the Zenodo repository at https://doi.org/10.5281/zenodo.45878. The horizontal dashed lines indicate the mean flux enhancement (1.0) and $3\sigma$ above the mean (where $\sigma$ is the standard deviation outside of flaring). The inset shows the flare peaking at 2010 May 25 01:31:02 UT, which is referred to as the “IF10” event in Kowalski et al. (2016). The inset shows time relative to peak

None of the flares in Fig. 1 are all that energetic. The flare with the largest peak flux enhancement of five would correspond to an energy of $3\times 10^{29}$ erg in the Johnson U bandpass, which is traditionally the optimal broadband filter ($\lambda _{{\rm{peak}}} \sim 3700$ Å, a full-width-at-half maximum of $\sim 700$ Å; Moffett 1974; Bessell and Murphy 2012)⁴ for stellar flare observations at ground-based observatories. This energy is nearly four times larger than the average U-band flare energy from Proxima Centauri (Walker 1981) and a factor of 100 smaller than the largest event on the Sun that has been observed at similar wavelengths (Neidig et al. 1994). Multi-wavelength scaling relations (Hawley and Pettersen 1991; Kretzschmar 2011; Schrijver et al. 2012; Osten and Wolk 2015; Namekata et al. 2017; Cliver et al. 2022b) facilitate comparison to the standard classification of solar flares according to their peak flux at Earth over the $\lambda = 1{-}8$ Å band, as measured by the X-ray Sensor (XRS) on the Geostationary Operational Environmental Satellite (GOES). The A1.0, B1.0, C1.0, M1.0, and X1.0-class flares correspond to a range of peak X-ray irradiances logarithmically spanning $10^{-8}$ to $10^{-4}$ W m$^{-2}$, while those at $\ge 10^{-3}$ W m$^{-2}$ are sometimes clipped (Hudson et al. 2024) and are given designations beginning with X10.0. The scaling relations imply that the large event on Proxima Centauri would correspond to approximately a GOES 1–8 Å C-class flare if it had occurred on the Sun. However, Prox Cen occasionally erupts with energies that are comparable to, or even much greater than, typical GOES X-class flares from the Sun (Güdel et al. 2002a; Howard et al. 2018). Prox Cen also hosts a near Earth-mass planet at $\sim 0.05$ au (Anglada-Escudé et al. 2016), and unlike our relatively safe distance from solar flares at 1 au, equivalent flare energies from this star would bathe the surrounding planet in 400x the flux of high-energy radiations. Higher-mass, M3-M4 stars are known to emit even higher-energy flares than Prox Cen. For example, the most luminous event known resulted in a remarkably fast $\Delta t \approx 35$ s rise to a peak flux change of $\Delta V = - 5$ mags ($L_{V,\rm{peak}} \approx 1.7 \times 10^{32}$ erg s$^{-1}$) from a young M4$+$M4 binary DG CVn (Caballero-García et al. 2015), for which the estimated values of the U-band energy and GOES class are $4.7 \times 10^{34}$ erg and X600,000, respectively (Osten et al. 2016; Youngblood et al. 2017). Pettersen (2016) describes a remarkable EV Lac flare with a peak magnitude change in the U-band of $-7.2$ mags ($L_{U,\rm{peak}} = 4.6 \times 10^{31}$ erg s$^{-1}$) and a U-band energy of $7.23 \times 10^{33}$ erg. Howard et al. (2019) highlight several flares with extreme amplitudes ($\Delta g^{\prime } \le -3$ mag) and energies ($E \approx 10^{35}{-}10^{36}$ erg) on M1-M4 stars in their optical survey.

The most common qualitative description of a flare light curve is a “FRED”: a fast-rise, exponential-decay. In Fig. 1, the flares exhibit a simple FRED shape at low time-resolution; however, at high-time resolution, much more variation in the temporal morphology is apparent, including two periods of the rise phase (1a, 1b), an extended peak (1c), two intervals of fast decay (2a, 2c), a stall between these two intervals (2b), and a gradual decay phase (3). These features are clear in some other events too (e.g., Fig. 2 of Kowalski et al. 2013, which is reproduced in the top left of Fig. 9 here), though the respective phases may have different durations, relative amplitudes, and integrated energies. High-time resolution observations of stellar flares were actually routine in the 1970s and 1980s using photometers, which have been superseded by high efficiency CCDs, culminating in the unprecedented precision from Kepler, K2, and TESS. The short timescale variations in stellar flares have actually long been recognized. The seminal study of Bopp and Moffett (1973) aptly summarized their findings as follows: “As the time resolution of observations has improved, the great complexity of the flare phenomenon has been revealed. The classical definition of a flare (i.e., rapid rise to maximum followed by a slower quasi-exponential decay) appears to be a gross oversimplification of the complex structures observed”.

3 Solar flares and the standard flare model

Only the very largest solar flare events could be detected in the optical if the Sun was at a similar distance as other stars. Recently, Kepler has provided the precision to detect a signal of 0.0170–0.0270% in optical enhancements that have been observed in Sun-as-a-star data (Woods et al. 2004; Kretzschmar 2011; Moore et al. 2014) and from slowly rotating G dwarfs (Maehara et al. 2012; Notsu et al. 2019; Okamoto et al. 2021). Nonetheless, it is generally thought that flares from other stars originate from the same or similar processes as solar flares. This is supported by the empirical Neupert effect (Neupert 1968), which was first reported in stellar flares in Hawley et al. (1995) and Güdel et al. (1996) (Sect. 7.7). In this section, we briefly review the standard solar flare model paradigm, whose phenomenology (e.g., footpoints, loops) and fundamental physical processes are widely adopted in the interpretation and analyses of stellar flares—either through direct application or through some sort of physical scaling to higher densities, magnetic fields, accelerated particle fluxes, etc... This section synthesizes decades of observations and theoretical work, and it presents our own (rather highly simplified) viewpoint of the entire process. Thus we give a non-exhaustive reference list. For more extensive reviews of solar flares, see Švestka (1976), Hudson (2007), Hudson (2011), Benz (2017), Hudson (2021), the entries within the Space Science Reviews Volume (Dennis et al. 2011), the books by Aschwanden (2005) and Tandberg-Hanssen and Emslie (2009), and the review by Shibata and Magara (2011). For modern, comprehensive reviews of the observations and modeling of solar flares, see Fletcher et al. (2011), Reeves (2022), Kerr (2022), and Kerr (2023).

In the standard flare paradigm, magnetic potential energy is transferred to the atmosphere, which responds by radiating away this energy as the flare. There are generally two categories of solar flares (Thalmann et al. 2019; Kazachenko 2023). Mass and magnetic field are ejected away from the Sun during eruptive flares. In contrast, magnetized plasma is not ejected in confined, or compact, flares or the eruption may be undetectable. The collective process of the mass eruption and the electromagnetic flare is called a “solar eruptive event” (SEE). The left panel of Fig. 2 illustrates the geometry of eruptive magnetic field above compact flare loops, with reconnection of magnetic fields in between (Shibata et al. 1995). This framework underlies most modern generalizations of flare-productive magnetic field topologies (e.g., Fig. 1 of Kazachenko et al. 2022), which we now describe in further detail.

An SEE begins in active regions where there are colliding, non-conjugate⁵ sunspots of opposite polarity (Chintzoglou et al. 2019; Toriumi and Wang 2019; Rempel et al. 2023). The demarcation that separates the bulk of the north and south polarities in an active region is termed the polarity inversion line (PIL). At the PIL, mixed negative and positive magnetic polarities appear as “salt-and-pepper” patterns in magnetograms where there is newly emerging magnetic flux. Ongoing magnetic flux emergence adds to a twisted, cool filament (flux rope) that is parallel to the PIL, and it is thus driven to an instability (Lin et al. 2001). The start of the SEE occurs when this filament becomes unstable and erupts through overarching magnetic field into interplanetary space, developing into a coronal mass ejection (CME). The CME produces a shock that accelerates a large number of solar energetic particles (SEPs; Reames 2021) away from the Sun; these particles may reach Earth within $\approx$ 10 min of the X-ray flare, peak about $\approx$ 6–12 h after, and persist for days. Another class of SEPs are prompt or impulsive-type SEPs, which originate from the flare site. The arrival of the CME itself within $\approx$ 18–50 h of the flare disturbs the terrestrial magnetic field (if the orientation of the CME magnetic field is opposite that of the Earth’s magnetic field upon arrival), induces DC electric fields and currents in the ground, and increases the particle flux into the poles and radiation belts. The ground-induced currents can damage transformers if power is not diverted away from high-load regions.

Returning to the Sun, the solar magnetic field lines that initially arch over a filament current channel (which is relatively low-lying in the corona; Sun et al. 2012; Rempel et al. 2023) are “stretched” during the eruption. Oppositely-directed magnetic field lines pinch together in a “current sheet” in the wake of the erupting filament (a simplified analogy to magnetic tension and retraction is often made to a rubber band building up elastic tension and releasing it by snapping back). The fields undergo magnetic reconnection at many so-called X-points, and the rapid speed of this process is thought to be facilitated by the tearing mode/plasmoid instability and MHD turbulence. The magnetic potential energy in the stretched field line is released as they shorten and retract (Sect. 4 of Longcope et al. 2018).

The fundamental physics of the conversion of magnetic energy into particle kinetic energy is described in the overview by Dahlin (2020). As the retraction of field proceeds to a more relaxed state, an Alfvén-ish-speed outflow (“jet”/“exhaust”) is directed from the reconnection region in the direction toward the stellar surface. Slow-mode MHD (“Petschek”) shocks (e.g., Longcope et al. 2016) and parallel electric fields (Egedal et al. 2015) form in the just-reconnected, bent field lines, enhancing the temperature and density in the reconnection outflow. Thus, the magnetic potential energy initially stored in the contorted magnetic field is converted into kinetic energy of bulk flows, the thermal energy of the plasma, and the kinetic energy of particles out of the ambient/thermal distribution.

Fig. 2

(Left) Illustration of several salient features and processes in the standard model of solar flares, reproduced from Shibata et al. (1995) with permission, copyright by AAS. (Right) Generalization of the chromospheric evaporation and condensation phenomena in radiative-hydrodynamic models of high-flux electron-beam heating at the footpoints of newly-reconnected flare loops (adapted from Kowalski and Allred 2018). Here, we also include an illustration of the radiative backwarming of the photosphere, which is predicted by some 1D models (Allred et al. 2006) of flare chromospheres that are optically thin at Balmer and Paschen continuum wavelengths (see Fisher et al. 2012, for the expected geometry of 3D radiative backwarming of the photosphere)

In large solar flares, a sizable fraction of the magnetic energy that is released seems to go into accelerating ambient particles (electrons, protons, ions) to mildly relativistic energies (Lin et al. 2003; Emslie et al. 2012; Warmuth and Mann 2016, 2020). The nonthermal (“beam”) distribution is a power-law extending from $\sim 10$ keV to tens of MeV for electrons and $\sim 1{-}1000$ MeV for protons. Recent observations have shown that the Sun is a prolific particle accelerator: $\approx 1{-}100$% of the pre-flare coronal electrons are accelerated (White et al. 2003; Kundu et al. 2009; Krucker et al. 2010; Oka et al. 2013; Krucker and Battaglia 2014; Fleishman et al. 2022; Kontar et al. 2023). The inferred nonthermal particle flux density into the lower atmosphere can attain a correspondingly extreme range, $\sim 10^{12}{-}10^{13}$ $\hbox {erg cm}^{-2}\,\hbox {s}^{-1}$, at some locations in solar flares (McClymont and Canfield 1986; Canfield et al. 1991; Wulser et al. 1992; Neidig et al. 1993a; Krucker et al. 2011).

Determining how and where the bulk of nonthermal particles is accelerated in the solar atmosphere is actually a very active research topic (for an overview, see Zharkova et al. 2011). Generally speaking, there are three classes of particle acceleration that may operate in the solar atmosphere (Aschwanden 2005, and see also Knuth and Glesener 2020 and Chen et al. 2020 for succinct overviews of the more modern classifications): first-order Fermi and shock (including shock-drift and diffusive shock), stochastic (“AC”, or resonant interaction with a spectrum of magnetic turbulence), and electric field (“DC”) acceleration. Many theoretical frameworks roughly fall under one of these types in one way or another, and all may operate to some degree over the volume and time-evolution of a flare.

In and around the volume of reconnection and its outflow jets, there have been exciting advances in producing power-laws from multiple first-order Fermi reflections (through curvature-drift motions of particles) among volume-filling plasmoids, which are contracting and coalescing magnetic “islands” (Shibata and Tanuma 2001; Drake et al. 2006; Oka et al. 2010; Karpen et al. 2012; Drake et al. 2013; Dahlin et al. 2014; Guidoni et al. 2016). These models have recently been extended to 3D (Dahlin et al. 2015) and expanded to much larger physical sizes (Drake et al. 2019; Arnold et al. 2021) than particle-in-cell simulations, demonstrating efficient acceleration of high-energy electrons into power-law tails that extend over many orders of magnitude. The Alfvenic-speed reconnection exhaust may develop turbulence, and stochastic acceleration within turbulence (Miller et al. 1996; Petrosian and Liu 2004; Petrosian 2012; Bian et al. 2012) in the presence of Coulomb collisions (Bian et al. 2014) is thought to be one of the possible acceleration mechanisms.⁶ Sub-Driecer (e.g., Holman 1985; Benka and Holman 1992; Alaoui et al. 2021) electric fields have also been investigated in accelerating a small fraction of the ambient particles, and super-Driecer electric fields (accelerating the entire ambient population) that result from the observed decay of magnetic field over a large loop-like volume have also been discussed (Fleishman et al. 2022), possibly related to a type of hybrid particle acceleration involving electric fields and turbulent reconnection (Vlahos and Isliker 2019). Despite the apparent lack of consensus among these impressive modeling and theoretical directions, a commonality is that a significant amount of energy is transferred to the atmosphere through shortening of magnetic field lines (facilitated by reconnection), and particles are efficiently accelerated over a large volume. For reviews of magnetic reconnection and particle acceleration, see also Uzdensky (2007), Cassak et al. (2008), Priest (2014), Cassak et al. (2017), Pontin and Priest (2022), Arber et al. (2015), Lazarian et al. (2020), Nishikawa et al. (2021), Ji et al. (2022).

We should mention that acceleration at a coronal termination shock where the reconnection outflow/exhaust collides with the lower-lying (previously reconnected) magnetic loops could also be important (e.g., Somov and Kosugi 1997), and plasmoid interaction with a fast-mode MHD shock at the collision interface could contribute to the accelerated particle flux into the lower atmosphere (Shibata and Tanuma 2001). Acceleration based on terrestrial aurorae or magnetotail sub-storm processes are also considered (Birn et al. 2017; Haerendel 2018), and some models argue for a mechanism where the bulk of electron acceleration takes place close to (Simnett and Haines 1990; Tsiklauri 2017) or within (Fletcher and Hudson 2008) the chromosphere. Whatever the case or cases may be, acceleration models should attempt to reconcile one way (e.g., Brown et al. 1998, and see also Section 3.6 in Aschwanden et al. 1996a) or another (Cheng et al. 2012) with energy-dependent time delays of hard X-rays (Aschwanden et al. 1995, 1996b, a; Aschwanden 1996; Qiu et al. 2012; Altyntsev et al. 2019; Knuth and Glesener 2020) and the general agreement between inferred path lengths and direct imaging measurements of flare loop sizes. Other critical tests are matching the properties of the white-light continuum radiation (Fletcher et al. 2007), reproducing the energy content in accelerated protons (e.g., $E \approx 7 \times 10^{32}$ erg above a proton kinetic energy of 30 MeV; Murphy et al. 1997), and generating strong nonthermal electron radio signatures at coronal heights (Chen et al. 2020).

In order to consider the rest of the flare process, let us assume that the nonthermal particle beams are injected into a region around the apex of a just-reconnected-and-retracted coronal loop. Initially, the particles are injected in relatively low-lying magnetic fields ($\lesssim 10$ Mm), resulting in the onset of the flare impulsive phase. Over time, larger (and weaker) magnetic fields reconnect and energize particle beams. The particle beams are injected into the newly reconnected loops with a distribution of pitch angles with respect to the magnetic field orientation. The $E \gtrsim 100$ keV electrons in the beams produce gyrosynchrotron radiation as they spiral (with Larmor radii of $\sim$ cm) in the magnetic fields. If the magnetic flux density [G] increases along the beam path (e.g., if the cross-sectional area of the loop correspondingly converges into the lower atmosphere), then the particles with small pitch angles with respect to the magnetic field direction will freely stream (“precipitate”) into the chromosphere. Those with larger initial pitch angles will gradually and adiabatically increase their pitch angle, which is a result of Faraday’s Law of Induction in the frame of the gyrating electrons. These particles can eventually reflect, or “mirror”, at the footpoints (which could occur in the transition region, high chromosphere, low chromosphere, or low corona).⁷ After some time, the mirrored particles too will scatter off enough ambient particles in the corona and precipitate out of the magnetic “trap” into the chromosphere. This is the so-called “trap$+$precipitation” paradigm (Melrose and Brown 1976). The highest energy particles take the longest to scatter and thus remain trapped for a time proportional to $\frac{E^{1.5}}{n_e}$ where E is the kinetic energy of the nonthermal electrons and $n_e$ is the ambient/thermal electron density that the beam heats (under the dilute beam assumption: $n_{{\rm{e}}} \gg n_{{\mathrm{e-beam}}}$). These times are rather short ($< 10$ s) for typical conditions and energies. In solar flares, proton/ion beams can certainly contribute to atmospheric heating (Procházka et al. 2018; Allred et al. 2020), nuclear excitation, a gamma-ray continuum spectrum, the 2.2 MeV neutron-capture Deuterium-formation line, the 511 keV positron-annihilation line, and neutrons that are directly detected at Earth (Murphy et al. 1997; Vilmer et al. 2011). For the sake of brevity, we further consider the effects of only accelerated electrons.

When the electron beam particles stream into the chromosphere, they encounter a wall of dense, partially ionized gas to which they rapidly lose their energy through Coulomb collisions with ambient electrons. At the same time, free-free radiative transitions in collisions with ambient protons occur in the chromosphere, producing hard X-ray ($E \gtrsim 25$ keV), nonthermal bremsstrahlung radiation that exhibits a power-law distribution. Hard X-ray light curves define the impulsive phase in solar flares. The standard model of hard X-ray emissivity that accounts for Coulomb energy loss as the electrons radiate bremsstrahlung is called the collisional (cold) thick target model (CTTM; Brown 1971, for reviews, see Brown et al. 2003; Holman et al. 2011; Kontar et al. 2011). Electron beam power-laws that are inferred from the CTTM and injected into radiative-hydrodynamic simulations are typically characterized by power-law indices of $\delta \approx 5$, a low-energy cutoff of $E_c \approx 15{-}25$ keV, and energy flux densities in the range of $\approx 10^{10}$ to $\approx 10^{11}$ $\hbox {erg cm}^{-2}\,\hbox {s}^{-1}$ (Carlsson et al. 2023). Often, the hard X-ray light curves consist of gradual variations with superimposed short bursts of ${\mathcal {O}}(0.5{-}5)$ s. In some interpretations, these timescales are sensibly related to the processes of prompt and delayed precipitation (described in the previous paragraph) of nonthermal electrons.

The relative displacement of the accelerated electrons from the slower ions in the corona results in a “return current electric field” (e.g., van den Oord 1990; Siversky and Zharkova 2009); this electric field drains energy from the beam and transfers it as a drifting velocity component to the ambient Maxwellian distribution of electrons. The drift is towards the loop apex. In steady state, the flux densities (J; el s$^{-1}$ cm$^2$) of the beam and return current drifting electrons are equal and opposite, thus preventing pulsar-strength, $B \approx 10^9$ G, magnetic fields from forming in the solar atmosphere. The resistivity of the background results in “Joule heating” of the coronal plasma ($\approx \eta e^2 J_{{\rm{beam}}}^2$; Holman 2012; Allred et al. 2020; Alaoui and Holman 2017). For non-dilute beams, the density of the electron beam is comparable to the ambient electron density, and interactions of plasma wave turbulence and parallel electric fields (e.g., Langmuir waves, electrostatic double layers) with the beam are additionally expected.

The nonthermal electron beam kinetic energy is thermalized in the low corona and mid-to-upper chromosphere; see right side of Fig. 2. Initially hydrogen primarily provides the radiative cooling that balances the chromospheric temperature increase. As the chromosphere increases from $T \approx 8000$ K to $T \approx 30{,}000$ K, hydrogen becomes fully ionized, and then helium I and helium II take over the radiative cooling in the range of $T \approx 30{,}000{-}80{,}000$ K. Thus, it is thought that non-equilibrium rates of helium and detailed radiative transfer are important for calculating accurate evolutionary states of the atmosphere. After helium is fully ionized, oxygen, carbon, and neon regulate the cooling at $T \approx 10^5$ K, which is the peak of the optically thin cooling curve (e.g., Cox and Tucker 1969; Rosner et al. 1978; Dere et al. 1997; Dorfi 1998). If enough heat is deposited to fully ionize these elements, then a temperature runaway, or explosion, ensues because the optically thin cooling decreases as a function of temperature up to $T \approx 20$ MK. Thermal conduction transports energy away from the region of maximum beam deposition. Thus, the temperature “bubble” expands, resembling a one-dimensional blast wave (which is more or less confined to the magnetic field direction) in an exponential atmosphere with pervasive beam heating. The increased pressures and their gradients result in a multitude of mass motions. Additionally, an increase in temperature of the very low corona (through thermal heat conduction) can also drive mass motions on top of the beam-generated mass motions. The sources and sinks of the equation for energy conservation are discussed further in Sect. 9.

The impulsively-generated upflows and downflows have been studied numerically as part of the standard solar flare paradigm for many decades (Livshits et al. 1981; Fisher et al. 1985a, b, c; Fisher 1989). They are respectively known as explosive chromospheric evaporation and chromospheric condensation (Fig. 2, right).⁸ The condensations and evaporations manifest as redshifts and blueshifts in cool and hot lines, respectively, with a transition temperature that has been constrained in solar flares (Milligan and Dennis 2009). The high-speed upflows fill the loops with $T \approx 10{-}30$ MK plasma, while the higher density downflows accrue/accrete mass over time and radiatively cool through $T \approx 10{,}000$ K (i.e., this is a non-adiabatic process). The chromosphere is compressed as in a snow plow effect while also experiencing continuous energy deposition directly by the beam. The condensations are a result of shock phenomena because they increase in density up to 10x the ambient density (whereas sound waves would ostensibly smooth out the density perturbations). The closest physical process to the temperature bubble and condensation that we have found (e.g., Kowalski et al. 2022) is that of a thermal wave (Section 6.6, pp. 671–672 of Zel’dovich and Raizer 1967, and see also Atzeni and Meyer-ter Vehn 2004 for similarities to supersonic ablative heating waves in laboratory implosion experiments), which are described as a “second-kind temperature wave” in Livshits et al. (1981).

The representative properties of the condensation and evaporation flows at any time in the evolution follow from a simple momentum balance equation:

$$\begin{aligned} \left( v \rho \Delta z \right) _{{\rm{cond}}} \approx - \left( v \rho \Delta z \right) _{{\rm{evap}}}, \end{aligned}$$

(2)

where $\Delta z$ is a physical depth range over which gas with a mass density of $\rho$ has a bulk velocity of v.

Here, we ignore the downward momentum of proton and electron beam particles (Ichimoto and Kurokawa 1984; Allred et al. 2015). Typical values of the condensation are $\Delta z \approx 30$ km (which scales inversely with surface gravity; Kowalski and Allred 2018) and $v \approx -50$ km s$^{-1}$. Typical evaporation flows are fully ionized and have $n_e \approx 10^{11}$ cm$^{-3}$ ($\rho \approx 3 \times 10^{-13}$ g cm$^{-3}$) with flow speeds of $v \approx 100-500$ km s$^{-1}$. The density of the condensation increases over time, which is compensated by an increase in the $\Delta z$ of the evaporation as chromospheric/transition region mass ablates into and fills the newly-reconnected magnetic loop. After the flows fill up a $\Delta z \approx 10$ Mm loop, then Eq. (2) implies a gas density of $10^{-9}$ g cm$^{-3}$ in the condensation, in agreement with numerical simulations.⁹ The coronal loops thus emit in soft X-rays, while the footpoints emit in optical and UV radiation. The condensation manifests as red-shifts in broad, chromospheric lines (Fe II, Mg II, H$\alpha$, Si II) that exhibit a red-wing asymmetry (RWA). The RWA evolves in velocity and intensity over $\approx 30$ s (Ichimoto and Kurokawa 1984; Zarro and Canfield 1989; Wulser et al. 1992; Graham et al. 2020).

Fisher (1989) calculated the analytic relationship between the beam energy flux deposited above the flare transition region ($F_{{\rm{evap}}}$), the peak downflow (condensation) speed ($v_{{\rm{peak}}}$) just below the flare transition region, and the pre-flare chromospheric mass density ($\rho _{{\rm{chrom}}}$) at which the flare transition region forms:

$$\begin{aligned} v_{{\rm{peak}}} \approx 0.6 \bigg ( \frac{F_{{\rm{evap}}}}{\rho _{{\rm{chrom}}}} \bigg )^{1/3} \end{aligned}$$

(3)

Typical upflows of 100–500 km s$^{-1}$ imply durations up to $\approx 30{-}90$s to fill a flare loop, at which time the confined, field-aligned flows from conjugate footpoints collide. Superhot (lower volume and lower density) $T \approx 30{-}50$ MK sources can be produced above the looptops of the cooler, $T \approx 20$ MK, plasma that is thought to be due to chromospheric evaporation; these are thought to be due to heating from slow-mode (Petschek) shocks in reconnection, as demonstrated in recent multi-dimensional MHD simulations (e.g., Longcope et al. 2016). However, some of such sources that are consistent with ultrahot 100–200 MK thermal spectra that appear early on are generally more consistent with nonthermal power-law spectra.

In the impulsive phase, bright sources are observed from the footpoints (Fig. 2, right) of flare loops in cooler emission lines and in the UV, optical, and infrared continuum radiation. A variety of source morphologies appear and could be affected by temporal and spatial resolution. However, one can typically separate the geometry into elongated “ribbons” and compact, brighter circular “kernels” (see, for example, the remarkable H$\alpha$ images in Kawate et al. 2016). Within a relatively short time, the footpoints light up quasi-sequentially along the PIL, forming flare ribbons (e.g., Qiu et al. 2017; Kazachenko et al. 2022). One ribbon typically develops in plage,¹⁰ and another ribbon develops in a penumbral or umbral region of sunspot. As the flare progresses, the distance between the two ribbons increases, which is the so-called perpendicular apparent motion. This is thought to be a signature of the reconnection occurring higher and higher in the corona, thus releasing energy from larger and larger loops, which is a hallmark of the famous “CSHKP” model of two-ribbon flares (Carmichael 1964; Sturrock 1966; Hirayama 1974; Kopp and Pneuman 1976). As a result, the brightest emission is generally seen at the “leading bright edge” of the ribbons as they spread apart (Fig. 2 left). It is thought that excitation from above, e.g. through particle beams, changes location rather than through cross-field thermal heat transport/diffusion in the low atmosphere. However, radiative backwarming from the X-ray, EUV, and FUV lines may play an important role in heating the surrounding atmosphere away from the brightest kernels (such as in the so-called “core-halo” morphology; Neidig et al. 1993a; Allred et al. 2006; Isobe et al. 2007; Fisher et al. 2012; Namekata et al. 2022a). After the excitation front passes, the impulsively-generated chromospheric condensation continues to propagate even in the absence of electron beam heating, but it eventually reaches pressure equilibrium with the lower atmosphere. The loops cool from tens of MK to 1 MK, first primarily through thermal conduction, then primarily through radiation, resulting in the appearance of the famous “post-flare” loops as seen in filtergrams such as AIA 171 Å (which are, really, post-impulsive phase loops; see Fig. 3 and Liu et al. 2013a). In cooler lines such as H$\alpha$ or He II 304, the very late phase of each flare loop is sometimes associated with the phenomena of coronal rain, which is a second phase of “condensation”. In this phase, dense material that was chromospheric before the flare drains out of the corona. New flare loops are continuously formed through the gradual decay phase of GOES soft X-ray solar flares (Warren 2006), but the details of the partition among various possible heating sources in the low atmosphere (thermal heat flux,¹¹ electron beams, proton beams, direct heating by waves) remain poorly understood, especially during stellar flares. We note that the origin of the long cooling times of post-flare loops on the Sun is still a very active area of research (e.g., Ryan et al. 2013; Liu et al. 2013b; Qiu and Longcope 2016; Bian et al. 2018; Zhu et al. 2018; Kerr et al. 2020; Reep et al. 2022; Allred et al. 2022; Ashfield and Longcope 2023).

Several notable and well-studied solar flares are the SOL1992-Jan-13T17:25 M2.0 flare (“the Masuda flare”; Masuda et al. 1994), the SOL2000-Jul-14T10:00 X5.7 flare (“the Bastille Day flare”,¹² e.g., Qiu et al. 2010), the SOL2003-Oct-28T11:30 X17 flare (one among the “Halloween storms”, e.g., Woods et al. 2004), the SOL2002-Jul-23T00:30 X4.8 (“double power-law”; Holman et al. 2003) flare, the SOL2014-Mar-29T17:48 X1.0 flare (NASA’s “best observed X-class flare”; e.g., Kleint et al. 2016), the SOL2014-Sep-10T17:45 X1.6 flare (e.g., Graham and Cauzzi 2015), and the SOL2017-Sep-10T16:06 X8.2 flare (e.g., Chen et al. 2020).

Fig. 3

SDO/AIA images around the EUV wavelength of 171 Å (Lemen et al. 2012), obtained from Helioviewer.org, showing thermal emission from flare loops at two times during a large X-class flare on the solar limb. The top panel shows the compact, low-lying thermal flare arcade (indicated by the arrow) during what is usually the peak of the nonthermal hard X-ray (not shown) luminosity from the footpoints. Later in the flare, the bright arcade volume exhibits an expansion horizontally and vertically—the expansion is, however, apparent in the sense that newly reconnected, larger magnetic loops participate in the flare as time progresses. The large arcades have traditionally been called “post-flare” loops, which shine brightly much longer after the end of the main hard X-ray impulsive phase (they are perhaps more accurately called “post-nonthermal-X-ray-footpoint” loops). It is generally thought that the flare loops in these images are the result of magnetic structures that have cooled down from much higher temperatures to $T \approx 1$ MK (e.g. Aschwanden and Alexander 2001). Note that a cusp-like (fork-shaped) dark region just above the bright loops in the top panel is where many models and observations suggest that the bulk of magnetic reconnection and particle acceleration occur (see Longcope et al. 2018; Chen et al. 2020). For illustrative comparisons of flare and active region loops, see Güdel (2004)

4 A survey of flare stars

The characteristics that lead to flaring in stellar atmospheres are generally thought to be some combination of the following: rapid rotation, an outer convective zone, and disorganized surface magnetic fields. The ratio of stellar rotation to the convective turnover time¹³ is called the Rossby number, or the inverse Coriolis number, which may be indicative of how much internal shear, differential rotation, and turbulent amplification of kinetic energy into magnetic energy contribute to the dynamo. The dynamo may operate in the shear flows around the tachocline interface, and/or throughout the turbulent convection zone, and/or in the near-surface shear layers. Stars with small Rossby numbers are found in the so-called “saturated regime” of quiescent magnetic activity (outside of large flares). The saturated regime is empirically characterized by nearly constant luminosity ratios of $L_X/L_{{\rm{bol}}} \approx 10^{-3}$ (Fig. 4) and $L_{H\alpha }/L_{{\rm{bol}}} \approx 10^{-4}$ (e.g., Pizzolato et al. 2003; Wright et al. 2011; Reiners et al. 2014; West et al. 2015; Wright and Drake 2016; Newton et al. 2017; Brun and Browning 2017; Wright et al. 2018), where $L_{{\rm{bol}}}$ is the quiescent bolometric luminosity, and the numerator is a luminosity (e.g., the soft X-ray luminosity, $L_X$; Schmitt and Liefke 2004; Wright et al. 2011) that is a proxy for magnetic heating in the corona or chromosphere.¹⁴ There is a rather large spread about the saturated value from star to star, and the Sun varies from $L_X/L_{{\rm{bol}}} \approx 5 \times 10^{-8}$ to $\approx 5 \times 10^{-7}$ over its 11 year sunspot cycle (Judge et al. 2003; Ayres 2015a). For comparison, the least active M dwarfs have ratios of $10^{-6}$ (Wright and Drake 2016; France et al. 2020; Brown et al. 2023; Engle 2024). At very small Rossby numbers, a super-saturated regime may occur, while stellar activity proxies follow a power-law with increasing Rossby numbers $\gtrsim 0.1$, relatively independent of spectral type¹⁵. For an introduction to dynamo theory in the stars and the Sun, see Brandenburg (2005), Brandenburg and Subramanian (2005), Charbonneau (2013), Brun and Browning (2017), and Schrijver et al. (2019). For a review of the literature, see Brun and Browning (2017) and Käpylä et al. (2023). Much recent progress has been made in modeling stellar dynamos in low-mass stars just below (Bice and Toomre 2020, 2022) and above (Browning 2008; Brown et al. 2020; Käpylä 2021) the fully convective transition, and with especially high-resolution simulations of the solar magnetic field and differential rotation (Hotta et al. 2022).

Fig. 4

Quiescent $R_X = L_X/L_{{\rm{bol}}}$ vs. the Rossby number for a large sample of stars. $R_X$ is “saturated” at $\approx 10^{-3}$ for stars with Rossby numbers $\lesssim 0.1$, and $R_X$ decreases according to a power-law for lower-activity stars. Image reproduced with permission from Wright et al. (2011), copyright by AAS

Magnetic fields buoyantly rise through the stellar photosphere; once in the corona, they undergo destabilization and reconnection to release energy into the atmosphere to produce flares. Flaring occurs at different times through stellar evolution as convective envelopes come and go with changes in internal structure, core fusion processes, and external factors such as tidal locking and synchronous rotation. An overview of the myriad of flare stars is best summarized in a color-magnitude diagram (CMD). Prolific and/or representative flare stars are shown on a Gaia CMD in Fig. 5. The data were obtained from Gaia Data Release 2 (Gaia Collaboration et al. 2018), except for a few cases that are only available through the Early Data Release 3 or through common filter transformations (van Leeuwen 2007) for some of the brightest stars. The background data are the Gaia photometry for a collection of stars compiled for the Palomar/Michigan State University spectroscopic survey (PMSU; Reid et al. 1995; Hawley et al. 1996; Gizis et al. 2002; Reid et al. 2002; Reid and Hawley 2005), downloaded from Neill Reid’s website.¹⁶ Many of these stars are discussed in the recent 10 pc Gaia sample (Reylé et al. 2021) and the CNS5 catalog (Golovin et al. 2023). Isochrones from the PARSEC v1.2S (Bressan et al. 2012) stellar evolution models are shown for solar age and metallicity and the age of the $\beta$ Pic moving group (Mamajek and Bell 2014). Characteristics¹⁷ of selected flare stars in the CMD are summarized in Table 1, and in some cases a reference to an example of a flare study featuring the respective star.

Fig. 5

Gaia color-magnitude diagram of absolute magnitude vs. color for notable flare stars (outside of flaring). Note that corrections for extinction and reddening from interstellar dust are not applied because most of these stars are in close proximity (within $\sim 20$ pc) of Earth. The linecolors and linestyles for the annotations are arbitrary. Open circles signify that the data were transformed from Hipparcos magnitudes and/or parallaxes. Several solar-metallicity isochrones from the PARSEC stellar evolutionary models are shown for reference

Table 1 Table of common and notable flare stars

Stellar classification according to binarity, evolutionary status, and spectral type is important for understanding the origin and release of magnetic energy. We describe general classifications of the types of stars that flare across the CMD (Fig. 5). In the following “early” spectral types refer to those of hotter effective temperatures (e.g., M0-M2), and “later” spectral types refer to the cooler effective temperatures (e.g., M7-M9); this jargon has no direct correspondence to time or age. Mid-type M dwarfs span the fully convective transition and roughly correspond to stars with spectral types $\approx$ M3-M6, though there is over an order of magnitude range of bolometric luminosities just within this sub-grouping.

• RS CVn. Though rapid rotation is often associated with stellar youth, tidal locking of binaries can lead to prolonged magnetic activity and enhanced flare rates for billions of years (Osten et al. 2012) in systems that would not normally produce energetic events. RS CVn systems are synchronously rotating, detached binaries consisting of an early-type main sequence star and an evolved, cooler G/K sub-giant or giant companion. In RS CVn systems, the flares are thought to originate from the cooler component. The cool star component has a large convection zone, which combined with very rapid rotation for its size, leads to large flares with rise and decay times on the order of days. The flares of RS CVns are among the highest energy ($10^{36}-10^{37}$ erg) and longest-lasting stellar flares. Flare rates from the EUV and X-rays indicate about one flare per day (Osten and Brown 1999), but optical enhancements are relatively rare due to the enormous background glare from the non-flaring source; peak U-band changes (Eq. (1)) are typically much less than a magnitude for energies of $E_U>10^{33}$ erg (rates of $\approx 0.17$ h$^{-1}$; Mathioudakis et al. 1992); the largest flares increase the U-band flux by a factor of $\approx 2$ (Doyle et al. 1990b). The best-studied RS CVn system is arguably HR 1099 (V711 Tau), which consists of a primary K1 IV (sub-giant) and a G5 V (main sequence) star with orbital separation of $\approx 3R_{{\rm{Sun}}}$ and a 2.8 day synchronized orbital and rotation period (Fekel 1983). The most extensive multi-wavelength study of the flaring and quiescence of RS CVns is presented in Osten and Brown (1999); Osten et al. (2000, 2002, 2003, 2004). A large catalog of the non-flaring properties of these stars is contained in Strassmeier et al. (1993).

• Algol-type flare stars are semi-detached binaries with a history of significant mass transfer from a cool evolved sub-giant onto a main sequence B-type star. Algol (B8V + G8III/K1IV/K2III) consists of eclipsing binary (EB) stars separated by 14.14$R_{{\rm{Sun}}}$, whose components have radii of 2.9 and 3.5$R_{{\rm{Sun}}}$ with a 2.87-day orbital period. Algol produces giant flares on occasion (van den Oord et al. 1989; Favata and Schmitt 1999), and serendipitous eclipses of the flares constrain the geometries and locations of the events on the cooler, less-massive subgiant component, Algol B (Schmitt and Favata 1999; Schmitt et al. 2003), which is also the source of the coronal radio emission (Lestrade et al. 1993). Algol B is one of the few stars besides the Sun whose magnetosphere has been spatially resolved in radio observations (Peterson et al. 2010). Detailed comparisons of Algol-type and RS CVn systems have been presented in Singh et al. (1996); Sarna et al. (1998); Drake (2003). They are otherwise quite similar, except the total optical (non-flaring) light from the system is dominated by the main sequence B star in Algol-like systems, whereas RS CVns consist of two cool stars (spectral type G or later) with one subgiant or giant star that dominates the total system light.¹⁸ The different locations on the CMD (Fig. 5) are clear.

• BY Dra-type stars are, here, informally referred to as closely orbiting (EB or non-EB), detached binaries consisting of two main-sequence, typically late K or early M, stars. Arguably, the best known BY Dra flare star is the eclipsing 0.8 day period binary YY Gem, which consists of nearly identical, and probably rotationally synchronized, dM0e stars with radii of $\approx 0.6R_{{\rm{Sun}}}$ and masses of $\approx 0.6M_{{\rm{Sun}}}$ separated by 3.83 $R_{{\rm{Sun}}}$ (Torres and Ribas 2002; Feiden and Chaboyer 2013; Butler et al. 2015; Kochukhov and Shulyak 2019). Other eclipsing BY Dra flare star systems are CM Dra and CU Cnc. BY Dra itself is composed of two active main-sequence K stars with semi-major axes of $7.4-8.4 \times R_{{\rm{Sun}}}$ and pseudo-synchronous rotation with a period of 5.9 d (Hełminiak et al. 2012). Formally, BY Dra is a term that describes the class of magnetically active spotted stars that exhibit periodic variation in their optical light curves, due to starspots and rotational modulation (General Catalogue of Variable Stars: Version GCVS 5.1; Samus’ et al. 2017)¹⁹); however, the BY Dra classification more uniquely describes (e.g., Reid and Hawley 2005) detached active binary systems that consist of main sequence late-type (K or M) stars that orbit very closely together,²⁰ like BY Dra itself.

• Pre-main sequence (PMS) stars, such as Young stellar objects (YSO’s) and T Tauri stars exhibit flaring emission (e.g., Flaccomio et al. 2018) in addition to the emission from accretion: the flow, shocks, and hotspots at the stellar surface all generate optical, UV, IR, and X-ray continuum and emission line radiation (e.g. Herczeg and Hillenbrand 2008) that can be both steady-state and transient. The optical, NUV, and FUV flares on T Tauri stars tend to be very energetic, $\gtrsim 10^{35}$ erg (Tofflemire et al. 2017; Hinton et al. 2022; Getman et al. 2023). Tofflemire et al. (2017) conducted a large optical photometric monitoring campaign on the T Tauri binary star DQ Tau and compared the color differences between accretion and flare variations. Getman et al. (2011) performed an extensive study of X-ray giant flares in the DQ Tau system. The Chandra Orion Ultradeep Project (COUP; Getman et al. 2005) detected many X-ray superflares from YSOs; the flares produced X-ray hardness ratios consistent with coronal temperatures of $T >10^8$ K (Getman et al. 2008). M-type stars take longer to contract onto the main sequence after losing their gaseous disks and are thus not actively accreting, but may exhibit inflated radii for tens of Myr (e.g., AU Mic) as they are still contracting. These are still PMS stars but are also often considered as dwarf, UV Ceti-type stars.

• UV Ceti-type stars are low-mass, M-type main sequence (dwarf) flare stars that exhibit frequent flaring and H$\alpha$ in emission in quiescence (Gershberg et al. 1999). There is a high correspondence between the dMe status and high flare rates (Kowalski et al. 2009). These stars are traditionally denoted as “dMe” or “MVe” stars, they exhibit optical rotational modulation (with periods typically in the 0.2–5 day range, but see West et al. 2015 and Newton et al. 2017), and many but not all are in the saturated regime ($L_x/L_{{\rm{bol}}} \approx 10^{-3}$) of quiescent X-ray luminosity. Many of the UV Cet/dMe-type stars are in the early stages of their main sequence evolution and are possibly only several hundred Myr or less in age. UV Ceti itself is a triple star system with two dM5.5e stars (A and B components) separated by 5.3 au (Benz et al. 1998); UV Ceti B is a binary BY Dra-type system (using our adopted definition above). This category includes some non-accreting PMS stars (AU Mic, AT Mic). A catalog of UV Ceti stars and their quiescent properties was compiled by Gershberg et al. (1999) and is provided on VizieR.²¹ Though the M-dwarfs are the most populous stars in the Galaxy (e.g., Reid and Hawley 2005), none are visible to the naked eye ($V \lesssim 6$) from Earth due to their low quiescent luminosities.

• Very low-mass (VLM) stars, ultracool dwarfs, Brown dwarfs. The spectral types M7 through early L span the edge of the core-hydrogen burning regime (i.e., the star/brown-dwarf boundary). The exoplanet host star TRAPPIST-1 is now a well-known, M8 flare star, due to its system of planets in or around the traditional habitable zone (Gillon et al. 2017; Agol et al. 2021). VB 10 is another prolific VLM flare star (Herbig 1956; Linsky et al. 1995; Fleming et al. 2000; Kanodia et al. 2022). The Güdel-Benz relationship (Güdel and Benz 1993) is a correlation between the quiescent radio and soft X-ray fluxes that holds over many orders of magnitudes but breaks down in the VLM regime (with the VLMs being radio overluminous compared to the X-rays). Nonetheless, flares from VLMs and L-type stars have shown dramatic optical continuum enhancements and broad and bright H$\alpha$ lines (Liebert et al. 1999; Schmidt et al. 2007; Gizis et al. 2013) that otherwise resemble the very energetic flares from the earlier M dwarf spectral types, as noted by Reid and Hawley (2005). White-light flare rates and magnitude variations have been rather well-characterized recently in time-domain surveys (e.g., ASAS-SN, Kepler, K2, NGTS; Schmidt et al. 2014, 2016; Gizis et al. 2017; Vida et al. 2017; Paudel et al. 2018, 2019; Jackman et al. 2019a; Paudel et al. 2020). The L-dwarf flare luminosities are remarkably super-bolometric and exceed magnitude changes of $-9$ in the V-band (Schmidt et al. 2014). X-ray superflares have been reported from L dwarfs (e.g., De Luca et al. 2020) as well.

• Weakly active and “inactive” M dwarfs (dM or MV) stars are “optically inactive” stars without H$\alpha$ in emission during quiescence. The inactive dM’s also flare on occasion (Paulson et al. 2006; Kowalski et al. 2009). These are thought to be descendants of MVe stars, which are inferred to hold on to their high levels of activity for several billion years (Hawley et al. 2000; West et al. 2008; Kiman et al. 2021). After several Gyr on the main sequence, the M-dwarf stars maintain a lower level of chromospheric emission in Ca II H & K (with sometimes H$\alpha$ showing a deeper absorption profile), Mg II h & k emission lines, transition region emission (e.g., C II), and faint X-rays. One notable example of a low-activity M dwarf is the $\approx 10$ Gyr-old, slowly rotating ($P \approx 150$ d) star Gl 699 (Barnard’s star; Walkowicz and Hawley 2009; Fontenla et al. 2016; Toledo-Padrón et al. 2019; France et al. 2020). Optical flares are rare (Hilton 2011; Hawley et al. 2014) on optically inactive stars, which do not exhibit much if any rotational modulation in Kepler. The visibility (contrast) of flares at shorter wavelengths is much greater, however, and thus the FUV and X-ray flares are rather prolific (Loyd and France 2014; France et al. 2016; Loyd et al. 2018a; Froning et al. 2019; France et al. 2020; Brown et al. 2023).

• Solar-type stars Flares from solar-type, G-dwarf stars have been observed rather infrequently, and mostly through serendipitous means in the FUV (Ayres et al. 1994), X-ray (Getman et al. 2008; Pye et al. 2015), and optical (Schaefer et al. 2000; Maehara et al. 2012). The detectable flares above the spatially integrated, background luminosity are very energetic, and in almost all cases they exceed broadband optical energies of $E > 10^{33}$ erg. Detailed flare rates in the EUV have been studied for the rapidly rotating, young, single G-type stars, EK Dra, 47 Cas, and $\kappa$ Ceti (Audard et al. 2000).

  The Kepler and K2 missions have transformed our knowledge about the white-light flare rates of G-type stars. The literature has been summarized in detail by Cliver et al. (2022b) and Hayakawa et al. (2023b), and we refer the reader to these works. Here, we provide only a brief synopsis. After the discovery of G-dwarf white-light superflares in Kepler (Maehara et al. 2012), Shibayama et al. (2013) and Notsu et al. (2013b) performed a comprehensive flare rate analysis of the 30-min cadence data (with some comparisons to the 1-min cadence data), and Candelaresi et al. (2014) extended analyses to lower mass K and M stars (see also Walkowicz et al. 2011, for analyses of 30-min cadence data of M star flares). Yang et al. (2018) also compares short and long cadence Kepler data. Notsu et al. (2013a), Nogami et al. (2014), and Karoff et al. (2016) began the detailed spectroscopic characterization of the chromospheric and photospheric properties of the solar-type superflare stars reported in Maehara et al. (2012) and Shibayama et al. (2013). Notsu et al. (2015a) and Notsu et al. (2015b) presented spectroscopic follow-up of 46 of the flare stars with 30-min cadence data in Shibayama et al. (2013). Maehara et al. (2015) analyzed 1-min cadence data of 23 solar-type superflare stars. The most prolific flare star (KIC 11551430; Table 1) was confirmed as a visual and spectroscopic binary in Notsu et al. (2019), who presented spectroscopic follow-up of 18 of these 23 stars, leveraging Gaia DR2 data as well to constrain evolutionary statuses on the subgiant branch or main sequence. Okamoto et al. (2021) re-analyzed all solar-type Kepler data and presented up-to-date flare rates. As a result of the group’s detailed follow-up, the calculated superflare occurrence rates on Sun-like stars ($T_{{\rm{eff}}} = 5600-6000$ K, $P_{{\rm{rot}}} = 20{-}40$ d) has decreased by $\approx$ an order of magnitude since Shibayama et al. (2013). These refined constraints provide interesting comparisons to extrapolated, multi-wavelength scaling relations from smaller solar flare energies and from large solar energetic particle events (see also, e.g., Figure 33 of Usoskin 2023). Wu et al. (2015) presented flare rates and power-law indices for 77 individual G stars that produce superflares, including KIC 10422252 which is the prolific flare star that is discussed in Shibayama et al. (2013). Davenport (2016), Van Doorsselaere et al. (2017), and Yang and Liu (2019) present Kepler flare catalogs, where the latter includes all spectral types including A-type stars (see below) and only long-cadence data. Lawson et al. (2019) discussed source contamination in the Davenport (2016) catalog in the higher mass stellar range. Recently, Aschwanden and Güdel (2021) analyzed the power-law distributions of the Kepler long-cadence data flare catalog of Yang and Liu (2019) and tested against self-organized criticality theory (Aschwanden 2014).

• Single active giants (Simon and Drake 1989) are classified into one of two general categories (see Ayres et al. 1998). The first type is helium core burning (post Helium flash) red clump active giants like $\beta$ Ceti. Giant flare events with rise times of $\sim 1$ day are observed from these stars: an impressive EUV light curve of $\beta$ Ceti is shown in Fig. 6 from Ayres et al. (2001). The active clump stars are rotating slower than evolutionary models predict, yet they have detectable magnetic fields (unlike the rest of this stellar population). The origin of the magnetic activity in these stars is not well understood, but it has been hypothesized that swallowing a companion planet or brown dwarf increases the shear in the stellar interior (Siess and Livio 1999), resulting in “magnetic rejuvenation” in old age. An alternative idea is that these stars, like some other active giant stars (Stepien 1993), are descendants of Ap stars with strong fossil fields (Tsvetkova et al. 2013). The second type of active giants occur in the Hertzsprung Gap, which is a short-lived phase of stellar evolution after rapidly rotating A and B-type stars leave the main sequence before rotationally breaking down and ascending the red giant branch. White-light flares have been widely reported in Kepler data of subgiants and giants (Notsu et al. 2019; Okamoto et al. 2021; Katsova et al. 2018; Kővári et al. 2020; Oláh et al. 2022), and some exceed $10^{38}$ erg from KIC 2852961, whose binary status is not yet clear.

• B5-F5 single stars have convective cores and radiative exteriors and are not expected to produce complex surface magnetic fields and flares. B stars and Ap stars are magnetic but the magnetism is thought to be a remnant of fossil fields, which are not disorganized enough to facilitate reconnection and impulsive energy release. The F2 star HR 120 (Mullan and Mathioudakis 2000) is a surprisingly early-type flare star that has a detailed EUV flare rate presented in Audard et al. (2000). White-light superflares from A-type stars have been reported in Kepler data (Balona 2012, 2019), and spectroscopic follow-up by Pedersen et al. (2017) found evidence for binarity in most of these sources. If originating from the companion, Balona (2012) argues that the energy requirements are yet rather large and unreasonable (see also Mullan 2009). The A7 star Altair is magnetically active, and its dynamo may operate locally in an equatorial zone that is cool enough for convection to occur due to the star’s very fast rotation (Robrade and Schmitt 2009). But to our knowledge, this “backyard star” is not known to produce white-light superflares. Algol-like systems consist of a B- or A-type star, but the flares originate from the cooler, evolved star.

Fig. 6

A small fraction of single, red clump giant stars like $\beta$ Ceti exhibit continuous variability that consists of week-long, giant EUV (60–180 Å) flares. A flare on $\mu$ Velorum (also shown) is even more unexpected in its evolutionary stage. Image reproduced with permission from Ayres et al. (2001), copyright by AAS

5 Flare rates and flare frequency distributions

The rates at which stars flare have been measured through continuous ground-based observations, continuous space-based observations, and coarse sampling with serendipitous detections. The coarse sampling may result in just one data point per flare; thus an accurate characterization of the flaring source (Sect. 4) is often desired through follow-up spectroscopy. The study of flare rates is important for a wide variety of wider applications, including coronal heating physics (Hudson 1991), characterizing serendipitous variability in time-domain and exoplanet surveys (Becker et al. 2004; Welsh et al. 2007; Hilton et al. 2010; Berger et al. 2013; Gezari et al. 2013; Hawley et al. 2016; Fuhrmeister et al. 2018; Mondrik et al. 2019; Howard et al. 2019; Chang et al. 2020; Jackman et al. 2021; Rodríguez Martínez et al. 2020; Koller et al. 2021; Webb et al. 2021; Wu et al. 2022), and studying the effects on exoplanet environments (e.g., Howard et al. 2018; Tilley et al. 2019). Much knowledge of stellar flares from low-mass stars is a result of monitoring a handful of nearby dMe stars, which reliably flare at a high rate. Large etendue²² time-domain capabilities of wide-field surveys such as the Sloan Digital Sky Survey, 2MASS, and ASAS-SN have been leveraged to characterize flare rates of inactive and early-type M dwarfs (Kowalski et al. 2009), at infrared wavelengths (Davenport et al. 2012), and of ultracool dwarfs (Schmidt et al. 2014, 2016, 2019).

5.1 Basic methods and early studies

Hilton et al. (2010) distinguish among flare rates, flare duty cycles (flaring fractions), and flare frequency distributions. Here, we focus on flare frequency distributions (FFDs), which refer to the energy dependence of the flare rates from continuous monitoring. Flaring fractions in sparsely sampled data are discussed in Sect. 5.3.

The seminal study of Lacy et al. (1976) presented an extensive flare rate analysis of several nearby active M dwarf (dMe) stars from hundreds of hours of ground-based monitoring in typical Johnson bandpasses. The high-time resolution data and observational details, along with some fundamental measurements such as equivalent durations, timescales, and peak amplitudes, can be found in Moffett (1974). In this section, we first present several basic measurements that are derived from flare light curves. The most fundamental calculation is the equivalent duration of a flare through a bandpass. The equivalent duration of a flare is the time that a star must spend in quiescence to produce the same energy as in the flare. It is a proxy of the flare energy, which is spectral type and luminosity class dependent, but can be readily obtained from relative photometry. The equivalent duration is (Gershberg 1972):

$$\begin{aligned} ED = \int _{t_{{\rm{start}}}}^{t_{{\rm{end}}}} {\frac{I(t)- I_q}{I_q}} \,dt = \int _{t_{{\rm{start}}}}^{t_{{\rm{end}}}} {I_f(t)} \,dt \end{aligned}$$

(4)

which has units of time and requires establishing flare start and end times. A small equivalent duration can result from a small flare energy and/or a large quiescent count flux ($I_q$) at Earth. Note that the symbol for intensity, I, is traditionally used, though this quantity is actually the spatially unresolved instrumental (I) count flux from the star²³. The equivalent duration is multiplied by quiescent luminosity through this bandpass to obtain the absolute flare energy.²⁴

The integrand of the equivalent duration is known as the “flare visibility” or “flare contrast” and is denoted as $I_f(t)$. The flux enhancement relative to the flare star’s quiescence is known as $I_f(t) + 1$ (Sect. 2). The numerator in the integrand is the “flare excess”, which can be denoted using a prime symbol to indicate a change. If the average flux is calculated from a continuum or pseudo-continuum region of a spectrum, then the flare-only flux is denoted as C$\lambda _c^{\prime }$, where $\lambda _c$ is a representative wavelength within the spectral window. The I quantity is formally the count flux of the flare star normalized to the count flux from a nearby non-variable star ($I(t) = C_{\rm{rel}} (t)$ in Eq. (1)) in the same exposure and same aperture radius. For a photon-counting instrument, such as a CCD, the theoretical value of I is

$$\begin{aligned} I(t) = \frac{\int {T(\lambda ) f_{\lambda } (\lambda , t) \lambda \, d\lambda }}{\int {T(\lambda ) f_{\lambda , 0} (\lambda ) \lambda \, d\lambda }} \end{aligned}$$

(5)

where $T(\lambda )$ is the total system response or effective area, including filter, atmospheric transmission (if applicable), and the detector gain; $f_{\lambda ,0}$ is the flux spectrum of a comparison (non-variable) star or stars, and $f_{\lambda }$ is the flux spectrum of the target star in units²⁵ in $\hbox {erg cm}^{-2}\, \hbox {s}^{-1} \text{\AA }^{-1}$. The ED is then multiplied by the bandpass (T) quiescent luminosity ($L_{q,T}$) to find the bandpass-integrated flare energy, $E_{T} = ED \times L_{q,T}$. For known and well-characterized bandpasses, these can be determined using zeropoint fluxes, $f_{{\rm{ZP}}}$ (Willmer 2018) and published apparent magnitudes²⁶, $m_T$, at low-levels of flare activity:

$$\begin{aligned} m_T = -2.5 \, \log _{10} \frac{\langle f_{q,\lambda } \rangle _T}{f_{\rm{ZP}, T}} \end{aligned}$$

(6)

where the zeropoint flux in the Vega mag system is determined by the numerical integration of a bandpass over the spectrum of the A0 V spectrophotometric standard star, Vega. If a flux-calibrated quiescent spectrum of the flare star is available, numerical integration over the total system response, $T(\lambda )$, including bandpass, alternatively yields the following:

$$\begin{aligned} L_{q,T} = \langle f_{q,\lambda } \rangle _T \, \Delta \lambda \, 4 \pi d^2 = \frac{\int {T(\lambda ) f_{q,\lambda } (\lambda ) \lambda } \,d\lambda }{\int {T (\lambda ) \lambda } \,d\lambda } \Delta \lambda \, 4 \pi d^2, \end{aligned}$$

(7)

where $\langle {f_{q,\lambda }} \rangle _T$ is the system-weighted flux (Sirianni et al. 2005), $f_{q,\lambda }(\lambda$) is the quiescent stellar spectrum at Earth in units of $\hbox {erg cm}^{-2}\, \hbox {s}^{-1} \text{\AA }^{-1}$, and d is the distance to the star. The width of the bandpass, $\Delta \lambda$, is usually taken to be the FWHM or the effective width, which vary from $\Delta \lambda \approx 600$ Å to 4000 Å, for the SDSS u-band and Kepler bands, respectively. In recent years, quasi²⁷-bolometric conversions from white-light filters have assumed a blackbody function for all flares, and all phases of flares (e.g., Shibayama et al. 2013). Other studies assume constant spectral energy distributions, with energies reported over the limits of the bandpass (e.g., Hawley et al. 2014). Several quiescent luminosities used for M dwarfs are given in Table 16 by Moffett (1974), but it should be noted that stars, especially active M dwarfs, have a rather large spread of U-band magnitudes reported in the literature (e.g., reported quiescent apparent U-band magnitudes for the flare star YZ CMi range between 13.70–13.85). Total system response characteristics through U bandpasses are also inherently uncertain (e.g. Maíz Apellániz 2006) and to some extent, intractable. Thus, non-negligible systematic uncertainties in absolute energies are to be expected from a combination of these various issues (see Hilton 2011, for empirical comparisons of flare energies in SDSS u and Johnson U from two different telescopes). In Appendix A, we show several of the broadband and narrowband filters that are commonly employed in optical stellar flare studies.

FFDs are presented as either downward cumulative, $Q(>E)$, or differential, n(E), distributions in log$_{10}$-log$_{10}$ space (see Audard et al. 2000, for comparisons of these methods). The most common practice is to fit a power-law (pareto) model to the unbinned, downward cumulative FFD:

$$\begin{aligned} Q(>E) = N \bigg (\frac{E}{E_0} \bigg )^{\beta } \end{aligned}$$

(8)

with power-law index, $\beta$, and N total flares greater than or equal to a fiducial low-energy limit, $E_0$ (e.g., the detection completeness limit) observed within a monitoring duration of $\Delta t$ (hr). The usual convention is that $\beta$ is negative. The differential FFD (# of flares per unit energy) is then

$$\begin{aligned} n(E) = - \frac{dQ}{dE} = N \frac{\alpha - 1}{E_0} \bigg (\frac{E}{E_0} \bigg )^{-\alpha } \end{aligned}$$

(9)

where $\beta = 1 - \alpha$ (again, the convention is that $\beta$ is negative and $\alpha$ is positive). Note that if $A \int _{E_0}^{\infty } n(E)dE = 1$ is solved for A, the converted differential distribution $A\ n(E)$ is a probability density function (PDF): given a flare occurred, $A\ n(E)dE$ is the probability for a flare to have an energy between E and $E+dE$.

Convenient analytic forms exist for the power-law index and its statistical uncertainty from a direct (i.e., unbinned) maximum likelihood (ML) analysis (Wall and Jenkins 2003; Clauset et al. 2009, Appendix B):

$$\begin{aligned} {\hat{\beta }}_{\text {ML}} = \frac{N}{\sum _{i=1}^N ln \frac{E_i}{E_0}} \end{aligned}$$

(10)

from which the result follows that $\sigma _{{\hat{\beta }}_{\text {ML}}} \approx {\hat{\beta }}_{\text {ML}}/\sqrt{N}$.

Alternatively, weighted least-squares fits of the power-law index (slope) and intercept of a line in log-log space has been sometimes employed in stellar flare studies. For example, Lacy et al. (1976) express the power-law for each flare star in their sample of active M dwarfs as

$$\begin{aligned} \log _{10} \nu (>E) = a + \beta \, \log _{10} \, E \end{aligned}$$

(11)

where $\nu (>E) = Q(>E)/\Delta t$ is the number of flares per hour greater than or equal to energy E (expressed in erg) through an optical bandpass T. (To increase statistical confidence, flares that are observed in several bandpasses were multiply-counted using empirical bandpass energy conversions with the effective observing time adjusted accordingly.) Of course, the counts in even well-populated cumulative distributions have statistical uncertainties that are asymmetric in logarithmic space, and the values in cumulative distributions are non-independent. Nevertheless, this method is employed as a sufficient approximation in some cases given the statistics that are possible with the relatively low number of flares in ground-based observing campaigns.

Power-laws are typically fit to FFDs over an intermediate range of energy that is carefully chosen (e.g., Silverberg et al. 2016) to exclude the high-energy and low-energy ends, which may have incomplete sampling that can affect the fits. Note that Clauset et al. (2009) present more sophisticated methods for uncertainty estimation for power-laws and have been employed in some more recent flare rate studies (Medina et al. 2020). The Markov Chain Monte Carlo approach is described in Davenport et al. (2016), and Davenport et al. (2019) use an appropriate uncertainty analysis for the counts (Gehrels 1986). Modifications to the basic power-law form (Eq. (11)) and completeness functions have been used to make adjustments to the the low-energy end (Rosner and Vaiana 1978; Aschwanden 2015; Davenport et al. 2012; Medina et al. 2020; Okamoto et al. 2021) and high-energy tail (Aschwanden and Güdel 2021) of FFDs. For reference, solar flare frequency distributions (Cliver et al. 2022b) exhibit power-law indices in (nonthermal) hard X-ray energy of $\approx 1.5$ and peak luminosity of $\approx 1.7$ (Crosby et al. 1993), whereas the thermal soft X-ray peak flux exhibits a larger power-law index of 2.0 (Veronig et al. 2002a; Aschwanden and Freeland 2012; Hudson et al. 2024).

We summarize the general range of values of $\alpha$ calculated from FFDs in the X-ray and the optical regimes. Ground-based optical and space-based XEUV (0.01–10 keV, which corresponds to 1.24–1240 Å, and where only EUV data, such as DS count rates, are available, a two-temperature, $T = 6+23$ MK, optically thin free-free model extrapolates to the full XEUV range to calculate energies) flare rates have generally given values from $\alpha \approx 1.4{-}2.2$ with statistical uncertainties of $\gtrsim 0.1$ (Audard et al. 2000). RS CVn systems exhibit power-law indices in the EUV on the lower end of this distribution (Osten and Brown 1999). Wu et al. (2015) summarize the power-laws from 77 G-type stars. Hawley et al. (2014) and Ilin et al. (2021) summarize values of $\alpha$ for low-mass stars in the optical, and Loyd et al. (2018a) present power-laws for low-mass stars of different activity levels in the FUV. A useful document that summarizes flare rates from low mass M stars of various spectral types is Osten (2016).²⁸ Among the active M dwarfs, the higher luminosity stars (in quiescence) exhibit flatter power-laws (Pettersen et al. 1984), and thus the sum of the energy in all flares is dominated by the occurrence of the highest-energy flares (Lacy et al. 1976). Alternatives to the power-law functional form (Eq. (9)) are summarized in Aschwanden (2011), Rosner and Vaiana (1978), and Sakurai (2022).

A few stars have had power-law indices measured independently several times; these stars show some interesting systematic differences in the reported power-law properties. In the M dwarf study of Lacy et al. (1976), the famous flare star AD Leo appeared to be an outlier from the others. Due to the relatively short observing time in the Lacy et al. (1976), its FFD was recalculated after more data were collected. Pettersen et al. (1984) found a value of $\alpha = 1.62 \pm 0.09$ from 85 U-band flare events (with 115 individual flare peaks) in 111 hr of monitoring. The most energetic flare was $10^{33}$ erg. The U-band FFD of Proxima Centauri has been quantified by Walker (1981) for flare energies between $5\times 10^{27}$ erg and $10^{30}$ erg, giving a similar value of $\alpha \approx 1.7$ to AD Leo but flatter than other stars of the same quiescent luminosity, such as CN Leo. Recently, Davenport et al. (2016) doubled the number of flares in the statistics for Prox Cen using white-light optical data from the MOST satellite and calculated a similar power-law slope to Walker (1981). Howard et al. (2018) found a steeper power-law for Proxima Centauri from the Evryscope survey and investigated the effect on the power-law index with and without an extreme superflare outlier.

Lacy et al. (1976) and Pettersen et al. (1984) found trends in power-law indices with effective temperature or quiescent U-band luminosity: larger luminosity (earlier-type) main-sequence stars exhibiting flatter power-laws (i.e., smaller $\alpha$, $\mid \beta \mid$). However, the study of Pettersen et al. (1984) did not include the results of Walker (1981) for the interesting case of Proxima Centauri, which shows a flatter power-law than expected; evidence for and against Proxima Centauri being less magnetically active than others of the same sub-type, such as CN Leo, are summarized in Reiners and Basri (2008). Audard et al. (2000) suggested that $\alpha$ derived from the XEUV may exhibit an opposite trend with spectral type, though this is reported with low confidence. They find an increasing rate of $E > 10^{32}$ erg flares corresponds with greater quiescent X-ray luminosity.

Other interesting quantities that have been calculated from flare star monitoring are the energy per flare (average flare energy) and the radiated flare energy per unit time (average luminosity due to flaring). These quantities are usually calculated through the U band. They correlate with the non-flaring, bandpass quiescent luminosity (Lacy et al. 1976) or non-flaring, bolometric luminosity of the star over several orders of magnitude. It should be noted that the number of detectable flares is anti-correlated with the quiescent luminosity, an effect largely due to brighter stellar background that raises the detection limit (e.g., Davenport et al. 2012). Nonetheless, including many small, unobserved flares in the total energy calculation does not change the trends. From modern space-based data, it seems that these relationships generally extend to active binaries, which comprise much higher luminosity stellar components (Osten et al. 2012). An extrapolation of these relationships to the Sun, however, dramatically over-predicts the flare energy release rate in optical solar flares by many orders of magnitude, as first noted by Lacy et al. (1976). From Kepler data (Sect. 5.2), a somewhat similar quantity of $L_f$ vs. $L_{{\rm{Kp}}}$ has been constructed by Lurie et al. (2015) and is making for powerful comparisons for normalized flaring efficiencies among stars of different spectral types (e.g., Davenport et al. 2019). Davenport (2016) investigates evidence for saturation in the relative flare luminosity given by $L_f/L_{{\rm{Kp}}}$ at small Rossby numbers, analogous to the quiescent $L_X/L_{{\rm{bol}}}$ saturation (Sect. 4).

There has been a severe paucity of low-activity M dwarf flare rates until Hawley et al. (2014) presented white-light 1-min cadence data from Kepler and some previously unpublished ground-based U-band data from the PhD work within Hilton (2011). These studies revealed that stars that are optically inactive in quiescence (Sect. 4) produce broadband optical flares, a result that had been quantified in sparsely sampled data from the Sloan Digital Sky Survey (Kowalski et al. 2009) and from serendipitous spectral detections (Paulson et al. 2006). Loyd et al. (2018a) compared FUV flare energies among stars with comparable quiescent bolometric luminosities and found that FFDs are shifted according to the quiescent flux in the FUV.

The different apparent behaviors at the high-energy ends of stellar flare FFDs are still largely unexplained. Does the power-law that is well-fit to a middle range of energies continue on indefinitely, or does it turn over at the high-energy end and follow a steeper power-law index (Doyle and Mathioudakis 1990; Osten et al. 2012) or another type of roll-over (e.g., Osten et al. 2004; Dal 2020; Aschwanden and Güdel 2021; Cliver et al. 2022b; Sakurai 2022) for some stars? What physical parameters determine the upper limits of flare energies (and peak luminosities) that are possible on a given star? Extending the power-law that is fit to the FFD of YZ CMi predicts a U-band flare to occur about once per month with an energy that is a factor of 5000–10,000 larger than its mean flare energy (Lacy et al. 1976; Kowalski et al. 2010). Dal and Evren (2011b) discuss an upper limit to the energy released in dMe flares is approached as the flares become longer in total duration. Shibata et al. (2013) and Aulanier et al. (2013) present theoretical grounds and semi-empirical relations for upper limits of flare energies. Empirical constraints from very long monitoring times have only been possible recently using Kepler data.

5.2 White-light flare rates & Kepler

The release of Kepler data (Borucki et al. 2010) has transformed knowledge of flare rate statistics in several important ways. First, the high precision and long monitoring times facilitate detection of relatively rare white-light flares from older, less active low-mass stars, from much larger luminosity stars such as solar-type stars, and from higher mass K and early-type M dwarf stars, and evolved stars. Second, it provided many more hours of monitoring per star: ground-based monitoring efforts resulted in $\approx 30{-}100$ hours per star, while the nominal baseline for Kepler was two contiguous months of observations per star. The flare statistics correspondingly improved, but the high precision also revealed comparable amounts of background flux variation due to starspots and rotation. Nonetheless, the Kepler field began as a relatively poorly known star field (due to its relatively faint magnitude limits), and systematic contamination of the flare sources from binaries, subgiants, and other non-flare sources was an issue until targeted, ground-based spectroscopic (and asteroseismic) analyses caught up and Gaia DR2 parallaxes became available. Additionally, the vast majority of data from the Kepler mission is at 30-min cadence, and most superflares span just two to four data points. The 1-min cadence data were made possible for select stars through Guest Observer programs (e.g., Hawley et al. 2014).

The dM4e star GJ 1243 has been observed for 11 months at 1-min cadence with Kepler, making it the longest observed flare star in the sky besides the Sun. The FFDs and flare properties of this star have been studied extensively by Ramsay et al. (2013), Hawley et al. (2014), Davenport et al. (2014), Davenport (2016), Silverberg et al. (2016). From the database of over 6000 flares, Silverberg et al. (2016) obtained a value of $\alpha = 2.008 \pm 0.002$, and Davenport et al. (2020) inferred $\alpha = 1.942 \pm 0.001$. From the first two months of data, Hawley et al. (2014) found $\alpha = 2.01$. They also noted that the turnover in the FFD at the low-energy end, which has often been attributed to a detection limit effect, may be astrophysical at some level since flares of these energies should be readily detectable above the precision of Kepler. The FFD of GJ 1243 is shown in Fig. 7(top) in comparison to several other M dwarfs spanning different spectral types and levels of quiescent magnetic activity (dM v. dMe). Comparable diagrams with compilations of FFDs of many individual stars are shown in Shakhovskaia (1989), Ramsay et al. (2013), and Dal (2020). Lurie et al. (2015) presented a detailed analysis of the GJ 1243 AB system, where they separated the Kepler light into the two stellar components and found values of $\alpha \approx 2$ for both stars, but a higher flare rate occurs on the slower rotator. The bottom panel of Fig. 7 shows the Kepler energies converted to energies in the U-bandpass. The flare rates of “less-active” M3–M5 stars are closer to the FFD of GJ 1243, whereas the very-active M3-M5 stars include the dMe stars from Lacy et al. (1976), such as YZ CMi, EV Lac, and AD Leo. A continuum of flare rates is clearly evident as the vertical offsets of the power-laws, while different power-law slopes may be a result of vast differences in statistics from ground and space-based monitoring (Hawley et al. 2014).

Fig. 7

Flare frequency distributions of several main-sequence M stars in the Kepler band (top) and in the U-band (bottom). Note the spread in flare rates in the bottom panel among the very active and less active stars of the same spectral types. Image reproduced with permission from Hawley et al. (2014), copyright by AAS

The spread in flare rates among active M-dwarfs of similar spectral types (Fig. 7) contrasts with the findings of Kowalski et al. (2009), who found that the average M4Ve-M5Ve flare frequency distribution from sparsely sampled data of field stars in the Sloan Digital Sky Survey (SDSS) Stripe 82 is sensibly in line with that of the well-known flare star YZ CMi (which is one of those among the “more active” stars in Fig. 7). Hilton (2011) began to investigate the chance alignment of FFDs in Kowalski et al. (2009) using Monte Carlo sampling of simulated flares at the SDSS riuzg cadence. The resolution of this issue (which is a non-trivial adjustment to the statistical rate calculation; E. Hilton, priv. communication 2013) will be important for accurately mapping flare rates of the Galaxy using sparsely sampled data from the LSST (Ivezić et al. 2019).

The immaculate flare statistics of Kepler facilitated new timing analyses of flare occurrence and correlation. The flare occurrence of GJ 1243 was investigated as a function of the phase of the longer-timescale, flux modulation that exhibits a $\approx 1$% amplitude and is attributed to stellar rotation. No significant correlation was found at either maximum or minimum flux levels. The relationship between flare occurrence and phase was further investigated in detail with larger samples of flares on GJ 1243 (Silverberg et al. 2016) and for a large sample of M dwarfs (Lurie et al. 2015; Doyle et al. 2018). Morris et al. (2018) reported a larger flare occurrence when the very low-mass (M8V) star TRAPPIST-1 is slightly brighter than average and when it goes through the increasing flux phase in the Kepler light curve. Dal and Evren (2011a) investigated the occurrence of different types of flares with phase from a large sample of ground-based observations of flares.

Hawley et al. (2014) finds that the intervals between successive flares (“waiting times”) on an active M dwarf (GJ 1243) over many rotation periods is consistent with a single exponential probability distribution, $p(\Delta t_{{\rm{next}}}; \tau _0)$ $\propto \frac{1}{\tau _0} e^{-\Delta t_{{\rm{next}}}/\tau _0 }$ (where p is the probability density and $\tau _0$ is the average interval between consecutive flares). Equivalently, the number of flares occurring in fixed time intervals, $\Delta t$, is Poisson-distributed, $p(N_{{\rm{flares}}}$; $r, \Delta t) = (r \Delta t)^{N_{{\rm{flares}}}} e^{-r \Delta t}/(N_{{\rm{flares}}}!)$ (e.g., Lacy et al. 1976). Stellar flares from active stars thus occur independently and continuously at a constant rate, $r = 1/\tau _0$ (note that some studies use $\lambda$ instead of r for the rate), except possibly on some short time-intervals ($\Delta t \lesssim 30$ min) within complex events (Hawley et al. 2014; Davenport et al. 2014). Correlations between waiting time and flare energy, as may be expected from magnetic energy buildup and release (e.g., Rosner and Vaiana 1978) in a single region on the star, are not clearly consistent with the stellar data (Hawley et al. 2014). Efforts to understand the flare occurrence as a function of light-curve phase are ongoing in the era of TESS data (Ikuta et al. 2023).

Davenport et al. (2014) constructed a canonical flare template from the normalized, classical (single-peaked) flares of GJ 1243 in short-cadence Kepler data. They modeled the decay phase with a sum of two exponentials and the rise phase with a polynomial. The template reflects a remarkable self-similarity in many flares on 1-min integrations. Further analysis in this work demonstrated that many complex flares can be decomposed into linear superpositions of several templates with adjustable parameters. They then injected synthetic flares over a long time baseline and simulated flare statistics, concluding that the number of complex flares observed in the data is not fully consistent with a chance superposition of many uncorrelated classical flares. This is tantalizing, quantitative evidence for a solar-like “sympathetic” (e.g., Lynch et al. 2016) flaring property of stars.

The rates of flares in the high-energy regime have been investigated in detail with Kepler data. Silverberg et al. (2016) discuss a deviation from a single power-law between $E_{{\rm{Kp}}} = 10^{32.5}{-}10^{34}$ erg flares in the GJ 1243 sample. (In the following, flare energies now refer to the energies integrated over a $T=9000$ K blackbody function, as reported in the respective studies, rather than an energy over a particular bandpass). The upper limit flare energy for white-light flares has been pretty well established from the long monitoring times from Kepler to be $E_{\max } \approx 10^{37}$ erg for subgiants, $E_{\max } \approx 10^{36}$ erg for rapidly rotating G-dwarfs, and $E_{\max } \approx 10^{35}$ erg for slower rotating G dwarfs (Notsu et al. 2019; Okamoto et al. 2021). For G, K, and M stars in the Kepler field, Candelaresi et al. (2014) discuss a correlation between a star’s superflare ($E> 5\times 10^{34}$ erg) rate and its light-curve amplitude modulation due to starspots. The ground-based Evryscope network has also provided new insights into $g^{\prime }$-band superflare rates: Howard et al. (2019) find that the rate of $E \ge 10^{33}$ erg superflares (by extrapolating the best-fit power-laws down over 0.5–1.5 orders of magnitude in energy for the K5–M2 stars) decrease from the averaged active K5–K7 stars to the averaged active M4 stars. Similar to many previous studies (Sect. 5.1), they report a larger average flare energy for earlier-type stars, which they attribute to the “size of the stellar convective region”. Pettersen (1989) review and analyze the time-averaged luminosity due to flaring (and quiescent magnetic activity proxies) in a sample of dKe and dMe stars, and they correlate it with the convective envelope volume. Recent stellar evolution models of dM flare stars in eclipsing binary systems (Feiden and Chaboyer 2013, 2014) and dK/dG stars (Matt et al. 2011) also show that early type stars, which have larger mean flare energies, have larger convective zone volumes.

5.3 Timing on long scales: flare rates and stellar age

By characterizing the FFDs across many types of stars in different evolutionary stages, astronomers begin to piece together how flare rates change over billions of years. With age-dating techniques such as open cluster membership/association and gyrochronology, quantitative evolutionary tracks in the flare history of stars are possible. At the youngest stellar ages, high flare rates and energies have been reported from PMS stars in nearby open clusters, such as Orion, in the X-ray and optical (Getman et al. 2005; Jackman et al. 2020). Signatures of accretion (e.g., in T Tauri systems) indicate a very young age range for some flare stars (e.g., Tofflemire et al. 2017). Membership among moving groups, such as the well-known $\beta$ Pic moving group (with the flare star members AU Mic and AT Mic AB) probe ages around $\sim$ 20 Myr for stars that are no longer actively accreting gaseous material from a circumstellar disk. At older ages ($\gtrsim 100$ Myr), flare rates and properties have been reported in open clusters, such as the Pleiades (Stelzer et al. 2000; Ilin et al. 2021). At the oldest ages, globular cluster and the Galactic bulge membership facilitates calibration of the flare-rate–age relationship, but such stars are located at very large distances. As stars spin down due to angular momentum loss over their main sequence evolution, flare rates are expected to correspondingly decrease. This correlation is supported by detections of stellar superflares from main sequence, rapidly rotating (periods < 5 d) G-type stars. However, binarity plays an important role in maintaining flaring among a population stars to old ages, compared to the Galactic disk populations. This may occur through a co-evolutionary process in M dwarf-white dwarf binaries (Morgan et al. 2016) or rotational tidal synchronization in RS CVn-like systems (Osten et al. 2012). Binarity may also directly affect flare occurrences and energies (e.g., Doyle and Mathioudakis 1990; Gao et al. 2008; Adams et al. 2011, and see references in Sect. 7.5), but it is in general a difficult property to assess without adaptive optics (e.g., Clarke et al. 2018) and high-resolution spectroscopy (e.g., Notsu et al. 2013a), which are observational techniques that are limited to nearby stars.

Following the statistical inference of quiescent magnetic activity lifetime analyses of West et al. (2008), the vertical distances above or below the Galactic plane have been used as proxies of stellar age for M dwarf stars in the field. The basic idea is that over time, stars experience repeated gravitational scatterings in the midplane of the Galaxy and are consequently more likely to be found at larger vertical displacements. Large time-domain surveys, such as the SDSS and soon the LSST, allow characterization of flare rates in multiple populations of stars in sparsely sampled data sets. The stellar flaring fraction and the flare fraction (duty cycle) are two quantities that are useful in mapping the flare rate of a galaxy and the star’s vertical distance from the galactic plane. A stellar flaring fraction is the fraction of stars that produce some number of flares, while flaring fraction, or duty cycle, is the fraction of sparsely sampled epochs that flare. UV and optical surveys have been leveraged to map these quantities with increasing vertical distances (Welsh et al. 2007; Kowalski et al. 2009; Hilton et al. 2010; Walkowicz et al. 2011). There is suggestive evidence that flaring fraction from sparsely sampled data indicates that flaring occurrence among a population decreases faster with stellar age than other measures of quiescent magnetic activity (Kowalski et al. 2009; Hilton et al. 2010). One may be tempted to speculate that this is one possible source for the spread in flare rates among the mid-type dMe stars in Fig. 7, but more investigation is needed. The effects of solar-like activity cycles (Crosby et al. 1993) might also sensibly contribute to the variation of flare rates within a population of stars of the same spectral type (but see Pettersen et al. 1984; Davenport et al. 2020).

Most recently, the improved methods and data samples for flare statistics, asteroseismology, gyrochronology, and open cluster calibrations using Kepler and K2 data have facilitated much more quantitative characterization of the flare rate as a function of stellar age (e.g., Ilin et al. 2021). Davenport et al. (2019) used Kepler field stars and gyrochronology (Mamajek and Hillenbrand 2008) to present the first stellar flaring fraction as a function of age and $g-i$ color, which maps to spectral type (Covey et al. 2007). Davenport et al. (2019) used MCMC non-linear least squares to constrain $\alpha (t)$ in an equation similar to Eq. (11) over a range of stellar masses, from G to M dwarfs, and ages older than $t \ge 10$ Myr. They found that $\alpha (t)$ remains approximately constant, flare rates indeed decrease over time for the M dwarfs, but that the flare rate decreases much less than for the G dwarfs. Notsu et al. (2019) and Okamoto et al. (2021) characterized the evolution of differential flare frequency of G-type stars in the Kepler field by splitting up into different rotational period bins from $<5$ d to 20–40 d. They further divided the $< 5$ d sample into smaller bins and found that flare frequency is approximately constant for periods < 3 d, for which the corresponding ages are $\lesssim 500$ Myr but are not as well constrained from gyrochronology calibrations. Johnstone et al. (2021) combine several state-of-the-art techniques to calculate the quiescent $L_X$ as a function of age and mass, and then they leverage a scaling relation (discussed in Sect. 5.1) from Audard et al. (2000) to map the rate of high-energy XEUV flares on GKM stars over ages spanning 2 Myr to 5 Gyr.

6 Light curve analyses

In addition to integrated flare energy in a particular bandpass (Sect. 5), a number of other strictly empirical quantities can be calculated from single-bandpass data with a moderately fast cadence. A useful decay timescale measurement is the time it takes for a light curve to reach 1/e of its maximum flux (e.g., Maehara et al. 2015; Namekata et al. 2017), which prevents ambiguities in determining the bona-fide end of a flare. This measure also avoids model-dependent uncertainties in fitting an exponential function to flares with events in the decay phase, which may also show lengthening timescales (Osten et al. 2005; Davenport et al. 2014). Here, we summarize a few relationships between decay times and other quantities. The total rise and total decay times of flares are not well-correlated (Moffett 1974; Lacy et al. 1976; Hawley et al. 2014), which possibly suggests some loss of “memory” about the details of the rise phase after the peak. Integrated energy and total duration exhibit a tighter correlation. Larger peak-flare amplitudes typically correspond to longer decay times and larger energy flares, but there is significant scatter. Namekata et al. (2017) find that the 1/e decay duration scales with total white-light energy to the 1/3 power and magnetic field to the $-5/3$ power for flares on the Sun, rapidly rotating G stars (see also Maehara et al. 2015), and M stars (see also Howard et al. 2019). Kővári et al. (2020) found that the 1/3 power scaling between duration and energy also extends to a late-type giant. For similar-duration ($\approx 5-10$ min) solar flares and white-light flares on solar-type stars, the energies of the latter are several orders of magnitude larger (cf. Figs. 8–9 of Namekata et al. 2017).

6.1 Classification of flares by optical broadband evolution

Several classification schemes of flares have been developed based on optical broadband durations and light curve shapes (“morphologies”). Bopp and Moffett (1973) and Moffett (1974) established some of the basic descriptive terminology of flare classification, such as complex, spike, slow, typical, or multipeak (a type of complex flare). More recently, a condensed scheme denotes flares as either classical or complex (Davenport et al. 2014); several classical and complex events are evident throughout a small section of the Kepler data of GJ 1243 in Fig. 8. Kowalski et al. (2019b) argued that flares with several events that comprise the rise phase can appear as relatively simple, “classical” flares if degraded to 1-min time-integrations. A classification of solar flares into “gradual”, “impulsive”, and “secondary” was proposed by Cliver et al. (1986) based on extremely high-time resolution, hard X-ray data. Hard X-ray data is rarely available for stellar flares, and impulsive and gradual classifications have thus been determined from broadband U- and optical light curves.

Fig. 8

A randomly chosen section of the 1-minute cadence Kepler light curve of GJ 1243 from Davenport et al. (2014) showing a variety of classical and complex flares superimposed on the modulation, which is due to starspots (Davenport et al. 2015). The data were obtained from https://github.com/jradavenport/GJ1243-Flares. Note, there was a single flare event (not shown) that was observed simultaneously in both Kepler and in the U-band (Hawley et al. 2014). For reference, this flare had a relative peak contrast of $\Delta F/F_{{\rm{median}}} = 0.014$ in the Kepler band, a peak contrast of $I_{f} = 0.9$ in the U band, and a total U-band energy of $10^{31}$ erg

A common classification utilizes a measure of the impulsiveness (also referred to as “impulse” or “impulsivity”), of a flare. Hawley and Pettersen (1991) quantified the impulsiveness as the fraction of energy emitted in the impulsive phase, which ends after the fast-decay phase(s) in broadband optical radiation. They found that about 70% of the energy was emitted in the impulsive phase during two²⁹ flares that are separated by $\sim 2$ orders of magnitude in total energy. These were both considered very impulsive events. Kowalski et al. (2011) classified three events as “gradual”, “impulsive”, and “traditional” in ultra-high-cadence (0.16 s) photometry data. Kowalski et al. (2013) calculated the impulsiveness index (${\mathcal {I}}$) as the peak value of $I_f$ divided by the FWHM ($t_{1/2}$) of the broadband (U) light curve to classify events into “impulsive flares (IF)”, “hybrid flares” (HF), and “gradual flares” (GF). Examples of high-energy, $E_U > 10^{32}$ erg, GF and IF events are shown in Fig. 9 (top). Lower-energy and lower-amplitude IF, HF, and GF events are shown on the same axes in Fig. 9 (bottom). Dal and Evren (2010) divide a large sample of dMe flares into fast and slow flares and suggest the differences are related to the position on the stellar disk. This hypothesis could be verified with spectral properties that are thought to change with position on the solar disk (Neidig et al. 1993b).

Fig. 9

Representative high-cadence optical light curves of M Ve flares. (Top left) The impulsive and gradual decay phases of a high-energy gradual flare (GF) event on EV Lac; reproduced with permission from Kowalski et al. (2013). The impulsive phase and gradual decay phases are indicated. At the peak (vertical black dashed line), GF-type events like this one have large Balmer jumps, which are still much smaller than predicted by optically thin hydrogen recombination theory alone (e.g., Sect. 7.2.1). The $t_{1/2} = 13$ min is the FWHM of the light curve, and the $U$-band impulsiveness index, ${\mathcal {I}} = I_f/t_{1/2}$, is 0.5 for this flare. (Top right) An example high-energy IF event ($t_{1/2} = 2.5$ min, ${\mathcal {I}} = 8.3$) on EQ Peg A, reproduced with permission from Kowalski et al. (2013). (Bottom) Examples of lower energy flares within the flare IF/HF/GF classification scheme based on the impulsiveness in the U-band or a narrow-band filter (NBF) with $\lambda _c = 3500$ Å. The HF and IF events have similar peak contrasts ($I_f = 2.7{-}3.7$) but differ significantly in the values of $t_{1/2}$ (14 s and 120 s). The IF and GF events are referred to as IF4 and GF2 in Kowalski et al. (2016); the HF event is referred as the HST-2 flare in Kowalski et al. (2019b). Note, the flare-only impulsive-phase spectrum of the IF event is shown in comparison to common photometry filters in Appendix A

The Kowalski et al. (2013) IF/HF/GF impulsiveness index is primarily a diagnostic of the broadband evolution during the impulsive phase of a flare. Thus, it is smaller for flares with a gradual decay phase that starts at a flux level that is a significant fraction of the peak flux. Among the IF-type events, there is a large range, 30–70%, in the percentage of the total energy that is emitted in the impulsive phase. Several spectral quantities (Balmer jump ratios, Balmer line-to-continuum ratios; Sect. 7.2) at the time of peak broadband flux correlate with this index; one possible interpretation is that the average impulsive-phase heating rate varies among events. If the heating rate is attributed to the flux of electron beams, then the inter-flare variation of these optical spectral quantities suggests (Kowalski et al. 2016) an important connection between large-scale flare evolution and electron acceleration, which is sometimes discussed in the context of nonthermal hard X-ray properties of solar flares (McClymont and Canfield 1986; Holman et al. 2011). This expands upon the foundational paradigm that stellar flare scaling relations (e.g., H$\gamma$ v. U-band energies), which extend over many orders of magnitude, mean that differences from flare to flare primarily result from the flaring surface area (e.g., Butler et al. 1988; Hawley and Pettersen 1991). The empirical differences in IF/HF/GF spectra thus eventually motivated the semi-empirical modeling work in Kowalski (2022), which suggests that differences among flares further result from scaling the relative areas of higher and lower concentrations of particle beams in the stellar atmosphere. Recent analyses of H$\alpha$ line profiles from different locations within solar flare ribbons have begun to quantitatively explore the implications for spatially unresolved stellar data (Namekata et al. 2022a).

The fraction of flares within each of the IF, HF, and GF types is not yet well-established. The reported numbers of flares that are IF/HF/GF (Kowalski et al. 2013) are biased toward the IF type because the larger signal-to-noise at peak facilitates spectral characterization. Note that rigorous distinctions cannot be made between an HF and GF classification for certain events (Kowalski et al. 2019b). It is also not clear how many IF events produce larger Balmer jumps (Kowalski et al. 2016) than expected according to the original Kowalski et al. (2013) IF/HF/GF classification scheme. GF events in the decay phase of highly-impulsive events also have shown some of the smallest Balmer jumps (in absorption).

Nonetheless, it is rather clear that the most gradual-type events in the light-curve evolution exhibit the largest Balmer jumps ratios and ratios of H$\gamma$/C4170$^{\prime }$ at the peak time in broadband optical flux, whereas the most impulsive-type events show the smallest Balmer jumps in emission and the largest fraction of optical energy that is radiated in the continuum (Sect. 7.1). This is possibly an interesting connection to the Type I and Type II white-light flare classification on the Sun (Sect. 7.2), but the stellar events that would, ostensibly, be most analogous to Type II white-light flares in the spectral properties are actually the stellar events that are most impulsive.

We should clarify that an impulsive-type flare exhibits an impulsive phase and a distinct gradual decay phase, while a gradual-type flare can show a gradual decay phase that is clearly distinct from its impulsive phase. Figure 8 (bottom 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 (Vievering et al. 2023), 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 Broomhall et al. (2019a) and Zimovets et al. (2021).

QPPs in stellar flares are reported on a variety of timescales, from 10 s to 30 min, and they occur in the impulsive and gradual phases. Mathioudakis et al. (2006) report a 10 s period during the u-band flare in the peak phase using a wavelet analysis (Torrence and Compo 1998) and discuss different interpretations, including one in which the reconnection and particle acceleration are modulated through MHD oscillations in a nearby large loop (Nakariakov et al. 2006). Anfinogentov et al. (2013) detect a damped 32 min oscillation during the decay phase of a U-band megaflare (Kowalski et al. 2010) on YZ CMi using Lomb–Scargle and auto-correlation analyses of the detrended light curve. Doyle et al. (2022) discuss shorter, 1–4 min, 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 (Inglis et al. 2015, 2016) is another powerful technique for evaluating the statistical significance of QPPs. Inglis et al. (2015) apply the AFINO method to the stellar event in Anfinogentov et al. (2013) and Kowalski et al. (2010) 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. (2016) 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 Balona et al. (2015). Mathioudakis et al. (2003) report remarkable oscillations with a 240 s period during the peak phase of a U-band flare on the RS CVn II Peg, and Mitra-Kraev et al. (2005b) 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 Broomhall et al. (2019b). QPPs of 320 s and 660 s periods were reported over the peak phase of a remarkably energetic flare on a PMS M3 star (Jackman et al. 2019b).

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 (Sect. 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: thermal (T~10⁴ K) response of impulsive footpoint heating

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. 2013 reproduced in the top-left panel of Fig. 9 and the event studied in Hawley et al. 1995). The optical and NUV radiation is thought to be thermal radiation that originates from photospheric and/or chromospheric heights. For reference, Table 2 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 2 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 Eq. (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 2) is shown in Fig. 10. The temporal sequence in the figure is representative of many stellar flares, although quantitative differences among events occur (Houdebine et al. 1991; Houdebine 2003; Kowalski et al. 2013, 2019b). 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 (Hawley and Pettersen 1991; Houdebine et al. 1993a; Houdebine 2003; García-Alvarez et al. 2002; Kowalski et al. 2013). The blue continuum flux is the fastest to rise to its maximum and then begin its decay, followed by the highest order Balmer lines (Hawley and Pettersen 1991). The Balmer H$\alpha$ line is the slowest to decay relative to its peak (see also Namekata et al. 2020). 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., Hawley and Pettersen 1991; García-Alvarez et al. 2002, and see also the IF4/F2 event in Kowalski et al. 2016). The neutral helium lines exhibit a very wide range of decay timescales: The blue He I lines decay quickly (Hawley and Pettersen 1991), but He I 10830 Å (studied in a different MVe flare) exhibits one of the most gradually evolving light curves (Schmidt et al. 2012).

Quantitative assessments of the time-scales of emission line and continuum fluxes have been investigated with a time-decrement (Kowalski et al. 2013, 2019b). 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. 10 consistent with a simple model consisting of a slowly decreasing average temperature over the flare footpoints on the star (Gurzadian 1984; Houdebine et al. 1991)? This explanation would be analogous to the cooling thermal loop sequence in solar flares (Aschwanden and Alexander 2001) 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. 10. 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 (Kunkel 1970; Drake and Ulrich 1980), vertical and transverse spatial inhomogeneities (Kowalski et al. 2015; Kowalski 2022), and non-LTE radiative transfer (Hawley and Fisher 1992; Allred et al. 2006) 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. 10.

Fig. 10

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 2. The emission line fluxes are integrated over wavelength, and all light curves are normalized to their respective peak fluxes. The C4170$^{\prime }$ 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. Image reproduced with permission from Kowalski et al. (2013), copyright by AAS

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. 10 (see also Hawley and Pettersen (1991); García-Alvarez et al. (2002); Kowalski et al. (2019b)). The observed flare continuum radiation from the NUV through the optical is collectively referred to as the white-light³⁰. 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 \ge 4000$ Å, flare-only continuum radiation in the spectra of the impulsive phase of stellar flares is a color temperature matching a $T \approx 8000-14{,}000$ K blackbody (Mochnacki and Zirin 1980; Hawley and Pettersen 1991; Katsova et al. 1991; Paulson et al. 2006; Kowalski et al. 2013; Gizis et al. 2013; Lalitha et al. 2013; Kowalski et al. 2016). 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. 11, which was taken at the peak of the event whose evolution is showcased in Fig. 10. Single-component blackbody function models in the range of $T \approx 9000{-}12{,}000$ K can be fit to this flare spectrum at $\lambda \ge 4000$ Å outside of the major emission lines (Kowalski et al. 2013). 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. 11.

Fig. 11

Flare spectrum at the peak of the large, impulsive-flare (IF) event on YZ CMi in Fig. 10. Several blackbody fits to optical ($\lambda \ge 4000$ Å) continuum regions are shown. The flare spectrum (black) has the pre-flare (dotted) spectrum subtracted. Image reproduced with permission from Kowalski et al. (2013), copyright by AAS

The broadband photometry observations from Hawley and Pettersen (1991) 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 Kowalski et al. (2013) 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 the bottom panel of Fig. 9 which has an energy of $E_U\approx 1.6 \times 10^{31}$ erg (Kowalski et al. 2016). 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. 11 is a highly impulsive event. We refer the reader also to the discussions in Kowalski et al. (2016) and Kowalski (2023) 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 2.

In echelle flare spectra with high-resolving power, Schmitt et al. (2008) and Fuhrmeister et al. (2008) 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. 12a). At lower resolving power, similarly blue spectral trends at wavelengths shortward of the Balmer limit have been noted in Kowalski et al. (2013) 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 \mathrm{\text{\AA }}}/f^{\prime }_{6010 \mathrm{\text{\AA }}} =$C4170$^{\prime }/$C6010$^{\prime }$ (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 (1970) and Moffett (1974), and it was clearly evident and extensively characterized in impulsive-phase flare spectra using modern instrumentation and high-time resolution in Kowalski et al. (2013, 2016). An impulsive phase spectrum during a $E _U \approx 10^{31}$ erg, gradual/hybrid-type flare on GJ 1243 is shown in Fig. 12(b) over a broader spectral range than in panel (a). This is one of two flares analyzed in Kowalski et al. (2019b) 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. 11, the Balmer jumps are evident as smaller flux excesses above the blackbody extrapolations to $\lambda \lesssim 3700$ Å.

Fig. 12

(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. Image reproduced with permission from Fuhrmeister et al. (2008), copyright by ESO. 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 (Kowalski et al. 2019b) 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. Image reproduced with permission from Kowalski et al. (2019b), copyright by AAS

Impulsive-phase continuum flux ratios, which are “colors” if they are converted to magnitude differences, are calculated from flare spectra in Kowalski et al. (2013), Kowalski et al. (2016), and Kowalski et al. (2019b) and are shown in Fig. 13. 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. 13 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$^{\prime }$) 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. 13 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 (Dhillon et al. 2007). 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 Kowalski et al. (2016). The flare color from relative photometry is the ratio of two values of

$$\begin{aligned} I_f(t) \langle f_{q,\lambda } \rangle _T \end{aligned}$$

(12)

in the notation of Eq. (7). Several flares were observed simultaneously with spectra over the impulsive to gradual decay phases, which facilitated unambiguous interpretation of the filter ratios (cf. Fig. 3 in Kowalski et al. 2016). 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. 2016), in line with the various blackbody fits to the peak flare spectrum in Fig. 11.

Fig. 13

Compilation of impulsive-phase continuum flux ratios (color-color diagrams) from spectra (Kowalski et al. 2013, 2016, 2019b). 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 (Hawley and Pettersen 1991) are estimated from extrapolations of blackbody curve fits over the $\lambda = 3800{-}4440$ Å range (Kowalski et al. 2013). The wavelength range of C4170$^{\prime }$ can be seen in the right panel of Fig. 12. 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; Kowalski and Allred 2018) 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 (Kowalski et al. 2016) 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

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., Kowalski et al. 2013; Kowalski 2023). 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 (Kowalski et al. 2013). For equal $\lambda \approx 3615$ Å flare-only flux levels during the fast rise and fast decay phases, the fast decay phase is more than $\Delta T \approx 1000$ K cooler. Mochnacki and Zirin (1980) observed similar trends for optical color temperatures using continuum data during several MVe flares at $\lambda = 4200{-}6900$ Å (e.g., see their Fig. 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 Kowalski et al. (2013); 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 Hawley and Pettersen 1991) and retrospectively to the deviations from blackbody fits to the photometry of that flare (bottom right panel of Fig. 11 in Hawley and Fisher 1992). See Figure 31 of Kowalski et al. (2013) 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$ Å (Kowalski et al. 2019b) 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. 13) that spans the blue-optical and red-optical regimes. On the other hand, Fuhrmeister et al. (2008) report hotter blackbody temperatures over the red-wavelengths (not shown) in the flare in Fig. 12a. 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 (Hawley et al. 2014).

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. 14 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. 14), and Balmer line-to-continuum ratios (namely, the H$\gamma$/C4170$^{\prime }$ ratios—see Appendix B) was established from simultaneous spectral and broadband photometry observations (see also Kowalski et al. 2016, 2019b).

Fig. 14

Energy of the hydrogen Balmer component vs. the optical continuum energy evolution in a large sample of MVe flares. 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 (Kowalski et al. 2016, 2019b) since this result have been consistent with these trends. Image reproduced with permission from Kowalski et al. 2013, copyright by AAS

Stellar flares in the NUV generally exhibit a rapid time evolution (Brasseur et al. 2019). 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 Brasseur et al. (2023). Kowalski et al. (2019b) 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 (Hawley and Pettersen 1991; Robinson et al. 1995). For a review of the timing of the IUE/NUV and IUE/FUV spectral integrations within the AD Leo Great Flare, see Kowalski (2022). 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% (Hawley et al. 2007; Kowalski et al. 2019b), 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. 12(a).

7.2.3 Broadband radiated energy budgets

Empirical radiative energy fractionation from the X-ray to the radio provides model-independent³¹ characterization of stellar flares. The energy in the optical and NUV continuum provides a large contribution to the fractionation. To date, Osten and Wolk (2015) 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 (Hawley and Pettersen 1991). 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. 10). The energy radiated in soft X-rays (0.04–2 keV) is estimated to be $\sim 11$x the H$\gamma$ energy (Butler et al. 1988) and to be comparable to that in the U-band, which is roughly 1/6 of the integrated 1200–8000 Å white-light energy (Hawley and Pettersen 1991; Hawley et al. 1995), and a factor of 1/10 of the bolometric flare energy (Osten and Wolk 2015). Hawley and Pettersen (1991) constrain $E_U = 0.64 E_{2000{-}3260\,\mathrm{\text{\AA }}}$, however the bandpass-averaged energy [erg $\mathrm{\text{\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}}}$ (Hawley et al. 2014). Lacy et al. (1976) determine $E_U = 1.2 E_B = 1.8 E_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 Fig. 4.2 in Kowalski 2012).

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., Hawley et al. 1995) 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. Audard et al. 1999, 2000) is sometimes that the broad-band flare energy is accurately accounted for with a scaled 2-temperature quiescent coronal model (Güdel et al. 1997). 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 ± 10–20 min window (Doyle et al. 1988a; Osten et al. 2005)), down to 1/3 (Hawley et al. 1995), to as large as a factor of $\gtrsim 10$ (Güdel et al. 2002a; Osten et al. 2016). Namekata et al. (2020) 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/3 E_{0.01{-}10 \rm{keV}} = 1/3 E_{{\rm{XEUV}}}$ in Osten and Wolk (2015), 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 Pettersen et al. 1984) corresponding to 1 flare per day with $E_{{\rm{XEUV}}} \gtrsim 3 \times 10^{32}$ erg (Audard et al. 2000) for non-simultaneous flares on AD Leo. Tristan et al. (2023) 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$. Schmitt et al. (2019) extend the comprehensive, multi-wavelength analysis of a $E \approx 10^{34}$ erg flare on the K0Ve star AB Dor from Lalitha et al. (2013); they also report comparable energies released in soft X-rays and in an NUV bandpass. This is consistent with a paradigm that a majority (Osten and Wolk 2015; Kuznetsov and Kolotkov 2021) 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 (Osten et al. 2005). New broadband observations have suggested much larger energy releases than expected³² in the NUV regime relative to the optical (Brasseur et al. 2023), while others are rather well-matched by current RHD models that invoke large heating rates (Osten et al. 2016; Kowalski et al. 2019b). 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 (Robinson et al. 2005) (but see the discussion of bright Galex sources in Sect. 7.3). In the NIR, Davenport et al. (2012) present model predictions for the magnitude enhancements and energies relative to the U band. Tofflemire et al. (2012) 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 Zhilyaev et al. (2007). 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. Hawley et al. (2003) 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 (Kunkel 1970; Mullan 1976; Mullan and Tarter 1977; Cram and Woods 1982; Doyle et al. 1989; Hawley and Fisher 1992; Maas et al. 2022; Namekata et al. 2022a; Brasseur et al. 2023), 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 (Robinson et al. 2005; Getman et al. 2023).

However, there are several important assumptions in broadband colorimetry that warrant spectroscopic verification or calibration (e.g., Kowalski et al. 2016). We discuss a few here. First, a correction to the contribution from emission lines in the optical bandpasses must be considered (Hawley and Pettersen 1991; Hawley and Fisher 1992). 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. 14, Fig. 34 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. Allred et al. (2006) 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! Kowalski et al. (2019b) 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 Hawley et al. (2003) 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 Zhilyaev et al. (2007).

The results from broadband flare colorimetry have largely motivated spectroscopic investigations (Giampapa 1983; Hawley and Pettersen 1991; Kowalski et al. 2013) of flares and detailed modeling (Hawley and Fisher 1992), 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 (Hawley and Pettersen 1991; Phillips et al. 1992; Hawley et al. 2003; Ayres 2015b; Froning et al. 2019; Kowalski 2022; Doyle et al. 2012, 2013). 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. (2003) 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 (Hawley and Pettersen 1991). 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 Phillips et al. (1992) to be consistent with radiative backheating of the temperature minimum region.

Remarkable M-star FUV flare spectra were analyzed in Hawley and Pettersen (1991), Redfield et al. (2002), Loyd et al. (2018b), Froning et al. (2019), MacGregor et al. (2021), and Feinstein et al. (2022). In nearly all cases, the continuum spectra are either very flat (in units of $\hbox {erg cm}^{-2}\, \hbox {s}^{-1} \text{\AA }^{-1}$; see also the spectra compiled in Butler et al. 1981; Bromage et al. 1986; Byrne et al. 1990; Phillips et al. 1992) 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, Hawley and Pettersen (1991) and MacGregor et al. (2021) had simultaneous multi-wavelength observations that were analyzed in detail, and Froning et al. (2019) and Loyd et al. (2018b) discuss redshifted line asymmetries in the FUV. Ayres (2015b) 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. 15); 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 (Fang and Ding 1995).

Fig. 15

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$ ± 25 Å, 1610 ± 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. Image reprodcued with permission from Ayres (2015b), copyright by AAS

Broadband comparisons of FUV and NUV responses during flares have been possible due to the time-tagged photon capability of the Galex mission. Robinson et al. (2005) analyzed a giant flare in both Galex bands, and Welsh et al. (2006) 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 (Morrissey et al. 2007; Fleming et al. 2022). More recently, the gPhoton tool (Million et al. 2016) 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 (Brasseur et al. 2019). 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 (Brasseur et al. 2023).

Potential signatures of proton beams are found in both the FUV and in the gamma rays. Loh et al. (2017) and Song and Paglione (2020) 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 (Orrall and Zirker 1976; Kerr et al. 2023). (Woodgate et al. 1992) 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 (Robinson et al. 1993, 2001; Feinstein et al. 2022).

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 Dulk (1985) 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; White et al. 2003). A large multi-wavelength flare from the dM3.5e star EV Lac was studied in Osten et al. (2005) and a light curve is shown Fig. 16. 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 Osten et al. (2005) 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 (Sect. 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

A remarkable radio flare on the dM3.5e EV Lac with the spectral indices calculated from observations at the two observed frequencies. Image reproduced with permission from Osten et al. (2005), copyright by AAS

Figure 17 shows three radio flares from the comprehensive multi-wavelength study of the RS CVn system HR 1099 (Osten et al. 2004). The data show low polarization during the flares and a similar spectral index evolution patterns among these long-duration flares. Güdel et al. (1996) and Smith et al. (2005) present radio analyses of stellar flares that will be discussed below (Sect. 7.7, Sect. 11) in multi-wavelength contexts.

Fig. 17

Flux (top row), polarization (middle row), and spectral indices (bottom row) during three gyrosynchrotron radio flares from the RS CVn HR 1099. Image reproduced with permission from Osten et al. (2004), copyright by AAS

The emission at lower frequencies results from a plethora of possible phenomena. Decimetric flares³³ can be very luminous and highly circularly polarized indicating coherent emission from beam-plasma instabilities (“plasma emission”) or electron cyclotron maser (Lang et al. 1983; Bastian et al. 1990; Osten and Bastian 2008; Villadsen and Hallinan 2019, see references and discussion in Bastian 1990; Güdel et al. 1996; Osten et al. 2005). Coherent bursts can occur at higher frequencies too, and a remarkable example is the 100% left-hand circularly polarized 30 mJy peak-flux, 3-min burst at 5 GHz discussed in Smith et al. (2005) 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 bursts³⁴ have been reported in very high time-resolution data from Arecibo (Osten and Bastian 2006). A surprising lack of Type II radio bursts has been constrained in the context of extensive ground-based monitoring of optical flaring (Crosley and Osten 2018). Many recent numerical, theoretical, and observational efforts have sought to explain the lack of Type II radio bursts from other stars (Alvarado-Gómez et al. 2018; Mullan and Paudel 2019; Alvarado-Gómez et al. 2020; Wood et al. 2021, see also Cliver et al. 2022a).

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 (MacGregor et al. 2018, 2020). Particularly remarkable is the extremely impulsive flare from Proxima Centauri with multi-wavelength observations (MacGregor et al. 2021). 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 ±20% (see MacGregor et al. 2020). 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, Beskin et al. (2017) 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 (Phillips et al. 1996; Bower et al. 2003; Salter et al. 2008, 2010; Brown and Brown 2006; Adams et al. 2011). In the case of the remarkable periastron event in the binary T Tauri V773 A (Massi et al. 2006), the flare was interpreted as synchrotron radiation. There is a notable observation of an RS CVn event at $\nu \le 100$ GHz showing a break and a rise toward higher frequencies (Beasley and Bastian 1998): 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 (Kaufmann et al. 2004). Colliding magnetospheres at periastron are also attributed to some flare-like optical and NUV variability in the PMS DQ Tau system (Salter et al. 2010; Getman et al. 2023), while the optical variability properties at periastron at other times follow the expectations of accretion alone (Mathieu et al. 1997; Tofflemire et al. 2017).

7.6 Soft X-ray and EUV: loops filled with evaporated, thermal $T>10^7$ K 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. 18.

Fig. 18

X-ray spectra from Chandra during quiescence (top) and flaring (bottom) in the active binary $\sigma ^2$ Coronae Borealis. The authors note the enhanced continuum flux increasing towards shorter wavelengths in the flare spectrum. Image reproduced with permission from Osten et al. (2003), copyright by AAS

An extensive review of stellar X-ray flares and analysis methods (prior to and including 2004) are provided by Güdel (2004) and Güdel (2009). A general review of X-ray spectroscopy of stars is contained in Güdel and Nazé (2009). Since 2004, there have been many important X-ray flare studies (Osten et al. 2005; Robrade and Schmitt 2005; Stelzer et al. 2006; Fuhrmeister et al. 2007; Nordon and Behar 2007, 2008; Pandey and Singh 2008; Huenemoerder et al. 2010; Liefke et al. 2010; Fuhrmeister et al. 2011; Pandey and Singh 2012; Lalitha et al. 2013; Pye et al. 2015; Namekata et al. 2020; Paudel et al. 2021; Karmakar et al. 2022). 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).

7.7 The Neupert effect in stellar flares

The backbone of the solar-stellar flare connection is the Neupert effect (Neupert 1968). 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 (1 s$^2$–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., Dennis and Zarro 1993; McTiernan et al. 1999; Veronig et al. 2002b, 2005; Warmuth et al. 2009). 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., Kane et al. 1985; Neidig 1989; Hudson et al. 1992; Neidig and Kane 1993; Fletcher et al. 2007; Krucker et al. 2011, 2015; Kleint et al. 2016). A remarkable example of the stellar Neupert effect is shown for a flare from Proxima Centauri (Güdel et al. 2002a) in Fig. 19.

Fig. 19

A remarkable example of the Neupert effect in a stellar flare in Proxima Centauri; figure reproduced with permission, from Güdel et al. (2002a). 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 Güdel et al. (2004) and Reale et al. (2004)

A stellar flare Neupert effect was first reported in Hawley et al. (1995) using the correlation between the extreme-ultraviolet (65–190 Å) flux and the integrated U-band energy during a long-duration, 3.5-h rise phase of a flare from the dM3e AD Leo. This was followed shortly thereafter in Güdel et al. (1996) 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 (Güdel et al. 2002b; Osten et al. 2004), dMe stars (Güdel et al. 2002a, 2004; Mitra-Kraev et al. 2005a; Wargelin et al. 2008; Fuhrmeister et al. 2011; Caballero-García et al. 2015; Tristan et al. 2023), dK-stars (Lalitha et al. 2013), T Tauri stars (Audard et al. 2007). There are also notable exceptions to the expected multi-wavelength relations from the Neupert effect in stellar flares (Doyle et al. 1988a; Osten et al. 2005); 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 Xia et al. 2021) has been modeled with a large, low-energy cutoff $E \gtrsim 100$ keV in the nonthermal electron beam (Warmuth et al. 2009), which would reasonably generate small amounts of collisional heating within and chromospheric evaporation into the corona (however, the explanation is still being investigated; Alaoui and Holman 2017).

The Neupert effect in Flare D in Güdel et al. (1996) was investigated in detail and was compared to a similar multi-wavelength relationship in a solar flare from Dennis and Zarro (1993). 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 McTiernan et al. 1999). Güdel et al. (1996) 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:

$$\begin{aligned} \frac{d}{dt} \Big ( 3 n_e k_B T V \Big ) = \alpha f_{{\rm{radio}}}(t) - n_e^2 V \psi (T) \end{aligned}$$

(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 cm$^{3}$ s$^{-1}$. The left-hand side of Eq. (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 Fisher and Hawley 1990). 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 Güdel et al. 1996).

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; Hawley and Fisher 1992). Kowalski et al. (2013) 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 (Laming 2004) 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 (Güdel et al. 1999), as have relative changes in high-FIP abundances such as Ne (Osten et al. 2005). An example of an abundance analysis of an RS CVn flare in Fig. 18 is shown in Fig. 20 where an abundance change in Fe was seen, but other expected FIP dependencies were not. Recently, Paudel et al. (2021) 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 (Warren 2014), but Laming (2021) 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.

Fig. 20

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 ^22$ Coronae Borealis. They note that there is a slight enhancement of Mg, Si, O, and Ne during the flare, and Fe increases by nearly twice. Image reproduced with permission from Osten et al. (2003), copyright by AAS

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., Cram and Woods 1982, 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., Allred et al. 2006). 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., Kowalski et al. 2015). 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 (Hawley et al. 1995, 2003; Zhilyaev et al. 2007; Kowalski et al. 2013), calculating LTE continuum spectra for a range of optical depths $\tau (\lambda ) = [0,\infty )$ (Kunkel 1970; Eason et al. 1992), and modeling the NLTE Balmer line decrements, optical depths, and continuum flux ratios (Drake and Ulrich 1980; Jevremovic et al. 1998a; García-Alvarez et al. 2002; Morchenko et al. 2015). The pioneering study of Kunkel (1970) 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. (2006) 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 (1970) 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 flares.³⁵

The higher quality spectral observations of the Great Flare of AD Leo (Hawley and Pettersen 1991) 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 (Hawley and Fisher 1992), 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 1992; Osten and Wolk 2015, 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 (Kowalski et al. 2013). 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. Cram and Woods (1982) 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 (Houdebine 1992; Mauas and Falchi 1996; Jevremovic et al. 1998b; Christian et al. 2003; Fuhrmeister et al. 2010; Schmidt et al. 2012). An example of this modeling technique from Fuhrmeister et al. (2010) is shown in Fig. 21, where the chromospheric temperature adjustments result in satisfactory fits to the emission lines in the giant CN Leo flare (Fig. 12(a)) from Fuhrmeister et al. (2008). Further, a depth-dependent filling factor and missing photospheric heating were inferred through this technique.

Fig. 21

Demonstration of semi-empirical temperature adjustments used to model the red- and blue-wavelength optical spectra of a giant flare from CN Leo. See also Fuhrmeister et al. (2011) and Fuhrmeister et al. (2005) for similar approaches for modeling flares from Proxima Centauri and the dM6 star AZ Cnc, respectively. Image reproduced with permission from Fuhrmeister et al. (2010), copyright by ESO

How to create these deep transition region locations and chromospheres self-consistently with a flare heating function, $Q_{{\rm{flare}}}$? One early hypothesis (Mullan and Tarter 1977) 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). Hawley and Fisher (1992) 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 (Hawley and Pettersen 1991), 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 (Hawley and Fisher 1994; Abbett 1998; Allred et al. 2005, 2006, 2015) 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 (Nizamov 2019).

Significant advances were made in modeling the temporal evolution of the atmospheric response in gas-dynamic simulations, starting in the 1980s (Livshits et al. 1981; McClymont and Canfield 1983; Fisher et al. 1985c, b, a; Emslie et al. 1992). An accurate treatment of energy deposition in a partially ionized, nonuniform chromosphere due to an injected nonthermal electron distribution was developed (Hawley and Fisher 1994) based on analytic formulae for Coulomb collisions (Emslie 1978, 1981a). With the inclusion of radiative-transfer, this effort resulted in a series of papers on radiative³⁶-hydrodynamic modeling of solar flares (Abbett and Hawley 1999; Allred et al. 2005) using a version of the RADYN code (Carlsson and Stein 1992, 1995, 1997) 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. (2006). 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 (Lin et al. 2002) hard X-ray spectroscopic imaging era (Dennis et al. 2022). A double-power law distribution with a low-energy cutoff of $E_c=37$ keV was obtained from the hard X-ray modeling of Holman et al. (2003), 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 (Hawley and Fisher 1992) and semi-empirical models (Fuhrmeister et al. 2010). 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., Hawley et al. 2003), 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 (Kowalski et al. 2010). The lower flux model was run for several minutes in Allred et al. (2006), and it well-reproduced several transition region flare lines from Hawley et al. (2003). 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 Allred et al. (2006) 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. (1981) 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 Katsova et al. 1997). Gan et al. (1992) 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$ (Ichimoto and Kurokawa 1984; Kowalski et al. 2022; Namekata et al. 2022a) and in more optically thin lines such as Fe II in IRIS NUV spectra (Kowalski et al. 2017b; Graham et al. 2020). 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 (Kowalski et al. 2015) and the RH code (Uitenbroek 2001). (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. (2006) 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 Kowalski et al. (2015) is shown in Fig. 22. 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) (Hudson 1972; Ricchiazzi and Canfield 1983; Aboudarham and Henoux 1986; Fang et al. 1993) 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 (Zharkova and Kobylinskii 1993) do not seem to predict the observed blue continuum properties and Balmer jump strengths (e.g., Fig. 11).

Fig. 22

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 (Kowalski et al. 2015). (Top panel) The cumulative of the contribution function to the emergent intensity, $I(\lambda =3550 \mathrm{\text{\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 in the online supplementary material

The formation of the optical and NUV continuum radiation in a F13 electron beam model is illustrated in Fig. 22. 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 Kowalski et al. 2016; Kowalski 2016, 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. Kowalski and Allred (2018) reported self-similar patterns in chromospheric condensation models (using a basic HI³⁷ 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 Kowalski and Allred (2018) is shown in Fig. 13 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.

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., Kowalski et al. 2016, 2017b). 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 (Kowalski et al. 2019b), and success has been achieved in an even broader-wavelength comparisons to a superflare event (Osten et al. 2016). 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}$ $\hbox {erg cm}^{-2}\,\hbox {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 Kowalski et al. (2017b) and Kowalski (2022). 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. Namekata et al. (2020) 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 (Fig. 23). The models with energy transport through only thermal conduction were largely discrepant.

Fig. 23

(Top) Model spectra of H$\alpha$ that were calculated with the RADYN code and analyzed in Namekata et al. (2020). The larger electron beam heating rates reproduce the broad, symmetric Balmer lines in an AD Leo flare. (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. Image reproduced with permission from [top] Namekata et al. (2020) and [bottom] from Notsu et al. (2024), copyright by the authors

Similarly to Namekata et al. (2020), the approaches in Kowalski et al. (2017b), Kowalski (2022), and Brasseur et al. (2023) employ semi-empirical forward modeling with the RADYN and RH codes. These studies use large low-energy cutoffs ($E_c \ge 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/[\mathrm{g\ cm}^{-2}] \gtrsim -2$ to gas temperatures hotter than 10,000 K (see also Kowalski 2023). 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., Kontar et al. 2012) that similarly results in a large fraction of the beam heating over $\log _{10} m_c = -2.2$ to $-1.2$ (Kowalski 2023) 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}$ $\hbox {erg cm}^{-2}\,\hbox {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 (Kowalski et al. 2010, 2013). 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. 13 as star symbols. Alternative explanations are summarized in Anfinogentov et al. (2013). Kowalski et al. (2013) 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 (Hawley and Fisher 1992; Aschwanden and Alexander 2001). 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.

Schmitt et al. (2008) used the Palermo-Harvard adaptive regridding code (Peres et al. 1982; Betta et al. 1997) to model the hydrodynamic response of the giant X-ray CN Leo flare in Fig. 12(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. 24. 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 Fuhrmeister et al. 2010, and Fig. 21 here). They successfully reproduced the X-ray luminosity evolution by multiplying the model by an area of $10^{19}$ cm$^{2}$. The model X-ray luminosity is shown in the right panel of Fig. 24. 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 (Reale et al. 2004) and the active giant HR 9024 as well (Testa et al. 2007; Argiroffi et al. 2019). Reale et al. (2004) model the second X-ray flare in the decay of the event in Fig. 19 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.

Fig. 24

(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. Images reproduced with permission from Schmitt et al. (2008), copyright by ESO

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 (whatever cm$^{-3}$), which is a function of space and time, can be written (in laboratory frame coordinates) as

$$\begin{aligned} \frac{\partial w}{\partial t} + \nabla \cdot (w {\textbf{v}}) = 0 \end{aligned}$$

(14)

in three spatial dimensions (if the time-derivative is 0, then the equation leads to the steady-state solution), where ${\textbf{v}}$ is the macroscopic gas velocity vector, $w{\textbf{v}}$ is the flux³⁸ of the quantity in w, and we have excluded a diffusion term, $-D \nabla w$ (Weinberg 2021), 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 Emslie et al. (1992) and Reep et al. (2022).

The conservation equations are solved simultaneously with an adaptive grid equation (Dorfi and Drury 1987) 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 (van Leer 1977, 1979) 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 Kneer and Nakagawa (1976), McClymont and Canfield (1983), Carlsson (1998), Abbett (1998) for general overviews.³⁹ We especially recommend the lucid and comprehensive lecture on the topic in Dorfi (1998).

The Eulerian (fixed space) forms of the conservation equations that are solved in RADYN flare modeling are written in Allred et al. (2015). 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:

$$\begin{aligned} \frac{\partial e}{\partial t} + \frac{\partial (e v_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 \end{aligned}$$

(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 Smith and Auer (1980) and Fisher et al. (1985c). The internal energy density is the sum of thermal and excitation energy densities:

$$\begin{aligned} e = \frac{3}{2} k_B T \sum _s n_s + \sum _{s,i} n_{s,i} \chi _{s,i} \end{aligned}$$

(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., Thorne and Blandford 2017):

$$\begin{aligned} \varepsilon = \rho \left( \frac{1}{2} v_z^2 + e^{\prime } + \Phi \right) \end{aligned}$$

(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

$$\begin{aligned} {\mathfrak {F}} = \rho \left( \frac{1}{2} v_z^2 + h + \Phi \right) v_z \end{aligned}$$

(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 Eq. (15) is readily obtained.⁴⁰

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 2014) within each hydrodynamic time-step. The solution to the intensity at each depth point is integrated over angle using a Gaussian quadrature numerical weighting (Chandrasekhar 1960), 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 Alfvén 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 Emslie (1978) and a unified model calculation presented in Allred et al. (2015). Additional formalism about the theory of energy loss on ionized species and neutrals is contained in Trubnikov (1965), Chambe and Henoux (1979), Emslie (1981a), Fang et al. (1993), Hawley and Fisher (1994), Mott and Massey (1949), Snyder and Scott (1949). 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.5 k_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., Benz 2002). 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 logarithm⁴¹$\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 (Ricchiazzi and Canfield 1983) $\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 Emslie (1978) and Hawley and Fisher (1994); Allred et al. (2015, 2020) 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 Emslie (1978),

$$\begin{aligned} \Big ( \frac{dE}{dz} \Big )_{{\rm{hyd}}}&= \Big ( \frac{dE}{dz} \Big )_{{\rm{hyd, I}}} + \Big ( \frac{dE}{dz} \Big )_{{\rm{hyd, N}}} \end{aligned}$$

(19)

$$\begin{aligned}&= -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}}} \end{aligned}$$

(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}$ erg keV$^{-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 in Eq. (19).

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 cold-target, Coulomb stopping depths of electrons and protons as functions of their initial kinetic energies, KE$_0$, are shown for a model atmosphere of a flaring M dwarf in Fig. 25.

Fig. 25

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 Eq. (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. 22), 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 (Allred et al. 2020)

The heating from a power-law distribution of electrons and protons is analytically solved in Emslie (1978) following Brown (1973) and Lin and Hudson (1976) in an integral from, $Q_{{\rm{beam}}}(m_c)$, which is a function of column mass over an atmosphere with arbitrary, constant ionization fraction. Hawley and Fisher (1994) 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 Kowalski et al. (2015, 2016).

An alternative method is to numerically solve the distribution function, $f(z,\mu ,E,t)$, of the nonthermal electrons (Leach and Petrosian 1981; McTiernan and Petrosian 1990), and evaluate the spatial gradient/divergence of the energy flux. This is the Fokker–Planck solution (Rosenbluth et al. 1957), 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 (Hamilton and Petrosian 1987; Mauas and Gómez 1997; Zharkova and Gordovskyy 2006; Allred et al. 2015). We refer the reader to the comprehensive description of the collisional Boltzmann transport equation in Trubnikov (1965) and the application of the Fokker–Planck solution to flare modeling in Allred et al. (2020). Other helpful references are Thorne and Blandford (2017) and Boyd and Sanderson (2003). 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:

$$\begin{aligned} Q_{{\rm{flare}}}(z) = Q_{{\rm{beam}}}(z) = \frac{d }{dz} \int _{\mu } \int _{E} \mu v E f(z,\mu ,E) dE d\mu \end{aligned}$$

(21)

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

10 Models of chromospheric and transition region emission lines

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 flare.⁴² 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. 1990). 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. 26 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.

Fig. 26

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. (1973) (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 Vidal et al. 1971), 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 (Barklem et al. 2000). The Balmer H$\gamma$ line is also less affected by thermal, $T=10^4$ K ion-dynamic effects (Stehlé 1994): we show an area-normalized (rather than peak-normalized) profile from Stehlé and Hutcheon (1999), 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-like⁴³ 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: Oks and Gershberg (2016) 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). Oks and Gershberg (2016) 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; Anstee and O’Mara 1995; Barklem and O’Mara 1998). 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${\textbf {+}}$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 (Barklem et al. 2000). 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.

• Microturbulence,⁴⁴ 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 Kowalski et al. (2017b) and Zhu et al. (2019), respectively. Non-Gaussian microturbulence is discussed often in the interpretation of solar flare transition region lines (Dudík et al. 2017), 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., Mathioudakis et al. 1999). Osten et al. (2005) 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., Bornmann 1987), 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). Hawley et al. (2007) 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

  $$\begin{aligned} KE = \frac{3}{2} \sigma _{{\rm{los}}}^2 \rho _{{\rm{chrom}}} \Delta z_{{\rm{chrom}}} A_{{\rm{flare}}} \end{aligned}$$

  (22)

  where $A_{{\rm{flare}}}$ is the projected area on the star (Sect. 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., Drake and Ulrich 1980). 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 Wood et al. (1996), Rathore and Carlsson (2015), and Hansteen et al. (2023). Essentially, opacity broadening follows from the Eddington-Barbier relation for a depth-dependent source function over optically thick wavelengths of a line.⁴⁵ 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. 26 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. (2022). Namekata et al. (2020) 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 (Vida et al. 2019). 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 Wu et al. 2022, for an overview). Many stellar flare observations, however, with echelle resolving power have rather long exposure times ($t_{{\rm{exp}}} \approx 15$ min), 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 (Namekata et al. 2020; Wu et al. 2022).

  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. Wu et al. (2022), Namekata et al. (2022b), and Namizaki et al. (2023) use Voigt/Lorentzian profiles to isolate a Gaussian fit to a red-shifted components in large stellar flares. The bottom panel of Fig. 23 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 (Notsu et al. 2024). Eason et al. (1992) and Fuhrmeister et al. (2008) 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. Fuhrmeister et al. (2011) presents detailed line profiles at high-time resolution during flares on Proxima Centauri with VLT/UVES. Vida et al. (2016) and Vida et al. (2019) discuss blueshifts observed across the Balmer series with some temporal resolution. Notsu et al. (2024) 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. 23(bottom)). Honda et al. (2018) report on high-time cadence monitoring of EV Lac and present persistent blue-wing asymmetries through the entire duration of a flare (see also Maehara et al. 2021). Namekata et al. (2022c) discuss a transient blueshifted H$\alpha$ absorption feature, while Namekata et al. (2022b) 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. (2013) analyze the redshift evolution of a broad H$\alpha$ line component in a flare on the K0Ve star AB Dor. Gunn et al. (1994b) 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 (Brosius and Inglis 2018).

  Flux-weighted line-center redshifts in Si IV and C IV during and after the impulsive phase of several dMe flares were reported in Hawley et al. (2003). Chromospheric-line, red-wing emission in H$\alpha$, Na I D, and He I 5876 Å was described in Fuhrmeister et al. (2005). These observations were tested against scaling relations (e.g., Eq. (3)) from analytic models of impulsive-phase chromospheric condensations (Fisher 1989). 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. Wu et al. (2022) 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. 31). Fuhrmeister et al. (2008) and Kanodia et al. (2022) report on red wing enhancements in He I 7065 Å and 10830 Å emission lines, respectively, in two flares. Fuhrmeister et al. (2020) show remarkably broad He I 10830 Å wings during a flare on the dM3e star EV Lac; correspondingly broad H$\alpha$ profiles were analyzed in Fuhrmeister et al. (2018) and consist of a spectrally unresolved red wing satellite component. Loyd et al. (2018b), Froning et al. (2019), and France et al. (2020) report redshifts in other transition region lines, such as C II, from low-mass stellar flares. Linsky and Wood (1994) and Ayres (2015b) 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 (Houdebine 1992; Houdebine et al. 1993b; Fuhrmeister et al. 2008; Kowalski 2022). Extreme blueshifted (Houdebine et al. 1990) and redshifted (Woodgate et al. 1992; Bookbinder et al. 1992) emission has been reported out to several 1000 km s$^{-1}$ in chromospheric and transition region lines. Koller et al. (2021) 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 (Rubio da Costa and Kleint 2017), or equally blueshifted and redshifted (Doyle et al. 1988a). 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$ Å (Hawley and Pettersen 1991; Namekata et al. 2020). 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). Doyle et al. (1988a) and Eason et al. (1992) 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 Vidal et al. (1973) profiles for electron densities of $n_e \approx 10^{15}$ cm$^{-3}$ (see also Houdebine et al. 1993a, b; Gunn et al. 1994a). The better fits with multiple Gaussian components in these studies support the hypothesis of some type⁴⁶ 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$ $\hbox {km s}^{-1}$ in the wings. Montes et al. (1999) 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 (Dimitrijević and Sahal-Bréchot 1992; Dimitrijevic and Sahal-Brechot 1993). Hawley and Pettersen (1991) and García-Alvarez et al. (2002) present much broader Balmer series lines in comparison to the Ca II K emission line (see Fig. 4 of García-Alvarez et al. 2002) in the impulsive phases of two large dMe flares (see also Lalitha et al. 2013).

We show two remarkable examples of the hydrogen broadening in stellar flare spectra in Fig. 27. 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 6 l 699 (Barnard’s star) that was studied in detail and modeled with the RADYN and MULTI codes in Paulson et al. (2006). The bottom panel displays a sequence of spectra during a flare from the dM4.5e YZ CMi at the same times⁴⁷ as shown for the H$\gamma$ spectra in Fig. A.1, top panel, of Vida et al. (2019). 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 France et al. (2020). A detailed analysis of the H9 Balmer line was presented for the large flare on CN Leo in Fuhrmeister et al. (2008). Broadened Paschen lines have been reported in the decay phase of a large flare from the very low mass star VB10 (Kanodia et al. 2022). Similar comparisons are, in principle, possible with the broadening of the high-order Paschen series and the Ca II infrared triplets (Neidig and Wiborg 1984), for which opacity broadening is perhaps less of a concern in comparisons to Ca II H and K.

Fig. 27

Dramatic broadening differences between hydrogen and metallic emission lines, such as the Ca II resonance lines, in flare spectra are considered evidence for the important role of pressure broadening of hydrogen. The pressure broadening of hydrogen constrains the charge densities (and optical depths) in the flare chromosphere. This conclusion is supported by detailed spectral synthesis with, e.g., the RADYN, RH, MULTI codes. The top figure shows these effects during a flare on the dM4 star Gl 699, reproduced from Paulson et al. (2006), and the bottom panel shows the coarse temporal evolution of the broadening differences during a flare (Vida et al. 2019) on the dM4.5e star YZ CMi

10.2.1 Recombination edge modeling

Accurate interpretation of stellar spectra from EUV through NIR wavelengths depends on modeling how the wings of the hydrogen lines overlap and merge near bound-free ionization limits (Hubeny et al. 1994). Accurate modeling, which should lack any sharp edge feature, opens up diagnostic power from the rich spectral region at wavelengths just longward of a recombination limit. For example, the Inglis–Teller relation is (Inglis and Teller 1939; Rutten 2003)

$$\begin{aligned} \log _{10} n_e = 23.2 - 7.5 \log _{10} n_{\max }^{{\rm{Balmer}}} \end{aligned}$$

(23)

which has commonly been used to determine electron density from the last discernible Balmer emission line near the Balmer recombination limit (Kurochka and Maslennikova 1970; Hawley and Pettersen 1991). Besides stellar flares, many types of astrophysical and laboratory phenomena also lack discontinuities in spectra around the expected location of hydrogen recombination edges (Herczeg and Hillenbrand 2008; Esteban et al. 2004; Meigs et al. 2013; Shull et al. 2012). In solar flares, features such as a “blue continuum bump”, a shift of the edge to redder wavelengths, or a completely flat continuum have been reported in observations (Neidig 1983; Donati-Falchi et al. 1985). The complete lack of any type of edge feature, and a lack of a jump in flux across the expected edge location, is ostensibly an indication of blackbody radiation.

However, the lack of an edge is also expected from recombination theory that includes non-ideal, level dissolution. The fundamental question pertaining to stellar flares is thus if there is enough instrumental blending of the wings and/or overlapping, linearly-superposed wing opacities to account for the lack of the observed discontinuities. A remarkable example of stellar flare spectra with no edge feature is shown in Fig. 12(a). Instead, these spectra (and others) show a series of faint, broad hydrogen lines that fade into a continuum. Other spectra, albeit at much lower resolving power, show some noisy features here, however (Hawley and Pettersen 1991). Ostensibly, the conservation of oscillator strength density should result in a smooth transition of observed flux from a pseudo-continuum of a confluence of high-order, highly broadened and blended emission lines to the onset of recombination radiation. A phenomenological description of the transfer of bound-bound opacity to bound-free opacity was developed by Däppen et al. (1987) using the occupational probability ($w_n$) formalism in the non-ideal equation-of-state of Hummer and Mihalas (1988). The “dissolved-level” bound-free Balmer continuum opacity longward of the Balmer edge has been recently incorporated into spectral models of stellar flares (Kowalski et al. 2017b; Kowalski 2022). However, self-consistent modeling of the pressure broadening of the bound-bound Balmer lines (Tremblay and Bergeron 2009) was not also included in the first models (Kowalski et al. 2015). This incorrectly supported the extremely dense chromospheric condensations in 1D RADYN models (Sect. 8.2) as realistic models of flare density enhancements at cool temperatures (Kowalski et al. 2016). Otherwise, the dissolved-level continuum opacities have largely reconciled observations with models by explaining the lack of an edge around $\lambda = 3646{-}3700$ Å and the range of observed pseudo-continuum properties in between the longer-wavelength, high-order Balmer lines at $\lambda = 3700{-}4000$ Å. The non-discontinuous, transition of dissolved-level opacities at $T \approx 10^{4}$ K have apparently reconciled the hydrogen recombination (Balmer jump) theory of the NUV with the blackbody-like theory of the optical continuum radiation.

The physical cause of the continuous, dissolved-level opacity that was parameterized in Däppen et al. (1987) for hydrogen is still not well understood. Experiments with Rydberg atoms support the idea that Landau–Zener transitions (e.g., Zwiebach 2022) at avoided level crossings (Zimmerman et al. 1979; Rubbmark et al. 1981; Pillet et al. 1984; Gallagher et al. 1989) facilitate a cascading ionization process. In this picture, an electron is photo-excited to a high-lying level n that is pressure-broadened to overlap with the broader, next higher-energy level $n+1$, thus effectively “dissolving” the level n. This description apparently explains the Inglis-Teller relationship through dynamic microfield ionization (Hummer and Mihalas 1988). There are actually many types of calculations (Littman et al. 1978; Damburg and Kolosov 1979; Bergeman 1984; Seidel et al. 1995; Benenti et al. 1999; Fisher and Maron 2002, 2003; Cho et al. 2022, and see the review in Sects. 86.2.2–86.2.3 in Drake 2006; Weisheit and Murillo 2006)⁴⁸ that may contribute to “lowering the ionization potential” or “continuum lowering” in a hydrogen plasma with microfield fluctuations and collisional perturbations.⁴⁹

Due to level-dissolution in dense regions of flare atmospheres that are not very optically thick in the Balmer continuum, applying a simple relationship, such as the Inglis–Teller, that connects the last-visible hydrogen line to a single electron density may be misleading. Kowalski et al. (2022) discuss a scenario wherein level-dissolution of hydrogen in a chromospheric condensation can generate an approximately grey opacity window at $\lambda \approx 3700{-}3750$ Å. Thus, the Balmer lines that are in emission in the emergent spectra actually originate from deeper layers with smaller ambient (thermal) electron densities and no apparent gas velocities. Due to the heterogeneous vertical stratification of the flare atmosphere, the high-order hydrogen Balmer series in emergent intensity spectra can be as broad as or less broad than the low-order Balmer lines (H$\gamma$, H$\beta$, H$\alpha$), which are very optically thick in some of the more recent RHD flare models. Some models of stellar flare spectra at the Balmer limit also incorporate spatial and temporal heterogeneity within the flare ribbons (Kowalski 2022). These seek to explain the inconsistencies in large $n_{\max }^{{\rm{Balmer}}} \gtrsim 15$ (Eq. (23)) being observed (indicative of low electron densities, $n_e= 10^{13}{-}10^{14}$ cm$^{-3}$) with the necessity to also produce the observed continuum properties from a location on the star with much larger electron density ($n_e \gtrsim 10^{15}$ cm$^{-3}$; $n_{\max }^{{\rm{Balmer}}} \approx 11$). Note that in occupational probability theory, however, there is not one single maximum upper level for all atoms in a plasma.

Figure 28 directly compares two different modeling techniques within the Balmer jump spectral region. The emergent intensity is calculated in LTE from a static slab with $dl = 400$ km, $T=10^4$ K, and $\rho =10^{-10}$ g cm$^{-3}$ ($n_e = 3.9 \times 10^{13}$ cm$^{-3}$). Donati-Falchi et al. (1985) modeled several solar flare spectra (with high resolving power) using Voigt profiles up to $n_{\max } = 2\times 10^4 T^{0.25} n_e^{-0.25} \approx 65-85$; this is a saddle-point estimate with Debye-screening of the proton field (Mihalas 1978). This method convolves the Balmer edge with the Voigt function corresponding to the line with $n_{\max }$. The “current (2022)” method uses the TB09$+$HM88 profiles (Tremblay and Bergeron 2009) and the dissolved level bound-free continuum opacity. Many other various approximations for the pressure broadening of hydrogen have been used in models of solar and stellar flares. The ambiguities were discussed extensively in Johns-Krull et al. (1997) and were resolved only in recent years. In Fig. 29, we summarize the wide variety of broadening prescriptions for the Balmer H11 line, whose integrated flux is sensitive to the amount of level dissolution in the densities of flare atmosphere models (Kowalski et al. 2015). For a description of each method presented in the figure, we refer the reader to Kowalski et al. (2022) where similar comparisons are shown for the Balmer H$\gamma$ line.

Fig. 28

Direct comparisons of two methods (Donati-Falchi et al. 1985; Kowalski et al. 2017b) for calculating the Balmer jump spectral region in solar and stellar flares. The emergent LTE intensity is calculated from a slab with $n_e = 3.9 \times 10^{13}$ cm$^{-3}$. The method from Donati-Falchi et al. (1985) predicts a “blue continuum bump”, whereas the current treatment with occupational probabilities predicts an extension of the Balmer continuum intensity to wavelengths longward of the “ideal” recombination limit (indicated by an $\infty$ symbol) at $\lambda = 3646$ Å

Fig. 29

A variety of approximations to opacity line profile functions for the Balmer series have been used in stellar flare modeling codes over many decades. This figure summarizes the differences in models of the line profile function, $\phi _{\lambda }$, for the H11 Balmer line (see also Johns-Krull et al. 1997, for a similar figure for the H10 line). The calculations use $n_e = n_p = 3.16 \times 10^{14}$ cm$^{-3}$, or about 10x larger than in Fig. 28, and thermal Doppler broadening at $T = 10^4$ K is included. The “ideal” profiles, $\phi _{VCS}$, are the VCS73 calculations from Lemke (1997), and do not account for line narrowing due to occupational probability / level dissolution. Several Voigt profiles are also plotted for a range of Lorentzian FWHM values, $\Gamma$. Currently, the most accurate theory corresponds to the $\phi _{{\mathrm{TB09+HM88}}}$ (solid black) curve

11 A comprehensive multi-wavelength stellar flare dataset

What fundamental physical parameters can one infer using a multi-wavelength dataset of a stellar (super)flare? A comprehensive analysis of a series of two consecutive superflares from the binary dM4e$+$dM4e system DG CVn was presented in Osten et al. (2016), who summarize the plasma and geometric properties that can be constrained in a flaring atmosphere. This culminates an extensive history of stellar flare research, and we briefly review the salient points to the algorithm.

11.1 Inferring nonthermal electron properties from stellar data

The nonthermal radiation from power-law electrons in stellar flares occurs as hard X-rays ($E \gtrsim 10$ keV) and as incoherent gyrosynchrotron radiation in the radio/microwaves at cm wavelengths. The radiative processes are thought to be similar to solar flares, in which the hard X-rays are emitted through collisional bremsstrahlung in the chromosphere, whereas the radio flux originates from predominantly higher energy electrons ($E \approx 100--300$ keV) trapped in magnetic fields; see the review in White et al. (2011). The distribution of electrons is parameterized in the hard X-rays and radio using two power-law indices, $\delta _x$ and $\delta _r$, respectively. From hard X-rays, the particle flux distribution is

$$\begin{aligned} F(E) = F_o \left( \frac{E}{E_c} \right) ^{-\delta _x} \end{aligned}$$

(24)

where E is the electron kinetic energy in keV, $E_c$ is the low-energy cutoff, and $F_o$ is the differential (specific) electron flux density (el s$^{-1}$ cm$^{-2}$ keV$^{-1}$) around the cutoff energy. Similarly, a nonthermal electron beam density distribution determined from the radio is written as

$$\begin{aligned} n(E) = n_o \left( \frac{E}{E_c} \right) ^{-\delta _r} \end{aligned}$$

(25)

where $n_o$ is the differential number density (el cm$^{-3}$ keV$^{-1}$) of electrons around energy $E_c$. If the same population of electrons is responsible for both X-rays and radio, then $\delta _x \approx \delta _r$ is expected for electrons with velocity v very close to the speed of light c; otherwise, $\delta _x \approx \delta _r - 0.5$. The different forms (i.e.,, flux density vs. number density) of these equations is because the injected power-law flux distribution is inferred through the collisional thick target model of the hard X-ray spectrum, whereas a power-law density distribution is determined through observations at optically thin microwave frequencies (see below).

Note that the integral of Eq. (24) multiplied by E gives the total beam energy flux density (erg cm$^{-2}$ s$^{-1}$), which is an important input for flare modeling (e.g., F10, F13, etc...). Integrating Eq. (25) over energy gives the total beam density, N, where N is the quantity that is used in the gyrosynchrotron expressions in Dulk (1985) if $E_c$ is set to 10 keV. In the studies of Smith et al. (2005), Osten et al. (2016), Dulk and Marsh (1982), Dulk (1985), Eq. (25) is integrated from $E_c$ to $\infty$ to replace our $n_o$ with the prefactor $N(t) \frac{\delta _r-1}{E_c}$ where N(t) is the total electron density (el cm$^{-3}$) above the cutoff (and the same is usually done for Eq. (24)). Equation 24 can be multiplied by a probability density function for the angular distribution, $M(\theta )$, to specify the directivity of particles with respect to the magnetic field (Allred et al. 2020). A highly-collimated Gaussian distribution is usually chosen for $M(\theta )$. Then, a “1.5D” solution to the Fokker–Planck equation (Sect. 9.2) accounts for the changes in the pitch angles of the particles.

The best diagnostic of accelerated particles in stellar flares is optically thin frequencies in the mm or radio. If the observations are in the radio at GHz frequencies, one must ensure that these frequencies are above the peak frequency ($\approx 10$ GHz) so as to be in the optically thin part of the gyrosynchrotron spectrum. An expression from Dulk (1985) converts from a power-law index in the optically thin flux spectrum, $S_{\nu } \propto \nu ^{\alpha _r}$, to the power-law index in the electron distribution, $\delta _r$:

$$\begin{aligned} \alpha _r = 1.22 - 0.90 \delta _r \end{aligned}$$

(26)

The power-law index in the optically thick part of the spectrum is much less sensitive to the power-law index of nonthermal electrons.

Following Smith et al. (2005), the kinetic energy of nonthermal particles can be constrained from the radio spectrum as follows. The power per volume in nonthermal electrons above the assumed low-energy cutoff value, $E_c$, is given by Eq. (9) of Osten et al. (2016). In this equation, the total number of nonthermal electrons, N(t)V(t), over the flaring volume V(t) is an unknown, but it is directly proportional to the optically thin spectral radio luminosity, $L_{\nu }$. By constraining the value of $\delta _r$ in this part of the spectrum and invoking some prior knowledge or assumption about the peak frequency, an integral over the entire radio spectrum can be calculated—including the optically thick part—to give the time-integrated energy of nonthermal electrons in the flare. This leaves one unknown, the magnetic field, B, in the gyrosynchrotron emissivity formula for the optically thin part of the spectrum (see below for further constraints on the magnetic field, but it must be assumed that the B derived from the radio originates from the same B that is derived from the soft X-rays). If the radio or mm data are available at only one frequency, and thus a spectral index cannot be determined from the data, then contours of nonthermal energies and plausible ranges of magnetic field strengths and spectral indices are reported, as in Fig. 30, which is reproduced from Osten et al. (2016) (see also MacGregor et al. 2021, where a spectral index was derived from mm observations).

Fig. 30

A range of possible nonthermal electron beam power-law indices and magnetic field strengths that can be constrained from gyrosynchrotron data of stellar flares. See also Smith et al. (2005) and MacGregor et al. (2021). Contours of constant electron beam kinetic energies add a third constraint. Cross sectional footpoint areas from optical data can constrain the beam flux density for forward modeling inputs (Osten et al. 2010, 2016). Image reproduced with permission from Osten et al. (2016), copyright by AAS

X-ray and optical data can constrain the physical parameters of flares further. From time-resolved X-ray spectra, one obtains the volume emission measure, $\mathcal {VEM}$. The $\mathcal {VEM}$ is $\int n_e^2 dV$, which is the amount of optically thin emitting plasma at a fixed X-ray temperature, $T_X$. Typical values of the $\mathcal {VEM}$ (in M-dwarf flares) are $10^{51}{-}10^{52}$ cm$^{-3}$. In a $\log _{10}$ $T_X$ vs. $\log _{10}\ \sqrt{\mathcal {VEM}}$ diagram, the peak $T_X$ may clearly occur before the peak $\mathcal {VEM}$ (Favata et al. 2000; Aschwanden et al. 2008; Liefke et al. 2010). Over the decay phase evolution of a flare, the slope in this diagram provides one of three inputs for the widely employed flare loop modeling of Reale et al. (1997). Combined with an exponential decay constant from the X-ray light curve and the maximum X-ray temperature, a semi-circular loop’s half-length, l, can be inferred. The loop sizes in giant X-ray flares are typically reported to be on the order of the stellar radius (e.g., Favata et al. 2000; Reale et al. 2004; Osten et al. 2010; Lalitha et al. 2013; Osten et al. 2016). Note, large 50–200 Mm loop lengths are also inferred under some interpretations of QPPs (Sect. 6.2) in optical broadband light curves (Anfinogentov et al. 2013; Doyle et al. 2022).

Optical and/or NUV data provide a way to relax the geometrical assumption of a single, giant monolithic loop (Sect. 12), while also providing constraints on the electron density and the magnetic field strength in the flare corona. The optical and NUV probe the footpoint radiation, from which a dimensionless ratio of $R_{{\rm{foot}}}/(2l)$ can be calculated, where $R_{{\rm{foot}}}$ is an equivalent radius of the area of an optical flare footpoint: $A_{{\rm{flare,opt}}}=2 \pi R_{{\rm{foot}}}^2$. By combining with an alternative expression for the $\mathcal {VEM}$ in terms of these quantities and using the value of the $\mathcal {VEM}$ from X-ray spectral fitting, a representative value of the coronal $n_e$ is readily obtained that is independent of the number of identical loops in a hypothetical stellar arcade (Eq. 25 of Osten et al. 2016). One also can assume that the magnetic pressure must be at least equal to the gas pressure of the evaporated MK plasma for it to be confined to the loops instead of cooling through expansion, as in a scenario analogous to a filament eruption or coronal mass ejection (Cully et al. 1994). The magnetic field confinement condition places a lower-limit on the coronal magnetic field of the flaring loops.

Mullan et al. (2006) homogeneously analyzed a large sample of many sizes of X-ray stellar flares from various spectral types and evolutionary stages using “Haisch’s simplified analysis” (HSA; Haisch 1983), which involve scaling relations (see also van den Oord and Mewe 1989) that infer coronal temperature, electron density, loop length, and confining magnetic field strength from an emission measure and light curve decay timescale. Typical X-ray temperatures derived from multi-component fits to $E=0.5{-}10$ keV spectra are similar to those in solar flares, $T_X \approx 15{-}30$ MK (e.g., Osten et al. 2005; Liefke et al. 2010; Namekata et al. 2020). At the low coronal temperature end, $T \approx 2-15$ MK, electron density variations, or lack thereof, are constrained with X-ray spectra of the f (forbidden) and i (intercombination) components of helium-like triplets (e.g., bottom right spectra in Fig. 19; Güdel et al. 2002a; Osten et al. 2003; Testa et al. 2004; Osten et al. 2005; Liefke et al. 2010). When significantly detected, the flare-enhanced electron densities are reported in the range of $10^{11}-10^{12}$ cm$^{-3}$ at these cooler coronal temperatures.

Osten et al. (2016) summarize the physical parameters in three of the largest dMe flares that have been observed with X-ray spectra (cf. their Table 5). The inferred temperatures are $T \approx 50$, 140, 290 MK; the electron densities are 3 $\times 10^{11}$, 3 $\times 10^{12}$ 3 $\times 10^{11}$ cm$^{-3}$; the confining magnetic field flux densities are 230, 1100, 580 G, respectively; and $\mathcal {VEM}$ ranges between $10^{54}{-}10^{55}$ cm$^{-3}$. Generally similar numbers are found in the largest RS CVn flares, which produce integrated X-ray energies and luminosities one to two orders of magnitude larger (Osten et al. 2007; Karmakar et al. 2023). YSO’s also exhibit such extreme volume emission measures and temperature ($40-80$ MK; Vievering et al. 2019). These superhot temperatures correspond with X-ray luminosities that are close to or in excess of the bolometric luminosity of the star⁵⁰. Note that in an X-ray superflare from the dM3.5e star EV Lac with comparable temperatures ($T_{X, \rm{peak}} \approx 70$ MK) that were derived from the hydrodynamic modeling analysis of Reale et al. (1997), Favata et al. (2000) argue that accounting for the $E > 10^{34}$ erg of energy release during the flare places further constraints on coronal magnetic fields in the $\approx 4$ kG range before the flare.

The peak stellar flare values of $\mathcal {VEM}$ and $T_{X}$ generally compare well with scaling relations over many orders of magnitude, from superflares to solar flares, in the famous “Hertzsprung-Russell (H-R) diagram of solar and stellar flares” (Yokoyama and Shibata 1998). The H-R diagram of solar and stellar flares was first theoretically explained using MHD scaling relations in Shibata and Yokoyama (1999). This diagram was investigated further in Aschwanden et al. (2008), and the tightest correlation was reported in the empirical relationship

$$\begin{aligned} \mathcal {VEM}(T_p) = 10^{50.8} \big ( \frac{T_p}{\mathrm{10\ MK}} \big )^{4.5 \pm 0.4} [\mathrm{cm^{-3}}] \end{aligned}$$

(27)

where $T_p$ is the peak flare X-ray temperature. Their relation for decay time and peak temperature is, in principle, useful for inferring physical parameters if only the decay phase is caught in the observations (Ayres 2015b); however, the empirical scatter about the trends is an order of magnitude in both quantities. Getman et al. (2008) and Osten et al. (2016) suggest that the hottest ($T_X > 100$ MK) stellar flares exhibit a “saturated”-like volume emission measure of $\approx 10^{55}$ cm$^{-3}$ and do not follow a solar extrapolation in the flare H-R diagram. Shibata and Yokoyama (1999) include two additional parameters (either loop length or magnetic field and ambient thermal electron density) in their scaling relation analogue (see also further alternatives summarized in Heinzel and Shibata 2018) to Eq. (27) to explain the apparent offset in volume emission measure from the solar extrapolation. Howard et al. (2022) include a much smaller flare from an active M-dwarf flare in the context of stellar superflares and solar flares. Osten et al. (2016) report on a stellar superflare (their “F2” event) with $T_X \approx 50$ MK that is within a factor of 2.5 in the predicted $\mathcal {VEM}$ given by Eq. (27).

Large stellar flares occasionally trigger Swift’s Burst Alert Telescope (BAT trigger numbers 172969, 310125, 331592, 596958, 625898), which is intended to automatically slew to gamma-ray bursts and other extra-galactic transients.⁵¹ The combination of X-ray data from the BAT and at softer X-ray energies from the Swift’s XRT facilitates classical model hypothesis testing among thermal and nonthermal models out to $E \approx 100$ keV. In superflares from II Peg and DG CVn, nonthermal bremsstrahlung models were found to be equally statistically robust fits as a superhot thermal model fit. However, in the superflare event from DG CVn (Osten et al. 2016) that the BAT spectra triggered on, the $T_X = 290$ MK thermal model was favored over a nonthermal model because the time-evolution of the X-rays was delayed with respect to the optical V-band data (in line with the thermal, empirical Neupert effect; Caballero-García et al. 2015), and a footpoint area from the optical gives (even by stellar flare standards) unrealistic collisional thick target electron beam energy and flux into the lower atmosphere. The collisional relaxation times were estimated to be rather short compared to the (classical) thermal conduction times, thus further supporting the thermal interpretation. Additional evidence for thermal interpretation of superflare X-ray emission from YSO’s is discussed in Vievering et al. (2019).

In soft X-ray flare spectra, features around the Fe 6.4 keV K$\alpha$ emission line and the 6.7 keV Fe XXV thermal emission line complex are also reported. These features have been interpreted as either electron beam excitation (Emslie et al. 1986) or non-LTE fluorescence of a larger-area photospheric ‘halo’ (Osten et al. 2007; Drake et al. 2008; Osten et al. 2010). However, there are some calibration issues in Swift data around these features that have been reported only recently (Pagani et al. 2011). The recalibrations greatly diminish the 6.4 keV feature while enhancing the 6.7 keV complex. Fe K$\alpha$ fluorescence was studied in a flare from an active giant star using a different observatory (Chandra) without instrumental problems and higher resolving power (Testa et al. 2008); see Sect. 12.2.

Table 3 highlights a few important multi-wavelength diagnostics in stellar flares that have been discussed in this review. Of course, the interpretations of these diagnostics are complicated by issues pertaining to heterogeneities of the unresolved nature of stellar flare sources (in time, and in all three spatial directions), and thus are generally model-dependent. In the final section of this review, we discuss assumptions for and inferences of stellar flare geometries and the associated implications for heterogeneity.

Table 3 A summary of flare diagnostics

12 Stellar flare geometries

Stellar flares are spatially unresolved point sources. Yet, from spectral data, model fitting, and sufficient time-resolution, there are many properties of the spatial dimensions of flaring sources that can be inferred. These reveal similarities to solar flare phenomenology at similar wavelengths and strongly reinforce the solar-stellar connectivity of flare physics. In this section, we review the methods and results of stellar flare geometry calculations.

Optical and NUV observations constrain the chromospheric/photospheric flare areas. A blackbody temperature, $T_{{\rm{flare}}}$, and the fraction (“filling factor”) of the visible stellar hemisphere that is flaring, $X_{{\rm{flare}}} = R_{{\rm{flare}}}^2/R_{{\rm{star}}}^2$, can be fit to spectral data or multi-band photometry through the following equation (Hawley et al. 1995, 2003):

$$\begin{aligned} f_{\lambda , \rm{Model}}^{\prime } = \Big (\pi B_{\lambda }(T_{{\rm{flare}}}) - S_{\lambda ,q}\Big ) X_{{\rm{flare}}} \frac{R^2_{{\rm{star}}}}{d^2} \end{aligned}$$

(28)

where $B_{\lambda }(T_{{\rm{flare}}})$ is the Planck function intensity model of the flare, $S_{\lambda ,q}$ is the surface flux of the quiescent star, and $f_{\lambda , \rm{Model}}^{\prime }$ is the ‘flare-only’ model flux at Earth above the atmosphere that is to be fit to the observed, calibrated [$\hbox {erg cm}^{-2}\, \hbox {s}^{-1} \text{\AA }^{-1}$] flare-only flux, $f_{\lambda }^{\prime }$. The ratio of flare-only fluxes at two wavelengths is a flare color and is the minimum information that is needed to solve for $T_{{\rm{flare}}}$ in Eq. (28) (through linearization and iteration or a lookup table). An m-component, linear superposition of RHD surface flux flare spectra, $S_{\lambda , \rm{RHD}, m}$, or a multi-thermal blackbody model are typically required to reproduce the optical and NUV continuum observations over a wavelength range spanning more than $\Delta \lambda \approx 1000$ Å. Additional terms are thus included in Eq. (28). For an m-component RHD model superposition, the solution is a linear least squares fit to an observed flare-only spectrum, $f_{\lambda }^{\prime }$, resulting in the estimated parameters, ${\hat{X}}_{\rm{flare},m}$ (see Osten et al. 2016; Kowalski et al. 2017b; Kowalski 2022). For further discussion, see Appendix C.

Constraints on $X_{{\rm{flare}}}$ from blackbody fitting to optical flare-only spectra result in inferences that a small fraction of the visible stellar hemisphere flares for even the largest, $E \approx 10^{34}$ erg, dMe events. For example, the broadband photometry in the Great Flare of AD Leo (Hawley and Pettersen 1991) constrain $X_{{\rm{flare}}}$ to be 0.5% (Hawley and Fisher 1992) in the first impulsive peak. Fits to the flare-only, optical continuum in the rise/peak spectrum of this event are consistently small (0.2%) but with a correspondingly larger blue-optical blackbody color temperature of $T_{{\rm{flare}}} \approx T_{BB} = 11,600$ K (Kowalski et al. 2013). An $m=2$ component RHD model of this spectrum gives a similar filling factor for a model with large electron-beam heating fluxes, which are combined with the second RHD model component with filling factors that are $\approx 2{-}10$x larger in order to account for most of the narrow emission line flux (Kowalski 2022, see also Cram and Woods 1982; Hawley and Fisher 1992; Kowalski et al. 2010). The corresponding projected area at the star is $A_{{\rm{flare}}} = X_{{\rm{flare}}} \pi R_{{\rm{star}}}^2$, which may instead consist of a series of K circular flare kernels. Inferences in smaller dMe flares from broadband photometry (Hawley et al. 2003) and optical spectroscopy (Kowalski et al. 2016)⁵² are consistent with very small areas $\approx 5{-}10\times 10^{16}$ cm$^{2}$, thus reinforcing the analogy to compact white-light kernels that are observed at the footpoints of impulsively heated loops in solar flares (e.g., Fletcher et al. 2007). During the most energetic dMe events, however, peak-phase flare areas are inferred to be several orders of magnitude larger than solar flare white-light and NUV areas (Krucker et al. 2011; Kowalski et al. 2017a).

In the gradual decay phase, and in the impulsive phase of gradual-type dMe flare events with large Balmer jump ratios, blackbody fits to broadband photometry may be misleading (Sect. 7.2.4). Blackbody fits to the gradual phase continuum give a rather inconsistent picture of the post-peak footpoint areal evolution from event to event; perhaps, heterogeneous RHD model superpositions in the gradual decay phase are much more important for accurate areal inferences than in the $10^4$ K blackbody-like rise and peak phase. A superposition of RHD model components in the gradual decay phase of a megaflare on the dM4.5e star YZ CMi is shown in Fig. 31 from Kowalski et al. (2017b). This model was previously developed to explain the NUV to V-band ratio in the DG CVn superflare event reported in Osten et al. (2016), and it also adequately modeled the NUV continua and Balmer jumps in much smaller, gradual-type events (Kowalski et al. 2019b). This “DG CVn superflare multi-thread (F13) model” consists of several RHD model spectra over the evolution of the atmosphere that is shown in Fig. 22, thus superposing snapshots from a range of optical depths, temperatures, and densities (see Sect. 12.1) in order to explain the gradual decay phase.

Fig. 31

Gradual-decay phase, flare-only spectrum of a megaflare on the MV4.5e star YZ CMi (Kowalski et al. 2010). Multi-thread RHD modeling of another superflare on DG CVn (Osten et al. 2016) adequately reproduces the Balmer jump, optical continuum shape, and some of the hydrogen Balmer emission line properties (but some important shortcomings remain). Each RHD model component is calculated self-consistently, but the linear superposition of the model components is semi-empirical; see discussion in Kowalski et al. (2017b) and Kowalski et al. (2019a)

Various size scales in solar and stellar flares are shown in Fig. 32.

Fig. 32

Sizes to scale. An image of the Sun at $\lambda = 6173$ Å from the Helioseismic and Magnetic Imager (HMI; Scherrer et al. 2012) on SDO, a representative size of an active mid-type (MV3e-MV4.5e) M star ($R_{{\rm{star}}} = 0.3 R_{{\rm{Sun}}}$), and a circular flare area ($\approx 5 \times 10^{18}$ cm$^{2}$ ) that corresponds to about 0.5% of the visible hemisphere of this hypothetical M star. For comparison, Earth and Jupiter are shown to scale. The insets each have a field-of-view of $\approx 10^{19}$ cm$^{2}$, and they show an image of the solar active region in the NW quadrant during the hard X-ray impulsive phase of the SOL2014-03-29T17:48 X1.0 flare (total intensity left, difference image right). For a detailed analysis of the HMI data of this solar flare, we refer the reader to Kleint et al. (2016). Jupiter image credit: NASA, ESA, and J. Nichols; Earth image credit: GOES/NOAA/NESDIS

12.1 A dominant loop or an arcade?

There are two general approaches to modeling spatially unresolved flare light curves and spectra (Liefke et al. 2010). The first approach is summarized in Karmakar et al. (2023), and references therein, where it is most recently applied to an RS CVn X-ray superflare. This modeling follows the analysis of Reale et al. (1997) and describes the rise and decay of a flare as result of a dominant monolithic loop (or a superposition of a few dominant loops that are heated and cooled simultaneously). The heating and cooling timescales in the loop are on the order of tens of minutes to hours (see also Reale et al. 2004; Testa et al. 2007; Argiroffi et al. 2019). In Reale et al. (2004), the initial X-ray event in the same flare in Fig. 19 was modeled with this method, and the subsequent X-ray flare in the event was modeled as an arcade of $\approx$ five loops that were sequentially heated in a nearby region. Coincidentally, a solar analogy inspired a similar hypothetical explanation for the NUV and optical footpoint ignition during the decay phase of a different M dwarf flare in Kowalski et al. 2017b. Longitudinal MHD oscillations in a monolithic, giant (200 Mm) loop have been proposed as an explanation for the triggering of quasi-periodic, U-band, light-curve oscillations in the decay phase of this megaflare event (Anfinogentov et al. 2013). The comparisons of particularly dense snapshots (e.g., $t=2.2$ s; Kowalski et al. 2015) from RHD models to observations also inherently make the assumption that a single dominant kernel at any one time contributes to the bulk of the flux received at Earth (but see Namekata et al. 2022a, for recent tests of this assumption).

An alternative hypothesis is known as “multi-thread modeling”, which was pioneered by Hori et al. (1997), Reeves and Warren (2002), and Warren (2006) to explain Sun-as-a-star soft X-ray solar flare light curves. These works addressed the large discrepancy between the observations of long decay times in GOES soft X-rays and the very short decay times in models of individual loops; multithread techniques have further been used to model a variety of solar flare phenomena (e.g., Reeves et al. 2007; Rubio da Costa et al. 2016; Reep and Toriumi 2017). Elements of this approach and of observations of sequential footpoint brightenings in the hard X-ray and optical wavelengths in solar flare ribbons (e.g., Kosovichev and Zharkova 2001; Qiu et al. 2010) have inspired stellar flare RHD model superpositions (Osten et al. 2016; Kowalski et al. 2017b) for more direct comparisons to spatially unresolved observations. Getman et al. (2011) argue that the rise phase of a giant X-ray flare from the PMS star DQ Tau consists of a multi-threaded source; see also Schmitt et al. (2008) in reference to the modeling approach showcased in Fig. 24. Doyle et al. (2022) suggest that a series of tens of loops contribute to the geometrical dimensions, as inferred from their QPP (Sect. 6.2) analysis, in two high-energy flares. Mathioudakis et al. (2006) discuss the loop triggering scenario of Emslie (1981b) as an alternative hypothesis for short-duration QPPs in white-light. Osten et al. (2005) hypothesize that an arcade of heterogeneous gyrosynchrotron-emitting loops may explain the microwave spectral indices in a large dMe flare.

Multi-thread modeling approaches seek to be more representative of a solar flare arcade geometry, following arguments previously made in stellar (Hawley and Fisher 1992) and solar modeling (Hawley and Fisher 1994). However, ad hoc geometrical assumptions and simplifications must be imposed due to the absence of appropriate (i.e., broad spectral coverage, high-time resolution, unsaturated data at comparable wavelengths) constraints from solar data. Recent analyses of IRIS solar data at high-time resolution, such as the “new area” calculations of Graham et al. (2020), and spatially resolved spectral line analyses (Namekata et al. 2022a), may be able to provide clearer guidance to stellar models of sequential footpoint ignition. Nonetheless, there may be important differences between solar and stellar flares that complicate a direct application of the thermal multi-thread modeling of Warren (2006): shorter-during footpoint pulses (e.g., Mathioudakis et al. 2006), greater importance of the gradual-phase white-light continuum radiation (Hawley and Pettersen 1991; Heinzel and Shibata 2018) and particle acceleration (Osten et al. 2005), and larger apparent loop height growth over the flare evolution ($\tau _c \propto l^2$; e.g., Osten et al. 2005, 2016) all are examples of properties that distinguish stellar flares.

Novel stellar flare applications of the two-ribbon solar flare reconnection model of Kopp and Poletto (1984) are described in Poletto et al. (1988) and Güdel et al. 1999 (see also Haisch 1983). Güdel et al. (1999) model the $\Delta t \approx 8000$ s rising phase of a flare on the RS CVn UX Ari in this framework. This flare exhibited a remarkable peak X-ray luminosity of $>10^{32}$ erg s$^{-1}$, and they infer the evolution of the reconnection heights and electron densities (which evolve through $\approx 10^{11}$ cm$^{-3}$ at peak). A long-duration decay phase is predicted by the model, and the paper discusses how the lack of nonthermal particles in the early rise phase is a potential limitation of the original assumptions.

12.2 Stellar core-halo

Multi-dimensional geometrical analysis of flare sources is possible by combining hydrodynamic modeling of loop lengths and high resolution spectral observations at Fe K$\alpha$. The technique is described in Testa et al. (2008), who infer an X-ray source height that irradiates a large area of the photosphere, from which the Fe K$\alpha$ is observed at a viewing angle $\phi$. The geometry is shown in Fig. 33. The combined constraints from the Palermo-Harvard code and the MOCASSIN 3D Monte Carlo radiation code (Ercolano et al. 2003, 2008; Drake and Ercolano 2007) are reproduced from their paper on the right of Fig. 33. Recently, Werner band H$_{2}$ flare emission has been reported in France et al. (2020) and interpreted as fluorescence due to O VI $\lambda = 1032$ Å irradiation of the temperature minimum region ($T \le 1500$ K).

Fig. 33

Geometrical inference from spectral analysis of the 6.4 keV (1.94 Å) Fe K$\alpha$ emission line during a flare on the giant star HR 9024, reproduced from Testa et al. (2008) with permission. The X-rays from hot flare loops fluoresce the neutral iron in the surrounding photosphere. The modeling leverages hydrodynamic modeling from Testa et al. (2007)

In 1D radiative-hydrodynamic modeling, XEUV backwarming of the photosphere is approximated using the treatments of Gan and Fang (1990), Hawley and Fisher (1992), Hawley and Fisher (1994) and Abbett and Hawley (1999) through exponential integrals and detailed radiative transfer (Allred et al. 2015; Wahlstrom and Carlsson 1994). A multi-dimensional extension of UV and optical radiative backwarming is illustrated in a cartoon of Figure 4 of Fisher et al. (2012), and the geometrical dilution of chromospheric radiative flux impinging on the photosphere is given by their Eq. 25.

The remarkable, state-of-the-art 3D MHD codes that have been recently developed with the primary goal to study quiet Sun magnetism are not ready to address some of the most pressing current problems in stellar flare physics. Flare models require accurate treatments of non-equilibrium physics and detailed radiative transfer for comparisons to chromospheric observations. An approximation to nonthermal particle heating can in principle be included as $Q_{{\rm{flare}}}$ in the energy equation in 3D simulations of small flares (Frogner et al. 2020). The geometrical grid spacing in 3D MHD models, however, is immobile and too coarse (e.g., see Sect. 8.2) to resolve the gradients and contribution function features for the optically thick spectral lines and continua. Section 4.6.1 of Bjørgen et al. (2019) discusses the issues facing a 3D MHD solar flare simulation that predict chromospheric emission lines in response to thermally conducted fluxes. The RADYN code (for example) was originally a “quiet” Sun code, but it had implemented the 1D adaptive grid scheme of Dorfi and Drury (1987) and was flexible enough to simulate particle beam heating in a non-solar gravity. These features proved fortuitous in applications to stellar flare models. The stellar community is in urgent need of collaborative development of multi-dimensional extensions of Carlsson’s and Dorfi’s works, such as Stökl and Dorfi (2007), in order to understand the comprehensive, radiative-hydrodynamic response of stellar atmospheres to the impulsive release of magnetic energy.

13 Conclusions and future outlook

In conclusion, we summarize six, non-exhaustive big-picture questions that are (in our view) at the forefront of stellar flare research.

• Stellar flares routinely attain energies and peak luminosities that are factors of $10^2$–$10^4$ greater than the largest solar flares. How are the flare energies connected to photospheric magnetic flux densities/geometries, particle acceleration properties, flaring footpoint areas, and energy transport processes in the atmospheres of other stars?

• Stellar superflares exceed the quiescent bolometric stellar luminosity. The rate of superflares, however, is not characterized with high statistical significance. What are the rates of events in the superflare regime for solar analogs and other low-mass cool stars as a function of stellar age?

• On the Sun, the brightness of chromospheric flaring sources is transient at any single location for (arguably) $\lesssim 10{-}20$ s. A comprehensive, self-consistent model (e.g., within frameworks similar to Reeves and Forbes (2005) and Rempel et al. (2023) pertaining to solar flare soft X-rays) has not yet been developed for optical or radio stellar flare light curves that last for tens of minutes to hours. What are the principal physical parameters that drive the observed temporal scales of stellar flare light curves of optically thick radiation in the rise and decay phases?

• Blue asymmetries in Balmer lines are often interpreted as evidence of mass eruptions from a star. Are there direct analogues (e.g., as in Mg II; Tei et al. 2018) on the Sun, and does chromospheric evaporation in stellar flares contribute to blueshifted spectral features in cool lines?

• Turbulent mass motions and other non-equilibrium macroscale processes (e.g., UV and X-ray radiative backwarming) in stellar flares remain poorly understood in three spatial dimensions. If gas-dynamic turbulence is important for explaining the observed spectral line shapes, what processes generate large Reynold’s numbers in flare atmospheres and, more specifically, in flare chromospheres?

• Optical and NUV stellar flare spectra are inconsistent with optically thin hydrogen recombination theory. Electron-beam and thermally-conducted heating fluxes that are within our expectations of the standard solar flare model paradigm do not explain the observed strength and spectral shape of the continuum radiation at NUV and optical wavelengths in M-dwarf flares. What is the source of stellar flare white-light continuum radiation, which is prominent in almost all of the Kepler (and TESS) broadband flare data?

Supplementary Information. A movie of Fig. 22 is provided in the supplementary online material of this article.

Appendices

Appendix A Optical filter reference

We include several appendices in this Living Review. This Appendix A serves as a quick-reference for optical broadband filter curves. Appendix B shows supplementary spectral constraints from optical dMe flare data, Appendix C reviews two continuum modeling approaches with slab geometries (blackbody fitting and optically thin hydrogen emissivity spectra) and the continuum modeling approach with parameterized chromospheric condensations (Kowalski and Allred 2018), and Appendix D reviews the salient physics (through a mostly qualitative exposition) of the pressure broadening of the spectral lines of hydrogen.

In Fig. 34, we show the optical broadband and narrowband filters that are commonly used in stellar flare studies. There are numerous databases that provide these, such as the SVO Filter Profile Service at http://svo2.cab.inta-csic.es/theory/fps/, and some are published alongside photometry calibration surveys (Bessell and Murphy 2012). For comparison, we show the low-resolving power, flare-only spectrum (from Fig. 3 of Kowalski et al. 2016) over the impulsive phase of the IF event from the bottom panel of Fig. 9. The spectral shape from the calibrated ULTRACAM narrowband filters at the flare peak time is in agreement with the spectrum.

Fig. 34

Quantum efficiency curves of common optical filters, normalized to peak transmission. An example of a terrestrial atmospheric transmission curve from the Sloan Digital Sky Survey is shown for an airmass of 1.0 at the Apache Point Observatory. One should note that the quantum efficiency of the CCD, telescope primary mirror, and other optical elements are not included. For reference, the impulsive phase flare spectrum from the IF4 event in Kowalski et al. (2016) is shown. The wavelengths that are affected by the dichroic are in a lighter shade of red. The ULTRACAM photometry (NBF3500, NBF4170, NBF6010) at the peak of this flare were calibrated independently, and they are scaled to account for the factor of ten shorter integration time. The flare data were obtained from the Zenodo repository at https://doi.org/10.5281/zenodo.45878

Appendix B A summary of supplemental model constraints in the optical

In this appendix, we supplement the impulsive phase flare colors from spectra, shown in Fig. 13, with the peak colors calculated from narrow-band continuum filters. Figure 35(top) is a summary of the peak-phase ULTRACAM filter (Appendix A) ratios from Kowalski et al. (2016). These flare colors were calculated at higher time resolution than from the spectra that were discussed in the main text (Sect. 7.2). Figure 35(middle) shows the evolution of the ULTRACAM filter ratios in the large IF1 event described in row (1a) of Table 2; see also Kowalski (2023). Average fluxes in the ULTRACAM narrowband filters are calculated from the spectra of the energetic IF3 event (also on YZ CMi) from Kowalski et al. (2013) and are shown for comparison in Fig. 35(middle). Note, the peak flare spectrum for this event is showcased in Fig. 11. In the bottom panel, we show Balmer jump ratios against the ratios of the H$\gamma$ line flux to the flare-only continuum flux around $\lambda = 4170$ Å at the time corresponding to the optical continuum peak flux. (The H$\gamma$ line flux is continuum-subtracted, integrated over the wavelengths of the line, and pre-flare subtracted; C4170$^{\prime }$ is the average the flare-only flux over $\lambda = 4155-4185$ Å). Figure 35(bottom) combines the impulsive-phase spectral measurements from Kowalski et al. (2013) with the same quantities obtained in the four flares (indicated on the top axis) with low-resolving-power, optical spectra that were analyzed in Kowalski et al. (2016) and Kowalski et al. (2019b). It is our hope that as more flares are observed, these figures will be updated to provide compelling and comprehensive constraints for future multi-dimensional flare models.

Fig. 35

(Top) Flare-peak colors from ULTRACAM narrowband filter ratios; similar to Fig. 13 and see text for details. (Middle) Evolution of flare colors in two large flares on YZ CMi. Time increases over the tracks in a clockwise direction. These data have small statistical uncertainties (which are shown) but non-negligible systematic errors. The uncertainties are correlated, which we show for representative 5% systematic absolute color calibration uncertainties in both flare colors located at $(x,y) \approx (2,2)$. The error-ellipses indicate representative joint probabilities: the outer is the joint 68% confidence interval, while the inner translates to the marginal 68% uncertainties. (Bottom) Flare-only, line-to-continuum flux ratios from the same peak spectra in Fig. 13 in the main text

Appendix C Slab model reference

1.1 Blackbody fitting to stellar flare optical data

In this appendix, we summarize the assumptions and techniques employed by fitting isothermal blackbody functions to spectra, broadband photometry, and narrowband photometry of flares. A blackbody “color” temperature refers to the temperature of a blackbody that is fit to a ratio of two fluxes (e.g., a flare color calculated from a ratio of two flare-only continuum fluxes; Sect. 7.2.1), where color is equivalently a difference of magnitudes in two bandpasses. The term also refers to the temperature of a blackbody that is fit to the absolute value of more than two continuum flux values (outside of prominent emission lines) over a limited spectral range ($\Delta \lambda \approx 1000-2000$ Å). Generally, a color temperature contains rather limited information about the actual nature of the radiation, such as whether it is optically thick, thin, an intermediate optical thickness, or which atomic processes contribute to the emergent radiative flux. A color temperature does not unambiguously constrain how the spectrum of continuum radiation extends to wavelengths outside of the range over which it is calculated. Within the total wavelength span that a flare color is calculated (e.g., the ratio of flare fluxes in narrowband filters around 4170 Å and 6010 Å), there may be emission lines that affect the calculated value, and the actual spectrum between two wavelengths may be best represented by a more complex continuum model. Therefore, flare color temperatures are equivalent color temperatures of blackbody functions.

Other measures of temperature include gas temperature, radiation (or brightness) temperature, and effective temperature (in addition to the color temperatures of optically thin recombination and bremsstrahlung spectra). None of these are possible to robustly measure directly from optical stellar flare observations and must be inferred from detailed modeling. In principle, flare effective temperatures could be measured, but there has not yet been a more complete broadband spectral observation of a flare than analyzed in Hawley and Pettersen 1991. To demonstrate the ambiguities in these measures of temperature, it is interesting to note that for an optically thin hydrogen recombination and free-free spectrum calculated from LTE level populations for a slab with a gas temperature of $T_{{\rm{gas}}}= 10{,}000$ K (Appendix C.2, Fig. 38), the radiation temperature at $\lambda = 4170$ Å is $T_{\rm{rad}} = 3600$ K. The emergent intensities are calculated for a slab thickness of $dl = 0.1$ km to ensure that $\tau (\lambda ) < 0.005$ across the optical regime. The color temperature measured from the ratio of intensities at $\lambda = 4170$ Å to 6010 Å is $T_{{\rm{FcolorR}}}=5300$ K. Thus, $T_{{\rm{gas}}}> T_{{\rm{FcolorR}}} > T_{\rm{rad}}$ in this example.

Another interesting ambiguity in inferred gas temperatures occurs when the emergent continuum flux is formed over significant ($\tau \approx 1$) optical depth in dynamic, non-LTE, non-isothermal flare atmospheres. In Fig. 36, we summarize the argument in Kowalski (2023), which is also closely related to the analyses in the appendices of Kowalski et al. (2017b). Two Planck functions with $T = 13{,}000$ K and $T = 14{,}000$ K are shown. If wavelength-dependent, continuum optical depths are large over an inhomogeneous temperature structure, then one could possibly, for example, measure radiation temperatures of $T_{\rm{rad}} = 13{,}000$ K at $\lambda = 4170$ Å and $T_{\rm{rad}} = 14{,}000$ K at $\lambda = 6010$ Å. The ratio of the corresponding emergent intensities gives a color temperature that is only $T_{{\rm{FcolorR}}} = 10,500$ K. Thus, $T_{\rm{rad}} > T_{{\rm{FcolorR}}}$ in this example. For semi-infinite atmospheres that exhibit a linear dependence between the optical depth and the (LTE or non-LTE) source function, the Eddington–Barbier approximation relates the radiation temperature to the gas temperature at the depth corresponding to $\tau (\lambda )=1$; however, dynamic flare atmosphere calculations are usually much more complicated over the large physical depth ranges of continuum formation (e.g., in chromospheric condensations and the stationary layers below; cf. Fig. 2(right) and Fig. 22), which results in differences between the radiation temperature and the gas temperatures over which the emergent continuum intensity forms (see Kowalski 2023, for details).

Fig. 36

Demonstration showing the relationship between a color temperature and wavelength-dependent radiation (brightness) temperatures at $\lambda = 4170$ Å and $\lambda = 6010$ Å

The temperatures and filling factors calculated from optical and NUV flare data are in general obtained through two methods:

• 1. From a ratio of flare-only fluxes in two bandpasses, one can solve for the temperature and filling factor ($T_{{\rm{flare}}}$, $X_{{\rm{flare}}}$) of a Planck function following Hawley et al. 1995 (i.e., by equating ratios of the RHS of Eq. (28) in Sect. 12 to the observed flare-only ratio). The ratio equation is a transcendental equation for $T_{{\rm{flare}}}$ and can be solved numerically using linearization and iteration (i.e., Newton–Raphson) or through a simple lookup table. The result for $T_{{\rm{flare}}}$ is used to infer $X_{{\rm{flare}}}$.

  As discussed in Kleint et al. (2016) and Castellanos Durán and Kleint (2020), by not accounting for the fraction of stellar flux ($X_{{\rm{flare}}}$) that no longer shines through with the quiescent brightness, such as in a blackbody flare (see Hawley et al. 1995) or a flare with a transient overlying absorbing layer, it is possible to calculate a color temperature of $\approx 8500$ K from flare-only irradiance observations of solar flares, as reported in Kretzschmar (2011). In Fig. 37, we demonstrate what happens to the calculated color temperature of an optically thick ($\tau \gg 1$) flare source if one does not use the proper equation in Hawley et al. (1995). For example, a 5800 K Planck function subtracted from a 7800 K Planck function at $\lambda _1=4170$ Å and $\lambda _2 = 6010$ Å gives a blackbody color temperature of 9000 K (c.f., where the 9000 K solid line and the 7500 K dashed line approximately intersect). However, for certain conditions, $S_{\lambda , \rm{flare}} \gg S_{\lambda , q}$ where $S_{\lambda ,\rm{flare}}$ is the model flare surface flux spectrum, one may ignore the background subtraction errors (in which case, the problem must still be solved numerically or through a simple lookup table). Using the actual quiet Sun photospheric intensity, which is larger than a Planck function at these wavelengths, a 7500 K Planck function is sufficient to produce a 9000 K blackbody shape in the subtraction residuals. The surface fluxes of M dwarfs are much fainter than on G dwarfs. M-dwarf surface fluxes are lower than blackbody functions given by their effective temperatures, and the errors indicated in Fig. 37 are consequently smaller. For optically thick emission lines, similar issues arise if the model surface fluxes are comparable to the pre-flare surface fluxes (thus, it is fortuitous that large electron beam flux models produce spectra that are extremely bright compared to the pre-flare surface flux spectrum!).

• 2. Second, a color temperature can be fit to a spectrum of observed fluxes. A Planck function is typically fit to the excess (flare-only) spectrum, or to flare-only multi-band photometry data (Hawley et al. 2003), with $n\lambda$ points ($n > 2$) and Gaussian uncertainties, thus deriving the maximum likelihood parameters ${\hat{T}}_{{\rm{flare}}}, {\hat{X}}_{{\rm{flare}}}$. This can be accomplished with non-linear least squares through linearization and iteration for unconstrained parameter fits and with the Levenberg-Marquardt algorithm (e.g., IDL’s mpcurvefit) for constrained parameter fits. The color temperatures (e.g., $T_{{\rm{BB}}}$) that were quoted from blackbody fitting in Kowalski et al. (2013) were obtained by fitting Planck functions directly to the observed flare-only spectral fluxes (i.e., assuming $B_{\lambda }(T_{{\rm{flare}}}) \gg S_{\lambda , q}$). Also in that paper, the excess RHD model flux spectra were similarly compared to a Planck function temperature without subtracting a pre-flare from the Planck function. In this approach, the color temperature is just a convenient parametrization for comparing models to flare-only spectral observations, since temperatures are easier to compare than color indices in different filters or at different wavelengths (as noted in the main text). Whatever one does to the data, one should also do analogously to the model, whether the model is a blackbody or an RHD spectrum, irrespective of spectral type of the star (from RHD models of M dwarf flares that are optically thin in the chromospheric continuum radiation, it is actually desirable to isolate the chromospheric flare radiation by subtracting the pre-flare surface flux because the model photospheric spectrum, which does not have proper molecular opacities, otherwise shines through and obfuscates comparisons to observations).

Fig. 37

Lines of constant blackbody flare temperatures ($T_{{\rm{flare}}}$) and lines of constant blackbody color temperatures that are calculated after subtracting a quiescent blackbody surface flux spectrum, $S_{\lambda ,q} = B_{\lambda }(T_q)$, from the flare model spectrum, $S_{\lambda , \rm{flare}} = B_{\lambda }(T_{{\rm{flare}}})$. The color temperature calculations, here shown for $\lambda _1=4170$ Å and $\lambda _2=6010$ Å, from stellar flare data usually account for or consider these issues in M dwarf flares (Hawley et al. 1995; Kowalski et al. 2016), but hotter flare color temperatures in higher-mass stars could be a result of “oversubtraction” effects

The main purpose of subtracting a pre-flare spectrum from the flare observations is to remove the relatively large, wavelength-dependent flux from the area on the star that is not flaring so that the flare emission is more accurately characterized. A method for determining an accurate subtraction of the preflare spectrum based on molecular bandpass veiling by the flare continuum is described in Kowalski et al. (2013). Alternatively, a very wide slit may be employed, but spectral purity and wavelength resolution are sacrificed; see Kowalski et al. (2016).

In summary, the first method uses two filters and assumes that the flare spectrum is inherently blackbody radiation emitting at a single temperature $T_{{\rm{flare}}}$, and the second method is used for calculating color temperatures ($T_{{\rm{BB}}}, T_{{\rm{FcolorR}}}$) directly from a spectrum of flare-only fluxes assuming that the effects in Fig. 36 are not important (although this assumption can readily be relaxed in the second method).

1.2 Hydrogen continuum emissivity models

Kunkel (1970) calculated a grid of optically thin, LTE hydrogen recombination models with free-free contributions. These are often used as a starting place for comparisons to flare observations (e.g., Fig. 13). The equations are not included in Kunkel (1970) but were restated in Kowalski (2012). Here, we update the calculations with free-free Gaunt factors from CHIANTI and include the dissolved-level Balmer continuum components (longward of the Balmer edge at $\lambda = 3646$ Å) in the spectra. These are shown in Fig. 38 for a range of temperatures and a gas density of $\rho = 10^{-8}$ g cm$^{-3}$. The Balmer jump ratios are calculated from these spectra and are shown in Fig. 39. All spectra are available as FITS tables for a range of gas densities in the supplemental online material hosted on Zenodo.

Fig. 38

Hydrogen LTE emissivity calculations with bound-free (recombination) emissivity from H I and H$^-$ and free-free (thermal bremsstrahlung) emissivity from H II. The spectral shapes here are identical to the shapes of intensity spectra, $I_{\lambda }$, that emerge from optically thin, homogeneous slabs. A blackbody function and a mono-energetic synchrotron spectrum is shown for comparison. The smooth transitions from the Balmer jump to the optical wavelengths are caused by the “dissolved level Balmer continuum” (Dappen et al. 1987), which extends to redder wavelengths for larger electron densities

Fig. 39

Balmer jump ratios calculated from the LTE emissivity models in Fig. 38 and from the same calculations of the LTE emissivity in a hydrogen gas with a density of $\rho = 10^{-9}$ g cm$^{-3}$

At constant slab temperature, a large range of flare colors result from adjusting the continuum optical depths (Morchenko et al. 2015). In Fig. 40, we show the LTE continuum spectra from RHD model parameterizations (Kowalski and Allred 2018); the colors of these spectra are calculated and are shown in Fig. 13 as the “KA18 CC Models”. A range of reference column masses ($m_{{\rm{ref}}}$) were chosen to increase the continuum optical depths within the model chromospheric condensations. Occupational probabilities ($w_n$) for the dissolved level Balmer continuum opacities are included as in Fig. 38.

Fig. 40

(Left) LTE continuum flux spectra calculated from the parameterized evolution of chromospheric condensations from RHD models. (Right) Examples of two of the parameterized atmospheres

Appendix D Pressure broadening of hydrogen lines in flares

Electric fields split the degenerate energy levels of hydrogen atoms, Rydberg atoms, and hydrogen-like ions by an amount directly proportional to the field magnitude ($\propto E^1$; e.g., Condon and Shortley 1963; Bethe and Salpeter 1957; Gallagher 2006; Goldman and Cassar 2006), which is commonly known as the linear Stark, or Stark–Lo Surdo⁵³ effect. The splitting increases for larger n, and $n=1$ experiences only second-order energy shifts, which are small. Within a partially ionized gas, the electric fields from ambient, thermal charges (electrons and protons) within a Debye radius broaden spectral lines. This explains the broad hydrogen absorption profiles in main sequence A stars and white dwarfs (combined with, of course, large optical depths in their photospheric layers).

The symmetric, atomic energy-level shifts that are linearly proportional to an external, uniform electric field ${\mathcal {E}}$ are described by the quantum numbers nqm, according to the following absolute value (in cgs units):

$$\begin{aligned} | \Delta E^{(1)} | = | \frac{3 a_0 e}{2 Z} n q {\mathcal {E}} | \end{aligned}$$

(D1)

In Eq. (D1), n is the principal quantum number, q is the difference of two parabolic quantum numbers, $k_1 - k_2$, that satisfy $n = k_1 + k_2 + \mid m \mid + 1$, $m = m_l$ is the magnetic quantum number, $a_0$ is the Bohr radius, e is the electron charge, and Z is the nuclear charge (where $Z=1$ for hydrogen).

Figure 41 shows the splitting of the energy levels of hydrogen in a static electric field. The energy shifts were calculated using Eq. (D1) and the quadratic (second-order) and cubic (third-order) perturbation terms for hydrogen (Condon and Shortley 1963; Bethe and Salpeter 1957; Drake 2006; Goldman and Cassar 2006). The second-order and third-order terms create asymmetric splittings but are only important for hydrogen at very high electric field strengths. The maximum possible wavelength shift for an allowed Balmer transition ($n_{{\rm{lower}}} = 2$) and $\Delta m \pm 1,0$ follows from Eq. (D1) and is proportional to (van Dien 1949):

$$\begin{aligned} \bigg ( \frac{n_{j}^2}{n_{j}^2 - 4} \bigg )^2 ( n_{j} (n_{j} - 1) + 2 ) n_e^{2/3} \end{aligned}$$

(D2)

where $\big ( \frac{n_j^2}{n_j^2 - 4} \big )^2$ is the appropriate conversion from dE to $d\lambda$, and $n_j$ is the upper level of the Balmer transition. This formula, which goes as $n_j^2$ in the limit of large $n_j$, is sometimes used as an equivalent Doppler width that is added in quadrature with the true Doppler width in a Voigt profile. However, this approximation, among others (Sutton 1978), predicts inaccurate scalings for the line shapes and widths within a series. Thus, we recommend that these approximations not be used in place of proper line broadening calculations, which we briefly summarize next.

Fig. 41

Energy level splitting for hydrogen in a static external electric field as a function of the electric field (in SI units). A principal level n is split into roughly $2n-1$ (including the unshifted) sub-states. For even n, there is an unshifted component. Every split sub-state ($k_1, k_2, m$) is a mix of all the l-states for that value of n. The maximum splitting of the $n=5$ level shown here corresponds to 0.016 eV, which would correspond to a maximum wavelength displacement of $\Delta \lambda \approx 24$ Å for the H$\gamma$ transition. Only levels that can result in a dipole transition with $n=2$ have been plotted. The transitions with $\Delta m = 0$ are known as the $\pi$ components, and $\Delta m \pm 1$ are the $\sigma$ components. The field strength that forms a “saddle point” (potential energy maximum) and destroys all levels $>n$ in hydrogen is shown as a function of n (if there is a saddle-point in the potential, then tunneling is allowed as well; Bethe and Salpeter 1957). The field ionization of hydrogen has been observed (see the photograph from Rausch v. Traubenberg et al. (1930) that is reproduced in Montgomery et al. 2022). The Inglis-Teller relation corresponds to the critical field that causes the “bluest” sub-state of n to intersect with the “reddest” sub-state of $n+1$, which occurs at a much smaller electric field for a given n than the saddle point estimate (the ratio of the critical field to the saddle-point field is given by the factor known as $K_n$ in the literature). Holtsmark probability distributions of electric field magnitudes are shown in the top panel to illustrate an occupational probability ($w_n$) calculation for $n=11$ and $n_e = 10^{15}$ cm$^{-3}$. A “normal field”, $F_0$, is the electric field magnitude at an average interparticle separation; the critical fields for $n_e =10^{14}$ and $10^{15}$ cm$^{-3}$ are indicated by vertical dotted lines in the top panel

In a partially ionized gas, the electric microfield amplitude and direction vary from atom to atom and are not constant in time even at one atom because the charge distribution changes. The microfield amplitude probability is described by a probability distribution function (e.g., in the top panel of Fig. 41). A brief phenomenological summary of the microphysics of the broadening of hydrogen spectral lines follows. See Barklem (2016), Gigosos (2014), Gomez et al. (2022), Hubeny and Mihalas (2014) for more comprehensive and rigorous reviews.⁵⁴ The densities of perturbing, ambient (thermal) protons ($n_p$) and electrons ($n_e$) affect the total electric pressure broadening of radiative (dipole) transitions in hydrogen differently. The slow-moving protons and ions produce quasi-static (meaning that the microfield changes slowly compared to radiative/orbital timescales; see Griem 1974; Hummer and Mihalas 1988) electric microfields that split the energy levels of hydrogen into sub-states (Fig. 41), which are superpositions of the various orbital angular momentum states, l, within each n (“manifold”). The microfield magnitude is given by a Hooper probability distribution (which is a modified Holtsmark distribution that accounts for Debye screening and plasma correlations; Nayfonov et al. 1999). The collisional perturbations from ambient electrons are more complicated. Classically, the individual elastic electron impacts on each atom are randomly occurring and well-separated in time, thus causing a damping/Lorentzian profile function with broadening FWHM $\Gamma$ that is proportional to $n_e$. However, on very short timescales, which correspond to large frequency displacements from line center (in the far line wings) in Fourier-conjugate space, the surrounding electrons are nearly still and thus produce a quasi-static (Holtsmark) broadening proportional to $n_e^{2/3}$. The transition in electron broadening from its impact to its quasistatic regimes occurs at wavelength shifts, of order $|\lambda - \lambda _{{\rm{rest}}}| \approx 20$ Å (for the H$\gamma$ line). Thus, large optical depths can make the effects of the quasi-static limit of electrons important in the wings. Additionally, lifetime broadening occurs if there are inelastic effects, such as collisional transitions to states of different principle quantum number n (e.g., Hummer and Mihalas 1988; Bransden and Joachain 2000), but these are generally small for well-separated levels (and are ignored in the “no-quenching” assumption in the VCS73 theory).

The energy-level perturbations from protons and from electron collisions are not independent random variables, and the line profiles thus cannot be modeled as convolutions of two probability distribution functions, such as a Lorentzian convolved with a Holtsmark (or Hooper) distribution. This is because the electron collisions with hydrogen are transient (short duration) and non-adiabatic⁵⁵ (Smith et al. 1969; Sobel’man et al. 1995; Gigosos 2014; Hubeny and Mihalas 2014): they can cause transitions and redistribution among the microfield-split sub-states of the same n. Quantum calculations (Eq. (10) of Tremblay and Bergeron 2009, which is reproduced from Vidal et al. 1970) that account for the perturbing electrons in both of their limiting regimes are crucial. If the dynamical aspects (Stehlé 1994) of the ions are ignored, the normalized line profile functions are then the ionic microfield averages of the electronic broadening profiles calculated at each value of the quasi-static (QS), ionic microfield $\mathcal {E_{{\rm{QS}}}}$. When the spacing between the upper levels $n_j$ becomes small and the broadening causes the states to have a significant probability of overlapping, level dissolution occurs and the bound-bound transitions become narrower (Seaton 1990; Tremblay and Bergeron 2009) and less Holtsmark-like in the wings. The extensions of the VCS73 profiles that account for the effects of level dissolution through occupational probabilities (Sect. 10.2.1) are known as the TB09$+$HM88 profiles. Cho et al. (2022) discuss comparisons of high-density analytic profiles to new numerical calculations without level dissolution effects (see also Stehlé and Jacquemot 1993; Gomez 2017). To our knowledge, the TB09$+$HM88 profiles are yet the only quantum calculations extended to large $n_{j}$ that also include a framework, albeit phenomenological, for ionization-limit lowering consistently with the bound-bound opacities and line shapes. This method is adequate for model spectra of the A0 V star Vega (e.g., Kowalski et al. 2017b), and thus by extension for flare models, given that there is an unfortunate paucity of echelle observations in this spectral region in stellar (Fuhrmeister et al. 2008, 2011) and solar (Donati-Falchi et al. 1985) flares.