MAGIC Observations of the Nearby Short Gamma-Ray Burst GRB 160821B^*

GRB 160821B triggered the Swift-BAT detector at 22:29:13 UT on 2016 August 21 (hereafter T₀; Siegel et al. 2016). With the reported burst duration T₉₀ = 0.48 s, the event is classified as a short GRB. The spectrum of the prompt emission in the keV–MeV range is described by a power law with an exponential high-energy cutoff, with photon index 0.11 ± 0.88, peak energy E_p = (46.3 ± 6.4) keV, and fluence S(15–150 keV) = (1.0 ± 0.1) × 10⁻⁷erg cm⁻² (Palmer et al. 2016). The prompt emission was also detected by Fermi-GBM at the same trigger time as Swift-BAT, with T₉₀(50–300 keV) ∼ 1 s (Stanbro & Meegan 2016; refined later to 1.088 ± 0.977 s; von Kienlin et al. 2020). The spectrum is fit with a cutoff power-law function ⁴⁵ with E_p = 92 ± 28 keV and fluence S(10–1000 keV) = (1.7 ± 0.2) × 10⁻⁷ erg cm⁻². A host galaxy was identified (Xu et al. 2016) with spectroscopic redshift z = 0.162 (Levan et al. 2016; Lamb et al. 2019; Troja et al. 2019), making this one of the nearest short GRBs to date. With this redshift, the isotropic energy is estimated to be E_iso ∼ 1.2 × 10⁴⁹ erg, which is toward the low end of the known distribution for short GRBs, but not unusual (Berger 2014).

At the time of the GBM trigger, the burst was near the border of the standard FoV of Fermi-LAT (<60°). No emission was detected by LAT in the energy range 0.3–3 GeV (Section 3.1 for more details).

Follow-up observations were performed in the radio, optical, and X-ray bands (see corresponding light curves in Figure 1). Swift X-Ray Telescope (XRT; Burrows et al. 2005) started observations 57 s after the BAT trigger. The X-ray light curve, retrieved from the public on-line repository (Evans et al. 2009), reveals complex behavior, with an initial plateau followed by a steep decay. After ∼ 10³ s, a more commonly observed type of decay is seen, with ∝ t^−0.8, where t is time since T₀ (Troja et al. 2019).

Zoom In Zoom Out Reset image size

The optical afterglow was first reported by the Nordic Optical Telescope (Xu et al. 2016), and confirmed by the William Herschel Telescope (Levan et al. 2016), Gran Telescopio Canarias (Jeong et al. 2016), and the Hubble Space Telescope (Troja et al. 2016). Observations at different epochs confirmed that the optical source was fading, with a reported magnitude r = 22.6 ± 0.1 mag at 0.95 hr after T₀. The GRB is located in the outskirts of the host spiral galaxy, at ∼15 kpc projected distance from its center (Lamb et al. 2019; Troja et al. 2019). In the radio band, VLA (6 GHz) detected a fading source consistent with an afterglow (Fong et al. 2016).

MAGIC is a system composed of two imaging air Cherenkov telescopes, both with a mirror diameter of 17 m. The system is located at 2200 m above sea level, at the Roque de Los Muchachos Observatory in La Palma, Canary Islands, Spain. The integral sensitivity of the system is 0.66% of the Crab Nebula flux above 220 GeV with a 50 hr observation (Aleksić et al. 2016).

MAGIC started observing GRB 160821B at the Swift-BAT position (R.A.: +18^h39^m57^s; decl.: +62^d23^m34^s) on 2016 August 21, 22:29:37 UT, 24 s after T₀. The observation started from a zenith angle of 34°, and continued until 4 hr after T₀ (August 22, 2:29 UTC), reaching a zenith angle of 55°. The level of the night sky background (NSB) light was relatively high, due to the presence of the Moon. The NSB quickly increased during the observations as the Moon rose, ranging from two to eight times the level during dark nights (Ahnen et al. 2017). The first ∼1.7 hr of the data (until ∼0:10 UTC) were strongly affected by clouds, while the remaining ∼2.2 hr were taken under better weather conditions. The atmospheric transmission was measured using an LIDAR facility installed at the MAGIC site (Fruck & Gaug 2015). In the first part of the observation, the transmission at a height of 9 km fluctuated between 40% and 70% (average ∼60%) relative to that in good weather conditions, which prevented the use of standard analysis procedures. On the other hand, this was above 85% during the latter part, where the data quality was good enough for standard analysis with corrections for the transmission. Since the data for the first hours are potentially crucial for clarifying the physics of GRBs, we applied a dedicated analysis to recover this data, as described below.

The data analysis was performed entirely with the MAGIC standard analysis package called MARS (Zanin et al. 2013). The telescope performance, such as the sensitivity, energy threshold, and systematic errors for observations performed under nominal conditions can be found in Aleksić et al. (2016). In this analysis, however, we used dedicated software configurations optimized for high NSB levels. The description and performance study of this method can be found in Ahnen et al. (2017). After calibration of the data, we applied a more stringent image cleaning procedure compared to standard ones to remove a larger amount of spurious signals. As a consequence, the energy threshold is increased. We used data from known, bright gamma-ray sources (Crab Nebula, Mrk 421) observed under similar conditions (NSB and atmospheric transmission) to optimize the analysis cuts, in accordance with expected changes in the shower image parameters. The best sensitivity was obtained with a cut corresponding to an energy of ∼0.8 TeV. Thus, we used this threshold in order to maximize the sensitivity and search for possible signals.

The data were selected from a region of interest (RoI) of 12° centered on the best GRB position provided by Swift-XRT. We analyzed the data up to 10,000 s after T₀, considering a zenith angle cut of 100° to avoid Earth limb photons. A nonstandard FoV limit of ∼70° (in contrast to the normal 60°) was chosen in order to account for all analyzable data, including the entire RoI of 12° at the earliest times, when GRB 160821B was at a boresight angle of ∼61°. The source went outside the Fermi-LAT FoV at T₀ + 2315 s. During this first interval, the source was mainly around the border of the FoV of the instrument due to a previous automatic re-pointing request in the direction of GRB 160821A (Longo et al. 2016a, 2016b). It later reentered the FoV from T₀ + 5285 s up to T₀ + 8050 s, following the standard survey mode, moving in the LAT FoV from 70° down to 10° with respect to the boresight.

No hints of a detection were registered, and upper limits were derived adopting a Bayesian approach. Pass8 Source data (Atwood et al. 2013) in the energy range 0.3–3 GeV were selected and analyzed with the unbinned likelihood method, similar to that for the Second LAT GRB Catalog (Ajello et al. 2019), using the FermiTools package version 1.2.1, ⁴⁶ and the corresponding instrument response functions P8R3_SOURCE_V2. The GRB was modeled as a pointlike source, having a fixed power-law spectrum with photon index −2. The spectral parameters of other sources were kept fixed to those derived from the 4FGL catalog (Abdollahi et al. 2020). For the background modeling, appropriate isotropic extragalactic and Galactic models iso_P8R3_SOURCE_V2_v1.txt, gll_iem_v07.fits ⁴⁷ were employed. Two main time intervals were considered during which the source was in the FoV, from T₀ to T₀ + 2315 s and from T₀ + 5285 s to T₀ + 8050 s. The resulting upper limits in the 0.3–3 GeV energy band are, respectively, 1.7 × 10⁻⁶ cm⁻² s⁻¹ and 9.0 × 10⁻⁷ cm⁻² s⁻¹.

The excess significance at the GRB position observed by MAGIC is 3.1σ (pre-trial). It is derived using the prescription of Equation (17) of Li & Ma (1983), based on the distribution of the squared angular distance (θ²) between the reconstructed source position and the nominal source position (Figure 2). We tested two analysis cuts used for MAGIC data, the cut described above and an alternative cut that is optimized for low-energy events, so that a trial factor of two should be considered, leading to a corrected significance value of 2.9σ.

Zoom In Zoom Out Reset image size

We also computed a significance sky map of the observation (Figure 3) and found a spot with high significance (4.7σ pre-trial) 0.05 deg away from the GRB position. The post-trial significance of seeing such a hot spot at any place in the sky map is 3.0σ (1232 trials). Since this hot spot is close to the GRB position, we discuss whether it can be a possible signal from the GRB that appears displaced from its actual position.

Zoom In Zoom Out Reset image size

The systematic error in the telescope pointing is typically <002 and maximally ∼003 even with strong wind gusts. Thus the offset of 005 cannot be attributed to the telescope pointing alone. We also checked in the 4FGL catalog (Abdollahi et al. 2020) that there are no previously known GeV gamma-ray sources within 1 deg around the spot that could be potential TeV emitters.

We considered possible shifts of the reconstructed source position for a weak source embedded in a background that is fluctuating at a comparable level. We performed a Monte Carlo study simulating the sky maps, and found that the centroid position of the hot spot can be spread over a larger area than that of the actual signal. The hot spot position is distributed as a two-dimensional Gaussian with a width 2.6 times larger than that of the signal. The probability of the reconstructed position of such weak sources falling outside the original 1σ contour of the point-spread function (PSF; 0.045 deg in radius) is 24%. Therefore we conclude that the 005 offset seen in the sky map is well explained by statistical fluctuations in the case of weak signals, and that the significance of 3.1σ (pre-trial) conservatively computed at the Swift-XRT position can be regarded as evidence of a signal from the GRB.

We note that in addition to the trial factor discussed above, follow-up observations of other GRBs in the MAGIC GRB program may be considered as further trials. Among the 69 GRBs followed up by MAGIC in stereoscopic mode since 2009 (Carosi et al. 2015; MAGIC Collaboration et al. 2021, in preparation), the only short GRBs other than GRB 160821B observed under acceptable conditions were 140930B, 160927A, and 180715A, all with delays longer than 5000 s, and none with measured redshifts. Properly accounting for such observations as trials is difficult and not discussed in this paper, as they are subject to hidden observational and analysis biases, implying unequal trial factors.

In order to estimate flux values, we divided the data into two sets according to the weather conditions during the observations. The first 1.7 hr are characterized by low atmospheric transmission (average ∼60%), while the remaining 2.2 hr had good weather conditions. The first 1.7 hr are further subdivided into two time bins, to better represent the results on a logarithmic timescale. The resulting bins in time since T₀ are 24–1216 s, 1258–6098 s, and 6134–14,130 s. The flux is estimated by integrating the signal above 0.5 TeV, the peak energy of the reconstructed gamma-rays when assuming a power-law spectrum with photon index −2, convolved with the effective area. Because of the low significance in the first two time bins, we calculated 95% confidence-level flux upper limits using the method described in Rolke et al. (2005), obtaining 1.1 × 10⁻¹¹ cm⁻² s⁻¹ and 5.4 × 10⁻¹² cm⁻² s⁻¹, respectively. For the third time bin, we can similarly derive a flux upper limit of 3.0 × 10⁻¹² cm⁻² s⁻¹. On the other hand, despite the limited significance, we can also derive the flux for the last time bin, assuming that the excess is a real signal, which gives 9.9 ± 4.8 × 10⁻¹³ cm⁻² s⁻¹.

In order to check for the possibility of an unknown, unrelated gamma-ray source at the GRB position, we carried out an additional observation about a year after the GRB (2017 September 11–14, T₀ + 3.3 × 10⁷ s) and obtained 7.6 hr of good-quality data. The result is a flux upper limit of 4.4 × 10⁻¹³ cm⁻² s⁻¹ (>0.5 TeV, 95% C.L.), which is about half of the value discussed above for the putative signal. If a steady source was present at the position, an observation of 7.6 hr (instead of 2.2 hr) should result in a flux measurement with a smaller error, 9.9 ± 2.6 × 10⁻¹³ cm⁻² s⁻¹. The confidence belts of the flux inferred earlier and the flux upper limit derived later marginally overlap at the 2σ level on both sides, so the hypothesis of a steady source is disfavored, although it does not exclude the possibility of a variable source that is unrelated to the GRB.

Because of the low significance, an unfolded spectral energy distribution could not be derived, even for the third time bin with data obtained during good weather. The error box shown in the right panel of Figure 4 indicates only the reconstructed flux for this bin, derived from the photon flux by integrating over the energy range 0.5–5 TeV and assuming a power-law spectrum with photon index −2 (horizontal edges of the box). The height of the box corresponds only to statistical errors for the photon flux, and does not account for systematic errors related to the assumed spectral index.

Zoom In Zoom Out Reset image size

Several distinct components can contribute to the radio to X-ray emission of short GRBs after their prompt emission. The main component is synchrotron radiation from electrons accelerated in external forward shocks, triggered by interactions between the relativistic jet and the ambient medium (hereafter, simply "afterglow" radiation). In some cases, another component can arise from a reverse shock propagating into the jet ejecta. Two additional components are unique to short GRBs. Often seen in X-rays is "extended emission," where a relatively shallow temporal decay during the first tens to hundreds of seconds is followed by a much steeper decay, widely thought to be related to long-lasting activity of the central engine (either a magnetar or a black hole resulting from an NS merger; Norris & Bonnell 2006; Lü et al. 2015). Finally, optical-infrared kilonova emission can occur on timescales of days, powered by freshly synthesized r-process elements ejected in NS mergers (Metzger 2019).

All four of the aforementioned components are actually observed in GRB 160821B. Hereafter, our modeling focuses on the afterglow component from the external forward shock. Thus, we only consider the X-ray data at t > 10³ s, excluding the extended emission that can be clearly seen at earlier times in Figure 1 (see also Zhang et al. 2018). The kilonova emission has been inferred to dominate the optical/nIR band from 1 day to 4 days after T₀ (Kasliwal et al. 2017; Jin et al. 2018; Lamb et al. 2019; Troja et al. 2019).

The broadband light curves are shown in Figure 4 (left panel). We adopt the X-ray light curve from Troja et al. (2019) and model the broadband emission as synchrotron emission from the external forward shock, considering the simplest case of impulsive energy injection. The modeling is performed with a numerical code that self-consistently solves the evolution of the electron distribution, accounting for continuous electron injection with a power-law energy distribution (${dN}/d\gamma \propto {\gamma }^{-p}$), synchrotron, synchrotron-self-Compton (SSC) and adiabatic losses, synchrotron self-absorption, and γγ pair production (for a description of the code, see MAGIC Collaboration et al. 2019b and references therein).

The broadband SED at t ∼ T₀ + 3 hr is shown in Figure 4 (right panel). The consistency between the X-ray and optical spectral indices (F_ν ∝ ν^−0.8) suggests that the X-ray and optical bands are located between the characteristic synchrotron frequency ν_m and the cooling frequency ν_c. The radio data at 6 and 10 GHz together with optical and X-ray data constrain ν_m to be located between the radio and optical bands. The radio emission from the forward shock is then expected to increase with time (see dashed green curve in the left panel of Figure 4), implying that the observed radio emission at early times is dominated by another component, most likely from the reverse shock (Lamb et al. 2019; Troja et al. 2019; Lamb 2020). To be consistent with the radio upper limits at later times, ν_m must cross the radio band. All together, these observations constrain its value to be ν_m ≳ 4 × 10¹² Hz at t ∼ 10⁴ s and ${F}_{{\nu }_{{\rm{m}}}}^{\mathrm{syn}}\sim 0.03$ mJy. The model parameter space is further constrained by the requirement ν_c > ν_X up to at least 4 days (from the observed lack of a clear temporal break in X-rays). Order-of-magnitude estimates for the model parameters can be inferred by solving the equations

$$\begin{eqnarray*}&&\begin{array}{l}{\nu }_{{\rm{m}}}{(t\sim {10}^{4}{\rm{s}})\sim 2\times {10}^{12}\mathrm{Hz}{\epsilon }_{{\rm{e}},-1}^{2}(p-2)}^{2}/{(p-1)}^{2}{\epsilon }_{{\rm{B}},-4}^{1/2}\,{E}_{{\rm{k}},50}^{1/2}\\ \qquad \qquad =\,4\times {10}^{12}\,\mathrm{Hz},\\ {F}_{{\nu }_{{\rm{m}}}}^{\mathrm{syn}}(t\sim {10}^{4}\,{\rm{s}})\sim 0.04\,\mathrm{mJy}\,{\epsilon }_{{\rm{B}},-4}^{1/2}\,{n}_{-1}^{1/2}\,{E}_{{\rm{k}},50}=0.03\,\mathrm{mJy},\mathrm{and}\\ {\nu }_{{\rm{c}}}^{\mathrm{syn}}(t\sim 1\,\mathrm{day})\sim 4\times {10}^{20}\,\mathrm{Hz}\,{\epsilon }_{{\rm{B}},-4}^{-3/2}\,{n}_{-1}^{-1}\,{E}_{{\rm{k}},50}^{-1/2}\gt 2.4\times {10}^{18}\,\mathrm{Hz}\end{array}\end{eqnarray*}$$

(see, e.g., Sari et al. 1998; Panaitescu & Kumar 2000; Granot & Sari 2002), where E_k is the initial, isotropic-equivalent kinetic energy, n is the density of the surrounding medium, _e and _B are the fraction of energy dissipated behind the shock in accelerated electrons and the magnetic field, respectively, and p is the power-law index of the injected electron energy distribution.

We find good agreement for values of the model parameters within the following ranges: Log (E_k/erg) = [50–51], Log (_e) =[ − 1; − 0.1], Log (_B) = [ − 5.5; − 0.8], Log (n/cm⁻³) =[ − 4.85; − 0.24], and p = [2.2; 2.35]. The inferred values are very similar to the values inferred by Troja et al. (2019).

There is degeneracy between the parameters, which can be understood as follows: since ${\nu }_{{\rm{m}}}\propto {\epsilon }_{{\rm{e}}}^{2}\sqrt{{\epsilon }_{{\rm{B}}}\,{E}_{{\rm{k}}}}$ and ${F}_{{\nu }_{{\rm{m}}}}\,\propto {E}_{{\rm{k}}}\sqrt{{\epsilon }_{{\rm{B}}}\,n}$ for a fixed value of _e, the other parameters must satisfy ${\epsilon }_{{\rm{B}}}\propto {E}_{{\rm{k}}}^{-1}$ and $n\propto {E}_{{\rm{k}}}^{-1}$. E_k < 10⁵⁰ erg would imply large values of _B and n, resulting in ν_c < ν_X.

The result of the modeling is compared with observations in Figure 4. The reverse shock and kilonova components (dotted–dashed and dotted curves in the left panel) are taken from Troja et al. (2019).

Note that in contrast to Troja et al. (2019) and our modeling here, Lamb et al. (2019) proposed a different, multi-zone interpretation for the afterglow, invoking emission from a narrow jet component, as well as a slower outflow component caused by energy injection from the central engine at late times. This different interpretation is mainly driven by a different analysis of the X-ray data, resulting in an X-ray light curve with evidence for a double peak.

Assuming that the TeV gamma-ray signal obtained from MAGIC observations of GRB 160821B is real, we discuss possible mechanisms for TeV emission in short GRBs and assess their viability in accounting for these observations.

4.2.1. Synchrotron-self-Compton Emission (SSC)

Considering the parameter space constrained from the observed synchrotron emission in the radio, optical, and X-ray bands, we estimate the associated SSC component, and compare the results with MAGIC observations. Given the wide range of values still allowed, the expected flux of the SSC emission can vary by a few orders of magnitude.

The energy range covered by MAGIC observations lies in the range ${\nu }_{{\rm{m}}}^{\mathrm{SSC}}\lt {\nu }_{\mathrm{MAGIC}}\lt {\nu }_{{\rm{c}}}^{\mathrm{SSC}}$ where the flux can be analytically estimated by

$$\begin{eqnarray}&&{F}^{\mathrm{SSC}}(\nu )=\tau \,{F}^{\mathrm{syn}}({\nu }_{{\rm{m}}}^{\mathrm{syn}}){\gamma }_{{\rm{m}}}^{2}\,{\left(\displaystyle \frac{\nu }{{\nu }_{{\rm{m}}}^{\mathrm{SSC}}}\right)}^{-(p-3)/2},\end{eqnarray} \tag{ 1 }$$

subject to corrections of the order of $\sim \mathrm{ln}(\nu /{\nu }_{{\rm{m}}}^{\mathrm{SSC}})$, where τ = σ_T R n/3 and ${\nu }_{{\rm{m}}}^{\mathrm{SSC}}\simeq {\gamma }_{{\rm{m}}}^{2}\,{\nu }_{{\rm{m}}}^{\mathrm{syn}}$ (Sari & Esin 2001).

The SED in the right panel of Figure 4 shows the inferred flux in the energy range 0.5–5 TeV under the assumption that the gamma-ray signal from GRB 160821B is real (thin red box). This TeV flux would imply a large amount of energy in the SSC component, with a Compton parameter Y > 1. A large TeV flux is obtained for large values of _e, since they allow for lower _B and higher n (see equations in Section 4.1). Adopting the following parameter values: Log (E_k/erg) = 51, Log (_e) = − 0.1, Log (_B) = − 5.5, Log (n/cm⁻³) = − 1.3, and p = 2.2, at the time of the last MAGIC observation, the external shock radius is R ∼ 2.4 × 10¹⁷ cm, the bulk Lorentz factor is Γ ∼ 16, and the characteristic Lorentz factor of the electrons is γ_m ∼ 3.8 × 10³. The expected SSC flux at 1 TeV is then F^SSC(1 TeV) ∼ 2 × 10⁻¹³ erg cm⁻² s⁻¹. The expected SSC spectrum accounting for extragalactic background light (EBL) attenuation (Domínguez et al. 2011) is shown as the red solid curve in Figure 4. Compared with the flux suggested from MAGIC observations at 0.1 day, the maximum flux expected in one-zone SSC models falls short by about an order of magnitude.

4.2.2. Proton Synchrotron Emission

4.2.3. Other Possibilities