The Bright and the Slow -- GRBs 100724B \& 160509A with high-energy cutoffs at $\lesssim100\;$MeV

TL;DR

This study investigates two bright GRBs with unexpectedly low high-energy cutoffs, fitting models to determine their physical parameters and revealing potential biases in gamma-ray burst detection.

Contribution

It provides the first detailed modeling of high-energy cutoffs in these GRBs using opacity-based models, estimating Lorentz factors and exploring detection biases.

Findings

High-energy cutoffs at 20-150 MeV are well fitted by opacity models.

Estimated Lorentz factors range from 100 to 400.

Lower cutoff energies may bias Fermi-LAT detection against certain GRBs.

Abstract

We analyze the prompt emission of GRB 100724B and GRB 160509A, two of the brightest Gamma-Ray Bursts (GRBs) observed by Fermi at $\lesssim{\rm MeV}$ energies but surprisingly faint at $\gtrsim100\;$MeV energies. Time-resolved spectroscopy reveals a sharp high-energy cutoff at energies $E_c\sim20-60\;$MeV for GRB~100724B and $E_c\sim80-150\;$MeV for GRB~160509A. We first characterize phenomenologically the cutoff and its time evolution. We then fit the data to two models where the high-energy cutoff arises from intrinsic opacity to pair production within the source ($\tau_{\gamma\gamma}$): (i) a Band spectrum with $\tau_{\gamma\gamma}$ from the internal-shocks motivated model of Granot et al. (2008), and (ii) the photospheric model of Gill & Thompson (2014). Alternative explanations for the cutoff, such as an intrinsic cutoff in the emitting electron energy distribution, appear to be less natural. Both models provide a good fit to the data with very reasonable physical parameters, providing a direct estimate of bulk Lorentz factors in the range $\Gamma\sim 100-400$, on the lower end of what is generally observed in Fermi GRBs. Surprisingly, their lower cutoff energies $E_c$ compared to other Fermi-LAT GRBs arise not predominantly from the lower Lorentz factors, but also at a comparable level from differences in variability time, luminosity, and high-energy photon index. Finally, particularly low $E_c$ values may prevent detection by Fermi-LAT, thus introducing a bias in the Fermi-LAT GRB sample against GRBs with low Lorentz factors or variability times.

Tables

Table 1. Table 1 : Summary of various spectral models used in this work.

Model	Description	Ref.
$f_{Band}$	Band function	eq. (E2)
$f_{BB}$	Band function plus blackbody	eq. (E4)
$f_{BHec}$	Band function with high-energy exponential cutoff	eq. (E3)
$f_{Bbkpo}$	Band function with broken power-law spectrum above the peak	eq. (E5)
$f_{Bgr}$	Band function with a gradual break in power-law spectrum above the peak	eq. (E6)
$f_{BG}$	Band function with high-energy spectral break due to $\gamma \gamma$ pair opacity	§5.1; (Granot et al., 2008)
$f_{GT}$	Spectrum from delayed pair breakdown model in a strongly magnetized jet	§5.2; (Gill & Thompson, 2014)
$f_{th}$	Quasi-thermal spectrum described by a power-law plus a Wein peak	eq. (8); (Gill & Thompson, 2014)

Table 2. Table 2 : Value of the -log(likelihood) S 𝑆 S for the Band model f B ​ a ​ n ​ d subscript 𝑓 𝐵 𝑎 𝑛 𝑑 f_{Band} , and the T ​ S 𝑇 𝑆 TS obtained respectively with a Band + Black body ( f B ​ B subscript 𝑓 𝐵 𝐵 f_{BB} ) and a Band with exponential cutoff ( f B ​ H ​ e ​ c subscript 𝑓 𝐵 𝐻 𝑒 𝑐 f_{BHec} ) as alternative hypotheses, for GRB 100724B. In parenthesis we report the significance of the improvement. In the p null subscript 𝑝 null p_{\rm null} columns we also report the null-hypothesis probability for the models for LLE data. Given the limited number of simulations used to measure p null , we cannot reliably measure probabilities smaller than 10 − 3 superscript 10 3 10^{-3} . Therefore, we report p < null 10 − 3 {}_{\rm null}<10^{-3} in these cases.

#	Time interval	$S$ of $f_{Band}$	p	pLLE	TS of $f_{BB}$		p	pLLE	TS of $f_{BHec}$		p	pLLE
1	2.11 - 9.39	1374.3			16.8	($3.5\sigma$)	0.3	0.002	27.0	($5.2\sigma$)	0.38	0.23
2	9.39 - 13.11	1079.6			74.6	($8.3\sigma$)	0.05		74.2	($8.6\sigma$)	0.26	0.52
3	13.11 - 14.61	604.9	0.05	0.002	12.8	($2.9\sigma$)	0.13	0.004	16.4	($4\sigma$)	0.13	0.05
4	14.61 - 19.14	1227.2			97.6	($9.5\sigma$)	0.26		127.2	($11.2\sigma$)	0.93	0.61
5	19.14 - 27.23	1448.3			12.0	($2.8\sigma$)	0.12		57.3	($7.6\sigma$)	0.88	0.17
6	27.23 - 62.52	2174.9	0.1	0.001	50.8	($6.7\sigma$)	0.25		27.1	($5.2\sigma$)	0.21	0.09
7	62.52 - 63.65	544.7	0.09	0.03	3.8	($3.1\sigma$)	0.4	0.001	7.8	($2.8\sigma$)	0.59	0.82
8	63.65 - 82.23	1888.2			51.8	($6.8\sigma$)	0.26		40.6	($6.4\sigma$)	0.21	0.36
9	82.23 - 148.41	2393.0	0.07	0.12	7.2	($1.9\sigma$)	0.58	0.15	3.2	($1.8\sigma$)	0.78	0.32

Table 3. Table 3 : Value of the -log(likelihood) S 𝑆 S for the Band model f B ​ a ​ n ​ d subscript 𝑓 𝐵 𝑎 𝑛 𝑑 f_{Band} , and the T ​ S 𝑇 𝑆 TS obtained with a Band with exponential cutoff ( f B ​ H ​ e ​ c subscript 𝑓 𝐵 𝐻 𝑒 𝑐 f_{BHec} ) as alternative hypothesis, for GRB 160509B. In parenthesis we report the significance of the improvement. We also report the null-hypothesis probability p null for each model, for LAT and LLE data.

#	Time interval	$S$ of $f_{Band}$	p	pLAT	TS of $f_{BHec}$		p	pLAT
1	9.712-11.045	970.5	0.23	0.12	7.1	($2.4\sigma$)	0.14	0.15
2	11.045-12.042	894.07	0.36	0.05	7.6	($2.5\sigma$)	0.59	0.31
3	12.042-14.449	1516.65			93.8	($9.6\sigma$)	0.23	0.33
4	14.449-17.783	1620.88			42.1	($6.4\sigma$)	0.62	0.48
5	17.783-18.480	847.16	0.08	0.003	21.2	($4.5\sigma$)	0.44	0.23
6	18.480-18.667	233.18			67.5	($8.1\sigma$)	0.59	0.77
7	18.667-19.044	527.27			43.9	($6.5\sigma$)	0.09	0.24
8	19.044-20.249	1105.23			63.1	($7.8\sigma$)	0.87	0.89
9	20.249-21.787	1124.36	0.56	0.09	5.1	($2.0\sigma$)	0.07	0.35
10	21.787-25.254	1460.61	0.15	0.23	6.3	($2.2\sigma$)	0.19	0.22

Table 4. Table 4 : Ratio of intrinsic parameters as described in eq. ( 12 ) with GRB 160509A as the reference system. In systems where no spectral cutoff was observed (marked with ∗ ), the maximum observed photon energy is quoted. ∗∗ Minimum bulk- Γ Γ \Gamma (or the actual inferred Γ Γ \Gamma when a cutoff is observed) calculated using a one-zone analytical model employing a more elaborate radiation field spectrum as compared to a simple power-law used in LS01. In all three cases, the inferred Γ min ∼ 2 ​ Γ B ​ G similar-to subscript Γ min 2 subscript Γ 𝐵 𝐺 \Gamma_{\rm min}\sim 2\Gamma_{BG} :  a Abdo et al. ( 2009b ) , b Ackermann et al. ( 2011 ) , c Ackermann et al. ( 2010 )

	080916C	090926A	090510	160509A
$z$	4.35	2.1062	0.903	1.17
$E_c^{\ast }$ (GeV)	${3.0}^{\ast }$	0.4	${30.5}^{\ast }$	0.08
${\Gamma }_{\min }^{\ast \ast }$	$887\pm {21}^a$	$720\pm {76}^b$	$1218\pm {61}^c$	-
${\Gamma }_{BG}$	451	319	628	363
$\beta$	$-2.21$	$-1.71$	$-1.85$	$-1.60$
$L_{0,52}\left(\mathrm{erg}/s\right)$	55.78	3.42	1.24	6.18
$t_{v,z}$ (ms)	374	48	6.3	23
${\psi }_{\left(1+z\right)}$	$-0.20$	$-0.18$	$0.02$
${\psi }_{\Gamma }$	$0.41$	$-0.58$	$0.81$
${\psi }_{\beta }$	$0.58$	$0.41$	$0.26$
${\psi }_{L_{0,52}}$	$-0.81$	$0.50$	$0.46$
${\psi }_{t_{v,z}}$	$1.03$	$0.63$	$-0.37$

Equations

$$f_{BG}(E)=f_{Band}(\alpha ,\beta ,E_p,K)\Lambda (\beta ,\frac{\Delta R}{R_0},\psi ,m,b),$$

$${\Gamma }_2={\left[3.969\frac{10^{4(2+\beta )}{L}_{0,52}}{(1+z{)}^{\beta }{f}^{\beta +1}}{\left(\frac{-\beta }{2}\right)}^{-\frac{5}{3}}\frac{33.4\mathrm{ms}}{t_{\mathrm{v}}}\right]}^{1/(2-2\beta )}.$$

$$R_{\mathrm{br}}\sim 2{\Gamma }_{\mathrm{br}}^{2}c\mathcal{R}{t}_{\mathrm{eng}}=5.4\times 10^{12}{\Gamma }_{\mathrm{br},3}^2{\mathcal{R}}_1{t}_{\mathrm{eng},0}\mathrm{cm}.$$

$${\ell }_{\mathrm{th},\mathrm{br}}^{\prime }$$

$$R_{\mathrm{diss}}$$

$${\Gamma }_{\mathrm{diss}}$$

$${\ell }_{\mathrm{tot}}^{\prime }={\ell }_{\mathrm{th}}^{\prime }+{\ell }_{\mathrm{heat}}^{\prime }=(1+{\xi }_{\mathrm{th}}){\ell }_{\mathrm{heat}}^{\prime }\lesssim {\ell }_{\gamma ,\mathrm{crit}}^{\prime }$$

$$f_{th}^{\prime}(E^{\prime})\propto\left\{\begin{array}[]{ll}{E^{\prime}}_{\rm pk}^{2-\alpha}\exp({-E^{\prime}_{\rm pk}/k_{B}T^{\prime}_{\rm th}}){E^{\prime}}^{\alpha}&E^{\prime}<E^{\prime}_{\rm pk}\\ {E^{\prime}}^{2}\exp(-E^{\prime}/k_{B}T^{\prime}_{\rm th})&E^{\prime}>E^{\prime}_{\rm pk}\end{array}\right.$$

$$R=\frac{3{\sigma }_T}{16\pi m_e{c}^3}\frac{L_{\mathrm{iso}}}{{\Gamma }^3{\ell }_{\mathrm{tot}}^{\prime }}.$$

$$\Gamma$$

$$E_{c,z}=$$

$$=$$

$${\psi }_{\xi ,i}=\frac{\log (E_{c,0}({\xi }_i)/E_{c,0})}{\log (E_{c,z,i}/E_{c,z,0})},$$

$$E_{\mathrm{syn},\max }=\frac{\Gamma }{(1+z)}\frac{9}{4{\alpha }_F}{m}_e{c}^2=0.159\Gamma (1+z{)}^{-1}\mathrm{GeV}.$$

$${\Gamma }_{\gamma \gamma ,\max }(E_c)=(1+z)\frac{E_c}{m_e{c}^2}=176\left(\frac{1+z}{3}\right)\left(\frac{E_c}{30\mathrm{MeV}}\right),$$

$$E_c\ge E_{sa}=\frac{\Gamma }{1+z}{m}_e{c}^2=51.1\left(\frac{3}{1+z}\right)\left(\frac{\Gamma }{300}\right)\mathrm{MeV},$$

$${\tau }_{T,\pm }\sim ({\epsilon }_1^{\prime }{)}^{2(1+\beta )}={\left[\frac{{\Gamma }_{\gamma \gamma ,\min }(E_c)}{{\Gamma }_{\gamma \gamma ,\max }(E_c)}\right]}^{2(-1-\beta )}\ll 1,$$

$${\tilde{\tau }}_{T,\pm }\sim ({\epsilon }_1^{\prime }{)}^{1+\beta }={\left[\frac{{\Gamma }_{\gamma \gamma ,\min }(E_c)}{{\Gamma }_{\gamma \gamma ,\max }(E_c)}\right]}^{2(1-\beta )}\gg 1.$$

$$\Gamma \approx \min \left[{\Gamma }_{\gamma \gamma ,\min }(E_c),{\Gamma }_{\gamma \gamma ,\max }(E_c)\right].$$

$$E_c\approx \max \left[E_{sa}(\Gamma ),E_1(\Gamma )\right].$$

$$\eta ({\vec{p}}_s,t)+R<113^{\circ }-PSF(\theta ,E_{\min }),$$

$$f_{PHec}(E)=K{\frac{E}{E_{piv}}}^{-\alpha }{e}^{-\frac{E}{E_c}},$$

$$f_{Band}(E)=K\{\begin{array}{ll}{E}^{\alpha }\exp \left(\frac{-E}{E_0}\right)& E<(\alpha -\beta )E_0\\ {\left[(\alpha -\beta )E_0\right]}^{\alpha -\beta }\exp (\beta -\alpha )E^{\beta }& E\ge (\alpha -\beta )E_0\end{array},$$

$$f_{BHec}(E)=f_{Band}(E)\times \exp \left(-\frac{E}{E_c}\right),$$

$$f_{BB}(E)=f_{Band}(E)+A\frac{E^2}{\exp \left(\frac{E}{kT}\right)-1}$$

$$f_{Bbkpo}(E)=\{\begin{array}{ll}K{E}^{\alpha }\exp \left(\frac{-E}{E_0}\right)& E<(\alpha -\beta )E_0\\ K{\left[(\alpha -\beta )E_0\right]}^{\alpha -\beta }\exp (\beta -\alpha )E^{\beta }& (\alpha -\beta )E_0\le E<E_c\\ K{\left[(\alpha -\beta )E_0\right]}^{\alpha -\beta }\exp (\beta -\alpha )E_c^{\beta -{\beta }_2}{E}^{{\beta }_2}& E\ge E_c\end{array}$$

$$f_{Bgr}(E)=f_{Band}(E)\times {\left[1+{\frac{E}{E_c}}^{n\Delta \beta }\right]}^{-1/n},$$

$$\Delta \beta =\frac{(\beta +1)(2-\beta )}{\beta -1},$$

Full text

keywords

1 Introduction

The $\gamma$-ray emission from Gamma-Ray Bursts (GRBs) is believed to originate within an ultra-relativistic jet, which is launched during the collapse of a massive star (for long duration GRBs that last $\gtrsim 2\,$s, MacFadyen & Woosley, 1999) and likely also during the merger of two compact objects (for short duration GRBs that last $\lesssim 2\,$s, Rezzolla et al., 2011). However, the mechanisms that produce the prompt emission of GRBs are still debated (see e.g. the recent review by Kumar & Zhang, 2015). An important question is the composition of the jet, which remains unresolved, and for which two scenarios have been proposed: a baryonic jet where particles are accelerated converting thermal energy into bulk motion (fireballs) (Rees & Meszaros, 1994), or a Poynting-flux-dominated jet (Lyutikov & Blackman, 2001). The composition of the jet in turn determines the dominant dissipation mechanism that converts the energy content of the jet into heat and accelerated particles that radiate the observed prompt emission. For example, in baryonic jets energy dissipation can be attributed to internal shocks (e.g. Rees & Meszaros, 1994; Morsony et al., 2010; Lopez-Camara et al., 2013), and/or collisional heating due to inelastic collisions between neutrons and protons (Beloborodov, 2010). On the other hand, in a Poynting-flux-dominated jet, where most of the energy is stored in the magnetic field, magnetic reconnection occurring in an outflow with a striped magnetic field structure or due to magnetohydrodynamic turbulence can dissipate magnetic energy and power the prompt emission (e.g., Thompson, 1994; Lyutikov & Blandford, 2003; Zhang & Yan, 2011).

In the context of fireball models, the dominant emission mechanism was thought to be synchrotron radiation, possibly also accompanied by synchrotron self-Compton. In particular, the highly-variable prompt emission has been attributed to synchrotron emission from particles accelerated in multiple internal shocks, i.e., shocks that occur when a faster shell ejected by the central engine collides with a slower shell within the outflow. Such a scenario has been used to explain the non-thermal spectrum that characterizes GRBs. The efficiency that internal shocks can achieve in converting energy into radiation appears to be insufficient to explain the luminosity of some GRBs (Lazzati et al., 1999; Kobayashi et al., 1997), unless the spread in Lorentz factor between the colliding shells is large (Kobayashi & Sari, 2001). Also, a non-negligible fraction of GRBs show spectra that are difficult to explain with pure synchrotron emission (Preece et al., 2002; Yu et al., 2015a; Burgess et al., 2015; Axelsson & Borgonovo, 2015). For this reason, some GRBs have been modeled with phenomenological models adding a thermal component to the non-thermal one (Ryde, 2005; Guiriec et al., 2011; Axelsson et al., 2012; Guiriec et al., 2013; Burgess et al., 2014; Yu et al., 2015b; Guiriec et al., 2015; Nappo et al., 2017).

Because of these issues with the so-called “standard” fireball paradigm, another class of fireball models has emerged, which we call for simplicity photospheric models (for example Ryde, 2004; Pe’er et al., 2005; Beloborodov, 2010; Vurm et al., 2011; Lazzati et al., 2013). In this class of models the spectrum of a GRB is explained as reprocessed quasi-thermal radiation coming from the photosphere, i.e. the surface where radiation and matter decouple, typically after the acceleration of the fireball has ended for thermal acceleration, or possibly during the acceleration phase for magnetic acceleration (which is slower than thermal acceleration). A thermal or quasi-thermal initial spectrum is reprocessed within the jet to produce the non-thermal spectrum commonly observed in GRBs. The differences between the various photospheric models lie in the mechanisms responsible for the reprocessing of the thermal spectrum, which in turn requires different ingredients: strongly-magnetized or non-magnetized jets, baryon-dominated or baryon-poor, or other factors.

Here we present the analysis of the prompt emission of GRB 100724B and GRB 160509A, both detected by the Fermi Gamma-ray Space Telescope instruments. These two GRBs are very bright at low energy, but they do not show any emission above 1 GeV during the prompt phase, which sets them apart from bursts of comparable low-energy fluence such as GRB 080916C, GRB 090902B and GRB 090926 (Ackermann et al., 2013b). Moreover, the high-energy emission above 1 GeV, widely thought to originate from a different mechanism than the prompt emission (for example, external shock), picks up after the prompt phase is finished. This gives us the rare possibility of studying the prompt emission without any contamination from the high-energy component. Both GRBs show a very evident spectral cutoff in the $10-200$ MeV energy range with respect to the extrapolation of the low-energy component. We interpret it as pair production opacity, which allows for a measurement of the bulk Lorentz factor of the jet. While other cases of sub-GeV cutoffs have been reported (Ackermann et al., 2013b; Tang et al., 2015), GRB 100724B and GRB 160509A are by far the two brightest ones, and allow for an in-depth analysis impossible in the other cases. We also perform a detailed time-resolved analysis and measure the time evolution of the bulk Lorentz factor in both GRBs. Our detailed analysis allows us also to verify the viability of specific physical models. We choose to consider one model related to the “standard” fireball picture and one photospheric model. In particular, among many possibilities, we choose the semi-phenomenological internal-shock model of Granot et al. (2008) featuring a detailed modeling of the pair production opacity, and the photospheric model of Gill & Thompson (2014). These models provide a natural explanation for the spectral cutoff, and we have readily available numerical codes which provide the spectra foreseen by the two scenarios as a function of physical parameters (see section 5 for more details).

In § 2 we present the Fermi observatory. We then present the main features of GRB 100724B (§ 3) and GRB 160509A (§ 4). In particular, we establish phenomenologically that the high-energy data cannot be modeled extrapolating the low-energy spectrum, requiring instead a high-energy cutoff in the $10-200\;$MeV energy range. Next, in § 5 we interpret such a feature in the context of physical models. We finally discuss our results (§ 6) and provide our conclusions (§ 7). Throughout this paper we will use the “Planck 2015” flat cosmology (Planck Collaboration et al., 2016), with $H_{0}=67.8$ km s*-1* Mpc*-1* and $\Omega_{m}=0.308$.

2 The Fermi observatory

Fermi orbits the Earth at an altitude of $\sim 565$ km. Its pointing is continuously changing in a pattern that allows its instruments to survey the entire sky approximately every 3 hours.

The Large Area Telescope (LAT) (Atwood et al., 2009) is a pair-conversion telescope operating in the energy range from around 20 MeV up to over 300 GeV. For this study we use the P8_TRANSIENT020E class of LAT data, and the corresponding instrument response function, and the LAT Low-Energy data (LLE), available on the Fermi Science Support Center (FSSC) website111http://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermille.html. When compared to P8_TRANSIENT020E data, LLE data feature a higher acceptance especially below 100 MeV, at the expense of a higher background contamination and a very limited spatial resolution. It is designed for the spectral analysis of short-duration transients such as GRBs and solar flares.

On board Fermi is also the Gamma-Ray Burst Monitor (GBM). It is comprised of 12 sodium iodide (NaI) detectors sensitive in the 8 keV $-$ $\sim$1 MeV energy range, and 2 bismuth germanate (BGO) detectors sensitive in the 200 keV $-$ $\sim$40 MeV energy range. The detectors are arranged to allow GBM to probe continuously all the sky not occulted by the Earth, with the exception of the time interval when the spacecraft is going through the South Atlantic Anomaly and data taking is suspended. In this work we use the GBM data and tools publicly available on the FSSC website.

3 GRB 100724B

3.1 Observations

The bright GRB 100724B triggered Fermi/GBM (Meegan et al., 2009) at 00:42:05.99 on 2010-07-24 (Bhat, 2010) ($t_{0}$ in the following). It was also detected by Fermi/LAT and a preliminary localization was reported (Tanaka et al., 2010). GRB 100724B was also detected by Konus-Wind (Golenetskii et al., 2010), AGILE (Marisaldi et al., 2010; Giuliani et al., 2010; Del Monte et al., 2011) and Suzaku (Uehara et al., 2010). This burst has the third greatest fluence to date at low energy ($<10\;$MeV) among all the LAT-detected GRBs, exceeded only by GRB 090902B and the record-breaking GRB 130427A (Ackermann et al., 2013b, 2014). The initial localization has been improved in Ackermann et al. (2013b). We use in this paper an even more refined localization, $R.A.=123.47\,^{\circ}$ and $Dec.=75.88\,^{\circ}$ (J2000), obtained as described in Appendix A.

Another burst, GRB 100724A, was detected few seconds later by Swift (Markwardt et al., 2010) in a position occulted by the Earth for Fermi. Therefore, even if Fermi/GBM is a non-imaging full-sky monitor, this second GRB was not observed by any GBM detector and does not therefore affect the analysis presented in this work. However, follow up efforts focused on this second GRB and therefore no multi-wavelength data are available for GRB 100724B.

The light curve of GRB 100724B is shown in Figure 1. During the main emission episode the signal in the LAT was exceptionally intense in the 30 MeV–100 MeV energy range, but nothing was detected above 100 MeV. There is a faint “precursor” peak before $t_{0}$, the main emission episode going from $t_{0}$ to $\sim t_{0}+150\,$s, and then a late soft peak starting at $\sim t_{0}+180$ s.

3.2 Spectral analysis of the prompt emission

We consider GBM detectors NaI 0 and 1, because they are the only two low-energy detectors seeing the GRB at an off-axis angle of less than $40\,^{\circ}$. Furthermore, we select the BGO detector closest to the GRB direction (BGO 0). We use Time-Tagged-Events data provided by the GBM team and publicly available on the FSSC website. We generate custom response matrices (rsp2 files, with one new response every time the spacecraft slew by 0.5 deg) using the public tool gbmrspgen222http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gbmrspgen.html and using our best localization for the source. We use NaI data in the energy range 8 keV $-$ $\sim$1 MeV but excluding the energy range $30-40$ keV, which contains the K-edge feature. We use BGO data from channel 2 to channel 125 (corresponding to the energy range $\sim 217$ keV - $38$ MeV). We also use LAT LLE data above 30 MeV.

We estimate the background for all GBM detectors and for LLE data by fitting off-pulse intervals with one polynomial function for each channel, and then interpolating such fit to the on-pulse interval (for details see Ackermann et al., 2013b). This way the time-varying background –and in particular the Earth Limb contribution – is naturally taken into account.

We first perform a time-integrated analysis using the same time interval used in Ackermann et al. (2013b), i.e., the GBM $t_{90}$ time interval, from $t_{0}+8.195\;$s to $t_{0}+122.882\;$s. We use the Multi-Mission Maximum Likelihood framework (3ML) for all spectral analysis performed in this paper (Vianello et al., 2015). We find very similar results: the spectrum can be successfully modeled by a Band function, a phenomenological model traditionally used in describing GRB spectra (Band et al., 1993), multiplied by an exponential cutoff. The formulas for the Band function and the Band with exponential cutoff function are reported in Appendix E, eqs. (E2) and (E3). The best fit parameters for the time-integrated analysis are $\alpha=-0.69\pm 0.02$, $\beta=-2.01^{-0.02}_{+0.03}$, $E_{p}=330\pm 10\;$keV and $E_{c}=48\pm 6\;$MeV. The fluence in the 1 keV –$\,10\;$GeV energy range is $(4.7\pm 0.3)\times 10^{-4}\;{\rm erg\;cm^{-2}}$. Given the very high signal-to-noise ratio of this time-integrated analysis we can rebin the spectra in order to have at least $30$ counts in each bin, and then we can use a standard $\chi^{2}$ test. We obtain $\chi^{2}=376.2$ for $391$ d.o.f., corresponding to a p-value of $\sim 0.7$.

In order to further study this feature, we then perform a time-resolved spectral analysis. The choice of the time intervals requires a trade-off. Choosing many time bins gives good time resolution but low sensitivity for detecting features, due to the decreased statistics in each spectrum. On the contrary, choosing few bins gives good sensitivity at the risk of smearing the time evolution of the parameters. In this paper we are mainly interested in the study of the cutoff, thus we choose to focus on LLE data, which cover the energy range where the cutoff is measured, and we decide our time bins based on the variability seen in the LLE light curve. In particular, we apply the Bayesian Blocks algorithm (Scargle et al., 2013) that finds the most probable segmentation of the observation into time intervals during which the photon arrival rate has no statistically significant variations, i.e., it is perceptibly constant. We have used the implementation provided in the tool gtburst of the Fermi Science Tools, using a probability of false positives of $p_{0}=0.01$. Applying it to LLE data we find the 9 intervals between 0 and $\sim 150\;$s shown in Figure 1.

We note that these intervals do not cover the faint and soft “precursor” peak that can be seen between $t_{0}-42\;$s and $t_{0}-13.3\;$s, nor the faint and soft late peak between $t_{0}+175\;$s and $t_{0}+235\;$s, since they do not show any LLE emission. For the precursor, we find that it is well described by the power-law with exponential cutoff model $f_{PHec}$ (eq. E1), with parameters $\alpha=-0.4^{-0.3}_{+0.4}$ and $E_{p}=130^{-50}_{+90}$ keV**, and a p-value computed as above of $p=0.62$. The faint and soft late peak is described again by $f_{PHec}$, with parameters $\alpha=-1.35^{-0.12}_{+0.13}$ and $E_{p}=142^{-35}_{+60}$ keV ($p=0.45$).**

We now focus on the main emission episode. We extract the spectra and compute the response matrices for each detector and each interval. Initially we consider a pool of commonly used phenomenological spectral models (summarized in Table 1) in order to characterize the spectra without having to assume a specific theoretical framework. We will consider two specific physical models later on (see section 5). Our phenomenological models are based on the Band model: a) the Band model itself $f_{Band}$ (eq. E2); b) a Band plus black-body model $f_{BB}$ (eq. E4), which was used for the modeling of this GRB in Guiriec et al. (2011); c) the Band model multiplied by an exponential cutoff $f_{BHec}$ (eq. E3); d) a Band model where the high-energy power law changes photon index abruptly at a cutoff energy ($f_{Bbkpo}$, eq. (E5)); and e) a Band model with a smooth spectral break, suggested on theoretical grounds in Granot et al. (2008) ($f_{Bgr}$, eq. (E6)). We also apply a group of alternative models, namely the log-parabolic spectral shape (Massaro et al., 2010), a broken power law, and the smoothly broken power law of Ryde (1998). However, they yield large residuals and in all time intervals considered here they describe the data significantly worse than the models based on the Band function. Therefore, we disregard them from now on. We also use a procedure to mitigate effects due to inter-calibration issues between the instruments. We take one instrument as reference (NaI 0), and then we introduce a multiplicative constant for every other detector. Such constant is left free to vary in the fit between 0.7 and 1.3, corresponding to an inter-calibration uncertainty of up to 30%. This “effective area correction” reduces the biases due to systematic errors in the total effective area of the instruments with respect to the reference one.

For each time interval we measure separately the significance of the black body in model $f_{BB}$ and of the exponential cutoff in model $f_{BHec}$ with respect to the Band model alone $f_{Band}$. We rely on the Likelihood Ratio Test, which uses as Test Statistic ($TS$) twice the difference in log-likelihood between the null hypothesis (the Band model in our case) and the alternative hypothesis (either $f_{BHec}$ or $f_{BB}$ in our case). The details of this procedure can be found in Appendix B. We also measure the p-value for a goodness-of-fit test using a procedure equivalent to the classic $\chi^{2}$ test but more appropriate for Poisson data. In particular, we follow the method proposed by Cousins (2013) based on Monte Carlo simulations. It is well known that the goodness-of-fit p-value $p$ can be misleading when the data points have very different uncertainties, because the points with smaller errors will dominate. This is the case in our situation, where GBM data provide a much larger statistic than LLE data. Therefore, we measure separately the null-hypothesis probability for the entire dataset ($p$) and for the LLE data alone ($p_{\rm LLE}$) in order to investigate whether a model is able to describe the data both at low and at high energy.

We report the results in Table 2. The TS of $f_{BHec}$ with respect to $f_{Band}$ and the corresponding significance of the improvement (9th column) is large ($>4\sigma$) for all intervals except for the two where the GRB is faint. The quality of the fit is good both overall and for LLE data in particular, as shown by the p-values $p$ and $p_{\rm LLE}$ (last two columns). On the other hand, the improvement obtained with $f_{BB}$ with respect to $f_{Band}$ is large ($>4\sigma$) for 4 intervals (6th column). The values for the overall null-hypothesis probability $p$ for $f_{BB}$ seem to indicate a good fit (7th column), however while the model describes well the low-energy data it does not describe well LLE data, as shown by $p_{\rm LLE}$ (8th column, see also fig. 2). Hence, we can conclude that while $f_{BB}$ models well the low-energy data – as already concluded in Guiriec et al. (2011) – it fails to describe well LLE data.

Summarizing, $f_{BHec}$ is a more parsimonious model than $f_{BB}$ (it has one less free parameter) and also provides a better description of LLE data in all intervals. It is therefore our model of choice. This result appears to be at odds with what is reported in Guiriec et al. (2011). We note however that these authors did not use LLE data, which is where the advantage of $f_{BHec}$ over $f_{BB}$ becomes evident, and used different time intervals for the time-resolved spectral analysis. They also used a different localization for the GRB, provided by the Fermi/GBM with a large localization error, which has an impact on the response matrices used for GBM data and therefore on the modeling of the spectrum (Connaughton et al., 2015; Burgess et al., 2018). This GRB has also been studied by Del Monte et al. (2011) using AGILE data. The spectrum they measured is much harder than what Fermi measured, and with a much larger flux. If the characteristics measured by AGILE were true, we would have detected with the LAT a large number of photons above 100 MeV which we do not see. We discuss in Appendix D a plausible motivation for this discrepancy.

The procedure described here considers only statistical uncertainties. A study of the effects of systematic uncertainties on the significance of the cutoff that are not neutralized by the use of the “effective area correction” is reported in Appendix C, and demonstrates that the improvement given by the cutoff is unlikely to be due to systematic uncertainties in the response of the instruments.

The significance of the cutoff with respect to the simple extrapolation of the low-energy spectrum being established, we compare $f_{BHec}$ with the two models with power-law shape after the cutoff ($f_{Bbkpo}$ and $f_{Bgr}$) to assess whether the spectrum is curved (exponential cutoff) or not (power law break) after the break. We find that $f_{Bbkpo}$ and $f_{Bgr}$ never provide a better fit as measured by $p_{null}$ with respect to the exponential shape despite having more parameters, which favors a curved spectrum above the cutoff.

In Figure 2 we show the best fit parameters for $f_{BHec}$ for the intervals of the main emission episodes. The parameters $\alpha$ and $\beta$ decrease during the first peak, increase in the second peak, and then decrease again. $E_{\rm p}$ shows a similar evolution. This tracking behavior is common in GRBs (Ford et al., 1995; Kaneko et al., 2006). The cutoff energy $E_{\rm c}$ increases slightly with time.

4 GRB 160509A

4.1 Observations

GRB 160509A triggered Fermi/GBM on 2016-05-09 at $t_{0}=$ 08:58:45.22 UTC. It was also localized on-board by Fermi/LAT (Longo et al., 2016), one of only 5 cases over 8 years of mission. This allowed for a quick follow up and localization by Swift (Kennea et al., 2016), which in turn allowed for a redshift measurement by Gemini/North of $z=1.17$ (Tanvir et al., 2016) when the afterglow was still bright. We adopt the position of the afterglow measured by Gemini North (R.A. $=311.7538\,^{\circ}$, Dec. $=76.1081\,^{\circ}$). The prompt emission (Figure 3) consists of a soft “precursor” peak between $t_{0}-5.0$ and $\sim t_{0}+5.5$ s, followed by a much brighter main episode which lasts until $t_{0}+40$ s. After a quiescent time, there is another very soft emission episode, visible only in the low-energy detectors, from $\sim t_{0}+300$ s until $\sim t_{0}+400$ s. This excess localizes in roughly the same direction as the main episode, although with large statistical uncertainty (GBM team, private communication), therefore it is likely to be associated with GRB 160509A. Similarly to GRB 100724B, during the main emission episode the LAT detected many photons associated with the GRB in the 30 –100 MeV energy range, but surprisingly few above 100 MeV (see last two panels in Figure 3).

4.2 Spectral analysis of the prompt emission

We use here the same technique and energy selections discussed in section 3.2. For the first and second episode we use data from GBM detectors NaI 0, NaI 3 and BGO 0, which are the detectors in the most favorable position to observe the GRB. Since the pointing of the Fermi satellite changed between the first two emission episodes and the third one, for the latter we used NaI 0, NaI 6, NaI 9 and BGO 1, which were the detectors closest to the direction of the GRB at that time. Contrary to the case of GRB 100724B, the Earth Limb was far from the direction of GRB 160509A during the prompt emission. For all intervals we hence use LAT LLE data from 30 MeV up to 100 MeV, and LAT standard data above 100 MeV.

The spectrum accumulated over the entire duration of the GRB, from $t_{0}$ to $t_{0}+400$ s, can be well described with the $f_{BHec}$ model. Thanks to the high signal-to-noise ratio we have many counts in each bin in the spectrum; thus we can assume Gaussian statistics and apply a normal $\chi^{2}$ test. We obtain $\chi^{2}=392.4$ with 381 d.o.f, corresponding to a p-value of $p=0.33$.. The cutoff is required with high significance ($>10\sigma$) with respect to the Band model alone. We measure a fluence of $(3.2\pm 0.5)\times 10^{-4}$ erg cm*-2* in the 1 keV – 10 GeV energy range, corresponding to an isotropic emitted energy of $E_{\rm iso}=(2.6\pm 0.4)\times 10^{54}$ erg. The contribution to this quantity by the precursor and the late emission episode is negligible.

The spectrum of the “precursor” peak is well described by a power law with exponential cutoff (eq. E1), with $\alpha=-1.03^{-0.06}_{+0.07}$, $E_{c}=410^{-80}_{+100}\;$keV, and $K=1.42^{-0.3}_{+0.4}\;$ph. cm*-2* s*-1*. This is very similar to the “precursor” peak in GRB 100724B.

The third, late episode is faint and soft as well. We divide it in two intervals, $297.45-358.9\;$s and $358.9-400.88\;$s from $t_{0}$. Their spectra are both well described by a Band model. The best fit parameters are respectively $\alpha=-1.21^{-0.07}_{+0.06}$, $\beta=-1.98^{-0.24}_{+0.10}$, $E_{p}=270^{-50}_{+90}\;$keV, and $\alpha=-1.27^{-0.04}_{+0.05}$, $\beta=-2.20^{-0.13}_{+0.09}$, $E_{p}=138^{-13}_{+15}\;$keV. Adding an exponential cutoff, or any other component like a thermal component, does not significantly improve the fit. This can of course either be intrinsic, or just due to the lack of sufficient statistics, especially at high energies.

We focus then on the main episode, much brighter than the other two. The Band model $f_{Band}$ overestimates the amount of LLE signal by a large amount and the improvement obtained by adding an exponential cutoff to the Band model is very large. We obtain $TS=278$ for $f_{BHec}$, corresponding to a significance of $16.6~{}\sigma$. The addition of a black body, instead, returns a lower $TS=120$. Moreover, the $f_{BB}$ model does not describe well LAT data, yielding very large residuals. This is reflected by the p-values returned by the $\chi^{2}$ test – which again we can apply in virtue of the very high statistics – that are respectively $0.25$ for the $f_{BHec}$ model and $\sim 10^{-6}$ for the $f_{BB}$ model.. We therefore do not consider $f_{BB}$ as a viable model for the time integrated analysis for this GRB and our current energy and interval selection.

As for GRB 100724B, we run the Bayesian Blocks algorithm on LLE data with the same setup to determine the time intervals for the time-resolved spectral analysis of the main episode. We show these intervals as the black lines in Figure 3.

In the 6 intervals where the GRB is bright we find again that the addition of an exponential cutoff to the Band spectrum improves the fit significantly ($>4\sigma$), as shown by the 5th column in table 3. There we also report the p-value for the goodness of fit test for the entire dataset ($p$) and for LAT and LLE data ($p_{\rm LAT}$), computed as described in section 3.2. It shows that the $f_{BHec}$ model provides a good description both overall and for LLE and LAT data in particular. In the case of this GRB, contrary to GRB 100724B, the addition of a black body to the Band spectrum (i.e., the $f_{BB}$ model) does not yield a significant improvement in most intervals. Moreover, residuals in LLE and LAT data are very significant. Therefore, we avoid computing the p-values for the goodness-of-fit test for this model, which is very computationally intensive, and do not consider it as a good model for this GRB given our selections.

As for GRB 100724B we tested whether the two models with a power-law shape after the cutoff ($f_{Bbkpo}$ and $f_{Bgr}$) provide a better fit than $f_{BHec}$. The $p_{null}$ obtained with these models is worse for all intervals despite their increased complexity with respect to $f_{BHec}$. Hence, as for GRB 100724B, the Band with exponential cutoff model $f_{BHec}$ provides a better fit with less parameters and is therefore statistically preferred, and the shape of the spectrum after the cutoff appears to be curved.

The best fit parameters as a function of time are shown in the left panel in Figure 4. The time evolution of $\alpha$ and $E_{\rm p}$ is similar to the case of GRB 100724B: there is a decreasing trend with some variability tracking the light curve in correspondence of the second peak, while the cutoff energy appears to increase slightly with time. Also $\beta$ is tracking the flux, becoming harder when the flux increases, but there is no decreasing trend. The luminosity as a function of time is shown in the right panel of Figure 4, computed in the energy range 1 keV – 10 GeV: the values of a few $10^{53}$ erg s*-1* are quite typical for long-duration GRBs (Yonetoku et al., 2004).

5 Interpretation and physical modeling

In the previous sections we have established phenomenologically that an high-energy exponential cutoff in GRB 100724B and GRB 160509A must be added to the extrapolation of the low-energy Band spectrum in order to successfully model LAT data. In this section we provide some possible interpretations.

Among many possibilities (for example Ryde, 2004; Pe’er et al., 2005; Beloborodov, 2010; Vurm et al., 2011; Lazzati et al., 2013), we consider two scenarios: i) the cutoff is due to pair-production opacity that attenuates a non-thermal spectrum (produced for example by synchrotron emission during internal shocks); or ii) a photospheric model where the cutoff arises due to the development of an electron-positron pair cascade in a highly magnetized, dissipative and baryon-poor outflow. In the first scenario we adopt a hybrid approach by considering the phenomenological Band model, traditionally used in modeling the non-thermal spectrum of GRBs, and we multiply it by a $\gamma$-$\gamma$ attenuation factor computed from first principles in Granot et al. (2008). It features a self-consistent semi-analytic calculation of the impulsive emission from a thin spherical ultra-relativistic shell (model $f_{BG}$, see eq. 1). The calculation accounts for the fact that, in impulsive relativistic sources, the timescale for significant variations in the properties of the radiation field within the source is comparable to the total duration of the emission episode, and therefore, the dependence of the opacity to pair production on space and time cannot be ignored. In the second scenario we instead adopt the photospheric model of Gill & Thompson (2014), which we will call $f_{GT}$ in the following, as this model produces spectra which are strikingly similar to the phenomenological $f_{BHec}$ model that is a good description of the data. We describe both models in some detail next.

5.1 Pair Opacity Break in Impulsive Relativistic Outflows - the $f_{BG}$ model

The model of Granot et al. (2008) features an expanding ultra-relativistic spherical thin shell. The emission, which may arise from internal shock heated electrons, is assumed to be isotropic in the shell’s comoving frame. Its comoving luminosity scales as a power law with dimensionless energy $\varepsilon^{\prime}=E^{\prime}_{\rm ph}/m_{e}c^{2}$, where $m_{e}$ is the electron mass and $c$ is the speed of light, and with radius $R$, $L^{\prime}_{\varepsilon^{\prime}}\propto(\varepsilon^{\prime})^{1+\beta}R^{b}$, where $\beta$ is the (high-energy) photon index. The shell’s Lorentz factor (LF) is also assumed to vary as a power law with radius, $\Gamma\propto R^{-m/2}$. The emission episode lasts between radii $R_{0}$ and $R_{0}+\Delta R$, where the fractional radial width $\Delta R/R_{0}$ determines how impulsive the emission is.

The optical depth $\tau_{\gamma\gamma}$ to pair production ($\gamma\gamma\to e^{+}e^{-}$) is calculated along the trajectory of each test photon that reaches the observer. Its contribution to the observed flux is attenuated by a factor of $\exp(-\tau_{\gamma\gamma})$, leading to a quasi-exponential (after adding contributions from different emission radii and angles) cutoff in the instantaneous spectrum. Depending on the value of $\Delta R/R_{0}$, the time integrated spectrum either features a smoothly broken power law cutoff ($\Delta R/R_{0}\sim 1$) or a quasi-exponential cutoff that asymptotes into a power law ($\Delta R/R_{0}\gtrsim 1$).

In order to enable a (semi-) analytic calculation, the effects of the pairs that are produced in this process are neglected. This is a reasonable approximation as long as their Thomson opacity is $\tau_{{\rm T},\pm}\ll 1$. Below we will examine the validity of this approximation.

In practice, we compute the attenuation factor333We note that here that $-\beta$ plays the role of the parameter called $\alpha$ in Granot et al. (2008). $\Lambda(\beta,\Delta R/R_{0},\psi,m,b)$ due to the $\gamma$-$\gamma$ opacity through a numerical code which implements the computation described in Granot et al. (2008). We then define $f_{BG}$ as:

[TABLE]

where $f_{Band}$ is the Band model (eq. E2), and $\psi=E/\varepsilon_{z}$ where $\varepsilon_{z}=(1+z)\varepsilon$ is the dimensionless photon energy in the source’s cosmological frame at redshift $z$, while $\varepsilon$ is the value measured at Earth. In order to reduce the number of free parameters to a manageable number, we fix $m=0$, which corresponds to a shell in coasting phase as expected from an internal shock scenario, and $b=0$, which corresponds to assuming a comoving spectral emissivity independent of radius. Therefore, we have 6 free parameters: $\alpha,\beta,E_{p},K,\Delta R/R_{0},\psi$.

The LF $\Gamma_{0}=100\Gamma_{2}$ can be estimated by using this relation: 444The numerical coefficient in the expression here for $\Gamma_{2}=\Gamma_{0}/100$ is larger by a factor of $\pi^{1/(1-\beta)}$ compared to Eq. (126) of Granot et al. (2008), correcting an error in the latter equation.

[TABLE]

Here $L_{0,52}=4\pi d_{L}^{2}(1+z)^{-\beta-2}F_{0}/(10^{52}\;{\rm erg~{}s^{-1}})$, where $d_{L}$ is the luminosity distance of the burst, $F_{0}$ is the (unabsorbed) energy flux ($\nu F_{\nu}$) foreseen by the high-energy power law of the Band model at 511 keV. The parameter $f=\psi/m_{e}c^{2}$ relates to the parameter $C_{2}=10^{-(6+4\beta)}\left(\frac{f}{(E_{c}/5.11~{}{\rm GeV})}\right)^{\beta+1}$ that appears in eq. (126) of Granot et al. (2008). Its exact value is not known a priori and can only be determined numerically. To that end, we obtain the value of $C_{2}$ from model fits to both GRBs considered in this work and, as expected in Granot et al. (2008), find that its value is of order unity (see Figure 5). We extract from the data $t_{\rm mv}$, which is the minimum variability time scale detected in the light curve, defined as the rise time of the shortest significant structures. Therefore, $t_{\rm mv}\approx(1+z)\Delta R/2c\Gamma_{0}^{2}=t_{0}\Delta R/R_{0}$ where $t_{0}=(1+z)R_{0}/2c\Gamma_{0}^{2}$ is the arrival time of the first photons to the observer. Since we chose to define the variability timescale as $t_{v}\equiv 2t_{0}$ when deriving eq. (2), we obtain that it can be expressed in terms of $t_{\rm mv}$ as follows: $t_{v}=2~{}t_{\rm mv}/(\Delta R/R_{0})$.

We determine $t_{\rm mv}$ using a wavelet analysis. Similar techniques have been already used by many authors to study the variability of GRBs (Walker et al., 2000; MacLachlan et al., 2013; Golkhou et al., 2015). In contrast to these authors, we adopt the Continuous Wavelet Transform (CWT) in place of the Discrete Wavelet Transform (DWT), as the CWT allows for a much better resolution in the spectrum (Torrence & Compo, 1998).

We start by obtaining a light curve with a bin size of $10^{-4}$ s of the entire time interval with a bright LLE emission (respectively $t_{0}+8.195$ – $t_{0}+122.882\;$s for GRB 100724B and $t_{0}+10$ – $t_{0}+25\;$s for GRB 160509A). Then we compute the wavelet power spectrum $W$ as a function of the time scale $\delta t$, as described in Torrence & Compo (1998), with the correction suggested by Liu et al. (2007). The result is shown in Figure 6 (dots). In order to measure the variance of the power spectrum due to the Poisson fluctuations of the background, we generate 10 thousand simulated background light curves with the same duration and binning as the original light curve, and a background rate estimated in an off-pulse interval, measuring the wavelet spectrum for each realization. We then plot the 99% containment interval for each time scale $\delta t$ (blue shaded region) centered on the median (dotted line). In the wavelet power spectrum, Poisson noise follows a power law $W\propto\delta t^{-1}$. This is evident for very short time scales, where the data are dominated by noise. The first time scale that deviates from the noise power law outside the 99% c.l. region represents our estimate of $t_{\rm mv}$. We obtain $t_{\rm mv}\sim 0.3$ s for GRB 100724B and $t_{\rm mv}\sim 0.05$ s for GRB 160509A. We also run the Bayesian Blocks algorithm on the GBM+LLE dataset and confirm that we find the shortest significant structures with a duration of respectively $\sim 0.6$ s and $\sim 0.1$ s, corresponding to $\sim 2~{}t_{\rm mv}$ as expected.

5.2 Delayed Pair-Breakdown in a High-$\sigma$ Relativistic Jet - the $f_{GT}$ model

The model presented by Gill & Thompson (2014) considers the breakout of a strongly magnetized, baryon-poor jet from the confining envelope of a Wolf-Rayet (WR) star at a breakout radius

[TABLE]

The outflow bulk-LF at breakout $\Gamma_{\rm br}=3\Gamma_{\rm br,3}$ is modest and ranges from $\sim 3-10$, and $t_{\rm eng}=1~{}t_{\rm eng,0}~{}{\rm s}$ represents the typical time scale over which the central engine (a black hole in this case) remains active. $\mathcal{R}=10\mathcal{R}_{1}$ is a factor that governs the geometry of the ouflow at deconfinement ($\mathcal{R}<1$ for ‘jet’ geometry and $\mathcal{R}>1$ for ‘pancake’ geometry). It was shown by Thompson & Gill (2014) that the advected quasi-thermal radiation field at breakout has a relatively flat spectrum ($E^{\prime}f^{\prime}_{\rm th}(E^{\prime})\propto E^{\prime 1+\alpha}$, with $\alpha\sim-1$) below the Wien peak at energy $E^{\prime}_{\rm pk,br}\simeq 0.1m_{e}c^{2}$ in the fluid-frame. The enthalpy density of the jet at breakout is dominated by the magnetic field, with compactness $\ell_{B,\rm br}^{\prime}\gtrsim\ell_{\rm th,br}^{\prime}$, where the advected quasi-thermal radiation field has compactness

[TABLE]

and where $\sigma_{T}$ is the Thomson cross-section and $E_{\gamma,{\rm iso}}=10^{54}E_{\gamma,54}$ erg is the total isotropic equivalent energy of the radiation field.

Post jet-breakout, the outflow is accelerated to high bulk-LF $\Gamma\sim 10^{2}-10^{3}$, and the radiation field compactness $\ell^{\prime}_{\rm th}(R)\propto R^{-4}$ and optical depth of the flow $\tau_{T,\pm}(R)\propto R^{-3}$ drop with radius. The thin baryonic layer, that was lifted from the WR envelope during breakout, suffers a corrugation instability (akin to a Rayleigh-Taylor instability) as it feels an effective gravity in its rest frame $g^{\prime}_{\rm eff}=-c^{2}d\Gamma/dR$ due to the acceleration of the outflow. This breaks the baryonic layer into multiple plumes, which lose radiation pressure support at a critical radiative compactness, $\ell^{\prime}_{\gamma,{\rm crit}}\sim m_{p}/Y_{e}m_{e}\simeq 4\times 10^{3}$, where $m_{p}$ is the proton mass and $Y_{e}\approx 0.5Y_{e,1/2}$ is the electron fraction in a long-GRB, and begin to lag behind as the magneto-fluid continues to accelerate. This differential motion between the two components leads to strong inhomogeneities in the magnetofluid and dissipation of the magnetic energy in the form of a turbulent cascade. The dissipation zone is radially localized at

[TABLE]

and the corresponding bulk-LF of the outflow is

[TABLE]

The Thomson depth of the pairs at the dissipation radius is $\tau_{T,\pm,{\rm diss}}\lesssim 10^{-4}$ and the dissipated magnetic energy with compactness $\ell^{\prime}_{\rm heat}$ goes into heating the pairs. The initially relativistically hot pairs inverse-Compton (IC) scatter the peak thermal photons to high energies. As the average energy of the pairs drops, due to pair production, the IC scattered peak gradually moves to lower energies, and finally merges with the thermal peak.

The total radiative compactness of the flow after dissipation can be written as

[TABLE]

where $\xi_{\rm th}=\ell^{\prime}_{\rm th}/\ell^{\prime}_{\rm heat}$ sets the heating compactness $\ell^{\prime}_{\rm heat}$ relative to the thermal compactness $\ell^{\prime}_{\rm th}$. The quasi-thermal soft seed photon spectrum is described as

[TABLE]

where $\alpha$ sets the spectral power-law index below $E^{\prime}_{\rm pk}=3k_{B}T^{\prime}_{\rm th}$, $T^{\prime}_{\rm th}$ is the temperature of the radiation field and $k_{B}$ is the Boltzmann constant. The free parameters of this model are: $\alpha$, $T^{\prime}_{\rm th}$, $\xi_{\rm th}$, and $\ell^{\prime}_{\rm tot}$. The comoving radiation spectrum is then formed using a one-zone time-dependent kinetic code that involves integro-differential equations for both the radiation and particle distributions in the frame of the outflow (see Gill & Thompson, 2014 for further details of the numerical scheme).

To reduce the number of independent model parameters, so as to make the fitting procedure computationally tractable, we set the low-energy power-law index $\alpha$ to that obtained from fitting the $f_{BHec}$ profile. In addition, an estimate of the average comoving radiation field compactness can be obtained from the burst luminosity, such that

[TABLE]

Demanding that $R>R_{\rm diss}$, the above equation can be iterated to determine the correct $\ell^{\prime}_{\rm tot}$ that satisfies the radius constraint. For example, taking the value for the luminosity $L_{\gamma,\rm iso}=2.32\times 10^{53}~{}{\rm erg~{}s}^{-1}$ in the third time-interval of GRB 160509A from Figure 4, we find an emission radius $R=6.7\times 10^{14}~{}{\rm cm}$ for $\Gamma=200$ and $\ell^{\prime}_{\rm tot}=70$. Ideally, $\ell^{\prime}_{\rm tot}$ should remain as an independent parameter since it depends strongly on radius. Redshift information is not available for GRB 100724B, which makes it less trivial to ascertain the correct dissipation radius and $\ell^{\prime}_{\rm tot}$. Therefore, for simplicity, we use the same radiative compactness $\ell^{\prime}_{\rm tot}=70$ for this burst as well, with the underlying assumption being that GRB 100724B had a similar intrinsic brightness as GRB 160509A. This appears to be a reasonable assumption, given the similarity between the two bursts.

5.3 Model Fitting and Results

Next, we fit the two physical models described in the previous subsections to the data. In order to compute the attenuation factor for $f_{BG}$ we have implemented the semi-analytical computation described in Granot et al. (2008) in a code that, in the spirit of reproducible research (Donoho et al., 2009), we make publicly available555https://github.com/giacomov/pyggop.

To fit the $f_{GT}$ model to the data we use templates that are produced by a numerical code which is very computing-intensive, requiring a medium-sized computer farm. We therefore release, in place of the code, the templates that can be used to reproduce our results. These templates are interpolated during the fit procedure to give the final results. As explained in section 5.2 the model has 3 parameters, including the low-energy photon index $\alpha$. The computer code returns differential photon flux as a function of dimensionless rest frame energy. Therefore, in order to fit the data in the observer frame, we need to multiply the dimensionless rest frame energy by $m_{e}c^{2}$ and by a scale factor $\eta$ so that $\Gamma=\eta(1+z)/(m_{e}c^{2})$, where $\Gamma$ is the bulk Lorentz factor. We also of course need a normalization for the model, for a total of 5 free parameters. In order to reduce the number of templates we need to generate, we fix $\alpha$ to the index measured with the $f_{BHec}$ model. Indeed, the fit would converge there anyway since $\alpha$ is the only parameter affecting the spectrum at low energies. This of course does not reduce the number of free parameters of the model, since $\alpha$ is still measured on the data, but it allows us to reduce the complexity of the problem.

In Figure 7 we show the best fit models for $f_{BHec}$, $f_{GT}$ and $f_{BG}$ for both GRBs. Even though the $f_{BG}$ model tends to predict a much higher flux at high photon energies, it is still fully compatible with the data due to the low statistics at high energies. This can be seen in the count spectra shown in Figures 8 and 9. On the other hand, the $f_{GT}$ model is very similar to $f_{BHec}$. The p-values for the goodness of fit test $p$ and $p_{LAT}$ for all three models, computed as described in section 3.2, are $\gtrsim 0.05$ for all intervals. We conclude that all 3 models appear to describe our data well. For high compactness the $f_{\rm GT}$ model features a fairly prominent pair annihilation line, visible for example in the best fit models for GRB 160509A (middle right panel of Figure 7). It is currently not detectable by Fermi/LAT as it is smeared out by energy dispersion effects in the detector, and indeed it is not apparent once the model is folded with the response of the instrument (Figure 9, blue dashed line).

The best fit parameters of the $f_{BG}$ model for both GRBs are shown in Figures 10 and 11. The values of $\Delta R/R_{0}$ are of order unity (typically $\sim 1-2$ and ranging from $\sim 0.5$ to $\sim 3.5$) for GRB 100724B, and somewhat larger, typically around $\sim 2-3$ (ranging between $\sim 1.5$ and $\sim 5.5$) for GRB 160509A. This is in reasonable agreement with the expectations of the internal shocks model (for which the physical setup of the Granot et al. (2008) model is particularly well suited), as are the typical inferred emission radii ($\sim 10^{13}-10^{15}\;$cm for GRB 100724B and $\sim 10^{13}-10^{14}\;$cm for GRB 160509A). The LF is relatively low when compared to what is inferred for other LAT-detected GRBs (Ackermann et al., 2013b), ranging from $\sim 140$ to $\sim 400$ for GRB 160509A, while for GRB 100724B it depends on the unknown redshift but for typical redshifts it is broadly similar (ranging from $\sim 70$ to $\sim 310$ for $1\leq z\leq 4$). This may account for the relatively low values of the cutoff energy $E_{c}$ (e.g. as inferred for the $f_{BHec}$ model and is shown in Figures 2 and 4), of $\sim 15-50\;$MeV for GRB 100724B (except in the last time bin, where $E_{c}\sim 200\;$MeV), and $\sim 100\;$MeV (ranging between $\sim 20\;$MeV and $\sim 400\;$MeV) for GRB 160509A. In turn, this may demonstrate the fact that slower GRBs tend to be dimmer in the LAT energy range, thus producing a selection effect in favor of faster GRBs in the LAT GRB sample. This effect would be more pronounced when not accounting for LAT-LLE only detections (with no photons detected above $\sim 100\;$MeV). Finally, the self consistency of the $f_{BG}$ (and Granot et al. 2008) model requires that $\tilde{\tau}_{T,\pm}<1$ and therefore $\Gamma_{0}/\Gamma_{\rm max}=\Gamma_{\rm\gamma\gamma,min}(E_{c})/\Gamma_{\rm\gamma\gamma,max}(E_{c})<1$. This is satisfied, at least marginally, in all time bins (with the possible exception of time bin 3 in GRB 100724B). We conclude that $f_{BG}$ is a viable interpretation for both GRBs.

The best fit parameters of the $f_{GT}$ model for both GRBs are shown in Figure 12.

The bulk-LF of the outflow in the $f_{GT}$ model was determined by laterally shifting the comoving spectrum in energy by the factor $\Gamma/(1+z)$ to fit the observed spectrum. We show the temporal evolution of $\Gamma$ for both GRBs in Figure 12, where the fluctuations in $\Gamma$ are correlated with fluctuations in the observed flux, and therefore luminosity, for both bursts. In the case of GRB 160509A, this behavior clearly coincides with the two broad peaks observed in the BGO and LLE emission. Since there is no redshift available for GRB 100724B, $\Gamma$ could only be determined in the engine-frame. Assuming a typical redshift of $z\sim 2$ we obtain $\Gamma\sim 90-270$ throughout the entire prompt emission phase. For GRB 160509A, $\Gamma$ varies by a factor $\sim 5$ during the prompt phase and peaks at $\Gamma\lesssim 500$. In light of the fact that the underlying numerical model is one-zone, such a high value for $\Gamma$ is typically found from one-zone estimates (e.g. Lithwick & Sari, 2001) as compared to that obtained from models of Granot et al. (2008) and Hascoët et al. (2012).

In the $f_{GT}$ model, as the outflow expands to larger radii, the comoving temperature of the quasi-thermal radiation field should drop due to adiabatic cooling. This behavior is clearly seen in the evolution of $\theta_{\rm th}$ in the case of GRB 160509A; the existence of a similar trend is less clear for GRB 100724B.

The appearance of a quasi-thermal spectrum at smaller radii and, consequently, larger $\theta_{\rm th}=k_{B}T_{\rm th}/m_{e}c^{2}$ is quite naturally explained in the $f_{GT}$ model. Such an emission, with no high-energy component, is expected to escape from optically thin regions of the outflow before any dissipation has occurred. Since it originates at smaller radii, it should arrive at the observer earlier than the main burst. Thus, we associate this quasi-thermal component to the precursors observed in both GRBs. In particular, the precursor of GRB 100724B can be described well with the quasi-thermal spectrum predicted by the GT model $f_{\rm th}$ (in the observer-frame) in eq. (8), yielding $\alpha=-0.79^{-0.17}_{+0.20}$ and $kT=43^{-5}_{+7}\;$keV. The precursor of GRB 160509A can be similary fit, with best fit parameters $\alpha=-1.13\pm 0.07$ and $kT=87\pm 21\;$keV.

5.4 Comparison to other Fermi/LAT GRBs

The most striking property of GRBs 100724B and 160509A is the clear need for a high-energy spectral cutoff in their prompt emission with respect to the extrapolation of the low-energy component. For GRB 100724B the cut-off energy in its time-resolved spectrum typically lies in the range $E_{c}\sim 20-60\;$MeV with high statistical significance, and in the case of GRB 160509A the cutoff typically appears at energies $E_{c}\sim 80-150\;$MeV.

In earlier LAT GRBs, for example GRB 080825C, there was marginal evidence for a cutoff at an energy of $E_{c}\sim 1.8\;$MeV (Abdo et al., 2009a), which if true does not have a good natural explanation. In GRB 090926A (Ackermann et al., 2011) there was a high-energy spectral cutoff at $E_{c}\sim 1.4\;$GeV in the time-integrated spectrum, and at $E_{c}\sim 0.4\;$GeV in one time bin of the time-resolved spectrum, which has been nicely interpreted as arising due to intrinsic opacity to pair production in the source, in which case it implies a bulk-LF of $\Gamma\sim 300-700$ for the prompt emission region, depending on the exact model assumptions about the emission.

The upper end of this range corresponds to a simple one-zone model in which the radiation in the outflow’s frame is uniform, isotropic and time-independent. In this case, if the photon number spectrum is described by a power-law for photon energies $\varepsilon_{\rm pk}<(\varepsilon\equiv E/m_{e}c^{2})<\varepsilon_{c}$, such that $f(\varepsilon)=f_{\rm pk}(\varepsilon/\varepsilon_{\rm pk})^{\beta}$, where $\varepsilon_{\rm pk}$ and $\varepsilon_{c}$ are, respectively, the dimensionless peak and cutoff energies, then an estimate of the bulk-$\Gamma$, which corresponds to the condition that $\tau_{\gamma\gamma}(\varepsilon>\varepsilon_{c})>1$, is given by Lithwick & Sari (2001, LS01 hereafter; eq. (5) therein),

[TABLE]

where the numerical values are for $\beta=-2.2=-11/5$, which is typical of the values measured for the prompt GRB. The above estimate also assumes a redshift of $z=2$, luminosity distance $d_{L}=4.8\times 10^{28}d_{L,28.7}$ cm, variability time $t_{\rm mv}=10^{-1}t_{\rm mv,-1}$ s, peak photon number flux $f_{\rm pk}=10^{7}f_{\rm pk,7}~{}{\rm cm}^{-2}~{}{\rm s}^{-1}~{}{\rm erg}^{-1}$ at a peak energy $E_{\rm pk}=250\varepsilon_{\rm pk,-0.3}$ keV, and cutoff energy $E_{c}=100\varepsilon_{c,2.3}$ MeV.

We use the analytic one-zone method of LS01 to calculate the bulk-$\Gamma$ for the case of GRB 160509A and compare it with $\Gamma$ obtained from the $f_{BG}$ and $f_{GT}$ model in Figure 13. The right panel of Figure 13 shows how the LS01 method of estimating bulk-$\Gamma$ generally yields a value that is a factor $1.5-2$ times higher than that given by a fully time-dependent model (which yields the lower limit on $\Gamma$ in the case of GRB 090926A) where the radiation field starts from zero at the emission onset and is calculated self-consistently as a function of time, space and direction (Granot et al., 2008; Hascoët et al., 2012; Gill & Granot, 2018). In GRB 110731A there is a (slightly marginal) detection of a cutoff at $E_{c}\sim 0.4\;$GeV, which similarly implies $\Gamma\sim 300-600$ if interpreted as due to intrinsic pair production in the source (Ackermann et al., 2013a).

In Figure 14 we compare the bulk-$\Gamma$ estimates obtained from the $f_{BG}$ and LS01 models for several GRBs. Since GRB 100724B lacks redshift information, we show the evolution of $\Gamma$ with redshift. Other Fermi/LAT detected GRBs that do not show any spectral cutoff, namely GRBs 080916C and 090510, are also shown. However, for these $\Gamma$ should be interpreted as a lower limit $\Gamma_{\gamma\gamma,\rm min}$.

There are large differences in the observed cutoff energies $E_{c}$ between different GRBs (see, e.g., Table 4). To understand which properties among the different GRBs are leading to different cutoff energies (or lack thereof, in which case the highest energy observed photon was used), we can express eq. (2) (which relies on Eq. (126) of Granot et al. 2008; see footnote 4) in terms of the intrinsic parameters, such that the cutoff energy in the central engine frame (quantities in this frame are expressed with a subscript $z$) is

[TABLE]

where $t_{v,z}=t_{v}/(1+z)$. We compare each given GRB (subscript ‘$i$’) with GRB 160509A (in particular the results of time-bin 6; subscript ‘[math]’) where we quantify the effect of a change in the parameter $\xi=\{\Gamma,~{}\beta,~{}L_{0,52},~{}t_{v,z},~{}z\}$ as follows (see Table 4),

[TABLE]

where $E_{c,z,0}(\xi_{i})$ is the value of $E_{c,z}$ for GRB 160509A obtained from eq. (5.4) by replacing the parameter $\xi$ by its value for GRB $i$ and keeping all other parameters fixed to their measured values for GRB 160509A. In GRBs where a cutoff was not observed, and instead only a lower limit on $E_{c}$ was derived, which was used to derive a lower limit on $\Gamma$ (GRBs 080916C and 090510), we use these lower limits for this comparison.

When examining the origin of the large differences in the observed cutoff energies $E_{c}$ between different GRBs, we find that the effect due to differences in redshifts between the different GRBs is sub-dominant, as can be seen in Table 4. This further implies that typically most of the differences between the observed break energies $E_{c}$ are intrinsic (i.e. the differences in $E_{c,z}$).

One might naively expect that the dominant intrinsic parameter that would account for the different $E_{c,z}$ and $E_{c}$ values would be the Lorentz factor, $\Gamma$, since in eq. (5.4) it appears with a power larger than that of $t_{v,z}/L_{0,52}$ by a factor of $2-2\beta\sim 6$ (for typical values of $\beta\sim-2$). However, out of the three GRBs we considered for comparison here, which have a disparate set of intrinsic properties, only in one of them, GRB 090510, does it appear to be dominant (by a factor of $\sim 2$) over the effects of the differences in $t_{v,z}$ and $L_{0,52}$, while in GRB 090926A it has a comparable effect and in GRB 080916C it is sub-dominant. This occurs since despite the larger power of $\Gamma$ in the expression for $E_{c}$, it varies by a smaller factor between different GRBs compared to $t_{v,z}$ and $L_{0,52}$. The dependence of $E_{c}$ and $E_{c,z}$ on $\beta$ is non-trivial (see eq. (5.4)), but the effect of its variation between different GRBs is typically comparable to that of the other physical parameters. The comparison of the intrinsic properties of different GRBs in Table 4 ultimately shows that it is likely that differences in many of these properties jointly contribute to the appearance of high-energy spectral cutoffs (or lack thereof). This obviously constitutes a broad set of possibilities, and it is not unlikely to find GRBs with spectral cutoffs where only the difference in bulk-$\Gamma$ makes the dominant contribution.

Before GRB 100724B (Ackermann et al., 2013b) there was no clear direct evidence for a high-energy spectral cutoff at an energy $E_{c}\ll 1\;$GeV. The prompt GRB spectrum of some GRBs is consistent with a Comptonized spectrum featuring a power-law with an exponential cutoff, but with a typical peak energy $E_{p}\lesssim 1\;$MeV, so this cannot really be considered as a high-energy cutoff, and most likely has a different physical origin. On the other hand, there was indirect evidence for a high-energy cutoff at tens of MeV from the extrapolation of the GRB spectrum and LAT upper-limits (Beniamini et al., 2011; Guetta et al., 2011; Fermi Large Area Telescope Team et al., 2012). Therefore, GRB 100724B shows the first clear-cut detection of a high-energy cutoff at well below a GeV. Other Fermi-detected GRBs were shown to have similar sub-GeV cutoffs by Tang et al. (2015).

6 Discussion

GRBs must be compact sources in order for their short variability time scales to be explained, and at the same time they must be optically thin above the pair production threshold in order to produce the observed high-energy $\gamma$-rays. This is the so-called “compactness” problem (Piran, 1999), which in the fireball model is solved by assuming a high bulk Lorentz factor $\Gamma$ for the emitting shells. However, $\Gamma$ cannot be arbitrarily large; thus the system will become optically thick above a certain energy, producing a spectral cutoff. The observed high-energy spectral cutoffs in both GRBs (100724B and 160509A) can therefore be interpreted this way. Both theoretical models considered in this paper ($f_{BG}$ and $f_{GT}$) produce a high-energy spectral cutoffs due to pair production opacity, and they describe the data well. Still, it is prudent to ask if instead the spectral cutoffs can be explained due to some intrinsic limitation of the emission process in producing high-energy photons, which would introduce a cutoff at an energy different (and possibly lower) than the cutoff expected from the pair-production opacity.

One such case, for example, would be if the underlying mechanism for the prompt emission was synchrotron. Then an exponential cutoff would be observed at fluid-frame energies $E^{\prime}>E^{\prime}_{\rm syn,max}=\gamma_{e,\rm max}^{2}(B^{\prime}/B_{Q})m_{e}c^{2}$, where $\gamma_{e,\rm max}m_{e}c^{2}$ is the maximum energy to which electrons are accelerated in the dissipation region, $B^{\prime}$ is the local magnetic field in the fluid frame, and $B_{Q}=m_{e}^{2}c^{3}/e\hbar$ is the quantum critical field. In this case, $\gamma_{e,\rm max}$ will be limited if the synchrotron cooling time of electrons, $t^{\prime}_{syn}=6\pi m_{e}c/(\sigma_{T}B^{\prime 2}\gamma_{e})$, is shorter than their acceleration time. The shortest viable length scale over which electrons are accelerated is given by their Larmor radius, $r_{L}^{\prime}=\gamma_{e}m_{e}c^{2}/(eB^{\prime})$, which yields an acceleration time $t_{\rm acc}^{\prime}=r_{L}^{\prime}/c$. Comparison of the two timescales gives the factor $\gamma_{e,\rm max}^{2}B^{\prime}=6\pi e/\sigma_{T}$, where $\sigma_{T}=(8\pi/3)(e^{2}/m_{e}c^{2})^{2}$ is the Thomson cross-section. From this we immediately find that the maximum energy of synchrotron photons is $E^{\prime}_{\rm syn,max}=(9/4\alpha_{F})m_{e}c^{2}$, where $\alpha_{F}=e^{2}/\hbar c=1/137$ is the fine-structure constant (e.g. Guilbert et al., 1983; de Jager & Harding, 1992; Piran & Nakar, 2010; Atwood et al., 2013). In the observer frame, this limiting energy translates to

[TABLE]

It is much higher than the cutoffs observed in the two GRBs discussed in this work. In other GRBs photons approach or even exceed $E_{\rm syn,max}$, which suggests very efficient electron acceleration. Since the efficiency of electron acceleration is not expected to drastically change between different GRBs, this would support a different origin for the high-energy cutoffs in GRBs 100724B and 160509A, in agreement with our interpretation of an intrinsic opacity to pair production origin.

In most GRBs there is no observed high-energy cutoff with respect to the extrapolation of the low-energy component, so the maximal observed photon energy $E_{\rm max}$ is used as a lower limit for any possible cutoff energy, which in turn sets a lower limit, $\Gamma_{\gamma\gamma,{\rm min}}$, on the LF of the emitting region, $\Gamma$, through the condition that $\tau_{\gamma\gamma}(E_{\rm max})<1$. In cases where a high-energy cutoff is actually observed at an energy $E_{c}$ and may be attributed to intrinsic pair production, then $\Gamma_{\gamma\gamma,{\rm min}}$ can serve as an actual estimate of $\Gamma$ when $E_{\rm max}$ is replaced by $E_{c}$, as shown in eq. (5.4). For the $f_{BG}$ model $\Gamma_{\gamma\gamma,{\rm min}}(E_{c})$ is given by eq. (2).

This is valid only as long as $\Gamma_{\gamma\gamma,{\rm min}}$ is lower than

[TABLE]

for which $E_{c}$ corresponds to a comoving photon energy of $E^{\prime}_{c}>m_{e}c^{2}$, so that photons near the cutoff can pair-produce with other photons of comparable energy. When the effects of pair cascades on the spectral cutoff are ignored, then for a given $\Gamma$ the cutoff energy always satisfies

[TABLE]

as long as the high-energy photon index is lower than $-1$ ($\beta<-1$ where $dN_{\rm ph}/dE\propto E^{\beta}$, as is almost always the case in GRB prompt spectra), so that $dN_{\rm ph}/d\log E=E(dN_{\rm ph}/dE)$ decreases with photon energy $E$. Here $E_{sa}$ is the minimal energy of photons that can “self-annihilate”, i.e. interact with other photons of the same energy. This occurs because of the following reason. Let us denote by $E_{1}$ the photon energy above which the optical depth to pair production is large, $\tau_{\gamma\gamma}(E_{1})=1$, and also denote by $\varepsilon=E/m_{e}c^{2}$ and $\varepsilon^{\prime}=\varepsilon(1+z)/\Gamma=E/E_{sa}$ the observer and comoving frame photon energies in units of the electron rest energy.

If $E_{1}>E_{sa}\Leftrightarrow\varepsilon^{\prime}_{1}>1\Leftrightarrow\Gamma_{\gamma\gamma,{\rm min}}(E_{c})<\Gamma_{\gamma\gamma,{\rm max}}(E_{c})$ then the Thomson opacity of the $e^{\pm}$ pairs that are produced is small (Gill & Granot, 2018),

[TABLE]

and can therefore be ignored for our purposes. Therefore, we call this the “thin” regime. In this regime there will be a cutoff at $E_{c}=E_{1}>E_{sa}\Leftrightarrow\varepsilon^{\prime}_{c}=\varepsilon^{\prime}_{1}>1$.

However, if $E_{1}<E_{sa}\Leftrightarrow\varepsilon^{\prime}_{1}<1\Leftrightarrow\Gamma_{\gamma\gamma,{\rm min}}(E_{c})>\Gamma_{\gamma\gamma,{\rm max}}(E_{c})$ then the situation changes. In this regime photons in the energy range $E_{1}<E<E_{sa}\Leftrightarrow\varepsilon^{\prime}_{1}<\varepsilon^{\prime}<1$ will on the one hand have a large initial optical depth to pair production, $\tau_{\gamma\gamma}(E)\sim(E/E_{1})^{-1-\beta}=(\varepsilon^{\prime}/\varepsilon^{\prime}_{1})^{-1-\beta}>1$. However, on the other hand since they can only pair-produce or annihilate with photons of energy $\varepsilon^{\prime}_{\rm an}\geq 1/\varepsilon^{\prime}>1>\varepsilon^{\prime}$, and there are fewer such photons, they all quickly annihilate, so that the optical depth $\tau_{\gamma\gamma}(E)$ rapidly drops below its initial value to well below unity, and most of the photons in this energy range remain with no photons that they can pair-produce with. Therefore, only photons of energy $E>E_{sa}\Leftrightarrow\varepsilon^{\prime}>1$ can fully annihilate, and there is a cutoff only at $E_{c}=E_{sa}>E_{1}\Leftrightarrow\varepsilon^{\prime}_{c}=1>\varepsilon^{\prime}_{1}$. The latter implies that in this regime $\Gamma=\Gamma_{\gamma\gamma,{\rm max}}(E_{c})$, i.e. the bulk-LF of the emission region is given by eq. (14). In this regime the Thomson optical depth of the pairs that are produced (neglecting their possible annihilation, hence the tilde in its notation) is large (Gill & Granot, 2018),

[TABLE]

Therefore, we call this the “thick” regime.

Since typically $\beta\sim-2$, $\tilde{\tau}_{T,\pm}$ in the thick regime scales as a fairly high power ($\sim 6$) of the ratio $\Gamma_{\gamma\gamma,{\rm min}}(E_{c})/\Gamma_{\gamma\gamma,{\rm max}}(E_{c})$, so even if $\Gamma_{\gamma\gamma,{\rm min}}(E_{c})$ exceeds $\Gamma_{\gamma\gamma,{\rm max}}(E_{c})$ only by a factor of a few one might already have $\tilde{\tau}_{T,\pm}\gtrsim 10^{2}$. In this case, the spectrum is modified due to Compton scattering by cold pairs, which also brings the cutoff energy below $E_{sa}$ (Gill & Granot, 2018). In addition, large $\tilde{\tau}_{T,\pm}$ decreases the radiative efficiency due to adiabatic losses and washes away much of the temporal variability, as photons must diffuse out of the emission region on a diffusion time larger than the dynamical time. On the other hand, in this regime pair annihilation can become important, and the expansion of the emitting region also dilutes its opacity and causes a non-isotropic photon distribution in the comoving frame, which further suppresses pair production.

Altogether the results for the LF and cutoff energy in the two regimes discussed above can be summarized by

[TABLE]

[TABLE]

6.1 $f_{GT}$ vs $f_{BG}$ Model

Spectrally, the one important place where the $f_{GT}$ model differs from the $f_{BG}$ model is where the cutoff energy lies with respect to $m_{e}c^{2}$. If the comoving radiation field compactness is high the $f_{GT}$ model always yields the comoving frame cutoff energy $E_{c}^{\prime}<m_{e}c^{2}$, such that $\Gamma_{\gamma\gamma,\rm min}(E_{c})>\Gamma_{\gamma\gamma,\rm max}(E_{c})$. As argued above, this consequently would yield a high pair Thomson depth and significantly alter the observed spectrum and temporal variability. However, since pair-production and pair-annihilation effects are self-consistently accounted for in the $f_{GT}$ model, the $e^{\pm}$ pair Thomson depth is always regulated to $\tau_{T,e^{\pm}}\sim 1-5$ in the dissipation region. What the $f_{GT}$ model does not account for is the pair opacity accumulated over the line of sight as the photon travels from its emission point to the observer while interacting with other photons en route. This is the essence of the $f_{BG}$ model. Still, this additional pair opacity effect will not significantly alter the spectrum obtained in the $f_{GT}$ model as the high-energy spectrum is already exponentially suppressed due to pair-production in the dissipation region. As a result, the cutoff energy cannot be made appreciably smaller due to additional pair-opacity effects.

When comparing the two spectral models, we find that the $f_{BG}$ model yields $\Gamma$ values that are on average comparable to that obtained from the one-zone $f_{GT}$ model. Both models have additional parameters, other than the ones used for fitting in this work, that introduce some degeneracy in the final outcome.

A potential test for both models is that photons above the cutoff energy $E_{c}$ are expected to arrive preferentially near the beginning of pulses in the lightcurve, as compared to near their peak or during their tails. However, this requires good photon statistics within a single spike of the lightcurve, which was not available so far with Fermi/LAT but may become possible in the future with the Cherenkov Telescope Array (CTA; e.g. Inoue et al. 2013).

6.2 Comparison with Other Work

In a recent work, using Swift X-ray data along with ground-based optical, infrared and radio data, Laskar et al. (2016) analyzed the afterglow emission and determined $\Gamma(t_{\rm dec})\approx 330$ for GRB 160509A at the deceleration time $t_{\rm dec}\approx 460\;$s for a constant density circumburst environment ($k=0$). In comparison, we find a mean bulk-LF of $\Gamma_{\rm BG}\sim 220$ from the $f_{BG}$ model and $\Gamma_{GT}\sim 215$ from the $f_{GT}$ model over the entire duration of the prompt phase. The apparent discrepancy between our results and that of Laskar et al. (2016) critically depends on the density profile of the circumburst medium, $\rho_{\rm ext}\propto R^{-k}$. They find a much lower $\Gamma(t_{\rm dec})=34$ (and $t_{\rm dec}\approx 170\;$s) for the wind case ($k=2$), where the density of the surrounding medium is determined by stellar winds from the progenitor star. The actual value of $\Gamma$ may be somewhere in the middle depending on the value for $k$, which is likely also intermediate ($0<k<2$), as suggested by some afterglow modelings (e.g., Kouveliotou et al., 2013), and is viable given the uncertain wind velocity and mass loss rate history at the massive star progenitor’s last years (which determine the density profile around the deceleration radius corresponding to the afterglow onset). In that case, the results of this work would be consistent with that obtained from the multi-wavelength afterglow analysis.

Moreover, the effective duration of the prompt emission in GRB 160509A is $\sim 20-30\;$s (see Fig. 3), i.e. much shorter than its $T_{90}\approx 370\;$s, which is dominated by a weak and very soft emission episode around $T\sim 300-400\;$s. This may suggest either an earlier deceleration time, $t_{\rm dec}\sim 30\;$s for a relativistic reverse shock (or “thick shell”), or alternatively if $t_{\rm dec}\gg 30\;$s the one would expect a Newtonian reverse shock (or “thin shell”), in which case the correspondingly weaker reverse shock would tend to imply a lower value for $\Gamma(t_{\rm dec})$, which could be consistent with the values we derived for models $f_{BG}$ and $f_{GT}$.

Finally, for strong internal shocks a good part of the outflow energy may reside in internal energy of the baryons just after the shells collide. It eventually transforms back to kinetic energy above the internal shock emission radius leading to a larger asymptotic $\Gamma(t_{\rm dec})$ compared to $\Gamma$ of the emitting plasma during the internal shocks themselves. A similar effect may arise in a Poynting-flux-dominated outflow if the emission occurs during the acceleration phase.

In the work of Tang et al. (2015), a total of eight GRBs (including GRB 100724B) that were detected by Fermi were found to have spectral cutoffs between tens of MeV and several 100 MeV. They derived the bulk-$\Gamma$ for these GRBs using a simple one-zone analytical model, akin to the LS01 model, and found that for majority of the cases $\Gamma>\Gamma_{\gamma\gamma,\rm max}$. This led them to estimate the actual bulk-LF by its maximum value given by $\Gamma_{\gamma\gamma,\rm max}$. As shown by Gill & Granot (2018), estimating $\Gamma$ in this way can lead to underestimating its true value by as much as an order of magnitude, since in this case the spectral break energy is modified by pair cascades. It is clear from Figure 14 that $\Gamma$ obtained from simple one-zone analytic models will exceed $\Gamma_{\gamma\gamma,\rm max}$ in the case of GRB 100724B (unless $z\ll 1$ which is unlikely), whereas the much more detailed and self-consistent model of Granot et al. (2008) generally yields $\Gamma\lesssim\Gamma_{\gamma\gamma,\rm max}$ for all redshifts.

Recently, a sub-photospheric dissipation model of Pe’er et al. (2006) was used in the work of Ahlgren et al. (2015) to fit the time-resolved spectra of GRB 100724B using the code developed in Pe’er & Waxman (2005). The underlying GRB model producing the prompt-phase spectrum has many similarities with the model of Gill & Thompson (2014), in particular the continuous and slow-heating of electrons which then Compton up-scatter the soft thermal emission to produce the spectrum above the peak. The major difference between the two models is that the model of Gill & Thompson (2014) assumes a Poynting-flux-dominated baryon-pure outflow whereas the model advanced in Pe’er et al. (2006) assumes a kinetic energy dominated baryonic jet. They also find a much larger mean $\Gamma\approx 443$ for GRB 100724B, while assuming a redshift $z=1$, in comparison to $\Gamma_{GT}\approx 180$ and $\Gamma_{BG}\approx 10^{2}$ using the two models considered in this work. More importantly, the model fit in Ahlgren et al. (2015) lacks a spectral cutoff at high energies as sharp as the model considered in this work. Consequently, it yields a poorer fit in the 1 MeV to 1 GeV energy range (compare the upper right panel in Figure 2 in Ahlgren et al. (2015) to Fig. 8). A spectral cutoff is naturally and self-consistently produced in both the $f_{GT}$ and $f_{BG}$ models.

7 Conclusions

We presented a detailed time-resolved analysis of two bright Fermi/LAT GRBs, GRB 100724B and GRB 160509A, that provide the clearest examples of sub-GeV high-energy cutoffs during the prompt emission. We characterized phenomenologically the high-energy cutoffs, which we measure respectively in the range 20-60 MeV and 80-150 MeV. We have shown through the fitting of two physical models that the observed cutoff can be interpreted as the result of intrinsic opacity to pair production at the source, while it appears to be too low to be explained as originating from the limitation of the particle acceleration process.

In particular, a semi-phenomenological model of an impulsive relativistic outflow with detailed $\gamma\gamma$ opacity computation presented in Granot et al. (2008) can describe the data well and self-consistently, yielding estimates for the emission onset radius $R_{0}\sim 10^{13}-10^{14}\;$cm and for the fractional size of the emission zone $\Delta R/R_{0}\sim 1-5$ that are consistent with the internal shocks model. The one-zone photospheric model of Gill & Thompson (2014) can also describe the data well. Moreover, it predicts a drop in the comoving temperature of the seed quasi-thermal radiation field which is clearly observed in GRB 160509A (but not in GRB 100724B as the details of the model depend on the redshift, which is lacking in this case). The estimate for the bulk Lorentz factors derived by using the model of Granot et al. (2008) are typically in the range $\Gamma_{0}\sim 100-300$, and they are on average comparable to those obtained from the one-zone photospheric model of Gill & Thompson (2014). These estimates are a factor of a few to several smaller than the lower limits derived for bright LAT GRBs and a factor of $\sim 2$ smaller than values inferred from high-energy cutoffs, which were generally obtained for LAT-detected GRBs from a one-zone analytical model (see for example Tang et al., 2015). Indeed such a factor of $\sim 2$ difference exists also when deriving $\Gamma$ for the same GRB using different models (see, e.g., Figures 13 and 14).

Because of opacity to intrinsic pair production, slower GRBs tend to be fainter in the LAT energy range and are therefore more difficult to detect. This may produce a selection bias against deriving lower Lorentz factors from the detection of high-energy cutoffs. We also note that our measurement for $\Gamma_{0}$ is in line with the upper limit estimated in Nava et al. (2017) for GRBs observed but not detected by the LAT.

We find that the differences in observed cutoff energies $E_{c}$ between different GRBs are predominantly intrinsic, and arise not only from the different Lorentz factor $\Gamma$ of their emission regions, but also from differences in other intrinsic parameters, namely their variability times $t_{v,z}$, isotropic equivalent luminosities $L_{0,52}$, and high-energy photon index $\beta$.

The two GRBs analyzed in this work have relatively low inferred Lorentz factors compared to other Fermi/LAT GRBs. They were still detected by Fermi/LAT despite their relatively low cutoff energies of $E_{c}\lesssim 100\;$MeV, since they are extremely bright at $\lesssim\;$MeV energies. This may introduce a bias in the Fermi/LAT GRB sample against GRBs with low Lorentz factors $\Gamma$, as well as short variability times (corresponding to small emission radii), as these would lead to low cutoff energies $E_{c}$, which would make them more difficult to detect with Fermi/LAT. $E_{c}$ also decreases as the isotropic equivalent luminosity ($L_{0,52}$) increases, so that highly luminous GRBs would require a higher Lorentz factor in order to be detected by Fermi/LAT. This may introduce an apparent positive correlation between the isotropic equivalent luminosity $L_{\rm iso}$ and $\Gamma$, such that $\Gamma\propto L_{\rm iso}^{1/(2-2\beta)}$ with all else being equal. A positive correlation between $\Gamma$ and $L_{\rm iso}$ has indeed been claimed in the literature (e.g. Lü et al., 2012). The possible apparent correlation we point out is not expected to be very tight, and is not expected to appear in the time-resolved spectroscopy of a single GRB (in which such a correlation would most likely be of intrinsic origin). This correlation may be modified by the fact that more luminous GRBs may be detected for a slightly lower $E_{c}$ with possible correlations with $\beta$ or $t_{v}$.

The Fermi/LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l’Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Études Spatiales in France. This work performed in part under DOE Contract DE- AC02-76SF00515. JG and RG acknowledge support from the Israeli Science Foundation under Grant No. 719/14. RG is supported by an Open University of Israel Research Fund.

Appendix A High-energy detection and localization of GRB 100724B

Shortly after the beginning of the prompt emission, due to the Fermi satellite’s orbital motion the location of GRB 100724B on the sky started to move closer and closer to the Earth Limb (EL), which is a very powerful source of gamma rays. Given the altitude of Fermi, the EL appears in LAT standard data as a curved band centered around $\eta=113\,^{\circ}$, where $\eta$ is the angle with the zenith of the spacecraft (Zenith angle). Thus, at a given energy $E$ there is potentially significant contamination from the EL for all regions of the field of view at a zenith angle $\eta\gtrsim 113\,^{\circ}-PSF(\theta,E)$, where $PSF(\theta,E)$ is the size of the PSF at off-axis angle $\theta$ and at the energy $E$ as measured for example by the 90% containment radius. Due to the orbital motion of the Fermi spacecraft, the zenith angle $\eta(\vec{p})$ of a fixed point in the sky $\vec{p}$ is continuously changing. The procedure suggested by the LAT team to limit the contamination from the Earth Limb proceeds as follows. Let us fix a Region Of Interest (ROI) in sky coordinates, centered around the position $\vec{p}_{s}$ of the source, and with a radius $R$: there are time intervals in which such ROI is clean from EL contamination (good time intervals, GTIs) and other intervals in which it is not (bad time intervals, BTIs). Let us fix a minimum energy for our analysis $E_{\rm min}$. The GTIs are all the time intervals in which:

[TABLE]

The computation of such GTIs is performed by the tool gtmktime, part of the Fermi Science Tools. The analysis is then performed only on the GTIs.

We first search for high-energy emission from GRB 100724B during the prompt emission in a circle of $15^{\circ}$ around the location provided by Fermi/GBM. Given the high zenith angle of the source, the condition in eq. (A1) does not return any GTI unless we increase $E_{\rm min}$ up to 300 MeV. There are only 2 photons above that energy, and the source is not detected. However, high-energy emission from GRBs has been proved to last much longer than the prompt emission (Ackermann et al., 2013b). We therefore consider a longer time interval, covering up to $15\;$ks after the beginning of the prompt emission, and adopt $E_{\rm min}=1\;$GeV to recover exposure. The condition in eq. (A1) returns now 10962 s of GTIs in the first $15\;$ks after $t_{0}$. In this case 6 events survive the cuts, among which a $\sim 10\;$GeV photon at $t_{0}+2239\;$s. Note that there are still no photons during the prompt emission with this selection. We detect and localize the source producing a Test Statistic (TS) map for this time interval, with a likelihood model containing the appropriate Galactic and isotropic templates provided by the LAT collaboration as well as all the point sources from the 3FGL catalog (Acero et al., 2015). We model GRB 100724B as a point source with a power-law spectrum, and we keep the parameters for all 3FGL sources fixed to their catalog values. No 3FGL source is detected in the small time window of our analysis, thus this choice is irrelevant for our final results.

This analysis results in a firm detection of the GRB with a significance of $\sim 7\sigma$. The best fit power law has a photon index of $\alpha=-1.8\pm 0.4$, a value typical for the high-energy emission of GRBs (Ackermann et al., 2013b), with an average flux of $8.26\times 10^{-10}$ erg cm*-2* s*-1* (1-100 GeV). Our localization is shown in Figure 15. The cross marks the best fit position, corresponding to $R.A.=123.47\,^{\circ}$ and $Dec.=75.88\,^{\circ}$ (J2000), while the white lines correspond to the 68% and 90% c.l. containment regions. We used this position in the following analysis. The most accurate localization available in the literature before this work was reported in the Fermi/LAT GRB Catalog (Ackermann et al., 2013b), and corresponds to R.A. = $119.59\,^{\circ}$, Dec. = $75.86\,^{\circ}$ (J2000) with a 68% containment radius of $0.88\,^{\circ}$ (dashed line in Figure 15). This localization was based on the detection during the prompt phase, obtained by relaxing the zenith angle cut using a threshold of $110\,^{\circ}$, and it is compatible with the one found in this work, but it features a much larger containment radius. The excess we find is not related to any known source, and lies within the GBM localization region for GRB 100724B. The source is not detected in any other time interval. We therefore identify it as the high-energy counterpart of the burst.

Appendix B Likelihood Ratio Test

The Likelihood Ratio Test (LRT) is a procedure for model selection. It is a statistical test that helps to select between two nested models $m_{0}$ and $m_{1}$ the one that best captures all significant features of the data. The model $m_{0}$ is the null-hypothesis, and represents the simpler model, while $m_{1}$ is the alternative hypothesis and is more complex than $m_{0}$. The Test Statistic (TS) is twice the difference in the log-likelihood $S$ between the two models. Wilks’ theorem states that under certain assumptions $TS$ is asymptotically distributed as a $\chi^{2}$ with $n$ degrees of freedom, where $n$ is the difference in number of parameters between $m_{0}$ and $m_{1}$. In the cases of interest here such theorem is not guaranteed to hold, as our setup violates some of its hypotheses (see Protassov et al., 2002, for details). Also, we might not be close enough to the asymptotic regime. We have then to rely on Monte Carlo simulations to calibrate the distribution of $TS$. In particular, we perform 1 million simulations of the null hypothesis $m_{0}$ using the responses of the instruments. Then, we fit each simulated dataset with both $m_{0}$ and $m_{1}$, recording the values of the log-likelihood (respectively $S_{0}$ and $S_{1}$). We can then study the distribution of $TS=2~{}(S_{0}-S_{1})$. For the comparison between $m_{0}=f_{Band}$ and $m_{1}=f_{BHec}$ we find that the $TS$ is distributed as $ֿ12*χ*^2_1$,asshowninpanel1ofFigure~{}\ref{fig:ts_bhec}.Instead,forthecomparisonbetween$m_0=f_Band$and$m_1 =f_BB$wefindthattheTSisdistributedas$*χ*^2_2$,asexpectedfromWilks^{\prime}theorem(panel2).Wealsorepeatedthissimulationexerciseforeachtimeintervalusedinouranalysis,using100,000simulations,andverifiedthatineachcaseweareinthesameregime.Wecanthenusethesetwodistributionstodeterminethesignificanceoftheblackbodyandthecutoff.\par$

Appendix C Systematic uncertainties

The effective area correction we used during the fit, by introducing a multiplicative constant free to vary independently for each detector except one (see text), can help neutralizing effects due to systematic errors in the inter-calibration between the different detectors, introduced for example by underestimating or overestimating the effective area of a detector. It cannot, however, account for distortion in the spectrum introduced by errors in the energy-by-energy or channel-by-channel measurement of the instrument response model. In this section we study the impact of such uncertainties on our results.

We first briefly describe LLE data and the procedure we used to estimate the amount of systematic uncertainty in the relative response. We then proceed to study how such uncertainty, coupled with the uncertainty on the Fermi/GBM response, can modify our results.

C.1 LLE class

The LAT Low energy (LLE) technique is an analysis method designed to study bright transient phenomena, such as GRBs and solar flares, in the 30 MeV–1 GeV energy range. The LAT collaboration developed this analysis using a different approach than the one used in the standard photon analysis, which is based on sophisticated classification procedures (a detailed description of the standard analysis can be found in Atwood et al., 2009; Ackermann et al., 2012). The idea behind LLE is to maximize the effective area below $\sim$ 1 GeV by relaxing the standard analysis requirement on background rejection. The basic LLE selection is based on a few simple requirements on the event topology in the three sub-detectors of the LAT namely: a tracker/converter (TKR) composed of 18 x–y silicon strip detector planes interleaved with tungsten foils; an 8.6 radiation length imaging calorimeter (CAL) made with CsI(Tl) scintillation crystals; and an Anti-coincidence Detector (ACD) composed of 89 plastic scintillator tiles that surrounds the TKR and serves to reject the cosmic-ray background.

First of all, an event passing the LLE selection must have at least one reconstructed track in the TKR and therefore an estimate of the direction of the incoming photon. Secondly, we require that the reconstructed energy of the event be nonzero. The trigger and data acquisition system of the LAT is programmed to select the most likely gamma-ray candidate events to telemeter to the ground. The on-board trigger collects information from all three subsystems and, if certain conditions are satisfied, the entire LAT is read out and the event is sent to the ground. We use the information provided by the on-board trigger in LLE to efficiently select events which are gamma-ray like. In order to reduce the number of photons originating from the Earth Limb in our LLE sample we also include a cut on the reconstructed event zenith angle (i.e. angle $<$100*∘*). Finally we explicitly include in the selection a cut on the region of interest, i.e. the position in the sky of the transient source we are observing. In other words, the localization of the source is embedded in the event selection and therefore for a given analysis the LLE data are tailored to a particular location in the sky.

The response of the detector for the LLE class is encoded in a response matrix, which is generated using a dedicated Monte Carlo simulation for each GRB, and is saved in the standard HEASARC RMF File Format666Described here: http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#Sec:RMF-format.. LLE data and the relative response are made available for any transient signal (GRB or Solar Flare) detected with a significance above 4$\sigma$ through the HEASARC web site777FERMILLE, at http://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermille.html.

C.1.1 Validation and systematic uncertainties

Discrepancies between the actual response of the LAT and the response matrix derived from simulations can cause systematic errors in spectral fitting. We investigated the systematic uncertainties tied to the LLE selection by following the procedure described in Abdo et al. (2009c). In particular, we compared Monte Carlo with flight data, using the Vela pulsar (PSR J0835–4510) as a calibration source. The pulsed nature of the gamma-ray emission from this source (Abdo et al., 2010) gives us an independent control on the residual charged particle background. In fact, off-pulse gamma-ray emission is almost entirely absent, and a sample of “pure photons” can be simply extracted from the on-pulse region, after the off-pulse background is subtracted. Considering all time intervals during which the Vela pulsar was observed at an incidence angle $\theta<80^{\circ}$, we estimate the discrepancy between the efficiency of the LLE selection criteria in the LAT data and in Monte Carlo to be $\sim$17% below 100 MeV, decreasing to $\sim$8% at higher energies, with an average value $\sim$9% (note that this average is weighted by the Vela spectrum).

Additionally, we performed a spectral analysis of the Vela pulsar, comparing LLE results with the standard Fermi likelihood analysis. The $>$100 MeV flux obtained from the LLE analysis is 16% less than the flux reported by Abdo et al. (2010). This discrepancy can be attributed to the fact that the selection criteria between LLE and the standard LAT likelihood analysis are rather different with the former being much looser.

Finally, we also studied the energy resolution using large samples of simulated events with the Fermi/LAT full simulator. No significant bias was found, and the energy resolution for LLE is estimated to be $\sim$ 40% at 30 MeV, $\sim$30% at 100 MeV and $<$ 15% for energies greater than 100 MeV.

C.2 Fermi/GBM detectors

The systematic uncertainties on the responses of the Fermi/GBM detectors have been studied before launch, giving a calibration uncertainty of $\sim$ 10% for NaI detectors and $\sim$ 20% for BGO detectors (Bissaldi et al., 2009). In-flight calibration efforts have been limited so far, but preliminary results show a systematic uncertainty of around 15% in overall flux measurements, obtained comparing Fermi/GBM with other instruments (von Kienlin et al., 2009).

C.3 Effect of systematic errors on the measurement of the cutoff

In this section we describe a method we have used to estimate the impact that systematic uncertainties in the response of the instruments have on the detection of the cutoff.

For all instruments used in this work, the response is encoded in a response matrix. A response matrix888As encoded in OGIP RSP files is a bidimensional histogram having true energy on the x-axis and detected energy (or channels) on the y-axis, and a value proportional to the probability for a given photon of true energy $E_{true}$ to be detected at energy $E_{obs}$. Such responses are generated using simulations of the detectors, which must account for the geometry of the observation, the physical characteristics of the detector, and the physics of all the processes involved. All of these have uncertainties that are very difficult to study and to model, contributing to errors in the final response matrix. However, we have presented in the previous sections calibration studies that have been performed to estimate the overall systematic uncertainties of the effective area at different energies.

The approach we decided to use is to deliberately distort the response matrix in many different ways while keeping the difference in effective area between the original and the distorted matrix within the systematic uncertainties. In particular, this is our procedure: i) we consider the original response, with M true-energy channels and N observed-energy channels, and we generate a M x N matrix of uncorrelated noise ii) we smooth the noise matrix with a “diagonal kernel” that introduces correlation among the elements in the matrix along the diagonal direction iii) we smooth again the noise matrix with a Gaussian kernel, which removes unphysical jumps between neighboring elements iv) we re-normalize the noise matrix to have only elements between $1-f_{sys}$ and $1+f_{sys}$, where $f_{sys}$ is the fractional systematic uncertainty for the detector under examination. We used respectively 0.1 for NaI detectors, 0.2 for the BGO detector and 0.15 for LLE. We then multiply it by the original matrix to give a distorted matrix with at most $f_{sys}$ fractional variations v) we renormalize the distorted matrix to have the same total effective area as the original one, since the bias in the total area is already encapsulated by the effective area correction factor used in the fit.

Using this procedure we generated 500,000 distorted matrices for the data. For each realization we performed a fit, but using the distorted matrices in place of the original ones. We fit both $f_{Band}$ and $f_{BHec}$, and record the value of the best fit parameters as well as the value for TS. Among all intervals for both GRBs, interval 3 for GRB 100724B turned out to be the interval most affected by changes in the response matrices, and it consequently has the widest TS distribution, which is shown in Figure 17. We show the distribution with a logarithmic y-axis, which emphasizes the extremes of the distribution, but it must be noted that the vast majority of realizations change the TS value by only a few units. Among all intervals for both bursts, there are very rare instances where TS changes by up to 15 units, corresponding to very specific cases where changes in effective area of all the instruments conspire to change more drastically the significance of the cutoff. We can then assume a reduction of 15 units in TS for all time intervals as the most pessimistic case. It still would translate in 4 intervals above $5~{}\sigma$ for GRB 100724B, and 5 intervals for GRB 150509A. However, our simulations show that this is extremely unlikely: out of 500,000 simulations per interval, only a handful resulted in such a big change. Much more likely, systematic uncertainties change the value of TS by just few units, which does not affect almost at all our results. We can conclude that the detection of the cutoff is not likely affected by systematic uncertainties in the response of our instruments.

Appendix D Comparison with AGILE observations of GRB 100724B

GRB 100724B was detected and studied by AGILE as well (Del Monte et al., 2011). However, their results are different from what we observed with the LAT. For example, they observe a signal between the trigger time and 90 s with photons up to 3.5 GeV. They report a photon fluence in the 22 MeV - 3.5 GeV energy band of $0.25\pm 0.05$ ph. cm*-2* s and an energy fluence in the same band of $(4.7\pm 0.9\times 10^{-5})$ erg cm*-2*. They also measure a photon index of $-2.04^{+0.31}_{-0.14}$. These values correspond to a source much brighter and harder than what we see. To demonstrate this, we have performed a simulation of a source with the flux measured by AGILE. We used the tool gtobssim, part of the public Fermi Science Tools, which takes into account all the aspects of the LAT as well as the real pointing history of the satellite. We then compared it with what we observe. In Figure 18 we show a counts map of the simulation (left panel) and what we see (right panel), both using events above 300 MeV where the Earth Limb contamination in the data is small. For the sake of this comparison we did not introduce any cut on the zenith angle. It is apparent that what we observe is incompatible with what has been measured by AGILE. According to the simulation we should have detected $N_{pred}=60$ photons from the source above 300 MeV. We instead observe $N_{obs}=4$ photons. Even assuming that they all come from the source, the Poisson probability of observing $N_{obs}$ when we expect $N_{pred}$ is $5\times 10^{-21}$, i.e., essentially zero.

Given the extremely soft spectrum that we measure, we note that a possible culprit for such discrepancy could be the energy dispersion. For example, let us consider the photon that AGILE observed during the prompt emission with a reconstructed energy at 3.5 GeV. The LAT has an effective area above 1 GeV which is several times the one of AGILE/GRID, thus it should have detected several photons above 1 GeV in the same time interval, while we detect none. The energy resolution in AGILE/GRID is such that photons with energies below 100 MeV have a probability of being reconstructed well above 1 GeV of few percent (see upper left panel in Figure 3 in Chen et al., 2013). Given the brightness and softness of this GRB, there are hundreds of photons below 100 MeV in the LAT, as clear from the LLE light curve in Figure 1. AGILE/GRID observed 57 events. Thus it is entirely possible that one photon with a true energy below 100 MeV has a reconstructed energy above 1 GeV in AGILE/GRID. In Del Monte et al. (2011) it is not clear if and how the energy dispersion has been accounted for. However, ignoring the energy dispersion altogether or inaccuracies in its treatment can have a huge impact on the analysis of this burst. We note that in our analysis of the prompt emission, the energy dispersion in all instruments is accounted for.

Appendix E Spectral models

In this appendix we report the expressions for the spectral models used in this paper. All formulae are for the differential photon flux in photons cm*-2* s*-1* keV*-1*.

E.1 Power law with exponential cutoff

[TABLE]

where $\alpha$ is the photon index, $E_{c}$ is the cutoff energy, $K$ is the differential flux at $E_{piv}$, and $E_{piv}$ is the pivot energy, which we keep fixed to 1.

E.2 Band model

This is the model from Band et al. (1993), which consists of a low-energy power law and a high-energy power law joined by an exponential function:

[TABLE]

where $\alpha$ is the low-energy photon index, $\beta$ is the high-energy photon index, $E_{0}$ is the break energy and $K$ the normalization constant. It is easy to show that the peak energy in the Spectral Energy Distribution $E^{2}f_{Band}(E)$ is $E_{p}=(2+\alpha)E_{0}$.

E.3 Band model with exponential cutoff

A Band model multiplied by an exponential function:

[TABLE]

where $E_{c}$ is the cutoff energy.

E.4 Band model plus black body

A Band model plus a Planck function:

[TABLE]

E.5 Band model with power law break

A Band model where the second branch is a broken power law instead of a power law:

[TABLE]

E.6 Band model with smooth transition to a power law

A Band model where the high-energy power law changes index smoothly:

[TABLE]

with $\Delta\beta$ fixed to:

[TABLE]

as expected from theoretical considerations (Granot et al., 2008).