© A. J. Maas et al. 2022

1 Introduction

M-dwarfs are prime targets for the search for habitable conditions, that is, rocky planets in the habitable zone of the stars. The small stellar radius and low effective temperatures of M-dwarfs permit not only the detection of habitable planets but also the atmospheric characterization of those in close-in orbits (e.g. Kaltenegger & Traub 2009, and references therein). M-dwarfs represent roughly 73% of all stars (Dole 1964; Henry et al. 2006, Henry et al. 2018). Late M-dwarfs, including ultra-cool M-dwarfs, have masses M_* ≤ 0.2 M_⊙ and effective temperatures T_eff ≤ 3000 K (Gillon et al. 2017). Ultra-cool dwarfs are stellar and sub-stellar objects with spectral type later than M7V (Kirkpatrick et al. 1995, 1997). Mid and early M-dwarfs are more massive, 0.2 M_⊙ < M_* ≤ 0.4 M_⊙, 0.4 M_⊙ < M_* ≤ 0.6 M_⊙ respectively, and thus have surface temperatures of up to 3900 K (Heath et al. 1999; Farihi et al. 2006). Ground-based targeted searches such as MEarth (Irwin et al. 2009), TRAPPIST (Gillon et al. 2013), ExTrA (Bonfils et al. 2015) and SPECULOOS (Sebastian et al. 2020), as well as the space missions Kepler (Koch et al. 2010), K2 (Howell et al. 2014), and TESS (Ricker et al. 2014), have boosted the discoveries of exoplanets around nearby M-dwarfs. Consequently, there is growing interest in putting constraints on the flare rates of the host stars and the respective flare energies, temperatures, and areas (e.g Fuhrmeister et al. 2008; Schmidt et al. 2016; Davenport 2016; Vida et al. 2017; Günther et al. 2020; Howard et al. 2020; Johnson et al. 2021) in order to understand how these can affect planetary atmospheres and the prospects of searching for habitable conditions.

Stellar flares are explosive magnetic reconnection events (Pettersen 1989) occurring in main sequence (MS) stars with convective envelopes. When opposing magnetic field lines approach each other, they are prone to reconnect and release the energy that was stored in the magnetic field in the form of accelerated particles and electromagnetic radiation. These stochastic events have a duration ranging from minutes to several hours (e.g. Benz 2017, and references therein). The flare occurrence frequency and the released energy are both dependent on the magnetic field of the stellar surface, which is known to decrease with stellar age for MS stars (Reiners et al. 2022). Due to the loss of angular momentum, the stellar dynamo quiets during the star’s lifetime (Skumanich 1972). For low-mass stars, this process takes longer than for solar type stars and thus, in comparison, low-mass stars experience a prolonged active time (West et al. 2008), they produce high UV and X-ray flux over longer periods (Audard et al. 2000; Hawley et al. 2014), which can alter the atmospheres of close-by planets and might be a serious threat to surface life (Segura et al. 2010; Tilley et al. 2019). Other studies suggest that additional UV flux from stellar flares might trigger and drive prebiotic chemistry on planetary surfaces (Ranjan et al. 2017; Rimmer et al. 2018). Several recent studies of stellar flares on nearby M-dwarfs with planetary companions address surface habitability (Ribas et al. 2016; O’Malley-James & Kaltenegger 2017; Meadows et al. 2018; Vida et al. 2019; Glazier et al. 2020) and use constant flare temperatures of T_eff = 9000−10000 K. However, a lower or higher flare temperature would impact the UV flux models for the star and therefore the flux incident on the planetary companions.

The primary output of flares are accelerated particles, which then precipitate back on the chromosphere of the star, emitting nonthermal radiation and subsequently heat the chromospheric plasma (Benz & Güdel 2010; Aschwanden et al. 2017; Benz 2017). The amount of heating differs locally and is strongest at the footpoints of the magnetic loop after reconnection. This local heating gives rise to thermal continuum emission and line emission in the optical, UV, and soft X-ray regime (Benz & Güdel 2010; Kowalski et al. 2013; Aschwanden et al. 2017). One example is the Hα line. At the footpoints of the magnetic fields, accelerated electrons and protons will recombine producing hydrogen emission lines, the so-called Ha ribbons (Kowalski et al. 2013). The main contribution to the total energy budget of a flare comes from the white-light continuum emission (Benz & Güdel 2010; Aschwanden et al. 2017). The total energy budget is frequently estimated by a black body of effective temperature T_eff = 9000−10000 K (Hawley & Pettersen 1991; Kretzschmar 2011; Kowalski et al. 2013; Shibayama et al. 2013; Davenport 2016; Paudel et al. 2018; Howard et al. 2018; Jackman et al. 2018, 2019; Günther et al. 2020). While this knowledge relies mainly on magnetically active M-dwarfs (dMe stars) with spectral types dM3e-dM4.5e from spectroscopic and photometric observations (e.g. Kowalski et al. 2013, and references therein), it is usually assumed that flares on ultra-cool M-dwarfs like TRAPPIST-1 with spectral type M7.5V behave similarly (Paudel et al. 2018; Günther et al. 2020; Glazier et al. 2020).

Multi-passband photometry is a useful tool for determining continuum black body temperatures of stellar flares (Hawley et al. 2003; Kowalski et al. 2013; Howard et al. 2020). The black body temperature of stellar flares is not constant and should be understood as an energy budget of the continuum emission. As flares are subject to different physical processes in different atmospheric layers, it is nontrivial to attribute a single temperature characterizing the flare. If flares are observed in several photometric filters, recent methods can be employed to retrieve black body temperatures (Howard et al. 2020). For typical active M-dwarf stars spectroscopic and photometric observations are consistent with a black body with temperatures of 9000−14 000 K at the peak and 5500−7000 K in the gradual phase (Kowalski et al. 2013). Howard et al. (2020) investigated a large sample of stars with photometric flares and spectral types of M0-M7V. These authors found the total temperature averaged over their sample to be $T_{\text{eff}}={11000}_{-2600}^{+3500}$ K and at peak $T_{\text{eff}}={14000}_{-3400}^{+8300}$ K.

For studying flares on low-mass stars in the context of exoplanet research, the most promising target is TRAPPIST-1. This target is a well-studied late M-dwarf with spectral type M7.5V (Gizis et al. 2000) and is known to host seven rocky exoplanets of which three lie within the habitable zone, which is the zone where the surface flux is just right for water to be liquid (Gillon et al. 2016, 2017; Luger et al. 2017). Because TRAPPIST-1 flares frequently, with a flare with an energy above 10²⁹ erg occurring every 1−2 days (Paudel et al. 2018), the impact of its flaring activity on the atmospheres of its orbiting exoplanets is being investigated (Vida et al. 2017; O’Malley-James & Kaltenegger 2017; Paudel et al. 2018; Glazier et al. 2020; Estrela et al. 2020). Paudel et al. (2018) analysed K2 data for TRAPPIST-1 and detected 39 flares, and a superflare rate of 4.4+2.0//−0.2 flares per year. A flare is considered a superflare if it reaches bolometric energies of 10³³−10³⁸ erg (Shibayama et al. 2013). Glazier et al. (2020) did not detect flares in a 2-yr survey with Everyscope, an array of small telescopes imaging the entire accessible sky all at once (Law et al. 2016; Ratzloff et al. 2019). Glazier et al. (2020) and Paudel et al. (2018) argue that the possible atmospheres of the planets in the habitable zone of TRAPPIST-1 are, with the current flare rate, not in danger of complete ozone depletion, but also suggest that the UV flux from TRAPPIST-1 is not enough to trigger and sustain abiogenesis. Tilley et al. (2019) showed that UV radiation and proton fluxes from frequent flaring can deplete the ozone layer of an Earth-like planet almost entirely over timescales of several million years. Estrela et al. (2020) examined the impact of the UV radiation from flares on the potentially habitable planets of TRAPPIST-1. These authors argue that UV-sensitive organisms could survive if the planet has an ozone layer or if their natural living environment is at least 8 m below surface ocean level. These studies highlight the importance of considering flares when assessing the habitability of a planet. Despite the growing interest in flaring processes of ultra-cool M-dwarfs and their subsequent UV radiation, the flare temperature of the continuum emission remain uncertain for late M-dwarfs. For a L-dwarf (L1V) flare, Gizis et al. (2013) found a black body spectrum consistent with an effective temperature of T_eff = 8000 ± 2000 K.

Here, we present an automated search for stellar flares on TRAPPIST-1 using MuSCAT1 (Multicolor simultaneous camera for studying atmospheres of transiting exoplanets; Narita et al. 2015) and MuSCAT2 (Narita et al. 2019) data and the open source package AltaiPony (Ilin 2021). For the first time using the simultaneous multi-color photometry of the MuSCAT instruments, we estimated effective black body temperatures of TRAPPIST-1 flares.

The paper is structured as follows: in Sect. 2, we present the instruments, observations and preparation of the photometric data. Section 3.1 describes the flare detection, and in Sect. 3, we describe the methods used to derive black body temperatures. The results are presented in Sect. 4 and discussed in Sect. 5. Finally, we summarize our findings in Sect. 6.

2 Observations

MuSCAT1 (M1) and MuSCAT2 (M2) are part of the Global Multi-Color Photometric Monitoring Network for Exoplanetary Transits. M1 is mounted on the 1.88m telescope at the Okayama Astro-Complex in Japan and M2 is located at the Teide Observatory in Tenerife, Spain and mounted on the 1.52 m Telescopio Carlos Sánchez (TCS). While M1 permits the use of three photometric filters g, r, and z_s at the same time, M2 observes in four band-passes simultaneously, namely g (400−550 nm), r (550−700 nm), i (700−820 nm), and z_s (820−920 nm) (Narita et al. 2015, 2019). M1 and M2 are equipped with 1024 ×1024 pixel CCDs with a pixel scale of 0.36 and 0.44 arcsec pixel⁻¹, resulting in a field of view of 6.1 × 6.1 arcmin and 7.4 × 7.4 arcmin, respectively. Both instruments are used for follow-up observations of transiting exoplanets by for example TESS for confirmation and/or validation, and studies of exoplanet atmospheres. The data are reduced and optimized with the M2 transit pipeline (Parviainen et al. 2020). The pipeline covers the reduction of generic (nontransit) photometry, transit analysis, and more specific TESS follow-up analysis. Moreover, the pipeline tags exposures where any of the pixels inside a photometry aperture for a star are close to the linearity limit of the CCD. The exposures where either the target star or any of the comparison stars may be saturated are excluded from the analysis. The uncertainty on the data points is inferred using a rolling median with sigma clipping equivalent to the flux error determined by the PDCSAP pipeline for Kepler and TESS data. We estimate the average white noise scatter during the fitting in the pipeline, assuming the uncertainties do not change significantly during the night.

For the analysis, we used 48 nights of M2 between 2018 and 2020, and 11 nights of M1 between 2016 and 2020, with an observation of 1−3 h per night for both instruments. Thus, we had a total observation time of four days per filter for the M2 instrument; we refer to Table A.1 for details.

We also consider the K2 data set of TRAPPIST-1 (Howell et al. 2014). Paudel et al. (2018) analyzed the short-cadence data of the K2 Campaign 12 (Gilliland et al. 2010) and report 39 flares. An overview of the data of TRAPPIST-1 for M1, M2 and K2 is given in Table A.1. In Figs. 1 and 2 both flare light curves of M2 are shown in all four photometric bands.

Fig. 1 Light curve in all four filters of MuSCAT2 (M2) of Flare 1: the flare is visible in the g, r, and i filter, whereas the S/N in the zs-filter is too low. To better distinguish between the four M2 channels, we add constant values to the flux. The time is given in days in reduced Julian date (RJD), RJD = BJD − 2454 833 days where BJD is the barycentric Julian date.

3 Methods

3.1 Flare detection

For the flare detection we relied on the AltaiPony package (Ilin 2021). The selection criteria were adopted from Chang et al. (2015), where every outlier with three subsequent data points exceeding the local σ-level by a factor of three, was considered a flare candidate. σ represents the local standard deviation. To discriminate the real flares among the candidates found by the algorithm, the standard method is a flare injection recovery (hereafter FLIR), which constrains the probability of a given candidate being recovered in a given light curve. However, as the M2 data for TRAPPIST-1 suffer from strong variation in the photometric scatter due to weather, moon and telescope issues, we would have to perform the FLIR on every single observation night independently. For this reason, we first detected the flares with AltaiPony and afterwards, instead of performing FLIR on all light curves, we visually selected the flare candidates and adopted FLIR for M2 example light curves for each pass-band. The results of the FLIR of the M2 example light curves were used to explore the lower detection limit in bolometric energy. For M2, we obtained a detection limit of E_bol = 0.74 × 10³⁰ erg, see Appendix A.

Fig. 2 Light curve in all four filters of MuSCAT2 (M2) of Flare 2: the flare is visible in the g- and r-filter, whereas the S/N in the i- and zs-filter are too low. Time and flux are as indicated in Fig. 1.

3.2 Flare temperatures

The color temperature of a stellar flare is defined as the effective continuum temperature associated with a black body inferred from its spectral properties (e.g. Howard et al. 2020, and references therein). To estimate the color temperatures for MuSCAT flares we applied two different methods. First, the MuSCAT instruments provide the opportunity with their different photometric filters to measure the spectral energy distribution (SED) of flare events, that is the energy emitted by the flare per unit area, time and wavelength. We attempted to extrapolate the SED leaving the effective temperature and the flaring area as free parameters. The flaring area parameter a is given as a fraction of stellar radius and is the part of the stellar surface where the continuum emission of the flare is emitted. We assume it to be circular. The applied black body model can be found in Eq. (C.6). To avoid confusion, we refer to the effective temperatures estimated from the SED as SED temperatures T_SED. Second, we followed the methodology from Howard et al. (2020) cf. Sect. 5, to infer the flare color temperature using only two different filters. The global flare temperature T_glob is defined as the black body effective temperature linked to the total amount of flux of the flare, whereas the peak temperature T_peak is the temperature associated with the average temperature inside the FWHM of the flare (Howard et al. 2020). The temperatures can be estimated for each multi-color data point yielding a time evolution in flare temperatures.

As we did not have any simultaneous flare observed with K2 at different bands, this analysis was done only for the MuSCAT instruments.

3.2.1 SED temperature method

1. We computed the flux of TRAPPIST-1 in each MuSCAT filter using the corresponding segments of the mega-MUSCLES semi-empirical model (Wilson et al. 2021) and calibrated the data to absolute fluxes using the Gaia DR2 catalog (Evans et al. 2018). For the calibration we followed Ilin et al. (2021, cf. Sect. 2.4.3). First, we integrated the mega-MUSCLES spectra over the Gaia G-band response function, and normalized the SED to the Gaia G-band flux for TRAPPIST-1. Integrating over the MuSCAT-band response curves, we obtained the fluxes in the MuSCAT2-bands.

2. We then fit the relative flare light curve with the flare template aflare1 (Davenport et al. 2014). From the best-fit flare profile (cf. Sect. 3.3), we calculated the percentage flux change of TRAPPIST-1 during the flare.

3. The flare SED was then obtained by multiplying the percentage flux change of the flare with the quiescent flux for any MuSCAT filter.

4. We fit the flare SED with a black body model leaving the effective temperature and the flaring area as free parameters, cf. Appendix C.1.

This task can be performed for both the total flare flux, that is, the flux arising from the combined rise, peak and decay phase, and the peak flare flux, which is the flux emitted at the peak of the flare.

3.2.2 Two-filter temperature method from Howard et al. (2020)

1. We computed the radiation spectrum of a black body with temperature T_eff as a function of wavelength λ. We multiplied all filters separately with the black body spectrum and integrated them over the M1/M2 wavelength range to obtain the passband-specific flux. To account for the filter sensitivity we normalized the fluxes by dividing each calculated passband-specific flux by the total filter throughput.

2. We took the ratios R_j between the flux observed in all different filters, namely M2 g, r, i and z_s where j ∈ [(g/r), (g/i), (g/z_s)]. By convention, we took every single ratio such that the flux of the filter with a higher central wavelength was in the denominator and the flux of the filter with a lower central wavelength is in the numerator. In total we had three ratio functions to probe the color temperatures [R_(g/r), R_(g/i), R_(g/zs)].

3. This process was repeated for effective black body temperatures T_eff = [50, 50 000] K with 5 K steps.

4. To compute the effective black body temperatures of the observed flares, we inverted the different ratio functions such that we had the temperature as a function of the respective filter ratio T_eff (R_j).

In Fig. 3, the ratio functions R_j are plotted for all different filter combinations for the MuSCAT filters as a function of effective black body temperature T_eff. The R_j(T_eff) functions indicate that the retrievable color temperature information is limited. Where R_j(T_eff) is getting asymptotic, we reach the R_j specific sensitivity limit; see Table C.1. Beyond these limits, we cannot give reasonable estimates of the color temperature. In comparison, the ratio functions approach the asymptote faster if the total covered wavelength range between the two filters is smaller.

Following Hawley et al. (2003); Howard et al. (2020); Castellanos Durán & Kleint (2020), we calculated the wavelength-specific flux F_λ, $$F_{\lambda }=\frac{A_{\text{flare}}{B}_{\lambda }\left(T_{\text{flare}}\right)}{d^2}$$(1)

where A_flare is the area of the flare, B_λ(T_flare) is the black body spectrum of the flare with color temperature T_flare, and d is the distance of the star. A_flare does not depend on λ (Günther et al. 2020; Howard et al. 2020), taking the ratios between two passband fluxes yields: $$\frac{F_{{\lambda }_1}}{F_{{\lambda }_2}}\approx \frac{B_{{\lambda }_1}\left(T_{\text{flare}}\right)}{B_{{\lambda }_2}\left(T_{\text{flare}}\right)}=R_j.$$(2)

With Eq. (2) we used the ratio functions to infer the temperatures for the observed flares. Following Howard et al. (2020), this could be done in three different ways. We first investigated the flare temperature epoch by epoch, and then the global flare temperature. Finally, and importantly, we examined the flare peak temperature. We employ the calculation of all three measures and probe their consistency with the SED temperatures from the total flux and the peak flux.

Fig. 3 Ratio R of integrated fluxes observed in the MuSCAT bandpasses uniquely determines the effective black body temperature of a flare. A black body of temperature Teff is separately multiplied by the spectral response functions of the MuSCAT filters to produce the ratios Rj. We obtain different sensitivities for each filter pair: R(g,zs) is theoretically more sensitive in temperatures below 5000 K while R(g,zs) may be sensitive to temperatures up to 50 000 K.

3.3 Equivalent duration

The equivalent duration (ED) is defined as the amount of time it would take the star in its quiescent state to produce the same amount of energy that is released in a flare event (Gershberg 1972; Hunt-Walker et al. 2012). Thus, it is the time integral of the dimensionless flux, $$\text{ED}={\displaystyle \int \frac{F_{flare}-F_{\text{q}}}{F_{\text{q}}}\text{d}t,}$$(3)

where F_q is the quiescent flux of the star and F_flare is the flux observed during the flare. To estimate the equivalent duration, we fit the flare template aflare1 (Davenport et al. 2014) to each selected flare and thus acquire a best-fit approximation using a Markov chain Monte Carlo (MCMC) approach (Foreman-Mackey et al. 2013). The analytic expressions for the flare template are given in Appendix B.1. We used 5000 steps in the MCMC chain with 40 walkers, and discard the first 1000 steps as the burn-in phase. To put constraints on the parameters for the flare template (amplitude, peak time and full width at half maximum (FWHM)) we define Gaussian prior probability distributions, see Appendix B.2. After the fitting process, we integrate the flare template with the best-fit parameters for each flare event over time and thus obtain the ED, making use of the trapezoidal sum for integration. In Fig. 4, a flare profile fit is shown for the r-band.

Stellar flares are often multi-peaked (Günther et al. 2020) or show quasi-periodic pulsations (QPP; Mathioudakis et al. 2003; Mitra-Kraev et al. 2005). We focused in this work on modeling classical single-peaked flare events and disregard periodic oscillations in the decay phase in our analysis because both would not affect our temperature estimation; cf. Sect. 5.5.

Fig. 4 Flux profile fit of Flare 2 with M2 r-filter data: the orange curve indicates the best-fit using the flare template aflare1. The black lines show 100 random samples from the Markov chain. The lower panel shows the residuals between the observation and the model. The corresponding posterior probability distribution of the MCMC-fit is shown in Fig. B.1. The lower panel reveals systematic errors in the residuals which might be due to superimposed flares or oscillations. This is further discussed in Sect. 5.

4 Results

We found no flares in the M1 data and two flares were observed in all four filters in the M2 data. Flare 1 had a bolometric energy of E_bol = 2.99 × 10³¹ erg and occurred on August 25 2020 and Flare 2 had a bolometric energy of E_bol = 1.51 × 10³¹ erg and occurred on October 15 2020. Both flares are shown in all passbands in Figs. 1 and 2, respectively.

4.1 SED temperatures

To infer the SED for our M2 flares, we used the methods described in Sect. 3.2.

The inferred SEDs for both flares are shown in Fig. 5. The uncertainties were propagated from the best-fit of the flare flux profiles and the uncertainties of the mega-MUSCLES spectrum, as well as the uncertainties on the Gaia flux. To find a black body temperature consistent with the empirically measured filter fluxes, we used our black body model (c.f. Eq. (C.6)) and employed a fit to the data using emcee (Foreman-Mackey et al. 2013). We used 40 walkers and 25 000 steps in the chain and discarded the first 5000 steps as the burn-in phase. The convergence was checked visually by plotting the chains against step number and by computing the integrated autocorrelation time (Sokal 1996; Goodman & Weare 2010; Foreman-Mackey et al. 2013). We report autocorrelation times of roughly 35 steps for each trial. The used prior probability distributions can be found in B.1. In Fig. 6, the combined posterior probability distributions for both flare samples are shown and reveal a satisfactory resolution of the expected degeneracy between the area parameter and the flare temperature. The posterior distributions overlap marginally for the flare area parameter but indicate that the flare temperatures are indeed different between the two flares. The best-fit temperatures are $T_{\text{SED}}={7940}_{-390}^{+430}$ K for Flare 1 and $T_{\text{SED}}={6030}_{-270}^{+300}$ K for Flare 2. The same task is performed for the peak temperatures yielding $T_{{\text{SED}}_{\text{p}}}={13650}_{-1250}^{+1550}$ K for Flare 1 and $T_{{\text{SED}}_{\text{p}}}={8300}_{-550}^{+700}$ for Flare 2. The best-fit area parameters to the SEDs are $a={0.266}_{-0.021}^{+0.023}$ for Flare 1 and $a={0.171}_{-0.017}^{+0.019}$ for Flare 2. And for the peak SED $a_{\text{p}}={0.536}_{-0.070}^{+0.062}$ for Flare 1 and $a_{\text{p}}={0.213}_{0.031}^{+0.027}$ for Flare 2. The coverage of the stellar surface is represented in Table 1 and indicates that the flares cover up to 28% of the surface of TRAPPIST-1 at the peak and in the combined rise, peak and gradual phases up to 8%.

Fig. 5 SED for the two observed flares in M2. The SEDs are fit with a model where the temperature Teff and the area parameter a are left as free parameters using emcee again. The flux in zs-filter is slightly inconsistent with the best-fit black bodies for both flares. The uncertainties of the model are indicated with shaded areas. The uncertainty bars in wavelength should be understood as the respective filter widths. The SED for Flare 2 is multiplied by a constant factor of 2.5 for a better visual comparison between the two SEDs. The purple and orange lines indicate 100 randomly drawn samples from the Markov chain and mark the uncertainty of the fitted dashed best-fit curves.

4.2 Flare temperature by epoch

The flare temperature was calculated epoch-by-epoch to acquire the time evolution of the color temperature, see Fig. 7. Thus, fluxes were computed for every time step and compared with each other as presented in Sect. 3.2.

Because the observations in different passbands yield different cadences, we had to bin the observed flux to compare the theoretical black body fluxes between two filters. The uncertainty in each time step was propagated from the uncertainty of the filter-specific fluxes to the ratio functions.

To obtain a smoothly varying model for our epoch-by-epoch color temperature, we made use of our best-fit flare flux profile model. We fit the flare template aflare1 to the flare flux profile for each flare and MuSCAT passband, obtaining a smooth flux model for the flare event. Subsequently, we performed the same steps as before to retrieve the color temperatures from these flux profiles, cf. Sect. 3.2. We used our color temperature model to compute the global flare temperature, cf. Sect. 4.3.

Fig. 6 Posterior probability distribution for both flare samples: The sample of Flare 1 is represented in blue and the samples of Flare 2 are in orange. The sampling is done in logarithmic parameter space such that both the temperature T and the area parameter a are given here in logarithmic scale.

4.3 Global vs peak flare temperature

If signal-to-noise ratio (S/N) is low, the temperature evolution is not a precise measure for the color temperature of the flare and other measures should be used (Namekata et al. 2020). A more sophisticated approach to determining the temperature of a given flare makes use of the total instead of the epoch-by-epoch flux. Using our best-fit approximation for the flare flux profile (see Sect. 3.3), we calculated the total flux in each filter by integrating the flare template aflare1 with our best-fit parameters obtained for each flare. We compared the acquired flux ratios for each filter pair and applied the ratio functions R_j to yield an estimate of the ratio-specific global flare temperature. The same can be done for the peak fluxes that is, the flux within the FWHM of the flare. Table 2 shows the global and the peak flare temperatures, respectively, for every possible filter combination in comparison to the derived SED temperatures. Also, we combined our results for all filters by weighting the temperatures by their uncertainties.

Fig. 7 Temperature evolution of Flare 1: using the ratio functions Rj, we infer for each time step a flare color temperature by comparing the observed flux ratios with the theoretical black body flux ratios in between two MuSCAT filters using the χ2 method. Different ratio functions Rj yield different temperature evolutions. The uncertainty range is indicated by black markers and is valid for all three temperature evolutions. The lines represent the smoothly varying model for the flare color temperature which is derived from the best-fit parameters of the MCMC fitting to the flare light curve using the flare template aflare1. The orange shaded area gives the overlaid uncertainty of all three models.

5 Discussion

Flare temperatures on active M-dwarf stars (dM3e-dM4.5e) were previously derived from spectroscopic and photometric observations (Kowalski et al. 2013; Johnson et al. 2021) – T_glob = 9000−10000 K and T_peak = 9000−14000 K. It is usually assumed that flares on ultra-cool dwarfs have similar temperatures (e.g. Kowalski et al. 2013; Paudel et al. 2018; Günther et al. 2020; Glazier et al. 2020, and references therein). In this work, we showed with real photometric data that this assumption cannot be transferred to TRAPPIST-1, an ultra-cool M-dwarf with spectral type M7.5V. Thus raises the question of whether or not flares on late M-dwarfs should in general be modeled with cooler temperatures than early M-dwarfs.

5.1 SED versus two-filter method

The results from both the SED and the two-filter method imply that the global and peak temperature for the observed TRAPPIST-1 flares are lower than expected from the empirical flare temperatures used in modeling the energy budget of a flare (e.g. Segura et al. 2010; Kowalski et al. 2013; Tilley et al. 2019; Howard et al. 2020, and references therein), except for the peak temperature of Flare 1, where the temperature from both applied methods is consistent with the literature. The two methods are consistent within a 2 − σ−confidence interval; see Table 2 for details. The color temperature evolutions are shown in Figs. 7 and C.4. For both flares, the different ratio functions are consistent within the uncertainties. Comparing both methods, the clear advantage of the SED method is the simultaneous use of all filter fluxes and the simultaneous fit of the area parameter. As a consequence, the SED method is not affected by the disadvantage that we do not have a blue filter in our data set. The uncertainties on the two-filter methods could be improved by MCMC sampling of the theoretical black body flux values rather than using a simple χ² reduction. However, this is beyond the scope of this paper as the two-filter method was only meant to prove the temperature consistency of our SED method.

For Flare 1, the peak temperature for both applied methods is between 9000 and 14 000 K and is therefore consistent with the literature. For Flare 2, the peak temperature lies within 2-σ from the literature value for both methods. The reason for these consistent peak temperatures could be twofold. The fitting of the temperatures becomes more difficult toward higher temperatures because of the lack of a bluer photometric filter. This is partly monitored by the larger uncertainties that we obtain for the peak SED temperatures in comparison to the global SED temperatures. On the other hand, if the flare has high peak SED temperatures but lower global SED temperatures, this could imply that the cooling mechanisms are different in the atmosphere TRAPPIST-1 in comparison to earlier M-dwarfs.

5.2 SED temperature uncertainties

For the SED method, the estimated uncertainties were inferred by Gaussian uncertainty propagation taking the uncertainty from the flare flux profile and the uncertainty on the mega-MUSCLES spectrum. The marginalized uncertainties on the temperature were estimated directly by taking the 84th and 16th percentiles of the MCMC samples. Observing Fig. 5 closely, the uncertainties on the z_s-flux are the largest, extending over the whole SED flux regime for Flare 2. The z_s-flux is more sensitive to the flux-profile-fitting component, because of the low flare S/N in the z_s-band. For Flare 2 the low S/N prevented an accurate fit to the flare flux profile in the z_s-band and therefore we obtain large uncertainties.

Strong emission lines can also contribute to uncertainties. For instance, the peak flare flux could be contaminated by Ha and H_β emission by up to 10% and 8.8% of the total flux (Kowalski et al. 2013), respectively. Taking the maximum of the Ha contribution into account, the r-filter flux was reduced for the global calculation and for the peak calculation by roughly 30% and 26%, respectively. Indeed, this correction for the Hα contribution results in higher flare temperatures because the r-filter flux is corrected toward lower flux values, meaning that the SED is more likely to be fit with higher temperatures, unless another line contamination is considered. This could be compensated by the H_β emission in the g-filter. We reduced the flux in g-filter by approximately 50% and 20% for the total and peak calculations, respectively. As a consequence of the flux reduction in the g-filter the overall shape of the SED is shallower for both cases such that it is more likely to fit with lower temperatures than that of the reported black body. Contamination from both lines is probably is present in our data and therefore they will cancel each other out to a certain degree. However, it is impossible to correct the SED for the real Hα and H_β contributions. By omitting line emission in our calculations our uncertainties on the inferred black body temperatures are increased, but also it does not indicate a clear bias since the strongest line emission are prone to cancel each other out to a certain degree in the fitting procedure.

5.3 Flare location

The position of the flare on the stellar surface has an impact on the received flux and therefore might alter the SED. Because of the rotation of TRAPPIST-1 with period P ∈ [1.4, 3.295] days (Gillon et al. 2016; Vida et al. 2017; Luger et al. 2017), the flare might move over the stellar surface during its emission. The observed flares have durations of less than 30 min and therefore their angular movement on the surface is Δα ∈ [4.5, 1.9]°. The projected angles are even smaller, such that the time was not sufficient for the observed flares to move from the limb to the center or vice versa. Thus, we assumed that the projected surface area of TRAPPIST-1 is constant. Applying simple laws of limb darkening (Claret 2000), it was evident that the flux at the limbs reduced already at a line-of-sight angle of roughly 60 degrees to half its value at the center. Because the M2 data are given in relative units we used the mega-MUSCLES Spectra for TRAPPIST-1 to derive the SED of the flare. Bluer colors are more affected by limb darkening; see Fig. C.5. As discussed in Sect. 5.4, the flare temperature determines the slope of the SED and therefore the reduction of the flux close to the limb can be neglected. If we assume that the flares occurred close to the limb, the shape of the SED would result in a higher temperature than if they occurred in the vicinity of the center; see Fig. C.5. Therefore, we assumed that our best-fit temperatures are consistent with lower temperatures and independent of the exact position on the stellar surface.

5.4 Black body model

Our model relies on two main characterizing parameters, the flare temperature T_eff and the flare area parameter a; see Appendix C.1. The black body distribution determines the shape and slope of the model. Here, a is the constant of proportionality between the observed flux and our model. Using the MuSCAT data directly, the flare area parameter would not be a measure for the real covered radius fraction, because the MuSCAT data are given in relative units. We therefore calibrated our observation to the Gaia absolute fluxes. Qualitatively, the flare area parameter shows the expected behavior. Flare 1 emits more flux and therefore has also a larger flare area parameter than Flare 2. Because the flare temperature is the parameter that determines the shape of the SED, our model can still determine T_eff satisfactorily even if no absolute calibration is done but only the relative MuSCAT calibration. We tested this behavior for both flares obtaining the same temperatures but different flaring area parameters. Also, we performed a model sensitivity test for Flare 1 by multiplying the flare flux by arbitrary constant factors; see Fig. C.6. Multiplying the SED with constant factors resulted in a change in the flare area parameter but left the temperatures unchanged.

5.5 Residuals in flare profiles

In Fig. 4 the fit shows residual systematic errors. This could be due to complex flare events with multiple peaks or QPPs. Any superimposed flare would increase flux. As discussed in the previous section, the amount of flux is regulated by the flare area parameter and therefore a superimposed flare would not contaminate our flare temperatures, unless it changes the shape of the SED. Additionally, we employed a sinusoidal fit to the residuals with period P = 29.02 ± 6.53 min. Typical durations for QPPs range from several seconds to tens of minutes for M-dwarf flares (e.g. Van Doorsselaere et al. 2016; Million et al. 2021). Because our period is consistent with the typical QPP periods we considered our systematic errors to be explained by a QPP. The flux oscillations do not imply the flare temperatures and therefore the main result of our paper remains valid. For this reason, we neglected the observed QPP for our calculations.

5.6 Flare temperatures per spectral type

Figure 8 shows the flare continuum temperatures for TRAPPIST-1 and other M- and L-dwarfs from different photometric and spectroscopic studies (Mochnacki & Zirin 1980; Hawley & Pettersen 1991; Kowalski et al. 2013; Gizis et al. 2013; Howard et al. 2020). The upper panel reveals a possible trend for low-mass stars. The hottest flares seem to get cooler toward lower masses, such that the higher right corner is not populated by flares. This could be explained by a selection bias. As Howard et al. (2020) pointed out in their energy-temperature relation, more energetic flares seem to be hotter; see Fig. 9. These higher energetic flares are more likely for earlier type M-dwarfs (Howard et al. 2019; Günther et al. 2020) and therefore hotter flare temperatures T_eff > 10 000 K could be caught more frequently for early M-dwarfs than for later types like TRAPPIST-1.

Toward the tail of the low-mass-stars, there are only three flare temperatures measured. All three of them have bolometric flare energies between 10³¹ and 10³² and are in the lower energy tail of the energy–temperature relation. Therefore, it is difficult to verify the trend with spectral type, as there are not yet enough measurements of ultra-cool dwarf flares. Later type stars turn fully convective at roughly spectral type M4V and have cooler atmospheres, i.e. surface temperatures of below T_eff = 3200 K (Chabrier & Baraffe 1997; Reid & Hawley 2013). If the trend with spectral type is physical it might be due to a different flare generation process for fully convective stars in comparison to stars with radiative cores. Also, the heating mechanisms of the lower chromosphere driven by the accelerated particles could be different in late M-dwarfs. Without further observation and modeling of ultra-cool M-dwarf chromospheres, this discussion remains speculative.

Figure 9 shows the energy-temperature relation including the two TRAPPIST-1 flares for bolometric energies. We fit a power law of the form log₁₀(T_eff) = α · E_bol + β with two coefficients a and β to the Howard et al. (2020) sample, including the two TRAPPIST-1 flares, and obtain for the total temperatures α_tot = (6.98 ± 0.03) × 10⁻²,β_tot = 1.5 ± 0.4 and for the peak temperatures α_peak = (9.67 ± 6) × 10⁻², β_peak = 0.8 ± 0.7. Both peak and global temperatures of the TRAPPIST-1 flares are consistent with the uncertainty of the relation. However, we argue that the correlation is complete only if the flaring area is included in the relation because we expect that the heating process is affected not only by the energy released from the magnetic field but also by the size of the area that is heated. This could explain the large observed scatter especially for the peak temperatures. To obtain a correlation between energy-temperature and area, we propose to measure area and temperature simultaneously and dependent on each other in future similar studies; for example, by taking the dominant cooling mechanism in the chromosphere we could model how the temperature is distributed at the flare footpoints and how the area evolves then with time in dependence of the temperature. However, this is beyond the scope of the present work.

From the density map in the lower panel of Fig. 8, it follows that cool flares T_eff < 8000 K are also present in early spectral types and are not exclusively occurring on late M-dwarfs. Nevertheless, it appears that they are less frequent than flare temperatures around T_eff = 10 000 K. Following the energy-temperature relation, surprisingly, lower flare temperatures are not denser on the map, as they should occur more frequently. This could be explained by selection effects due to observational limitation for example.

Figure 10 shows a histogram of measured total temperatures from the same sample as in Fig. 8. The distribution reveals the peak at 10 000 K as expected. Our TRAPPIST-1 flares with 6030 and 7940 K, respectively, seem less likely but not at all unlikely. The histogram shows that the numbers of flares in a given temperature range increase, until 11 000 K is reached. Afterwards, there appears to be a cut. Between 11 000 and 12 000 K the observed numbers decrease strongly. It appears that flares with T_eff = 10 000−11 000 K are the most likely observed ones and flares with temperatures below 10 000 K appear to be more likely than flares with T_eff > 11 000 K. This could be explained by the empirical energy-temperature relation mentioned above. Flares with higher energy are less likely to occur and therefore also flares with higher temperatures are less frequent.

Taking the energy-temperature relation and the flare frequency distribution (FFD; see following section), we predict that we would have to observe TRAPPIST-1 for at least 15 days to find a flare with temperatures above 10 000 K. Taking this into account, we conclude that the flares observed on TRAPPIST-1 show lower temperatures than expected when compared to the average temperatures from previous studies (Hawley et al. 2003; Shibayama et al. 2013; Kowalski et al. 2013; Howard et al. 2020). However, if the flare frequency distribution and the energy-temperature relation are considered it is not surprising that we observed lower flare temperatures for TRAPPIST-1.

Fig. 8 Spectral type and stellar effective temperature versus total flare temperature for different photometric and spectroscopic studies: The upper panel uses the stellar effective temperature if available from the TESS input catalogue (Stassun et al. 2019). For the Kowalski sample, we use conservative uncertainties of ±500 K. The lower panel shows a density map of flare temperatures per spectral type that is consistent with the upper stellar effective temperatures. The flare temperature axis is represented in log-scale. The two TRAPPIST-1 flares are overplotted and the blue horizontal lines represent the average temperatures from the two biggest temperature samples.

Fig. 9 Energy-temperature relation in log-log representation: in the upper panel the global flare temperatures are plotted versus the bolometric energy, whereas in the lower panel the peak temperatures are used. The flares from Howard et al. (2020) are indicated in blue and our two TRAPPIST-1 flares are shown in orange. The blue dashed line indicates a power-law fit to the data.

Fig. 10 Histogram of total flare temperatures from the spectroscopic and photometric observations shown in Fig. 8 with the peak of the distribution close to 10 000 K.

5.7 Flare frequency distribution and habitability

One goal of our analysis was to update the FFD of TRAPPIST-1 (Paudel et al. 2018; Glazier et al. 2020) and show the implications of a different flare temperature. The FFD is the cumulative rate of flares per time unit (Gershberg 1972) and is frequently described by a power law (Davenport 2016; Paudel et al. 2018; Günther et al. 2020; Glazier et al. 2020). In log-log space, this power law has a linear appearance, $${\log }_{10}v=\alpha {\log }_{10}{E}_{\text{bol}}+\beta ,$$(4)

where ν is the occurrence rate of the flares in a given time unit, E_bol is the bolometric energy α indicates the slope and β is the y-intercept in log-log space. The bolometric energies were calculated as presented in Appendix B.3 and we adopted a value of 7000 ± 500 K for the flare temperature, which is consistent with the results of this paper.

First, we took the flares found by Paudel et al. (2018) in the K2 light curve and computed the bolometric energies by adopting their EDs following Sect. B.3.

The updated FFD for the combined information of M1, M2 and K2 is shown in Fig. 11. We observe a truncated power law with a turnoff at roughly 10³⁰ erg. This is a result of flares with lower energies close to the K2 detection limit as discussed for Kepler by Hawley et al. (2014). The theoretical detection limit for M2, shown by the vertical dashed line, is given by the g-band; see Sect. 3.1. The orange squares indicate the two M2 flares, while the black points represent the K2 flares (Paudel et al. 2018).

Using the power law ansatz for the flare frequency Eq. (4), we fit the flare frequency with the MCMC approach of Wheatland (2005) implemented in AltaiPony. The best-fit to the data is indicated by the red solid line, as well as 100 samples of the 2500 MC trials. By extrapolating the FFD to larger bolometric energies, we explored the FFD into the regime of low flare occurrence, that is higher energies. We analyzed potential danger to the atmospheres of the TRAPPIST-1-companions and possible abiogenesis. The blue dash-dotted line indicates the best-fit if the bolometric energies were calculated with a black body temperature of 10 000 K. Comparing the two best-fit solutions with the different adopted temperatures, we obtained α_7000K 1.45 ± 0.20, α_{10 000 K} = 1.46 ± 0.24. The higher temperature led to a slightly steeper slope. However, the impact is not significant because both slopes are consistent with each other. In the abiogenesis zone of an MS star, specific prebiotic chemical reactions are possible that enable the synthesis of ribonucleic acid (RNA; Rimmer et al. 2018; Günther et al. 2020). The flare frequency ν needed to trigger and sustain prebiotic chemistry is given by Günther et al. (2020); $$v\ge 25.5{\text{days}}^{-1}\left(\frac{{10}^{34}\text{erg}}{E_{\text{U}}}\right){\left(\frac{R_{\text{star}}}{R_{\odot }}\right)}^2{\left(\frac{T_{\text{star}}}{T_{\odot }}\right)}^4$$(5)

where E_U is the U-band energy, R_star the stellar radius and T_star the stellar surface temperature. In Fig. 11, the abiogenesis zone is plotted in green, adopting the surface temperature and radius of TRAPPIST-1 of T = 2648 ± 26 K (Wilson et al. 2021) and R = 0.1192 ± 0.0013 R_⊙ (Agol et al. 2021). Both the red and blue lines do not intersect with the abiogenesis zone, indicating, that the high energetic flux from TRAPPIST-1 flares is likely to be insufficient to trigger and sustain prebiotic chemistry. Glazier et al. (2020) came to the same conclusion.

The yellow shaded area in Fig. 11 marks the ozone depletion zone around TRAPPIST-1, where frequent high-energy flares may erode exoplanet atmospheres and lead to subsequent surface sterilization. Tilley et al. (2019) modeled the impact of M-dwarf flares on the atmosphere of an Earth analog planet. These authors argued that Earth’s ozone layer should erode if flares with bolometric energies ≥10³⁴ erg hit the atmosphere at a frequency of ν₃₄ ≥ 0.4 day⁻¹. Günther et al. (2020) gave a lower limit for the frequency ν₃₄ ≥ 0.1 day⁻¹. We adopted the more conservative approach of Günther et al. (2020), which also allowed us to remain consistent with the analysis of Glazier et al. (2020). As the best-fit did not intersect with the ozone-depletion zone, we conclude that the UV flux is insufficient to deplete a complete Earth-like ozone layer.

Fig. 11 Flare frequency distribution for K2 and M2 data. The FFD for TRAPPIST-1 is updated by our flares indicated in orange. The black points are adopted from (Paudel et al. 2018). The best fit to the power law is chosen from 2500 MC trials, here only 100 random samples are shown to mark the uncertainty of the model. Also overplotted are the abiogenesis zone in green (Rimmer et al. 2018; Günther et al. 2020) and the ozone-depletion zone in yellow (Tilley et al. 2019; Günther et al. 2020) for TRAPPIST-1.

6 Conclusion

We used 59 nights of multi-color photometric observations of MuSCAT1 and MuSCAT2 to search for flares on TRAPPIST-1. We found two flares in the MuSCAT2 light curves. Moreover, we discovered that the black body temperatures associated with the total emitted flux are likely to be cooler than previously suggested for ultra-cool M-dwarfs. We inferred temperatures of $T_{\text{SED}}={7940}_{-390}^{+430}$ K for Flare 1 and $T_{\text{SED}}={6030}_{-270}^{+300}$ K for Flare 2. We obtained peak temperatures $T_{\text{SEDp}}=13{620}_{-1220}^{+1520}$ K for Flare 1 and $T_{\text{SEDp}}={8300}_{-550}^{+660}$ for Flare 2. The observed peak black body temperatures for TRAPPIST-1 also show marginally cooler temperatures for Flare 2, while they are consistent with the literature for Flare 1. It could be that different cooling mechanisms in the stellar atmosphere are responsible for this behavior and that cooling is more efficient toward later spectral types. This remains speculative without further spectroscopic measurements and modeling. Lower black body temperatures for flares lead to different UV surface fluxes and therefore have an impact on the habitability estimation on exoplanets around ultra-cool M-dwarfs, even though we showed the bolometric energies remain marginally affected. While using the 9000−10 000 K assumption one has to keep in mind, that flares with lower temperatures are more likely to occur and thus, we suggest using a flare frequency temperature distribution to account for the different flare temperatures in the future.

Furthermore, we conclude that, based on our data, it is not possible to verify whether the lower observed temperatures on TRAPPIST-1 are intrinsic to ultra-cool dwarfs or are the result of lower energies generated for these late-type stars, making it more likely to observe lower flare temperatures as suggested by the energy-temperature relation. Our prediction of roughly 15 days of observation to observe a flare with T_eff > 10 000 K is a considerable task for future work and would also show whether or not the energy-temperature relation can be expanded to the low-mass end of M-dwarfs.

Our SED method proves to be suitable for estimating flare temperatures and areas. It is a good trade-off between previous multi-color photometric- and spectroscopic approaches because of its efficiency and low cost. Even though spectroscopic methods enable the resolution of the specific line emission, we show that our model allows us to infer precise temperatures and areas associated with stellar flares.

Further multi-color observations are needed to put constraints on other M-dwarf flares in order to verify whether or not the observed behavior is indeed a general characteristic of ultra-cool M-dwarfs. The reported trend of lower flare temperatures toward later spectral types could be investigated by a thorough observation campaign of several late M-dwarfs and early L-dwarfs to populate the low-mass tail in Fig. 8. M-dwarfs with a surface temperature of ~3000 K are particularly suitable targets as they show enhanced activity in comparison to other M-dwarfs and earlier type stars (Günther et al. 2020). We recommend temperature measurements of similar stellar types to TRAPPIST-1 with bluer filters. Also, we encourage the spectroscopic observation of ultra-cool M-dwarf flares in order to infer Hα and Hβ contributions, for example, to improve the results.

Acknowledgements

A.J.M. would like to thank the anonymous referee for his/her thoughtful review of the work. This article is based on observations made with the MuSCAT2 instrument, developed by ABC, at Telescopio Carlos Sanchez operated on the island of Tenerife by the IAC in the Spanish Observatorio del Teide. This article was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. E.I. acknowledges support from the German National Scholarship Foundation. E.E.-B. acknowledges financial support from the European Union and the State Agency of Investigation of the Spanish Ministry of Science and Innovation (MICINN) under the grant PRE2020-093107 of the Pre-Doc Program for the Training of Doctors (FPI-SO) through FSE funds. N.N. acknowledges the support by JSPS KAKENHI Grant Number JP18H05439, JST CREST Grant Number JPMJCR1761. M.M. is supported by Grant-in-Aid for JSPS Fellows, Grant Number JP20J21872. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia, Innovación y Universidades through projects PID2019-109522GB-C53. The authors want to thank Maximilian Günther and Ward Howard for helpful discussions based on their expertise on stellar flares on M-dwarfs.

Appendix A Flare detection limits of the MuSCAT 2 instrument

To infer the detection limits, we generated a mock light curve, for each M2 filter by adopting normalized flux and adding Gaussian noise with a mean equal to the median of the respective noise level of all M2 observations for TRAPPIST-1. The typical observation time is calculated from Table A.1 by adapting the total observation time per filter and their averaged cadences. The results from FLIR are shown in Figure A.1, where the injection of 20 000 classical single-peaked events with the Davenport et al. (2014) template aflare I were visualized for each filter. The FLIR was executed iteratively, i.e. only one flare was injected per iteration in the light curve to avoid overlapping of injected flares. In Fig. A.1, the dark regions mark low recovery probability, while the yellow areas indicate high recovery probability¹.

We defined the detection limit as the region where the FLIR delivers 100% recovery probability, for the first time, because there we were certain that any flare event with these properties would be detected. This approach may seem too conservative but as the light curve is artificial, it seems appropriate to remain as conservative as possible. The detection limits in terms of bolometric energy, FWHM, and amplitude of the flare are listed in Table A.2. The ranges for FWHM and amplitude implied that the recovery probability was calculated in FWHM and amplitude bins, such that in the given bin with the border of the given ranges, we have 100% recovery probability. The bolometric energy was then calculated using the central value of these ranges by adopting the equations presented in Section B.3, where the uncertainty was computed by Gaussian error propagation and is constrained by the width of the FWHM and the amplitude for a given filter.

We also had a limit for the detection of long-duration and hence high-energy flares. If the FWHM of the flare exceeds the duration of the observation (1.5−2 h), the flare was not completely resolvable and the fitting process might not converge as described in Section 3.3. However, this limit was difficult to quantify because of its multidimensionality. It depends not only on the FWHM of the flare but also on its peak time and was therefore not quantified in this work. Moreover, the detection of such highly energetic flares was unlikely in our data set. By looking at the FFD of TRAPPIST-1 we identify that flares with bolometric energy of E_bol > 10³⁴ occur only once in 100 days. As our total observation time is roughly 5 days with the combined MuSCAT observation we neglected the detection limit of long-duration flares.

Fig. A.1 Flare injection recovery via heat maps for each MuSCAT 2 (M2) filter, Panels (a),(b),(c), and (d) correspond to the M2 g, r, i and zs photometric filters. The detection limits in amplitude and FWHM of the flares are for all bands visible and marked by the dark regions where the recovery probability is low.

Appendix B Flare energetics

Appendix B.1 Analytic flare template

Davenport et al. (2014) used for the generation of their flare template aflare1 flares that occurred on GJ1243 (M4.0Ve). The duration of used flares for the generation of the template was in the range of 20–75 min. Therefore, both of our observed flares with durations of over t > 30 min fit well in the selected flare range. The analytic formula for the flare flux profiles relies on a combination of polynomial and exponential fits. The rise phase, i.e., −1t_1/2 ≤ t ≤ 0, is described by a polynomial function (Davenport et al. 2014): $$\begin{array}{lll}{F}_{\text{rise}}\hfill & =\hfill & 1+1.941\left(\pm 0.008\right)t_{1/2}-0.175\left(\pm 0.032\right)t_{1/2}^2\hfill \\ \hfill & \hfill & -2.246\left(\pm 0.039\right)t_{1/2}^3-1.125\left(\pm 0.016\right)t_{1/2}^4\hfill \end{array}$$(B.1)

with t_1/2 = x − t_peak/FWHM being the full width in time at the maximum of the flare and x the observational time (Kowalski et al. 2013). The decay phase, i.e., t > 0, is described by (Davenport et al. 2014): $$\begin{array}{lll}{F}_{\text{decay}}\hfill & =\hfill & 0.6890\left(\pm 0.0008\right)e^{-1.600\left(\pm 0.003\right)t_{1/2}}+\hfill \\ \hfill & \hfill & 0.3030\left(\pm 0.0009\right)e^{-0.2783\left(\pm 0.0007\right)t_{1/2}}\hfill \end{array}$$(B.2)

Therefore the full template reads: $$F_{\text{flare}}=\left(F_{\text{rise}}+F_{\text{decay}}\right)\cdot \left|ampl\right|$$(B.3)

where ampl is the amplitude in the flux of the flare.

Appendix B.2 Gaussian prior probability

For the MCMC fitting of the flare flux profiles, we needed to define prior probability distributions for the three parameters of the flare template aflare1. The estimate of the amplitude am_rec from AltaiPony was taken as the mean of the respective Gaussian prior and 0.05 in terms of normalized flux as a conservative standard deviation. The peak time t_peak was calculated as the argmax of the normalized flux for a given flare event. We therefore chose the standard deviation of the peak time to be one data point in each observation. The FWHM_est was only approximated very roughly by dividing the duration of the flare by the arbitrary value of four. For this reason the standard deviation for the FWHM was chosen to be 1/24 days and the prior was therefore rather uninformative. For an overview of the adapted Gaussian priors see Table B.1.

Appendix B.3 Bolometric flare energies

In this chapter, we present how the flare energies were calculated. We followed the bolometric flare energy calculation of (Shibayama et al. 2013; Günther et al. 2020), where the stellar luminosity and the best-fit flare flux profile are used.

Fig. B.1 Corner plot of MCMC flare profile fit of Flare 2 in the r passband. Shown are the correlations between the three parameters of aflare1 - peak time, full width at half maximum and amplitude. The dashed lines indicate 15th, 50th 84th quantiles.

First, we calculated the bandpass independent area A_flare(t) of the flare assuming the flare luminosity can be modeled as a black body with constant effective temperature of T_flare = 9000 ± 500K, $$A_{\text{flare}}\left(t\right)=\frac{F_{\text{flare}}\left(t\right)+F_{\text{q}}}{F_{\text{q}}}\pi R_{\text{star}}^2\frac{{\displaystyle \int R_{\lambda }{B}_{\lambda }\left(T_{\text{eff}}\right)d\lambda }}{{\displaystyle \int R_{\lambda }{B}_{\lambda }\left(T_{\text{flare}}\right)d\lambda }}$$(B.4)

where denotes the instrument response function, β_λ the black body radiation, T_eff the effective temperature of the star in its quiescent state, and R_star refers to the stellar radius. The effective flare temperature T_flare was chosen as a lower conservative limit in order to remain consistent with previous studies (Davenport 2016; Paudel et al. 2018; Günther et al. 2020) and references therein, but it was not a correct assumption, because the flare temperature not only changes with time but also depends on the flare phase. Also, we show in this work that the temperature associated with the whole flare event is lower than frequently adopted for ultra-cool M-dwarfs.

Having A_flare(t), we can compute the bolometric flare energy E_flare of the flare, using its corresponding luminosity, $$E_{\text{flare}}={\displaystyle \int L_{\text{flare}}\left(t\right)dt={\sigma }_{\text{sb}}{T}_{\text{flare}}^4{\displaystyle \int A_{\text{flare}}\left(t\right)dt}}$$(B.5)

We calculate the uncertainties of the calculation using Gauss uncertainty propagation. Integrations of this section were performed via the trapezoidal sum rule.

Appendix B.4 Filter-specific flare energies

We calculated not only the bolometric but also the filter-specific energy, in order to infer the detection limits. We obtained the filter-specific energies by adopting a conversion constant from the bolometric energies. Assuming the flare has a bolometric luminosity L_flare(t) and filter specific luminosity L_λflare(t), then: $$L_{\text{flare}}\left(t\right)={\sigma }_{\text{sb}}{T}_{\text{flare}}^4{A}_{\text{flare}}\left(t\right)$$(B.6) $$L_{\lambda \text{,flare}}\left(t\right)=L_{\text{flare}}\left(t\right)\cdot c$$(B.7)

where the flaring area A_flare(t) is not passband-specific, T_flare is the effective blackbody temperature from the flare and c denotes the dimensionless conversion constant. By using L = Fπd², we obtain: $$c=\frac{L_{\lambda ,\text{flare}}}{L_{\text{flare}}}=\frac{F_{\lambda ,\text{flare}}}{F_{\text{flare}}}$$(B.8)

with F_flare denoting the flux from the flare event. Because the bolometric energy associated with the flares was defined as in equation (B.5), we derived the passband-specific energy using equation (B.8): $$E_{\text{flare},\lambda }={\displaystyle \int L_{\lambda ,\text{flare}}dt=c{\displaystyle \int L_{\text{flare}}dt}}$$(B.9)

The calculated conversion constants for M1/M2 filters are given in Table B.2.

Appendix C Flare temperatures

Appendix C.1 Flare black body model

In this section, we show our used black body model for the fitting process described in Section 3.2. We start at the apparent brightness of the observer, $$l_{\lambda }\frac{L_{\lambda }}{4\pi d^2},$$(C.1)

with L_λ the luminosity and d the distance of the source. As the flare occurs on the visible side of the star, we integrate only the visible area from the star. Thus, we obtain $$L_{\lambda }=4{\pi }^2{R}_{*}^2{B}_{\lambda }\left(T_{*}\right),$$(C.2)

with R_* the radius of the star and B_λ(T_*) the black body spectrum of temperature T_*. Substituting back in to the apparent brightness, $$l_{\lambda }=\frac{R_{*}^2}{d^2}\cdot \frac{1}{\exp \frac{hc}{k_b{T}_{*}\lambda }-1}\cdot \frac{2h{c}^2}{{\lambda }^5}\cdot \pi .$$(C.3)

Fig. C.1 The semi-empirical mega-MUSCLES spectra for TRAPPIST-1 include emission and absorption lines. In comparison the black body temperature model for TRAPPIST-1 Teff = 2648 K is shown.

Under the assumption that a flare with an effective temperature of T_flare covers a circle with radius R_flare on the surface of the star, we can perform the same calculation. $$l_{flare}=\frac{R_{flare}^2}{d^2}\cdot B_{\lambda }\left(T_{\text{flare}}\right)\cdot \pi =\frac{{\left(a\cdot R_{*}\right)}^2}{d^2}\cdot B_{\lambda }\left(T_{\text{flare}}\right)\cdot \pi =\tilde{a}\cdot l_{\lambda }\left(T_{\text{flare}}\right),$$(C.4)

with ã = a² and a being the fraction of covered surface radius by the flare. In the optically thick case, we obtain the total apparent brightness, $$l_{\text{total}}=\tilde{a}{l}_{\lambda }\left(T_{\text{flare}}\right)+\left(1-\tilde{a}\right)\cdot l_{\lambda }\left(T_{*}\right).$$(C.5)

As we obtain the apparent flare brightness directly as a result of our fitting pipeline, we fit the obtained brightness independent of the stellar brightness. Thus, we obtain the black body model for a flare with effective temperature T_flare and fractional area parameter ã, $$l_{\text{flare}}=\tilde{a}\cdot \left(l_{\lambda }\left(T_{\text{flare}}\right)-l_{\lambda }\left(T_{*}\right)\right).$$(C.6)

Fig. C.2 The SED of the peak fluxes for the two observed M2 flares. The SEDs fit with a model where the temperature Teff was left as a free parameter using again emcee. The uncertainties of the model are indicated with 100 random samples from the chain.

Fig. C.3 Posterior probability distribution for both peak flare samples: The sample of flare 1 is represented in blue and the samples of flare 2 are in orange. The sampling is done in logarithmic parameter space such that both the temperature T and the area parameter a are given on logarithmic scale.

Fig. C.4 Color-temperature evolution for both observed TRAPPIST-1 flare with MuSCAT 2. In the left column, all temperature evolutions for the different filter ratio functions Rj(Teff) are shown, and in the right column, the same is done for flare 2. The orange line indicates the model drawn from the best-fit to the flare flux profile, with respective model uncertainties represented by the orange shaded areas.

Fig. C.5 Center-to-limb variation for the MuSCAT filters. The limb darkening is modeled with quadratic law (Claret 2000).

Fig. C.6 Model sensitivity test for the flare area parameter: three SEDs are shown for Flare 1. The blue SED indicates the SED with absolute calibration to Gaia and the orange and green curves represent fluxes with altered SEDs.

References

All Tables

All Figures

	Fig. 1 Light curve in all four filters of MuSCAT2 (M2) of Flare 1: the flare is visible in the g, r, and i filter, whereas the S/N in the zs-filter is too low. To better distinguish between the four M2 channels, we add constant values to the flux. The time is given in days in reduced Julian date (RJD), RJD = BJD − 2454 833 days where BJD is the barycentric Julian date.
In the text

	Fig. 2 Light curve in all four filters of MuSCAT2 (M2) of Flare 2: the flare is visible in the g- and r-filter, whereas the S/N in the i- and zs-filter are too low. Time and flux are as indicated in Fig. 1.
In the text

Fig. 3 Ratio R of integrated fluxes observed in the MuSCAT bandpasses uniquely determines the effective black body temperature of a flare. A black body of temperature Teff is separately multiplied by the spectral response functions of the MuSCAT filters to produce the ratios Rj. We obtain different sensitivities for each filter pair: R(g,zs) is theoretically more sensitive in temperatures below 5000 K while R(g,zs) may be sensitive to temperatures up to 50 000 K.

In the text

Fig. 4 Flux profile fit of Flare 2 with M2 r-filter data: the orange curve indicates the best-fit using the flare template aflare1. The black lines show 100 random samples from the Markov chain. The lower panel shows the residuals between the observation and the model. The corresponding posterior probability distribution of the MCMC-fit is shown in Fig. B.1. The lower panel reveals systematic errors in the residuals which might be due to superimposed flares or oscillations. This is further discussed in Sect. 5.

In the text

Fig. 5 SED for the two observed flares in M2. The SEDs are fit with a model where the temperature Teff and the area parameter a are left as free parameters using emcee again. The flux in zs-filter is slightly inconsistent with the best-fit black bodies for both flares. The uncertainties of the model are indicated with shaded areas. The uncertainty bars in wavelength should be understood as the respective filter widths. The SED for Flare 2 is multiplied by a constant factor of 2.5 for a better visual comparison between the two SEDs. The purple and orange lines indicate 100 randomly drawn samples from the Markov chain and mark the uncertainty of the fitted dashed best-fit curves.

In the text

	Fig. 6 Posterior probability distribution for both flare samples: The sample of Flare 1 is represented in blue and the samples of Flare 2 are in orange. The sampling is done in logarithmic parameter space such that both the temperature T and the area parameter a are given here in logarithmic scale.
In the text

Fig. 7 Temperature evolution of Flare 1: using the ratio functions Rj, we infer for each time step a flare color temperature by comparing the observed flux ratios with the theoretical black body flux ratios in between two MuSCAT filters using the χ2 method. Different ratio functions Rj yield different temperature evolutions. The uncertainty range is indicated by black markers and is valid for all three temperature evolutions. The lines represent the smoothly varying model for the flare color temperature which is derived from the best-fit parameters of the MCMC fitting to the flare light curve using the flare template aflare1. The orange shaded area gives the overlaid uncertainty of all three models.

In the text

Fig. 8 Spectral type and stellar effective temperature versus total flare temperature for different photometric and spectroscopic studies: The upper panel uses the stellar effective temperature if available from the TESS input catalogue (Stassun et al. 2019). For the Kowalski sample, we use conservative uncertainties of ±500 K. The lower panel shows a density map of flare temperatures per spectral type that is consistent with the upper stellar effective temperatures. The flare temperature axis is represented in log-scale. The two TRAPPIST-1 flares are overplotted and the blue horizontal lines represent the average temperatures from the two biggest temperature samples.

In the text

	Fig. 9 Energy-temperature relation in log-log representation: in the upper panel the global flare temperatures are plotted versus the bolometric energy, whereas in the lower panel the peak temperatures are used. The flares from Howard et al. (2020) are indicated in blue and our two TRAPPIST-1 flares are shown in orange. The blue dashed line indicates a power-law fit to the data.
In the text

	Fig. 10 Histogram of total flare temperatures from the spectroscopic and photometric observations shown in Fig. 8 with the peak of the distribution close to 10 000 K.
In the text

Fig. 11 Flare frequency distribution for K2 and M2 data. The FFD for TRAPPIST-1 is updated by our flares indicated in orange. The black points are adopted from (Paudel et al. 2018). The best fit to the power law is chosen from 2500 MC trials, here only 100 random samples are shown to mark the uncertainty of the model. Also overplotted are the abiogenesis zone in green (Rimmer et al. 2018; Günther et al. 2020) and the ozone-depletion zone in yellow (Tilley et al. 2019; Günther et al. 2020) for TRAPPIST-1.

In the text

	Fig. A.1 Flare injection recovery via heat maps for each MuSCAT 2 (M2) filter, Panels (a),(b),(c), and (d) correspond to the M2 g, r, i and zs photometric filters. The detection limits in amplitude and FWHM of the flares are for all bands visible and marked by the dark regions where the recovery probability is low.
In the text

	Fig. B.1 Corner plot of MCMC flare profile fit of Flare 2 in the r passband. Shown are the correlations between the three parameters of aflare1 - peak time, full width at half maximum and amplitude. The dashed lines indicate 15th, 50th 84th quantiles.
In the text

	Fig. C.1 The semi-empirical mega-MUSCLES spectra for TRAPPIST-1 include emission and absorption lines. In comparison the black body temperature model for TRAPPIST-1 Teff = 2648 K is shown.
In the text

	Fig. C.2 The SED of the peak fluxes for the two observed M2 flares. The SEDs fit with a model where the temperature Teff was left as a free parameter using again emcee. The uncertainties of the model are indicated with 100 random samples from the chain.
In the text

	Fig. C.3 Posterior probability distribution for both peak flare samples: The sample of flare 1 is represented in blue and the samples of flare 2 are in orange. The sampling is done in logarithmic parameter space such that both the temperature T and the area parameter a are given on logarithmic scale.
In the text

Fig. C.4 Color-temperature evolution for both observed TRAPPIST-1 flare with MuSCAT 2. In the left column, all temperature evolutions for the different filter ratio functions Rj(Teff) are shown, and in the right column, the same is done for flare 2. The orange line indicates the model drawn from the best-fit to the flare flux profile, with respective model uncertainties represented by the orange shaded areas.

In the text

	Fig. C.5 Center-to-limb variation for the MuSCAT filters. The limb darkening is modeled with quadratic law (Claret 2000).
In the text

	Fig. C.6 Model sensitivity test for the flare area parameter: three SEDs are shown for Flare 1. The blue SED indicates the SED with absolute calibration to Gaia and the orange and green curves represent fluxes with altered SEDs.
In the text