Searching For Gamma-ray Emission from Stellar Flares

Flares from magnetically active dwarf stars should produce relativistic particles capable of creating gamma-rays. So far, the only isolated main sequence star besides the Sun to have been detected in gamma-rays is TVLM 513-46546. Detecting gamma-ray flares from more dwarf stars can improve our understanding of their magnetospheric properties, and could also indicate a diminished likelihood of their planets' habitability. In this work, we stack data from the Fermi Gamma-ray Space Telescope during a large number of events identified from optical and X-ray flare surveys. We report an upper limit of gamma-ray emission from the population of flare stars. Stacking results towards control positions are consistent with a non-detection. We compare these results to observed Solar gamma-ray flares and against a model of emission from neutral pion decay. The upper limit is consistent with Solar flares when scaled to the flare energies and distances of the target stars. As with Solar flares, the neutral pion decay mechanism for gamma-ray production is also consistent with these results.


INTRODUCTION
The Sun is a -ray source both while quiescent (Abdo et al. 2011;Orlando & Strong 2008) and when it flares (Ajello et al. 2021;Ackermann et al. 2017;Ajello et al. 2014).The first Fermi-LAT Solar Flare Catalog (FLSF, Ajello et al. 2021) provides detailed insights into -ray Solar flares, such as their emission mechanism and associations with coronal mass ejections (CMEs).Extreme Solar events can severely impact the Earth's atmosphere and even create geomagnetic storms that impact human activities (Varela et al. 2022).Being more magnetically active than the Sun, red and ultra-cool dwarf stars can generate flares that are orders of magnitude more energetic than Solar flares.The extremely energetic stellar surface activity can lead to proton acceleration events such as CMEs.Stellar winds, CMEs and accompanying high energy photons can alter the atmosphere and magnetosphere of surrounding planets (Hazra et al. 2022;Chen et al. 2021).Some of these stellar flares can be categorised as superor even mega-flares as they reach output energies over 10 36 erg s −1 .Such events are predicted to emit -rays (Ohm & Hoischen 2018).These stars also have much higher flare frequencies compared to the Sun.With their large abundance in the Galaxy, and their prevalence for hosting planets, detecting -ray flare events from these stars would imply a greatly diminished number of habitable worlds.
Many recent ground-and space-based missions and facilities, searching for transient events at many wavelengths, are able to capture stellar flares and create large flare catalogs.These include ★ E-mail: yuzhesong@swin.edu.auTESS (Günther et al. 2020), Kepler (Hawley et al. 2014), Evryscope (Howard et al. 2019), MAXI (Tsuboi et al. 2016), the Australian Square Kilometer Array Pathfinder (Rigney et al. 2022), and the Deeper, Wider, Faster Program (Dobie et al. 2023).In Song & Paglione (2020), a 4 detection was achieved by phase-folding the lightcurve of TVLM 513-46546, an unusually active, nearby, and rapidly rotating radio dwarf star.However, a residual photon counts stacking search for -ray emission from 97 nearby flare stars did not result in a detection.With the advent of extensive flare surveys, we are now able to select a much larger sample with hundreds of stars and thousands of flares to search for stellar -rays.Further, development of -ray stacking techniques using joint likelihoods (Principe et al. 2021;Song et al. 2023;Paliya et al. 2020;Ajello et al. 2020;Paliya et al. 2019) indicates improved sensitivity in mining Fermi-LAT data for sub-threshold signals.In this work, we attempt a stacking method utilising joint likelihoods, combined with a windowing scheme that isolates the -ray data around the flare times, to search for -rays from 1505 flares from 200 flaring dwarf stars.

Sample Selection
We utilised optical flare data from the Evryflare (Howard et al. 2019) and TESS Flare (Günther et al. 2020) surveys.These optical surveys have high cadences of two minutes and cover most of the sky to provide accurate flare times for a very large number of uniformly distributed stars.For example, the Evryscope flare catalog observed 575 flares from 284 stars, and the TESS survey observed 8695 flares from 1228 stars.We also used the all-sky survey of X-ray flares by the Monitor of All-sky X-ray Image (MAXI) project (Tsuboi et al. 2016), which detected 21 giant flares from 13 active stars.These X-ray flares, while smaller in number compared to the optical flare surveys, are extremely energetic and may therefore be more likely contributors of -ray emission within the stacking survey.Only stars with Galactic latitude > 20 • were selected in this study to avoid the complicated -ray background caused by the Galactic plane.These surveys were able to identify, in many cases, multiple flares from any individual star during the observation (Table 1).Distance cuts were made so that the estimated flare fluxes could be within the range of fluxes that we can conceivably probe with the stacking methods (as high as ∼ 10 −11 cm −2 s −1 , detailed in § 3.2).For the Evryflare survey, stars within 50 pc were included, which resulted in 106 stars and 244 flares.For the TESS Flare survey, stars within 25 pc were included, which resulted in 86 stars and 1225 flares.For the MAXI flare survey, 8 stars satisfy the Galactic latitude cut, which included a total of 16 flares.No distance cut is applied to the MAXI flare stars given the low number of stars and extremely high flare energy.The 8 stars have distances ranging from 6.5 pc to 88 pc.

Fermi-LAT Data Analysis
The analysis in this work utilised the third revision of the PASS8 data, P8R3, released on Nov 26, 2018, of -LAT, along with the 10-year source catalog (the 4FGL-DR2, Abdollahi et al. 2020), hereafter the 4FGL.The latest Galactic interstellar emission model gll_iem_v07, and isotropic background model iso_P8R3_SOURCE_V3_v1 (Abdo et al. 2009) were used.Data analysis in this work was performed with Fermipy (Wood et al. 2017) version 1.0.1 1 based on the Fermi Science tools Anaconda distribution (Anaconda 2020) version 2.0.8 2 .
For each star, the analysis was performed both on the full mission elapsed time (MET), and within a specified time window (or windows) around the flare(s).The full MET analysis covers MJD = 54678.05to MJD = 59453.00.The window began just before the peak of the optical or X-ray flare, and extended for many hours beyond the peak, based on the prevalence of Solar -ray flares to have such long durations.For Evryflare and MAXI targets, we used Fermi-LAT data from 15 minutes before the start of each observed flare, and for 24 hours after the peak time of the flare.For TESS targets, however, in many cases there were a lot of recorded flares within the span of a day.All flares that happened within one day of one another were observed in the same window.A flare window started 15 minutes before the peak time of the first flare, and ended 24 hours after the peak time of the last flare identified in this window.If a target star had more than one flare window, then all the flare windows were combined using the tool gtselect.
Each ROI was a 21.2 • ×21.2 • square and centered at the location of a star, which corresponded to a ∼ 15 • radius region of interest (ROI).The data were also filtered using a zenith angle cut of 90 • to avoid bright emission from the Earth.Good time intervals were chosen with conditions DATA_QUAL==1 && LAT_CONFIG==1.In this analysis, the data were binned uniformly in 37 logarithmically-spaced energy bins, between 300 MeV and 100 GeV.The standard binned likelihood analysis process was performed on each ROI within the observation time periods described above, containing all the observed flare times of each star or the full MET.The star was added as a source with a power law (PL) spectrum to the center of the ROI model.To quantify the significance of any detection, the Test Statistic (TS) of each star was obtained through this process, defined as TS = 2 log L/L 0 , where L is the likelihood from the best-fit model containing the star, and L 0 is the likelihood of the model for the null hypothesis.

Stacking Analysis
The binned likelihood analysis of Fermi-LAT data returned a model of the ROI which was then used for the stacking analysis.A point source with a PL spectrum was placed at the location of the star, and the flux and spectral index were fixed.Only the background isotropic and Galactic diffuse normalisations were free to be fit by Fermipy's fit function.This process returned a log likelihood for the ROI given these spectral parameters for the central source.We repeated this process over a grid of flux and spectral index values.TS maps in flux-index parameter space for each ROI were made by subtracting the log likelihood of the grid point representing the null, which was chosen to be at the lowest flux in the grid (7.2 × 10 −12 cm −2 s −1 ) and softest spectral index of −4, then multiplying this result by 2. The stacked TS map was then made by summing the individual TS maps.
Test sources located at supposedly empty sky locations should be also analysed to serve as a set of controls to compare with the results of the stars.There are two ways of choosing the control comparison: first is to choose an "off" time window for each star with the same exposure and perform the same analysis as mentioned above; and second is to choose a random location within each ROI.Since we cannot guarantee that these relatively short time windows have the identical exposure in the "off" time, we perform the control comparison with the latter method.Test sources for this analysis were placed at a random location within the ROI five degrees away from the star, within the same flare window, and subject to the same analysis procedures described above.

Results
The individual source analysis was first applied to all 200 flare stars using data covering the entireMET.None of the stars have a TS value larger than 25, which traditionally indicates a significant detection.Stacking the fullMET data of these stars also returned results consistent with a non-detection.
When analysing the flare windows, no star had central sources with TS values larger than 25.An overwhelming number of the stars in fact had TS values ∼ 0. Overall, the TS distributions of both the stars and the test sources are similar to each other as well as to the theoretical null, which is proportional to a  2 distribution for two degrees of freedom (Fig 1).
We also report the stacked parameter space TS maps of the stars and the controls in the bottom two panels of Fig 1 .Since no individual target sources or control field test sources are detected, which is evident in the TS distributions in the top panel of Fig. 1, we only explore the spectral parameter space up to the point source detection sensitivity of the LAT.Any higher, and a source should be detected.Out of the 200 stars, 53 converged in the stacking analysis, containing a total of 298 flares.The spectral parameter stack for the flare stars, even though the total TS never reaches 25, peaks at a photon index of −2.5 +1 −1.5 .The monotonic increase in likelihood with flux is simply indicative of a flux upper limit.This result is consistent with the PL spectral index values observed for -ray solar flares (Ajello et al. 2021).The stacked TS map for the test sources, however, appears to be a non-detection with its maximum TS value in the corner of the parameter space at the highest flux and the softest spectral index.The difference in TS values of the test sources and flares, Δ(TS) = 7, also indicates that the existence of flare -ray emission is only marginally favored.To obtain the 90% confidence level upper limit of the flux of the stack, the flux of the stack is increased until , where L max is the maximum likelihood of the null hypothesis, and L max | 90 is the maximum likelihood of the model with increased flux, or Δ(TS) = 2.71 (Huber et al. 2012;Cowan 1997).Setting the PL index to be −2.5, we estimated the upper limit flux of the stacked flares to be 1.9 × 10 −10 cm −2 s −1 .Although this flux is about an order of magnitude below the point source sensitivity of the LAT, this analysis remarkably may constrain the PL spectral index of the stellar flares.
It bears repeating that the TS distributions for the flares and control field test sources are not only indistinguishable from each other, but also from the null.The null distribution is not an isolated peak at TS = 0, but the  2 , which means it stacks up to a nonzero value.We do indeed observe a non-zero, but insignificant, peak in TS for test sources in the bottom left panel of Fig. 1.We also note that the control field stack is spectrally distinguishable from the flares.This behaviour is known from other analyses of control field stacking (Song et al. 2023;Paliya et al. 2020) and implies that they trace or resemble the diffuse -ray background.

DISCUSSION
Compared to Song & Paglione (2020), the most obvious change made in this work is the much larger sample size.In Song & Paglione (2020), we examined flare stars that had been detected in radio and/or X-ray surveys, and utilized Fermi-LAT data from the entire MET.Flares from these stars might emit in -rays, but examining the full MET can dilute the signal and decrease the sensitivity.In this work, we choose optical and X-ray surveys that provide the time of every flare, which allows us to isolate each one individually and avoid signal dilution.More importantly, the analysis methods have been significantly updated and are more sensitive.Rather than stacking residual photon counts, the stacking analysis is now performed on the likelihood profile in spectral parameter space of each flare, which proves to be more sensitive.

Flare Frequency Distribution
Since the Sun is the only star with individual flares observed in -rays, it is our best template to understand any stacked flare signal in this study.To establish the appropriate context in order to compare stellar and Solar flares, we first examine the flare frequency distributions (FFDs) of all the flares investigated in this work.The FFD describes the rate of flares above a given energy , and typically follows a power law: where  and  are free parameters to be fitted.The power law index  is often used as an indication of the magnetic activity of the stars (Shakhovskaia 1989;Paudel et al. 2018).In examining the FFDs, we can illustrate the differences and similarities between the Solar and stellar flares, and justify the scaling of the Solar flares in the following analysis.We estimate the total flare energy with the GOES soft X-ray (SXR) observations of solar flares, catalogued by Plutino et al. (2023), who provide a detailed list of SXR flares between 1986 -20203 .The flare times in the SXR catalog are matched to the FLSF catalog.All SXR flares that fall between the estimated start and end time of a FLSF flare, are counted towards that FLSF flare.The summed integrated flux of all SXR flares within a FLSF flare, multiplied by 4AU 2 , is the total flare energy.The caveat of this estimation is that the total energy of a solar flare released in SXR only serves as a lower limit.A potentially more accurate estimate of flare energy is from the proton energy, given that the flare energy should be 20 times the proton kinetic energy (Ohm & Hoischen 2018).However, it is beyond the scope of this work to correlate SXR luminosity to total flare energy output.For the five FLSF flares that are studied in Aschwanden et al. ( 2017), we estimate their flare energies as 20 times the Solar energetic particle energy.
The FFDs of the all the flares in this study, as well as all Solar flares in FLFS, are produced using Altaipony4 (Ilin 2021;Ilin et al. 2019;Davenport 2016), in which  and  are estimated using an MCMC power law fitting method (Wheatland 2004).All the FFDs are presented in Fig 2 .Given the same flare energy, the stars in this survey flare far more frequently than the Sun, and their flare energies are also much higher compared even to -ray Solar flares.The -ray Solar flare FFD is generally very low and steep in comparison.
The vast majority of the flares in the sample can be thought of as extremely high energy versions of solar flares.While the flare energy varies, the power law slope  of the Solar flares FFD are comparable to those of the TESS and Evryscope flares, indicating they have similar physical origins.In contrast, the Solar and stellar flares have very different  values compared to the MAXI flares, which is expected as these X-ray megaflares are associated with young stars or RS CVn systems.Regardless, these X-ray megaflares are still not detected in the stack, and contribute to the significance of the stack as much as the rest of the sample.

Comparison With Solar Flares
Having justified the common origin and scalability of stellar and Solar flares in the previous subsection, we now estimate the expected stellar -ray signal based on an examination of the Solar flares in the FLSF observed between 300 MeV and 10 GeV by Fermi-LAT.Our stacking results place a sensitive upper limit on the average -ray flux from these flares and constrain the PL index of their -ray emission.
As a simple comparison, we first scaled the -ray fluxes of the 26 PL solar flares to a distance of 25 pc, the average distance of the 53 stars being stacked.Even the most energetic solar flare only has a -ray flux of 1.42 × 10 −16 cm −2 s −1 , which is many orders of magnitude below the range of parameter space being examined in this work.These scaling results are plotted as the black data points in Fig. 3.
Assuming that flare flux depends linearly on total flare energy (which we substantiate in the next subsection), we can further scale the -ray flux of the FLSF solar flares using the estimated Solar flare energy described in § 3.1, and where  stellar = 2.3 × 10 33 erg is the median flare energy of the sample stellar flares,  solar is the flare energy of any given Solar flare, and  = 25 pc is the average distance to the target stars.The range of these scaled Solar flare -ray fluxes is 1×10 −16 to 3×10 −11 cm −2 s −1 , which is below the LAT sensitivity limit of ∼ 10 −9 cm −2 s −1 , but overlaps with the fluxes probed by our stacking method.

Emission Modeling
In this section, we explore how our sensitive flux upper limit constrains the flare physics for these stars.We use Naima (Zabalza 2015;Kafexhiu et al. 2014) 5 , a Python package that computes radiation from non-thermal particle populations and also does MCMC fitting to spectra.Solar -ray flares appear to be well described by the decay of neutral pions which are created when non-thermal protons strike the Solar atmosphere.Proton populations with PL spectral indices Γ p ranging from −6 to −3.2 yield spectra consistent with Solar -ray flares (Ajello et al. 2021).We use input proton spectra with Γ p = −6 and −3, and three different values for the target density of the stellar atmosphere: 10 8 , 10 10 , and 10 12 cm −3 .These reflect the proton densities theoretically estimated for flare stars by Ohm & Hoischen (2018), and for the Solar disk model of Seckel et al. (1991), given the average depth for proton absorption.For any given combination of Γ p and atmospheric density, the proton spectrum normalisation can be determined from the total proton kinetic energy, which is 5% of the flare energy.
The model results for each proton density are plotted as the overlapping shaded areas in Fig 3 .The lower and upper boundaries of each area are for Γ p = −6 and −3, respectively.The -ray flux of Solar flares as a function of their SXR flare energy estimation, and the upper limit from the stacking results are also plotted.The pion decay models are consistent with both the Solar flares as well as the upper limit from the stellar flare stack.This further indicates that emission mechanism of stellar superflares is likely similar to that of Solar flares.Additionally, the results shown in Fig. 3 indicate that the energy conversion from flare to proton kinetic energy should not be more efficient than 5%.If more flare energy were converted to proton energy, creating the same level of -ray flux would require a less energetic flare, which would start to contradict the upper limit.A recent study by Kimura et al. (2023) indicated that no more than 0.1% of total flare energy output is converted into non-thermal protons, which is also consistent with our upper limit.
We note that the predicted photon spectral indices for Solar flares, using a power-law with exponential cutoff (PLEC) model, ranges from 3.5 to 4.5 (Kafexhiu et al. 2018).Due to the low detection significance of our results, we do not use PLEC models, only power-law models.For this reason, we cannot directly compare the upper limit to these model predictions.However, this range agrees with those -ray Solar flares from the FLSF catalog modeled with the PLEC spectrum.

Prospect of TeV Observations
Ohm & Hoischen (2018) suggested that TeV emission should be present from stellar flares.Very high energy observatories, such as SHALON, recently claimed detection of TeV emission from the direction from M dwarfs (Sinitsyna et al. 2019).The Cherenkov Telescope Array6 (CTA) should be sensitive enough to observe TeV stellar flares, and in fact, modeling from Ohm & Hoischen (2018) suggests superflares from DG CVn will be detectable by CTA.If Target of Opportunity (ToO) observation is adopted by the CTA consortium as part of the observing plans, it can be taken advantage of to detect TeV flares.Combined with the large number of flares anticipated (Clarke et al. 2024;Kowalski et al. 2009) from the Vera C. Rubin Observatory (Ivezić et al. 2019), and broker software such as Fink7 (Möller et al. 2021), CTA can quickly slew towards the flaring star for followup TeV observation.Targets for the ToO observations can be triggered follow-ups from ground based all-sky monitoring missions.Evryscope introduces a pipeline for low-latency transient detection which is suitable for detecting superflares (Corbett et al. 2023).These triggered events could potentially be used for lowlatency follow-up with CTA.Possibly included in these triggered follow-ups, TRAPPIST-1 would be an interesting source to focus on.At a distance of 12 pc and predicted to have 4 +1.9 −0.2 superflares per year (Glazier et al. 2020), it should be at least as detectable in TeV as DG CVn.Detecting GeV and TeV -ray emission around this planet-hosting star can help further understand habitability of exoplanets.Additionally, recent JWST observations of TRAPPIST-1 of transits during flares (Howard et al. 2023) could be useful for atmospheric characterization efforts.Multi-wavelength observations of TRAPPIST-1 planetary transits during flares could potentially be helpful towards these efforts.

CONCLUSION
In this work, we used sensitive stacking methods to search for any potential -ray emission associated with energetic stellar flares.Stacking the LAT data using the fullMET shows no detection, while gating the data around the flare times returns a sensitive upper limit of the flare -ray emission.The same analysis on empty test locations as a control returns a null result.Modeling this upper limit with Naima and comparing it with Solar -ray flares indicates that the common emission mechanism is likely neutral pion decay generated in the stellar atmosphere during flare events.To remain consistent with the stellar flare upper limit, the proton acceleration efficiency should not exceed 5%.

Figure 1 .
Figure 1.Top: Distributions of TS values for: stacked flare windows of each star (blue), test sources (green), and  2 /2 distribution with 2 d.o.f.Bottom left: parameter space stacking for the 298 test sources.Bottom right: parameter space stacking of 53 stars containing 298 observed flares with a TS peak at photon index of −2.5 and an upper limit in flux at 1.9 × 10 −10 cm −2 s −1 .The blue contour in the figure indicates TS = 2.71, which is traditionally used to indicate the 95% upper limit.The distribution of the PL index of all 26 FLSF cataloged solar flares with a PL spectral model is displayed on the right.

Figure 2 .
Figure 2. Flare frequency distributions of all flares investigated in this work and the FLSF Solar flares.Flare energies of the Solar flares are estimated as described in the text above using integrated SXR flux.Flare energies of the stellar flares are the values from the respective flare catalog.

Figure 3 .
Figure 3. Gamma-ray flux as a function of flare energy for the Solar flares and the stellar flares.Black data points are FLSF Solar flares plotted with SXR solar flare energy estimation; red data points are FLSF Solar flares that have SEP energy estimates from Aschwanden et al. (2017).The blue upper limit is the stellar -ray flux.Its flare energy value and uncertainty are the average and range of all the flare energies used in the stack.The overlapping green, yellow and purple shaded areas represent the results from the pion decay modeling implemented with Naima.Each shaded region represents a different target proton density, and their vertical limits depend on proton index: −3, bounded from above, and −6, bounded from below.

Table 1 .
Samples Drawn from Each Flare Survey Instrument Wavelength Max Dist.No. Stars No. Flares Reference