Rising Near-ultraviolet Spectra in Stellar Megaflares

Stellar flares are the most dramatic examples of variability that a cool star experiences while on the main sequence. Current understanding (A. F. Kowalski 2024) is that stellar flares are produced as a result of magnetic reconnection in the tenuous outer atmosphere of the star, which happens as a result of the jostling and colliding of magnetic fields that are rooted in the lower stellar atmosphere. Flares involve all layers of the stellar atmosphere, with a variety of physical processes, from plasma heating to particle acceleration to mass motions (A. O. Benz & M. Güdel 2010; Y. Notsu et al. 2024). These three factors result in flare emissions being produced across the electromagnetic spectrum. While flares on the Sun are the best studied due to the Sun's proximity, flaring itself is a product of magnetic activity and can be seen on a variety of stars with outer convective envelopes, from pre-main-sequence to young solar-type to evolved stars (H. Yang & J. Liu 2019; S. Okamoto et al. 2021). Apart from solar studies, M dwarfs are the most studied for stellar flares, due to a combination of population statistics (M dwarfs are the most common type of star in the nearby solar neighborhood; J. J. Bochanski et al. 2010) and flare statistics (M dwarfs tend to have high flaring rates; S. Candelaresi et al. 2014; H. Yang & J. Liu 2019).

Flare events can have prodigious releases of energy of up to 10³⁶ erg, or more than 10,000 times bigger than the largest solar flares (e.g., H. Maehara et al. 2012; R. A. Osten et al. 2016). In particular, the flare radiation in the ultraviolet (UV) can have significant impacts on nearby potentially habitable planets, which may be either helpful (e.g., P. B. Rimmer et al. 2018) or harmful (e.g., A. Segura et al. 2010; M. A. Tilley et al. 2019). The UV spectral region is key to understanding the impacts that flares have on exoplanet habitability by examining their effects on the atmosphere under a variety of assumptions about atmosphere composition and characteristics (A. Segura 2018).

In early studies of M dwarf flares, S. L. Hawley & G. H. Fisher (1992) established that the blue optical flare spectral energy distribution (SED) was consistent with a blackbody of temperature around 10⁴ K. This model was based largely on far-UV (FUV) and optical spectra in the impulsive phase of a flare first studied in S. L. Hawley & B. R. Pettersen (1991), and it has remained lore in the stellar flare community. The original explanation of the broadband increase, reprocessing of upper atmospheric (coronal) radiation, was demonstrated in subsequent papers (J. C. Allred et al. 2006) to be negligible compared to the energy deposited in the lower atmosphere by electron beams. More recent high-cadence, low-resolution flare spectral atlases have shown the predominance in the blue optical for a continuum component consistent with such a blackbody, as well as hydrogen Balmer jump emission and higher-order hydrogen Balmer lines (A. F. Kowalski et al. 2010, 2013). Despite the evidence for multiple components, the functional form of a single-temperature blackbody has persisted in explaining flare continuum emission from the UV to optical wavelengths (e.g., R. O. P. Loyd et al. 2018a; M. N. Günther et al. 2020; W. S. Howard et al. 2020; S. Okamoto et al. 2021; V. L. Berger et al. 2024; A. D. Feinstein et al. 2024, among many others). Spectroscopic (C. S. Froning et al. 2019; W. S. Howard et al. 2023) and broadband (C. E. Brasseur et al. 2023; R. R. Paudel et al. 2024) investigations have started to consider multicomponent models from the FUV to the near-infrared (NIR), but spectra that span the Balmer jump can be helpful in limiting the parameter space and breaking degeneracies (A. F. Kowalski et al. 2019).

The optical and UV portions of the flare radiation trace emission produced in the lower stellar atmosphere. It was perplexing given the generally good agreement between solar and stellar flares in terms of scaling relations for plasma heating (e.g., K. Shibata & T. Yokoyama 1999, 2002) that the large optical continuum enhancements seen in M dwarf flares could not be reproduced using solar flare models (J. C. Allred et al. 2006). Current state-of-the-art radiative-hydrodynamic (RHD) models of solar and stellar flares (J. C. Allred et al. 2005; A. F. Kowalski et al. 2017b) have demonstrated an ability to reproduce continuum flare features, as well as the strongest line emission, in the near-UV (NUV) through blue optical (see also C. E. Brasseur et al. 2023) through the action of accelerated particles impacting the lower stellar atmosphere. Solar studies generally utilize hard X-ray observations of nonthermal bremsstrahlung emission and radio observations to diagnose the presence and action of accelerated particles (G. D. Holman et al. 2011; E. P. Kontar et al. 2011; S. M. White et al. 2011). The former is not available for stellar studies due to sensitivity issues, and radio observations in the microwave regime have not had the frequency coverage to constrain the properties of the accelerated particles in stellar flares (R. A. Osten et al. 2005). Earlier studies (M. Güdel et al. 1996) have suggested a disconnect in particle acceleration relative to plasma heating for the stellar flares studied versus solar flares. Also, the large gap between solar and stellar optical continuum flares in energy-duration diagrams can be attributed to physical scaling relations that imply larger magnetic fields (and flaring volumes) in active stars (H. Maehara et al. 2015; K. Namekata et al. 2017; see also studies of X-ray superflares, such as F. Favata et al. 2000 and R. A. Osten et al. 2016). Confirmation that the UV spectral region can diagnose the characteristics of magnetic energy and accelerated particles in stellar flares would open new regimes for exploring solar–stellar connections.

Stellar flares have been studied at UV wavelengths predominantly in the FUV (e.g., S. L. Hawley et al. 2003; K. France et al. 2016; R. O. P. Loyd et al. 2018a) for many decades, and current generations of planet-transit-hunting telescopes like the Transiting Exoplanet Survey Telescope (TESS) and Kepler routinely pick up the optical counterparts of these events (J. R. A. Davenport 2016; M. N. Günther et al. 2020). Despite this seemingly advanced topic, recent discoveries have thrown doubt on how well we understand the stellar flare process. Some of this arises from the relative lack of observations in the NUV wavelength region, despite this being where such a 10⁴ K blackbody would peak. Surprisingly, in a recent multiwavelength flare study, only 25% of the NUV flares on M dwarfs had an optical counterpart (R. R. Paudel et al. 2024), and some NUV flares clearly show a much greater response than a T ≈ 9000−10⁴ K blackbody would predict from the available optical constraints (J. A. G. Jackman 2022; C. E. Brasseur et al. 2023; J. A. G. Jackman et al. 2023; see also A. F. Kowalski et al. 2019). These results motivate renewed scrutiny of the NUV spectral region in stellar flares and form the basis of a Hubble Space Telescope (HST) Treasury program designed to use the NUV as a "fulcrum" to bridge the more-often-studied FUV and the very well-studied optical regions (HST Guest Observer 17464 "From High-Energy Particle Beam Heating in Stars to Ozone Destruction in Planets: NUV Spectra as the Fulcrum for a Comprehensive Understanding of Flaring M Dwarf Systems"). In the following sections, we describe early results from this Treasury program, highlighting two remarkable events and what they can tell us about the physical processes at work in stellar flares. Section 2 describes the data, Section 3 describes the light curve and spectral analysis, Section 4 discusses the application of RHD models to understand the continuum and line emission, Section 5 discusses the results, and Section 6 concludes. Appendix A describes the data reduction, and Appendices B and C present light curves and flare spectra that supplement the figures in the main text.

CR Dra is a low-mass, M dwarf star that is well known to flare in the NUV (B. Y. Welsh et al. 2006; C. Million et al. 2016; J. A. G. Jackman et al. 2024). It is a spectroscopic binary with a trigonometric parallax distance of 20.15 ± 0.17 pc (Gaia Collaboration et al. 2018, 2023), consisting of an M0 star and a lower-mass companion (V. S. Tamazian et al. 2008). The combined spectral type is M1e (I. N. Reid et al. 1995; S. L. Hawley et al. 1996), and there is a single source in Gaia with a G-band absolute magnitude that is well above the main sequence at the color of other M1 stars (A. F. Kowalski 2024). ¹⁰

The Cosmic Origins Spectrograph (COS) on HST observed two extraordinary ("mega" ¹¹ ) flares on CR Dra during the observations for Guest Observer (GO) Treasury Program 17464. The NUV G230L grating was employed with the 2950 Å central wavelength, giving spectral coverage at λ ≈ 1683–2082 Å in Stripe A (NUVA) and λ ≈ 2773−3170 Å in Stripe B (NUVB). The linear dispersion is 0.39 Å pixel⁻¹, and the resolving power, R, varies from 2100 to 3900. The photons are time tagged, which allows for spectra to be extracted over any time interval. TESS (G. R. Ricker et al. 2015) simultaneously observed CR Dra at 20 s and 120 s cadences in Sectors 76 and 78 in Cycle 6. Appendix A describes the data processing in detail.

The broadband HST light curves of the two megaflare events are shown in Figure 1. These light curves exclude the Mg ii h and k emission lines and the pixels within each stripe that are affected by detector shadow (Appendix A.1). The light curves for NUVA and NUVB are given as average flare-only flux densities over the respective wavelength ranges, $\left\langle {f}_{\lambda }^{{\prime} }\right\rangle$, where the prime symbol indicates that a preflare flux is subtracted. Flare Event 1 consists of many extraordinary variations superimposed on two major peaks. The HST observations stop during the flare due to the end of the orbit. This time corresponds to the end of the fast decay phase in the Δt = 20 s TESS light curve (Appendix B). The energy in the TESS bandpass is E_TESS = 7.5 × 10³³ erg, and the peak luminosity is 6.1 × 10³⁰ erg s⁻¹, making this an extremely energetic event. Flare Event 2 is over 10× more energetic, having E_TESS ≈ 10³⁵ erg. Flare Event 2 is the most luminous NUV flare that has been reported (to our knowledge) on an M dwarf star; the spectral luminosity in NUVA is over a factor of 10 larger than the most luminous flare spectra compiled from the International Ultraviolet Explorer (IUE)/Short-Wavelength Prime Spectrograph in K. J. H. Phillips et al. (1992). The rise phase increases the NUVB stellar flux (≈1.24 × 10⁻¹⁴ erg cm⁻² s⁻¹ Å⁻¹) by a factor of ≳200, but the buffer filled as HST stopped observations due to an overlight condition. After this time, the broadband optical/NIR response continues to increase, and the TESS luminosity reaches 8.8 × 10³¹ erg s⁻¹ at peak (Appendix B). Flare Event 2 is about 25 times more impulsive than Flare Event 1 in the TESS band, where we follow A. F. Kowalski et al. (2013) and define the impulsiveness index of a light curve as the flare-only peak flux divided by the full width at half-maximum (FWHM) duration, t_1/2.

Zoom In Zoom Out Reset image size

Within each flare event, the NUVA light curves are more impulsive and show peak spectral flux densities that are larger than at the longer wavelengths of NUVB (Figure 1). The impulsiveness index is 2–2.3× larger and the peak fluxes are ≈1.4–1.5× higher in NUVA during the two major peaks in Flare Event 1. By the end of the HST observations of Flare Event 2, the peak flux in NUVA is a factor of ≈1.7 larger than NUVB. These properties are qualitatively similar to broadband Galaxy Evolution Explorer (GALEX)/FUV (1350–1750 Å) and GALEX/NUV (1750–2800 Å) filter ratios during stellar flares (R. D. Robinson et al. 2005; B. Y. Welsh et al. 2006; C. Million et al. 2016; S. W. Fleming et al. 2022; V. L. Berger et al. 2024). The impulsiveness differences between NUVA and NUVB are also qualitatively similar to S. L. Hawley et al. (2003), who found a power-law relationship between the time evolution of FUV continuum regions (averaged over λ = 1420−1452 Å and 1675−1710 Å) and the time evolution of the U-band (λ ≈ 3200−4000 Å) flux. Over the entire impulsive phase light curve of Flare Event 1, we calculate the power-law index, α, from ${\mathrm{log}}_{10}\left\langle {f}_{\lambda ,\mathrm{NUVA}}^{{\prime} }\right\rangle =\beta +\alpha {\mathrm{log}}_{10}\left\langle {f}_{\lambda ,\mathrm{NUVB}}^{{\prime} }\right\rangle$ as α ≈ 1.5. Over the rise phase of Flare Event 2, we calculate α ≈ 1.3. These values are similarly >1, like the value of α = 1.65 that is reported in the S. L. Hawley et al. (2003) data. A nonlinear relationship means that the flux ratios vary as functions of the fluxes, which has remained unexplained in U-band and FUV comparisons. However, previous interpretations of broadband UV flare data (e.g., S. L. Hawley et al. 2003; U. Mitra-Kraev et al. 2005; R. D. Robinson et al. 2005; B. Y. Welsh et al. 2006; V. L. Berger et al. 2024) have been severely limited to simple models, such as emission line slabs or single-temperature blackbody functions, due to the lack of an NUV spectral fulcrum connecting the blue edge of the U band to the FUV.

For the first time, we present UV spectral observations at λ ≲ 3200 Å that show that there are two blackbody color temperatures in the NUV with T_NUVA > T_NUVB. The spectra and the fitting results from the second major peak in Flare Event 1 and from the rise phase of Flare Event 2 are shown in Figure 2. Continuum radiation dominates at both λ ≲ 3200 Å (in our NUVB stripe) and ≲2100 Å (in our NUVA stripe, which some may consider as part of the FUV or mid-UV) in the impulsive phases of these flares. The inferred NUVB continua are strikingly flat within this range; we use the nonlinear least-squares Levenberg–Marquardt algorithm and fit blackbody temperatures of $\hat{T}\,\approx \,\rm{9000\unicode{x02013}10,000}$ K. The NUVA spectra are overall brighter than the NUVB spectra, as noted in the continuum light curves (Figure 1), and they exhibit an ascent toward shorter wavelengths. We independently fit blackbody temperatures to NUVA and find that ${\hat{T}}_{\mathrm{NUVA}}$ ranges from ≈16,000 K to nearly 18,000 K; see Appendix C for several other flare spectra and the corresponding fits. A two-temperature blackbody fit to NUVA and NUVB results in a high-temperature component that is only a few hundred degrees hotter than the fit to only NUVA. A single-temperature fit to the ratio, $\left\langle {f}_{\lambda ,\mathrm{NUVA}}^{{\prime} }\right\rangle /\left\langle {f}_{\lambda ,\mathrm{NUVB}}^{{\prime} }\right\rangle$, of the wavelength-averaged flux densities in the spectrum of Flare Event 2 gives a blackbody color temperature of ${\hat{T}}_{\mathrm{NUV}}=14,790\pm 430$K. This is not representative of the continuum shape within either NUV stripe, thus demonstrating that one filter ratio spanning a large wavelength range across the UV does not accurately characterize the flare SED.

Zoom In Zoom Out Reset image size

There are probably many faint emission lines (J. W. Cook & G. E. Brueckner 1979; J. G. Doyle & J. W. Cook 1992) that blend together and form a pseudocontinuum on top of the bona fide continuum radiation. However, the contribution appears to be energetically small in the spectra of these flares. We calculate that the continuum fits comprise 97%–99% of the wavelength-integrated fluxes within NUVA and NUVB. This is expected based on knowledge from an echelle flare spectrum of the impulsive phase of an M dwarf flare at λ = 3300–3800 Å, which shows a very small contribution from emission lines in comparison to a T ≈ 11,000 K blackbody continuum model fit (B. Fuhrmeister et al. 2008). We have identified regions in the HST spectra (indicated by the gray bands in Figure 2) that correspond to the low points between obvious emission lines (e.g., J. G. Doyle & J. W. Cook 1992). These regions best represent the bona fide continuum fluxes, which are used in the fitting procedure. Notably, the best-fit blackbody temperatures are largely robust to the choice of fitting windows outside of the emission lines of Al iii (1854.716, 1862.790 Å) in NUVA (indicated by vertical dashed lines) and Mg ii h and k in NUVB. The other two vertical dashed lines in NUVA indicate Si ii λ1808.0 and λ1816.9, which are among the brightest in the preflare and gradual phase spectra of solar flares (J. G. Doyle & J. W. Cook 1992; P. J. A. Simões et al. 2019).

The gradual temporal evolution of the Mg ii emission line fluxes further justifies attributing the dominant source of radiative energy to continuum radiation. We integrate over the wavelengths of the resonance lines and the nearby subordinate lines, which also brighten in the flares, and we show the continuum-subtracted flux evolution at the 20 s cadence of TESS in Figure 1. In Flare Event 1, there is a striking lack of similarity between the Mg ii emission line response and the two main wavelength-integrated peak fluxes, which are dominated by the continuum radiation. The Mg ii light curve rises and peaks more gradually, and it remains significantly elevated above preflare levels while still decaying at the start of the next HST orbit (Appendix B).

Previous RHD modeling (A. F. Kowalski 2022) of the multiwavelength observations of the impulsive phase of the Great Flare of AD Leo (S. L. Hawley & B. R. Pettersen 1991; S. L. Hawley & G. H. Fisher 1992) provides a starting point for explaining the new HST spectra in terms of collisional heating from nonthermal, power-law electron beams. A. F. Kowalski (2022) find that several combinations of models with large electron beam energy flux densities are possible explanations for the hydrogen Balmer spectra and broad-wavelength constraints, including the rather flat shape of the FUV continuum, in the Great Flare. We use the publicly available output (on Zenodo; A. F. Kowalski et al. 2024a) from a grid (A. F. Kowalski et al. 2024b) of RADYN (M. Carlsson & R. F. Stein 1992, 1995, 1997; J. C. Allred et al. 2015; M. Carlsson et al. 2023) stellar flare models and perform a complete parameter space search for radiative surface flux continuum spectra that are able to explain the NUVA and NUVB shapes in Flare Event 2 (Figure 2, bottom). The shapes of the model spectra are determined by the ratios of the detailed continuum fluxes at λ = 1800 and 2069 Å in NUVA and at λ = 2830 and 3080 Å in NUVB. These are compared to the ratios of flare-only fluxes in the data at averages over λ = 1778.7–1805 Å and 2065.16–2080.14 Å in NUVA and over λ = 2820–2850 Å and 3011–3170 Å in NUVB.

Only the most energetic electron beam heating model (cF13-500-3) with an energy flux density of 10¹³ erg cm⁻² s⁻¹, a low-energy cutoff of E_c = 500 keV, and a number-flux power-law index of δ = 3 above this low-energy cutoff can account for the rising slope of Flare Event 2 within NUVA. This continuum model is shown in Figure 3. However, the model fails to fully account for the flux at the longer wavelengths in NUVB. We thus follow A. F. Kowalski (2022) and fit linear superpositions of every combination of two RHD model spectra to the four continuum regions in the data,

$$\begin{eqnarray}&&{f}_{\lambda ,\mathrm{obs}}^{{\prime} }=\displaystyle \frac{{R}_{\mathrm{star}}^{2}}{{d}^{2}}\left({\hat{X}}_{1}{F}_{\lambda ,1}^{{\prime} }+{\hat{X}}_{2}{F}_{\lambda ,2}^{{\prime} }\right),\end{eqnarray} \tag{ 1 }$$

where ${\hat{X}}_{1}$ is the best-fit filling factor of RHD model component 1, ${F}_{\lambda ,1}^{{\prime} }$ is the radiative surface flux model of component 1 with the preflare subtracted, R_star = 4 × 10¹⁰ cm is an approximate radius of the flare star, and d is the distance to CR Dra (see A. F. Kowalski 2022 for other details). The fitting is performed to 98% of the observed flux densities, which allows the models to better account for the continuum distribution that underlies a pseudocontinuum pedestal of faint emission lines within the noise. The filling factor of each component is directly proportional to the area of the respective flaring source at the star. We report ${\hat{X}}_{\mathrm{rel}}={\hat{X}}_{2}/{\hat{X}}_{1}$, and we leave out the factor ${\hat{X}}_{1}\tfrac{{R}_{\mathrm{star}}^{2}}{{d}^{2}}$ from the labels of the model combinations in the figure legends.

Zoom In Zoom Out Reset image size

As in the modeling of the AD Leo Great Flare, a combination of the cF13-500-3 with an electron beam flux density of 2 × 10¹² erg cm⁻² s⁻¹, a smaller low-energy cutoff (E_c = 37 keV), and a hard power-law index, δ = 2.5, gives an adequate explanation for the overall properties of the full UV continuum range in the new data. For Flare Event 2 (Figure 3), the best-fit filling factors, $\hat{X}$, of each model imply a ratio of the area of the low flux density model (m2F12-37-2.5) to the area of the high flux density model (cF13-500-3) of ${\hat{X}}_{\mathrm{rel}}={\hat{X}}_{2}/{\hat{X}}_{1}={\hat{X}}_{\mathrm{low}}/{\hat{X}}_{\mathrm{high}}\approx 3.4$. The same combination of models results in a value of ${\hat{X}}_{\mathrm{rel}}=7.3$ for the spectrum of Flare Event 1 (Appendix C), which suggests that the combination is a rather flexible model for the impulsive phase of very energetic flares. Section 5 discusses interpretations of the two-component RHD model fitting.

Though the two-component RHD model of the NUVB is still not flat enough, detailed radiative transfer modeling of the Mg ii emission lines lends credence to the general approach of superposing these RHD models. We use the RH code (H. Uitenbroek 2001) to calculate the spectra with the 10-level quintessential model Mg ii atom (J. Leenaarts et al. 2013) and the updated quadratic Stark damping (Y. Zhu et al. 2019) from the STARK-B database (M. S. Dimitrijević & S. Sahal-Bréchot 1995; S. Sahal-Bréchot et al. 2011). ¹² Figure 4 shows representative spectra of Mg ii that are calculated from snapshots of the cF13-500-3 and m2F12-37-2.5 models. The model spectra are convolved with the G230L line-spread function and binned to the linear dispersion (0.39 Å pixel⁻¹). A modest scaling of the total model shows general agreement with the observations of the early rise phase of Flare Event 2. The Mg ii h and k emission lines at high spectral resolution are notoriously difficult to reproduce in solar flare models (F. Rubio da Costa & L. Kleint 2017; Y. Zhu et al. 2019; A. Sainz Dalda & B. De Pontieu 2023; G. S. Kerr et al. 2024a), but the models here appear to have a rather reasonable amount of radiative energy loss at optically thick lines relative to the radiative energy loss through the continuum (see K. Namekata et al. 2020; A. F. Kowalski 2022; A. F. Kowalski et al. 2024b for comparisons of the model grid to hydrogen Balmer lines).

Zoom In Zoom Out Reset image size

Until now, there have been only a handful of NUV spectra during the impulsive phase of stellar flares (R. D. Robinson et al. 1995; S. L. Hawley et al. 2007; A. F. Kowalski et al. 2019; J. A. G. Jackman et al. 2024). The AD Leo Great Flare spectra in the NUV from the IUE/Long-Wavelength Prime spectrograph are saturated in the impulsive phase; thus, spectra have been merged from different phases of the flare with heterogeneous exposure times. This merged flare spectrum, shown in Figure 3, has been widely used as a template flare SED in exoplanet photochemistry modeling (A. Segura et al. 2010; O. Venot et al. 2016; M. A. Tilley et al. 2019; R. J. Ridgway et al. 2023) and assessments of UV-dependent pathways to the origin of life (S. Ranjan et al. 2017; M. Z. Armas-Vazquez et al. 2023). In the fiducial Great Flare template, the NUV is more luminous than the FUV. The emission line fluxes are larger relative to the continuum radiation than in the CR Dra flare spectra, as expected for a spectrum that was obtained from the gradual decay phase (S. L. Hawley & B. R. Pettersen 1991). Based on spectra and photometry in the optical and U band of other flares (S. L. Hawley & G. H. Fisher 1992; S. L. Hawley et al. 1995, 2003; B. Fuhrmeister et al. 2008; A. F. Kowalski et al. 2013), it was expected that with sufficient signal-to-noise and a large enough flare, the spectral flux density in the impulsive phase would be seen to clearly peak and turn over at λ ≈ 3000 Å, like a T ≈ 10⁴ K blackbody or blackbody-like spectrum within the hydrogen Balmer continuum.

Now, we see that the inferred continuum radiation in the short-wavelength NUV rises well above the continuum fluxes in the longer-wavelength UV around λ = 3000 Å in two megaflare events on CR Dra. We conclude that there is a large amount of energy from (at least some) stellar flares in the impulsive phase that is radiated in the short-wavelength NUV continuum—and, by a sensible extrapolation, in the FUV as well. This energy has previously been unaccounted for through the AD Leo Great Flare template (Figure 3), blackbody extrapolations from optical photometry (e.g., T. Shibayama et al. 2013; M. N. Günther et al. 2020), and RHD extrapolations based on Balmer jump strengths (A. F. Kowalski et al. 2019; A. F. Kowalski 2022). In Figure 5, we show a T = 10⁴ K blackbody anchored to the TESS flare-only flux (blue star symbol) averaged over UTC 2024 May 17 11:45:00.85–11:47:00.85 during Flare Event 2. The HST/COS spectra in this figure correspond to the same time interval, and the observed fluxes within NUVA and NUVB are 6.1× and 2.6× larger, respectively, than the extrapolation of the T = 10⁴ K blackbody.

Zoom In Zoom Out Reset image size

Hotter blackbody color temperatures, ranging from T ≈ 15,000 K to T ≈ 40,000 ± 10,000 K, have been recently reported in several stellar flare events with spectra at shorter wavelengths in the FUV, either at λ ≲ 1700 Å or at λ ≲ 1400 Å (R. O. P. Loyd et al. 2018b; C. S. Froning et al. 2019; M. A. MacGregor et al. 2021; A. D. Feinstein et al. 2022). We suggest that the increasing values of f_λ within the NUVA range represent the long-wavelength extension of this phenomenon, which is the most luminous component of the flare radiation in our data. Even with two-temperature blackbody fits (Section 3), we do not find temperatures larger than 20,000 K within the rise phase spectra of Flare Event 2. Thus, there is no evidence in our spectra that is suggestive of extremely large temperatures near T = 30,000–50,000 K, which have been claimed in the optical (W. S. Howard et al. 2020) and the FUV (R. D. Robinson et al. 2005; C. S. Froning et al. 2019) in other flares. As is well known, the FUV has many bound–free edges and bright resonance lines (C ii, Si iv, C iv) from metals that make interpretation difficult, especially in broadband data; the quiescent FUV spectrum of an M1 star has comparably large blackbody temperatures as in flares (A. D. Feinstein et al. 2022). The NUVA, however, is predominantly affected by H i bound–free opacities at the expected temperatures of flare chromospheres (e.g., A. García-Gil et al. 2005; A. F. Kowalski 2025, in preparation). Thus, the NUVA wavelength range provides straightforward comparisons to currently available model predictions. The simultaneous spectral constraints within the NUVB range connect the continuum radiation that rises into the FUV to the flare flux at the blue edge of the U band, which has traditionally been used in stellar flare studies over many decades (e.g., C. H. Lacy et al. 1976). In the smaller flares of S. L. Hawley et al. (2003; see flare 8 in their Figure 10), the FUV continuum fluxes are clearly fainter than the U-band fluxes. However, we note that the FUV is well above their T ≈ 9500 K blackbody fit; less energetic and less luminous flares (e.g., similar to those in A. F. Kowalski et al. 2019) will be analyzed from the HST Treasury program in future work.

Two phenomenological blackbody color temperatures are necessary to separately account for the blue optical (λ = 4000–4800 Å) and red optical (λ = 5900–7000 Å) continuum spectra in stellar flares (A. F. Kowalski et al. 2013, 2016); see, for example, the IF3 event peak spectrum in Figure 32 of A. F. Kowalski et al. (2013) and the accompanying discussion. The peak flare-only optical spectrum of this flare is shown in Figure 5. In this event, a cooler blackbody can be fit to the redder wavelengths closer to the head of the Paschen continuum, while blue optical wavelengths at λ = 4000−4800 Å rise above this; a hotter blackbody model temperature is needed (A. F. Kowalski et al. 2013). We have shown that a similar empirical phenomenon occurs across the NUV beginning at wavelengths just shortward of the U band. This potential connection between blackbody modeling in different wavelength regimes suggests that the NUV and optical continua in the impulsive phases of stellar flares are similarly affected by a common opacity source. H i bound–free opacity, as in our RHD models with very large electron beam heating rates (A. F. Kowalski et al. 2017b), is dominant in both regimes, while other opacity sources such as H⁻ are not. A best-fit combination of the RHD models (Section 3) is shown in Figure 5 and predicts generally similar properties as the IF3 optical peak spectrum: a small Balmer jump and decreasing color temperatures ¹³ (not shown) toward the wavelengths of the TESS bandpass (λ ≳ 6000 Å). Although simultaneous optical spectra are not available for Flare Event 2 to constrain the Balmer jump and optical emission lines, the observed TESS flare-only flux in Figure 5 suggests that the extrapolation of the RHD model to the NIR is relatively robust. The cF13-500-3 component in the two-model fit extends almost precisely to the TESS data constraint, while the two-component fit to the NUV predicts a TESS filter-weighted flux that is about 2× too large. RHD models with smaller surface fluxes than the m2F12-37-2.5 and cooler blackbody curves with T ≈ 4200 K, which better account for the NUVB shape in two-component blackbody fitting (Section 3), predict optical fluxes far above the TESS constraint. Also, much smaller model surface fluxes require filling factors that are larger than the visible stellar atmosphere, which is probably unrealistic. Detailed investigation into the NUV and TESS relationship in all of the flares in the Treasury program will be the subject of a second paper (R. A. Osten et al. 2024, in preparation).

Other sources of large heating rates that are similar to the E_c = 500 keV model may result in alternative explanations of the CR Dra megaflare spectra. We leverage the new and unique capabilities of the FP (J. C. Allred et al. 2020) and FP+RADYN codes (J. C. Allred et al. 2022; G. S. Kerr et al. 2023, 2024a) to calculate the transport of proton beam energy and momentum in an M dwarf atmosphere. As generally expected from the equations for energy loss in a cold target (e.g., A. G. Emslie 1978; S. L. Hawley & G. H. Fisher 1994; see Figure 25 of A. F. Kowalski 2024), we calculate that a proton beam with E_c = 10 MeV, a power-law index of δ = 3, and an energy flux density of 10¹³ erg cm⁻² s⁻¹ above the cutoff produces a thermal response in the low stellar atmosphere that is similar to the fully relativistic, high-energy electron beam heating model (cF13-500-3). Thus, a superposition of a high-energy proton beam simulation with a lower-energy electron beam simulation (e.g., the m2F12-37-2.5 model) could plausibly serve as an explanation for the megaflare spectra. It is interesting that the lower atmospheric impact sites of accelerated protons in solar flares are significantly displaced from the spatial locations of hard X-ray footpoint sources, which are associated with nonthermal electron energy deposition (G. J. Hurford et al. 2006). High-energy (E ≫ 1 MeV) protons/ions produce a variety of unique nonthermal radiation signatures (R. J. Murphy et al. 1997; G. H. Share et al. 2004; N. Vilmer et al. 2011), but these are not currently detectable at stellar distances (Y. Song et al. 2024). G. D. Fleishman et al. (2022) discuss potential signatures in gyrosynchrotron radiation, which is often prominent in stellar flares, that might help to clarify the roles of high-energy ion beams.

Among the current generation of RHD models of stellar flares with electron beam heating (J. C. Allred et al. 2006; A. F. Kowalski et al. 2017b, 2024b; K. Namekata et al. 2020), only the largest heating rates in the low chromosphere and temperature minimum region reproduce the continuum radiation in the NUVA spectra of both flare events on CR Dra. Much smaller electron beam flux densities (≈10¹⁰–10¹¹ erg cm⁻² s⁻¹) are widely used to model solar flares (e.g., J. C. Allred et al. 2005; L. Fletcher & H. S. Hudson 2008; D. Kuridze et al. 2015; F. Rubio da Costa et al. 2016; M. Carlsson et al. 2023; G. S. Kerr et al. 2024a, 2024b; P. J. A. Simões et al. 2024b), and the low-energy cutoffs are typically constrained to be E_c < 40 keV (e.g., A. M. Veronig et al. 2005; L. Kleint et al. 2016; A. Warmuth & G. Mann 2016). In Figure 3, we demonstrate that an electron beam with an energy flux density of 5 × 10¹¹ erg cm⁻² s⁻¹ and a low-energy cutoff of 25 keV produces a Balmer continuum spectrum that has the opposite slope as the observations. As expected, a simple optically thin, static slab model calculation in local thermodynamic equilibrium at T = 10,000 K (dashed red line; W. E. Kunkel 1970; A. F. Kowalski 2024) is even worse. This clearly demonstrates that the continuum radiation in stellar flares is not as simple as a static, uniform, optically thin slab at T = 10,000 K as assumed by some (P. J. A. Simões et al. 2024a) in the solar community. Other possible sources of very bright NUVA fluxes, such as a hot T ≳ 10⁷ K optically thin thermal source, are not plausible hypotheses (S. L. Hawley & G. H. Fisher 1992). C. S. Froning et al. (2019) use a semiempirical hot, dense chromospheric condensation to model large blackbody color temperatures in the FUV. However, the densities and temperatures in this model component are not self-consistently produced in current RHD models.

The larger inferred low-energy cutoffs and flux densities of electron beams in M dwarf flares may effectively originate during magnetic reconnection (see discussion in A. F. Kowalski et al. 2024b) and particle acceleration (e.g., R. J. Hamilton & V. Petrosian 1992). However, there has been little work on the origin of accelerated particles in stellar atmospheric conditions not typically found in the solar corona. Alternatively, A. F. Kowalski (2023) suggests that similar beam distributions to those with large E_c values could result from time-dependent transport effects (E. P. Kontar et al. 2012; see also M. Karlický & E. P. Kontar 2012; R. Pechhacker & D. Tsiklauri 2014; L. F. Ziebell & P. H. Yoon 2022; V. Annenkov & E. Volchok 2023) in strong magnetic fields, which may saturate the low coronae of dMe stars (S. M. White et al. 1994) and inhibit ion-acoustic beam instabilities (K. W. Lee & J. Büchner 2011; T. C. Li et al. 2012). The beam density of the cF13-500-3 model during coronal propagation is ≈4 × 10⁸ cm⁻³, which is well below the ambient densities determined from X-ray spectra of dMe stars outside of flaring times (R. A. Osten et al. 2006; C. Liefke et al. 2008). Thus, the Buneman instability and the availability of electrons are not concerns. Gyrosynchrotron radiation at optically thin radio frequencies (R. A. Osten et al. 2016; I. I. Tristan et al. 2024, in preparation), ALMA flare emission (which has been observed to be temporally correlated with a highly impulsive response in the FUV; M. A. MacGregor et al. 2021), and linear polarization in U-band bursts (G. Beskin et al. 2017) could help to place constraints on the number of relativistic electrons and the strengths of magnetic fields in stellar flares (A. M. MacGregor et al. 2020).

We interpret the two-component electron beam model fits (Section 4) to represent heterogeneous flaring sources on stars in analogy to chromospheric sources in solar flares. Notably, the ratio of areas, ${\hat{X}}_{\mathrm{low}}/{\hat{X}}_{\mathrm{high}}={\hat{A}}_{\mathrm{low}}/{\hat{A}}_{\mathrm{high}}\approx 3$, that we infer from the average rise phase spectra of Flare Event 2 is a very similar value to the ratio of the area of bright ribbons to the area of bright kernels in some solar flares (e.g., A. F. Kowalski 2022). Thus, the areas heated by high-energy electrons with high flux density (cF13-500-3) could represent concentrated heating in kernels, which are surrounded by extended ribbons that are heated by lower-energy electron beams (m2F12-37-2.5). ¹⁴ The impulsiveness calculations from the light curves in NUVA and NUVB (Section 3; Figure 1) could further support a model framework wherein two rather distinct heating components with different temporal evolution contribute to the total flare continuum flux; within the solar analogy, the kernels are transient while the ribbons are persistent. The inferred flare areas, ${\hat{A}}_{\mathrm{high}}\approx 2\times {10}^{19}$cm² and ${\hat{A}}_{\mathrm{low}}\approx 6.5\times {10}^{19}$ cm², are larger in the CR Dra stellar megaflares than the total areas directly measured in solar flare kernels and ribbons, which typically range from 10¹⁶ to 10¹⁸ cm² but sometimes extend up to ≈2 × 10¹⁹ cm² at certain wavelengths (D. F. Neidig et al. 1994; L. Fletcher et al. 2007). As additional comparison, the inferred flare areas in the Great Flare of AD Leo with the mF13-500-3 and m2F12-37-2.5 models (A. F. Kowalski 2022) are factors of 10 and 2 smaller than the respective RHD components (cF13-500-3 and m2F12-37-2.5) that are fit to Flare Event 2 on CR Dra. The CR Dra megaflare areas are yet only a small fraction of the visible stellar hemisphere, $\hat{X}\approx 0.01$, which implies compact sources.

The spatial resolution of the solar data constrains how spectral and temporal heterogeneity across chromospheric flare sources contributes to Sun-as-a-star signals (K. Namekata et al. 2022; A. G. M. Pietrow et al. 2024). K. Namekata et al. (2022) demonstrate that the brightest and broadest Hα spectra differ markedly from spatially integrated spectra, which more closely resemble spectra from the diffuse regions with weaker emission line intensity. Ostensibly, this supports superposing RHD model spectra to represent heterogeneous stellar flare sources. Similar comparisons of solar data from the Interface Region Imaging Spectrograph (IRIS) to other flares in the HST Treasury program that are similar in energy to large solar flares would be valuable for constraining the number of heterogeneities in UV sources. For example, it would be valuable to determine if any very bright kernels in solar flares (e.g., D. B. Jess et al. 2008; S. Krucker et al. 2011) exhibit a brighter and more impulsive FUV continuum source than expected, and, further, if these sources also have a gradual response in Mg ii h and k. H. Tian et al. (2015) discuss a nearly simultaneous evolution of the intensities of Mg ii, the NUV continuum, and the FUV continuum in the impulsive phase of a solar flare. In another flare, the FUV continuum intensity in the bright ribbons appears to exhibit a more prominent spike phase than the NUV continuum intensity, which peaks after the wavelength-integrated Mg ii k light curve. To our knowledge, these differences have not been previously discussed and interpreted in terms of solar flare heating mechanisms, perhaps because it is generally thought that irradiation of the low atmosphere by transition region lines adequately explains the FUV continuum response (M. E. Machado & J. C. Henoux 1982; J. G. Doyle & K. J. H. Phillips 1992; H. Tian et al. 2015; P. R. Young et al. 2015). Solar Dynamics Observatory/Atmospheric Imaging Assembly 1700 Å (J. R. Lemen et al. 2012) data may facilitate closer comparisons than the IRIS/FUV spectra to the rising NUVA spectra in stellar flares, but these data lack spectral information and saturate in the impulsive phase of large solar flares.

The relatively gradual response of the Mg ii emission lines in the new NUV spectra (Figure 1, top) may indicate that the most prominent sources of chromospheric heating and continuum formation in stellar flares differ from bright kernels (e.g., H. Tian et al. 2015; L. Kleint et al. 2016) in solar flares. A. F. Kowalski et al. (2019) compare the Mg ii line and NUV continuum (λ ≈ 2650 Å) fluxes in two lower-energy flares. Two events in their sample show longer t_1/2 values, and one shows a delay of the peak Mg ii line flux at 60 s cadence. S. L. Hawley & B. R. Pettersen (1991) have previously discussed similarities to the Ca ii H and K line fluxes, which show delayed light-curve peaks, in the gradual decay phase of the 1985 April 12 Great Flare of AD Leo. If Mg ii h and k respond very gradually and peak after the continuum and hydrogen Balmer lines like the resonance lines of Ca ii H and K are widely observed to do in other dMe events (D. García-Alvarez et al. 2002; A. F. Kowalski et al. 2013), then the light-curve evolution in the top panel of Figure 1 is expected. To date, no self-consistent explanation exists for the delayed peaks of the Mg ii and Ca ii resonance lines in stellar flares. However, contributions from white-light "postflare" loops (P. Heinzel & K. Shibata 2018; K. E. Yang et al. 2023) and radiative backheating of large regions of the chromosphere (S. L. Hawley & G. H. Fisher 1992; S. L. Hawley et al. 2003; G. H. Fisher et al. 2012) surrounding the extended beam-heated ribbons are tantalizing possibilities.

To further argue that extreme heating to T ≳ 10,000 K occurs in the deep stellar atmosphere where large continuum optical depths can develop (e.g., A. F. Kowalski 2023), we show a suggestive resemblance of the HST/COS SED of Flare Event 2 to the observed UV spectrum of Vega (R. C. Bohlin 2014; R. C. Bohlin et al. 2014) in Figure 5. A. F. Kowalski et al. (2011) find that the newly formed flare optical spectra of secondary events in the decay phase of an M dwarf megaflare resemble the spectrum of Vega. A. F. Kowalski et al. (2017b) modeled these flare spectra also with heating from an electron beam with E_c = 500 keV, but the energy flux density is a factor of 5 smaller than the cF13-500-3 model. In Figure 5, the UV spectrum of Vega bears a few interesting similarities to the flare spectra of CR Dra: the overall shape of the Vega spectrum in the NUVB range is relatively flat, and there is a spectral break that is followed by a steeper increase at shorter wavelengths starting at λ ≈ 2400 Å. The break in Vega's spectrum is due to a confluence of Fe ii bound–bound transitions (A. García-Gil et al. 2005); we speculate that spectral syntheses of the relativistic electron-beam-heated model atmospheres with Fe ii opacity may produce flatter NUVB shapes that are closer to the observations. This additional source of NUV opacity may improve the RHD model fits, which currently overpredict the observed TESS flux in Flare Event 2.

The NUV flare spectra have important applications beyond stellar flare physics. The spectra improve the empirical and phenomenological (blackbody) interpretations of broadband UV flare data (R. D. Robinson et al. 2005; C. E. Brasseur et al. 2023; I. I. Tristan et al. 2023; V. L. Berger et al. 2024; R. R. Paudel et al. 2024), which are more readily obtained than spectra. There are several upcoming missions that will have one or two broadband UV imaging capabilities with the stated goals of characterizing stellar activity. The Star-Planet Activity Research CubeSat (D. R. Ardila et al. 2022) will have two bandpasses at λ = 1530–1710 Å and λ = 2580−3080 Å, the Quick Ultra-Violet Kilonova surveyor (J. Krtička et al. 2024; N. Werner et al. 2024) will have two bands at λ = 1400–1900 Å and λ = 2600–3600 Å, ULTRASAT (Y. Shvartzvald et al. 2024) will have one band at λ = 2200–2900 Å, and the Ultraviolet Explorer (S. R. Kulkarni et al. 2021) will have two bandpasses at λ = 1390–1900 Å and λ = 2030–2700 Å. With just one or two bandpasses, it is not possible to constrain a spectral superposition with two RHD models or two blackbody temperatures; thus, additional constraints will be necessary for accurate modeling of stellar flaring activity. Other, related applications of the new NUV flare spectra pertain to exoplanet habitability experiments (X. C. Abrevaya et al. 2020) and photochemical simulations. The FUV radiation and short-wavelength NUV photons photolyze atmospheric carbon dioxide, methane, oxygen (O₂), and water, among other gases. Thus, the correct account of the flux in this wavelength region is fundamental to predict the effect of flares in planetary atmospheric chemistry and prebiotic chemistry (A. Segura et al. 2010; F. Tian et al. 2014; R. O. P. Loyd et al. 2018a; E. W. Schwieterman et al. 2019; H. Chen et al. 2021). For example, more FUV would produce more O₃ in an O₂-rich atmosphere, while in a CO₂-rich atmosphere, more CO and O₂ would be created and CH₄ photolysis could lead to haze formation.

It has generally been assumed that a T ≈ 10⁴ K blackbody model is a sufficiently accurate representation for M dwarf flare radiation, but there has never been a precise test of the predicted peak flux with NUV spectra. We present such spectra in the impulsive phase of two megaflare events, which provide the first constraints on what happens spanning Δλ ≈ 1500 Å across the NUV. The long-wavelength NUV spectra have blackbody color temperatures of T ≈ 10⁴ K, in line with general expectations from longer wavelengths within the U band (B. Fuhrmeister et al. 2008). However, the unassuming (due to a lack of commonly studied resonance lines) wavelength range from 1700 to 2100 Å varies more impulsively and exhibits larger peak continuum fluxes and hotter blackbody color temperatures. These unprecedented findings suggest that the peak of the white-light continuum radiation is located at λ ≲ 1700 Å and a break to a rising spectrum occurs between λ ≈ 2100 Å and ≈2700 Å. For the first time, we have shown time-resolved evolution of the Mg ii emission line flux in the impulsive phase of a stellar flare. The wavelength-integrated Mg ii emission line flux evolution is almost completely decoupled from the time evolution of the continuum radiation, which dominates the integrated energy in the NUV.

These spectral properties constrain a remarkable amount of impulsive phase heating to the deep stellar atmosphere. Optically thin hydrogen recombination radiation has been used to model solar flare continuum intensities in the IRIS/NUV (e.g., P. Heinzel & L. Kleint 2014; L. Kleint et al. 2016; A. F. Kowalski et al. 2017a), which overlaps with our NUVB range, and continuum intensities around 2000 Å(M. Dominique et al. 2018), which overlaps with our NUVA range. Considered together, the new short- and long-wavelength NUV spectra provide the most compelling evidence to date that this paradigm does not explain flaring M dwarf continuum radiation (even with a modestly heated upper photosphere; D. F. Neidig et al. 1993; J. C. Allred et al. 2006; L. Kleint et al. 2016; A. F. Kowalski et al. 2017a), an idea that was pioneered by W. E. Kunkel (1970) and S. L. Hawley & G. H. Fisher (1992).

Among all current stellar flare RHD models, the one that best reproduces the continuum radiation in the new spectra at λ = 1700–2100 Å has very large heating rates from relativistic electron beams in the deep atmosphere, which results in hot blackbody-like (i.e., originating from large and wavelength-dependent continuum optical depths) thermal radiation. A superposition of two model radiative flux spectra, which may plausibly represent bright spatially extended ribbons and compact kernels, is a semiempirical explanation for two blackbody temperatures across the NUV. However, the particle beam energies in each model component are much larger than the canonical values in solar flares. The consequences of large flux densities of relativistic electron beams in stellar flares should be taken seriously in the absence of evidence of high-energy proton beams or another physical explanation for the NUV spectra.

A.F.K. acknowledges funding support from HST GO 17464 and helpful comments on the paper from Dr. John P. Wisniewski. Y.N. acknowledges funding from NASA ADAP 80NSSC21K0632, NASA TESS Cycle 6 80NSSC24K0493, NASA NICER Cycle 6 80NSSC24K1194, and HST GO 17464. We thank the Program Coordinator Bill Januszewski for help planning the HST observations. A.F.K. thanks Dr. Joel C. Allred for updating the RADYN code to the RADYN+FP code, thus allowing the proton beam simulation; Dr. Mats Carlsson for the use of the RADYN code; Dr. Suzanne L. Hawley for helpful conversations about near-ultraviolet flare spectra; Dr. Enoto Teruaki for helpful discussions about an early draft; Clara Brasseur for help on proposal preparation and discussions; Dr. Ward Howard for discussions on the timestamps in TESS headers; Elena Mamonova for helpful discussions; and Dr. Adina Feinstein for helpful discussions about FUV spectra. We make use of Paul Tol's color tables (https://personal.sron.nl/~pault/).

Some of the data presented in this article were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute. The specific HST/COS and TESS observations analyzed can be accessed via doi:10.17909/dbr7-3f98. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support to MAST for these data is provided by the NASA Office of Space Science via grant NAG5-7584 and by other grants and contracts. Funding for the TESS mission is provided by the NASA Explorer Program. The NICER data used in this study were obtained through the GO program Cycle 6 ID: 7055. NICER analysis software and data calibration were provided by the NASA NICER mission and the Astrophysics Explorers Program. This research also makes use of observations from the Las Cumbres Observatory (LCO) global telescope network. We used the LCO observation time allocated to the University of Colorado.

Software: astropy (Astropy Collaboration et al. 2022).

A.1. HST/COS Data

A.2. TESS Data

A.3. NICER Data

A.4. Las Cumbres Observatory Data

In this Appendix, we show supplementary light-curve data of Flare Event 1 (Figure 6) and Flare Event 2 (Figure 7), and we describe the calculation of energies from the TESS SAP_FLUX data. The additional light-curve data provide important contextual information for the HST/COS light curves that are discussed in the main text.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The supplementary light curves indicate that the two NUV flare events are associated with luminous soft X-ray flares. The E = 0.3–4.0 keV X-ray light-curve data (Figure 6; Appendix A.3) during Flare Event 1 show a typical slow rise and late peak after the TESS peak. The E = 0.3–4.0 keV X-ray light curve of Flare Event 2 (Figure 7) extends just past the peak of the TESS light curve. After about 1 hr into the flare, the X-ray flux is decaying but is still elevated by just as much as near the peak of the TESS light curve. The X-ray peak delays with respect to the NUV light-curve peaks are qualitatively consistent with the Neupert effect (W. M. Neupert 1968), which is widely reported in solar (e.g., B. R. Dennis & D. M. Zarro 1993; A. Veronig et al. 2002) and stellar (S. L. Hawley et al. 1995; M. Güdel et al. 1996, 2002; R. A. Osten et al. 2004, 2016; B. Fuhrmeister et al. 2011; S. Lalitha et al. 2013; M. D. Caballero-Garcia et al. 2015; I. I. Tristan et al. 2023) flares. Detailed analyses of the X-ray data will be presented in a future paper (Y. Notsu et al. 2025, in preparation).

The TESS data show the red optical and NIR extension of the white-light response in each flare, and they are presumably dominated by continuum emission (S. L. Hawley & B. R. Pettersen 1991; S. J. Schmidt et al. 2007; B. Fuhrmeister et al. 2008; A. F. Kowalski et al. 2013). The TESS light curve of Flare Event 1 (Figure 6) shows that the HST orbit ends at the termination of the fast decay phase. At the start of the next orbit, the Mg ii line emission is still elevated and continues to decay. The V-band light-curve data (Appendix A.4) reach a peak ${\rm{\Delta }}f/{f}_{\mathrm{preflare}}\approx 0.27$ and have a slightly faster decay than the TESS light curve, which peaks at ${\rm{\Delta }}f/{f}_{\mathrm{preflare}}\approx 0.05$. We note similar low-amplitude variations in the decay phases of the V- and TESS-band light curves. The HST stopped observing about 160 s into the fast UV rise of Flare Event 2 (Figure 7), after which the flaring TESS flux doubled before reaching a peak. Presumably, the UV would have also continued to rise dramatically after the end of the HST observations. TESS data at Δt = 120 s cadence are available for both flare events, but we do not show them here.

TESS SAP_FLUX, f, is in data units. We convert to ${\rm{\Delta }}f/{f}_{\mathrm{preflare}}$, where ${\rm{\Delta }}f=f(t)-{f}_{\mathrm{preflare}}$. Then we convert to energy as follows. The HST/Faint Object Spectrograph (FOS) quiescent spectrum (J.-C. Augereau & H. Beust 2006; I. I. Tristan et al. 2023) of AU Mic is a representative template of an M1e star, which we scale to the quiescent V-band magnitude of CR Dra. I. N. Reid et al. (1995) and I. N. Reid & S. L. Hawley (2005) ¹⁷ report V = 9.97 for CR Dra, but the SIMBAD value is significantly brighter (V = 9.46; J. R. Ducati 2002). The yearly averaged V and I_c magnitudes are derived from the Kamogata/Kiso/Kyoto Wide-field Survey (KWS) data. ¹⁸ The data points with large photometry errors (>0.08 mag for the V band; >0.05 mag for the I_c band) are removed, and the data within each observing season are averaged. At the time of the 2024 observations of CR Dra, the Johnson V-band magnitude was ≈10.1 (the KWS V- and I_c -band data were calibrated using V and I_c magnitudes of nearby stars in the Hipparcos main catalog; M. A. C. Perryman et al. 1997). We thus scale to a V-band filter-weighted flux density (M. Sirianni et al. 2005) corresponding to V = 10.1 using a Johnson V-band zero-point (C. N. A. Willmer 2018) and the photon-counting sensitivity curve from J. Maíz Apellániz (2006). Following J. R. A. Davenport et al. (2012), we extend the HST/FOS spectrum to the TESS bandpass by scaling a dM1e optical template (J. J. Bochanski et al. 2007) from the Sloan Digital Sky Survey and a spectrum of the M1 star HD 42581 from the NASA Infrared Telescope Facility catalog of cool star templates (M. C. Cushing et al. 2005; J. T. Rayner et al. 2009). The photon-counting sensitivity curve of the TESS bandpass gives the quiescent flux of CR Dra to be 6.1 × 10⁻¹³ erg cm⁻² s⁻¹ Å⁻¹. The flare equivalent durations (R. E. Gershberg 1972), the distance (d = 20.3 pc; Gaia Collaboration et al. 2018) to CR Dra, and the TESS bandpass FWHM (3982 Å) are used to calculate flare energies, E_TESS, and peak luminosities. An approximate transformation to energies calculated for a T = 9000 K blackbody flare, as often done in the literature, is to multiply our bandpass energies by a factor of ≈5; for a T = 10,000 K blackbody, the factor is ≈6. The flare-only flux of Flare Event 2 (Figure 5) is calculated by multiplying the value of ${\rm{\Delta }}f/{f}_{\mathrm{preflare}}=0.09$ by the TESS bandpass quiescent flux above.

Several calculated quantities from the TESS light curves are summarized in Table 1.

Table 1. Calculated Quantities from TESS (Δt = 20 s) Light Curves

Flare	${f}_{\mathrm{preflare}}$	Peak ${\rm{\Delta }}f/{f}_{\mathrm{preflare}}$	t1/2	Impulsiveness a	Equivalent Duration	Total Duration	e−1 Duration
	TESS DN		(s)	(10−5 s−1)	(s)	(s)	(s)
Flare Event 1	84354	0.051	911 b	4.5–5.6	62	5500	880
Flare Event 2	89163	0.734	650	113	757	17,000	720

Notes. t_1/2 is the FWHM of the light curve, and e⁻¹ duration is the elapsed time between the peak flux time and the time when the flux reaches a value of peak flux divided by e (H. Maehara et al. 2015; K. Namekata et al. 2017).

^a Impulsiveness index calculated as peak I_f divided by t_1/2. The total durations are roughly determined according to the time when the gradually decaying flux returns to the level of the preflare flux, f_preflare, within the scatter. ^b Including both peaks in Flare Event 1 gives a full width at 0.45 of the maximum of 1122 s.

Download table as:  ASCIITypeset image

This Appendix supplements Figures 2 and 3 with other flare spectra and their corresponding blackbody and RHD model fits at several interesting times during Flare Events 1 and 2. These demonstrate how the best-fit blackbody temperatures and RHD filling factors differ compared to the integrated spectra over the second major peak in Flare Event 1 and over the entire fast rise phase interval that was observed in Flare Event 2. For example, the spectrum from the last 30 s before HST stopped in Flare Event 2 is shown in the bottom panels of Figure 8 and Figure 9. The blackbody fit to NUVA has a color temperature of nearly ${\hat{T}}_{{\rm{NUVA}}}\,\approx \,\rm{19,000}$ K, and the value of ${\hat{X}}_{\mathrm{rel}}$ is ≈1.2. The satisfactory fits of the two-component RHD model in the second major peak of Flare Event 1 (Figure 9 (top)), when the NUVA flux was relatively fainter than in Flare Event 2, demonstrates the versatility of this approach. The middle panel of Figure 9 demonstrates the fitting over the time interval corresponding exactly to the last six data points of TESS before HST stopped observations (${\hat{X}}_{\mathrm{rel}}=4.5;$ Figure 5). The spectrum that is summed over the impulsive phase of Flare Event 1, consisting of two major peaks and their respective fast decay phases, is shown in Figure 8 (middle). This flare spectrum has the highest signal-to-noise ratios, but the relative fluxes in NUVA and NUVB vary considerably over the course of the impulsive phase (Figure 1). The flare spectrum over the first fast rise phase of Flare Event 1 (Figure 8) does not have any flux contributing from a previously decaying peak; the NUVB is remarkably flat, and the NUVA has a rather large blackbody color temperature of ${\hat{T}}_{\mathrm{NUVA}}\approx$ 16,600 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size