Giant white-light flares on fully convective stars occur at high latitudes Free

ABSTRACT

White-light flares are magnetically driven localized brightenings on the surfaces of stars. Their temporal, spectral, and statistical properties present a treasury of physical information about stellar magnetic fields. The spatial distributions of magnetic spots and associated flaring regions help constrain dynamo theories. Moreover, flares are thought to crucially affect the habitability of exoplanets that orbit these stars. Measuring the location of flares on stars other than the Sun is challenging due to the lack of spatial resolution. Here we present four fully convective stars observed with the Transiting Exoplanet Survey Satellite that displayed large, long-duration flares in white-light which were modulated in brightness by the stars’ fast rotation. This allowed us to determine the loci of these flares directly from the light curves. All four flares occurred at latitudes between 55° and 81°, far higher than typical solar flare latitudes. Our findings are evidence that strong magnetic fields tend to emerge close to the stellar rotational poles for fully convective stars, and suggest that the impact of flares on the habitability of exoplanets around small stars could be weaker than previously thought.

1 INTRODUCTION

Cool stars, i.e. stars with an outer convective envelope, are known to display a variety of magnetically driven phenomena, such as starspots, flares, and coronal mass ejections, collectively called magnetic activity (Schrijver & Zwaan 2000). Sun-like stars possess a radiative core, and their magnetic activity is understood to be driven by a magnetic dynamo of the α − Ω type (Parker 1955; Roberts 1972). For low-mass stars that are fully convective (spectral type M4V and later, Stassun et al. 2011), a type of dynamo must be operating that does not require the presence of a radiative core. Surprisingly, the relation between rotation and magnetic activity for stars on either side of the fully-convective boundary is the same in both X-ray (Wright & Drake 2016; Wright et al. 2018) and Hα (Newton et al. 2017) observations, suggesting that the dynamo operating in fully convective stars has a similar efficiency as the one in our Sun. One important difference between Sun-like stars and low-mass M dwarfs is that Sun-like stars spin-down quickly, while M dwarfs remain fast rotators much longer (Irwin et al. 2011; Newton et al. 2016). Fast rotation fuels enhanced magnetic activity and therefore frequent flaring, which is exploited in studies of stellar flares in M dwarf stars down to the bottom of the main sequence (Stelzer et al. 2006; Robrade, Poppenhaeger & Schmitt 2010; Schmidt et al. 2016; Paudel et al. 2018).

On the Sun, emerging magnetic fields drive flares, and stronger fields drive more energetic flares, making them highly localized markers of dynamic magnetic field concentrations (Priest & Forbes 2002). In optical light curves, flares manifest as sudden brightenings that last between a few minutes and several hours. The onset of a flare is characterized by a fast rise in stellar flux which peaks in a sharp maximum, followed by an exponential, and then a more gradual, decay back toward the quiescent state.

The location of flares on the stellar surface constrains the surface field strength and topology, and tests predictions from stellar dynamo theories. On the Sun, the positions of large flares follow the cyclic latitude variations of sunspots (Zhang et al. 2007). The latter is known as the butterfly diagram, and is crucial to the understanding of the solar dynamo (Gnevyshev 1977). Fast-rotating young solar-like stars and low mass stars have been predicted to form magnetic spots close to the rotational poles (Schuessler & Solanki 1992; Granzer et al. 2000; Gastine et al. 2013; Yadav et al. 2015), and this prediction has been confirmed by spectroscopic and spectropolarimetric observations (e.g. Strassmeier & Rice 1998; Marsden et al. 2006; Järvinen et al. 2007, 2008; Ișık et al. 2018; Barnes et al. 2015, 2017). If these spots are co-spatial with large flares, they may be classified as intensely active regions with a dynamic local field configuration.

Differences in flaring region loci can also, at least partially, explain the variety seen in flare light curve morphologies. The convex shape of a flare light curve in X-rays on V773 Tau (Skinner et al. 1997) could be well described by rotational modulation. Later, effects of rotation and eclipses have been put forward as an explanation for unusual flare light curves multiple times (Stelzer et al. 1999; Montmerle et al. 2000; Johnstone et al. 2012). However, only few examples of unambiguously located stellar flares exist to date (Schmitt & Favata 1999; Wolter et al. 2008; Peterson et al. 2010).

Here, we present observations of four giant flares that occurred on fast-rotating fully convective dwarfs detected in photometry from the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015). These flares are remarkable because they lasted longer than a full rotation period of the star so that the characteristic flare light curve was modulated by the flaring footpoint’s movement in and out of view on the stellar surface. In Section 2, we describe our selection of both the stellar and the flare sample, and explain how we determined stellar properties for the promising candidates. In Section 3, we introduce a model of large flares with rotationally modulated light curves, show how we can use it to infer the latitude of the flaring regions, and discuss the assumptions made in the model. We fit our model to the observations, and present the inferred flaring region latitudes in Section 4. We briefly discuss quasiperiodic pulsations as an alternative model, place our results in the context of current research on stellar surface magnetic fields and dynamo thoeries for fully convective stars, and explore their implications for exoplanet habitability in Section 5. We conclude with a short summary and an outlook in Section 6.

2 DATA

We started our systematic search for long-duration superflares on fully convective stars with the entire TESS archive. First, we used a colour and brightness cut to restrict our search to fully convective nearby stars (Section 2.1). Second, we searched their 2-min cadence light curves for flares, and selected those candidate events where we could detect a rotation period that was shorter than the visible flare decay (Section 2.2). Our final sample consists of four giant flares on four stars with spectral types M5-M7 and rotation periods between 2.7 and 8.4 h. Two stars in our sample lacked measurements of projected rotation velocity in the literature. For those we obtained high-resolution spectra (Section 2.3), so that we could determine the properties of all stars, including the inclination of the stellar rotation axis (Section 2.4).

2.1 Stellar sample and TESS light curves

We compiled a sample of nearby fully convective dwarf stars of spectral type M5 and later that were observed with TESS in the first two years of observations (Cycles 1 and 2). We selected bright stars with TESS magnitude <15, which we crossmatched within 5 arcsec with nearby stars (parallax >20 arcsec) in Gaia (Gaia Collaboration 2016) DR2 (Gaia Collaboration 2018). We required good G and RP colour measurements, i.e. band flux divided by its error (flux_over_error in the Gaia catalog) >10, and applied a colour cut to select for fully convective stars. We applied 1.27 < G − RP < 1.75, based on the relation between G − RP band colour and spectral type (Kiman et al. 2019) derived from spectroscopically identified M and L dwarfs in the Sloan Digital Sky Survey (York et al. 2000) (SDSS). In total, we analysed 3175 targets that were observed with TESS.

2.2 Flare sample

We searched for stars that displayed both long duration flares in their light curves, and clear rotational modulation outside of flares with a period shorter than the visible flare decay. Modulation of the flare itself was not a selection criterion. The typical TESS light curve spans one Sector, or 27 days, and we used only stars with 2-min cadence data to ensure that we could fit the shapes of the superflares in sufficient detail. Our sample covers TESS Cycles 1 and 2, Sectors 1-26, encompassing both hemispheres above and below the ecliptic, that is, about 73 per cent of the whole sky (Ricker et al. 2015). We used the flare-finding algorithm AltaiPony (Ilin 2021) to detect flares on 6291 light curves that were obtained for the 3175 sample stars.

From the 17 646 automated flare detections, we selected those with relative amplitudes above 20 per cent. For each light curve with at least one such flare, we derived the Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) using lightkurve (Lightkurve Collaboration et al. 2018) in order to determine the stellar rotation period as the signal with the highest peak and SNR > 2.5. For the promising cases, we compared the periodograms derived from the TESS light curves that were systematicaly corrected with the Presearch Data Conditioning pipeline (PDCSAP_FLUX, Jenkins et al. 2016), and those extracted from the raw flux (SAP_FLUX). We selected candidates where the measured flare duration was at least 15 per cent of the detected rotation period. The duration was chosen to ensure that we did not exclude the superflares that were split into multiple detections by AltaiPony due to rotational modulation. We arrived at this threshold by injecting model superflare events into simulated light curves with typical Gaussian noise levels, and recovering their properties using AltaiPony (see also Section 5.1).

2.2.1 Manual inspection of flare events

We identified 343 candidate events for superflares which we inspected visually. Out of those, in 19 cases the detected rotation period was a fluke in the periodogram where multiple peaks could not be separated into rotational and non-rotational signal. Some candidates were instrumental false positives (19), others were contaminations from Solar system objects that mimicked flare events (6), both of which could be disambiguated from the Target Pixel Files from which the light curves were extracted. The majority, however, were flares with total duration below one rotational period (291), partly because aliases in the periodogram suggested faster rotation, but mostly because the conservative 15 per cent relative flare duration threshold allows selection of short flares. Finally, four flares lasted for multiple rotation periods but exhibited a complex flare morphology that indicated the presence of more than two superimposed flares in each event (Table 1).

We observed four long duration superflares, i.e. flares with TESS band energies above 10³³ erg, on fully convective M dwarf stars in the red-optical regime (‘white-light’) with TESS: TIC 277539431 (2MASS J10551532-7356091), TIC 452922110 (2MASS J08055713 + 0417035), TIC 44984200 (2MASS J08380224-5855583), and TIC 237880881 (2MASS J01180670-6258591), hereafter referred to as TIC 277, TIC 452, TIC 449, and TIC 237.

2.3 High-resolution spectra

We collected high-resolution spectra for two of our sample stars that did not have literature values for their projected rotation velocity vsin i available. TIC 449 and TIC 277 were observed with the high-resolution spectrograph (HRS, Crause et al. 2014) at the Southern African Large Telescope (SALT, Buckley, Swart & Meiring 2006) on February 8 and 2020 February 9, respectively. The HRS provides spectra covering a wavelength range of |$5500\!-\!8900\,$|Å in its red arm with a spectral resolution of about |$40\, 000$| in our chosen medium-resolution mode. Since our targets have low effective temperatures, the bulk of the stellar emission was observed redwards of |$7000\,$|Å. The data were reduced with the PEPSI data reduction software (Strassmeier, Ilyin & Steffen 2018). The reduction followed the standard steps of bias overscan, detection, and subtraction, scattered light extraction from the inter-order space and subtraction, definition of échelle orders, optimal extraction of spectral orders, wave-length calibration, and a self-consistent continuum fit to the full two-dimensional (2D) image of extracted orders. The one-dimensional spectra resulting from combining the échelle orders were used to infer the rotational line broadening as described in Section 2.4.5.

2.4 Stellar properties

Our goal was to simultaneously constrain the latitude (θ_f) of the superflare on the stellar surface together with its longitude, and underlying shape. In order to infer θ_f from the flare light curves we required reliable information about the stars on which they occurred (Table 2). Rotation periods (P_rot) enter both the stellar inclination calculation, and the flare modulation model itself. Spectral types are needed to determine the luminosity of the stars, and from these the flare energies (E_f) and flaring region sizes (ω). Combining P_rot with stellar radius (R_*), and projected rotational velocity vsin i, we could calculate the inclination i of the stellar rotation axes. Eventually, knowing i allows us to break the partial degeneracy between i and θ_f in the model fit.

2.4.1 Rotation periods

We calculated the uncertainty on P_rot by combining the results from sine fits to all available TESS light curves for each target assuming Gaussian uncertainties on the flux, and derived the posterior distributions for P_rot using the Python Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013). The modulation was stable (Fig. 1), i.e. it did not show any sign of evolution in the phase-folded light curves of 1-3 consequent Sectors for each star. The resulting values for P_rot were then used in the inclination calculation, and the modulation model fits.

2.4.2 Spectral types

We adopted spectral types from previous work (Kraus et al. 2014; Kirkpatrick et al. 2016; Phan-Bao et al. 2017) that used low-resolution optical or IR spectroscopy for TIC 237, TIC 277, and TIC 449. For TIC 452, we compared a spectrum obtained from the FAST spectrogragh at the Fred Lawrence Whipple Observatory to the Bochanski et al. (2007) composite spectra. For TIC 277 and TIC 449, spectral types were also consistent with the literature values when comparing the SALT spectra to those in Bochanski et al. (2007). We selected the nearest stars within 5 arcsec of our targets in the Two Micron All-Sky Survey (2MASS, Skrutskie et al. 2006) and Gaia DR2 databases in order to obtain optical and infrared photometry and distances from parallaxes.

2.4.3 Stellar quiescent flux

The absolute released flare flux determines the modelled size of the flaring region in the photosphere. However, the light curves provided by TESS did not have accurately calibrated absolute fluxes. We calculated absolute quiescent fluxes in the TESS band using stellar spectrophotometry, transforming the flux in Gaia G band to the TESS band. We first constructed representative spectra for each spectral type that covered the TESS-Gaia range (⁠|$\sim 3500-11\, 000\,$|Å) by combining spectroscopic templates (Bochanski et al. 2007; Schmidt et al. 2014a) with the IR prism spectra from the Spex Prism Library (Burgasser et al. 2004; Burgasser & McElwain 2006; Burgasser, Cruz & Kirkpatrick 2007; Burgasser et al. 2008; Kirkpatrick et al. 2010; Burgasser 2014). We then integrated over the Gaia G band filter response curve, and normalized the spectral energy distribution (SED) to the Gaia value. To obtain fluxes in TESS, we integrated over the response curve for TESS. The major uncertainty in the spectrophotometric fluxes was in the spectral type used to select the SED, so we assigned uncertainties to the fluxes with a difference of one spectral subclass based on typical differences in the literature, which is higher than uncertainties quoted within individual works. We then converted flux to luminosity using the Gaia distances.

2.4.4 Stellar radii

We calculated stellar radii based on the relationship between absolute K_S magnitude and radius for M dwarfs (Mann et al. 2015). The absolute K_S magnitude was calculated from the apparent K_S magnitude in 2MASS and the geometric Gaia DR2 distance of the stars (Bailer-Jones et al. 2018). The uncertainties in the reported radii are based on a combination of the scatter in the magnitude–radius relationship (∼3 per cent) and the uncertainties on distance and K_S magnitude.

Systematic errors in the model fits may arise in the R_* calculation, which would affect the derived inclination, and relative size of the flaring region. Radii that are inflated by ∼12 − 14 per cent compared to stellar structure models have been reported for both components of the rapidly rotating late M dwarf equal mass wide binary GJ 65 (Kervella et al. 2016) which they attributed to youth and the suppression of convection due to strong magnetic fields. However, we determine the radii of our targets from absolute K_S magnitudes, with their distance being well known from Gaia DR2 data (Mann et al. 2015). This method is self-correcting for inflation because an M dwarf that is abnormally large will also have an absolute brightness that will be larger than for a non-inflated M dwarf.

2.4.5 Projected rotation velocity vsin i

The inclination i of the stellar rotation axis is partially degenerate with the flare latitude θ_f. This ambiguity needs to be broken by an independent measurement of i, which we achieve by measuring the projected stellar rotation velocity vsin i.

For two stars (TIC 452 and TIC 237) we used existing vsin i values (Kesseli et al. 2018; Kraus et al. 2014). The other two, TIC 449 and TIC 277, did not have existing vsin i measurements in the literature. For those, we used the collected SALT spectra (see Section 2.3) to determine the projected rotation velocity.

For the vsin i inference we used a technique that employs cross correlation functions (CCFs) of rotationally broadened template spectra to calibrate a relation between the FWHM of the CCF and vsin i (Reiners, Joshi & Goldman 2012). We calibrated the normalized CCFs using high-resolution optical CARMENES (Calar Alto high-Resolution search for M dwarfs with Exo-earths with Near-infrared and optical Echelle Spectrographs, Quirrenbach et al. 2010) spectra of CN Leo (M6), LP 731-058 (M6.5), and Teegarden’s Star (M7). This was done to reduce bias in the vsin i values derived with individual templates, and to account for uncertainty in the spectral types of the SALT targets. The template stars rotate at vsin i < 3 km s⁻¹ (Reiners et al. 2018). We broadened these three template spectra using the eniric (Neal & Figueira 2019) software assuming a limb darkening coefficient ϵ = 0.6 (Claret 1998) in the range from 1 to 60 km s⁻¹ taking 1 km s⁻¹ steps. Then we measured the FWHM of the CCFs to derive a calibration function for vsin i from the broadened templates.

In the selection of spectral regions, we started with suggestions for suitable wavelengths for low-mass stars in Reiners et al. (2018) in the 7410–8700 Å range that overlapped with our SALT spectra. We found the region 7938–7955 Å surrounding the Rb 7948 Å line to be best suited for the chosen technique, and the spectra collected for our target stars (Fig. 2). This line is strong, shows isolated line wings, and is largely unaffected by telluric features. These properties make it sensitive to changes in vsin i in the 10−50 km s⁻¹ range, while in the other inspected regions the lines were blended, or much weaker, or both. We fit a third order polynomial to the calibration function, and used the fit to propagate the uncertainty on FWHM to the resulting vsin i values (Fig. 3). We found uncertainties on vsin i in the 1.2−1.9 km s⁻¹ range for individual pairs of templates and spectra, and no apparent trends with the spectral type of the templates, so that we treated the results as independent measurements. Averaging over the best-fit results for the pairs we arrived at uncertainties of 0.8 km s⁻¹ for TIC 449 and 1.0 km s⁻¹ for TIC 277 with Gaussian error propagation.

2.4.6 Stellar inclinations

We used the stellar rotation period P_rot and radius R_* along with the vsin i to infer the stellar inclination angle i. To account for the statistical dependence between equatorial velocity (determined from P_rot and R_*) and vsin i, we used the framework of Masuda & Winn (2020). We used an MCMC sampler (Goodman & Weare 2010) to produce a posterior probability density for cos i given the measurements. We required 0 < cos i < 1, and assumed Gaussian priors on vsin i, P_rot and R_*. The posteriors for i were non-Gaussian for the stars determined to be close to equator-on, but could be approximated by a double Gaussian fit which we subsequently used as an informative prior in the model fits (Fig. 4). We give the median and 68 per cent confidence interval as the best-fitting values and errors in Table 2, using the convention that i < 90°.

3 METHODS

The presence (or absence) of a periodic modulation in the flare flux observed in the optical light curves in our sample can be modeled by a bright flaring region on the stellar surface that periodically, partially or fully, rotates in and out of view (Section 3.1, Fig. 5). We discuss the simplifications made in the design of the model (Section 3.2), such as assuming a |$10\, 000\,$|K flare temperature, a spherical cap as the shape of the flaring region, and a universal flare time evolution parametrization. We also motivate why we neglect differential rotation in our analysis, and why we assume that multi-peaked flares are co-spatial. Our assumptions do not pose considerable limitations on our conclusions. They may, however, partly explain the cases where the model could not capture all features of the flare light curve, such as deviations from the flare template near the peak.

3.1 The flare modulation model

Our (3+3+1)-parameter model was fit using an MCMC method (Fig. 6). The first three parameters are the flare peak time t_f, the corresponding longitude ϕ_f, and the latitude θ_f, which define the timing of the flare and its location on the stellar surface. The other three are the relative flare amplitude A, FWHM_i, and FWHM_g, which determine the intrinsic shape of the light curve. The last parameter is the inclination i of the stellar rotation axis, which is constrained by independent observations. The geometric modulation is determined by i, θ_f and ϕ_f that we combined to give the distance of the flaring region from the limb. The model light curve was constructed in three steps:

• Determine the size of the flaring region from the flare temperature T_f, and A, given the stellar radius R_* and luminosity L_*.

• Split the flaring region into spatial elements (θ, ϕ) to calculate the spatially resolved rotational modulation of the flare flux.

• Apply the modulation to the flare light curve and sum the spatial elements to obtain the model luminosity L_model.

3.1.1 Size of the flaring region

3.1.2 Rotational modulation and model luminosity

The flaring region with radius ω/2 centred on (θ_f, ϕ_f) was represented by an ensemble of N evenly distributed spatial elements. It was assumed to emit uniformly, so each spatial element emitted the same fraction F_f(t) = L_f(t)/N of L_f(t) at any given time t during the flare. Before and after the flare the region emitted the average quiescent stellar flux. To model the maximum flare flux F_f(t), we used the adapted D14 model.

For ease of calculation, we converted t to phase |$\hat{t}$| in units of radian using the rotation rate of the star |$\hat{t} = 2\pi (t - t_{0})/P_{\text{rot}}$| where t₀ marked the start time of the light curve that we selected for the model fit, and P_rot was the rotation period of the star.

3.1.3 Fitting the model to the observed light curve

First, we masked the flare and fit the dominating starspot modulation and long-term trends in the light curve using Gaussian Process (GP, Rasmussen & Williams 2006) regression with celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018). celerite2 is a software that provides a library of covariance functions, or kernels, to use for GP regression on one-dimensional data sets like time series, and more specifically, light curves. Following Foreman-Mackey et al. (2017), we chose a mixture of two kernels, called SHOTerm, which are derived starting from the dynamics of stochastically driven damped simple harmonic oscillators (SHOs) as a general model for stellar variability. The first is a RotationTerm, which is itself a mixture of two SHOTerm kernels with periods P_rot and P_rot/2, and used to model the rotational variability. The second is a single SHOTerm that covers all other variability that may be present in the light curve. We obtained posterior distributions on the parameters using the same wide Gaussian priors (σ = 2) for all fits except for the prior on the rotation period in the RotationTerm kernel. There, we used the period P_rot obtained from the Lomb-Scargle periodogram as a prior (see Sections 2.4.1 and 2.2). We then subtracted the prediction for the light curve in the time interval of the masked flare to obtain the residual flare light curve to which we subsequently fit our model.

In a second step, we constructed a likelihood function assuming that the uncertainty in flux (quadratic addition of uncertainty from GP fit and PDCSAP_FLUX_ERR, termed ‘original uncertainty’ from here on) was Gaussian. Our priors were uniform over the physically permissible ranges (>0 for the flare parameters, [0°, 360°] for the longitude of the flare, and [ − 90°, 90°] for θ_f) for all parameters but i, for which we adopted empirical priors (see Section 2.4.6). Using the MCMC method (emcee, Foreman-Mackey et al. 2013) we determined the posterior distributions of the flare peak time t_f; the intrinsic flare properties FWHM_i, FWHM_g, and A; the location of the flaring region, θ_f and ϕ_f; and an updated distribution of the inclination i.

3.2 Model assumptions and their motivation

3.2.1 Flare emission in the optical is a 10 000 K blackbody

We assumed that the observed red-optical emission of our flares stems from a blackbody component at effective temperature |$T_f\approx 10\, 000$| K, as is common in both the Sun and other low-mass stars (Hawley & Fisher 1992; Kretzschmar 2011; Kowalski et al. 2013; Kowalski, Mathioudakis & Hawley 2018). Hotter temperatures have occasionally been reported for M dwarf superflares in the past (Robinson et al. 2005; Froning et al. 2019). Typically, the short, impulsive phase can be hotter than |$\sim 10\, 000$| K, but cools rapidly to |$\lessapprox 10\, 000$| K where it emits during the longer, gradual phase (Fuhrmeister et al. 2008; Howard et al. 2020). In our model, T_f primarily affects the size of the flaring region. A large region will be more smoothly modulated than a small bright kernel. If an initially hotter temperature with subsequent rapid cooling was present, our model would overestimate the size of the flaring region in the first few min of the flare decay, but the protracted gradual phase that characterizes long-duration superflares would be largely unaffected.

3.2.2 The flaring region is a spherical cap

We also assumed that the flaring region was a simple spherical cap on the stellar surface. White-light flares on the Sun take more complex shapes in general, such as pairs of non-circular footpoints, or crescent-shaped ribbons (Hudson, Wolfson & Metcalf 2006). The morphologies of active regions that are associated with the strongest flares are also the most complex (Sammis, Tang & Zirin 2000). In particular, non-circular shapes may cause departures from our model when the region moves across the stellar limb, but less so when the region is fully in view. Moreover, effects like saturation of flare loops might violate the assumption that the flaring region is flat (Heinzel & Shibata 2018; Jejčič, Kleint & Heinzel 2018). However, unlike in solar observations, we can neglect limb darkening effects (see Appendix  B).

3.2.3 Flares follow the same time evolution parametrization

Flares often display a complex time evolution that does not follow the D14 empirical template. The majority of the largest flares are complex. Based on the ∼6800 flaring events that Günther et al. (2020) detected in TESS M dwarf light curves, out of 10 per cent of the flares with the highest amplitudes (A > 0.24) more than 52 per cent were best fit with superpostions of two or more D14 flare templates, while in the other 90 per cent this fraction was only ∼20 per cent.

Some flares cannot be resolved into superpositions of multiple classical flares (‘unusual’ flares). While we expect about one out of four flares to be complex (Hawley et al. 2014; Davenport et al. 2014; Günther et al. 2020), unusual flares are rare with about 1 per cent incidence rate (Davenport et al. 2014).

The D14 model was derived from the flares found in Kepler light curves of GJ 1243, a frequently flaring dM4e star. We based our flare model on the empirical parametrization derived for this star. However, if the time evolution of the flare and the geometry of the flaring region deviate from our assumptions, we expect this to mostly affect the impulsive phase of the flare. To avoid overfitting the gradual decay phase of the flare which reveals more about the rotational modulation than the impulsive phase due to its relatively long duration, we made two adjustments. First, we chose to decouple the time scales for the impulsive and gradual phases. This is consistent with models that attribute different physical processes to the impulsive and gradual phases of flares (Benz & Güdel 2010). Second, we increased the uncertainties on the flux above A/2 in each of the flares to 20 per cent of the flare-only flux, as shown in the bottom panels in Fig. 6. The number was an ad hoc choice aiming to permit relatively strong variability in the amplitude while preserving the overall shape of the flare.

3.2.4 Differential rotation effects are negligible

Differential rotation may in principle cause systematic errors in the inclination estimates. Assuming the predicted differential rotation (Ω_equator − Ω_pole)/Ω_bulk ∼ 0.02 (Yadav et al. 2015) for fully convective stars, derived inclinations would increase by about 0.5°−5°. Assuming complete degeneracy between i and θ_f, the flare latitudes would move poleward accordingly. The poleward shift would be at the high end of the estimate for flare latitutes that were initially close to the pole, and at the low end for flares near the equator. However, fully convective stars do not appear to reach such strong differential rotation, instead observations indicate typical values about one order of magnitude lower (Barnes et al. 2017). Moreover, the degeneracy is only partial, because one can, for instance, distinguish a pole-on from an equator-on flare for a relatively broad range of inclinations, but less so for smaller differences in latitude (compare, for instance, the left- and rightmost columns in Fig. 9). Overall, as long as differential rotation is Sun-like, i.e. rotation is slower at the pole than at the equator, the real flare latitudes will be a bit closer to the poles than what we derive in this work.

Another effect of differential rotation would cause the rotation period derived from the light curve to differ from the rotation period at the latitude of the flare. The input period would be either too long or too short, resulting in a poor fit if an otherwise strong rotational modulation of the flare light curve is present. However, we observed no lag between between the minima of modulation and the spot induced variability.

3.2.5 Flare footpoints are co-spatial with the loci of magnetic field emergence

One caveat concerns the interpretation of the flare latitude as the locus of magnetic field emergence (see Section 1, and references therein): We cannot exclude that the region of white-light emission footpoints could have wandered across the stellar surface in the time between emergence and flare eruption. The distance between where the magnetic field concentration is at the time of flux emergence and at the time of white-light emission will depend on the duration of the build-up phase of the flare, and the conditions of the sub-surface plasma. Our stars show stable starspot configurations over timescales of multiple weeks to months which suggests that large scale spot evolution is very slow, but we cannot exclude more dynamic conditions on smaller scales.

4 RESULTS

Polar spots have been observed in fast-rotating young Sun-like stars (Ișık et al. 2018). However, on the present-day Sun, sunspots, and the flares they are associated with, occur in a belt around the equator, and almost never at latitudes above 30°. We found all our selected flares at significantly higher latitudes around and above 55°.

We found four candidates with long-duration superflares, and shorter measurable rotation period P_rot: TIC 277 (Sector 12), TIC 449 (Sector 10), TIC 452 (Sector 7), and TIC 237 (Sector 1). They rotate so quickly (P_rot = 2.7 − 8.4 h) that the observed white-light superflares last longer than a full rotation period on each star. Three of the four stars display strong flux modulation during the superflare decay in phase with the stellar rotation. This can most straightforwardly be explained by rotational modulation, with the flaring region rotating in and out of view on the stellar surface.

In Fig. 6 we show all investigated superflare light curves. The shape of the modulation contains information on the latitude of the flaring region: the combination of flare latitude and the inclination angle of the stellar rotation axis determine how long a flaring region is visible as the star rotates. Similarly, unless the star is viewed strictly pole-on, the absence of flare modulation also contains latitude information, indicating that the flare occurs close enough to the rotational pole for the modulation to be indistinguishable from the noise, and that the pole is visible all the time; this is exploited for TIC 237 which displays much more subtle but still measurable rotational flare modulation.

We found that all four stars displayed their giant flares at high latitudes ≥55°. In Table 2, we give the results for the marginalized posterior distributions of the model parameters, and in Fig. 7, we illustrate the inferred active regions loci at flare peak time on the stellar surface. The full posterior distributions can be viewed in detail in Appendix  C.

4.1 Model setups: single and double-flare, original and increased uncertainty around the flare peaks

The stars TIC 277 and TIC 452 both showed single modulated flares. The flares we detected on TIC 449 and TIC 237 were complex because they were better fit with a superposition of two flares than with a single flare (see Section 3.2.3).

TIC 237 displayed two flares in close sequence, both without apparent rotational modulation. We fitted this light curve under the assumption that both flares occurred at the same location, which is common on the Sun (Uchida 1968; Wheatland & Craig 2006; Török et al. 2011) and proposed for other stars as well (Vida et al. 2016; Gizis et al. 2017). TIC 449 displayed excess flux indicating a secondary flare during the second rotation of the main event that was, however, less clear than in the case of TIC 237. A model with a single underlying flare for TIC 449 is unable to reproduce the shape of the light curve through rotational modulation, in particular around the secondary peak. A model with two co-spatial flare events reproduces the light curve better.

Since flares often display complex evolution in the impulsive phase and the flare peak, we accomodated possible deviations from the flare template by increasing the uncertainty to 20 per cent of the excess flux around the peak (above A/2, see Section 3.2.3) for all flares except for TIC 277 where the flare peak most likely occurred behind the limb, and is therefore not observed (bottom panels in Fig. 6). This is motivated by the fact that our model is most sensitive to the (absence of) rotational modulation in the gradual decay phase. Deviations in the short impulsive phase, which our model cannot capture, however, can cause over-fitting of the template because they are more likely intrinsic to the flare than to effects of geometrical foreshortening.

We could improve or achieve consistent modelling of the modulated decay phase of the flare in all cases when we increased the uncertainties. In Appendix  A, we show the alternative fits with original uncertainties from the TESS light curves (Figs A1, A2 (a) and (b), A3). While we preferred the fits shown in Fig. 6, the inferred latitudes were overall high even in the solutions with higher amplitudes and stronger systematic deviations in residual flux (>48°, see Table A1).

4.2 Star by star: model fits, alternative and disfavoured setups

4.2.1 TIC 277

TIC 277 (WISEA J105515.71-735611.3) is an M7 dwarf 13.7 pc away that rotates with a 4.56-h period and is seen nearly equator-on. The superflare occurred at |$\theta _f=\left(80.9^{+0.5}_{-0.6}\right)^\circ$| so that the flaring region was only briefly hidden from view during each rotation.

The flare peak occurred while the region was facing away from the observer (Fig. 6 (a)) at phase ϕ_f = 0.64 ± 0.01, and was a clear single flare event, so that we did not fit the model with alternative setups. Note that the center of the flaring region becomes visible around ϕ ∼ 0.70, while other parts of it stay behind the limb longer (at lower latitudes) or shorter (at higher latitudes).

The uncertainty on θ_f is small despite the relatively high uncertainty on inclination (Fig. 4). This is due to the sensitivity of the modulation model to flares that occur simultaneously close to the observed limb and to the rotational pole (see posterior probability distributions in Fig. C1). For the flare in TIC 277, the flare light curve morphology alone constrains the latitude to very high values even without knowing the inclination.

One possible caveat is the flare temperature. If the flare is visible throughout almost the entire rotation, which is the case for TIC 277, i can be considered a lower limit on θ_f, unless the flaring region is very large. The flaring region size depends on the flare temperature, which we here assumed to be |$10\, 000$| K. If T_f were only 6500 K, the observed flare could be placed at about 70° latitude. However, this temperature would be more typical for less energetic solar flares (see Appendix in Castellanos Durán & Kleint 2020, and references therein) than for stellar superflares, which are usually hotter (see Section 3.2.1). Similarly, a higher temperature could place the flaring region even closer to the poles.

4.2.2 TIC 452

TIC 452 (2MASS J08055713 + 0417035) is an M7 dwarf that has literature vsin i = (29.7 ± 1.5) km s⁻¹ (Kesseli et al. 2018), and a Gaia DR2 distance of 22 pc. A full rotation period of this star is 4.22 h, and the rotation axis is viewed at an inclination angle of about 49°. Its chromospheric activity places TIC 452 among active fast-rotating low-mass stars with saturated Hα luminosity (Newton et al. 2017). Both the single flare fit with original (Fig. A1) and the preferred fit with increased uncertainties above A/2 (Fig. 6 (b)) consistently implied a flare latitude of ≈63°.

In the remaining two stars, the simple model had to be adjusted to a more complex flare morphology.

4.2.3 TIC 449

TIC 449 (SCR J0838-5855) is an M6 dwarf 11.1 pc away rotating at P_rot = 2.71 h. The flare on TIC 449 was visible for almost five rotation periods. While the flux appeared modulated, the flaring region never fully rotated out of view. The light curve strongly suggests the presence of a secondary flare. However, the burst in the second maximum could also stem from other complex flaring behaviour in the wake of the large event. We therefore fitted both a single- and a co-spatial two-flare solution to the light curve with increased and original uncertainties above A/2 (see Section 4.1). The first and second events are indexed accordingly in Tables 2 and A1.

The fit with the least systematic residual flux in the important gradual decay phase was achieved in the two-flare fit with increased peak uncertainties (Fig. 6 (c)) which indicated that the flares occurred at high latitudes at (72 ± 1)°. Because of the contamination by complex behaviour in the secondary maximum the results of the fits were inconsistent with each other, but yielded high latitudes regardless: (86 ± 1)° for the two-flare fit with original uncertainties (Fig. A2 (a)), (83 ± 1)° for the single flare fit with original uncertainties (Fig. A2 (b)), and (70 ± 1)° for the single flare with increased uncertainties (Fig. A2 (c)). While a secondary flare is a likely scenario, other conceivable explanations include a flare region shape that deviates substantially from our circular simplification, or a cluster of >2 flares.

4.2.4 TIC 237

Finally, TIC 237 (2MASS J01180670-6258591) was our earliest type star, an M5 dwarf at a distance of 46 pc with literature vsin i = (14.4 ± 2.6) km s⁻¹ (Kraus et al. 2014). TIC 237 is a probable member of the Tucana-Horologium association (Ujjwal et al. 2020), a young nearby moving group with age estimates spanning 4−45 Myr (Ujjwal et al. 2020; Bell, Mamajek & Naylor 2015; Kraus et al. 2014).

The systematic deviation of the fit with original uncertainties, which implied θ_f ∼ 48°, was stronger in the decay phase (Fig. A3) than that with increased uncertainties on both peaks (Fig. 6 (d)), which suggests a flaring region latitude of |$\left(55^{+5}_{-6}\right)^\circ$|⁠.

In TIC 237, the flare modulation was very subtle yet measurable, in particular, in the decay phase of the second event. We show a close-up in Fig. 8, in which we placed the flaring region at 25° higher and lower latitutes relative to the best-fit solution, respectively. The observed modulation appears as an intermediate between stonger modulation at lower latitudes, and vanishing modulation closer to the pole. We cannot rule out that at least a fraction of the observed modulation is intrinsic to the flare. If all variation was intrinsic, the flare would reside on the pole. In a less extreme case, if there was a mixture of intrinsic variation and rotational modulation that gives the appearance of more rotational modulation than is actually present, the derived latitude would still be biased to lower values. The opposite case, where rotational modulation is masked by intrinsic variation, is less likely to occur: While variation that mimicks rotational modulation can manifest as modulation with many different flare latitudes and longitudes, a variation that masks rotational modulation has much fewer degrees of freedom. We therefore suggest to give more credence to the lower error quoted for TIC 237 than the upper error. The same reasoning applies, albeit to a lesser degree, to the other three flares.

5 DISCUSSION

The consistently high latitudes of the superflares studied in this work are significant, even though the sample size is small. If superflares were uniformly distributed across the stellar surface, the probability of finding all four superflares above 55° would be about 0.1 per cent, i.e. very unlikely, provided that our detection method is not biased towards finding preferentially high latitude flares (Section 5.1). If flares occurred at all latitudes θ_f equally, and all inclinations i were equally likely to be observed, purely geometric considerations imply that as many flares would be observed in the equatorial strip [ − 26.4°, 26.4°] as in the complementing northern and southern polar caps taken together.

Superflares that last longer than a full rotation period of the star on fully convective M dwarfs are rare, with about one event per star per 120 yr of continuous observation based on our sample, and distinct from the quasiperiodic pulsations phenomenon (Section 5.2). Yet each event is an opportunity to probe the magnetic nature (Section 5.3) and the space weather conditions (Section 5.4) on these cool, ubiquitous stars.

5.1 Observing polar versus equatorial rotationally modulated superflares

To test if our selection mechanism disfavours equatorial flares, we analysed the recovery efficiency of our pipeline on a grid of simulated flares that we injected into light curves with typical noise properties. To produce the synthetic light curves, we used our model to produce flares on a uniform grid of latitudes θ_f and inclinations i. We fixed the stellar properties and flare FWHM values (FWHM_i ≡ FWHM_g), and varied the underlying flare amplitude A and peak time t_f to capture the effects of photometric noise levels and the visibility of the flare peak on the detection probability and recovered properties. We show an example of 25 flares from this grid of model flares that occur at different combinations of i and θ_f in Fig. 9. The majority of flares shown in this example would have been recovered, except for events that remain attenuated for the entire duration of the event because they occurred close the limb (red scatter). These can be both polar flares seen at high inclinations (Fig. 9 (a)), and equatorial flares seen at very low inclinations (Fig.  9 (y)). Overall, our automated flare detection method (Section 2.2) did not show any bias towards polar over equatorial flare locations.

5.2 Quasiperiodic pulsations are ruled out

Rotational modulation is the preferred scenario for our observations because the period of the modulation matches the rotation period of the star. A different type of flare modulation, magnetohydrodynamic oscillations of plasma following strong flares, has been seen on the Sun (Aschwanden et al. 1999; Nakariakov et al. 1999) and some other stars (López-Santiago et al. 2016; Pillitteri et al. 2014). Such quasiperiodic pulsations have been mainly observed in coronal X-ray emission or in chromospherically active individual spectral lines (Srivastava et al. 2008; Nakariakov & Ofman 2001; Mitra-Kraev et al. 2005) with periods ranging from a few and up to about 20 min (Inglis & Dennis 2012). White-light flare oscillations have been measured in a number of Sun-like and low-mass stars, but the detected periods were predominantly between 5 and 45 min, and none exceeded 100 min (Pugh et al. 2016), remaining shorter than the shortest modulation period we found. In contrast, the flare modulation we report here is of the order of several h, and is observed in white-light.

5.3 Stellar dynamos, magnetic field emergence, and polar spots

Large polar spots have long been suspected to exist on fully convective fast rotators, based on indirect methods such as (Zeeman) Doppler Imaging (DI or ZDI), which map spectral modulations to stellar surfaces (Morin et al. 2008b, 2010; Barnes et al. 2017; Kuzmychov, Berdyugina & Harrington 2017). Moreover, a number of simulations predict that strong magnetic fields emerge near the poles in fully convective fast rotators (Gastine et al. 2013; Yadav et al. 2015; Weber & Browning 2016).

The body of knowledge on magnetic field configurations of rapidly rotating, fully convective stars is still small, not least because the targets are faint and atomic lines that are typically employed for these methods in earlier type stars fade out and the spectrum grows dominated by molecular signatures. However, magnetic fields of up to 7 kG have been reported from fully convective fast rotators (Shulyak et al. 2017), far above the equipartition values suggested by solar-type dynamo models (Pevtsov et al. 2003). The average value is about 5 kG as derived from molecular spectropolarimetry in a sample of nine M1 to M7 dwarfs (Afram & Berdyugina 2019). Among the later M5-M7 dwarfs, cool spots appeared to occupy the majority of the stellar atmosphere, and at a greater depth than in the earlier types.

In ZDI maps of fully convective dwarfs with spectral types ≥ M4, a wide range of magnetic filling factors was observed (See et al. 2019). This can be interpreted as manifestation of a bimodal dynamo with either strong bipolar or weak multipolar fields (Morin et al. 2010). The strong bipolar branch is associated with fast rotation, a stable spot configuration, and little differential rotation (Gastine et al. 2013), suggesting that the stars in our sample belong on this branch, and rotate as near-solid bodies.

Studies on individual stellar systems indicate the presence of both branches. V374 Peg (M4V, P_rot = 0.446 d, Newton et al. 2016), showed no dominating polar spot, and weak contrast small spots on its surface at most (Morin et al. 2008a). GJ 65 A and B, two rapidly rotating fully convective dwarfs with similar properties, revealed very different spot configurations (Barnes et al. 2017). While GJ 65 A appeared to have spots in a band around 30−40°, but also larger polar spots at 50 − 85°, observations of GJ 65 B are more consistent with extended spots at intermediate latitudes (50−56°) and no spots above 70°. GJ 791.2A (M4.5V, P_rot = 0.257 d, Newton et al. 2016) showed spots at all latitudes, but circumpolar spots were larger than those closer to the equator. The spot structure of LP 944-20 (M9V) could be reconstructured with a dominating large spot at 79°, and marginal detections of smaller spots at lower latitudes (Barnes et al. 2015).

Finally, Berdyugina et al. (2017) reconstructed a magnetic field loop from giant flare-like bursts of the fast-rotating (2.84 h, Hallinan et al. 2008) M8.5 dwarf (Deshpande et al. 2012) LSR J1835 + 3259, one of which appeared to be modulated by rotation. Lacking a precise stellar rotation axis inclination, they could not derive the latitude of the loop’s footpoints. However, they recovered the geometry of the loops, and found that if the bursting region occupied ∼3 per cent of the stellar surface  (Hallinan et al. 2015) – comparable to the sizes of flaring regions in our sample – its magnetic field strength could reach 10 kG.

In this context, our results are further evidence that in fully convective, rapidly rotating dwarfs, strong magnetic field concentrations may indeed emerge at latitudes that are closer to the rotational pole than to the equator, where they create active regions that can erupt in giant flares.

5.4 Implications for stellar space weather

Our findings are relevant to exoplanets; several thousand of such planets outside of the Solar System have been discovered up to now (Perryman 2018). Planets residing in the habitable zones around their host stars, that is, orbital distances that allow for liquid water on the planetary surface, are of prime interest. The habitable zone is located close to the host star in the case of fully convective M dwarfs, enabling efficient detection of temperate planets, in comparison to more massive host stars. In order to characterize habitability of temperate planetary surfaces, the longevity and composition of exoplanetary atmospheres needs to be taken into account.

Numerical simulations have shown that flares and coronal mass ejections, which are directly associated with flares in the solar case (Webb & Howard 2012), can erode or alter exoplanetary atmospheres, and cause loss of exoplanetary oceans (Shields, Ballard & Johnson 2016; Tilley et al. 2019). Our results suggest that the largest flares on such fully convective stars preferably occur at high latitudes. This decreases the flare energy that is received at the planetary habitable zone orbit, but will less affect X-ray and (E)UV emission (Section 5.4.1). Coronal mass ejections are thought to accompany stellar superflares (Drake et al. 2013; Odert et al. 2020). However, in contrast to flares, they are only indirectly and rarely observed in stars to date (Moschou et al. 2019). If they can escape the strong magnetic field of the star (Alvarado-Gómez et al. 2018) and propagate at large angles to the planetary orbital plane they would pose a lesser risk to their planets’ environments than predicted under the assumption of a Sun-like low latitude distribution (Tilley et al. 2019, Section 2.4, and our Section 5.4.2). These considerations assume that the planetary orbital plane and the stellar equator are roughly aligned, which is not well understood for late-type systems (Section 5.4.3).

5.4.1 The impact of flares decreases with increasing flaring region latitude

With the present results, we have the opportunity to estimate the relative flux received by exoplanets, that are potentially hosted by these or similar stars, as a function of flare latitude. All else equal and assuming spin-orbit alignment, placing the flares in our sample at 45° or 80° latitude, the optical flux emitted in the direction of a distant observer in the orbital plane decreases by ∼30 and ∼85 per cent compared to the same event located near the equator. The exception is the first flare in the light curve of TIC 237 that, when placed at 80°, is visible for a larger fraction of the stellar rotation period due to the relatively large size of the active region. We show the typical relative attenuation of optical flux as a function of flare latitude in Fig. 10.

This comparison assumes that compact bright kernels in the stellar photosphere contribute the majority of the blackbody radiation to the flare events. This is true for most solar flares, but it is possible for superflares with particularly large flare region sizes (⪆2.5 per cent of the total stellar surface) to saturate the density in the flare loop arcades in the gradual decay phase such that they become optically thick and contribute significantly to, and may even become the dominant source of, white-light emission (Heinzel & Shibata 2018; Jejčič et al. 2018). Flare emission from a volume with considerable height, as it is proposed by Heinzel & Shibata (2018), instead of a flat slab on the stellar surface may explain the residual flux in our model fits in the impulsive phase because geometrical foreshortening cannot be applied to an extended 3D loop. Yet the good overall agreement between our model and data in the gradual phase suggests that the flat circular kernel is a satisfactory approximation for the purpose of accurate latitude inference.

However, the more energetic (E)UV and X-ray emission is produced at higher altitudes on the Sun. Spectropolarimetric observations of a large flare on a late M dwarf suggest a similar stratification for the stellar case (Berdyugina et al. 2017). Furthermore, stellar coronae are optically thin, meaning that the X-ray photons originating from them are not subject to the self-shadowing effects that optical emission experiences near the stellar limb. The X-ray flux received by the planets may therefore be much less attenuated than the optical continuum appears in Fig. 10. Ultimately, an accurate assessment of flare radiation impact on habitable zone exoplanets will have to include the combined effects of flaring region location, strong magnetic fields, high density, and large pressure gradients in the compact atmospheres of late type stars and brown dwarfs on the 3D emission distribution of superflares.

5.4.2 Coronal mass ejections and energetic particle events

The majority of X-class flares on the Sun are associated with coronal mass ejections (CMEs, Compagnino, Romano & Zuccarello 2017). A similar correlation is expected for stellar superflares (Odert et al. 2020), but direct CME observations are missing to date (Crosley & Osten 2018), perhaps due to differences in the spectral signature between solar and stellar CMEs (Alvarado-Gómez et al. 2019, 2020). Solar energetic particle (EP) events typically carry a few per cent of the energy of the CME they are associated with (Aschwanden et al. 2017). Extrapolating the power law relation between X-ray flux and solar EP fluence to stellar superflares is poorly constrained so that stellar events may contain one to five orders of magnitude more energy than the largest EP events on the Sun (Herbst et al. 2019). If stellar EP events associated with the most powerful flares are ejected from the star at high latitudes, as suggested by our observed flares, they may typically not collide with an exoplanet’s atmosphere, rendering their effect on habitability smaller than for equatorial flares. However, such particles coming from CMEs and EPs might also be deflected towards the equatorial plane by a strong stellar dipole field (Kay et al. 2019). The impact of particle events also hinges upon their ability to leave the stellar magnetosphere. In M dwarfs, there is evidence that CMEs and EPs may be suppressed by the large overarching dipole field that these stars often possess (Alvarado-Gómez et al. 2018; Fraschetti et al. 2019).

5.4.3 Key uncertainty: Spin-orbit alignment

All of the above assumes that exoplanetary orbital and their host stars’ rotation axes are aligned. While there is statistical evidence for spin-orbit alignment in Kepler exoplanet hosts (Mazeh et al. 2015), cases of misaligned systems suggest that alignment is less common than one would expect from protoplanetary disk formation (Winn & Fabrycky 2015; Bourrier et al. 2018; Louden et al. 2021; Hjorth et al. 2021). For small planets around M dwarfs, however, only few measurements exist, all of which indicate spin-orbit alignment (Hirano et al. 2020; Addison et al. 2020; Palle et al. 2020; Stefansson et al. 2020). Spin-orbit alignment in low-mass star systems with Earth-like planets is therefore not yet well understood, which ultimately limits our ability to assess the influence of flare emission direction on exoplanetary space weather.

6 SUMMARY AND CONCLUSIONS

In a systematic analysis of fully convective stars observed with TESS, we detected four stars that displayed giant flares which were modulated in brightness by the stars’ rapid rotation. The exceptional morphology of the modulation allowed us to directly localize these flares between 55° and 81° latitude on the stellar surface. Our findings are evidence that strong magnetic fields tend to emerge close to the rotational poles of fast-rotating fully convective stars, and suggest a reduced impact of these flares on exoplanet habitability.

Our search encompassed the first two years of TESS observations, during which one hemisphere per year was observed. By summer 2021, TESS will have completed its second full-sky scan that will give us the opportunity to put the conclusions of this work to test with more detections of long duration superflares.

The strength of our method lies in its complementarity to spectroscopy and spectropolarimetry not only with respect to the type of scientific information used to infer the stars’ magnetic properties, but also regarding the time scales probed. Flares trigger sudden and drastic changes in the local magnetic field configuration. Space-based time series photometry can capture its effects with a few min or less resolution. In contrast, ground based (Z)DI reconstructs larger scale changes on time scales of h and longer.

This work is a step towards understanding spatially constrained flare evolution on stars other than the Sun. Together with multiwavelength flare studies (Kowalski et al. 2013; Guarcello et al. 2019; Namekata et al. 2020; Maehara et al. 2020) that provide insight into the propagation of flare events through the stellar atmosphere we may soon see the first fully empirical spatio-temporal flare reconstructions on low mass stars. Such models would help us understand the so far elusive relation between stellar flares, coronal mass ejections, and energetic particle events, and better assess the impact of stellar magnetic activity on the space weather environment of exoplanet systems.

ACKNOWLEDGEMENTS

The authors would like to thank Uwe Wolter for his thorough reading and constructive criticism that markedly improved the quality of this work. EI acknowledges support from the German National Scholarship Foundation. KP acknowledges support from the German Leibniz Community under grant P67/2018. JRAD acknowledges support from the DIRAC Institute in the Department of Astronomy at the University of Washington. The DIRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences, and the Washington Research Foundation. This paper includes data collected with the TESS mission, obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This study was enabled by the TESS guest observer program G022164 led by JSP, ERN, JRAD, SJS, and Adam Burgasser. This research has benefitted from the SpeX Prism Library, maintained by Adam Burgasser. We are thankful to Andrew Vanderburg for providing the code used to calculate the inclination probability distributions.

Some of the observations reported in this paper were obtained with the Southern African Large Telescope (SALT). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Based on data from the CARMENES data archive at CAB (INTA-CSIC). This research also made use of the Python packages numpy, pandas astropy, and the cross-match service provided by CDS, Strasbourg.

DATA AVAILABILITY

TESS light curves are publicly available through the Mikulski Archive for Space Telescopes (https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html). The Gaia data used in this work are publicly available through the Gaia archive (https://gea.esac.esa.int/archive/). Data from the 2MASS Point Source Catalogs is accessible through the Centre de Donnees astronomiques de Strasbourg (CDS) via the VizieR Service (http://vizier.u-strasbg.fr/cgi-bin/VizieR?-source = II/246]). CARMENES spectra are publicly available through the CARMENES data archive (http://carmenes.cab.inta-csic.es/gto/jsp/reinersetal2018.jsp). SpeX Prism Library is available through http://www.browndwarfs.org/spexprism. Reduced SALT spectra ara available via Zenodo, under DOI:10.5281/zenodo.4332142. Low-resolution spectra obtained from the FAST spectrogragh at the Fred Lawrence Whipple Observatory (used for spectral typing of TIC 452) are available online at: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi%3A10.7910%2FDVN%2FL5WNU4&version=1.0. Python code including the flare modulation model, model fitting scripts, Gaussian Process regression scripts, scripts that were used to produce the figures and tables in this manuscript is publicly available through GitHub (https://github.com/ekaterinailin/MalachiteMountains). The code used to calculate the inclination probability distributions was provided by A. Vanderburg for use in this work and is available by request to the authors and with permission from A. Vanderburg.

REFERENCES

APPENDIX A: ALTERNATIVE AND DISFAVOURED FITS TO THE LIGHT CURVES

In our analysis we used different model setups to fit the light curves. The results we achieved using setups that we ruled out or disfavoured are summarized in Table A1, and illustrated in Figs A1–A3 for TIC 452, TIC 449, and TIC 237, respectively.

APPENDIX B: LIMB DARKENING

The geometric flare modulation model does not account for limb darkening effects that can affect the shape of the modulated flare light curve. The effect increases with flare region size: The larger the fraction of the stellar disk that covered by the flaring region the more do differences in observed stellar flux across the disk take effect.

In Fig. B1, we illustrate a scenario that represents the flare with the largest flaring region in our sample (first flare on TIC 237), which is comparable to some of the largest region sizes derived for flares on fully convective dwarfs in the literature (Schmidt et al. 2014b; Xin et al. 2021). The effect is negligible compared to the noise in the light curves in our sample.

APPENDIX C: MCMC MODEL FITS

In Figs C1–C4, we show the posterior distributions for the MCMC fits with the best-fit results quoted in Table 2. The upper left density plot in each figure shows the partial degeneracy between θ_f and i, which is best seen in TIC 277 (Fig. C1), where the empirical prior for i was non-Gaussian. In the distributions for TIC 277, TIC 449, and TIC 237, there are some correlations worth noting:

The flare parameters A₁, FWHM_{i, 1}, and FWHM_{g, 1} of the first event in TIC 449 (Fig. C3) are sensitive to the uncertainty on the flare peak time t_{f, 1}, but the flare latitude is not affected. The tail in the distribution of FWHM_{g, 2} suggests that the secondary flare could be fit with a single exponential in the decay phase.

In TIC 277, inclination and flare amplitude are correlated. Higher inclination implies stronger geometrical foreshortening. The distribution suggests that there are three particularly favorable solutions, depending on A, which seems counterintuitive at first. But since higher A implies a larger flaring region, we can interpret the solutions as combinations of different flare region sizes viewed at certain inclinations near the pole that result in comparable light curve morphologies.

Finally, in TIC 237, θ_f has a tail towards lower latitudes that correlates with higher flare amplitudes A₁ and A₂. However, as we can see in Fig. 8, rotational modulation in the decay phase of the event would become more and more pronounced with decreasing latitude, thereby breaking the partial degeneracy in A₁, A₂ and θ_f.