# Searching For Gamma-ray Emission from Stellar Flares

Yuzhe Song(宋宇哲),1,2,3,4,5 Timothy A. D. Paglione,3,4,5 and Ekaterina Ilin6
1Center for Astrophysics & Supercomputing, Swinburne University of Technology,
Hawthorn, VIC 3122, Australia
2OzGrav - ARC Centre for Excellence of Gravitational Wave Discovery, Hawthorn,
VIC 3122, Australia
3Department of Physics, the Graduate Center, City University of New York, 365
Fifth Ave., New York, NY 10016, USA
4Department of Earth & Physical Sciences, York College, City University of New
York, 94-20 Guy R. Brewer Blvd., Jamaica, NY 11451, USA
5Department of Astrophysics, American Museum of Natural History, Central Park
West at 79th Street, New York, NY 10024, USA
6ASTRON, Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4,
Dwingeloo, 7991 PD, The Netherlands
E-mail<EMAIL_ADDRESS>
(Accepted 2024 May 24. Received 2024 May 02; in original form 2024 January 23)

###### Abstract

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.

###### keywords:

stars: flare – methods: data analysis – gamma-rays: stars – stars: low mass

††pubyear: 2024††pagerange: Searching For Gamma-ray Emission from Stellar
Flares–References

##  1 Introduction

The Sun is a $\gamma$-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 $\gamma$-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 super- or even mega-flares as they reach
output energies over $10^{36}$ erg s-1. Such events are predicted to emit
$\gamma$-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 $\gamma$-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 TESS (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\sigma$ 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 $\gamma$-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
$\gamma$-rays. Further, development of $\gamma$-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 $\gamma$-ray data around the flare times, to search for
$\gamma$-rays from 1505 flares from 200 flaring dwarf stars.

##  2 Observations and Results

###  2.1 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
$\gamma$-ray emission within the stacking survey. Only stars with Galactic
latitude $>20$° were selected in this study to avoid the complicated
$\gamma$-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 $\sim 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.

Table 1: Samples Drawn from Each Flare Survey  Instrument | Wavelength | Max Dist. | No. Stars | No. Flares | Reference
---|---|---|---|---|---
Evryscope | Optical | 50 pc | 106 | 264 | Howard et al. (2019)
TESS | Optical | 25 pc | 86 | 1225 | Günther et al. (2020)
MAXI | X-ray | 88 pc | 8 | 16 | Tsuboi et al. (2016)

###  2.2 Fermi-LAT Data Analysis

The analysis in this work utilised the third revision of the PASS8 data, P8R3,
released on Nov 26, 2018, of $Fermi$-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.1111https://fermipy.readthedocs.io/en/latest/ based on the Fermi Science
tools Anaconda distribution (Anaconda, 2020) version
2.0.8222https://fermi.gsfc.nasa.gov/ssc/data/analysis/.

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.05 to 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 $\gamma$-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\degr\times 21.2\degr$ square and centered at the location
of a star, which corresponded to a $\sim 15\degr$ radius region of interest
(ROI). The data were also filtered using a zenith angle cut of $90\degr$ 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{{\cal L}/{\cal
L}_{0}}$, where ${\cal L}$ is the likelihood from the best-fit model
containing the star, and ${\cal L}_{0}$ is the likelihood of the model for the
null hypothesis.

###  2.3 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\times 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.

###  2.4 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
$\sim 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 $\chi^{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 $\gamma$-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, $\Delta$(TS) $=7$,
also indicates that the existence of flare $\gamma$-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 ${\cal L}_{\rm
max}-{\cal L}_{\rm max}|_{90}=2.71/2$, where ${\cal L}_{\rm max}$ is the
maximum likelihood of the null hypothesis, and ${\cal L}_{\rm max}|_{90}$ is
the maximum likelihood of the model with increased flux, or $\Delta$(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\times
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.

Figure 1: Top: Distributions of TS values for: stacked flare windows of each
star (blue), test sources (green), and $\chi^{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\times
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.

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
$\chi^{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 $\gamma$-ray background.

##  3 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
$\gamma$-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.

###  3.1 Flare Frequency Distribution

Since the Sun is the only star with individual flares observed in
$\gamma$-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 $E$, and typically follows a power law:

$f(>E)=\frac{\beta}{\alpha-1}E^{-\alpha+1},$ (1)

where $\alpha$ and $\beta$ are free parameters to be fitted. The power law
index $\alpha$ 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 -
2020333https://github.com/nplutino/FlareList. 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 $4\pi\textrm{AU}^{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
Altaipony444https://altaipony.readthedocs.io/en/latest/#id10 (Ilin, 2021; Ilin
et al., 2019; Davenport, 2016), in which $\alpha$ and $\beta$ 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 $\gamma$-ray Solar flares. The $\gamma$-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 $\alpha$ 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 $\alpha$ 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.

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.

###  3.2 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 $\gamma$-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 $\gamma$-ray flux from these flares and
constrain the PL index of their $\gamma$-ray emission.

As a simple comparison, we first scaled the $\gamma$-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 $\gamma$-ray
flux of $1.42\times 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 $\gamma$-ray
flux of the FLSF solar flares using the estimated Solar flare energy described
in § 3.1, and

$\frac{E_{\rm stellar}}{E_{\rm solar}}\times(\frac{\rm AU}{\rm d})^{2},$ (2)

where $E_{\rm stellar}=2.3\times 10^{33}$ erg is the median flare energy of
the sample stellar flares, $E_{\rm solar}$ is the flare energy of any given
Solar flare, and $d=25$ pc is the average distance to the target stars. The
range of these scaled Solar flare $\gamma$-ray fluxes is $1\times 10^{-16}$ to
$3\times 10^{-11}$ cm-2 s-1, which is below the LAT sensitivity limit of $\sim
10^{-9}$ cm-2 s-1, but overlaps with the fluxes probed by our stacking method.

###  3.3 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)555https://naima.readthedocs.io/en/latest/index.html, a Python package
that computes radiation from non-thermal particle populations and also does
MCMC fitting to spectra. Solar $\gamma$-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
$\Gamma_{\rm p}$ ranging from $-6$ to $-3.2$ yield spectra consistent with
Solar $\gamma$-ray flares (Ajello et al., 2021). We use input proton spectra
with $\Gamma_{\rm 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
$\Gamma_{\rm 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
$\Gamma_{\rm p}=-6$ and $-3$, respectively. The $\gamma$-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
$\gamma$-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.

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 $\gamma$-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.

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 $\gamma$-ray Solar flares from the FLSF catalog
modeled with the PLEC spectrum.

###  3.4 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 Array666https://www.cta-
observatory.org/ (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 Fink777https://fink-broker.org/ (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 low-latency 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 $\gamma$-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.

##  4 Conclusion

In this work, we used sensitive stacking methods to search for any potential
$\gamma$-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 $\gamma$-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
$\gamma$-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%.

## acknowledgements

This work was supported in part by NSF grant AST-1831412 and the PSC-CUNY
award # 63785-00 51. The first author is also supported by the Australian
Research Council Centre of Excellence for Gravitational Wave Discovery
(OzGrav), through project number CE170100004. This project made use of
computational systems and network services at the American Museum of Natural
History supported by the National Science Foundation via Campus
Cyberinfrastructure Grant Awards # 1827153 (CC$*$ Networking Infrastructure:
High Performance Research Data Infrastructure at the American Museum of
Natural History).

## Data Availability

Due to the size of output files for Fermi-LAT data analysis, the authors can
share configuration files used to process all LAT data. Post processing
scripts, including the stacking analysis, will be hosted on Github as well.
Before made available on Github, they will be available upon request provide
this work is properly cited.

## References

  * Abdo et al. (2009) Abdo A. A., et al., 2009, ApJ, 703, 1249
  * Abdo et al. (2011) Abdo A. A., et al., 2011, ApJ, 734, 116
  * Abdollahi et al. (2020) Abdollahi S., et al., 2020, ApJS, 247, 33
  * Ackermann et al. (2017) Ackermann M., et al., 2017, ApJ, 835, 219
  * Ajello et al. (2014) Ajello M., et al., 2014, ApJ, 789, 20
  * Ajello et al. (2020) Ajello M., Di Mauro M., Paliya V. S., Garrappa S., 2020, ApJ, 894, 88
  * Ajello et al. (2021) Ajello M., et al., 2021, ApJS, 252, 13
  * Anaconda (2020) Anaconda 2020, Anaconda Software Distribution, https://docs.anaconda.com/
  * Aschwanden et al. (2017) Aschwanden M. J., et al., 2017, ApJ, 836, 17
  * Chen et al. (2021) Chen H., Zhan Z., Youngblood A., Wolf E. T., Feinstein A. D., Horton D. E., 2021, Nature Astronomy, 5, 298
  * Clarke et al. (2024) Clarke R. W., et al., 2024, arXiv e-prints, p. arXiv:2402.06002
  * Corbett et al. (2023) Corbett H., et al., 2023, ApJS, 265, 63
  * Cowan (1997) Cowan G., 1997, Statistical data analysis
  * Davenport (2016) Davenport J. R. A., 2016, ApJ, 829, 23
  * Dobie et al. (2023) Dobie D., et al., 2023, MNRAS, 519, 4684
  * Glazier et al. (2020) Glazier A. L., Howard W. S., Corbett H., Law N. M., Ratzloff J. K., Fors O., del Ser D., 2020, ApJ, 900, 27
  * Günther et al. (2020) Günther M. N., et al., 2020, AJ, 159, 60
  * Hawley et al. (2014) Hawley S. L., Davenport J. R. A., Kowalski A. F., Wisniewski J. P., Hebb L., Deitrick R., Hilton E. J., 2014, ApJ, 797, 121
  * Hazra et al. (2022) Hazra G., Vidotto A. A., Carolan S., Villarreal D’Angelo C., Manchester W., 2022, MNRAS, 509, 5858
  * Howard et al. (2019) Howard W. S., Corbett H., Law N. M., Ratzloff J. K., Glazier A., Fors O., del Ser D., Haislip J., 2019, ApJ, 881, 9
  * Howard et al. (2023) Howard W. S., et al., 2023, ApJ, 959, 64
  * Huber et al. (2012) Huber B., Farnier C., Manalaysay A., Straumann U., Walter R., 2012, A&A, 547, A102
  * Ilin (2021) Ilin E., 2021, The Journal of Open Source Software, 6, 2845
  * Ilin et al. (2019) Ilin E., Schmidt S. J., Davenport J. R. A., Strassmeier K. G., 2019, A&A, 622, A133
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Kafexhiu et al. (2014) Kafexhiu E., Aharonian F., Taylor A. M., Vila G. S., 2014, Phys. Rev. D, 90, 123014
  * Kafexhiu et al. (2018) Kafexhiu E., Romoli C., Taylor A. M., Aharonian F., 2018, ApJ, 864, 148
  * Kimura et al. (2023) Kimura S. S., Takasao S., Tomida K., 2023, ApJ, 944, 192
  * Kowalski et al. (2009) Kowalski A. F., Hawley S. L., Hilton E. J., Becker A. C., West A. A., Bochanski J. J., Sesar B., 2009, AJ, 138, 633
  * Möller et al. (2021) Möller A., et al., 2021, MNRAS, 501, 3272
  * Ohm & Hoischen (2018) Ohm S., Hoischen C., 2018, MNRAS, 474, 1335
  * Orlando & Strong (2008) Orlando E., Strong A. W., 2008, A&A, 480, 847
  * Paliya et al. (2019) Paliya V. S., Domínguez A., Ajello M., Franckowiak A., Hartmann D., 2019, ApJ, 882, L3
  * Paliya et al. (2020) Paliya V. S., Ajello M., Cao H. M., Giroletti M., Kaur A., Madejski G., Lott B., Hartmann D., 2020, ApJ, 897, 177
  * Paudel et al. (2018) Paudel I., Bar-Tal A., Levy G. J., Rotbart N., Ephrath J. E., Cohen S., 2018, Agricultural Water Management, 204, 126
  * Plutino et al. (2023) Plutino N., Berrilli F., Del Moro D., Giovannelli L., 2023, Advances in Space Research, 71, 2048
  * Principe et al. (2021) Principe G., Di Venere L., Orienti M., Migliori G., D’Ammando F., Mazziotta M. N., Giroletti M., 2021, MNRAS, 507, 4564
  * Rigney et al. (2022) Rigney J., et al., 2022, MNRAS, 516, 540
  * Seckel et al. (1991) Seckel D., Stanev T., Gaisser T. K., 1991, ApJ, 382, 652
  * Shakhovskaia (1989) Shakhovskaia N. I., 1989, Sol. Phys., 121, 375
  * Sinitsyna et al. (2019) Sinitsyna V. G., Sinitsyna V. Y., Stozhkov Y. I., 2019, Advances in Space Research, 64, 2585
  * Song & Paglione (2020) Song Y., Paglione T. A. D., 2020, ApJ, 900, 185
  * Song et al. (2023) Song Y., Paglione T. A. D., Tan J., Lee-Georgescu C., Herrera D., 2023, MNRAS, 524, 5854
  * Tsuboi et al. (2016) Tsuboi Y., et al., 2016, PASJ, 68, 90
  * Varela et al. (2022) Varela J., Brun A. S., Strugarek A., Réville V., Zarka P., Pantellini F., 2022, A&A, 659, A10
  * Wheatland (2004) Wheatland M. S., 2004, ApJ, 609, 1134
  * Wood et al. (2017) Wood M., Caputo R., Charles E., Di Mauro M., Magill J., Perkins J. S., Fermi-LAT Collaboration 2017, in 35th International Cosmic Ray Conference (ICRC2017). p. 824 (arXiv:1707.09551), doi:10.22323/1.301.0824
  * Zabalza (2015) Zabalza V., 2015, in 34th International Cosmic Ray Conference (ICRC2015). p. 922 (arXiv:1509.03319), doi:10.22323/1.236.0922