# Kepler Flares I. Active and Inactive M dwarfs

Suzanne L. Hawley11affiliation: Corresponding author<EMAIL_ADDRESS>22affiliation: Department of Astronomy, University of Washington, Box 351580,
Seattle, WA 98195 , James R. A. Davenport22affiliation: Department of
Astronomy, University of Washington, Box 351580, Seattle, WA 98195 , Adam F.
Kowalski22affiliation: Department of Astronomy, University of Washington, Box
351580, Seattle, WA 98195 33affiliation: NASA Goddard Space Flight Center,
Code 671, Greenbelt, MD 20771 , John P. Wisniewski22affiliation: Department of
Astronomy, University of Washington, Box 351580, Seattle, WA 98195
44affiliation: HL Dodge Department of Physics & Astronomy, University of
Oklahoma, 440 W. Brooks Street, Norman, OK 73019 ,
Leslie Hebb55affiliation: Department of Physics, Hobart and William Smith
Colleges, 300 Pulteney Street, Geneva, NY 14456 , Russell
Deitrick22affiliation: Department of Astronomy, University of Washington, Box
351580, Seattle, WA 98195 , Eric J. Hilton22affiliation: Department of
Astronomy, University of Washington, Box 351580, Seattle, WA 98195
66affiliation: Universe Sandbox, 911 E. Pike Street #333, Seattle, WA 98122

###### Abstract

We analyzed Kepler short-cadence M dwarf observations. Spectra from the ARC
3.5m telescope identify magnetically active (H$\alpha$ in emission) stars. The
active stars are of mid-M spectral type, have numerous flares, and well-
defined rotational modulation due to starspots. The inactive stars are of
early-M type, exhibit less starspot signature, and have fewer flares. A Kepler
to U-band energy scaling allows comparison of the Kepler flare frequency
distributions with previous ground-based data. M dwarfs span a large range of
flare frequency and energy, blurring the distinction between active and
inactive stars designated solely by the presence of H$\alpha$. We analyzed
classical and complex (multiple peak) flares on GJ 1243, finding strong
correlations between flare energy, amplitude, duration and decay time, with
only a weak dependence on rise time. Complex flares last longer and have
higher energy at the same amplitude, and higher energy flares are more likely
to be complex. A power law fits the energy distribution for flares with log
$E_{K_{p}}>$ 31 ergs, but the predicted number of low energy flares far
exceeds the number observed, at energies where flares are still easily
detectable, indicating that the power law distribution may flatten at low
energy. There is no correlation of flare occurrence or energy with starspot
phase; the flare waiting time distribution is consistent with flares occurring
randomly in time; and the energies of consecutive flares are uncorrelated.
These observations support a scenario where many independent active regions on
the stellar surface are contributing to the observed flare rate.

##  1\. Introduction

The Kepler satellite (Borucki et al., 2010) has ushered in a new era of
stellar photometric investigation, enabling light curve analysis with
unprecedented precision: approaching 10 ppm in bright targets (V=9-10) and 100
ppm even in faint targets (V=13-14). Light curves of this precision, over
periods of months to years with minimal interruption, provide a unique view of
stellar variability in optical light. Here, we are interested in the
variability that results from flares on low mass stars.

The advantages of Kepler over ground-based observing of stellar flares are
many. Since flares occur at unpredictable intervals, a long time baseline
uninterrupted by daylight, weather or the vagaries of telescope scheduling
greatly facilitates the determination of the flare frequency distribution
(FFD, the cumulative number of flares above a given energy that occur per unit
time). FFDs on low mass stars have been extensively examined from the ground,
beginning with the pioneering work of Gershberg (1972) and the seminal paper
by Lacy et al. (1976, hereafter LME) which presented results for 8 stars using
more than 400 hours (spread over two years) of ground-based monitoring
observations by Moffett (1974). These early studies showed that flares occur
more frequently but with lower energy on mid-M dwarfs (M3-5) compared to
earlier type, M0-M2 stars. However, there are significant selection effects at
both extremes of the distribution due to a) the relatively higher quiescent
luminosity of the earlier type stars, which makes it more difficult to detect
low energy flares; and b) the limited duration of the monitoring observations
which makes it less likely to observe the relatively rare, high energy flares.
Recently, Hilton (2011) compiled several hundred hours of ground-based
monitoring of M dwarfs and extended the investigation of FFDs to the regimes
of very low mass (VLM) dwarfs of type M6-M8 and of inactive, early type M
dwarfs. His results indicate that the trend of higher frequency but lower
energy flares extends to the latest M dwarfs, and that even “inactive” (no
H$\alpha$ emission in the quiescent spectrum) early-type M stars still flare,
although at much lower frequency than active stars of the same spectral type.
Audard et al. (2000) presented FFDs for several M dwarfs using X-ray data from
the EUVE satellite and obtained similar results to the optical studies. HST
investigations in the NUV (Robinson et al., 1995) and FUV (Robinson et al.,
2001) revealed smooth extension of the FFDs to the microflare regime for both
early and late M dwarfs. In the first reported long-duration optical
monitoring observations from space, Hunt-Walker et al. (2012) used data from
the MOST satellite on the active M3 star AD Leo to determine that the FFD from
a week of nearly continuous observation matched those previously found from
the ground.

With a power-law probability distribution for flare energies

$N(E)dE\propto E^{-\alpha}dE,$ (1)

a linear fit to the cumulative FFD gives a slope $-\alpha+1$. There is
considerable interest in the solar literature concerning the value of
$\alpha$, since a flare energy distribution with $\alpha>2$ extending to low
flare energies would potentially provide enough energy from flares to heat the
solar corona (Parker, 1988; Hudson, 1991; Hannah et al., 2011; Schrijver et
al., 2012). The M dwarf studies described above typically find flatter power
law distributions (smaller $\alpha$) than for the Sun, at least for ground-
based data. Güdel (2004) reviewed X-ray results and noted several studies that
obtain $\alpha>2$. Ramsay et al. (2013) recently presented FFDs from Kepler
data for two M dwarfs (including the active M4 star GJ 1243 which is contained
in our sample), and found steeper power-law slopes, also implying $\alpha>2$.
By contrast, the Kepler study of an early L dwarf examined by Gizis et al.
(2013) showed a relatively shallow power-law slope in its FFD.

The assurance that every flare has been observed over a long time period of
consecutive observation (approximately one month with Kepler) also provides
excellent statistics for the study of flare timing, for example the waiting
time between flares. The waiting time distribution can be used to address
whether flares occur independently, or are triggered by other flares
(sympathetic flaring), and also whether the phenomenon of precursor flares
(small flares immediately followed by large flares) is real. LME found that
the flares in their sample were consistent with random occurrence, but other
studies have presented evidence for sympathetic flaring on low mass stars
(Pazzani & Rodono, 1981; Doyle et al., 1990; Panagi & Andrews, 1995), RS CVn
stars (Osten & Brown, 1999), and on the Sun (Pearce & Harrison, 1990).
Wheatland (2010) investigated flares within a single solar active region, and
found that they followed a piecewise Poisson distribution, indicating random
occurrence, but with different rates as the magnetic topology of the region
evolved.

In addition, the precision of the Kepler light curves provides an
unprecedented view of periodic variations attributed to starspots rotating
across the visible hemisphere. Starspots on M dwarfs have previously been
difficult to measure due to the very small amplitudes ($<0.01$ mag) that are
typically observed (Irwin et al., 2011). However, variations of this amplitude
are easily measured with Kepler. It is well known on the Sun that flares are
associated with active regions (sunspots). The Kepler data now make it
possible to look for a correlation between starspot phase and flare occurrence
in low mass stars.

Finally, since Kepler monitors many stars simultaneously, it is possible to
compare these various measures of flare energy, frequency and timing in stars
of different mass, temperature and evolutionary state as well as in otherwise
similar stars that exhibit very different levels of magnetic activity. A
series of papers reporting spectropolarimetric observations of M dwarfs (Morin
et al., 2008; Donati et al., 2008; Morin et al., 2010) describes the change in
the large-scale magnetic field structure between early M dwarfs that retain a
radiative core, presumably have a solar-type dynamo and exhibit weak, mostly
toroidal field; and mid M dwarfs, which are fully convective and have a
turbulent dynamo, and exhibit stronger, more stable, poloidal fields. As these
authors point out, the large scale fields they infer from their
spectropolarimetry have relatively low field strength and most of the magnetic
energy in the active stars, with measured field strengths of several kG
(Johns-Krull & Valenti, 1996), must be contained in smaller, spatially
confined active regions that are unresolved with global polarization
measurements. Reiners & Basri (2009) reported Stokes V and I measurements of
early and mid M dwarfs and find that the large scale fields contribute at most
15% of the total magnetic flux. Investigation of flare activity in both active
and inactive stars of early and mid M dwarfs may provide additional insight
into the relative importance of the large and small scale fields, especially
at spectral types near M4, where stars transition to a fully convective
interior and turbulent dynamo (see discussion in Stassun et al. (2011)).

We have carried out a long term program of flare monitoring with Kepler in the
short (one-minute) cadence mode, with the ultimate goal of obtaining a deeper
understanding of the flare mechanism, how and when flares occur and how flare
properties are correlated with the underlying stars. This is the first in a
series of papers describing our results. In this paper (Paper 1), we focus on
the lowest mass stars in our Kepler program, active and inactive M dwarfs of
spectral type M1-M5. Davenport et al. (2014, hereafter Paper 2) examines flare
light curve morphology using a large sample of flares on the active M dwarf GJ
1243. Lurie et al (2014 in prep., hereafter Paper 3), investigates the active
binary system GJ 1245AB using the Kepler pixel files to separate the two stars
and analyze the individual light curves. Additional papers are also planned.

Section 2 of this paper describes the M dwarf sample and the process we use to
identify flares. Section 3 presents the results for flare frequencies and
energies, and the comparison between the active and inactive stars. Sections 4
and 5 concentrate on the active single star GJ 1243, and describe the
relationships between basic properties such as flare amplitude, duration and
energy, as well as our analysis of flare timing and correlation with starspot
phase. We offer some conclusions and discussion of future work in Section 6.

##  2\. Kepler M Dwarf Flare Sample

Figure 1.— Moderate resolution spectra obtained with the ARC 3.5m Telescope
and DIS spectrograph for the M dwarfs in the Kepler flare sample.

Figure 2.— Light curves for the active M dwarf GJ 1243 (top) and the active
binary GJ 1245AB (bottom). A flare with energy log $E_{Kp}\sim$30.5 is marked
with a vertical blue bar, and a flare with energy log $E_{Kp}\sim$31 is marked
with a vertical red bar; these are discussed in §4.1. Although GJ 1245AB is
brighter, it exhibits more noise likely due to variations in the amount of
total light from the binary falling within the Kepler photometric aperture
(see text).

Figure 3.— Light curves for the inactive M dwarfs GJ 4099 (M1), GJ 4113 (M2)
and GJ4083 (M3). Note the increased noise for GJ 4083, which is more than a
magnitude fainter than the other two inactive stars. Flares that were
identified in these light curve segments are marked with vertical red lines.

We chose the most active and brightest late type M dwarfs in the Kepler field,
GJ 1243 (dM4e) and GJ 1245AB (dM5e+dM5e) and the latest type inactive (no
H$\alpha$ emission) M dwarfs (GJ 4083, dM3; GJ 4099, dM1; and GJ 4113, dM2)
that were bright enough to give good signal in one-minute cadence data. To our
knowledge, all of the objects are single field stars, except for the GJ 1245AB
system which is comprised of two nearly identical active mid-M dwarfs
separated by 7”, and a third component (GJ 1245C), an M8 companion to GJ
1245A, separated by 0$\farcs$6 (Schroeder et al., 2000). The entire system is
measured together in the Kepler pipeline photometry. Since GJ 1245C is
approximately 3 magnitudes fainter than the GJ 1245AB system, we do not expect
that it contributes significantly to the combined light observed with Kepler.
In particular, a 3 magnitude flare in the Kepler bandpass on an M8 dwarf will
occur less than once every 100 hours according to results from Davenport et
al. (2012) and Hilton (2011), and such flares would only appear as brief,
short-lived, low-energy excursions in the GJ 1245AB light curve. We therefore
treat GJ1245AB as a binary system and analyze the light curve in terms of the
A and B components. In Paper 3 of this series we analyze the GJ 1245 system in
more detail and separate the components using recalibrated light curves from
the Kepler pixel files.

Table 1 gives the basic properties of the sample targets, including our
determination of rotation periods which we describe below. We obtained new
spectral types using moderate resolution spectra obtained with the DIS
spectrograph on the ARC 3.5m telescope at Apache Point Observatory. Figure 1
shows the six spectra (GJ 1245 A and B were resolved so that we were able to
take spectra of the individual components). The spectral typing followed the
procedure of Covey et al. (2007) using the Hammer software described therein.
Distances for the inactive stars in Table 1 were determined using published
$R-I$ colors transformed to SDSS $r-i$ (Davenport et al., 2006) and then
employing the $M_{r}$ vs $r-i$ photometric parallax relation given in
Bochanski et al. (2010), with the updated values given in the erratum to that
paper (Bochanski et al., 2012). The active stars have measured parallaxes as
reported in Reid et al. (1995).

###  2.1. Kepler Light Curves

Our Kepler Cycle 2 GO program 20016 comprised two months of one-minute cadence
data for each of the objects in the sample. Note that the one-minute cadence
is important to obtain well-sampled flare light curves, since many flares on M
dwarfs have a duration shorter than 30 minutes (e.g. Moffett, 1974) and nearly
all are shorter than the 90 minutes (3 observations of 30 minutes each)
required in e.g. the Walkowicz et al. (2011) discussion of stellar flares in
Kepler long cadence data.

We used the most recent publicly available version of the Kepler photometry,
which gives the Pre-search Data Conditioning Simple Aperture Photometry
(PDCSAP) flux (Smith et al., 2012). The PDCSAP fluxes are the result of a
systematic reprocessing of the entire Kepler database using a Bayesian
approach to remove systematics from the short and long term light curves.
However, some month to month corrections still needed to be manually applied.
We used a linear least-squares fit to each month of data to remove low-order
flux variations, and normalized the median flux levels between months.

Figure 2 illustrates representative 3-day segments from the light curves for
the active stars, with the fractional flux $\Delta F/F$ given by:

$\Delta F/F=(F_{i}-F_{o})/F_{o}$ (2)

where $F_{i}$ is the measured flux sampled at each time $i$, and $F_{o}$ is
the normalized median flux as described above.

The light curves display a simple, well-defined, periodic flux variation in GJ
1243 and a more complicated flux variation for the binary GJ 1245AB, due to
Kepler measuring the combined light from both stars. The periodic flux
variations are ascribed to starspots moving in and out of view on the visible
hemisphere as the star rotates. These starspot variations are expected on
magnetically active stars, and have been observed commonly on active G and K
stars even from the ground, with flux variations up to 10%. M dwarfs exhibit
much smaller flux variations due to starspots; even the very active stars only
show variations of $\sim$1%. If much of the star is covered with strong
magnetic fields, as discussed in Johns-Krull & Valenti (1996), then numerous
small (or a few large) spots may always be present on the surface, and thus
only small variations as the star rotates would be expected.

Many flares are evident even in these short data segments. For reference,
flares with energies in the Kepler bandpass (see §2.3 for discussion of energy
measurements) of log $E_{Kp}$ = 30.5 and 31 are indicated with vertical blue
and red bars respectively. These particular flare energies will be discussed
further in §4.1 below.

The increased noise in the GJ 1245AB light curve compared to GJ 1243 is
unusual since GJ 1245AB is a brighter target. However, inspection of the pixel
files for GJ 1245AB reveal that it is partially resolved, with much of the
light from the A component in particular falling outside the aperture chosen
for the pipeline photometry. Additional noise results from varying amounts of
light from the two stars falling within the Kepler photometric aperture due to
spacecraft jitter and other guiding issues. These issues also affect the flux
normalization, flare identification and flare energies assigned to GJ 1245AB,
adding uncertainty to the results obtained from the combined system.

Figure 3 shows 11-day segments of the inactive star light curves, indicating
far less flux variation due to starspots (note the flux scale is a factor of
ten reduced compared to Figure 2) and a few small flare events, marked with
red vertical bars. It is notable that even in these inactive (no H$\alpha$
emission) stars, we still see evidence for some magnetic activity from
starspots and flares.

Figures 4 and 5 give the corresponding Lomb-Scargle periodograms for the
2-month Kepler dataset for each target. In Figure 4, GJ 1243 (P=0.59 days) and
GJ 1245AB (P = 0.26 days for component A and P = 0.71 days for component B)
have very strong peaks in their periodograms indicating clearly periodic
variations which we attribute to rotational modulation from stable starspots
on the surfaces of these stars. Our periods for these stars are reported in
Table 1, with uncertainties found from comparing values obtained when
analyzing separate months of data (the formal uncertainties on the individual
fits are much smaller). Our periods agree with previous determinations for GJ
1243 (Irwin et al., 2011; McQuillan et al., 2013) and GJ 1245B (Morin et al.,
2010). To our knowledge this is the first measurement of the period of GJ
1245A. Note that the combined Kepler pipeline light curve for GJ1245AB does
not allow us to independently identify which period belongs to which star, but
our additional pixel-level investigation in Paper 3 verifies that the periods
are ascribed to the correct components, and provides the uncertainties quoted
on those periods.

Table 1Properties of the Kepler M dwarf flare sample

Star | KIC | Sp Type | $m_{Kp}$ | Active | EW H$\alpha$ (Å) | Dist (pc) | Period† (d)
---|---|---|---|---|---|---|---
GJ 4099 | 004142913 | M1 | 10.9 | n | … | 16.2⋆ | 2.0:
GJ 4113 | 004470937 | M2 | 10.9 | n | … | 14.2⋆ | …
GJ 4083 | 010647081 | M3 | 12.3 | n | … | 18.2⋆ | 3.3:
GJ 1243 | 009726699 | M4e | 12.7 | y | 4.0 | 12.05 | 0.5927 $\pm$ 0.0002
GJ 1245A | 008451881 | M5e | 11.7+ | y | 2.5 | 4.55 | 0.2632 $\pm$ 0.0001
GJ 1245B | 008451881 | M5e | 11.7+ | y | 3.8 | 4.55 | 0.709 $\pm$ 0.001

$\star$$\star$footnotetext: Distances determined via photometric parallax

++footnotetext: combined magnitude for both components

$\dagger$$\dagger$footnotetext: periods determined in this paper, see §2.1

The inactive stars GJ 4083 and GJ 4099 (see Figure 5) have several peaks of
about the same weak power (an order of magnitude weaker than in the active
stars) in their periodograms. These may be due to small transient starspots
appearing on the stellar surface, possibly at different latitudes which may
indicate differential rotation. We report the period of the strongest peak in
Table 1 for GJ 4083 and GJ 4099 (as indicated by red vertical lines in Figure
5), but the periods are very uncertain, and should only be regarded as
suggestive. The final inactive star, GJ 4113, does not show any significant
peak in its periodogram, and we are unable to determine its rotation period.
We also note that our analysis of two months of short cadence data is
sensitive only to periods less than 10 days. Previous examination of long
cadence Kepler data revealed marginal evidence for longer periods in the
inactive stars in our sample (McQuillan et al., 2013).

Figure 4.— Lomb-Scargle periodograms obtained from the Kepler light curves for
the active star GJ 1243 (top) and the active binary GJ 1245AB (bottom). The
insets show expanded views of the strongest peaks.

Figure 5.— Lomb-Scargle periodograms obtained from the Kepler light curves for
the inactive stars GJ 4099 (M1, top), GJ 4113 (M2, middle) and GJ 4083 (M3,
bottom). The peaks representing the periods reported in Table 1 for GJ 4099
and GJ 4083 are indicated. We were unable to determine a period for GJ 4113.

###  2.2. Flare Identification

In order to identify flares, we used a modified version of the flare-finding
approach described in Hilton (2011) and Hunt-Walker et al. (2012). A full
description of our flare finding methodology is given in Paper 2, and we
briefly summarize the algorithm here. Each month of Kepler data was treated
independently. The light curve was first de-trended to remove starspot
modulations using a variable-span smoothing algorithm, similar to the
“Supersmoother” method (Friedman, 1984). The local mean and standard deviation
($\sigma$) of the flux were then computed at each time sampled, and candidate
flares were flagged as excursions of two or more consecutive points more than
2.5$\sigma$ above the local mean flux. Each month of data was independently
inspected by at least two people, and both had to agree that a candidate flare
was real for it to be included in the sample. A flare type (classical or
complex) is assigned based on whether there is only one peak in the light
curve (classical) or more than one peak (complex), and these types were also
manually verified. Finally, the assigned start and end times for each flare
were examined, and in some cases modified, during the inspection process.

The shape of a typical flare light curve (the light curve temporal morphology)
is discussed in detail for GJ 1243 in Paper 2. Here we simply show in Figure 6
an idealized, classical (only one peak) flare with the usual fast rise and
fast decay (impulsive phase) followed by a slower exponential decay (gradual
phase) light curve, to illustrate the measurements we obtain directly from the
Kepler data for each flare – the amplitude, rise and decay times, duration and
equivalent duration. The amplitude of the flare is defined by the point in
time with the highest flux value (the flare peak), and is expressed as a
fraction of the local mean flux,

$Amplitude=(F_{peak}-F_{localmean})/F_{localmean}$ (3)

similar to the fractional flux $\Delta F/F$ defined for the light curves
(§2.1, Equation 2) but using the local mean flux rather than the global median
flux as the reference. The amplitude of the flare is thus measured relative to
the current state of the underlying star, including effects from starspots,
and represents the excess emission above the local mean flux. For example, the
idealized flare shown with an amplitude of 1 in Figure 6 corresponds to a peak
flux enhancement that is twice as large as the local mean flux. The duration
of the flare is the difference between the start time (the point in time when
the flare flux begins to deviate from the local mean flux) and the end time
(when the flare flux returns to the local mean flux). The start and end times
for each vetted flare were obtained during the manual inspection step,
including any point where at least two people agreed it belonged in the flare.
The rise time is then the time between flare start and flare peak, the decay
time is the time between flare peak and flare end, and the duration is by
definition the sum of the rise and decay times. The measured rise and decay
times and durations are given in integer minutes due to the one-minute
sampling in the Kepler data. Finally, the equivalent duration is the area
under the flare light curve, i.e. the fractional flux integrated over the
flare duration as illustrated by the shaded area in Figure 6. It is directly
related to the flare energy as described in the next section.

We note that the time sampling and the presence of noise in the light curve
means that the detection of the first point in the initial flare rise, and the
last point in the flare decay are typically underestimated compared to an
analytical flare model. However, comparison with the model for classical
flares discussed in Paper 2 shows that the main effect is an underestimation
of flare decay time and duration by about 20%, and this is independent of the
size of the flare. Further, since the empirical measurements effectively cut
off only the last part of the long exponential tail, there is almost no effect
on the energy of the flare. We have chosen to use the empirical measurements
in this paper, particularly so that complex and classical flares can be
uniformly compared.

The quarter of Kepler data used, the total length of the one-minute cadence
data that were analyzed, and the number of flares found for each target are
listed in Table 2. The quiescent Kepler luminosity and the range of flare
energies found for each object in the sample are also given in Table 2 and
described in the next section.

Figure 6.— Idealized classical flare light curve showing measured quantities –
amplitude, rise and decay times, duration and equivalent duration.

###  2.3. Determination of Flare Energies

The flare energy is calculated using the equivalent duration (Gershberg, 1972;
Hunt-Walker et al., 2012), which is the area under the flare light curve as
described above and shown in Figure 6. The equivalent duration is measured in
units of time; multiplying it by the quiescent luminosity of the star in the
Kepler bandpass gives the Kepler flare energy, which we denote as $E_{Kp}$ .

Obtaining the quiescent luminosity in the Kepler bandpass, $L_{Kp}$ for each
of our targets required several steps. First, we calculate the quiescent
Kepler flux (at Earth) using the measured apparent magnitude $m_{Kp}$ and the
Kepler zero-point magnitude. To find the zero-point, we convolved
spectrophotometric data for GJ 1243 from Kowalski et al. (2013) with the
Kepler filter curve to obtain the specific flux in the Kepler bandpass,
6.66$\times 10^{-14}$ ergs sec-1 cm-2 Å-1. With $m_{Kp}$ = 12.7 for GJ 1243,
we solved for the Kepler zero-point magnitude ${Kp_{0}}$ = -20.24. The
resulting Kepler fluxes for each target were multiplied by 4$\pi{d^{2}}$ and
the Kepler filter FWHM (4000Å) to give the Kepler quiescent luminosities
$L_{Kp}$ shown in Table 2. As a consistency check, we convolved the spectra
shown in Figure 1 with the Kepler filter curve to obtain measured Kepler
fluxes, and found that these were $\sim$30-80% of the values obtained from the
measured apparent magnitudes using the method outlined above. This was as
expected due to cloud cover and slit losses during the non-photometric
conditions when the spectral data were obtained.

A conservative estimate leads to errors of $\pm$0.2 in the final values of log
$E_{Kp}$, due to uncertainties in the Kepler zero-point, the measured Kepler
magnitudes, the distance estimates and the measurement of the flare equivalent
durations from the Kepler photometry.

The final step in the data reduction is the transformation between the flare
energy in the Kepler bandpass and in a broad-band filter such as the Johnson
U-band, traditionally used in flare studies because flares are quite blue
compared to the underlying stellar light. However, because equivalent
durations change in different bandpasses, we cannot simply multiply the Kepler
equivalent duration by the quiescent U band luminosity to obtain $E_{U}$.
Instead we determined the relationship between $E_{U}$ and $E_{Kp}$ using
simultaneous data obtained for a flare on GJ 1243 with both Kepler and in the
$U$-band with the NMSU 1-meter telescope at Apache Point Observatory. The
flare light curves are shown in Figure 7. The equivalent durations differ by
nearly two orders of magnitude, with $ED_{U}$ = 422 seconds, and $ED_{Kp}$ =
4.4 seconds. GJ 1243 does not have a published $U$ magnitude, so we adopted
U-B = 0.93 (the value for YZ CMi, a famous nearby flare star of similar
spectral type), to obtain log $L_{U}$ = 28.50 erg s-1. The flare energies are
then log $E_{Kp}$ = 31.31 and log $E_{U}$ = 31.10, or $E_{U}$ = 0.65$E_{Kp}$.

The energy scaling between the U-band and Kepler filter may vary between
flares for several reasons. The strong blue emission which characterizes
ground-based observations of flares is emitted primarily during the flare
impulsive phase, which varies in strength and duration (see e.g. Kowalski et
al., 2013). Flares with a very strong and long-lived impulsive phase will emit
relatively more energy in the U-band. During the flare gradual phase, Kowalski
et al. (2013) have recently discussed the existence of red continuum emission
in some flares which will contribute to the emission in the Kepler filter if
present. Flares that have been observed with many filters simultaneously are
rare, but a very large flare observed on the dM3e star AD Leo simultaneously
in the Johnson UBVR filters was reported in Hawley & Pettersen (1991). Their
Table 6 indicates that the sum of the B+V+R filter energies covers the Kepler
bandpass from $\sim$4000-8000Å, and gives $E_{U}$ = 0.4($E_{B}+E_{V}+E_{R}$)
$\sim$ 0.4($E_{Kp}$). Flare energy scaling relations have also been determined
by LME (and confirmed by Hilton 2011) for a few flares observed simultaneously
with U, B and V filters. LME found that $E_{U}=1.2E_{B}$ and $E_{U}=1.8E_{V}$.
Again using the sum of the B, V and R filters to approximate the Kepler
filter, these relations give $E_{U}=0.6E_{Kp}$ if $E_{R}=0.5E_{V}$ or
$E_{U}=0.4E_{Kp}$ if $E_{R}=2E_{V}$ (the value found in Hawley & Pettersen
(1991) for the large AD Leo flare). Clearly, the presence and strength of a
red continuum component in a flare will strongly affect the resulting energy
scaling.

For our purposes, we have chosen to adopt the scaling $E_{U}=0.65E_{Kp}$
determined from the direct measurement. However, variations in the scaling
factor from 0.4 to 0.65 lead to a range of $\sim$0.2 in log $E_{U}$, and we
display this range in our discussion in the next section.

##  3\. Flare Frequency on Active and Inactive M Dwarfs

Figure 8 illustrates our observed flare frequency distributions (FFDs), which
are diagrams of cumulative flare frequency (log number of flares per day with
energy greater than E) as a function of flare energy (log E).111This diagram
is similar to those shown in LME but with the energy as the independent
variable, following more recent studies (Audard et al., 2000; Hilton, 2011;
Ramsay et al., 2013). The results in the Kepler bandpass for the five objects
in our sample are given in the top panel, while the bottom panel shows the
Kepler data transformed into U-band energies and compared with previous
ground-based studies.

It is common to characterize the cumulative flare frequency as a power law in
energy, for energies above some limiting energy set by the completeness of the
sample. The low energy turnover seen most clearly in the active stars may be
partially a result of incompleteness in the detection of small flares and is
excluded from the fit. (In §4.1 below, we further discuss this turnover at low
energies, using the large sample of flares observed on GJ 1243.) GJ 1243 has a
very high flaring rate, comparable to those seen on the most active, well-
observed nearby flare stars such as AD Leo and YZ CMi. The power-law fit
(green solid line) to the GJ 1243 data shown in the top panel of Figure 8 was
obtained from a least squares fit to the high energy flares (using Poisson
uncertainties so that the very highest energy flares do not receive undue
weight), increasing the low energy limit until the slope of the power law did
not change within the uncertainty in the fit. This resulted in a limit of log
$E_{Kp}$ = 31, meaning that flares of this energy and above are used to fit
the power law prediction for flare frequency.

The GJ 1245AB flares are at somewhat lower energy, which is expected since the
stars are of slightly later spectral type, as discussed in LME. It is
surprising that the flare rate is reduced compared to GJ 1243 since flares
from both stars in the GJ 1245AB system are being measured. In addition, GJ
1245A is rotating significantly faster than GJ 1243. A simple activity-
rotation relation would predict stronger activity and more flares on GJ 1245A.
However, much of the light from GJ 1245A was excluded from the Kepler
photometric aperture, the threshold for detecting flares is higher because the
background flux comes from (parts of) both stars, and the signal-to-noise
ratio in the GJ 1245AB data is lower, as described in §2.1. All of these
factors add uncertainty to the analysis of the GJ 1245AB data from the
standard Kepler processing pipeline. While the combined data for GJ 1245AB
qualitatively indicate a somewhat steeper power law than for GJ 1243, we defer
quantitative discussion of the individual FFDs to Paper 3, where we
investigate this system in more detail, using individual light curves for each
component obtained from the Kepler pixel files.

The inactive stars show significantly lower flare frequency, but it is notable
that they still exhibit some energetic flares. Evidently, ”inactive” as
measured by the lack of H$\alpha$ emission does not mean that these stars have
no magnetic field. Rather, it appears that the field is weaker and/or less
organized into substantial active regions where magnetic heating can produce a
strong, persistent chromosphere. Despite the lack of a strong chromosphere,
magnetic energy is still being stored and released periodically in flares
which appear similar to the flares seen on the active stars. The FFDs for the
inactive stars suffer from poor sampling (small numbers), especially for the
mid-M star GJ 4083; linear least squares fits are shown, which are useful
mainly to characterize the general region populated by these stars in the FFD.
The power-law slopes and resulting values for $\alpha$ (the exponent in the
flare energy probability distribution, see Equation 1) for the four single
stars in our sample are given in Table 3.

The bottom panel of Figure 8 shows the Kepler data transformed into U-band
energies and the comparison with previous ground-based studies. The GJ 1243
FFD is illustrated with the green band; the inactive early type stars have
been combined into a single FFD (orange band, slope given in Table 3); and the
inactive mid-M star (GJ 4083) is shown with the blue band. The width of the
bands indicates the uncertainty due to the energy scaling as discussed in
§2.3. The dashed lines are power-law fits to hundreds of hours of monitoring
data in the U band on a large sample of M dwarfs, as reported in Hilton (2011,
see Figures 4.14 and 4.19). The light green dashed line corresponds to the
very active mid-M stars (including YZ CMi and AD Leo), the dark green dashed
line to a set of less active mid-M stars, the red dashed line to the inactive
early-M stars, and the blue dashed line to the inactive mid-M stars. The
slopes of these power-law fits to the ground-based data and resulting values
for $\alpha$ are also given in Table 3 for comparison.

The Kepler FFD for GJ 1243 falls between the very active and less active mid-M
relations. However, it is apparent that the Kepler data exhibit a steeper
power law, with fewer flares predicted at high energy. Steeper power law fits
to the FFD from Kepler data have also been found by Ramsay et al. (2013) who
compared results for two M dwarfs (including GJ 1243) with numerous individual
M dwarfs measured from the ground. One possible explanation is that the long
consecutive monitoring periods in the Kepler data provide better sampling of
the (rarer) high energy flares, and thus allow a more accurate estimate of the
power law fit in the high energy regime.

The results for the inactive early M stars in the Kepler sample, GJ 4099 (M1)
and GJ 4113 (M2) (orange band) exhibit flare frequencies similar to those
found by the extension of the ground-based relation (red dashed line) to
higher energies; the longer duration of continuous monitoring with Kepler
resulted in the detection of more flares at higher energy than previously were
seen on this class of stars. The small numbers in both the Kepler and ground-
based samples ($\sim$ 20 flares in each) mean that both fits suffer from
incompleteness, but it is reassuring that they occupy a similar region in the
diagram and plausibly result from measuring the same underlying distribution
using instruments with different sensitivity.

The inactive M3 star GJ 4083 is quite anomalous, with only 2 measured flares,
and a flare frequency distribution (blue band) more than a magnitude lower
than predicted by the ground-based relation for inactive mid-M stars (blue
dashed line) which itself only includes 3 flares, albeit in a much shorter
monitoring period. The flare rates on inactive mid-M stars apparently can vary
by almost two orders of magnitude, although the small numbers of flares
observed mean that this class of stars warrants further investigation. We note
that Paulson et al. (2006) and France et al. (2012) have observed
serendipitous flares on the inactive dM4 stars Gliese 699 and Gliese 876
respectively, although there are no measurements of the flare frequency
distributions for those stars. GJ 4083 is extremely quiet in its flaring rate
although it did show a small (erratic) signature of rotational modulation. It
may be a good candidate for inclusion in planet search programs that aim to
target stars with very low levels of magnetic activity.

Another measure that we can compute is the fractional energy released in
flares compared to the total energy output of the star. These ratios are
included in Table 2 as $f_{E}$, and clearly show that a higher fraction of the
stellar luminosity is released as flare energy in the active stars, in
agreement with Pettersen (1988) and Kowalski et al. (2009). The fraction of
time spent in a flaring state, $f_{flare}$, is also given in Table 2. The
active stars spend $\sim$30% of their time flaring at a level detectable with
Kepler, while the flaring time for inactive stars varies between
$\sim$0.01-1%. These fractions are higher than those found by Kowalski et al.
(2009) in their investigation of flares observed in SDSS imaging data,
probably due to the much higher flare threshold (0.7 mag enhancement in the
SDSS $u$ filter) adopted in that study.

In summary, our Kepler results broadly agree with ground-based data, although
we find a steeper power-law fit for the active mid-M star, perhaps due to
better sampling of high energy events. We suggest that the separation of stars
into only two groups (active and inactive) based on their H$\alpha$ emission
is inadequate to describe the full range of flare activity observed. The very
active and less active stars span more than a magnitude in flare frequency,
and the (relatively) inactive stars cover an additional 1-2 magnitudes at
lower frequency. This continuum in flaring activity across the entire range of
early-mid M dwarfs has important implications for modeling flare occurrence in
large surveys such as LSST (e.g. Hilton, 2011), and for predicting the effects
of flares on planets orbiting low mass stars.

Figure 7.— Simultaneous Kepler and U-band (from the NMSU 1m telescope at APO)
observations of a flare on GJ 1243. The start and stop times used to compute
the equivalent durations are indicated by vertical dashed lines.

Figure 8.— The top panel shows the cumulative flare frequency distributions,
together with least squares power-law fits, for the active and inactive M
stars in the Kepler sample, using the native Kepler flare energies. The bottom
panel uses the Kepler to U-band energy scaling from §2.2 to compare the Kepler
results with ground-based data. The mean power law relations for nearby M
dwarfs from Hilton (2011) are shown with dashed lines. See text for
discussion. Table 2Kepler Flare Energies and Statistics

Star | # Flares | Quiescent log $L_{Kp}$ | Range (log $E_{Kp}$) | Quarter | duration | $f_{flare}^{\dagger}$ | $f_{E}^{\star}$
---|---|---|---|---|---|---|---
GJ 4099 | 12 | 31.64 | 30.6–32.7 | Q8 | 42.9d | 0.015 | 7.2$\times 10^{-6}$
GJ 4113 | 8 | 31.53 | 31.1–32.1 | Q7 | 55.7d | 0.0056 | 1.9$\times 10^{-6}$
GJ 4083 | 2 | 31.18 | 30.7–31.2 | Q7 | 58.0d | 0.0005 | 2.7$\times 10^{-7}$
GJ 1243 | 833 | 30.67 | 28.3–33.1 | Q6 | 56.4d | 0.36 | 5.1$\times 10^{-3}$
GJ 1245AB | 450 | 30.22 | 29.1–32.5 | Q8 | 45.3d | 0.27 | 4.2$\times 10^{-3}$

$\dagger$$\dagger$footnotetext: $f_{flare}$ is the fraction of time spent in a
flaring state.

$\star$$\star$footnotetext: $f_{E}$ is the fraction of the star’s total energy
output released in flares.

Table 3Power law fits to FFDs from Kepler and the Ground

Star/Sp Types | source | power law slope | $\alpha$
---|---|---|---
GJ 4099 - M1 | Kepler† | -0.52 | 1.52
GJ 4113 - M2 | Kepler | -0.83 | 1.83
Inactive M1-M2 | Kepler | -0.72 | 1.72
Inactive M0–M2 | Ground⋆ | -1.06 | 2.06
GJ 4083 - M3 | Kepler | -0.67 | 1.67
Inactive M3–M5 | Ground | -0.25 | 1.25
GJ 1243 - M4 | Kepler | -1.01 | 2.01
Very Active M3–M5 | Ground | -0.53 | 1.53
Less Active M3–M5 | Ground | -0.67 | 1.67
GJ 1245AB - M5+M5 | Kepler | -1.32 | 2.32

$\dagger$$\dagger$footnotetext: Kepler results from this paper, see Figure 8.

$\star$$\star$footnotetext: Ground-based results from Hilton (2011), Table 4.3
(”Corrected” values for the slope $\beta$).

##  4\. Basic Flare Properties

The exquisite precision and very long duration of consecutive observations
provided by the Kepler data allow us to analyze some basic properties of
flares with unprecedented clarity. We have chosen to concentrate on the sample
of nearly 1000 flares from two months of monitoring on the single star GJ 1243
in our analysis. We examine correlations among the basic observed parameters -
amplitude, duration and energy - and discuss the completeness of the
observations at low energy. We then investigate correlations of rise time and
decay time with the flare duration.

###  4.1. Flare Amplitude, Duration and Energy

Figure 9.— The observed distributions of classical (black) and complex (red)
flares in amplitude, duration and log energy from two months of monitoring the
active M4 star GJ 1243 with Kepler. The full sample (classical and complex) is
shown for the energy diagram (light blue) and fit with a power law for
energies above log $E_{Kp}=31$. Figure 10.— Relationships between amplitude,
duration and energy for classical and complex flares on GJ 1243. Larger flare
energy is strongly correlated with longer duration, and is also correlated
with higher amplitude. Most high energy flares are complex, but there are
still some classical flares even near the highest energies. The amplitude and
duration relationship shows more scatter but still exhibits a significant
correlation, with longer flares typically having higher amplitude.

Figure 9 shows histograms of amplitude, duration and energy for the flares
classified as classical (only one peak) and complex (more than one peak). The
strong cutoff in the number of complex flares of short duration is not real,
but reflects the one minute time resolution of the Kepler data which makes it
difficult to confidently identify multiple peaks in flares shorter than $\sim$
15 minutes. The figure shows that classical and complex flares have similar
amplitude distributions but complex flares tend to last longer and be more
energetic.

If all flares exhibited an identical light curve, such as the classical flare
shown in Figure 6, then we would expect these measured parameters to be highly
correlated, subject only to measurement uncertainties. In fact, Figure 10
illustrates that both the classical and complex flares do show strong
correlations between amplitude, duration and energy. Note that the one-minute
cadence of the observations leads to the artificial integer binning of the
duration measurements; we expect that the durations would form a smooth
distribution as do the amplitudes and energies if measured with better time
resolution. In Paper 2, we analyze the light curve morphology using an even
larger (11 month) sample of classical flares on GJ 1243 and show that there is
a remarkable similarity in classical light curve shape that can be tied to
physical emission processes during flares. Briefly, the rise phase represents
a period of rapid impulsive heating characterized by strong blackbody-like
emission, while the decay phase first shows an initial rapid decline in the
blackbody-like emission followed by a longer slow decline representing gradual
cooling from enhanced Balmer and other chromospheric emission. The regions of
blackbody and Balmer emission may also be spatially distinct (Kowalski et al.,
2012). Kowalski et al. (2013) describes the spectroscopic evidence leading to
this interpretation, and the Kepler photometric data presented here and in
Paper 2 provide additional evidence that classical flares are well described
by such a model.

A power law fit to the energy histogram (above log $E_{Kp}$= 31, found from
the FFD, see Figure 8), recovers a slope of -0.95, very close to the slope in
the cumulative FFD (-1.01) as expected. The observations deviate from the
power law below log $E_{Kp}\sim 31$ which corresponds to duration about twenty
minutes and amplitude about 0.005 in fractional flux units. A typical flare of
this energy is indicated in Figure 2. However, Figure 2 also shows that even
flares of duration ten minutes (corresponding to amplitude 0.002 and log
$E_{Kp}\sim$ 30.5) are easily identified in the Kepler data, and we do not
believe we have missed a significant number of such flares, certainly not more
than half of them as would be suggested by Figure 9. One possibility is that
the lower energy flares do not follow the same power law distribution as the
high energy flares (for a discussion of the solar case, see Hannah et al.
(2011)). Alternatively, we may not always identify small flares as individual
events because they overlap in time with another flare and the entire event is
counted in our sample as a single complex flare of higher energy and longer
duration. Disentangling the individual flares in complex events is a difficult
problem with data that have no spatial resolution. In Paper 2, we present
results from a Monte Carlo simulation to investigate the incidence of
overlapping flares in the 11 month sample for GJ 1243.

The complex flares tend to have lower amplitude but longer duration than
classical flares at the same energy, and longer duration and higher energy in
general, in keeping with a model where they represent a superposition of
classical flares. The percentage of flares that are complex increases from 30%
for flares with duration between 20-30 minutes to 60% for flares with duration
more than an hour. Although we still see some high energy flares that exhibit
classical light curves, most high energy (long duration) flares are complex.
They may occur in active regions with complicated magnetic field structures
that are capable of storing more energy to be released as flares, and that are
therefore more prone to sequential heating of new flaring areas within the
same active region, analogous to the appearance of transient optical/UV foot-
points during flare evolution as seen in high spatial resolution solar
observations (Kosovichev & Zharkova, 2001; Wang et al., 2007). Sympathetic
flaring may also occur in nearby active regions, as has also been seen on the
Sun (Pearce & Harrison, 1990). However, longer duration flares are also more
likely to randomly overlap with another flare from an independent active
region. Paper 2 addresses the makeup of complex flares using an empirical
(classical) flare light curve to decompose the observed complex events, in
order to determine the relative importance of these two scenarios, which may
both play a role in forming flares with multiple peaks.

###  4.2. Flare Rise and Decay Times

Figure 11 shows the rise and decay times as a function of flare duration
(energy). The rise time is nearly always shorter than the decay time, and has
little correlation with the duration, while the decay time is strongly
correlated with the duration. The flare evolution typically includes a fast
rise phase and a longer decay phase, so it is natural that the decay phase
dominates the duration. Since the duration is also strongly correlated with
the flare energy, it is clear that most of the energy must be emitted by the
flare during the decay phase. Complex flares have longer rise times than
classical flares of the same duration, perhaps indicating that multiple events
are adding to produce the strongest peak in the light curve, where the rise
time is measured. The ratio of the rise and decay times shows a correlation
with flare duration such that longer duration flares spend relatively more
time in the decay phase, particularly for the classical flares. The ratio is
not constant, indicating that the rise phase and decay phase timescales do not
simply scale together in a linear fashion. In fact, the canonical flare
morphology that we find in Paper 2 is more complicated, with two distinct
decay phases (see discussion in §4.1). The longer duration flares are mostly
complex, and show little correlation between rise and decay timescales, likely
due to multiple heating events contributing to both phases.

Figure 11.— The flare rise time is not well correlated with the duration
(amplitude, energy) while the decay time is very strongly correlated. Flares
clearly spend the majority of their evolution in the decay phase, which
dominates the energy budget. Complex flares typically have longer rise times
at the same duration. The ratio of the rise and decay times decreases with
duration for classical flares, but complex flares, especially those of long
duration, show more complicated behavior (see text).

##  5\. Flare Timing

There are several aspects of flare timing that we can investigate with the
Kepler data on GJ1243, which has both numerous flares and well-defined
starspot modulation. We first look at the correlation of flare occurrence (and
energy) with starspot phase. Do flares preferentially occur when the star is
dimmest (largest surface coverage of dark spots on visible hemisphere)? Next,
we examine the waiting time (time between consecutive flares) distribution. Do
flares tend to cluster together, for example because a new active region has
appeared on the surface? Also, does the flare energy correlate with the energy
of the previous flare or the waiting time between flares? This might have
implications for the flare energy buildup and storage mechanisms.

###  5.1. Correlation with Starspot Modulation

Figure 12 shows the GJ 1243 light curve folded on the 0.59 day period
determined from the photometric modulation which we attribute to starspots on
the stellar surface. The flares appear as positive flux excursions above the
stable, periodic variation. The number of flares as a function of phase is
shown in the histogram and is consistent with no trend, with reduced
$\chi^{2}=1.1$. The flare energy is also randomly distributed, and shows no
significant correlation with starspot phase.

We are carrying out extensive models of the starspot distribution on GJ 1243
(Davenport et al, in prep.) which suggest that an asymmetric polar spot (or
more likely spot group) covering a significant fraction of the pole can
explain the primary feature of the observed rotational modulation in the light
curve. Since a large fraction of the polar spot group will always be in view,
we do not expect a strong correlation of flare occurrence with phase, if the
flares are associated with this spot group. Alternatively, the flares may be
coming from many, relatively small, active regions that contain most of the
magnetic flux from the star (Reiners & Basri, 2009), while the polar spot is
associated with the weaker dipolar field (Morin et al., 2008). These many
small spots would not contribute significant modulation compared to that from
the large polar spot group, and thus, again, we would not expect a correlation
of flare activity with starspot phase. We cannot distinguish between these
scenarios with our existing data.

Figure 12.— The top panel shows one month of GJ 1243 short cadence data folded
in phase on the starspot period (0.59 days). The red line indicates the median
value of $\Delta F/F$ at each phase, while the excursions above the median are
the flares. The middle panel illustrates that the number of flares per phase
bin is very flat, consistent with no trend with reduced $\chi^{2}=1.1$. The
bottom panel shows similarly that there is no correlation of flare energy with
phase.

###  5.2. Flare Waiting Time Distribution

The flare waiting time is often used in solar studies, and defined as the time
interval between the start of one flare and the start of the next flare. We
use the entire two-month sample of flares on GJ 1243 in our analysis,
including both classical and complex flares. However, note that complex flares
are counted as a single flare event, so the times between individual peaks in
complex flares are not considered in our analysis.

The top panel of Figure 13 shows the distribution of waiting times in the GJ
1243 sample. The decline in the number of flares observed with waiting times
below 30 minutes depends on the flare duration distribution, since a new
independent flare cannot be identified until the previous flare duration is
reached. Flares that occur while another flare is in progress are considered
part of a single complex flare. The bottom panel gives the cumulative
distribution, with the longest waiting time in the sample being almost nine
hours; a few flares have waiting times of five hours or more; and all flares
have waiting times more than a few minutes (a limit set by the time resolution
of the Kepler data and the duration of the shortest flares in the sample). The
cumulative distribution is well fit with a single exponential function (with
slope -0.375) for waiting times between 30 minutes and five hours, while the
extension of the fit to longer waiting times is adequate given the small
numbers. The data are consistent with a model where flares occur with a
single, random, Poisson distribution of waiting time for the entire extent of
our observations. In contrast, Wheatland (2010) investigated the evolution of
a single solar active region over a weeklong period, and found a piecewise
Poisson distribution of waiting times, with an initial flaring period having a
shallow slope (longer time between flares) followed by an active flaring
period with short waiting times and steeper slope, and then a decay period
again with longer waiting times and shallower slope.

Our interpretation is that GJ 1243 maintains a steady state of flaring over
long time periods (months) compared to both its rotation period (fraction of a
day) and compared to the behavior of active regions on the Sun (weeks). This
is consistent with a scenario where there are numerous small spot regions
(clustered in a large polar spot group, and/or spread out over the surface)
that are constantly emerging and decaying. At any one time, a large number of
regions are contributing to the overall flare rate so that while some may be
in the buildup or decay phases seen on the Sun, many others are in the
strongly flaring phase. The sum of all this flaring activity thus remains
relatively constant with time. We are investigating the flare frequency over
the entire eleven months of short cadence data on GJ 1243 in the Kepler
database and indeed it appears quite similar on monthly timescales (Wisniewski
et al., in prep.), which supports this steady state scenario of flare activity
on GJ 1243.

Figure 13.— Flare waiting times (time between successive flares) measured for
GJ 1243 are shown in the histogram (top panel) and as a cumulative
distribution (bottom panel). A single exponential fits the data well above a
wait time of 30 minutes. Flares at smaller wait times are underrepresented,
likely because many of them are counted as part of complex events.

###  5.3. Flare Energy Distribution

We define $\Delta$ log E as the energy difference between subsequent flares,
i.e. $\Delta$ log E = log Ei+1 \- log Ei, which is positive if a flare is more
energetic than the previous one. The distribution of $\Delta$ log E for the
full sample of GJ 1243 flares, including both classical and complex flares
with log E $>$ 31, is shown as the histogram in the top panel of Figure 14.
The flare energies appear to be randomly distributed, with no preferred scale
or sign for $\Delta$ log E. The smooth blue line represents a model
distribution of energy differences obtained by choosing flares randomly from a
power-law flare energy distribution with $\alpha=1.95$, found from the fit to
the energy histogram of all events in §4.1. The random model is consistent
with the observations of flares with log E $>$ 31, where the power law is
applicable.

The distribution of $\Delta$ log E with waiting time is shown in the bottom
panel of Figure 14. Again there is no correlation of waiting time with the
energy of the subsequent flare. We do expect, and observe, a region in the
lower left of the diagram at negative $\Delta$ log E (first flare is higher
energy) and small waiting time that is relatively unpopulated. This is because
larger flares typically have longer duration, and hence necessarily a longer
waiting time before the next flare can be unambiguously identified.

These results support a scenario where the individual flares we have measured
are randomly occurring at times and energies that do not depend on the
previous flare history. Thus, each event may be occurring in a separate active
region, not triggered by previous events. We do not find any evidence of
sympathetic flaring between individual events, as suggested for example by
Panagi & Andrews (1995). We stress again that we treat complex flares as
single events in this analysis. Within a complex flare, the sub-peaks that we
see may sometimes be triggered from the initial reconnection and energy
release, see e.g. Anfinogentov et al. (2013). We discuss the makeup of complex
flares in more detail in Paper 2.

Figure 14.— The top panel shows the distribution of $\Delta$ log E between
successive flares (black histogram). The model (blue line) gives the expected
distribution if flares are drawn randomly from the power law fit to the GJ
1243 flare energies, and is consistent with the observed distribution. The
bottom panel shows that the energy difference between successive flares is not
correlated with the waiting time. Only flares with energy log $E_{Kp}$ $>$ 31
are included.

##  6\. Summary

We have analyzed two months of Kepler short-cadence data on several low mass M
dwarfs, spanning a range of activity and spectral type. We obtained additional
spectroscopic data from the ground to assign spectral types, and used
simultaneous Kepler and U-band photometry of one flare to obtain an energy-
scaling relationship between Kepler and U-band flare energies. We report
periods for five stars; these agree with published periods for the two stars
that have previous determinations.

The flare frequency distributions (FFDs) from Kepler data probe higher flare
energies with better statistics than are available from the ground. We suggest
that this leads to steeper power-law fits to the FFDs compared to those found
from ground-based data, indicating that the power law energy distribution of
flares may approach $\alpha=2$, with implications for coronal heating.
However, we also find evidence that the power law fit does not continue to
lower energies, which are still above our detection limit. The missing low
energy flares may be subsumed into complex events, or there may be a real
turnover in the power law energy distribution of flares. We investigate these
possibilities further in Paper 2 of this series.

Comparison of the FFDs for inactive and active stars of early and mid-M
spectral type shows that there is a continuum of flare frequencies and
energies, and that the artificial separation into ”active” and ”inactive”
classes based solely on the presence or absence of H$\alpha$ emission in the
quiescent spectrum does not capture the full complexity of magnetic (flaring)
activity in M dwarfs. Flares from inactive stars, though relatively rare, can
be of high energy, and the large number of early-type, inactive M dwarfs in
the Galaxy means that these flares could be an important source of transient
phenomena in time domain surveys. In addition, flares from inactive stars may
still play an important role in the atmospheric chemistry of orbiting planets.

We also examined the properties of basic flare parameters for a sample of
nearly 1000 flares on the active M4 star GJ 1243, and found that flare
amplitude, duration and energy are all strongly correlated. High energy flares
typically have both long duration and high amplitude. Complex flares (with
more than one peak) have longer duration and higher energy at the same
amplitude. The fraction of flares that are complex also rises with flare
duration. These results support a scenario where complex flares represent the
superposition of classical flares, possibly from flares in separate regions
that randomly overlap in time and/or from triggered flaring within the same or
nearby active regions. Paper 2 continues the exploration of the makeup of
complex flares using a classical flare model to decompose the complex flare
light curves.

The precision and long monitoring period of the Kepler data also allowed us to
investigate the correlation of flare occurrence and energy with starspot phase
on GJ 1243. We found no significant correlation, indicating that the spot or
spot group that contributes to the periodic photometric modulation does not
represent the only, or even the primary, spatial location where flares occur.
Instead, the random occurrence of flares with phase indicates that flaring is
distributed across the stellar surface as it rotates, possibly in many small
active regions. These results are consistent with the large scale magnetic
structure model of Morin et al. (2008), which accounts for only a small
fraction of the total magnetic energy in active mid-M stars; most of the
energy is likely found in a large number of smaller active regions that are
not detectable in global polarization measurements but contribute to the large
magnetic field measurements of several kG with large filling factor (Johns-
Krull & Valenti, 1996).

Finally, we considered the flare waiting time distribution, and found that it
is well described with a single exponential, indicating a random Poisson
distribution of flare waiting times above 30 minutes. The evolution of a
single active region on the Sun shows a piecewise Poisson distribution, with
periods of weak flaring as the region builds and decays, and strong flaring in
between (Wheatland, 2010). We interpret the steady-state flaring distribution
that persists for the entire two months on GJ 1243 as additional evidence for
the contribution of many active regions in a strongly flaring state, such that
a few regions that are in the buildup and decay phases do not contribute
significantly to the overall rate. We do not find evidence for sympathetic
flaring, but note that complex events are treated as single flares in our
analysis. Sympathetic flaring within a given complex flare may still be
occurring.

The difference in energy between successive flares was also investigated, and
the results were consistent with a model where flares are randomly chosen from
a power-law energy distribution with $\alpha$ appropriate for GJ 1243. This
does not support the notion of independent precursor flares, where a large
flare is immediately preceded by a smaller flare. The waiting time between
flares also showed no correlation with the energy of the next flare,
indicating that the flares are occurring in independent regions, or that the
magnetic energy stored in a given region is much larger than is being emitted
in flares.

In summary, our data support a model where many small-scale active regions are
distributed across the stellar surface, contributing only marginally to the
rotational modulation of the star, and flaring independently at random
intervals, with no connection in time or energy between successive flares.

Additional studies using our Kepler short cadence data will investigate the
temporal morphology of the classical flare light curve and the makeup of
complex flares (Paper 2); the flaring behavior of the individual components of
the GJ 1245AB system (Paper 3); the flare frequency and multi-wavelength
behavior of flares in the full 11-month Kepler dataset, together with ground-
based data, on GJ 1243 (Wisniewski et al., in prep.); and the comparison of
flares on stars of types G, K and M (Hawley et al., in prep.). The Kepler
satellite has delivered stellar light curves of unprecedented precision and
duration and will provide a rich data source for flare studies for the
foreseeable future.

This work was supported by Kepler Cycle 2 GO grant NNX11AB71G and Cycle 3 GO
grant NNX12AC79G . JRAD acknowledges support from NASA ADP grant NNX09AC77G.
SLH, JRAD and LH acknowledge support from NSF grant AST13-11678. EJH, AFK, and
SLH acknowledge support from NSF grant AST08-07205. Observations reported here
were obtained with the Apache Point Observatory 3.5-meter telescope, which is
owned and operated by the Astrophysical Research Consortium. This paper
includes data collected by the Kepler mission. Funding for the Kepler mission
is provided by the NASA Science Mission directorate. Some of the data
presented in this paper were obtained from the Mikulski Archive for Space
Telescopes (MAST). STScI is operated by the Association of Universities for
Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST
for non-HST data is provided by the NASA Office of Space Science via grant
NNX13AC07G and by other grants and contracts.

## References

  * Anfinogentov et al. (2013) Anfinogentov, S., Nakariakov, V. M., Mathioudakis, M., Van Doorsselaere, T., & Kowalski, A. F. 2013, ApJ, 773, 156
  * Audard et al. (2000) Audard, M., Güdel, M., Drake, J. J., & Kashyap, V. L. 2000, ApJ, 541, 396
  * Bochanski et al. (2010) Bochanski, J. J., Hawley, S. L., Covey, K. R., West, A. A., Reid, I. N., Golimowski, D. A., & Ivezić, Ž. 2010, AJ, 139, 2679
  * Bochanski et al. (2012) —. 2012, AJ, 143, 152
  * Borucki et al. (2010) Borucki, W. J., et al. 2010, Science, 327, 977
  * Covey et al. (2007) Covey, K. R., et al. 2007, AJ, 134, 2398
  * Davenport et al. (2012) Davenport, J. R. A., Becker, A. C., Kowalski, A. F., Hawley, S. L., Schmidt, S. J., Hilton, E. J., Sesar, B., & Cutri, R. 2012, ApJ, 748, 58
  * Davenport et al. (2014) Davenport, J. R. A., Hawley, S. L., Hebb, L., Wisniewski, J. P., Kowalski, A. F., & et al. 2014, ApJ, submitted
  * Davenport et al. (2006) Davenport, J. R. A., West, A. A., Matthiesen, C. K., Schmieding, M., & Kobelski, A. 2006, PASP, 118, 1679
  * Donati et al. (2008) Donati, J.-F., et al. 2008, MNRAS, 390, 545
  * Doyle et al. (1990) Doyle, J. G., Butler, C. J., van den Oord, G. H. J., & Kiang, T. 1990, A&A, 232, 83
  * France et al. (2012) France, K., Linsky, J. L., Tian, F., Froning, C. S., & Roberge, A. 2012, ApJ, 750, L32
  * Friedman (1984) Friedman, J. H. 1984, A Variable Span Smoother, Tech. Rep. 5, Department of Statistics, Stanford University
  * Gershberg (1972) Gershberg, R. E. 1972, Ap&SS, 19, 75
  * Gizis et al. (2013) Gizis, J. E., Burgasser, A. J., Berger, E., Williams, P. K. G., Vrba, F. J., Cruz, K. L., & Metchev, S. 2013, ApJ, 779, 172
  * Güdel (2004) Güdel, M. 2004, A&A Rev., 12, 71
  * Hannah et al. (2011) Hannah, I. G., Hudson, H. S., Battaglia, M., Christe, S., Kašparová, J., Krucker, S., Kundu, M. R., & Veronig, A. 2011, Space Sci. Rev., 159, 263
  * Hawley & Pettersen (1991) Hawley, S. L., & Pettersen, B. R. 1991, ApJ, 378, 725
  * Hilton (2011) Hilton, E. J. 2011, PhD thesis, University of Washington
  * Hudson (1991) Hudson, H. S. 1991, Sol. Phys., 133, 357
  * Hunt-Walker et al. (2012) Hunt-Walker, N. M., Hilton, E. J., Kowalski, A. F., Hawley, S. L., & Matthews, J. M. 2012, PASP, 124, 545
  * Irwin et al. (2011) Irwin, J., Berta, Z. K., Burke, C. J., Charbonneau, D., Nutzman, P., West, A. A., & Falco, E. E. 2011, ApJ, 727, 56
  * Johns-Krull & Valenti (1996) Johns-Krull, C. M., & Valenti, J. A. 1996, ApJ, 459, L95
  * Kosovichev & Zharkova (2001) Kosovichev, A. G., & Zharkova, V. V. 2001, ApJ, 550, L105
  * 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
  * Kowalski et al. (2012) Kowalski, A. F., Hawley, S. L., Holtzman, J. A., Wisniewski, J. P., & Hilton, E. J. 2012, Sol. Phys., 277, 21
  * Kowalski et al. (2013) Kowalski, A. F., Hawley, S. L., Wisniewski, J. P., Osten, R. A., Hilton, E. J., Holtzman, J. A., Schmidt, S. J., & Davenport, J. R. A. 2013, ApJS, 207, 15
  * Lacy et al. (1976) Lacy, C. H., Moffett, T. J., & Evans, D. S. 1976, ApJS, 30, 85
  * McQuillan et al. (2013) McQuillan, A., Aigrain, S., & Mazeh, T. 2013, MNRAS, 432, 1203
  * Moffett (1974) Moffett, T. J. 1974, ApJS, 29, 1
  * Morin et al. (2008) Morin, J., et al. 2008, MNRAS, 390, 567
  * Morin et al. (2010) Morin, J., Donati, J.-F., Petit, P., Delfosse, X., Forveille, T., & Jardine, M. M. 2010, MNRAS, 407, 2269
  * Osten & Brown (1999) Osten, R. A., & Brown, A. 1999, ApJ, 515, 746
  * Panagi & Andrews (1995) Panagi, P. M., & Andrews, A. D. 1995, MNRAS, 277, 423
  * Parker (1988) Parker, E. N. 1988, ApJ, 330, 474
  * Paulson et al. (2006) Paulson, D. B., Allred, J. C., Anderson, R. B., Hawley, S. L., Cochran, W. D., & Yelda, S. 2006, PASP, 118, 227
  * Pazzani & Rodono (1981) Pazzani, V., & Rodono, M. 1981, Ap&SS, 77, 347
  * Pearce & Harrison (1990) Pearce, G., & Harrison, R. A. 1990, A&A, 228, 513
  * Pettersen (1988) Pettersen, B. R. 1988, in Astrophysics and Space Science Library, Vol. 143, Activity in Cool Star Envelopes, ed. O. Havnes, J. E. Solheim, B. R. Pettersen, & J. H. M. M. Schmitt , 49–60
  * Ramsay et al. (2013) Ramsay, G., Doyle, J. G., Hakala, P., Garcia-Alvarez, D., Brooks, A., Barclay, T., & Still, M. 2013, MNRAS, 434, 2451
  * Reid et al. (1995) Reid, I. N., Hawley, S. L., & Gizis, J. E. 1995, AJ, 110, 1838
  * Reiners & Basri (2009) Reiners, A., & Basri, G. 2009, A&A, 496, 787
  * Robinson et al. (1995) Robinson, R. D., Carpenter, K. G., Percival, J. W., & Bookbinder, J. A. 1995, ApJ, 451, 795
  * Robinson et al. (2001) Robinson, R. D., Linsky, J. L., Woodgate, B. E., & Timothy, J. G. 2001, ApJ, 554, 368
  * Schrijver et al. (2012) Schrijver, C. J., et al. 2012, Journal of Geophysical Research (Space Physics), 117, 8103
  * Schroeder et al. (2000) Schroeder, D. J., et al. 2000, AJ, 119, 906
  * Smith et al. (2012) Smith, J. C., et al. 2012, PASP, 124, 1000
  * Stassun et al. (2011) Stassun, K. G., et al. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West, 505
  * Walkowicz et al. (2011) Walkowicz, L. M., et al. 2011, AJ, 141, 50
  * Wang et al. (2007) Wang, L., Fang, C., & Ming-DeDing. 2007, Chinese J. Astron. Astrophys., 7, 721
  * Wheatland (2010) Wheatland, M. S. 2010, ApJ, 710, 1324


[a]Felix Riehn

# Sibyll★: ad-hoc modifications for an improved description of muon data in
extensive air showers

Ralph Engel  Anatoli Fedynitch

###### Abstract

Current simulations of air showers produced by ultra-high energy cosmic rays
(UHECRs) do not satisfactorily describe recent experimental data, particularly
when looking at the muonic shower component relative to the electromagnetic
one. Discrepancies can be seen in both average values and on an individual
shower-by-shower basis. It is thought that the muonic part of the air showers
isn’t accurately represented in simulations, despite various attempts to boost
the number of muons within standard hadronic interaction physics. In this
study, we investigate whether modifying the final state of events created with
Sibyll 2.3d in air shower simulations can achieve a more consistent
description of the muon content observed in experimental data. We create
several scenarios where we separately increase the production of baryons,
$\rho^{0}$, and strange particles to examine their impact on realistic air
shower simulations. Our results suggest that these ad-hoc modifications can
improve the simulations, providing a closer match to the observed muon content
in air showers. One side-effect of the increased muon production in the
considered model versions is a smaller difference in the predicted total muon
numbers for proton and iron showers. However, more research is needed to find
out whether any of these adjustments offers a realistic solution to the
mismatches seen in data, and to identify the precise physical process causing
these changes in the model. We hope that these modified model versions will
also help to develop improved machine-learning analyses of air shower data and
to estimate sys. uncertainties related to shortcomings of hadronic interaction
models.

##  1 Introduction

This study addresses the ’Muon Puzzle’ [1, 2], a key challenge in the
interpretation of extensive air shower (EAS) data. Specifically, we focus on
the discrepancy between observed and predicted muon numbers in EAS, along with
the associated inconsistency between the depth of shower maximum ($\langle
X_{\rm max}\rangle$) and the ground-level signal [3, 4]. In pursuit of gaining
better understanding of this problem, we explore ad-hoc modifications of the
hadronic interaction model Sibyll 2.3d [5, 6, 7]. While numerous standard [8,
9, 10, 11, 12] and exotic [13, 14, 15, 16, 17] mechanisms have been proposed
to address this problem, our study specifically focuses on quantitatively
evaluating several conventional mechanisms. For this purpose, we employ a
customized version of Sibyll, designed for direct use in realistic EAS
simulations. All modifications of the model satisfy the constraint that, on an
event-by-event level, all fundamental constraints such as energy-momentum and
quantum numbers of relevance are conserved. Furthermore, existing experimental
data from particle physics experiments are considered as guiding input,
limiting some of the changes to energy and phase space regions for which no
such data exist.

##  2 Sibyll★

We construct the custom models by modifying events generated with Sibyll 2.3d.
Once the initial event generation is complete, we start by letting all
hadronic resonances with a shorter lifetime than that of K${}^{0}_{\rm s}$
decay, except for $\pi^{0}$. We then go through the event’s particle list to
identify appropriate pairs or triples of pions, only considering the five
nearest neighbors for every pion. If the combined invariant mass of the
selected pions is sufficient and the sampling criterion is fulfilled, we
exchange them with a pair of new particles that together have the same total
momentum, invariant mass, and charge. We calculate the final momenta from the
invariant mass, the new particles’ masses, and a minor transverse momentum
sourced from an exponential distribution in transverse mass. In these altered
events, we maintain conservation of energy, momentum, and charge, although
(iso)spin conservation is not upheld.

The acceptance rate of these particle exchanges hinges on the total center-of-
mass (CM) energy ($\sqrt{s}$) and the fraction of longitudinal momentum
$x_{\rm F}$ (defined as $p_{z}/p_{z,\mathrm{max}}$, with momenta $p$ in the CM
frame) of the initial particle. The probability of exchanging particles is
parameterized by

$P_{i}~{}=~{}P_{i,0}\cdot|x_{\rm F}|^{\epsilon_{i}}\cdot f(\sqrt{s},E_{\rm
thr})\ .$ (1)

The emphasis given to forward or central particles depends on the chosen value
for the exponent $\epsilon_{i}$ in the $x_{\rm F}$-dependence. If
$\epsilon_{i}=0$, all particles receive equal weight, preserving the original
distribution’s shape in longitudinal phase space. As $\epsilon_{i}$ approaches
1, the forward part of the $x_{\rm F}$-spectrum is enhanced.

The energy dependence of $f(\sqrt{s},E_{\rm thr})$ is logarithmic. It is set
such that the rate is precisely zero below a threshold energy $E_{\rm thr}$
and reaches the nominal $P_{0,i}$ at lab energies of $10^{19}\,$eV
($1.37\times 10^{5}\,$GeV in CM frame). The threshold energy $E_{\rm thr}$ is
set at $5\,$GeV.

This parameterization of energy dependence represents a very gradual change in
particle production, from no change at low energies, where fixed target
experiments effectively limit the entire phase space, to LHC energies, where
only the central region is well constrained, up to the UHECR energy scale
(with no lab experiment constraints). Overall, this spans five orders of
magnitude in energy, and the modification scales from zero to one. As an
alternative, we allow a more drastic increase of the exchange rate towards
high energies (above $13\,$TeV in CM frame) represented as $P_{i}\to
P_{i,0}+P_{i,\mathrm{HE}}\cdot f_{\rm HE}(\sqrt{s},13\,\mathrm{TeV})$. This
mode depicts a swift transition to new physics beyond the LHC scale. Lastly,
we apply this algorithm for all projectiles or only for mesons.

Figure 1: Fraction of projectile energy carried by $\rho^{0}$, anti protons
and charged kaons in $\pi^{-}$C collisions [18, 19]. Lines are Sibyll 2.3d,
Sibyll 2.1 and different variants of Sibyll★. Figure 2: Multiplicities of
$\pi^{+}$, anti protons and charged kaons in pp collisions [20, 21, 22]. Lines
are Sibyll 2.3d, Sibyll 2.1 and different variants of Sibyll★.

Using the previously discussed algorithm, we construct different Sibyll 2.3d
variants with an aim to enhance muon production in EAS. We select three
distinct modifications known for their efficacy in muon production: $\rho^{0}$
production, baryon anti-baryon pair-production, and kaon production
enhancement [8, 9, 10, 11, 17]. These variants are denoted as S★($\rho^{0}$),
S★($\bar{p}$), and S★(K±,0).

In the $\rho^{0}$ variant, $\pi^{0}$ are directly substituted with $\rho^{0}$.
For the baryon pair and kaon pair variant, charge-neutral combinations of two
or three pions are replaced with $p\bar{p}$ or $n\bar{n}$ pairs, and K+K- or
K0 $\bar{\mathrm{K}}^{0}$ pairs respectively. We adjust the parameters for
each variant to align the energy in the corresponding component with NA61
measurements [19, 18] (refer to Fig.1 and Fig. 2).

Producing $\rho^{0}$ significantly impacts muon production in EAS because
$\rho$ mesons can form directly from the pion projectile. This mechanism
doesn’t hold for proton projectiles, hence for the $\rho^{0}$ variant, only
meson projectile interactions are modified. However, the modifications for
baryon pair and kaon production apply to any projectile and include a rapid
increase at energies beyond the LHC. The parameters are adjusted so that the
energy fraction carried by all hadrons except neutral pions at $10\,$EeV is
approximately the same across all variants ($\approx 0.82$).

Furthermore, we create a fourth variant using both $\rho^{0}$ and baryon pair
production (S★(mix)). In this scenario, we choose a more moderate increase of
$\rho$ production and dismiss the rapid increase of baryon production at high
energies. The parameters of the different variants are detailed in Tab. 1.

Table 1: Parameters in different variants of Sibyll★.  Label | $P_{i,0}$ | forward weight $\epsilon_{i}$ | projectiles | $P_{i,\mathrm{HE}}$
---|---|---|---|---
$\rho^{0}$ | 0.8 | 0.3 | mesons | -
$\bar{p}$ | 0.5 | 0.7 | all | 0.25
$K^{\pm,0}$ | 0.5 | 0.8 | all | 0.3
$\rho$-mix | 0.8 | 0.4 | mesons | -
$\bar{p}$-mix | 0.5 | 0.7 | all | -

Note that the mechanism driving the effectiveness of both strangeness
enhancement and increased baryon pair production in raising the muon number is
fundamentally the same: quantum number conservation (strangeness and baryon
number, respectively). The key difference is that baryon production works
effectively at all energies (no nucleon decay), while strangeness is only
conserved in EAS at high energies, where kaon decay is negligible. Note also
that since hyperons do not decay into kaons, enhanced hyperon production is
part of the variant with enhanced baryon production.

##  3 EAS predictions

Figure 3: $N_{\mu}$ and $X_{\rm max}$ for proton showers at $67^{\circ}$
across Sibyll variants. The left panel shows a substantial increase, up to
35%, in the number of muons for the mixed and $\rho$ variant. However, the
variation on the shower maximum between Sibyll 2.3d and its variants is less
than $7\,$g$/$cm2. The grey line represents the required increase in muon
count to align with the data from the Pierre Auger Observatory.

In our EAS predictions, we carried out full air shower simulations using
CORSIKA111CORSIKA v7.7420 [23]; The magnetic field and observation level are
set to the values of the site of the Pierre Auger Observatory. Low-energy
hadronic interactions were modeled using FLUKA 2021.2.9 [24]. for each S★
variant. We compared the resulting average depth of shower maximum ($\langle
X_{\rm max}\rangle$) and the average number of muons at ground level ($\langle
N_{\mu}\rangle$) with simulations using the unmodified Sibyll 2.3d. These
simulations were set to mirror conditions at the Pierre Auger Observatory
site. The results for a primary proton with energy $10\,$EeV and an incident
zenith angle of $67^{\circ}$ are presented in Fig. 3. These results show a
notable increase in the number of muons, while $\langle X_{\rm max}\rangle$
remains largely unaffected.

Figure 4: Comparison of the predicted values for $\langle X_{\rm max}\rangle$
and $\langle N_{\mu}\rangle$ from various Sibyll★ variants to the measurements
obtained from the Pierre Auger Observatory. [25]

In the analysis shown in Fig. 4, we compare the predicted values of $\langle
X_{\rm max}\rangle$ and $\langle\ln N_{\mu}\rangle$ from proton, helium,
nitrogen, and iron primaries with the measurement from the Pierre Auger
Observatory [25]. Our findings indicate that only the $\rho^{0}$ and the mixed
$\rho^{0}$ baryon-pair variant yield a sufficient number of muons to match the
levels shown in the data.

The model predictions in Fig. 4 for p, He, N and Fe all fall on a line. The
reason is that per the superposition model, $\langle X_{\rm max}\rangle$ and
$\langle\ln N_{\mu}\rangle$, have the same dependence on the primary mass
(linear in $\ln A$), e.g. $\langle\ln N_{\mu}\rangle(A,E)~{}=~{}(1-\beta)\ln
A\,+\,\langle\ln N_{\mu}\rangle(1,E)$, where $\beta$ is the exponent in the
energy dependence of the number of muons for proton primaries, that is
$\langle\ln N_{\mu}\rangle(1,E)~{}=~{}\beta\ln E$. In simplified models à la
Heitler-Matthews we have

$N_{\mu}=A\cdot\left(\frac{E}{A\cdot E_{\rm dec}}\right)^{\beta}$ (2)

with $\beta$ being computed as $\ln N_{\rm ch}/\ln N_{\rm tot}$, where $N_{\rm
ch}$ and $N_{\rm tot}$ denote the multiplicities of charged and all pions,
respectively [26]. More broadly, $\beta$ is associated with the fraction of
hadrons that carry enough kinetic energy to re-interact rather than decay. The
slopes in Fig. 4 decreases from Sibyll 2.1 towards the variants with the
highest levels of muon production. This behavior is common to all variants S★.
The greater a fraction of energy is kept in hadrons, the larger the increase
in the number of muons and the larger $\beta$, and consequently, the less
pronounced is the separation of the primary masses in $\langle
N_{\mu}\rangle$.

Figure 5:  Ratio between the average number of muons in iron and proton
showers (left panel) and the slope $\beta$ (right panel). Both are shown as a
function of the distance from the shower axis. The inset numbers in the left
panel show the ratio of the total numbers of muons between iron and proton
showers and the ratio of the muon densities at $1000\,$m in parentheses.

The data clearly favor a larger muon content at high energy leading to a
reduction in the mass resolution. However Fig. 5 shows that for specific
experiments the situation may be not as dramatic, as the separation between
primary masses (respectively $\beta$) varies with the distance from the shower
axis. The reason is that muons at different lateral distances are dominated by
different phases of the shower development [6, 27, 28, 29, 30].

##  4 Inclusive fluxes

We employed the MCEq code [31] for a comparative analysis of the atmospheric
muon and neutrino fluxes as predicted by the S★ versions against the original
Sibyll 2.3d. Considering the modification (S${}^{\bigstar}(\rho^{0})$) impacts
secondary pion interactions and (S${}^{\bigstar}(\bar{p})$) affects the
production of secondary baryons, both of which influence the shower
development in deeper atmospheric layers, we didn’t foresee changes to
inclusive fluxes, which primarily depend on the air shower’s early stages. Our
findings confirm no significant effects. However, since $K^{\pm}$ decays into
muons and neutrinos, for the S★(K${{}^{\pm},0}$) model, we identified a minor
increase in muon fluxes of around 5% at tens of TeV and PeV energies, and up
to a 20% increase in atmospheric neutrino fluxes.

##  5 Discussion

In this study, we explored custom variants of the Sibyll 2.3d model, aiming to
boost muon production in extensive air showers. Our focus centered on three
modifications: $\rho^{0}$ production, baryon anti-baryon pair-production, and
kaon production enhancement. Significantly, the $\rho^{0}$ and mixed
$\rho^{0}$ baryon-pair enhancements effectively aligned with the observed muon
production data from the Pierre Auger Observatory at $10\,$EeV. In contrast,
elevating strangeness or baryon production proved insufficient, even when
amplified to extreme levels. However, our results don’t definitively eliminate
these scenarios, especially considering that our implementation does not
permit the production of leading strangeness (e.g.
$p\to\Lambda^{0}\to\mathrm{K}^{+}\,+\,$n).

Our simulations show that these modifications increase the number of muons,
while the average depth of shower maximum remains largely unaffected. Another
observation is that the degree of primary mass separation through muon
measurement is predicted to be smaller in the modified models than in the
original Sibyll 2.3d and other interaction models as along as the total number
of muons is concerned. This is also expected within the Heitler-Matthews model
[26] as more energy is kept in the hadronic channel in each interaction,
increasing the exponent $\beta$. On the other hand, the muon density at
sufficient distance from the shower core (for example, 1000 m) is a very well
suited observable for composition studies.

Moreover, despite certain modifications leading to increased muon fluxes,
inclusive fluxes - primarily dependent on early air shower stages - remained
largely unchanged. Only the S★(K${{}^{\pm},0}$) model showed a minor increase
in muon fluxes and a noticeable increase in atmospheric neutrino fluxes, due
to the decay of $K^{\pm}$ into muons and neutrinos.

The provided variants of Sibyll 2.3d can be used to train machine-learning
models like deep neural networks to have a better description of the
measurements and to estimate systematic uncertainties stemming from
shortcomings of modeling hadronic multiparticle production.

### Acknowledgements

The authors acknowledge many fruitful discussions with colleagues of the
Pierre Auger and IceCube Collaborations. They have used computing resources
provided by the Academia Sinica Grid Computing Center (ASGC), supported by the
Institute of Physics at Academia Sinica. FR and RE are supported in part by
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No. 101065027 and by BMBF grant No.
05A2023VK4, respectively.

## References

  * [1] J. Albrecht et al., _The Muon Puzzle in cosmic-ray induced air showers and its connection to the Large Hadron Collider_ , _Astrophys. Space Sci._ 367 (2022) 27 [2105.06148].
  * [2] J.C. Arteaga-Velzquez, _An update on the combined analysis of muon data from nine air-shower experiments at cosmic-ray energies above 1 PeV_ , _PoS_ ICRC (2023) 466.
  * [3] Pierre Auger collaboration, _Deep-learning based reconstruction of the shower maximum $X_{max}$ using the water-Cherenkov detectors of the Pierre Auger Observatory_, _JINST_ 16 (2021) P07019 [2101.02946].
  * [4] Pierre Auger collaboration, _Probing hadronic interaction models with the hybrid data of the Pierre Auger Observatory_ , in _21st International Symposium on Very High Energy Cosmic Ray Interactions_ , 9, 2022 [2209.00744].
  * [5] E.-J. Ahn, R. Engel, T.K. Gaisser, P. Lipari and T. Stanev, _Cosmic ray interaction event generator SIBYLL $2.1$_, _Phys. Rev. D_ 80 (2009) 094003.
  * [6] F. Riehn, R. Engel, A. Fedynitch, T.K. Gaisser and T. Stanev, _Hadronic interaction model SIBYLL 2.3d and extensive air showers_ , _Phys. Rev. D_ 102 (2020) 063002.
  * [7] A. Fedynitch, F. Riehn, R. Engel, T.K. Gaisser and T. Stanev, _Hadronic interaction model SIBYLL $2.3\mathrm{c}$ and inclusive lepton fluxes_, _Phys. Rev. D_ 100 (2019) 103018.
  * [8] S. Ostapchenko, _QGSJET-II: physics, recent improvements, and results for air showers_ , _EPJ Web Conf._ 52 (2013) 02001.
  * [9] H.-J. Drescher, _Remnant break-up and muon production in cosmic ray air showers_ , _Phys. Rev. D_ 77 (2007) 056003.
  * [10] T. Pierog and K. Werner, _Muon production in extended air shower simulations_ , _Phys. Rev. Lett._ 101 (2008) 171101.
  * [11] P.K.F. Grieder, _The effect of $n\bar{n}$-production on particle spectra in vertically incident and inclined showers derived from simulations_, _Proc. of 13th Int. Cosmic Ray Conf., Denver_ 4 (1973) 2467\.
  * [12] S. Baur, H. Dembinski, M. Perlin, T. Pierog, R. Ulrich and K. Werner, _Core-corona effect in hadron collisions and muon production in air showers_ , _Phys. Rev. D_ 107 (2023) 094031 [1902.09265].
  * [13] G.R. Farrar and J.D. Allen, _A new physical phenomenon in ultra-high energy collisions_ , _EPJ Web Conf._ 53 (2013) 07007 [1307.2322].
  * [14] J. Alvarez-Muniz, L. Cazon, R. Conceicao, J.D. de Deus, C. Pajares and M. Pimenta, _Muon production and string percolation effects in cosmic rays at the highest energies_ , 1209.6474.
  * [15] M. Rybczyński and Z. Włodarczyk, _Puzzle of muons in extensive air showers_ , _Int. J. Mod. Phys. D_ 28 (2019) 1950097 [1902.09988].
  * [16] L.A. Anchordoqui, H. Goldberg and T.J. Weiler, _Strange fireball as an explanation of the muon excess in Auger data_ , _Phys. Rev. D_ 95 (2017) 063005 [1612.07328].
  * [17] L.A. Anchordoqui, C.G. Canal, F. Kling, S.J. Sciutto and J.F. Soriano, _An explanation of the muon puzzle of ultrahigh-energy cosmic rays and the role of the Forward Physics Facility for model improvement_ , _JHEAp_ 34 (2022) 19 [2202.03095].
  * [18] NA61/SHINE collaboration, _Measurement of meson resonance production in $\pi^{-}+$ C interactions at SPS energies_, _Eur. Phys. J._ C77 (2017) 626 [1705.08206].
  * [19] NA61/SHINE collaboration, _Measurement of hadron production in $\pi$–C interactions at 158 and 350 GeV/c with NA61/SHINE at the CERN SPS_, _Phys. Rev. D_ 107 (2023) 062004 [2209.10561].
  * [20] CMS collaboration, _Measurement of charged pion, kaon, and proton production in proton-proton collisions at sqrt(s) = 13 TeV_ , 1706.10194.
  * [21] CMS Collaboration collaboration, _Study of the inclusive production of charged pions, kaons, and protons in $pp$ collisions at $\sqrt{s}=0.9$, 2.76, and 7 TeV_, _Eur.Phys.J._ C72 (2012) 2164 [1207.4724].
  * [22] E. Albini, P. Capiluppi, G. Giacomelli and A.M. Rossi, _Mean charged hadron multiplicities in high-energy collisions_ , _Nuovo Cim._ A32 (1976) 101.
  * [23] D. Heck, J. Knapp, J.N. Capdevielle, G. Schatz and T. Thouw, _Corsika: a monte carlo code to simulate extensive air showers_ , _Wissenschaftliche Berichte, Forschungszentrum Karlsruhe_ FZKA 6019 (1998) .
  * [24] T.T. Böhlen, F. Cerutti, M.P.W. Chin, A. Fassò, A. Ferrari, P.G. Ortega et al., _The FLUKA Code: Developments and Challenges for High Energy and Medical Applications_ , _Nuclear Data Sheets_ 120 (2014) 211.
  * [25] Pierre Auger collaboration, _Measurement of the Fluctuations in the Number of Muons in Extensive Air Showers with the Pierre Auger Observatory_ , _Phys. Rev. Lett._ 126 (2021) 152002 [2102.07797].
  * [26] J. Matthews, _A Heitler model of extensive air showers_ , _Astropart. Phys._ 22 (2005) 387.
  * [27] L. Cazon, R. Conceição and F. Riehn, _Universality of the muon component of extensive air showers_ , _JCAP_ 03 (2023) 022 [2210.13407].
  * [28] S. Müller, R. Engel, T. Pierog and M. Roth, _Impact of muon detection thresholds on the separability of primary cosmic rays_ , _Astroparticle Physics_ 97 (2018) 174.
  * [29] L. Cazon, R. Conceição, M. Pimenta and E. Santos, _A model for the transport of muons in extensive air showers_ , _Astropart. Phys._ 36 (2012) 211 [1201.5294].
  * [30] I.C. Maris et al., _Influence of low energy hadronic interactions on air- shower simulations_ , _Nucl. Phys. Proc. Suppl._ 196 (2009) 86\.
  * [31] A. Fedynitch, R. Engel, T.K. Gaisser, F. Riehn and T. Stanev, _Calculation of conventional and prompt lepton fluxes at very high energy_ , _EPJ Web Conf._ 99 (2015) 08001 [1503.00544].