<figure><img src="image_1.png" /></figure>

Bright Gamma-Ray Flares Observed in GRB 131108A

M. Ajello<sup>1,5</sup>, M. Arimoto<sup>2,5</sup>, K. Asano<sup>3</sup>, M. Axelsson<sup>4,5</sup>, L. Baldini<sup>6</sup>, G. Barbiellini<sup>7,8</sup>, D. Bastieri<sup>9,10</sup>, R. Bellazzini<sup>11</sup>,

A. Berretta<sup>12</sup>, E. Bissaldi<sup>13,14</sup>, R. D. Blandford<sup>15</sup>, R. Bonino<sup>16,17</sup>, E. Bottacini<sup>15,18</sup>, J. Bregeon<sup>19</sup>, P. Bruel<sup>20</sup>, R. Buehler<sup>21</sup>,

S. Buson<sup>22</sup>, R. A. Cameron<sup>15</sup>, R. Caputo<sup>23</sup>, P. A. Caraveo<sup>24</sup>, E. Cavazzuti<sup>25</sup>, S. Chen<sup>9,18</sup>, G. Chiaro<sup>24</sup>, S. Ciprin<sup>26,27</sup>,

D. Costantin<sup>28</sup>, S. Cutini<sup>29</sup>, F. D'Ammando<sup>30</sup>, P. de la Torre Luque<sup>13</sup>, F. de Palma<sup>16</sup>, N. Di Lalla<sup>6</sup>, L. Di Venere<sup>13,14</sup>,

F. Fana Dirirsa<sup>31</sup>, S. J. Fegan<sup>20</sup>, A. Franckowiak<sup>21</sup>, Y. Fukazawa<sup>32</sup>, P. Fusco<sup>13,14</sup>, F. Gargano<sup>14</sup>, D. Gasparrini<sup>26,27</sup>,

N. Giglietto<sup>13,14</sup>, F. Giordano<sup>13,14</sup>, M. Giroletti<sup>30</sup>, D. Green<sup>33</sup>, I. A. Grenier<sup>34</sup>, M.-H. Grondin<sup>35</sup>, S. Guiriec<sup>23,36</sup>, E. Hays<sup>23</sup>,

D. Horan<sup>20</sup>, G. Jóhannesson<sup>37,38</sup>, D. Kocevski<sup>23</sup>, M. Kovac'ević<sup>29</sup>, M. Kuss<sup>11</sup>, S. Larsson<sup>5,39,40</sup>, L. Latronico<sup>16</sup>, J. Li<sup>21</sup>,

I. Liodakis<sup>15</sup>, F. Longo<sup>7,8</sup>, F. Loparco<sup>13,14</sup>, M. N. Lovellette<sup>41</sup>, P. Lubrano<sup>29</sup>, S. Maldera<sup>16</sup>, A. Manfreda<sup>6</sup>, G. Martí-Devesa<sup>42</sup>,

M. N. Mazziotta<sup>14</sup>, J. E. McEnery<sup>23,43</sup>, I. Mereu<sup>12,29</sup>, P. F. Michelson<sup>15</sup>, T. Mizuno<sup>44</sup>, M. E. Monzani<sup>15</sup>, E. Moretti<sup>45</sup>,

A. Morselli<sup>26</sup>, I. V. Moskalenko<sup>15</sup>, M. Negro<sup>16,17</sup>, E. Nuss<sup>19</sup>, M. Ohno<sup>46</sup>, N. Omodei<sup>15</sup>, M. Orienti<sup>30</sup>, E. Orlando<sup>15,47</sup>,

M. Palatiello<sup>7,8</sup>, V. S. Paliya<sup>21</sup>, D. Paneque<sup>33</sup>, Z. Pei<sup>10</sup>, M. Persic<sup>7,48</sup>, M. Pesce-Rollins<sup>11</sup>, V. Petrosian<sup>15</sup>, F. Piron<sup>19</sup>, H. Poon<sup>32</sup>,

T. A. Porter<sup>15</sup>, G. Principe<sup>30</sup>, J. L. Racusin<sup>23</sup>, S. Rainò<sup>13,14</sup>, R. Rando<sup>9,10</sup>, B. Rani<sup>23</sup>, M. Razzano<sup>11,59</sup>, A. Reimer<sup>15,42</sup>,

O. Reimer<sup>42</sup>, D. Serini<sup>13</sup>, C. Sgrò<sup>11</sup>, E. J. Siskind<sup>49</sup>, G. Spandre<sup>11</sup>, P. Spinelli<sup>13,14</sup>, H. Tajima<sup>15,50</sup>, K. Takagi<sup>32</sup>, D. Tak<sup>51,23</sup>,

D. F. Torres<sup>52,53</sup>, J. Valverde<sup>20</sup>, K. Wood<sup>54</sup>, R. Yamazaki<sup>55</sup>, M. Yassine<sup>7,8</sup>, S. Zhu<sup>56</sup>, Z. Lucas Uhm<sup>57</sup>, and Bing Zhang<sup>58</sup>

<sup>1</sup>Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634-0978, USA

<sup>2</sup>Faculty of Mathematics and Physics, Institute of Science and Engineering, Kanazawa University, Kakuma, Kanazawa, Ishikawa 920-1192, Japan

<sup>3</sup>Institute for Cosmic-Ray Research, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba, 277-8582, Japan

<sup>4</sup>Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden

<sup>5</sup>Department of Physics, KTH Royal Institute of Technology, AlbaNova, SE-106 91 Stockholm, Sweden

<sup>6</sup>Università di Pisa and Istituto Nazionale di Fisica Nucleare, Sezione di Pisa I-56127 Pisa, Italy

<sup>7</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, I-34127 Trieste, Italy

<sup>8</sup>Dipartimento di Fisica, Università di Trieste, I-34127 Trieste, Italy

<sup>9</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Padova, I-35131 Padova, Italy

<sup>10</sup>Dipartimento di Fisica e Astronomia "G. Galilei," Università di Padova, I-35131 Padova, Italy

<sup>11</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Pisa, I-56127 Pisa, Italy

<sup>12</sup>Dipartimento di Fisica, Università degli Studi di Perugia, I-06123 Perugia, Italy

<sup>13</sup>Dipartimento di Fisica "M. Merlin" dell'Università e del Politecnico di Bari, I-70126 Bari, Italy

<sup>14</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Bari, I-70126 Bari, Italy

<sup>15</sup>W.W. Hansen Experimental Physics Laboratory, Kavli Institute for Particle Astrophysics and Cosmology, Department of Physics and SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94305, USA

<sup>16</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Torino, I-10125 Torino, Italy

<sup>17</sup>Dipartimento di Fisica, Università degli Studi di Torino, I-10125 Torino, Italy

<sup>18</sup>Department of Physics and Astronomy, University of Padova, Vicolo Osservatorio 3, I-35122 Padova, Italy

<sup>19</sup>Laboratoire Univers et Particules de Montpellier, Université Montpellier, CNRS/IN2P3, F-34095 Montpellier, France

<sup>20</sup>Laboratoire Leprince-Ringuet, École polytechnique, CNRS/IN2P3, F-91128 Palaiseau, France

<sup>21</sup>Deutsches Elektronen Synchrotron DESY, D-15738 Zeuthen, Germany

<sup>22</sup>Institut für Theoretische Physik und Astrophysik, Universität Würzburg, D-97074 Würzburg, Germany

<sup>23</sup>NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>24</sup>INAF-Istituto di Astrofisica Spaziale e Fisica Cosmica Milano, via E. Bassini 15, I-20133 Milano, Italy

<sup>25</sup>Italian Space Agency, Via del Politecnico snc, I-00133 Roma, Italy

<sup>26</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Roma "Tor Vergata," I-00133 Roma, Italy

<sup>27</sup>Space Science Data Center—Agenzia Spaziale Italiana, Via del Politecnico, snc, I-00133, Roma, Italy

<sup>28</sup>University of Padua, Department of Statistical Science, Via 8 Febbraio, 2, I-35122, Padova, Italy

<sup>29</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Perugia, I-06123 Perugia, Italy

<sup>30</sup>INAF Istituto di Radioastronomia, I-40129 Bologna, Italy

<sup>31</sup>Department of Physics, University of Johannesburg, P.O. Box 524, Auckland Park 2006, South Africa

<sup>32</sup>Department of Physical Sciences, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan

<sup>33</sup>Max-Planck-Institut für Physik, D-80805 München, Germany

<sup>34</sup>AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, F-91191 Gif-sur-Yvette, France

<sup>35</sup>Centre d'Études Nucléaires de Bordeaux Gradignan, IN2P3/CNRS, Université Bordeaux 1, BP120, F-33175 Gradignan Cedex, France

<sup>36</sup>The George Washington University, Department of Physics, 725 21st St., NW, Washington, DC 20052, USA

<sup>37</sup>Science Institute, University of Iceland, IS-107 Reykjavik, Iceland

<sup>38</sup>Nordita, Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

<sup>39</sup>The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden

<sup>40</sup>School of Education, Health and Social Studies, Natural Science, Dalarna University, SE-791 88 Falun, Sweden

<sup>41</sup>Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA

<sup>42</sup>Institut für Astro- und Teilchenphysik, Leopold-Franzens-Universität Innsbruck, A-6020 Innsbruck, Austria

<sup>43</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA

<sup>44</sup>Hiroshima Astrophysical Science Center, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan

<sup>45</sup>Institut de Física d'Altes Energies (IFAE), Edifici Cn, Universitat Autònoma de Barcelona (UAB), E-08193 Bellaterra (Barcelona, Spain)

<sup>46</sup>MPA Research Group for Physical Geodesy and Geodynamics, H-1585 Budapest, Hungary

<sup>47</sup>Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, and Università di Trieste, I-34127 Trieste, Italy

<sup>48</sup>Osservatorio Astronomico di Trieste, Istituto Nazionale di Astrofisica, I-34143 Trieste, Italy

<sup>49</sup>NYCB Real-Time Computing Inc., Lattingtown, NY 11560-1025, USA

<sup>50</sup>Solar-Terrestrial Environment Laboratory, Nagoya University, Nagoya 464-8601, Japan

<sup>51</sup>Department of Physics, University of Maryland, College Park, MD 20742, USA; donggeun.tak@gmail.com

<sup>52</sup>Institute of Space Sciences (CSICIEEC), Campus UAB, Carrer de Magrans s/n, E-08193 Barcelona, Spain
Received 2019 September 25; revised 2019 November 11; accepted 2019 November 11; published 2019 November 27

Abstract

GRB 131108A is a bright long gamma-ray burst (GRB) detected by the Large Area Telescope and the Gamma-ray Burst Monitor on board the *Fermi Gamma-ray Space Telescope*. Dedicated temporal and spectral analyses reveal three $\gamma$-ray flares dominating above 100 MeV, which are not directly related to the prompt emission in the Gamma-ray Burst Monitor band (10 keV–10 MeV). The high-energy light curve of GRB 131108A (100 MeV–10 GeV) shows an unusual evolution: a steep decay, followed by three flares with an underlying emission, and then a long-lasting decay phase. The detailed analysis of the $\gamma$-ray flares finds that the three flares are 6–20 times brighter than the underlying emission and are similar to each other. The fluence of each flare, (1.6 $\sim$ 2.0) $\times 10^{-6}$ erg cm$^{-2}$, is comparable to that of emission during the steep decay phase, $1.7 \times 10^{-6}$ erg cm$^{-2}$. The total fluence from three $\gamma$-ray flares is $5.3 \times 10^{-6}$ erg cm$^{-2}$. The three $\gamma$-ray flares show properties similar to the usual X-ray flares that are sharp flux increases, occurring in $\sim$50% of afterglows, in some cases well after the prompt emission. Also, the temporal and spectral indices during the early steep decay phase and the decaying phase of each flare show the consistency with a relation of the curvature effect ($\hat{\alpha} = 2 + \hat{\beta}$), which is the first observational evidence of the high-latitude emission in the GeV energy band.

*Unified Astronomy Thesaurus concepts: Gamma-ray bursts* (629)

1. Introduction

Gamma-ray bursts (GRBs), the most luminous electro-
magnetic events in the universe, show two emission phases: the
prompt emission and the afterglow. The prompt emission, short
and spiky pulses, dominates in the keV–MeV energy range
with multiple spectral components (Guiriec et al. 2015, and
references therein). On the other hand, the light curve and
spectrum of the afterglow, emission from the interaction
between an outgoing blast wave from the central engine and a
circumburst medium (Mészáros & Rees 1997; Sari et al. 1998),
are characterized by a series of broken power laws, sometimes
accompanying bright flares (Nousek et al. 2006; Zhang et al.
2006). The afterglow is observed in a broad energy band from
radio to $\gamma$-ray. The flares are commonly observed in the X-ray
band (e.g., Romano et al. 2006), but rarely in the optical band
(e.g., Roming et al. 2006). The X-ray flares have been
explained as a result of the late-time activities of the central
engine (e.g., Fan & Wei 2005; Falcone et al. 2006, 2007;
Zhang et al. 2006; Chincarini et al. 2007; Galli & Piro 2007;
Lazzati & Perna 2007).

Due to the curvature effect of a spherical, relativistic jet
producing an X-ray flare, the decay phase of the X-ray flare
evolves in a certain way. This effect was first discussed by
Fenimore et al. (1996), and Kumar & Panaitescu (2000)
characterized the evolution of the temporal decaying index of
the X-ray flare ($\hat{\alpha}$) as a function of corresponding spectral index
($\hat{\beta}$), $\hat{\alpha} = 2 + \hat{\beta}$ in convention of $F_{\nu} \propto t^{-\hat{\alpha}\nu^{-\hat{\beta}}}$.
This relation has
been identified in many X-ray flares (e.g., Liang et al. 2006;
Chincarini et al. 2007; Jia et al. 2016; Uhm & Zhang 2016).

The *Fermi Gamma-ray Space Telescope* (*Fermi*) has
observed numerous GRBs and helped to uncover the exotic
high-energy evolution of GRBs. The high-energy emission
(>100 MeV) of GRBs observed by the Large Area Telescope
(LAT) on board *Fermi* shares common features: delayed onset

and lasting longer compared to keV–MeV emission, requiring
additional spectral components, and a power-law decaying
light curve (Ackermann et al. 2013; Ajello et al. 2019).
These GeV features can be interpreted as the early afterglow
emission (e.g., Ghisellini et al. 2010; Kumar & Barniol
Duran 2010; Tak et al. 2019). Abdo et al. (2011) reported
the GeV emission during vigorous X-ray flaring activities, but a
flare above the underlying afterglow emission in the GeV
energy band has not been reported before this work.

In this work, we first report the three bright $\gamma$-ray flares observed in GRB 131108A, which are $\sim$6–20 times brighter than the underlying light curve (Figure 1). We will compare spectral and temporal properties of three $\gamma$-ray flares and the X-ray flares. The broadband spectral analysis and the correlation test between the low- and high-energy bands will be described.

2. Observations

At 20:41:55.76 UTC on 2013 November 8 ($T_0$), LAT triggered on a bright high-energy emission from GRB 131108A (Racusin et al. 2013), which is simultaneously observed by the Gamma-ray Burst Monitor (GBM) on board *Fermi* (Younes 2013). The duration of the burst ($T_{90}^{60}$), is 18.2 s, but the high-energy emission lasts $\sim T_0 + 600$ s. With the observation of the bright afterglow of GRB 131108A by various instruments such as *Swift* (Chester & Stroh 2013; Stroh & Kennea 2013), *AGILE* (Giuliani et al. 2013), the accurate location and redshift of GRB 131108A were reported as (R.A., decl.) = (156.50, 9.66) with an uncertainty of 3"/6 in radius (Stroh & Kennea 2013) and $z \sim 2.40$ (de Ugarte Postigo et al. 2013), respectively. *Swift/XRT* started to observe the afterglow of GRB 131108A 3.9 ks ($\sim 1$ hr) after the GBM trigger, where the X-ray light curve decays smoothly in time (Stroh & Kennea 2013).

We perform a time-resolved analysis of LAT data in energy
range of 100 MeV–10 GeV with the *Fermi Science Tools*

<footnote>$^{53}$ Institució Catalana de Recerca i Estudis Avançats (ICREA), E-08010 Barcelona, Spain</footnote>

<footnote>$^{54}$ Praxis Inc., Alexandria, VA 22303, Resident at Naval Research Laboratory, Washington, DC 20375, USA</footnote>

<footnote>$^{55}$ Department of Physics and Mathematics, Aoyama Gakuin University, Sagamihara, Kanagawa, 252-5258, Japan</footnote>

<footnote>$^{56}$ Albert-Einstein-Institut, Max-Planck-Institut für Gravitationsphysik, D-30167 Hannover, Germany</footnote>

<footnote>$^{57}$ Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea</footnote>

<footnote>$^{58}$ Department of Physics and Astronomy, University of Nevada, Las Vegas, NV 89154, USA</footnote>

<footnote>$^{59}$ Funded by contract FIRB-2012-RBFR12PM1F from the Italian Ministry of Education, University and Research (MIUR).</footnote>

<footnote>$^{60}$ A duration where a GRB emits from 5% of its total counts to 95%.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** GBM and LAT light curves of GRB 131108A. The energy fluxes in the LAT energy band (100 MeV to 10 GeV) and in the GBM energy band (10 keV–1 MeV) are plotted in blue and orange, respectively. They are calculated from the best-fit model for each time interval in the spectral analysis with each instrument. The solid black line shows the fit of the LAT light curve consisting of five components: a simple power law (dotted gray line), a broken power law (dotted gray line), and three Gaussian functions (dotted red, green, and blue lines).</figcaption></figure>

(v11r5p3). We use "Transient020E" class events with the standard cuts. Photons from within the 15° region of interest around GRB 131108A are considered, and the maximum zenith angle is set to 100°. We fit a background rate from 3FGL sources (Acero et al. 2015), the galactic diffuse emission (gll_iem_v06), and the isotropic diffuse emission (iso_P8R2_TRANSIENT020_V6_v06).<fnref n="61" /> The LAT events observed in GRB 131108A are binned. For four sequential LAT events, we perform an unbinned likelihood analysis, and compute a test statistic (TS)<fnref n="62" /> for the burst. If the resultant TS is lower than 9 (equivalent to $3\sigma$), we add the next event to the bin and compute the TS again. Once we have the bin with the TS $\ge 9$, we collect the following four events, and repeat this procedure. As a result, each bin contains at least four LAT events, resulting in a TS $\ge 9$. For each of these bins, we perform an unbinned maximum likelihood fit on the energy spectrum with a simple power-law (PL) model.

The high-energy light curve of GRB 131108A shows an unusual evolution compared to other bright LAT GRBs (Figure 1); rather, it resembles the canonical X-ray early afterglow light curve though compressed to earlier and shorter timescales (seconds to tens of seconds compared to hundreds to thousands of seconds; Nousek et al. 2006; Zhang et al. 2006). We find the best description of the LAT light curve by fitting several models and their parameters with the maximum likelihood method. The light curve is well-fitted with five components ($\chi^2/\text{dof} = 20.8/34$): a simple power law, a broken power law, and three Gaussian functions for the three bright pulses above an underlying emission. Note that single-component models such as a simple power law ($\chi^2/\text{dof} = 130.8/47$) or a broken power law ($\chi^2/\text{dof} = 90.8/45$) are not a good model for this light curve. The model of each pulse can be replaced with the Norris function (Norris et al. 1996) ($\chi^2/\text{dof} = 20.5/31$), a broken power law ($\chi^2/\text{dof} = 20.2/31$), or a smoothly broken power law (Liang et al. 2006) ($\chi^2/\text{dof} = 19.4/28$). However, the Gaussian function is the best-fit model considering its statistics and the number of free parameters. The best-fit parameters for the three Gaussian

<table><caption>**Table 1**
The Physical Properties of Three Flares</caption><thead><tr><th></th><th>Peak Flux<sup>a</sup><br>(×10<sup>-6</sup> erg cm<sup>-2</sup> s<sup>-1</sup>)</th><th>Peak Time<br>(s)</th><th>FWHM<sup>b</sup><br>(s)</th><th>Fluence<sup>a</sup><br>(×10<sup>-6</sup> erg cm<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Flare 1</td><td>2.6 ± 1.3</td><td>3.5 ± 0.1</td><td>0.8 ± 0.2</td><td>2.0 ± 0.9</td></tr><tr><td>Flare 2</td><td>1.4 ± 0.6</td><td>6.0 ± 0.1</td><td>0.9 ± 0.2</td><td>1.6 ± 0.6</td></tr><tr><td>Flare 3</td><td>2.4 ± 1.0</td><td>8.5 ± 0.1</td><td>0.9 ± 0.2</td><td>1.6 ± 0.6</td></tr></tbody></table>

**Notes.**

functions are listed in Table 1. Note that there are hints of more than three flares but other fluctuations are insignificant, which are composed of one or two flux point. The fluence of each pulse is (2.0 ± 0.8), (1.6 ± 0.6), and (1.7 ± 0.7) × 10⁻⁶ erg cm⁻², totally (5.3 ± 1.2) × 10⁻⁶ erg cm⁻². The fluence of each pulse is comparable to that of emission during the early steep decay phase, (1.7 ± 0.4) × 10⁻⁶ erg cm⁻². The decaying index of the later segment of the broken power law is 1.6 ± 0.2, consistent with other Fermi-LAT GRBs (Ajello et al. 2019).

Considering the LAT light curve and its best-fit model, we define three time periods (Figure 2): the early steep decay period (time period 1; $T_0+0.3$ s–$T_0+2$ s), the three γ-ray unusual pulses with the underlying emission period (time period 2; $T_0+2$ s–$T_0+10$ s), and the long-lasting shallow decay period (time period 3; $T_0+10$ s–$T_0+20$ s). Note that the third time period can be extended until the end of the LAT emission, but stops at the end of the prompt emission for the joint-fit purpose. The evolution of the first and last periods is commonly seen in the LAT GRBs, but the phenomena of the second time period are noteworthy.

For the three time periods, we perform a broadband spectral analysis with GBM and LAT data in energy band from 10 keV to 10 GeV. Of the 12 NaI and 2 BGO detectors that make up GBM, 4 NaI detectors (0, 3, 6, and 7; 10 keV–1 MeV), and 2 BGO detectors (0 and 1; 200 keV–40 MeV) show a considerably high count rate above the background level, so that we

<footnote>a In 100 MeV–10 GeV.</footnote>

<footnote>b Full width at half maximum.</footnote>

<footnote n="61">⁶¹ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html</footnote>

<footnote n="62">⁶² The detection significance of the source above the background.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Count-rate curves and spectral energy distributions (SEDs). The top panel shows the scaled count-rate curves in different energy bands. The lower two panels show the joint-fit SEDs in energy band from 10 keV to 10 GeV. The color coding of the shaded region in the top panel and the spectrum in the middle panel indicates three time periods: pink (time period 1; 0.3–2 s), violet (time period 2; 2–10 s), and orange (time period3; 10–20 s). The bottom panel shows SEDs for three GeV flares (red, blue, and green). Each solid curve represents the best-fit spectral shape (thick) with 1 confidence level contour (shaded region) derived from the errors on the fit parameters.</figcaption></figure>

constitute a set of data from these detectors. In addition to the
GBM data, LAT Low Energy (LLE; 30 MeV–100 MeV) and
LAT (Transient020E; 100 MeV–10 GeV) data are used. The
background rate of GBM and LLE data is estimated by making
use of the analysis package, RMfit (version43pr2), by fitting a
time interval combined before and after the prompt emission

phase of the burst with a polynomial function. The LAT
background is estimated with “gtbkg”^63 provided by the Fermi

<footnote n="63">63 This tool generates a background spectrum file, which contains the total background rate from 3FGL sources, the galactic diffuse source, and the isotropic diffuse source.</footnote>
<table><caption>**Table 2**
Joint-fit Spectral Analysis Parameters</caption><thead><tr><th>Time Period</th><th>Model</th><th>α</th><th>β</th><th>Γ</th><th>E<sub>p,low</sub><br>(keV)</th><th>E<sub>p,high</sub><br>(MeV)</th><th>E<sub>f</sub><br>(MeV)</th><th>PG-stat</th><th>dof</th><th>BIC</th></tr></thead><tbody><tr><td>1<br>(0.3–2 s)</td><td>Band</td><td>-0.50<sup>+0.06</sup><sub>-0.06</sub></td><td>-2.07<sup>+0.02</sup><sub>-0.02</sub></td><td></td><td>320.6<sup>+22.3</sup><sub>-20.1</sub></td><td></td><td></td><td>745</td><td>695</td><td>771</td></tr><tr><td>2<br>(2–10 s)</td><td>Band</td><td>-0.48<sup>+0.07</sup><sub>-0.06</sub></td><td>-2.01<sup>+0.02</sup><sub>-0.02</sub></td><td></td><td>305.3<sup>+21.5</sup><sub>-19.8</sub></td><td></td><td>350.3<sup>+143.1</sup><sub>-87.6</sub></td><td>722</td><td>694</td><td>754</td></tr><tr><td></td><td>CPL + CPL</td><td>-0.73<sup>+0.02</sup><sub>-0.02</sub></td><td>-2.24<sup>+0.01</sup><sub>-0.01</sub></td><td></td><td>389.7<sup>+14.2</sup><sub>-13.4</sub></td><td></td><td></td><td>880</td><td>695</td><td>906</td></tr><tr><td></td><td>CPL + CPL</td><td>-0.50<sup>+0.07</sup><sub>-0.07</sub></td><td></td><td>1.69<sup>+0.02</sup><sub>-0.03</sub></td><td>394.2<sup>+13.3</sup><sub>-12.4</sub></td><td>126.2<sup>+25.6</sup><sub>-19.5</sub></td><td></td><td>821</td><td>693</td><td>860</td></tr><tr><td>3<br>(10–20 s)</td><td>Band</td><td>-0.94<sup>+0.05</sup><sub>-0.05</sub></td><td>-2.15<sup>+0.02</sup><sub>-0.02</sub></td><td></td><td>168.0<sup>+12.5</sup><sub>-11.2</sub></td><td></td><td></td><td>741</td><td>695</td><td>767</td></tr><tr><td></td><td>CPL + CPL</td><td>-0.92<sup>+0.13</sup><sub>-0.09</sub></td><td></td><td>1.68<sup>+0.07</sup><sub>-0.11</sub></td><td>198.9<sup>+13.7</sup><sub>-12.0</sub></td><td>68.7<sup>+21.5</sup><sub>-15.5</sub></td><td></td><td>742</td><td>693</td><td>781</td></tr><tr><td>Flare 1</td><td>CPL + CPL</td><td>-0.5<sup>fixed</sup><sub>b</sub></td><td></td><td>1.64<sup>+0.04</sup><sub>-0.04</sub></td><td>447.9<sup>+19.9</sup><sub>-18.7</sub></td><td>131.3<sup>+63.3</sup><sub>-35.7</sub></td><td></td><td>819</td><td>692<sup>c</sup></td><td>852</td></tr><tr><td>Flare 2</td><td>CPL + CPL</td><td>-0.5<sup>fixed</sup><sub>b</sub></td><td></td><td>1.67<sup>+0.04</sup><sub>-0.04</sub></td><td>396.0<sup>+18.1</sup><sub>-16.9</sub></td><td>193.1<sup>+126.7</sup><sub>-62.0</sub></td><td></td><td>763</td><td>692<sup>c</sup></td><td>796</td></tr><tr><td>Flare 3</td><td>CPL + CPL</td><td>-0.5<sup>fixed</sup><sub>b</sub></td><td></td><td>1.76<sup>+0.03</sup><sub>-0.03</sub></td><td>229.0<sup>+17.5</sup><sub>-15.7</sub></td><td>284.9<sup>+419.4</sup><sub>-124.4</sub></td><td></td><td>716</td><td>692<sup>c</sup></td><td>749</td></tr></tbody></table>

**Notes.**

a The cutoff energy $E_c$ is fixed to 50 MeV.

b The photon index for the low-energy CPL component is fixed to -0.5, which is the photon index of the best-fit model in the second time interval.

c Note that the change in the dof results from decrease in the number of the energy bin of LAT data.

**Science Tools.** The energy spectrum of each time period is fitted with various models by using the maximum likelihood method with Xspec (12.9.1) (Arnaud 1996). We use a Poisson data with Gaussian background statistic (PG-stat) for the parameter estimation, and then use the Bayesian information criterion (BIC; Schwarz 1978) for comparing the likelihood of fit and selecting the best-fit model. The best-fit model is a model with the lowest BIC value.

In the first time period characterized by the short bright emission commonly observed in the broad energy band from 10 keV to 10 GeV, the best-fit model for this time period is the Band function (Band et al. 1993) with a high-energy cutoff<sup>64</sup> (Table 2). The decrease in BIC as a result of adding the high-energy cutoff to the Band function is ~17 units, implying that the high-energy cutoff is strongly required. One alternative model is a combination of two spectral components, Band and a power law with an exponential cutoff (CPL).<sup>65</sup> This Band + CPL model describes the data slightly better (lower PG-stat), but the statistical improvement is not high enough to compensate the increase of a free parameter, making BIC higher than the best-fit model.

The second time interval where we found the unusual pulses shows a high count rate only in the low- (10 keV–1 MeV) and high-energy regimes (100 MeV–10 GeV) (upper panel in Figure 2). The observed data is best explained by a two-component model, CPL + CPL<sup>66</sup> (orange in Figure 2), which is preferred over a single-component model such as Band (Table 2). The two CPLs have distinct peak energies, $E_{p,low}$ ~400 keV and $E_{p,high}$ ~130 MeV, respectively (Table 2).

When any one of the CPL components is replaced with the Band function, $\beta$ becomes very soft so that the high-energy segment of Band is indistinguishable from the exponential cutoff. Therefore, the combination of Band and CPL is not necessary. The two-component scenario for GRB 131108A is also reported by Giuliani et al. (2014), who analyzed the AGILE (350 keV–30 GeV) data and reached the conclusion that the extrapolation of the low-energy spectral component could not explain the high-energy emission, and an additional spectral component with a peak energy at few MeV is required. The CPL dominating in the low-energy band has $\alpha \simeq -0.5$ consistent with that of the best-fit model in the first time period (Table 2), implying that the low-energy emission of the first and second time periods may be continuous. Given the best-fit model, the LAT emission is described by the high-energy CPL component. In addition, we perform a time-resolved spectral analysis for time intervals during the three high-energy pulses, and two distinct spectral components are again observed (Table 2 and Figure 2). The fluence of this high-energy component during the second time period is $5.9_{-0.8}^{+0.5} \times 10^{-6}$ erg cm$^{-2}$, comparable to the sum of fluence from three $\gamma$-ray pulses, $(5.2 \pm 1.2) \times 10^{-6}$ erg cm$^{-2}$ (Table 1). Most of the LAT emission during the second time period can be dominated by the three $\gamma$-ray pulses, and thus the high-energy CPL component may represent the spectral shape of the three $\gamma$-ray pulses.

During the third time period, short-soft pulses in the low-energy band (<500 keV) are observed. The best-fit model in this time period is the Band function (green in Figure 2). A CPL + CPL model does not give a better result, which requires two more parameters but resulting in the similar statistics (Table 2). After $T_0 + 20$ s, the LAT spectrum is well-described by a power law with a photon index $\Gamma = 2.8 \pm 0.3$.

Figure 1 and the upper panel of Figure 2 show that the low- (keV to few MeV) and high-energy (100 MeV–10 GeV) light curves evolve differently, and the broadband spectral analysis reveals the presence of the two spectral components. We check the correlation between the low- and high-energy light curves with the discrete correlation function (Edelson & Krolik 1988), which compares the variability of two light curves and estimates the time lag and the respective cross-correlation coefficient (e.g., Rani et al. 2009). For this purpose, we performed a time-resolved spectral analysis for the time bin

64

$$ \frac{dN}{dE} = \begin{cases} N_0 \left( \frac{E}{100 \text{ keV}} \right)^\alpha \exp \left( \frac{E(\alpha + 2)}{E_p} \right) & \text{if } E \le \frac{\alpha - \beta}{2 + \alpha} E_p, \\ N_0 \left( \frac{E}{100 \text{ keV}} \right)^\beta \left( \frac{E_p}{100 \text{ keV} 2 + \alpha} \right)^{\alpha - \beta} \exp(\beta - \alpha) & \text{if } \frac{\alpha - \beta}{2 + \alpha} E_p < E \le E_c, \\ N_0 \left( \frac{E}{100 \text{ keV}} \right)^\beta \left( \frac{E_p}{100 \text{ keV} 2 + \alpha} \right)^{\alpha - \beta} \exp(\beta - \alpha) \exp \left( \frac{E_c - E}{E_f} \right) & \text{if } E > E_c, \end{cases} \quad (1) $$

where $\alpha$ and $\beta$ are the low- and high-energy photon indices, respectively, $E_p$ is the peak energy of the Band function, $E_c$ is the cutoff energy which is fixed to 50 MeV, and $E_f$ is the e-folding energy for the high-energy cutoff.

$$ 65 \frac{dN}{dE} = N_0 \left( \frac{E}{100 \text{ keV}} \right)^\alpha \exp \left( - \frac{E(\alpha + 2)}{E_p} \right). $$

$$ 66 \frac{dN}{dE} = N_{0,\text{low}} \left( \frac{E}{100 \text{ keV}} \right)^\alpha \exp \left( - \frac{E(\alpha + 2)}{E_{p,\text{low}}} \right) + N_{0,\text{high}} \left( \frac{E}{100 \text{ keV}} \right)^\Gamma \exp \left( - \frac{E(-\Gamma + 2)}{E_{p,\text{high}}} \right). $$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Test of the curvature effect for early steep decay emission and three γ-ray flares. The upper panels show the light curve of the early steep decay emission and the three γ-ray flares separately after removing the $t_0$ effect. The lower panel shows the evolution of the photon index. The data points corresponding to the decay phase of each flare are in red, blue, and green. The early steep decay phase is marked in black. The solid lines in the upper panels show the expected light curve derived from the relation of the curvature effect, $\hat{\alpha} = 2 + \hat{\beta}$. These theoretical light curves are computed from the linear fit of the photon indices of the decay phase (solid line in the lower panel).</figcaption></figure>

used for the LAT light curve (Figure 1) with the GBM data, and computed flux in the range of 10 keV–10 MeV with the best-fit model for each time interval (orange points in Figure 1). When the entire interval is considered, the correlation between the two light curves, 10 keV–10 MeV and 100 MeV–10 GeV, is evident (correlation coefficient peak = 0.8 ± 0.1). However, if only the light curves in the second time period is considered, the correlation analysis does not suggest any correlation between the two light curves (correlation coefficient peak ~0.1).

Considering the temporal and spectral features, the γ-ray pulses invoke a distinct origin from the prompt emission of the low-energy band as well as the LAT extended emission. From now on, the individual γ-ray pulse is called a “γ-ray flare.”

### 3. Discussion

First of all, we compare and test the well-known properties of X-ray flares to the observation of the γ-ray flares in GRB 131108A.

A flux variation of the X-ray flares, a flux ratio of a flare to an underlying emission ($\delta F_\nu/F_\nu$), ranges from 6 (e.g., GRB 050406; Romano et al. 2006) to higher than 100 (e.g., GRB 050202B; Burrows et al. 2005). The γ-ray flares are about 6–20 times brighter than the underlying emission, which slowly changes in time, $F_\nu \sim 2.7 \times 10^{-7}$ erg cm$^{-2}$ s$^{-1}$ on average (see Figure 1), and thus the flux variation of the γ-ray flares is subnormal, compared to X-ray flares (Chincarini et al. 2007). The duration of the X-ray flares varies from few hours to a day (Chincarini et al. 2007; Swenson & Roming 2014), and there is an empirical relation between the onset time and the duration of the X-ray flares, $\delta t/t \sim 0.1$ (e.g., Chincarini et al. 2007, 2010; Swenson & Roming 2014). In case of the γ-ray flares, they last only a few seconds (Table 1), much shorter than the X-ray flares (Chincarini et al. 2007). Also, the γ-ray flares are observed in $<T_0 + 10$ s, which is earlier than any X-ray flares (Chincarini et al. 2007). Combining these two unusual features, surprisingly, the temporal characteristics of the γ-ray flares are not in conflict with the empirical relation. A comparison between the flux variability and the temporal variability of the γ-ray flares ($\delta F_\nu/F_\nu$ versus $\delta t/t$) shows that the γ-ray flares are consistent with X-ray flares (Chincarini et al. 2007). Furthermore, this comparison implies that the γ-ray flares are not

related to the fluctuations of the external shock as previously discussed for the X-ray flares (Ioka et al. 2005; Zhang et al. 2006).

The steep decay of the X-ray flares is regarded as a result of the curvature effect, which is identified by testing the relation, $\hat{\alpha} = 2 + \hat{\beta}$ (e.g., Liang et al. 2006; Chincarini et al. 2007; Jia et al. 2016; Uhm & Zhang 2016). It is possible that the decay phase of the γ-ray flares also show evidence of the curvature effect. Before testing the relation, we should remove the so-called $t_0$ effect (Zhang et al. 2006; Kobayashi & Zhang 2007). Each flare is attributed to the late-time activity of the central engine and thus has its own onset time ($t_0$). Since the shape of a light curve in the logarithmic space is very sensitive to the choice of $t_0$, the intrinsic light curve of the flare can only be provided if the light curve is shifted to the true $t_0$. Due to the underlying emission, however, the true onset of the γ-ray flares is ambiguous. Therefore, we properly choose the onset of each flare as the time when the flux of the flare is 1/100 of its peak. Figure 3 shows the γ-ray flares after shifting them to the proper $t_0$ for each flare. For these light curves, we test the curvature effect relation. After selecting the data points corresponding to the decaying phase, we fit the measured photon indices ($\Gamma$) with a linear function, $\Gamma = f(t-t_0)$ (solid line in lower panel of Figure 3). Next, the photon index is converted to the spectral index, $\hat{\beta} = \Gamma - 1$. We then apply the HLE relation and get the temporal index as a function of time, $\hat{\alpha} = f(t-t_0) + 1$. Finally, the light curve expected by the curvature effect is described by a function of time, $F_\nu = F_{\nu,0} (t-t_0)^{f(t-t_0)+1}$ (solid line in upper panel of Figure 3). We fit this function with the observed, shifted flux points and conclude that the decay phases of all three γ-ray flares are consistent with the expectation by the curvature effect (Figure 3). Also, we find the spectral softening during the decay phase of the flares, which is the well-known phenomenon identified in the X-ray flares (Chincarini et al. 2007; Falcone et al. 2007).

The X-ray flares are likely attributed to internal shocks, where accelerated electrons at the shocks radiate via the synchrotron process. On the other hand, the γ-ray flares with $E_p \sim 130$ MeV may originate from the synchrotron self-Compton (SSC) process from the same population of electrons that might have produced X-ray flares. In principle, there could be two possible cases for the inverse Compton process: SSC from the
internal emission region and the external inverse Compton (EIC) from the external shock region (Wang et al. 2006; Fan et al. 2008). The observation of the high-latitude emission in the flares disfavors the EIC origin and supports the SSC origin. Assuming the typical electron Lorentz factor, $\gamma_e \gtrsim 10^3$, the peak energy of the seed photon should be $\lesssim 0.13$ keV, which is far below the GBM energy band. Another possibility for the origin of the $\gamma$-ray flares is the residual collision in the internal dissipation process (Li & Waxman 2008). In this case, the $\gamma$-ray flares can be interpreted as the SSC counterpart of the optical emission produced by the residual collision at large radii. Note that there were no X-ray and optical observations during the prompt emission phase of this burst, so that these hypotheses cannot be tested.

The very first steep decay emission in the first time period corresponds to the tail of the first bright broadband pulse. This decay emission is also consistent with the curvature effect (the first panel in Figure 3).

The underlying emission in the second time period can be interpreted as the emission during the development of the forward shock (e.g., Maxham et al. 2011), and the long-lasting decay emission (the third time period) can be the continuous emission from the fully developed forward shock when the total energy is not noticeably increased by the additional energy injection (e.g., Mészáros & Rees 1997; Sari et al. 1998).

The observation of GRB 131108A uncovers a new phenomenon in the high-energy GRB light curve. Even though the three $\gamma$-ray flares were observed in the prompt phase of the burst, they showed the temporal and spectral properties similar to those of the X-ray flares. Also, we found the evidence of the curvature effect in the GeV energy band for the first time.

The Fermi-LAT Collaboration acknowledges support for LAT development, operation, and data analysis from NASA and DOE (United States), CEA/Irfu and IN2P3/CNRS (France), ASI and INFN (Italy), MEXT, KEK, and JAXA (Japan), and the K.A. Wallenberg Foundation, the Swedish Research Council and the National Space Board (Sweden). Science analysis support in the operations phase from INAF (Italy) and CNES (France) is also gratefully acknowledged. This work performed in part under DOE Contract DE-AC02-76SF00515.

Software: XSPEC (v12.9.1; Arnaud 1996), Fermi Science Tools (v11r5p3).

## ORCID iDs

M. Ajello https://orcid.org/0000-0002-6584-1703

M. Arimoto https://orcid.org/0000-0003-1250-7872

M. Axelsson https://orcid.org/0000-0003-4378-8785

F. D'Ammando https://orcid.org/0000-0001-7618-7527

A. Franckowiak https://orcid.org/0000-0002-5605-2219

M. Kuss https://orcid.org/0000-0003-1212-9998

S. Larsson https://orcid.org/0000-0003-0716-107X

F. Longo https://orcid.org/0000-0003-2501-2270

S. Maldera https://orcid.org/0000-0002-0698-4421

T. Mizuno https://orcid.org/0000-0001-7263-0296

I. V. Moskalenko https://orcid.org/0000-0001-6141-458X

V. S. Paliya https://orcid.org/0000-0001-7774-5308

J. L. Racusin https://orcid.org/0000-0002-4744-9898

R. Yamazaki https://orcid.org/0000-0002-1251-7889

## References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011, ApJL, 734, L27

Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23

Ackermann, M., Ajello, M., Asano, K., et al. 2013, ApJS, 209, 11

Ajello, M., Arimoto, M., Axelsson, M., et al. 2019, ApJ, 878, 52

Arnaud, K. A. 1996, in ASP Conf. Ser. 101, Astronomical Data Analysis Software and Systems V, ed. V. Systems, G. H. Jacoby, & J. Barnes (San Francisco, CA: ASP), 17

Band, D., Matteson, J., Ford, L., et al. 1993, ApJ, 413, 281

Burrows, D. N., Romano, P., Falcone, A., et al. 2005, Sci, 309, 1833

Chester, M. M., & Stroh, M. C. 2013, GCN, 1, 15476

Chincarini, G., Mao, J., Margutti, R., et al. 2010, MNRAS, 406, 2113

Chincarini, G., Moretti, A., Romano, P., et al. 2007, ApJ, 671, 1903

de Ugarte Postigo, A., Thoene, C. C., Gorosabel, J., et al. 2013, GCN, 15470, 1

Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646

Falcone, A. D., Burrows, D. N., Lazzati, D., et al. 2006, ApJ, 641, 1010

Falcone, A. D., Morris, D., Racusin, J., et al. 2007, ApJ, 671, 1921

Fan, Y.-Z., Piran, T., Narayan, R., & Wei, D.-M. 2008, MNRAS, 384, 1483

Fan, Y. Z., & Wei, D. M. 2005, MNRAS, 364, L42

Fenimore, E. E., Madras, C. D., & Nayakshin, S. 1996, ApJ, 473, 998

Galli, A., & Piro, L. 2007, A&A, 475, 421

Ghisellini, G., Ghirlanda, G., Nava, L., & Celotti, A. 2010, MNRAS, 403, 926

Giuliani, A., Longo, F., Verrecchia, F., et al. 2013, GCN, 1, 15479

Giuliani, A., Mereghetti, S., Marisaldi, M., et al. 2014, arXiv:1407.0238

Guiriec, S., Kouveliotou, C., Daigne, F., et al. 2015, ApJ, 807, 148

Ioka, K., Kobayashi, S., & Zhang, B. 2005, ApJ, 631, 429

Jia, L.-W., Uhm, Z. L., & Zhang, B. 2016, ApJS, 225, 17

Kobayashi, S., & Zhang, B. 2007, ApJ, 655, 973

Kumar, P., & Barniol Duran, R. 2010, MNRAS, 409, 226

Kumar, P., & Panaitescu, A. 2000, ApJL, 541, L51

Lazzati, D., & Perna, R. 2007, MNRAS, 375, L46

Li, Z., & Waxman, E. 2008, ApJL, 674, L65

Liang, E. W., Zhang, B., O'Brien, P. T., et al. 2006, ApJ, 646, 351

Maxham, A., Zhang, B.-B., & Zhang, B. 2011, MNRAS, 415, 77

Mészáros, P., & Rees, M. J. 1997, ApJ, 476, 232

Norris, J. P., Nemiroff, R. J., Bonnell, J. T., et al. 1996, ApJ, 459, 393

Nousek, J. A., Kouveliotou, C., Grupe, D., et al. 2006, ApJ, 642, 389

Racusin, J. L., Zhu, S., Kocevski, D., et al. 2013, GCN, 1, 15464

Rani, B., Wiita, P. J., & Gupta, A. C. 2009, ApJ, 696, 2170

Romano, P., Moretti, A., Banat, P. L., et al. 2006, A&A, 450, 59

Roming, P. W. A., Vanden Berk, D., Pal'shin, V., et al. 2006, ApJ, 651, 985

Sari, R., Piran, T., & Narayan, R. 1998, ApJL, 497, L17

Schwarz, G. 1978, AnSta, 6, 461

Stroh, M. C., & Kennea, J. A. 2013, GCN, 1, 15474

Swenson, C. A., & Roming, P. W. A. 2014, ApJ, 788, 30

Tak, D., Omodei, N., Uhm, Z. L., et al. 2019, ApJ, 883, 134

Uhm, Z. L., & Zhang, B. 2016, ApJL, 824, L16

Wang, X.-Y., Li, Z., & Mészáros, P. 2006, ApJL, 641, L89

Younes, G. 2013, GCN, 1, 15477

Zhang, B., Fan, Y. Z., Dyks, J., et al. 2006, ApJ, 642, 354