© The Authors 2026

1. Introduction

Exoplanet science is rapidly advancing towards the discovery of Earth-like planets (i.e. in terms of size and mass), situated within the liquid-water habitable zone (Owen 1980). Since the early 2010s, studies by Gaidos & Mann (2013) and Dressing & Charbonneau (2015) have shown that planets in the liquid-water habitable zone are more common around M dwarfs. Consequently, M dwarfs have emerged as the most prevalent type of star known to host Earth-like planets. It is estimated that approximately 40% of M dwarfs could harbour one or more planets within their habitable zone (Bonfils et al. 2013; Mulders et al. 2015a,b). The relative abundance of M stars makes them prime targets in the search for habitable worlds.

The question of whether M stars can sustain habitable planets, however, remains open, partly due to their strong chromospheric activity and frequent flare production (Segura et al. 2010). While flaring activity may alter and erode planetary atmospheres (Lammer et al. 2007; Segura et al. 2010; Tilley et al. 2019), it may also provide the necessary ultraviolet (UV) radiation for prebiotic chemistry processes (Ranjan et al. 2017; Rimmer et al. 2021).

The presence of self-sustained magnetic fields in low-mass stars significantly impacts their upper atmospheric structure and high-energy radiative environments. Various heating processes, including shock dissipation of acoustic waves generated in stellar surface convection zones (Narain & Ulmschneider 1996) and magnetic reconnection (Klimchuk 2006), cause these low-mass stars to exhibit significant temperature inversions in their outer atmospheres. The increased power and frequency of M dwarf flares may be attributed to their turbulent magnetic dynamos and strong magnetic fields. It remains uncertain whether the mechanisms that produce M dwarf flares are truly analogous to those of the Sun. It is believed that stars below approximately 0.35 M_⊙ become fully convective, lacking a radiative zone and consequently a tachocline (Chabrier & Baraffe 1997). This transition to full convection may involve a shift from a Solar-type magnetic dynamo governed by rotation to a turbulent dynamo that is potentially less dependent on rotation (Reiners et al. 2022). Despite extensive research on stellar chromospheres and coronae, the nature of these layers, the underlying processes that generate them, and their evolution remain poorly understood, particularly for M-dwarf stars (Pineda et al. 2021).

Flares, manifesting as broadband emission ranging from the near-UV (NUV) through optical wavelengths and, in some cases, extending into the far-UV (FUV) and near-infrared (NIR), have been extensively studied in recent years. This progress is largely due to space-based dedicated photometric monitoring missions such as Kepler (Borucki et al. 2010), its extension K2 (Howell et al. 2014), and the ongoing Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015).

Flares and associated energetic particle events have the potential to alter the chemistry and climate of exoplanetary atmospheres, particularly those with Earth-like characteristics (Lammer et al. 2007; Segura et al. 2010; Venot et al. 2016; Tilley et al. 2019; Herbst et al. 2019, 2024). However, these analyses were primarily based on observations of a limited number of well-characterised M dwarf flares at FUV wavelengths and extrapolations from solar and M dwarf flare activity observations at optical wavelengths, since there are significant challenges associated with direct UV and X-ray observations. Below the hydrogen ionisation edge at 912 Å, stellar emission is absorbed by the interstellar medium, making UV and X-ray light observable only from space.

The question of whether M dwarfs can host habitable planets has attracted significant attention recently. For a detailed review of the numerous factors influencing the habitability of planets orbiting G, K, and M dwarfs, see Airapetian et al. (2019), Engle (2024).

M dwarf flares are, on average, more energetic than solar flares (Hawley & Pettersen 1991; Davenport 2016; Namekata et al. 2017). Typical white-light energies of active M dwarf (dMe) flares are 10³⁰–10³² erg (Hawley et al. 2014) with occasional events exceeding 10³⁴ erg in the U band (Kowalski et al. 2015). The areal extent of flare footpoints cannot be measured directly, and hard X-ray emission is far too faint except during the largest dMe flares (Osten et al. 2010). Although Kowalski et al. (2015) did not determine if non-thermal (NT) electron fluxes on M dwarfs are indeed as high as 10¹³ erg cm⁻² s⁻¹, typical values of ∼10¹¹ erg cm⁻² s⁻¹ in strong solar flares cannot explain the properties of white-light emission during dMe flares. Thus, it is plausible to assume that the energy flux is higher.

Radiative-hydrodynamic (RHD, Allred et al. 2005; Kowalski et al. 2017) models grounded in physical principles have become essential for understanding the dynamic processes in stellar atmospheres. By coupling radiation and hydrodynamics, these models allow for self-consistent modelling that captures the complex interplay between energy transport, atmospheric structure, and observable phenomena such as flares. Validated against both solar and stellar observations, RHD models provide a robust framework for interpreting a wide range of stellar activity and its signatures across different spectral regimes when subsequently used in radiative transfer calculations to produce observable signatures (Froning et al. 2019; Kowalski 2022).

Spectroscopic studies (Froning et al. 2019; Howard et al. 2023) alongside broadband analyses (Paudel et al. 2024; Brasseur et al. 2023) have begun incorporating multi-component models spanning from FUV to NIR. The multi-wavelength approach of fitting the combinations of the RHD models, combining TESS optical data and the FUV observations from the Cosmic Origins Spectrograph on the Hubble Space Telescope (HST-COS), provides a comprehensive view of magnetic activity in young M dwarfs.

In this paper, we discuss a newly developed Young M dwarf flare (YMDF) model to study exoplanetary atmospheres in the vicinity of young and active M dwarfs. The paper is organised as follows: in Section 2 we describe our methods for developing the model components, as well as the flare sample from TESS and FUV observations used for evaluating synthetic light curves and spectra. In Section 3 we present the results obtained using the YMDF model and analyse the synthetic stellar activity datasets. A comparison of our proposed model and previous modelling efforts is discussed in Section 4.

2. Methodology

2.1. Stellar sample

In Mamonova et al. (2025) (hereafter, Paper I), we assembled a sample of 493 young and field M-K dwarfs and analysed their activity patterns, validating 86 714 flare events in TESS and Kepler data and 52 events in HST-COS FUV data. In this work, we utilised a portion of this sample, specifically the observational data of stars exhibiting flaring activity in both FUV and TESS bandpasses (see Table 1).

2.2. Stellar atmosphere models during a flare event

Stellar flares are believed to originate from processes similar to those observed on the Sun. In solar flares, NT electrons generate hard X-ray sources that are co-temporal and co-spatial with white-light radiation at the foot-points of reconnected magnetic loops. Thus, electron beam models, analogous to the ones used in this study, provide a physically grounded representation of the energy input in flares. The beam approach with RHD simulations models key stellar atmospheric properties such as opacities, flows, heating rates, and non-equilibrium radiation. RADYN assumes an impulsive flare phase driven by a downward-propagating charged particle beam accelerated in the corona, rapidly heating and expanding the lower atmosphere (i.e. chromospheric evaporation, Allred et al. 2015). This allows flare evolution modelling from impulsive to decay phases. Kowalski et al. (2024, 2017) successfully reproduced observed flare properties, including optical and NUV continuum and Balmer emission lines, validating the electron beam framework as a robust stellar flare model foundation.

To model stellar flare heating, a 1D time-dependent model employs a power-law distribution of accelerated electrons (Kowalski et al. 2024), coupled with the RADYN RHD code (Carlsson & Stein 1994, 1995, 1997). RADYN self-consistently solves the equations of mass, momentum, and internal energy conservation, as well as non-equilibrium level populations, charge conservation, and radiative transfer (Carlsson & Stein 1994). The non-local thermodynamic equilibrium (non-LTE) radiative transfer problem is addressed within RADYN using the accelerated Λ-iteration techniques developed by Scharmer (1981) and Scharmer & Carlsson (1985). The electron beam enables significant heating at deeper atmospheric levels. By adjusting the parameters of the electron distribution, such models achieve optical blackbody temperatures of approximately 10 000 K in the emergent radiative flux spectra, consistent with observational data. Parameters that describe the electron beam power-law distribution, such as the low-energy cutoff (E_c), power-law index (δ), and energy flux, govern the energy deposition and heating depth during flares.

Kowalski et al. (2024) suggested time-dependent stellar flare models of deep atmospheric heating in M dwarfs. In their grid the constant group includes models with large E_c and constant flux injections reaching 10¹³ erg cm⁻² s⁻¹, producing intense heating and continuum emission. The auxiliary group expands the grid with models that better reproduce observed features than the main group, comprised of impulsive beam models. The auxiliary models can have a lower E_c and a hard power-law index δ = 2.5, representing a relatively low flux but energetic beam, reaching peak flux near 10¹² erg cm⁻² s⁻¹ in a pulsed ramp up and down injection, analogous to solar flare hard X-ray bursts characterised by rapid increases followed by gradual decreases in intensity.

Kowalski et al. (2025) found that only the most energetic electron beam heating model can reproduce the rising slope of flares they observed within the NUVA wavelength range, and in their framework, a lower E_c provided an adequate explanation for the overall properties of the full UV continuum. They studied a recent super-flare in CR Dra, an M dwarf with properties broadly similar to AU Mic, although CR Dra has an uncertain age and a rotation period twice as fast. Thus, we adopted their approach: a linear combination of the two atmospheric models during a flare event, each scaled by its respective coefficient. The first atmospheric model (from the auxiliary group), m2F12-37-2.5, is characterised by a low flux density, an electron beam of 2 × 10¹² erg s⁻¹ cm⁻², a smaller low-energy cutoff E_c = 37 keV, and a hard power-law index δ = 2.5. The second, high flux density model (from the constant group), cF13-500-3, has an energy flux density of 10¹³ erg s⁻¹ cm⁻², a low-energy cutoff of E_c = 500 keV, and a number-flux power-law index of δ = 3 above this low-energy cutoff. The physical origin of these two models can be attributed to fainter ribbons (associated with the low-flux m2F12-37-2.5 model) and bright kernel or kernels (associated with the high-flux cF13-500-3 model) in the flaring portion of the stellar surface, analogous to those observed in solar flare studies (e.g. (Kowalski 2022, Fig. 6)).

To obtain the flux density of the models during flare events, we used the RH 1.5D code (Pereira & Uitenbroek 2015), which solves the non-LTE radiative transfer problem. The code employs the multi-level accelerated Λ-iteration method developed by Uitenbroek (2001), building on the approach of Rybicki & Hummer (1991) to handle partial frequency redistribution (PRD) effects. To efficiently solve the complex radiative transfer and atomic rate equations, RH 1.5D applies advanced numerical techniques that speed up convergence and improve solution stability, along with simplifying the system of the rate equations that govern atomic level populations, enabling the treatment of spectral line blends.

In addition to supporting overlapping radiative transitions, the RH 1.5D code solves the rate equations for multiple atoms simultaneously, consistently treating any overlapping transitions (Pereira & Uitenbroek 2015). It provides flexibility in the treatment of atoms by allowing transitions and continua to be computed either in non-LTE for active atoms (i.e. populations of atoms treated in detail and updated according to non-LTE radiative transfer) or under the assumption of LTE for passive atoms, treating them as opacity sources.

To run RH 1.5D, spectral setups have to be determined. We used a six-level hydrogen atom as the only non-LTE atom, and used the LTE conditions for Ca II, Mg II, Si, Al, Fe, He, N, Na, S atoms (hereafter, the H6 spectral setup). In addition, the maximum number of PRD iterations per main iteration was set to ten, and the convergence limit of PRD iterations in each main iteration is 1.0 × 10⁻². We also created a wavelength table from 10 to 55 000 Å in order to calculate F_ν for most frequencies. For further analysis, we also included C I, C II, and C III atoms in LTE as background opacities, in addition to the H6 spectral setup’s species, which we refer to as the H6CIII spectral setup hereafter.

2.3. YMDF model components

We developed the YMDF model² to reproduce the time-dependent spectra of flaring young M dwarfs. This module predicts the spectral evolution of individual flares utilising publicly available radiative outputs from the grid of RADYN stellar flare simulations (Kowalski et al. 2024; Kowalski 2022). These build upon foundational work by Carlsson & Stein (1992, 1995), Allred et al. (2015) and Carlsson et al. (2023).

To calculate the flux integrated over the TESS wavelength range of the early M dwarfs in our sample, we used the linear two-component model (Kowalski 2022; Kowalski et al. 2024, 2025):

$$\begin{array}{c}\hfill f_{\mathrm{obs}}^{\prime }=\frac{R_{★}^2}{d^2}({\widehat{X}}_1{F}_1+{\widehat{X}}_2{F}_2),\end{array}$$(1)

where R_★ is the radius of the flaring star observed by TESS, and d is the star distance from the Gaia astrometric space mission recent data release (GAIA DR3, see Fouesneau et al. 2023). The fluxes F₁ and F₂ are radiative surface fluxes derived by the corresponding models with the pre-flare flux subtracted. ${\widehat{X}}_1$ and ${\widehat{X}}_2$ give the best-fit filling factors of both RHD model components, representing the contribution of each model flux in the resulting flare flux. These filling factors represent a physical area comprising multiple flare-related processes occurring simultaneously. A constant electron beam in the cF13-500-3 model and a pulsed electron beam in the m2F12-37-2.5 model were treated as time-resolved phenomena with durations of several seconds. These components were averaged over a specific time interval in the models (0–9.8 s and 0–4.2 s for the m2F12-37-2.5 and cF13-500-3 models, respectively) to reflect multiple such events within the flare region, thereby capturing the spatio-temporal complexity of the flare evolution. To obtain the spectra at the peak of the flare, we co-added time-resolved spectra from the RH 1.5D. For this, we ran the RH 1.5D in the H6 and H6CIII spectral setups separately. From the flare flux at the peak, we obtained a temporal evolution of the flare employing several temporal flare models widely recognised in the literature, and we further discuss them in Sect. 2.4.

Due to the linearity of the integration, the model enables two key applications: i) generating synthetic flux densities over extended wavelength ranges at the peak of the flare, and ii) calculating integrated fluxes within specific spectral band passes (FUV and TESS). We calculated F₁ and F₂ flux densities integrated over the TESS wavelength range (i.e. ∫R_res(λ)F_λ, 1dλ and ∫R_res(λ)F_λ, 2dλ, where R_res(λ) is the TESS response function) from m2F12-37-2.5 and cF13-500-3 atmosphere models using RH 1.5D (see for details Paper I, Eq. (3)). To obtain the best fit of ${\widehat{X}}_1$ and ${\widehat{X}}_2$, we performed a systematic parameter search for ${\widehat{X}}_1$ and ${\widehat{X}}_2$ using Markov chain Monte Carlo (MCMC) sampling with the ‘emcee’ package (Foreman-Mackey et al. 2013). For FUV spectral fitting, we retained the established methodology while omitting instrumental response corrections (see Paper I, Appendix A.5).

2.4. Synthetic distributions

Given that the equivalent duration (ED) serves as a proxy for flare energy, the derived model coefficients can be linked to ED values, thereby enabling amplitude estimations for flares of specified EDs. The ED quantifies flare energy relative to the star’s quiescent luminosity through the integral relationship:

$$\begin{array}{c}\hfill ED={\displaystyle \int (\frac{F_{\mathrm{flare}}(t)}{F_{\mathrm{quiescent}}}-1)}dt,\end{array}$$(2)

where F_quiescent denotes the quiescent stellar flux and F_flare(t) represents the time-dependent flux during the flare event with duration t. The ED serves as a direct proxy for the specific flare energy when scaled by the star’s quiescent luminosity in a given wavelength range.

The temporal evolution of the flare flux, F_flare(t), can be reconstructed using several established models from the literature. Our Python implementation includes three approaches: the temporal flare model described by Tovar Mendoza et al. (2022), the framework introduced by Davenport et al. (2014), and a modified Gaussian rise and exponential decay model based on the work of Feinstein et al. (2024). Hereafter, we refer to these as the Mendoza, Davenport, and Feinstein models, respectively.

These single-flare time-resolved models require three principal parameters: amplitude, duration, and peak times. We first implemented the single-flare model with the amplitude equal to unity and arbitrarily chose the peak time. We sampled duration randomly from the distribution explained in Yang et al. (2023) and Zhang et al. (2024). They identified a distinct power-law correlation between flare duration D and released energy E, characterised by D ∝ E^γ, for events exceeding 10-minute durations. Below this temporal threshold, the relationship bifurcates, revealing that short-duration flares (D < 10 min) can exhibit energy release comparable to their longer-duration counterparts. For the M stars, the power-law index was constrained to γ = 0.3 (Yang et al. 2023, Fig. 14). The resulting normalised flux can be scaled by the peak flux amplitude derived from the YMDF model. By adding the quiescent stellar flux, the complete temporal evolution of the surface flux is produced. Integrating the model spectra over a specified wavelength range allows fitting of the resulting light curve to observations, such as the TESS white-light curve or the light curve obtained by integrating over the HST-COS G130M grating wavelengths. Subsequently, the model spectra can be compared to spectrally resolved HST-COS FUV data at specific time points, such as pre-flare or flare peak.

In Paper I, we proposed that a piecewise power law provides a more accurate representation of observed flare energy and EDs distributions in young M dwarfs. The YMDF model implementation enables both the sampling of EDs from a broken power law as well as single power-law distributions. For every flare in the simulated distribution, EDs are used to determine the coefficients of the m2F12-37-2.5 and cF13-500-3 atmosphere models, enabling the simulation of the flare’s spectral-temporal profile. This results in complete synthetic spectra of a star evolving for prolonged periods of time. These spectra can subsequently serve as inputs for exoplanetary atmospheric studies, providing critical data for chemistry kinetic models across a broad wavelength range from 10 to 50 000 Å.

3. Results

3.1. Light curve fitting and analysis.

We analysed the flare data in the TESS bandpass in our sample of 11 young M-K dwarfs by fitting the two-component model we described in Section 2.2. The TESS bandpass (∼6000–10 000 Å) centred on the Cousins I band, is largely devoid of prominent atomic and ionic emission lines, being dominated by the stellar continuum and a few broad molecular features (Ricker et al. 2015). This red-optical to NIR coverage was chosen for the TESS survey to maximise sensitivity to small planets transiting cool, red stars, but it lacks the rich array of diagnostic lines characteristic of the FUV, resulting in a fundamentally different spectral information content. Fitting the simplified model based on the H6 spectral setup yields a ratio between the low energy m2F12-37-2.5 and the high energy cF13-500-3 stellar atmosphere models of R = 5.03 ± 0.01. Both coefficients, ${\widehat{X}}_1$ and ${\widehat{X}}_2$, from Eq. (1) were free parameters while using MCMC sampling to fit to the observed data. Figure 1 presents the results of fitting ${\widehat{X}}_1$ and ${\widehat{X}}_2$ to observational data from TESS (left panel). We further tested the sensitivity of these wavelength ranges to the nomenclature of lines, using the H6CIII spectral setup in the RH 1.5D calculations. This inclusion shifted the mean ratio to R = 8.74 ± 0.02 for TESS range flares.

Fig. 1. Coefficients of the fitted models for flares in young M dwarfs. The values of ${\widehat{X}}_1$ and ${\widehat{X}}_2$ are plotted on the y- and x-axes, respectively. Each point represents an individual flare event and is colour-coded by the flare energy released in the TESS (left) or the FUV bandpass (right). Model fitting was performed using multiple setups: H6 and H6CIII for the TESS data, and COScut1 and COScut2 for the HST-COS data. Linear regressions fitted to the flare events for each star are shown as dashed coloured lines. To guide the eye, in the panels, the H6, H6CIII, COScut1 and COScut2 ratios are shown as dotted red, dotted black, dash-dotted grey, and dotted grey lines, respectively.

We proceeded to analyse flare data in the FUV in the same way as for TESS data. The 1040–1350 Å wavelength range is densely populated with spectral lines, making it a particularly rich region for stellar diagnostics. Among the most significant features in the FUV interval are strong emission lines (C II, C III, Si III, Si IV, N V), and the prominent and ubiquitous hydrogen Ly-α line at 1215.67 Å (Linsky 2017; France et al. 2013). Fitting the model based on the simplified H6 spectral setup is therefore a challenging task. To address these difficulties, we implemented the Kurucz line list (Kurucz & Bell 1995) in the RH 1.5D calculations in the range of 1083.990–1359.275 Å for species stated in Feinstein et al. (2022), Table A1) and the H6 spectral setup species. Additionally, we adjusted the wavelength window for integration of the flux density, resulting in two regimes. The first covers λ∈ [1060, 1360] Å, excluding the interval λ∈ [1210, 1225] Å to mask the Lyα emission. The second spans λ∈ [1170, 1430] Å. The HST-COS G130M grating configurations are discussed further in Appendix A.4. Hereafter, we refer to these spectral setups as COScut1 and COScut2, respectively. We used the time-tagged spectra (the detailed methodology can be found in Paper I) and integrated flux in FUV-B and FUV-A sectors for every timestamp. In these analyses, we produced light curves with 10 sec cadence to reach better precision in the fitting. As before, both coefficients, ${\widehat{X}}_1$ and ${\widehat{X}}_2$, in Eq. (1) were treated as free parameters in the fitting while performing the MCMC sampling. We plotted them in Fig. 1 (right panel) for the stars with a flare count ≥2 in FUV and the corresponding spectral setups.

For each star in the sample, we determined the best-fit coefficients, using TESS photometry data of flare events with ED > 15 s and all available data in FUV. The number of flares used in the analyses is stated in Table 1. The ratio between the fitted coefficients is indicated in each panel of Fig. 1 for all stars in the sample, and, for the corresponding spectral setups, the mean values of the ratios ${\widehat{X}}_1$/${\widehat{X}}_2$ are included in Table 2.

While the component ratio is sensitive to the spectral setup, our experimenting further with incorporating additional emission lines in the COScut1 and COScut2 setups moves the initial ratio values R∼ 11–13 closer to R∼ 7–8, consistent with TESS results and prior NUV analyses (see e.g. Kowalski et al. 2025). While including an exhaustive spectral line list in the radiative transfer RH1.5 calculations should, in principle, produce more reliable results, such calculations are computationally demanding and often suffer convergence issues. Overall, while the H6 model may be considered oversimplified, and the H6CIII model presents a ratio more consistent with previous studies, both models yield a satisfactory fit to the time-integrated flux exhibited by flares in TESS and FUV light curves, as we demonstrate below in Sect. 3.2.

3.2. Temporal evolution of individual flares

Previously, we used bulk flare data to fit model coefficients and determine the median ratio between low- and high-energy electron beam models. We now focus on individual flare light curves of different energies, temporal evolution, and duration, observed in TESS and FUV for stars in our sample. The morphology found in the observed flares varies largely, both in the TESS (Davenport et al. 2014; Yudovich et al. 2025) and FUV (Loyd et al. 2018a,b; Froning et al. 2019; Duvvuri et al. 2025) bandpasses. We plotted the temporal evolution of 16 flares in our sample in Fig. A.1 and the observational data for these flares are stated in Table A.1 for TESS and Table A.2 for FUV flares.

We further selected several observed flares and compared their fits with the flux obtained from models using inverse calculations. First, we determined the coefficients of the model components m2F12-37-2.5 and cF13-500-3 by utilising the observed ED and the ratios derived from our model setups. This involves integrating the contributions of each component over the corresponding wavelength ranges for TESS/FUV (F₁ and F₂), as well as integrating the pre-flare flux F₀ within the model framework, and substituting these values into the ED equation (see Eq. (2), full derivation can be found in Appendix A.3). Second, the resulting coefficients allowed calculation of the peak flare flux (Eq. (1)). Third, we employed a temporal model (the Mendoza model here) to scale the peak flux across all time stamps, thus capturing the temporal evolution of the spectrum. Throughout this work, we refer to these temporally and spectrally resolved model outputs as the ‘inverse models’.

Results for three TESS flares and three FUV flares are shown in Fig. 2 in the upper and lower rows, respectively. Observed data (red spheres with black error bars) are plotted along with the temporal flare models, illustrated here by coloured lines: Mendoza (light blue), Davenport (light coral), and Feinstein (teal). Inverse models of the H6 and H6CIII setups, derived using EDs from observed data, are shown as solid blue and green lines with circle markers, respectively. It is evident that the temporal evolution of flares in our sample is not fully captured by these models: the Mendoza model, characterised by rounded peaks and slow rise and decay, underperforms; the Davenport model exhibits a somewhat sharper peak but still fails to reproduce the full amplitude; and the Feinstein model reproduces peak amplitudes in many cases but inadequately models the slow decay. However, the simplified models can closely reproduce the emission observed by TESS and by COS-HST. We calculated time-integrated fluxes in HST-COS FUV (erg per square centimeter) and TESS (electron counts), reporting these values for all models and observed data in Table A.3 (TESS) and Table A.4 (FUV). We found that the Mendoza and Davenport models produced time-integrated fluxes similar to the observed values, whereas the Feinstein model underestimated this value. The simplified inverse models H6 and H6CIII yielded results comparable to those obtained by directly fitting coefficients to observed data and, in some cases, outperformed the Mendoza model, which most closely reproduces the time-integrated observed flux, in the entire dataset.

Fig. 2. Light curves of several flares observed in TESS and HST-COS. Upper row: TESS flux light curves (e− s−1) for flares denoted as AU Mic T2, HIP 107345 T1, and 2MASS J02365171-5203036 T1, as dashed lines with red spheres for 20-second observation stamps and black error bars with capped ends. Observation details are provided in Table A.1 in Appendix A.2. Note that errors can be smaller than the sphere size. Three temporal flare model MCMC fitting results are shown: Mendoza (light blue), Davenport (light coral), and Feinstein (teal) models fitted to the observed data. The H6 and H6CIII models, which were inversely calculated using observed EDs for these particular flares, are shown as solid blue and green lines with data-point spheres. In the upper row, the H6 model is hidden by the H6CIII due to overlap as their flux values are very similar. Bottom row: FUV HST-COS flux observations for flares denoted as AU Mic F1, AU Mic F2, and Karmn J07446+035 F1, plotted with consistent styles and model fits, as the upper row. The details of observations for these flares can be found in Table A.2 in Appendix A.4.

3.3. Model comparison to spectrally resolved data.

RHD models of chromospheres predict the formation of a hot, dense condensation layer as a result of electron beam heating. When the energy flux of the electron beam is sufficiently high, the flaring regions reach continuum optical depths of approximately τ ∼ 1, producing spectra characterised by temperatures around T ∼ 10⁴ K and exhibiting a small Balmer jump. These features are consistent with the optical flare signatures observed in dMe stars (e.g. Hawley & Fisher 1992; Kowalski et al. 2013, 2016). Froning et al. (2019) employed these parametrisations to explore the modifications required in RHD models to reproduce the hot, blue FUV continuum emission detected during the GJ 674 flare. While RHD models with very high electron beam energy fluxes (10¹³ erg cm⁻² s⁻¹) can generate the ∼10⁴ K thermal component observed in previous flares, they are unable to account for the blue FUV continuum, as the model spectra turn over in NUV (Kowalski et al. 2015).

In Fig. 3, we plotted a comparison of the surface flux as a function of the wavelengths at the peak of the flare, highlighting the performance of the two stellar atmosphere models. The m2F12-37-2.5 and cF13-500-3 models, shown as dark and light-green lines, respectively, were plotted assuming that the entire stellar surface is flaring. Although this scenario is unphysical, it enables a direct visual comparison of the spectral energy distributions (SEDs) predicted by each model across a broad wavelength range. We plotted here one example of the resulting spectra at the peak of the flare for AU Mic. The fitting procedure was performed for a median size flare, yielding ED_model = 43.55 s and a total energy release of 4.13 × 10³³ erg in the TESS bandpass. The corresponding fit is plotted as the solid red line. The obtained coefficients are used to produce the model spectrum at the peak of the flare in the H6 setup covering the entire model range (λ∈ [10;55 000] Å, shown as the dashed orange line). We determined ${\widehat{X}}_1$ and ${\widehat{X}}_2$ of the inverse H6 model from observed ED and the R = 5.03. It is shown as a dash-dotted blue line, with the pre-flare state as a dash-dotted light-blue line. Blackbody spectral models with temperatures of 3850 K, 6000 K, 8000 K, and 10 000 K are included for reference. The panchromatic flux of AU Mic (Feinstein et al. 2022) is shown in light grey, empirically benchmarking for the model predictions. The additional panels focus on the FUV and TESS wavelength ranges, where the differences in model components’ behaviour and the correspondence to the observed fluxes are most pronounced. For this medium-sized flare, the right inset (TESS bandpass) shows a modest increase in the continuum, whereas the left inset (FUV bandpass) reveals a pronounced continuum enhancement.

Fig. 3. Observed and modelled surface flux density. The m2F12-37-2.5 (dark green) and cF13-500-3 (light green) stellar atmosphere models are plotted at coefficients ${\widehat{X}}_1=1$ and ${\widehat{X}}_2=1$ (the theoretical assumption of the entire stellar surface exhibiting flaring activity). For AU Mic during a flare, the solid red and dashed orange lines show the fitted model spectrum at flare maximum using coefficients from the example fit (the flare data are stated in Table A.1 in Appendix A.2). The inverse H6 model at flare maximum, pre-flare, and the ‘fiducial flare’ model from Loyd et al. (2018b) with a comparable ED is plotted as dash-dotted blue, dash-dotted light-blue, and dashed purple lines, respectively. Blackbody curves for 3850 K, 6000 K, 8000 K, and 10 000 K are included in solid, dotted, dashed, and dash-dotted black lines. The panchromatic AU Mic flux from Feinstein et al. (2022) is plotted in light grey. Insets highlight the FUV (lower left) and TESS (lower right) wavelength ranges.

We further analysed the spectra of the sample stars by fitting the derived model-dependent parameters for several spectral setups. Feinstein et al. (2022) identified the continuum regions in the AU Mic spectrum, specifying these areas as being free of emission features. We adopted this approach to extract these regions from each time-tagged observed flux, subsequently re-binning the flux into 0.05 Å bins and further smoothing the signal by applying a rolling average over every 20 values, reaching an effective resolution of 1 Å. We experimented with various methods of binning data and found this to be the optimal approach for effectively revealing the continuum. The H6 and H6CIII spectral setups have a native resolution of 1–10 Å, while the COScut setups’ native resolution is 0.5–20 Å. All setups treat key spectral lines with enhanced resolution; therefore, the applied binning and smoothing provide an effective intermediate resolution.

In Fig. 4, we compare the obtained continuum spectra at the peak of the FUV flare in AU Mic. In order to extract the observed continuum flux, we re-binned the flux at the timestamp of 17660.5 s from the start of observation for the flare peak and at the timestamp of 16700.0 s for the quiescent continuum (plotted here as green spheres and grey diamonds, respectively). To interpret the continuum enhancement, we plotted the fitted COScut1 setup at the flare peak. For additional context, we included comparisons with the inverse models H6 with R = 5.03 and H6CIII with R = 8.74 calculated from the observed flare ED.

Fig. 4. Flare spectra in FUV HST-COS. Green circles and grey diamonds indicate the continuum flux at the flare peak and pre-flare, respectively, with uncertainties displayed as corresponding light-coloured lines. Only the upper error bars are shown for visibility, as the errors are symmetric. The best-fit model to the continuum fluxes, derived using the COScut1 spectral setup, is plotted as a light-red line. The dashed blue, dash-dotted green, and dash-dotted light-green lines correspond to the H6 and H6CIII setups at the peak, and pre-flare H6, respectively. The light-grey, light-orange, and dashed purple lines represent the sum of the extracted 1D spectra from all exposures in the orbit, the smoothed panchromatic spectrum of AU Mic from Feinstein et al. (2022), and the ‘fiducial flare’ model from Loyd et al. (2018b) at a similar ED, respectively.

We note that the COScut1, H6 and H6CIII setups all result in a good-to-eye fit in FUV-A segment (1225–1360 Å) and partially in FUV-B (1130–1210 Å), indicating that these models can potentially reproduce the continuum rise during the flare. However, the continuum rise observed in the blue short-wavelength region (∼1110–1120 Å) of the FUV-B segment remains unaccounted, and the rise in the H6CIII peak and pre-flare models is shifted to wavelengths shorter than ∼1110 Å. This rise observed in the blue FUV-B continuum has been detected and analysed in detail by Feinstein et al. (2022). We continue this discussion in Sect. 4.1.

Our attempt to fit the pre-flare flux to the observed FUV data was unsuccessful for two stars, both of which exhibit comparatively weak flare events, making accurate fitting particularly challenging. First, in 2MASS J18141047-3247344 (HD 319139), a spectroscopic binary consisting of K5 and K7 dwarfs (Nefs et al. 2012), the YMDF H6 model underestimated the pre-flare quiescent flux by two orders of magnitude. This discrepancy may be attributed to the generally lower levels of magnetic activity observed in K dwarfs compared to M dwarfs, as well as the complex nature of the binary system; the fit fails as expected since the model assumes a single star, while in reality, the system consists of two unresolved stars. In 2MASS J11173700-7704381, the youngest star in our sample at ∼5 Myr old, we observed a similar issue, with the YMDF model significantly underestimating the pre-flare quiescent flux. To better understand how the model could replicate the flares in a range of M-K dwarfs, additional observations in the FUV and NUV are needed. Such data would help us more accurately characterise the pre-flare emission and flare properties of young, low-mass stars, providing further constraints for theoretical models.

Interestingly, for Karmn J07446+035 (YZ CMi), we achieved a successful fit in FUV-A (see Fig. A.2, left panel, in Appendix A), despite it being the dimmest and oldest star in our sample. Similarly to AU Mic, in YZ CMi, we observe a continuum rise in the blue end of the FUV-B segment (approximately ≥1130 Å) comparable to features reported in AU Mic flares (Feinstein et al. 2022; Redfield et al. 2002), as we discuss further in Sect. 4.1. YZ CMi, an M4.5 dwarf, is approximately 50 Myr old, and is thought to be near the transition line between partially and fully convective regimes (see e.g. Chabrier & Baraffe 1997, 2000). Our model may slightly overestimate the pre-flare quiescent flux, but overall the agreement is good. YZ CMi is well known for its high level of magnetic activity, with flares observed across multiple wavelength regimes over many years (Roques 1961; Kowalski et al. 2017; Namizaki et al. 2023). This extensive observational history makes it a valuable benchmark for testing flare models in active M dwarfs. These results indicate that the YMDF model, though presently tested on young early M dwarfs, may also be applicable to stars with properties beyond those in our current sample.

3.4. FFDs in simulated populations.

In order to simulate star flare activity, the YMDF module: i) simulates the ED distribution based on a single or piece-wise power law; ii) produces randomly chosen durations from a bifurcated distribution (Yang et al. 2023; Zhang et al. 2024) based on given ED; iii) determines flare start times according to the probability of event occurrence described by a Poisson distribution; and iv) injects obtained events at corresponding time stamps according to chosen cadence to the quiescent spectra, resulting in temporally and spectrally resolved datasets that represent flare activity for a chosen period. For practical application to flare frequency analysis, we chose a period of 51 days to facilitate direct comparison with the approximate duration of the observation of AU Mic in TESS sectors 1 and 27, with excluded gaps. Thus, we generated synthetic light curves and spectra compatible with observed flare frequency distributions (FFDs). The model’s output includes time-resolved spectra spanning the 10–50 000 Å range with a temporal resolution of up to 20 sec cadence, sufficient to resolve the impulsive phases of M dwarf flares, as Hawley et al. (2014) similarly demonstrated for Kepler data with a 60 sec cadence.

We proceeded to analyse the ED distributions in the synthetic stellar activity datasets generated by the YMDF module described above. We performed two complementary statistical tests to evaluate how well our synthetic population represents the observed data. The first is a Kolmogorov–Smirnov (K–S) test, which assesses the null hypothesis that both samples are drawn from the same parent distribution. The second is the Pearson R test (Pearson 1905), which measures the strength of the underlying linear correlation between the variables. Considered together, the results of these tests provide a quantitative basis for identifying the approach that offers the best agreement with the observations.

In Fig. 5, we present the results of an example of three simulation runs, each producing 1000 ED distributions, and examine the corresponding FFDs. As discussed in Paper I (i.e. Fig. 5, left panel), we employed a broken power law with the empirically derived slopes α₁ = 1.39 and α₂ = 1.8. In addition, we modelled the single power-law populations based on α₁ and α₂, respectively. We applied the K–S test to assess whether the simulated and observed distributions originate from the same population, and only included distributions with p > 0.01 in the K–S test in the figure. To assess the robustness of our approach, we repeated the test 100 times. The results across runs were consistent: typically, up to 900 populations that successfully met our condition were generated using a single power law (α₁), while approximately 200–300 populations were generated using the broken power law (α₁ and α₂), and at most 2 populations with α₂ alone pass the K–S test.

Fig. 5. Cumulative FFDs (scatter) of simulated ED distributions for broken power-law relation (blue spheres) with α1 = 1.39 and α2 = 1.8 found from the observed distribution of AU Mic white-light flares. The intrinsic observed AU Mic’s FFD is plotted as red spheres (adopted from Paper I, Fig. 5, left panel). Simulated ED distributions for single power laws with α = 1.39 and α = 1.8 are shown as light-blue and lime spheres, respectively. The resulting distributions are from one example run with 1000 random samples, plotted only if the sample meets p > 0.01 in the K–S test. The solid grey and dashed blue guides correspond to power-law coefficients α = 2.0 and α = 1.5, respectively, for a range of β coefficients.

We evaluated the linear correlation between the intrinsic stellar ED distribution and the simulated distributions by calculating the Pearson R coefficients. The results are shown in Fig. 6, indicating that both the broken power law and the single power law with α₂ provide samples that correlate with the observed intrinsic distribution better than the distributions governed by the single power law with α₁. For AU Mic, the peak of the distribution with α₁ is shifted further from unity compared to the other two distributions. We analysed all stars in the sample using the same methodology and found our conclusions to be consistent throughout, likely due to the homogeneity of the young star sample. For some stars, the peak of the Pearson R distribution with α₁ can be as low as 0.6.

Fig. 6. Pearson R coefficients indicating the linear correlation with the intrinsic observed AU Mic ED distribution, shown for simulated ED distributions generated using a broken power law with observed slopes α1 = 1.39 and α2 = 1.8 (green bars), and for single power-law distributions using α1 (blue bars) and α2 (red bars), respectively. Results are based on an example run with 1000 random samples.

When combined, these results indicate that only the broken power-law approach meets the criteria of both statistical tests and provides the closest approximation to the observed distributions. Thus, by employing the YMDF model, we reconfirm the conclusions drawn in Paper I.

4. Discussion

4.1. Comparison of stellar activity models

Loyd et al. (2018b) made an impactful contribution to the study of M dwarf activity by introducing the ‘fiducial flare’ model (Loyd 2022). Utilising FUV observations of several M dwarfs, they constructed an idealised flare template that captures both the spectral and temporal characteristics of typical stellar flares observed in these stars. This model is constructed from high-cadence FUV observations and is characterised by a box-car form, with rapid rise, long plateau, and slower decay in flux. Spectrally, the model includes both continuum and emission line enhancements, reflecting the observed increases in chromospheric and transition region lines such as C II, Si III, and Si IV during flares. By quantifying these temporal and spectral components, the ‘fiducial flare’ model provides a reproducible template for assessing the atmospheric response of exoplanets to stellar flaring, facilitating systematic studies of M dwarf activity (see e.g. Louca et al. 2023; Amaral et al. 2025). The ‘fiducial flare’ model has since become a widely used tool to investigate the effects of both rare, high-energy and frequent, low-energy flares on the chemical and physical environments of planets orbiting M dwarfs.

Following the routines in Loyd (2022), we generated a synthetic ‘fiducial flare’ stellar activity dataset for AU Mic spanning 51 days. Generating such a dataset requires access to the panchromatic spectrum of the star in its quiescent state. Although several M dwarfs from the MUSCLES survey, with their high-resolution SED, were used to refine the ‘fiducial flare’ model, these stars are generally older than 800 Myr and significantly less active than those in our sample. For AU Mic, however, a panchromatic spectrum (Feinstein et al. 2022) is available and has been used throughout this study. This enabled a direct comparison between the YMDF and the ‘fiducial flare’ models. For further context, in Fig. 3, which presents the comparison in the observed and modelled surface flux density, we plotted the ‘fiducial flare’ model as a dashed purple line for a flare with a similar equivalent duration of the flare we used for fitting (ED_TESS = 41.22 s and 43.55 s, respectively). While the YMDF and ‘fiducial flare’ models broadly reproduce the observed TESS-band emission in agreement with each other, their predictions diverge in the FUV and NUV. Here, the ‘fiducial flare’ model shows reduced energy release in these ranges compared to the YMDF model.

Previously, our range of models provided satisfactory fits to the FUV-A continuum; however, notable discrepancies remain in reproducing the continuum rise observed in the blue portion of the FUV-B segment. Thermal bremsstrahlung (TB) is a principal emission mechanism responsible for the soft X-ray emission observed in solar flares (Shibata 1996; McTiernan et al. 2019). Feinstein et al. (2022) modelled the continuum rise during AU Mic FUV flares using blackbody profiles spanning 9000–11 000 K. They found that fits converged only for electron number densities unrepresentative of typical chromospheric conditions (∼10²² cm⁻³). From COS imaging, they confirmed that count rate enhancement is restricted to the spectral trace region on the detector, supporting that it is an astrophysical feature, and they concluded that while TB alone cannot explain the FUV excess, additional emission processes are likely involved.

The TB emission during AU Mic flares predicts significant continuum enhancement at extreme-UV wavelengths, exceeding the brightness of bound-bound emission lines and recombination continua by several orders of magnitude. However, no such enhancement was observed with EUVE during AU Mic flares (Cully et al. 1993; Monsignori Fossi et al. 1996) nor in other M dwarf FUV flare observations (Froning et al. 2019; Loyd et al. 2018a). In our sample, in Karmn J07446+035, we also observe a continuum rise in the blue end of the FUV-B segment (see Fig. A.3, left panel), while the rest of the sample’s stars were not observed in this particular spectral region. Therefore, the precise physical origin of the FUV-B blue continuum rise necessitates further investigation to account for unresolved processes, possibly involving TB emission, and thus not yet incorporated into our models.

We fitted the ‘fiducial flare’ model to the observed FUV continuum. This model also does not reproduce the sharp rise in FUV-B continuum. Overall, it tends to slightly overestimate the observed flux in FUV when compared to flares of similar ED_FUV. Generally, in this model, the FUV-B segment is elevated compared to FUV-A, as illustrated in Fig. 4 (observed ED = 335.65 s; ED_fiducial = 340.99 s). For smaller flares, the ‘fiducial flare’ model matches the observed continuum at the flare peak more closely (see Fig. A.2 in Appendix A, right panel; observed ED = 156.04 s, ED_fiducial = 154.92 s). In contrast, for the largest flare in our AU Mic sample, the ‘fiducial flare’ model predicts a flux that considerably overestimates the observations (see Fig. A.2 in Appendix A, left panel; observed ED = 639.31 s, ED_fiducial = 638.01 s).

A major challenge in modelling flare spectra in the ‘fiducial flare’ model is the requirement for a panchromatic spectra across the 10–50 000 Å wavelength range. The ‘fiducial flare’ model inadequately addresses gaps in spectral data (see e.g. Fig. 3, lower right inset, 6000–7000 Å). For AU Mic, this limitation results in an underestimation of flare energy in the TESS bandpass, yielding a modelled FUV energy release smaller than expected (the figure’s right inset). Therefore, while assuming continuous panchromatic spectral coverage, the ED predicted by the ‘fiducial flare’ model would be expected to exceed the observed values. Conversely, in the FUV domain, where the panchromatic AU Mic spectra has full coverage, at similar EDs to observed, the ‘fiducial flare’ tends to overestimate both the flux density and the time-integrated flux (see Table A.2) compared to observed values, thus affirming that the model over-predicts the observed flux density at certain EDs. Since FFDs are predominantly derived from statistically robust, extensive TESS datasets, these modelling limitations would contribute to an excess of small flares and a deficit of mid-size flares, which we explore further in Sect. 4.3.

4.2. Temporal evolution of a flare in stellar activity models

The morphology of flare events exhibits significant variability, reflecting complex and poorly understood underlying physical processes (Davenport et al. 2014; Howell et al. 2014). In our sample, we observe complex flare morphologies across different spectral bands, as illustrated for 16 flares exhibiting energies in a broad range (ED_TESS∈ [11.8–437.5], ED_FUV∈ [8.69–639.31]) in the TESS and HST-COS FUV observations (Fig. A.1 in Appendix A.5). Various parametrisations have been employed to describe observed flare profiles, including composite Gaussian-exponential models (Tovar Mendoza et al. 2022). The analyses have distinguished impulsive and gradual phases in white-light flares (Kowalski et al. 2013) and identified quasi-periodic oscillations (Pugh et al. 2015). High-cadence observations have revealed flare peak rollovers, motivating continuous models without discontinuous breaks between rise and decay phases (Kowalski et al. 2016; Howard & MacGregor 2022).

We further analysed the temporal evolution of individual flares by integrating the time-resolved model spectra over TESS and FUV (excluding the detector gap) wavelength ranges into corresponding light curves of AU Mic (see Fig. 2) at EDs in TESS and FUV similar to observed. The fit for TESS in the AU Mic T2 panel reveals that the temporal evolution of the ‘fiducial flare’ lies outside the panel axis (upper row, left panel) as at the similar ED the ‘fiducial flare’ model produces fluxes ∼2.35–2.55 × 10⁵ e⁻ s⁻¹, which is much less than the observed quiescent flux of 2.79 × 10⁵ e⁻ s⁻¹. Hence, the ‘fiducial flare’ is omitted in this panel to enhance visual clarity. We attribute the lower-than-quiescent flux to a deficiency in the ‘fiducial flare’ spectrum within the 6000–7000 Å region and to the model’s approach in handling an incomplete panchromatic spectrum (previously discussed in Sect. 4.1), while applying scaling for the box-car temporal evolution. Additionally, the time-integrated flux value for the ‘fiducial flare’ is much lower than observed (see Table A.3 in Appendix A.6). However, in the AU Mic F1 and AU Mic F2 panels, we successfully plotted the ‘fiducial flare’ temporal evolutions considering EDs similar to the observed. The box-car temporal form of the ‘fiducial flare’ fails to reproduce the observed temporal evolution of these flares in detail, but also overestimates the time-integrated flux (see Table A.4 in Appendix A.5).

To our knowledge, the literature lacks extensive and comprehensive statistical studies that thoroughly characterise the detailed morphology of stellar flares to allow realistically reproducing not only frequency but also temporal profile variability in flares. Consequently, in modelling temporally resolved spectral variations of stars, it is necessary to adopt simplified approaches to maintain computational feasibility. However, the YMDF module currently allows the selection among the three aforementioned temporal flare models, and future development may include the implementation of random model selection and/or combination to enhance flexibility.

4.3. FFDs in synthetic datasets

To assess the performance of flare models in reproducing observed flare activity, in Figure 7, we present cumulative FFDs of EDs in TESS for AU Mic, alongside simulations from the YMDF and ‘fiducial flare’ models. The YMDF model (panel b) closely matches observations (panel a), while the ‘fiducial flare’ model (panel c) exhibits discrepancies, notably over-predicting small flares and under-predicting medium-sized events. The panel (d) further compares ED distributions, demonstrating that a synthetic YMDF population following a broken power law with slopes α₁ and α₂ aligns well with the data. For the largest flares, however, both models are likely to yield similar results, accurately reflecting the rarity and energetic dominance of these events. Extending the analysis to the FUV range (Fig. 8), the YMDF model exhibits an excess of small flares relative to observed FFDs, which can be attributed to the use of the broken power law in the module with α₁ and α₂ equal to those observed in TESS. Nevertheless, it produces EDs closely matching the range of the observed values. In contrast, the ‘fiducial flare’ model shows a systematic shift towards higher EDs, indicating a tendency to produce more energetic flares than observed in FUV.

Fig. 7. Cumulative FFDs of simulated and observed EDs in the TESS range. Left: Corresponding FFDs in the TESS bandpass, including the observed data (panel ‘a’.), a representative synthetic distribution generated by the YMDF model (panel ‘b’), and a representative synthetic distribution generated by the ‘fiducial flare’ model (panel ‘c’). The synthetic distributions are produced for the period of observations corresponding to the observed periods for AU Mic in TESS (∼51 days). The solid grey and dashed blue guides denote reference power-law slopes of α = 2.0 and α = 1.5, respectively. Panel ‘d’: Synthetic distributions, produced for the period of observations corresponding to the observed periods for AU Mic in TESS. Blue, red, and green bars indicate distributions produced by the ‘fiducial flare’ model, the YMDF model, and the observed ED distributions across two TESS sectors (1 and 27) for Au Mic, respectively.

Fig. 8. Cumulative FFDs in the HST-COS FUV bandpass. Left: Observed data (panel ‘a’), a representative synthetic distribution generated by the YMDF model (panel ‘b’), and a representative synthetic distribution generated by the ‘fiducial flare’ model (panel ‘c’). These synthetic distributions correspond to the AU Mic observation period in FUV (∼11.7 hours). Panel ‘d’: Synthetic distributions produced over the AU Mic FUV observation period. The bars are coloured consistently with those in the preceding figure, except for corresponding to AU Mic observation period in FUV.

In Paper I, we found that the slopes of the broken power law differ between the FUV and TESS ranges; specifically, the high-energy slope (α₂) in the FUV closely matches the low-energy slope (α₁) in the TESS range. This correspondence, also prominent in the synthetic ED distribution, suggests that a single broken power law, with α₁ ∼ 1.5 and α₂ ∼ 2.0, could adequately characterise the entire FFD, at least for young, active M stars with rotation periods shorter than 10 days. If supported by further observations, this unified description could offer a practical way to interpret flare activity across a wide range of energies in young, active stars.

These results support the proposition that young M dwarfs’ flare activity follows a unified broken power law, with FUV observations probing lower-energy flares below TESS’s detection threshold, and TESS capturing the larger, more energetic flares. Such a unified framework facilitates the interpretation of multi-wavelength flare activity and advances understanding of stellar magnetic phenomena in young active stars. Accurately capturing mid-sized flares is significant for exoplanet habitability through their influence on atmospheric chemistry and escape, and we argue that it is better achieved by the YMDF model.

4.4. Observational and modelling challenges for M dwarf flare spectra

The panchromatic spectra of M dwarfs, required by previous models such as the ‘fiducial flare’, remain challenging to obtain due to their intrinsic faintness, particularly in the UV and X-ray regimes, necessitating long exposure times with sensitive instruments. Coordinating simultaneous observations across the full spectral range is logistically complex and compounded by their high variability during flares, which renders non-simultaneous data unreliable for representing true SED (see e.g. France et al. 2016; Froning et al. 2019). However, the YMDF model has the potential to generate representative panchromatic spectra for stars lacking comprehensive observational data.

The YMDF model relies on the Kowalski et al. (2024) RHD model grid, which demonstrated the ability to reproduce the behaviour of hydrogen lines, particularly in the NUV regime as shown in previous studies. Kowalski (2022) successfully modelled the impulsive phase of AD Leo’s Great Flare using RADYN simulations of electron beam heating, supporting similar approaches for fitting flare flux to the continuum while also reproducing observed hydrogen line profiles such as Hγ. García Soto et al. (2025) conducted a forward modelling of observed Hα and Hβ spectra from a small flare on an M4.5 dwarf using a comprehensive grid search across 43 RADYN models with electron fluxes ranging from 10¹⁰ to 10¹³ erg s⁻¹ cm⁻². Their results favour the less energetic mF11-17-3 model, consistent with the flare occurring on a less active, likely older and fainter M dwarf compared to our sample stars. Similarly, Notsu et al. (2025) compared observed Hα and Hβ line profiles from AU Mic flares with the aforementioned RHD model grid, finding that electron beam heating parameters reproducing NUV/optical continuum emissions also broadly reproduce the pronounced line broadening near flare peaks. The inferred electron beam flux densities of 10¹²–10¹³ erg cm⁻² s⁻¹ exceed those typical of solar flares, aligning well with our findings.

The most prominent line feature, the Balmer series, and particularly the Balmer jump at ∼3650 Å, shows that flare flux can exceed blackbody model predictions (Kowalski et al. 2013, 2019). Using a setup similar to ours, Kowalski et al. (2025, Fig. 5) demonstrate that, despite the lack of simultaneous optical spectra to constrain the Balmer jump and optical emission lines, the observed TESS flare-only flux supports the robustness of extrapolating the RHD model into the NIR.

From solar studies, elevated continuum emission extends into the NUV (∼2000–3000 Å), indicating a common origin with the optical continuum (Joshi et al. 2021). In the NUV, emission lines (notably Mg II and Fe II) contribute more significantly, making up 20–50% of the NUV flux (Hawley et al. 2007; Jackman et al. 2024), with higher-energy flares being more continuum-dominated. For example, Kowalski et al. (2019) measured line contributions of ∼40% in the 2510–2841 Å range for two GJ 1243 flares observed with HST-COS. However, our model setups, COScut1 and COScut2, incorporate all major lines as background opacities. Comparing these model setups with the H6 setup, both in light curves (Fig. 2) and spectra (Figs. 4, A.2, A.3), reveals comparable time-integrated fluxes, implying similarity in ED_FUV, and the fitting of these models to the continuum fluxes is yielding similar results. This suggests that the simplified H6 approach is at least capable of adequately capturing the flare energy budget, while addressing the challenge of modelling variable spectra in young M dwarfs.

The YMDF model in the H6 setup primarily captures the continuum enhancement observed during stellar flares, along with the prominent hydrogen lines evolution as shown in previous work. Multiple studies have demonstrated that simple models that focus on the continuum rise and hydrogen recombination can provide a robust match to observed flare spectra, even when more detailed line emission modelling is not feasible. For example, Hawley & Pettersen (1991) found that roughly one-third of the line emission occurs during the impulsive phase, whereas the majority of the continuum emission is released during this same period. Because the continuum dominates the overall energy output, most of the total flare energy is also emitted during the impulsive phase, at least in the 1200–8000 Å wavelength range. It is worth noting that radiation emitted within this wavelength range is particularly relevant for studies of exoplanet atmospheres. Photons at these wavelengths drive the majority of photodissociation reactions (Ranjan et al. 2020), which can ultimately lead to atmospheric depletion.

Studies of M-star flares indicate that FUV emission can precede white-light emission, implying that it may originate from the initial heating, compression, and evaporation of plasma at the flare footpoints (Froning et al. 2019; Hawley et al. 2003; MacGregor et al. 2021). Furthermore, flare observations of M dwarfs have reported temperatures reaching up to 20 000 K and in some cases as high as 40 000 K in both the optical and FUV regimes (see e.g. Howard et al. 2019; Froning et al. 2019). Although our model configuration is capable of reproducing the high energies released in both the FUV and NUV ranges, it is currently unable to replicate the observed differences in the onset times of white-light flares and FUV emission. This discrepancy highlights a limitation in the present modelling approach. Further investigation is required to identify and accurately parametrise the physical processes responsible for these timing differences.

4.5. Summary and conclusions

In this study, we tested several spectral setups of the combinations of the two RHD models with high- and low-energy electron beams. We found that the straightforward H6 non-LTE hydrogen atom model with the inclusion of several important atoms in LTE can satisfactorily reproduce the observed continuum rise in the TESS and partly in the FUV range for flares spanning a wide range of EDs. This approach models stellar flare continua for these two regions, and in NUV and optical according to previous literature (see e.g. Kowalski et al. 2025; Notsu et al. 2025; García Soto et al. 2025). Improvements should include expanding the atomic and molecular set, and enhancing the treatment of atmospheric structure and time-dependent effects. Comparison with high-cadence, multi-wavelength data will further constrain models and improve their realism.

The YMDF model demonstrated reasonable performance in the probed wavelength regimes while producing temporally resolved spectra, with the exception of the blue continuum rise in the FUV-B segment. The models calculated with the H6 and H6CII spectral setups reproduced the observed overall energy release in flares, in both the optical and FUV components, separately. While temporal flare models cannot fully capture the complex morphology of flares, they are still effective in reproducing the energy release across a range of flare events, particularly the Mendoza model. Our results indicate that FFDs of young M dwarfs follow a common statistical law, well-described by a two-coefficient broken power law. The YMDF module supports synthetic population generation based on these findings, with flexibility to modify the setup as needed.

Future work will apply YMDF to simulate the effects of flares, especially mid-size, on planetary atmospheres through detailed chemical kinetics modelling under high-energy irradiation. This will clarify the impacts on atmospheric composition and loss processes. Additionally, the YMDF framework will be potentially extended to model stellar activity beyond young early-type M dwarfs.

Data availability

The data used for this study are publicly available in the Mikulski Archive for Space Telescopes (MAST³). The code, produced by this study, is made publicly accessible in GitHub (https://github.com/cepylka/YMDF).

Acknowledgments

We acknowledge financial support from the Research Council of Norway (RCN), through its Centres of Excellence funding scheme, projects number 332523 (PHAB, Centre for Planetary Habitability) and number 262622 (RoCS, Rosseland Centre for Solar Physics). We thank James Davenport (University of Washington, USA) and Adalyn Gibson (University of Colorado Boulder, USA) for the fruitful discussion. We thank an anonymous referee for helpful comments and critiques, which have improved the quality of this manuscript.

References

Appendix A: Detailed methodology

A.1. Radiative transfer code RH1.5D

RH1.5D allows us to calculate synthetic disk-centre intensity I_ν, but in order to calculate monochromatic flux according to

$$\begin{array}{c}\hfill F_{\nu }(z)=2\pi {\displaystyle {\int }_0^1}{I}_{\nu }\mu d\mu ,\end{array}$$(A.1)

where μ = cos θ, one should calculate I_ν(μ). The previous version of the RH radiative transfer code allowed calculating that for μ with weights, so flux can be calculated as follows:

$$\begin{array}{c}\hfill F_{\nu }=2\pi {\displaystyle {\int }_{-1}^1}{I}_{\nu }\mu d\mu \mathrm{a}\mathrm{s}{F}_{\nu }=2\pi \frac{b-a}{2}({\omega }_1{I}_{\nu }(1){\mu }_1+{\omega }_2{I}_{\nu }(2){\mu }_2+...)\end{array}$$(A.2)

We calculate I_ν(μ) in order to obtain the μ-dependent intensity using the following approach. Specifically, we employed parameters that can be derived with RH1.5D, such as source function S_ν and total opacity χ:

$$\begin{array}{c}\hfill I_0(\mu )=\frac{1}{\mu }{\displaystyle {\int }_0^{\infty }}(S_{\nu }({\tau }_{\nu }^{\prime }){exp}^{-{\tau }_{\nu }^{\prime }/\mu }d{\tau }_{\nu }^{\prime }.\end{array}$$(A.3)

A.2. TESS data and fitting

The TESS electron rate F_e⁻/s converts to physical flux through a two-step photometric calibration process (Sullivan et al. 2015). The instrumental magnitude is derived from electron counts:

$$\begin{array}{c}\hfill T_{\mathrm{mag}}=-2.5{log}_{10}(F_{e^{-}/\mathrm{s}})+20.44\end{array}$$(A.4)

The TESS photometric system is anchored to the Vega-relative Cousins I band, with its zero-point flux density defined as

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{Jy}}=2416\times {10}^{-0.4T_{\mathrm{mag}}}\end{array}$$(A.5)

The 2416 Jy value originates from the Vega flux density at the central wavelength of TESS’s bandpass, λ_eff = 7865 Å, calibrated through synthetic photometry of Vega’s spectrum (Sullivan et al. 2015). To approximate the bandpass-integrated flux F_{erg, s⁻¹, cm⁻²}, the effective wavelength serves as a spectral anchor for unit conversion:

$$\begin{array}{c}\hfill F_{\mathrm{erg},{\mathrm{s}}^{-1},{\mathrm{cm}}^{-2}}\approx F_{\nu }^{\mathrm{Jy}}\times {10}^{-23}\times \frac{c}{{\lambda }_{\mathrm{eff}}^2}\times \mathrm{\Delta }{\lambda }_{\mathrm{TESS}}\end{array}$$(A.6)

where c is the speed of light and Δλ_TESS ≈ 4000 Å represents the approximate bandpass width. This approximation assumes wavelength-independent flux density across the bandpass. The final quantity represents the total flux integrated over TESS’s 6000–10000 Å sensitivity range.

We convert all model fluxes (in cgs units) using a simple approximation for the T0 Vega following (Sullivan et al. 2015). It is convenient to define a star’s TESS magnitude, assuming Vega has a flux density of F_λ = 4.03×10⁻⁶, erg s⁻¹ cm⁻² Å⁻¹ at T = 0. Then, using Eq. A.4, we convert the magnitudes to electron counts per second.

For several flares in our TESS sample, we additionally provide the observational data and state that in Table A.1.

A.3. YMDF inverse calculations of the model setup coefficients

For producing a synthetic temporally and spectrally resolved flare dataset, we determine the coefficients of the model components, m2F12-37-2.5 and cF13-500-3, by utilising the observed EDs and the ratios derived from our model setups. This is accomplished by integrating the contributions of each component over a broad wavelength range F₁ and F₂, as well as by integrating the pre-flare flux F₀ within the model framework, and substituting them into the ED equation (see Eq. 2):

$$\begin{array}{c}\hfill ED={\displaystyle \int (\frac{F_{\mathrm{model}}(t)·[cx(F_1-F_0)+x(F_2-F_0)]+F_0}{F_0}-1)}dt,\end{array}$$(A.7)

where c is the ratio specific to the model setup, and F_model(t) is a temporal flare model. The full derivation is presented below:

$$\begin{array}{cc}\hfill \mathit{ED}& ={\displaystyle \int \left(\frac{F_{\mathrm{model}}(t)·[cx(F_1-F_0)+x(F_2-F_0)]}{F_0}\right)}dt\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill \mathit{ED}& =x{\displaystyle \int \left(\frac{F_{\mathrm{model}}(t)·[c(F_1-F_0)+(F_2-F_0)]}{F_0}\right)}dt\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill A& =\frac{c(F_1-F_0)+(F_2-F_0)}{F_0},S={\displaystyle \int F_{\mathrm{model}}}(t)dt\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill \mathit{ED}& =x·A·S,x=\frac{\mathit{ED}}{A·S}\hfill \end{array}$$(A.11)

This approach enables inverse calculations of the coefficient for the second model component (x), the first model component is calculated by multiplying it by the ratio (c ⋅ x), and is used in the YMDF module to produce the simulated ED distribution.

A.4. COS data handling

We worked with the archival data retrieved by Paper I from the HST-COS G130M grating, programme IDs: 11533 (PI Green, 2MASS J18141047-3247344, observation ID "lb3q01060"); 14784 (PI Shkolnik, 2MASS J01521830-5950168, observation ID ldab08030; 2MASS J02001277-0840516, observation ID "ldab54030"; 2MASS J03315564-4359135, observation ID "ldab10010"; 2MASS J22025453-6440441, observation ID "ldab06030"; 2MASS J02365171-5203036, observation ID "ldab02010"); 15955 (PI Richey-Yowell, HIP-107345, observation ID "le5104030"; HIP 1993, observation ID "le5108030"); 16164 (PI Cauley, AU Mic, observation IDs "lebb01010", "lebb02010", "lebb03010"); 16482 (PI Roman-Duval, 2MASS J11173700-7704381, observation ID "lein2d010"); 17428 (PI France, AU Mic, observation IDs "lf71z1010", "lf7101010"; Karmn J07446+035, observation IDs "lf7109010", "lf71z9010"). For several flares in our FUV sample, we additionally provide the observational data and state that in Table A.2.

For AU Mic and Karmn J07446+035, the configuration of G130M grating provided coverage from approximately λ∈ [1060,1360] Å with a detector gap λ∈ [1210,1225] Å, masking the bright Lyα emission feature to avoid detector saturation. For the rest of the stars in the sample, the wavelength coverage of this configuration is λ∈ [1170,1430] Å and includes strong emission lines of C II, C III, Si III, Si IV, and N V formed in the stellar transition region, along with several weaker lines, including some coronal iron lines. Lyα and O I also fall within the bandpass; however, they are generally obscured due to contamination from geocoronal airglow.

All data were processed with the CALCOS pipeline, and time-resolved analysis was carried out using the corrtag products with a 10 s sampling interval. During the CALCOS calibration, regions affected by bright airglow, specifically around Lyα and O I, are automatically flagged and masked. We applied the same masking when integrating fluxes for light curve fitting in both the observed data and the fitted configurations (COScut1 and COScut2), ensuring consistency between the calibrated data and our analysis. The same masking procedure was also adopted when constructing light curves for the inverse models (H6 and H6CII), and in constructing the ‘fiducial flare’ light curves.

In Figs. 4, A.2, and A.3, the model spectra are shown with masking applied for the COScut1/COScut2 setups, but without masking in the inverse models. For direct model-observation comparisons, only continuum regions were used in the corresponding spectral images. In the integrated flux analysis (Fig. 2), identical masks were applied to both the observed and model data in each panel to ensure consistency.

A.5. Electron beam model components

We reported the ratios between the coefficients in the H6 and H6CIII setups derived from the TESS and FUV bandpasses (as shown in Fig. 1 and listed in Table 2). The ratios obtained from the HST-COS spectra (COScut1 and COScut2 setups) are systematically larger. We suggested that this issue can be resolved by including a more extensive list of spectral lines in the RH1.5 calculations. Additionally, this indicates differing behaviour among the model components. For broadband energy estimates based on model spectra, the coefficient ratio corresponding to the TESS bandpass provides an accurate description of the flare energetics across a wide wavelength range. In contrast, the HST-COS-based ratios tend to overestimate the relative contribution of the m2F12-37-2.5 component, likely due to the narrower spectral coverage and the dominance of this energetic beam model in the FUV regime. A plausible interpretation is that the HST-COS G130M grating samples the spectral region where the m2F12-37-2.5 model predicts the strongest flux, reaching values of 10⁶-10⁹ erg cm⁻² s⁻¹ under the idealised assumption of the entire stellar surface participating in the flare (see Fig. 3). The comparison further suggests a physical distinction between the spectral domains. While in the TESS bandpass, reproducing the observed continuum rise requires the inclusion of a persistent, high-energy electron beam (represented by the cF13-500-3 component), and in the FUV, the emission is better explained by the impulsive, lower-flux beam dominating m2F12-37-2.5 model, the optical and FUV ranges better represent different processes in flares, namely ribbons and kernels formation on the stellar surface, as we mentioned in Sect. 2.2.

A.6. Temporal flare models: further considerations

Fig. A.1. Morphology of several stellar flares observed in TESS (upper row) and HST-COS FUV (lower row). Upper row: Light curves obtained with TESS. The details of the observations for these flares can be found in Table A.1. In the panels, the labels 2MASS...036 and 2MASS...441 stand for 2MASS J02365171-5203036 and 2MASS J22025453-6440441, respectively. Bottom row: Light curves produced by integrating observed flux over the wavelength range in FUV with HST-COS. The details of the observations for these flares can be found in Table A.2. In the panels, the labels 2MASS...168, Karmn...+035, 2MASS...036, and 2MASS...441 stand for 2MASS J01521830-5950168, Karmn J07446+035, 2MASS J02365171-5203036, and 2MASS J22025453-6440441, respectively.

As previously discussed in Section 2.4, the Mendoza model provides a better match to the total time-integrated flux of the flare light curve, accurately reproducing the overall energy output. In practical terms, this means that when integrating the flare emission over time, the total energetic budget is very well recovered. However, it typically underestimates the flare amplitude (see Fig. 2), which indicates that the model has difficulty reproducing the sharp, high-intensity peak of some observed events. Therefore, as we adopted the simpler Gaussian rise–exponential decay temporal flare model proposed by Feinstein et al. (2022), we suggest that with a modification allowing the rise phase to be described by a Student’s t-distribution having a non-integer fractional degrees of freedom of 2.05, this model would successfully reproduce the flare amplitude. This modification accounts for sharper peaking and better captures the observed flare rise morphology, particularly for flares with fast impulsive phases. Additionally, we set the decay parameter to γ = −0.09, which corresponds to a slower decay in time and produces a better match to the gradual phase in some flares.

However, the complex morphology of observed flares suggests that this component should be modified to reflect the evolution of individual flares. This modification is not yet implemented in the YMDF module, partly because our approach cannot currently be grounded in statistical studies of flare temporal profiles, which are absent and remain a subject for future work. Of all the approaches evaluated, this parametrisation yielded the best agreement with the observed flare profiles. Although the fitted coefficients are expected to be larger when using the Feinstein model compared to the Mendoza and Davenport models, the resulting ratios of the obtained coefficients remain very similar between the two models (for 2MASS J02365171−5203036 (COScut2), ${\widehat{X}}_1/{\widehat{X}}_2\sim$ 11.30 and for AU Mic (COScut1), the corresponding ratios ∼ 13.54, respectively).

For several stars, namely 2MASS J02365171-5203036, 2MASS J18141047-3247344, and Karmn J07446+035, the Feinstein model enabled us to successfully fit the YMDF H6 model to the peak of the flare observed in the FUV HST-COS spectra (see e.g. Fig. 2 for AU Mic F1 and Karmn J07446+035 F1 panels in the bottom row), while for TESS fits it can overestimate the peak (see the same figure, AU Mic T2 and HIP 107345 T1 panels in the upper row). Although the modified Feinstein model yields a superior fit to the flare peak and enables a more accurate estimation of the energy at the most energetic phase, including for relatively small flares observed in the FUV with HST-COS comparing to the TESS observations, it does not capture the full energy release within a specific wavelength range, for which the Mendoza model performs significantly better (see Tables A.3 and A.4).

Fig. A.2. AU Mic spectra at the peak of flares observed with HST-COS. In both panels, green spheres denote continuum flux measurements at the flare peaks, while grey diamonds indicate the quiescent continuum levels. Associated uncertainties are depicted as light-coloured lines – green for the flare peak and grey for quiescence – with only upper error bars shown for clarity, given the symmetry of the errors. The light-red lines represent the best-fit obtained using the COScut1 setup. The dash-dotted green and light-green lines correspond to the H6CIII models with ${\widehat{X}}_1$/${\widehat{X}}_2$ = 8.74 derived from TESS photometry for the flare peak and pre-flare states, respectively, while the dashed blue lines show the simplified H6 models with ${\widehat{X}}_1$/${\widehat{X}}_2$ = 5.03. The background light-grey spectra are the xdsum products for each respective orbit, representing the summed, extracted 1D spectra from all exposures within an orbit and providing high signal-to-noise reference spectra for comparison. Light-orange lines display smoothed panchromatic AU Mic spectra from Feinstein et al. (2022), and dashed purple lines depict the ‘fiducial flare’ model from Loyd et al. (2018b) at similar equivalent durations, included for context.

Fig. A.3. Spectra at the peak of flares observed with HST-COS for Karmn J07446+035 (left) and 2MASS J02365171-5203036 (right). For both stars, green spheres indicate continuum flux measurements at the flare peak, while grey diamonds represent the quiescent continuum levels. Associated uncertainties are shown as light-coloured lines – green for flare peak and grey for quiescence – with only upper error bars displayed for clarity, as the errors are symmetric. The light-red lines correspond to the best-fit continuum models obtained using the COScut1 (left panel) and COScut2 (right panel) setups. The dash-dotted green and light-green lines represent the inverse H6CIII models for the flare peak and pre-flare states, respectively, while the dashed blue lines show the simplified inverse H6 models. The background light-grey spectra are the xdsum spectrum.

All Tables

All Figures

Fig. 1. Coefficients of the fitted models for flares in young M dwarfs. The values of ${\widehat{X}}_1$ and ${\widehat{X}}_2$ are plotted on the y- and x-axes, respectively. Each point represents an individual flare event and is colour-coded by the flare energy released in the TESS (left) or the FUV bandpass (right). Model fitting was performed using multiple setups: H6 and H6CIII for the TESS data, and COScut1 and COScut2 for the HST-COS data. Linear regressions fitted to the flare events for each star are shown as dashed coloured lines. To guide the eye, in the panels, the H6, H6CIII, COScut1 and COScut2 ratios are shown as dotted red, dotted black, dash-dotted grey, and dotted grey lines, respectively.

In the text

Fig. 2. Light curves of several flares observed in TESS and HST-COS. Upper row: TESS flux light curves (e− s−1) for flares denoted as AU Mic T2, HIP 107345 T1, and 2MASS J02365171-5203036 T1, as dashed lines with red spheres for 20-second observation stamps and black error bars with capped ends. Observation details are provided in Table A.1 in Appendix A.2. Note that errors can be smaller than the sphere size. Three temporal flare model MCMC fitting results are shown: Mendoza (light blue), Davenport (light coral), and Feinstein (teal) models fitted to the observed data. The H6 and H6CIII models, which were inversely calculated using observed EDs for these particular flares, are shown as solid blue and green lines with data-point spheres. In the upper row, the H6 model is hidden by the H6CIII due to overlap as their flux values are very similar. Bottom row: FUV HST-COS flux observations for flares denoted as AU Mic F1, AU Mic F2, and Karmn J07446+035 F1, plotted with consistent styles and model fits, as the upper row. The details of observations for these flares can be found in Table A.2 in Appendix A.4.

In the text

Fig. 3. Observed and modelled surface flux density. The m2F12-37-2.5 (dark green) and cF13-500-3 (light green) stellar atmosphere models are plotted at coefficients ${\widehat{X}}_1=1$ and ${\widehat{X}}_2=1$ (the theoretical assumption of the entire stellar surface exhibiting flaring activity). For AU Mic during a flare, the solid red and dashed orange lines show the fitted model spectrum at flare maximum using coefficients from the example fit (the flare data are stated in Table A.1 in Appendix A.2). The inverse H6 model at flare maximum, pre-flare, and the ‘fiducial flare’ model from Loyd et al. (2018b) with a comparable ED is plotted as dash-dotted blue, dash-dotted light-blue, and dashed purple lines, respectively. Blackbody curves for 3850 K, 6000 K, 8000 K, and 10 000 K are included in solid, dotted, dashed, and dash-dotted black lines. The panchromatic AU Mic flux from Feinstein et al. (2022) is plotted in light grey. Insets highlight the FUV (lower left) and TESS (lower right) wavelength ranges.

In the text

Fig. 4. Flare spectra in FUV HST-COS. Green circles and grey diamonds indicate the continuum flux at the flare peak and pre-flare, respectively, with uncertainties displayed as corresponding light-coloured lines. Only the upper error bars are shown for visibility, as the errors are symmetric. The best-fit model to the continuum fluxes, derived using the COScut1 spectral setup, is plotted as a light-red line. The dashed blue, dash-dotted green, and dash-dotted light-green lines correspond to the H6 and H6CIII setups at the peak, and pre-flare H6, respectively. The light-grey, light-orange, and dashed purple lines represent the sum of the extracted 1D spectra from all exposures in the orbit, the smoothed panchromatic spectrum of AU Mic from Feinstein et al. (2022), and the ‘fiducial flare’ model from Loyd et al. (2018b) at a similar ED, respectively.

In the text

Fig. 5. Cumulative FFDs (scatter) of simulated ED distributions for broken power-law relation (blue spheres) with α1 = 1.39 and α2 = 1.8 found from the observed distribution of AU Mic white-light flares. The intrinsic observed AU Mic’s FFD is plotted as red spheres (adopted from Paper I, Fig. 5, left panel). Simulated ED distributions for single power laws with α = 1.39 and α = 1.8 are shown as light-blue and lime spheres, respectively. The resulting distributions are from one example run with 1000 random samples, plotted only if the sample meets p > 0.01 in the K–S test. The solid grey and dashed blue guides correspond to power-law coefficients α = 2.0 and α = 1.5, respectively, for a range of β coefficients.

In the text

Fig. 6. Pearson R coefficients indicating the linear correlation with the intrinsic observed AU Mic ED distribution, shown for simulated ED distributions generated using a broken power law with observed slopes α1 = 1.39 and α2 = 1.8 (green bars), and for single power-law distributions using α1 (blue bars) and α2 (red bars), respectively. Results are based on an example run with 1000 random samples.

In the text

Fig. 7. Cumulative FFDs of simulated and observed EDs in the TESS range. Left: Corresponding FFDs in the TESS bandpass, including the observed data (panel ‘a’.), a representative synthetic distribution generated by the YMDF model (panel ‘b’), and a representative synthetic distribution generated by the ‘fiducial flare’ model (panel ‘c’). The synthetic distributions are produced for the period of observations corresponding to the observed periods for AU Mic in TESS (∼51 days). The solid grey and dashed blue guides denote reference power-law slopes of α = 2.0 and α = 1.5, respectively. Panel ‘d’: Synthetic distributions, produced for the period of observations corresponding to the observed periods for AU Mic in TESS. Blue, red, and green bars indicate distributions produced by the ‘fiducial flare’ model, the YMDF model, and the observed ED distributions across two TESS sectors (1 and 27) for Au Mic, respectively.

In the text

Fig. 8. Cumulative FFDs in the HST-COS FUV bandpass. Left: Observed data (panel ‘a’), a representative synthetic distribution generated by the YMDF model (panel ‘b’), and a representative synthetic distribution generated by the ‘fiducial flare’ model (panel ‘c’). These synthetic distributions correspond to the AU Mic observation period in FUV (∼11.7 hours). Panel ‘d’: Synthetic distributions produced over the AU Mic FUV observation period. The bars are coloured consistently with those in the preceding figure, except for corresponding to AU Mic observation period in FUV.

In the text

Fig. A.1. Morphology of several stellar flares observed in TESS (upper row) and HST-COS FUV (lower row). Upper row: Light curves obtained with TESS. The details of the observations for these flares can be found in Table A.1. In the panels, the labels 2MASS...036 and 2MASS...441 stand for 2MASS J02365171-5203036 and 2MASS J22025453-6440441, respectively. Bottom row: Light curves produced by integrating observed flux over the wavelength range in FUV with HST-COS. The details of the observations for these flares can be found in Table A.2. In the panels, the labels 2MASS...168, Karmn...+035, 2MASS...036, and 2MASS...441 stand for 2MASS J01521830-5950168, Karmn J07446+035, 2MASS J02365171-5203036, and 2MASS J22025453-6440441, respectively.

In the text

Fig. A.2. AU Mic spectra at the peak of flares observed with HST-COS. In both panels, green spheres denote continuum flux measurements at the flare peaks, while grey diamonds indicate the quiescent continuum levels. Associated uncertainties are depicted as light-coloured lines – green for the flare peak and grey for quiescence – with only upper error bars shown for clarity, given the symmetry of the errors. The light-red lines represent the best-fit obtained using the COScut1 setup. The dash-dotted green and light-green lines correspond to the H6CIII models with ${\widehat{X}}_1$/${\widehat{X}}_2$ = 8.74 derived from TESS photometry for the flare peak and pre-flare states, respectively, while the dashed blue lines show the simplified H6 models with ${\widehat{X}}_1$/${\widehat{X}}_2$ = 5.03. The background light-grey spectra are the xdsum products for each respective orbit, representing the summed, extracted 1D spectra from all exposures within an orbit and providing high signal-to-noise reference spectra for comparison. Light-orange lines display smoothed panchromatic AU Mic spectra from Feinstein et al. (2022), and dashed purple lines depict the ‘fiducial flare’ model from Loyd et al. (2018b) at similar equivalent durations, included for context.

In the text

Fig. A.3. Spectra at the peak of flares observed with HST-COS for Karmn J07446+035 (left) and 2MASS J02365171-5203036 (right). For both stars, green spheres indicate continuum flux measurements at the flare peak, while grey diamonds represent the quiescent continuum levels. Associated uncertainties are shown as light-coloured lines – green for flare peak and grey for quiescence – with only upper error bars displayed for clarity, as the errors are symmetric. The light-red lines correspond to the best-fit continuum models obtained using the COScut1 (left panel) and COScut2 (right panel) setups. The dash-dotted green and light-green lines represent the inverse H6CIII models for the flare peak and pre-flare states, respectively, while the dashed blue lines show the simplified inverse H6 models. The background light-grey spectra are the xdsum spectrum.

In the text