Gamma-Ray Bursts Afterglow Physics and the VHE Domain

Abstract

Afterglow radiation in gamma-ray bursts (GRB), extending from the radio band to GeV energies, is produced as a result of the interaction between the relativistic jet and the ambient medium. Although in general the origin of the emission is robustly identified as synchrotron radiation from the shock-accelerated electrons, many aspects remain poorly constrained, such as the role of inverse Compton emission, the particle acceleration mechanism, the properties of the environment and of the GRB jet itself. The extension of the afterglow emission into the TeV band has been discussed and theorized for years, but has eluded for a long time the observations. Recently, the Cherenkov telescopes, MAGIC and H.E.S.S., have unequivocally proven that afterglow radiation is also produced above 100 GeV, up to at least a few TeV. The accessibility of the TeV spectral window will largely improve with the upcoming facility CTA (the Cherenkov Telescope Array). In this review article, we first revise the current model for afterglow emission in GRBs, its limitations and open issues. Then, we describe the recent detections of very high energy emission from GRBs and the origin of this radiation. Implications on the understanding of afterglow radiation and constraints on the physics of the involved processes will be deeply investigated, demonstrating how future observations, especially by the CTA Observatory, are expected to give a key contribution in improving our comprehension of such elusive sources.

1. Introduction

Gamma-ray bursts (GRBs) are observed as transient sources of radiation displaying a distinctive pattern that consists of two different phases. The first phase is dominated by emission in the keV-MeV energy range, lasting from fractions of a second to several minutes, and reaching isotropic equivalent peak luminosities in the range L∼${10}^{49}$–${10}^{53}$ erg s${}^{-1}$. The bimodal distribution of the prompt emission duration reveals that there are two classes of GRBs, called short and long depending on whether the prompt emission lasts shorter or longer than 2 s [1,2]. The second emission phase, called the afterglow, follows the prompt with a delay of tens of seconds, and is detected on a very wide range of frequencies, from $\gamma$-rays to the radio band. The afterglow flux decays smoothly as a power-law in time for weeks or months, and the typical frequency of the radiation moves in time from the X-ray to the radio band. Since 2019, the detection of a few long GRBs between 0.3 and 3 TeV on time-scales from tens of seconds to a few days proved for the first time that GRBs can also be sources of radiation in the TeV band, where they can convey a sizable fraction (20–50%) of the total energy emitted during the afterglow phase [3,4,5].

All this prompt/afterglow emission is identified with radiation produced as a result of the launch of an ultra-relativistic ($\Gamma$∼100–1000) jet from a newly born compact object. The ejecta first undergoes internal dissipation (through mechanisms such as shocks between different parts of the outflow [6] or magnetic reconnection episodes [7,8]). In a second moment, the ejecta undergoes external dissipation [9], triggered by interactions with the ambient medium (e.g., the interstellar medium or the wind of the progenitor’s star [10]). The two different dissipation processes occur at different typical distances from the central engine (R∼${10}^{13-14}$ cm and R∼${10}^{15-20}$ cm) and generate two well distinguished emission phases, identified as the prompt and afterglow emission, respectively.

For long GRBs, it is widely believed that the involved energetics and time-scales and the successful launch of a relativistic jet can find justification in the collapsar model [11,12]. In this model, the core of a massive star collapses into a black hole and the accretion from the surrounding disk powers the launch of two opposite, collimated (${\theta }_{jet}$∼5–${10}^{\circ }$) outflows. A similar scenario also applies to short GRBs, where the black hole originates from the merger of two neutron stars (as recently proven by the association of a short GRB with a gravitational wave signal [13]) or a neutron star and a black-hole. An alternative model [14,15,16,17] considers a millisecond magnetar (i.e., a rapidly rotating neutron star) as the progenitor of long GRBs (or at least a fraction of them). This model has the advantage of more naturally explaining the detections of late time activity (${10}^2$–${10}^3$ s after the prompt onset) in the form of X-ray flares and plateaus, observed in about one third of the population.

Beside the nature of the progenitor’s star, another quite pressing open issue in GRB physics concerns the composition of the jet itself, i.e., the nature of the dominant energy stored in the outflow, which can be either magnetic (in the form of Poynting flux [7,14]) or kinetic (i.e., bulk motion of the matter). This uncertainty reflects an uncertainty on the mechanism extracting energy from the jet (i.e., the process converting part of the jet energy into random energy of the particles), which is identified with internal shocks in the latter case, and magnetic reconnection events in the case of a Poynting flux dominated outflows [14]. While internal shocks in a matter-dominated jet have been considered the mainstream model for a long time, tensions between some model predictions and observations have moved the attention in the last decade to a family of models based on magnetic jets [18,19,20]. In particular, internal shocks are not an efficient mechanism [21,22], and this is in contrast with the evidence that only a relatively small fraction (10–50%) of energy is still in the blast-wave during the afterglow phase, meaning that most of it must have been dissipated and radiated away during the prompt. It must be noted, however, that the estimate of the energy content of the blast during the afterglow is indirect, and contingent upon a proper modeling of the afterglow emission [23]. Investigations that took advantage of GeV emission detected by LAT (the Large Area Telescope onboard the Fermi satellite), reached the conclusion that the blast energy is usually underestimated by studies relying on X-ray emission, and inferred a prompt emission efficiency between 1–10% [24], which is still consistent with internal shocks. The nature and efficiency of the dissipation mechanism in the prompt phase are still matters of intense debate. In any case, the radiation is expected to be produced by the accelerated electrons, which efficiently lose energy via synchrotron cooling [6,25]. Inconsistencies between the expected synchrotron spectrum and the observed spectral shape of the prompt emission [25,26] have also called into question the nature of the radiative process. Recent works have performed major advances towards the comprehension of the radiative mechanism responsible for the prompt emission, supporting the synchrotron interpretation [27,28,29,30,31].

The nature of the afterglow emission is much better understood, at least on its general grounds. The interaction between the jet and the external medium triggers the formation of a forward shock running into the external medium and a reverse shock running into the ejecta [32,33,34]. These shocks are responsible for the acceleration of particles and for the deceleration of the outflow, eventually down to non-relativistic velocities [35,36,37]. The observed radiation is the result of synchrotron radiation from electrons accelerated at the forward shock [38]. A contribution from the reverse shock may also be relevant, typically in the radio and optical band [34,39]. Shock formation and particle acceleration in ultra-relativistic shocks are still not completely understood. Very important progresses have been achieved in the last decade on the theoretical side (see [40] for a recent review), especially thanks to numerical particle-in-cell (PIC) simulations [41,42]. Ultra-relativistic shocks in a weakly magnetized medium are found to be efficient particle accelerators, with ${\xi }_e\sim 1\%$ of the electrons being accelerated into a power-law distribution with spectral index $p\sim 2.5$, carrying about ${ϵ}_e\sim 10\%$ of the shock-dissipated energy. A strong magnetic turbulence is self-generated by the accelerated particles counter-streaming in the upstream, ahead of the shock, at a level of magnetization ${ϵ}_B=$ 0.01–0.1 [42]. PIC simulations, however, are currently probing time-scales that are orders of magnitude smaller than the dynamical time-scale of the blast-wave. This implies that results from simulations can only be extrapolated to the relevant time-scales, introducing a certain degree of uncertainty and caution in using the results as inputs for the modeling of GRB afterglows. What is still poorly understood, even though dedicated simulations are starting to give important clues [43], is how the micro-turbulence generated in the shock vicinity evolves (decays) with time. This is particularly important for a proper interpretation of the observations, since it is likely that the particles produce synchrotron and synchrotron-self Compton (SSC) photons in a region of decayed micro-turbulence, and hence feel a magnetic field with ${ϵ}_B<<0.01$.

The afterglow emission, its spectral shape from radio to $\gamma$-rays, and its temporal evolution from seconds to months, contain a wealth of (convoluted) information on blast dynamics, particle acceleration, magnetic turbulence generation and decay, and external density in the progenitor’s surroundings (up to a parsec scale). Nevertheless, since all these ingredients are poorly constrained from theoretical grounds, they enter the afterglow physics as free, unconstrained model parameters. The large degeneracy among different parameters and the small number of observables as compared to model variables are limiting our possibility to extract valuable and robust information from the modeling of the observed afterglow radiation. To go beyond the state-of-the-art, additional efforts are necessary both on the observational and theoretical sides.

An interesting opportunity has recently opened on the observational side, thanks to the discovery that GRBs can be sources of TeV radiation associated with the afterglow phase [3,4]. The characterization of the TeV spectra and light curves offers new observables to further constrain the unknown physics of the afterglow emission. These observations are expected to impact on our current understanding of the environment where GRBs explode (and hence on the nature of their progenitors), of the physics of ultra-relativistic shocks, and of the properties of the jet (e.g., bulk Lorentz factor and energy content). Constraining the jet properties is mandatory for a correct estimate of the prompt mechanism efficiency and then for determining its nature. It is then evident how the opening of this completely new energy window in GRBs is expected to boost the studies in a field that has many connections both with the general understanding of the GRB phenomenon and with topics of general interest, such as star formation and evolution, the last stages of massive stars and their environments and plasma physics under extreme conditions.

Given the impressive amount of new information that VHE observations are going to bring to the field, it is important to revise what is the state-of-the-art, what are the main issues and how we can benefit from the few existing and the upcoming observations in the TeV domain. This review revisits the present understanding of afterglow radiation, the discovery of very-high energy (VHE, $>100$ GeV) emission from GRBs and future prospects for the detection of GRBs at VHE with the next generation of Cherenkov telescopes.

For recent and complete reviews on GRB’s phenomenology and theoretical interpretation before the TeV era see [44,45]. An overview focused on high-energy emission (0.1–100 GeV) observations and interpretation can be found in [46].

This review is organized as follows. Section 2 presents an overview of the afterglow external shock model, revisiting our common understanding and phenomenological description of (i) the dynamics of the blast-wave, (ii) shock formation, particle acceleration and self-generation of turbulent magnetic field in the shock proximity, and (iii) the main processes shaping the radiative output, on the whole electromagnetic spectrum, from radio to very-high energy $\gamma$-rays. In Section 3, we propose a discussion of the main open issues of the afterglow model, outlining which observations are at odds with model predictions, which observed features are missing in the basic scenario and what are the present limitations that prevent us from extracting valuable information from the modeling of multi-wavelength afterglow radiation. In Section 4, we describe the recent discovery that GRBs can be bright TeV emitters. Each GRB with a firm (or with a hint of) detection by MAGIC or H.E.S.S. is discussed in detail. We present multi-wavelength observations and review the proposed interpretations of the detected emission. In Section 5, we compare the general properties of the detected GRBs both among each other and with the general population. We discuss how the TeV emission can help to solve some of the most important issues of the afterglow model. Finally, in Section 6, we discuss the prospects for future studies of TeV emission from GRBs with the next generation of Cherenkov telescopes and their expected impact on GRB physics.

2. The Afterglow Model

Afterglow emission refers to all the broad-band radiation observed from a GRB on longer timescales (minutes to months) as compared to the initial prompt radiation detected in hard X-rays [38,47,48]. Its temporal evolution is usually well described with simple decaying power-laws, in contrast with the short-time (<seconds) variability that characterizes the prompt emission [49,50,51,52]. These major differences place the emission region of afterglow radiation at larger radii (>10${}^{15}$ cm), pinpointing its origin in the processes triggered by the interaction between the jet and the circumburst medium.

The expansion of the relativistic jet into the external medium is expected to drive two different shocks: the forward shock, running into the external medium, and the reverse shock, running into the jet. The shocked ejecta and the shocked external medium, separated by the contact discontinuity, are both sources of synchrotron radiation from the accelerated electrons [53]. Most of the detected radiation is interpreted as emission from ambient particles energized by the forward shock. Spectra and lightcurves are then shaped by the environment where the GRB explodes, which in turn is strictly connected to the nature of the progenitor. The other player that shapes the properties of afterglow radiation is particle acceleration at relativistic shocks, which is thought to proceed via diffusive shock acceleration, but for which the details of the underlying physics still remain poorly constrained. Moreover, the overall luminosity of the afterglow radiation depends on the energy content of the blast-wave. Such an amount is determined by how efficiently the prompt mechanism has dissipated and released part of the initial explosion energy. Following these considerations, it is evident how the study of afterglow radiation impacts on the general understanding of the GRB phenomenon: the progenitor and its environment, the nature and efficiency of the mechanisms responsible for prompt emission, the properties of the jet, and the micro-physics of relativistic shocks.

In this section, the physics involved in the afterglow scenario is presented, with a particular focus on the forward shock emission and on the radiative output expected at VHE. This section is organized as follows: we revisit the physics of the jet dynamical evolution in its interaction with the ambient medium (Section 2.1), the particle acceleration mechanism (Section 2.2) and the resulting radiative output and its spectral shape (Section 2.3).

2.1. Jet Dynamics

After the reverse shock has crossed the ejecta, the dynamics of the blast-wave enters a self-similar regime ([35], BM76 hereafter). In a thin shell approximation, the reverse shock crossing time corresponds to the time when the blast-wave starts decelerating. The deceleration of the jet, caused by the collision with the external medium, becomes significant at the radius $R_{dec}$, where the energy transferred to the mass m collected from the external medium (∼$m\left(R_{dec}\right)c^2{\Gamma }_0^2$) is comparable to the initial energy ($E_0=M_0{\Gamma }_0{c}^2$) carried by the jet. This deceleration radius is typically of the order of $R_{dec}$∼${10}^{15}$–${10}^{16}$ cm, depending on the density of the external medium, the ejecta mass $M_0$ and initial bulk Lorentz factor ${\Gamma }_0$. Before reaching this radius, the ejecta expands with constant velocity (coasting phase).

Most analytic estimates of the afterglow evolution with the purpose of modeling data are developed for the deceleration phase, where the self-similar BM76 solution for adiabatic blast-waves is adopted [38,47,54]. Since VHE emission can be detected at quite early times (a few tens of seconds), we are also interested in the description of the coasting phase and in a proper treatment of the transition between coasting and deceleration.

In the following, to derive the evolution of the bulk Lorentz factor, we adopt the approach proposed by [55]. This method allows us to describe the hydrodynamics of a relativistic blast-wave expanding into a medium with an arbitrary density profile $\rho \left(R\right)$ and composition (i.e., enriched by pairs), and the transition from the free expansion of the ejecta to the deceleration phase, taking into account the role of radiative and adiabatic losses. The internal structure is neglected (homogeneous shell approximation), and the Lorentz factor $\Gamma$ considered is the one of the fluids just behind the shock front. In the deceleration phase, the self-similar solutions derived in BM76 are recovered by this method, both for the adiabatic and the fully radiative cases, and for constant and wind-like density profiles of the external medium. The presented approach also allows us to introduce a time-varying radiative efficiency, either resulting from a change with time of ${ϵ}_{\mathrm{e}}$ or a change in the radiative efficiency of the electrons. Equations reported here are valid after the reverse shock has crossed the ejecta. Corrections to the hydrodynamics before the reverse-shock crossing time can be found in [55].

2.1.1. Equation Describing the Evolution of the Bulk Lorentz Factor

The aim is to derive an equation describing the change $d\Gamma$ of the bulk Lorentz factor of the fluid just behind the shock in response to the collision with a mass $dm\left(R\right)=4\pi R^2\rho \left(R\right)dR$ encountered when the shock front moves from a distance R to $R+dR$ and with $\rho$ being the mass density. The change in $\Gamma$ is determined by dissipation of the bulk kinetic energy, the conversion of internal energy back into bulk motion, and injection of energy into the blast-wave. The latter is sometimes invoked to explain plateau phases in the X-ray early afterglow or to explain flux rebrightenings [56,57,58]. The following treatment neglects energy injection, which, however, can be easily incorporated in this kind of approach.

To write the equation for energy conservation, from which $d\Gamma /dR$ can be derived, we first need to recall how the energy density transforms under Lorentz transformations. In the following, we denote quantities measured in the frame comoving with the shocked fluid (comoving frame), with a prime, to distinguish them from quantities measured in the frame of the progenitor star (rest frame, without a prime).

The energy density in the comoving frame is $u^{\prime }=u_{\mathrm{int}}^{\prime }+{\rho }^{\prime }{c}^2$, where $u_{\mathrm{int}}^{\prime }$ is the comoving internal energy and ${\rho }^{\prime }$ is the comoving mass density. Applying Lorentz transformations, $u=(u^{\prime }+p^{\prime }){\Gamma }^2-p^{\prime }$, where $p^{\prime }$ is the pressure and is related to the internal energy density by the equation of state, and $p^{\prime }=(\widehat{\gamma }-1)u_{\mathrm{int}}^{\prime }=(\widehat{\gamma }-1)(u^{\prime }-{\rho }^{\prime }{c}^2)$, where $\widehat{\gamma }$ is the adiabatic index of the shocked plasma. The energy density is then given by: $u=u_{\mathrm{int}}^{\prime }(\widehat{\gamma }{\Gamma }^2-\widehat{\gamma }+1)+{\rho }^{\prime }{c}^2{\Gamma }^2$, which shows how the internal energy and rest mass density transform. The total energy in the progenitor frame will be $E=uV=u{V}^{\prime }/\Gamma$, where V is the shell volume in the progenitor frame, and can be expressed as:

(1)$E=\Gamma M{c}^2+{\Gamma }_{\mathrm{eff}}{E}_{\mathrm{int}}^{\prime },$

where:

(2)${\Gamma }_{\mathit{eff}}\equiv \frac{\widehat{\gamma }{\Gamma }^2-\widehat{\gamma }+1}{\Gamma },$

which properly describes the Lorentz transformation of the internal energy. Here, $M=M_0+m={\rho }^{\prime }{V}^{\prime }$ is the sum of the ejecta mass $M_0=E_0/{\Gamma }_0{c}^2$ and of the swept-up mass $m\left(R\right)$, and $E_{int}^{\prime }=(u^{\prime }-{\rho }^{\prime }{c}^2)V^{\prime }$ is the comoving internal energy. The adiabatic index can be parameterized as $\widehat{\gamma }=(4+{\Gamma }^{-1})/3$ to obtain the expected limits $\widehat{\gamma }\simeq 4/3$ for $\Gamma \gg 1$ and $\widehat{\gamma }\simeq 5/3$ for $\Gamma \to 1$. The majority of analytical treatments use $\Gamma$ instead of ${\Gamma }_{\mathit{eff}}$, which implies an error up to a factor of $4/3$ in the ultra-relativistic limit [55].

The blast-wave energy E in Equation (1) can change due to (i), the rest mass energy $dm{c}^2$ collected from the medium, (ii) radiative losses $d{E}_{rad}={\Gamma }_{\mathit{eff}}d{E}_{rad}^{\prime }$ and (iii) injection of energy. Ignoring possible episodes of energy injections into the blast-wave, the equation of energy conservation in the progenitor frame is:

(3)$d\left[\Gamma (M_0+m)c^2+{\Gamma }_{\mathit{eff}}{E}_{int}^{\prime }\right]=dm{c}^2+{\Gamma }_{\mathit{eff}}d{E}_{rad}^{\prime }.$

The overall change in the comoving internal energy $d{E}_{int}^{\prime }$ results from the sum of three contributions:

(4)$d{E}_{int}^{\prime }=d{E}_{sh}^{\prime }+d{E}_{ad}^{\prime }+d{E}_{rad}^{\prime }.$

The first contribution, $d{E}_{sh}^{\prime }=(\Gamma -1)dm{c}^2$, is the random kinetic energy produced at the shock as a result of the interaction with an element $dm$ of circum-burst material: as pointed out by BM76, in the post-shock frame, the average kinetic energy per unit mass $d{E}_{sh}^{\prime }/dm$ is constant across the shock, and equal to $(\Gamma -1)c^2$. The second term in Equation (4), $d{E}_{ad}^{\prime }$, is the internal energy lost due to adiabatic expansion, that leads to a conversion of random energy back to bulk kinetic energy. The third term, $d{E}_{rad}^{\prime }$, accounts for radiative losses.

From Equation (3), it follows that the variation of the Lorentz factor is:

(5)$\frac{d\Gamma }{dR}=-\frac{({\Gamma }_{\mathit{eff}}+1)(\Gamma -1)c^2\frac{dm}{dR}+{\Gamma }_{\mathit{eff}}\frac{d{E}_{ad}^{\prime }}{dR}}{(M_0+m)c^2+E_{int}^{\prime }\frac{d{\Gamma }_{\mathit{eff}}}{d\Gamma }},$

from which the evolution $\Gamma \left(R\right)$ of the bulk Lorentz factor of the fluid just behind the shock as a function of the shock front radius can be derived.

The term ${\Gamma }_{\mathrm{eff}}d{E}_{ad}^{\prime }/dR$, accounting for adiabatic losses, allows us to describe the re-acceleration of the fireball: this contribution, usually neglected, becomes important only when the density decreases faster than $\rho \propto R^{-3}$. To evaluate Equation (5), it is necessary to first specify $d{E}_{ad}^{\prime }$ and $E_{int}^{\prime }$.

2.1.2. Internal Energy and Adiabatic Losses

In specific cases, the adiabatic losses and the internal energy content can be expressed in an analytic form. The following treatment to estimate adiabatic losses and the internal energy content of the blast-wave assumes that, right behind the shock, the freshly shocked electrons instantaneously radiate a fraction ${ϵ}_{\mathrm{rad}}$ of their internal energy and then they cool only due to adiabatic losses [55]. By assuming that the accelerated electrons promptly radiate at the shock, and then they evolve adiabatically, one is implicitly considering either a fast cooling regime or quasi-adiabatic regime, in which case the radiative losses do not affect the shell dynamics.

Defining ${ϵ}_{\mathrm{e}}$ as the fraction of energy $d{E}_{sh}^{\prime }$ dissipated by the shock and gained by the leptons, the mean random Lorentz factor of post-shock leptons becomes (for a more detailed discussion see Section 2.2):

(6)${\gamma }_{acc,e}-1=(\Gamma -1){ϵ}_e/{\mu }_e.$

Here, ${\mu }_e={\rho }_e/\rho$ is the ratio between the mass density ${\rho }_e$ of shocked electrons and positrons (simply “electrons” from now on) and the total mass density of the shocked matter $\rho$. In the absence of electron-positron pairs ${\mu }_e=m_e/(m_e+m_p)\simeq m_e/m_p$.

Leptons then radiate a fraction ${ϵ}_{\mathrm{rad}}$ of their internal energy, i.e., the energy lost to radiation is $d{E}_{rad}^{\prime }=-{ϵ}_{rad}{ϵ}_{e}d{E}_{sh}^{\prime }=-ϵd{E}_{sh}^{\prime }$, with $ϵ\equiv$${ϵ}_{\mathrm{rad}}$${ϵ}_{\mathrm{e}}$ being the overall fraction of the shock-dissipated energy that goes into radiation. After radiating a fraction ${ϵ}_{\mathrm{rad}}$ of their internal energy, the mean random Lorentz factor of the freshly shocked electrons decreases down to:

(7)${\gamma }_{rad,e}-1=(1-{ϵ}_{rad})({\gamma }_{acc,e}-1)=(1-{ϵ}_{rad})(\Gamma -1)\frac{{ϵ}_e}{{\mu }_e}.$

The assumption of instantaneous radiative losses is verified in the fast cooling regime (${ϵ}_{rad}\sim 1$), which is required (but not sufficient) to have $ϵ\sim 1$ (i.e., a fully radiative blast-wave). In the opposite case ${ϵ}_{rad}\ll 1$, the evolution is nearly adiabatic ($ϵ\ll 1$), regardless of the value of ${ϵ}_{\mathrm{e}}$, and the details of the radiative cooling processes are likely to be unimportant for the shell dynamics. The case with intermediate values of ${ϵ}_{\mathrm{rad}}$ and $ϵ$ is harder to treat analytically, since the electrons shocked at radius R may continue to emit copiously at larger distances as well, affecting the blast-wave dynamics.

A similar treatment can be adopted for protons: if protons gain a fraction ${ϵ}_{\mathrm{p}}$ of the energy dissipated by the shock (with ${ϵ}_p=1-{ϵ}_e-{ϵ}_B$), their mean post-shock Lorentz factor will be:

(8)${\gamma }_{acc,p}-1=(\Gamma -1)\frac{{ϵ}_p}{{\mu }_p},$

where ${\mu }_p={\rho }_p/\rho$ is the ratio between the mass density of shocked protons ${\rho }_p$ and the total shocked mass density $\rho$. In the standard case, when pairs are absent, ${\mu }_p\simeq 1$. Since the proton radiative losses are negligible, the shocked protons will lose their energy only due to adiabatic cooling.

Adiabatic losses can be computed starting from $d{E}_{\mathrm{int}}^{\prime }=-p^{\prime }d{V}^{\prime }$, where $p^{\prime }$ is the pressure in the comoving frame. For N particles with Lorentz factor $\gamma$, the internal energy density is:

(9)$u_{\mathrm{int}}^{\prime }=\frac{N(\gamma -1)m{c}^2}{V^{\prime }}.$

The radial change of the Lorentz factor, as a result of expansion losses, is:

(10)${\left(\frac{d(\gamma -1)}{dr}\right)}_{ad}=-(\widehat{\gamma }-1)(\gamma -1)\frac{dln{V}^{\prime }}{dr}.$

To estimate the adiabatic losses, let us assume that the shell comoving volume scales as $V^{\prime }\propto R^3/\Gamma$, corresponding to a shell thickness in the progenitor frame $\sim R/{\Gamma }^2$. This scaling is correct for both relativistic and non-relativistic shocks in the decelerating phase (BM76). For re-accelerating relativistic shocks, Shapiro [59] demonstrated that the thickness of the region containing most of the blast-wave energy is still $\sim R/{\Gamma }^2$. For the sake of simplicity, changes in the comoving volume due to a time-varying adiabatic index or radiative efficiency are neglected. If the scaling $V^{\prime }\propto R^3/\Gamma$ is assumed, the equation can be further developed analytically, and reads:

(11)${\left(\frac{d(\gamma -1)}{dr}\right)}_{ad}=-(\widehat{\gamma }-1)(\gamma -1)\left(\frac{3}{r}-\frac{dln\Gamma }{dr}\right).$

The comoving Lorentz factor at radius R, for a particle injected with $\gamma \left(r\right)$ when the shock radius was r, will be

(12)${\gamma }_{\mathrm{ad}}(R,r)-1={\left(\frac{r}{R}\right)}^{3(\widehat{\gamma }-1)}{\left[\frac{\Gamma \left(R\right)}{\Gamma \left(r\right)}\right]}^{(\widehat{\gamma }-1)}(\gamma \left(r\right)-1).$

where $\gamma \left(r\right)$ is given by ${\gamma }_{rad,e}\left(r\right)$ (Equation (7)) for leptons, and by ${\gamma }_{acc,p}\left(r\right)$ (Equation (8)) for protons.

Considering the proton and lepton energy densities separately, the comoving internal energy at radius R will be:

(13)$E_{int}^{\prime }\left(R\right)=4\pi c^2{\int }_0^{R}dr{r}^2\left\{{\rho }_p\left(r\right)[{\gamma }_{ad,p}(R,r)-1]+{\rho }_e\left(r\right)[{\gamma }_{ad,e}(R,r)-1]\right\}.$

With the help of Equation (12), one can explicitly find $E_{int}^{\prime }\left(R\right)$ and insert it in Equation (5).

The other term needed in Equation (5) is $d{E}_{ad}^{\prime }/dR$. First, we have derived ${(d\gamma /dR)}_{ad}$ for a single particle. Now integrating over the total number of particles, again considering protons and leptons separately, one obtains:

(14)$\frac{d{E}_{ad}^{\prime }\left(R\right)}{dR}=-4\pi c^2(\widehat{\gamma }-1)\left(\frac{3}{R}-\frac{dlog\Gamma }{dR}\right){\int }_0^{R}dr{r}^2\left\{{\rho }_p\left(r\right)({\gamma }_{ad,p}-1)+{\rho }_e\left(r\right)({\gamma }_{ad,e}-1)\right\}.$

In Equations (13) and (14), it is assumed that only the swept-up matter is subject to adiabatic cooling, i.e., that the ejecta particles are cold.

As long as the shocked particles remain relativistic, the equations for the comoving internal energy and for the adiabatic expansion losses assume simpler forms:

(15)$\begin{array}{cc}\hfill E_{int}^{\prime }\left(R\right)& =4\pi c^2{\int }_0^{R}dr{r}^2\frac{r}{R}{\left[\frac{\Gamma \left(R\right)}{\Gamma \left(r\right)}\right]}^{1/3}\Gamma \left(r\right)\left\{\frac{{ϵ}_p}{{\mu }_p}{\rho }_p+(1-{ϵ}_{rad})\frac{{ϵ}_e}{{\mu }_e}{\rho }_e\right\}\hfill \\ & =4\pi c^2{\int }_0^{R}dr{r}^2\frac{r}{R}{\left[\frac{\Gamma \left(R\right)}{\Gamma \left(r\right)}\right]}^{1/3}\Gamma \left(r\right)\rho \left(r\right)\left({ϵ}_p+{ϵ}_e-ϵ\right)\hfill \\ & \frac{d{E}_{ad}^{\prime }\left(R\right)}{dR}=-E_{int}^{\prime }\left(R\right)\left(\frac{1}{R}-\frac{1}{3}\frac{dlog\Gamma }{dR}\right).\hfill \end{array}$

In the absence of significant magnetic field amplification, ${ϵ}_p+{ϵ}_e\simeq 1$ so that ${ϵ}_p+{ϵ}_e-ϵ\simeq 1-ϵ$, and the radiative processes of the blast-wave are entirely captured by the single efficiency parameter $ϵ$. In the fast cooling regime ${ϵ}_{rad}\sim 1$ and $ϵ\simeq {ϵ}_e$. In this case, the term ${ϵ}_p+{ϵ}_e-ϵ$ reduces to ${ϵ}_{\mathrm{p}}$, meaning that, regardless of the amount of energy gained by the electrons, in the fast cooling regime the adiabatic losses are dominated by the protons, since the electrons lose all their energy to radiation.

Evaluating these expressions for adiabatic blast-waves in a power-law density profile $\rho \propto R^{-s}$, one obtains:

(16)$E_{int}^{\prime }\left(R\right)=\frac{9-3s}{9-2s}\Gamma m{c}^2,\frac{d{E}_{ad}^{\prime }\left(R\right)}{dR}=-\frac{(9-s)(3-s)}{2(9-2s)}\frac{\Gamma m{c}^2}{R},$

where $\Gamma \propto R^{-(3-s)/2}$ as in the adiabatic BM76 solution has been used.

In the fully radiative regime $ϵ=1$, which implies $E_{int}^{\prime }=0$ and $d{E}_{ad}^{\prime }=0$, Equation (5) reduces to:

(17)$\frac{d\Gamma }{dm}=-\frac{({\Gamma }_{\mathit{eff}}+1)(\Gamma -1)}{M_0+m},$

which describes the evolution of a momentum-conserving (rather than pressure-driven) snowplow. Replacing ${\Gamma }_{\mathit{eff}}\to \Gamma$, the solution of this equation coincides with the result by BM76.

Since the model is based on the homogeneous shell approximation, the adiabatic solution does not recover the correct normalization of the BM76 solution. In this treatment, the total energy of a relativistic decelerating adiabatic blast wave in a power-law density profile $\rho \left(R\right)\propto R^{-s}$ is

(18)$E_0\simeq {\Gamma }_{\mathit{eff}}{E}_{int}^{\prime }\simeq \frac{4}{3}\Gamma E_{int}^{\prime }\simeq \frac{12-4s}{9-2s}{\Gamma }^{2}m{c}^2,$

so that the BM76 normalization can be recovered by multiplying the density of external matter in Equations (5), (13) and (14) by the factor $(9-2s)/(17-4s)$. To smoothly interpolate between the adiabatic regime and the radiative regime, the following correction factor should be adopted:

(19)$C_{BM76,ϵ}\equiv ϵ+\frac{9-2s}{17-4s}(1-ϵ).$

No analytic model properly captures the transition between an adiabatic relativistic blast-wave and the momentum-conserving snowplow, as $ϵ$ increases from zero to unity. The simple interpolation in Equation (19) joins the fully adiabatic BM76 solution with the fully radiative momentum-conserving snowplow.

In summary, Equations (5), (13) and (14), complemented with the correction in Equation (19) (which should by applied to every occurrence of external density and external matter) completely determine the evolution of the shell Lorentz factor $\Gamma$ as a function of the shock radius R.

2.2. Relativistic Shock Acceleration

The spectral shape of the afterglow emission is well described by power-laws over a wide energy range (from radio to GeV-TeV). This is the clear manifestation of the presence of an electron population that has been accelerated in a power-law energy distribution. In GRB afterglows, the main candidate to explain the accelerated non-thermal particles is a Fermi-like mechanism that operates with similar general principles as the non-relativistic diffusive shock acceleration: particles are scattered back and forth across the shock front by magnetic turbulence and gain energy at each shock crossing. The particles themselves are thought to be responsible for triggering the magnetic instability that produces the turbulent field governing their acceleration. The outcome of this acceleration process is determined by the composition of the ambient medium (electron-proton plasma in the case of GRB forward shocks), the fluid Lorentz factor (${\Gamma }_{GRB}>>1$, decreasing to non-relativistic velocity only after several weeks or months) and the magnetization $\sigma$ (i.e., the ratio between Poynting and kinetic flux in the pre-shocked fluid, $\sigma =B^2/(4\pi m_{p}n{c}^2)$), with B being the magnetic field strength. For GRB forward shocks, the magnetization is low, around ${10}^{-9}$ in the interstellar medium and in any case below ${10}^{-5}$ even for a magnetized circumstellar wind.

In this section, we summarize the present understanding of particle acceleration and magnetic field generation in electron-proton, ultra-relativistic, weakly magnetized shocks. The statements and considerations reported in this section refer specifically to this case (which is the one relevant for forward external shocks in GRBs) and might not be valid for magnetized plasma and/or mildly-relativistic flows and/or electron-positron plasma.

In general, the information that one would extract from theoretical/numerical investigations and compare with observations are: (i) the spectral shape of the emitting electrons (i.e., the minimum and maximum Lorentz factor ${\gamma }_{min/max}$ and the spectral index p), (ii) the acceleration efficiency (i.e., the fraction of electrons ${\xi }_e$ and the fraction of energy ${ϵ}_e$ in the non-thermal population) and (iii) the strength of the self-generated magnetic field, usually quantified in terms of fraction ${ϵ}_B$ of the shock-dissipated energy conveyed in the magnetic field. In particular, in order to compare with observations, the relevant ${ϵ}_B$ is the one in the downstream, in the region where radiative cooling takes place and the emission is produced.

After revisiting the state-of-the-art of the theoretical understanding (for recent reviews, see [40,60]), we discuss how particle acceleration and magnetic field amplification are incorporated in GRB afterglow modeling, and then we comment on the constraints on the above-mentioned parameters as inferred from the comparison between the model and the observations.

2.2.1. Inputs from Theoretical Investigations

Analytical approaches and Monte Carlo simulations generally rely on the assumption that electromagnetic waves, providing the scattering centers to regulate and govern the acceleration, are present on both sites of the shock, with a given strength and spectrum, so that the Fermi mechanism can operate. The particle distribution is then evolved under some assumption (such as diffusion in pitch angle) on the scattering process, and considering a test-particle approximation (i.e., the high-energy particles do not modify the properties of the waves).

The main success of these approaches is the verification that under these conditions power-law spectra are indeed produced and the predicted spectral index is in very good agreement with observations of afterglow radiation from GRBs [61]. The spectral index has been calculated for different assumptions on the equation of state, diffusion prescription and for a wide range of shock velocities [61]. A quasi-universal value $p\simeq 2.2-2.3$ is found in the ultra-relativistic limit. Figure 1 shows a comparison between analytical and numerical results as a function of the shock velocity for three different types of shocks (see [61] for details). In the ultra-relativistic limit ($\gamma \beta \gg 1$), the estimates of the spectral slopes converge to a universal value $p=s-2\sim 2.2$.

The investigation of relativistic shocks is complemented by particle-in-cell (PIC) simulations, where the non-linear coupling between particles and self-generated magnetic turbulence is captured from first principles.

The limitations of this technique are imposed by the computation time: for accuracy and stability, PIC simulations need to resolve the electron plasma skin depth $c/{\omega }_{pe}$ of the upcoming electrons (where ${\omega }_{pe}=\sqrt{4\pi e^2{n}_e/m_e}$ is the plasma oscillation frequency of the upstream plasma, $n_e$ is the proper density and $m_e$ is the electron mass), which is orders of magnitudes smaller than the scales of astrophysical interest. It is then difficult to follow the evolution on time-scales and length-scales relevant for astrophysics. Low dimensionality (1D or 2D instead of 3D) and small ion-to-electron mass ratios are additional limitations imposed by the computation time. Results of PIC simulations need then to be extrapolated to bridge the gap between the micro-physical scales and the scales of interest. With these caveats in mind, we summarize here the main achievements.

PIC simulations have shown that magnetic turbulence can be efficiently (${ϵ}_B\sim 0.01-0.1$) generated by the accelerated particles streaming ahead of the shock (in the so-called precursor region), where they generate strong magnetic waves, which in turn scatter the particles back and fourth across the shock. In particular, in the weakly magnetised shocks discussed in this section, the dominant plasma instability is thought to be the so-called Weibel (or current filamentation) instability [62], generated by the counter-streaming of the accelerated particles against the background plasma in the precursor region [42,63]. PIC simulations have demonstrated that as long as the fluid is ultra-relativistic ($\Gamma >5$), the main parameter governing the acceleration is the magnetization $\sigma$, i.e., the efficiency of the process is insensitive to $\Gamma$, as the precursor decelerates the incoming background plasma.

An example of downstream particle spectra derived by PIC simulations is shown in Figure 2 ([42]). The ion and electron spectra are shown for a 2D simulation with $\Gamma =15$, as an ion-to-electron mass ratio $m_i/m_e=25$, and $\sigma ={10}^{-5}$. The temporal evolution is followed up to $t=2500{\omega }_{pi}^{-1}$. The formation of a non-thermal tail is clearly visible.

The downstream non-thermal population is found to include around $\xi \simeq 3\%$ of the electrons, carrying ${ϵ}_e\simeq 10\%$ of the energy. The spectral index is around $p\sim 2.5$. The acceleration proceeds similarly for electrons and ions, since they enter the shock in equipartition (i.e., their relativistic inertia is comparable) as a result of efficient pre-heating in the self-excited turbulence in the precursor.

The maximum energy ${\gamma }_{max}$ increases proportionally to $t^{1/2}$ (see inset in Figure 2), slower than the commonly adopted Bhom rate [64], in which case ${\gamma }_{max}\propto t$. Extrapolating the ${\gamma }_{max}$ behavior to the relevant time-scales and considering that synchrotron cooling will limit the acceleration for high-energy particles, the electron maximal Lorentz factor is found to reach values ${\gamma }_{max}\sim {10}^7$ in the early phase of GRB afterglows, corresponding to synchrotron photon energies around 1 GeV, which is roughly consistent with observations. All these results on the particle spectrum are obtained on time-scales that are too short for the supra-thermal particles to reach a steady-state and their extrapolation to longer time-scales is not trivial.

A still debated open question (because it is computationally demanding) is how the magnetization evolves downstream. PIC simulations have found values of ${ϵ}_B\sim 0.1-0.01$ in the vicinity of the shock front. How this turbulence evolves on longer time-scales is still a matter of debate. The turbulence is expected to decay rapidly, on time-scales orders of magnitude shorter than the synchrotron cooling time. Magnetization is then predicted to be very different close to the shock and in the region where particle cooling takes place. Electrons would then cool in a region of weak magnetic field [65,66]. These considerations suggest that it might not be correct to define a single magnetization ${ϵ}_B$ in GRB modeling, infer its value from observations and compare with predictions from PIC simulations referring to the magnetization near the shock front. Magnetization values inferred from observations most likely probe a region downstream, far from the front shock (see Section 3.3 for a discussion).

Theoretical efforts are fundamental to provide physically motivated inputs for the phenomenological parameters included in the afterglow model. The large number of unknown model parameters, coupled with a limited number of constraints provided by observations, implies that constraints from the theory are of paramount importance for a correct interpretation of the emission in GRBs and for grasping the origin of their non thermal emission, from radio to TeV energies. On the other hand, despite the huge progresses in the theoretical understanding of relativistic acceleration, the theory is not quite yet to the point of providing robust inputs for modeling observations. It is then clear how the two approaches must be combined to gain knowledge on the micro-physics of acceleration and magnetic field generation, on the one hand, and on the origin of radiative processes and macro-physics of the emitting region (bulk Lorentz factor and energy content) of the sources, on the other hand.

2.2.2. Description of Shocks in GRB Afterglow Modeling

The theory of relativistic shock acceleration is applied to the GRB afterglow by introducing several unknown parameters in the model. These are the fractions ${ϵ}_e$ and ${ϵ}_B$ of dissipated energy gained by the accelerated particles and amplified magnetic field, the spectral index p of the accelerated particle spectrum and the fraction ${\xi }_e$ of particles, which efficiently enter the Fermi mechanism and populate the non-thermal distribution.

Recalling that the shock-dissipated energy (in the comoving frame) is given by $d{E}_{sh}^{\prime }=(\Gamma -1)dm{c}^2$ (see Section 2.1), the corresponding energy density is $u_{\mathrm{sh}}^{\prime }=(\Gamma -1){\rho }^{\prime }{c}^2$. From shock jump conditions, the density in the comoving frame ${\rho }^{\prime }$ is related to the density of the unshocked medium (measured in the rest frame) by the equation:

(20)${\rho }^{\prime }=4\Gamma \rho$

which is valid in both the ultra-relativistic and non-relativistic limits (see e.g., [35]). In the GRB afterglow scenario, it is usually assumed that pairs are unimportant and then the density of protons and electrons is the same: $n_p=n_e=n$. This implies that the mass is dominated by protons: $\rho =n{m}_p$. In this case, the available energy density that will be distributed to the accelerated particles (electrons and protons) and to the magnetic field can be expressed as:

(21)$u_{\mathrm{sh}}^{\prime }=4\Gamma (\Gamma -1)n{m}_p{c}^2$

A fraction ${ϵ}_B$ of this energy will be conveyed to the magnetic field:

(22)$u_B^{\prime }={ϵ}_B{u}_{\mathrm{sh}}^{\prime }={ϵ}_{B}4\Gamma (\Gamma -1)n{m}_p{c}^2$

from which it follows that the magnetic field strength $B^{\prime }$ is:

(23)$B^{\prime }=\sqrt{32\pi {ϵ}_B{m}_p{c}^{2}n(\Gamma -1)\Gamma }$

Similarly, for the accelerated electrons:

(24)$u_e^{\prime }={ϵ}_e{u}_{sh}^{\prime }={ϵ}_{e}4\Gamma (\Gamma -1)n{m}_p{c}^2=⟨\gamma ⟩m_e{c}^{2}4\Gamma {\xi }_{e}n$

where $⟨\gamma ⟩$ is the average random Lorentz of the accelerated electrons:

(25)$⟨\gamma ⟩=\frac{{ϵ}_e}{{\xi }_e}\frac{m_p}{m_e}(\Gamma -1)$

The accelerated non-thermal electrons are assumed to have a power-law spectrum as a result of shock acceleration. Their energy distribution can be described by a power-law $N(\gamma )d\gamma \propto {\gamma }^{-p}d\gamma$ for ${\gamma }_{min}\le \gamma \le {\gamma }_{max}$ where ${\gamma }_{min}$ is the minimum Lorentz factor of the injected electrons and ${\gamma }_{max}$ is the maximum Lorentz factor at which electrons can be accelerated. To derive the relation between ${\gamma }_{\min }$, ${\gamma }_{\max }$ and the model parameters, we consider the definition of the average Lorentz factor $⟨\gamma ⟩$:

(26)$⟨\gamma ⟩=\frac{{\int }_{{\gamma }_{min}}^{{\gamma }_{max}}N(\gamma )\gamma d\gamma }{{\int }_{{\gamma }_{min}}^{{\gamma }_{max}}N(\gamma )d\gamma }$

and solve the integrals. Equations (25) and (26) leads to (for $p\ne 1$):

(27)$\begin{array}{c}\left[\frac{ln\left(\frac{{\gamma }_{max}}{{\gamma }_{min}}\right)}{\frac{1}{{\gamma }_{min}}-\frac{1}{{\gamma }_{max}}}\right]=\frac{{ϵ}_e}{{\xi }_e}\frac{m_p}{m_e}(\Gamma -1)\mathrm{if}p=2\\ \left[\frac{p-1}{p-2}\frac{{\gamma }_{min}^{-p+2}-{\gamma }_{max}^{-p+2}}{{\gamma }_{min}^{-p+1}-{\gamma }_{max}^{-p+1}}\right]=\frac{{ϵ}_e}{{\xi }_e}\frac{m_p}{m_e}(\Gamma -1)\mathrm{if}p\ne 2\end{array}$

A simplified equation for ${\gamma }_{min}$ can be obtained assuming that ${\gamma }_{max}^{-p+2}\ll {\gamma }_{min}^{-p+2}$:

(28)${\gamma }_{min}=\frac{{ϵ}_e}{{\xi }_e}\frac{m_p}{m_e}\frac{p-2}{p-1}(\Gamma -1)$

Since p is expected to be $2<p<3$, this condition is verified for ${\gamma }_{max}\gg {\gamma }_{min}$. The minimum Lorentz factor is then not treated as a free parameter of the model, as it is calculated from Equation (28) as a function of the free parameters ${ϵ}_e,{\xi }_e$ and p. Concerning the prescription for the value of ${\gamma }_{max}$ (for details see Section 3.5), it usually relies on the condition that radiative losses between acceleration episodes are equal to the energy gains, where energy gains proceed at the Bhom rate. As we mentioned in the previous section, PIC simulations, however, have shown that this might not be the case.

A similar treatment can be adopted also for protons simply substituting ${ϵ}_e$ with ${ϵ}_p$, $m_e$ with $m_p$ and assuming a power-law energy distribution with spectral index q. As a result, the minimum Lorentz factor for protons can be derived as:

(29)$\begin{array}{c}\left[\frac{ln\left(\frac{{\gamma }_{max,p}}{{\gamma }_{min,p}}\right)}{\frac{1}{{\gamma }_{min,p}}-\frac{1}{{\gamma }_{max,p}}}\right]=\frac{{ϵ}_p}{{\xi }_p}(\Gamma -1)\mathrm{if}q=2\\ \left[\frac{q-1}{q-2}\frac{{\gamma }_{min,p}^{-q+2}-{\gamma }_{ma{x}_p}^{-q+2}}{{\gamma }_{min,p}^{-q+1}-{\gamma }_{max,p}^{-q+1}}\right]=\frac{{ϵ}_p}{{\xi }_p}(\Gamma -1)\mathrm{if}q\ne 2\end{array}$

Solving the equations assuming that ${\gamma }_{max,p}\gg {\gamma }_{min,p}$ leads to:

(30)${\gamma }_{min,p}=\frac{{ϵ}_p}{{\xi }_p}\frac{q-2}{q-1}(\Gamma -1)$

The equations for ${\gamma }_{min}$ and Equation (23) for B, coupled with the description of the blast-wave dynamics described in Section 2.1, provides all the necessary equations to derive the radiative output for a jet with energy E and initial bulk Lorentz factor ${\Gamma }_0$ expanding in a medium with density $n\left(R\right)$. The derivation of the radiative output is detailed in Section 2.3. To conclude the discussion about particle acceleration, in the next section we anticipate which constraints can be inferred on the physics of particle acceleration from multi-wavelength observations, once the afterglow model is adopted.

2.2.3. Constraints to the Acceleration Mechanism Provided by Observations

Assuming that accelerated particles have a power-law spectrum ($d{N}^{\mathrm{acc}}/d\gamma \propto {\gamma }^{-p}$) and the cooling is dominated by synchrotron radiation, the spectral slope p can be inferred from observations of the synchrotron spectrum and/or from the temporal decay of the lightcurves if observations are performed at frequencies higher than the typical frequency ${\nu }_{\mathrm{m}}$ of photons emitted by electrons with the Lorentz factor ${\gamma }_{min}$ (this is correct in both cases of fast and slow cooling regime). The estimated value of p from the afterglow modeling are spread on a wide range, from $p\sim 2$ to $p\sim 3$, suggesting that the spectrum of injected particles does not seem to have a typical slope, at odds with theoretical predictions. The determination of p, however, suffers from the uncertainties on the spectral index inferred from optical and X-ray observations, where the observed spectra are subject to unknown dust and metal absorption. A derivation of p from the decay rate of the lightcurves is also subject to the correct identification of the spectra regime, and partially also to the assumption on the density profile of the external medium, which is often unconstrained (see Section 3.2).

The typical value of ${ϵ}_{\mathrm{e}}$ inferred from the afterglow modeling is around 0.1, meaning that 10% of the shock-dissipated energy is gained by the electrons, spanning from 0.01 to large values, such as 0.8. Although this seems a large uncertainty, ${ϵ}_e$ is perhaps the most well-constrained parameter of the model, and is in good agreement with the values predicted by numerical investigations [42]. For the fraction ${ϵ}_{\mathrm{B}}$, on the contrary, the inferred values varies in a very wide range, typically from ${10}^{-5}$ to ${10}^{-1}$ [67,68,69]. Recent studies that incorporate Fermi-LAT GeV observations [24,70] have demonstrated that the typical values estimated for ${ϵ}_B$ can be even smaller, in the range ∼${10}^{-7}$–${10}^{-2}$. These values are needed in order to model GeV radiation self-consistently with radiation detected at lower frequencies, with repercussions on the estimates of the other parameters, such as n and E. These small values of the ${ϵ}_B$ needed to model the radiation have been tentatively interpreted as the sign of turbulence decay in the downstream [65,66]. As a consequence, even though the turbulence is strong (${ϵ}_B\simeq 0.1$) in the vicinity of the shock where the particle is accelerated, it becomes weaker at larger distances, in the region where particles cool (see Section 3.3). Small values of ${ϵ}_B$ are confirmed by the modeling of recent TeV detections of afterglow radiation from GRBs ([71,72], see Section 4).

Another parameter that one would like to constrain from observations is the fraction of particles ${\xi }_e$ that are injected into the Fermi process. In the vast majority of the studies, this parameter is not included (i.e., it is implicitly assumed that all the electrons are accelerated, ${\xi }_e=1$). This parameter is indeed difficult to constrain, as it is degenerate with all the other parameters [73].

Observations so far have not been able to identify the location of a high-energy cutoff in the synchrotron spectrum that would reveal the maximal energy of the synchrotron photons and then the maximum energy ${\gamma }_{max}$ of the accelerated electrons. Observations by Fermi-LAT are in general consistent with a single power-law extending up to at least 1 GeV. Photons with energies in excess of 1 GeV have been detected from several GRBs, the record holder for Fermi-LAT being a 95 GeV photon [74]. These photons cannot be safely associated to synchrotron radiation on the basis of spectral analysis, as their paucity makes it difficult to assess from spectral analysis whether they are consistent with the power-law extrapolation of the synchrotron spectrum or if they are indicative of the rising of a distinct spectral component. In any case, the Fermi-LAT detections are suggesting that synchrotron photons should be produced at least up to a few GeV. This is consistent with the limit commonly invoked for particle acceleration: if the acceleration proceeds at the Bhom rate ($t_{acc}\simeq r_L/c$) with $r_L=E/eB$ being the Larmor radius, and is limited by synchrotron cooling ($t_{syn}\simeq 6\pi m_{e}c/{\sigma }_T{B}^2\gamma$) then ${\gamma }_{max}\sim {10}^7-{10}^8$ can be reached. Even though this does not necessarily imply that acceleration must proceed at the Bhom limit, the value of ${\gamma }_{max}$ inferred from the detection of GeV photons is quite large and barely consistent with what is found by PIC simulations. Whether or not the observations are in tension with the present derivation of ${\gamma }_{max}$ from PIC simulations and theoretical arguments, strongly depends on a clear identification of the origin of photons in the GeV-TeV energy range. Present and future observations with Imaging Atmospheric Cherenkov Telescopes (IACTs) are the main candidates to shed light on this issue.

2.3. Derivation of the Radiative Output

The expected radiative output can be estimated by means of analytical approximations, which provide prescriptions for the location of the synchrotron self-absorption frequency ${\nu }_{sa}$, the characteristic frequency ${\nu }_m$ emitted by electrons with Lorentz factor ${\gamma }_{min}$, the cooling frequency ${\nu }_c$ emitted by electrons with Lorentz factor ${\gamma }_c$ and the overall synchrotron flux [38,47]. In these approaches, the synchrotron spectrum is in general approximated with power-laws connected by sharp breaks, but more sophisticated analytical approximations of numerically derived synchrotron spectra have also been proposed [54]. The associated SSC component in the Thomson regime [48] and corrections to be applied to the synchrotron and SSC spectra to account for the effects of the Klein–Nishina [75] cross section (see Section 2.3.1) are also available in the literature. These prescriptions are usually developed for the deceleration phase, when the Blandford–McKee solution [35] for the blast-wave dynamics applies, i.e., as long as the blast-wave is still relativistic. These models take as input parameters the kinetic energy content of the blast-wave $E_k$, the external density $n\left(R\right)=n_0{R}^{-s}$ (with $s=0$ or $s=2$), the fraction of shock-dissipated energy gained by electrons (${ϵ}_e$) and by the amplified magnetic field (${ϵ}_B$), and the spectral slope of the accelerated electrons p. During the deceleration phase, the initial bulk Lorentz factor ${\Gamma }_0$ does not play any role, but its value determines the radius (or time) at which the deceleration begins.

An alternative approach to estimate the expected spectra and their evolution in time consists in numerically solving the differential equation describing the evolution of the particle spectra and estimating the associated emission [76,77,78]. In this section, we describe a radiative code that simultaneously solves the time evolution of the electron and photon distribution. The code has been adopted, e.g., for the modeling of GRB 190114C presented in [71].

The temporal evolution of the particle distribution is described by the differential equation:

(31)$\frac{\partial N(\gamma ,t^{\prime })}{\partial t^{\prime }}=\frac{\partial }{\partial \gamma }\left[\dot{\gamma }N(\gamma ,t^{\prime })\right]+Q(\gamma ),$

where $\dot{\gamma }=\frac{\partial \gamma }{\partial t^{\prime }}$ is the rate of change of the Lorentz factor $\gamma$ of an electron caused by adiabatic, synchrotron and SSC losses and energy gains by synchrotron self-absorption. In the SSC mechanism, the synchrotron photons produced by electrons in the emission region act as seed photons that are up-scattered at higher energies by the same population of electrons. Such a scenario will generate a very high energy spectral component, which is the target of searches by IACTs such as MAGIC¹ and H.E.S.S.². In principle, up-scattering of an external population of seed photons can also be considered and included in the cooling term, but here we will ignore this mechanism (external Compton) and focus only on SSC. The source term $Q(\gamma ,t^{\prime })=Q^{acc}(\gamma ,t^{\prime })+Q^{pp}(\gamma ,t^{\prime })$ describes the injection of freshly accelerated particles ($Q^{acc}(\gamma ,t^{\prime })=d{N}^{acc}/d\gamma d{t}^{\prime }$) and the injection of pairs $Q^{pp}(\gamma ,t^{\prime })$ produced by photon–photon annihilation.

In the next sections, we explicitly define each one of the terms included in Equation (31) and how to estimate the synchrotron and SSC emission. To solve the equation, an implicit finite difference scheme based on the discretization method proposed by [79] can be adopted.

2.3.1. Synchrotron and SSC Cooling

The synchrotron power emitted by an electron with Lorentz factor $\gamma$ depends on the pitch angle, i.e., the angle between the electron velocity and the magnetic field line. In the following, we assume that the electrons have an isotropic pitch angle distribution and we use equations that are averaged over the pitch angle (e.g., [80]). The synchrotron cooling rate of an electron with Lorentz factor $\gamma$ is given by:

(32)${\dot{\gamma }}_{syn}\equiv \frac{d\gamma }{d{t}^{\prime }}{|}_{syn}=-\frac{{\sigma }_T{\gamma }^2{B}^{\prime 2}}{6\pi m_{e}c}$

The cross section for the inverse Compton mechanism is constant and equal to the Thomson cross section (${\sigma }_T$) as long as the photon energy in the frame of the electron is smaller than the rest mass electron energy $m_e{c}^2$. For higher photon energies, the cross section decreases as a function of the energy and is described by the Klein-Nishina (KN) cross section. To estimate SSC losses, we adopt the formulation proposed in [81], which is valid for both regimes. Defining the SSC kernel as:

(33)$K(\gamma ,{\nu }^{\prime },\widetilde{{\nu }^{\prime }})=\left\{\begin{array}{c}\frac{\epsilon }{\tilde{\epsilon }}-\frac{1}{4{\gamma }^2}\frac{\tilde{\epsilon }}{4{\gamma }^2}<\epsilon <\tilde{\epsilon }\hfill \\ 2qlnq+(1+2q)(1-q)+\frac{1}{2}(1-q)\frac{{(4\gamma \tilde{\epsilon }q)}^2}{(1+4\gamma \tilde{\epsilon }q)}\tilde{\epsilon }<\epsilon <\frac{4{\gamma }^2\tilde{\epsilon }}{1+4\gamma \tilde{\epsilon }},\hfill \end{array}\right.$

where:

(34)$\tilde{\epsilon }=\frac{h\widetilde{{\nu }^{\prime }}}{m_e{c}^2}\epsilon =\frac{h{\nu }^{\prime }}{m_e{c}^2}q=\frac{\epsilon }{4\gamma \tilde{\epsilon }(\gamma -\epsilon )}.$

$\tilde{\epsilon }$ and $\epsilon$ are the energies of the photons (normalized to the rest mass electron energy) before and after the scattering process, respectively. The two terms of Equation (33) account, respectively, for the down-scattering (i.e., $\epsilon <\tilde{\epsilon }$) and the up-scattering (i.e., $\epsilon >\tilde{\epsilon }$) process. The energy loss term for the SSC can now be calculated with the equation:

(35)${\dot{\gamma }}_{SSC}=\frac{d\gamma }{d{t}^{\prime }}{|}_{SSC}=-\frac{3h{\sigma }_t}{4m_{e}c{\gamma }^2}\int d{\nu }^{\prime }{\nu }^{\prime }\int \frac{d\widetilde{{\nu }^{\prime }}}{\widetilde{{\nu }^{\prime }}}{n}_{\widetilde{{\nu }^{\prime }}}\left(t^{\prime }\right)K(\gamma ,{\nu }^{\prime },\widetilde{{\nu }^{\prime }}).$

2.3.2. Adiabatic Cooling

As discussed in Section 2.1, particles lose their energy adiabatically due to the spreading of the emission region. This energy loss term should be inserted in the kinetic equation governing the evolution of the particle distribution. To derive the adiabatic losses, we rewrite Equation (10) as a function of energy losses $d\gamma$ in a comoving time $d{t}^{\prime }$:

(36)${\dot{\gamma }}_{ad}=\frac{d\gamma }{d{t}^{\prime }}{|}_{ad}=-\frac{\gamma {\beta }^2}{3}\frac{dln{V}^{\prime }}{d{t}^{\prime }},$

with $\beta$ being the random velocity of particles in a unit of c. The comoving volume $V^{\prime }$ of the emission region can be estimated considering that the contact discontinuity is moving away from the shock at a velocity $c/3$. After a time $t^{\prime }=\int dR/\Gamma \left(R\right)c$, the comoving volume is:

(37)$V^{\prime }=4\pi R^2\frac{c{t}^{\prime }}{3},$

and:

(38)${\dot{\gamma }}_{ad}=\frac{d\gamma }{d{t}^{\prime }}{|}_{ad}=-\frac{\gamma {\beta }^2}{3}\left(\frac{2\Gamma c}{R}+\frac{1}{t^{\prime }}\right).$

2.3.3. Synchrotron Self-Absorption (SSA)

Electrons can re-absorb low energy photons before they escape from the source region. The absorption coefficient ${\alpha }_{\nu }$ can be expressed as [80]:

(39)${\alpha }_{\nu }=-\frac{1}{8\pi {\nu }^{\prime 2}{m}_e}\int d\gamma P^{\prime }(\gamma ,{\nu }^{\prime }){\gamma }^2\frac{\partial }{\partial \gamma }\left[\frac{N(\gamma )}{{\gamma }^2}\right]$

valid for any radiation mechanism at the emission frequency ${\nu }^{\prime }$, with $P^{\prime }(\gamma ,{\nu }^{\prime })$ being the specific power of electrons with Lorentz factor $\gamma$ at frequency ${\nu }^{\prime }$ and assuming $h{\nu }^{\prime }\ll \gamma m_e{c}^2$. Thus, the SSA mechanism will mostly affect the low frequency range. This results in a modification of the lower frequency tail of the synchrotron spectrum as:

(40)$P_{\nu }^{\prime }\propto \left\{\begin{array}{c}{\nu }^{\prime 5/2}{\nu }_i^{\prime }<{\nu }^{\prime }<{\nu }_{SSA}^{\prime }\hfill \\ {\nu }^{\prime 2}{\nu }^{\prime }<{\nu }_{SSA}^{\prime }<{\nu }_i^{\prime },\hfill \end{array}\right.$

assuming a power-law distribution of electrons, with ${\nu }_i^{\prime }=min({\nu }_m^{\prime },{\nu }_{cool}^{\prime })$ and ${\nu }_{SSA}^{\prime }$ the frequency below which the synchrotron flux is self-absorbed and the source becomes optically thick.

2.3.4. Synchrotron and Inverse Compton Emission

Following [82], the synchrotron spectrum emitted by an electron with Lorentz factor $\gamma$, averaged over an isotropic pitch angle distribution, is:

(41)$P_{{\nu }^{\prime }}^{\prime syn}({\nu }^{\prime },\gamma )=\frac{2\sqrt{3}{e}^3{B}^{\prime }}{m_e{c}^2}{x}^2\left[K_{4/3}\left(x\right)K_{1/3}\left(x\right)-0.6x(K_{4/3}^2\left(x\right)-K_{1/3}^2\left(x\right))\right],$

where $x\equiv {\nu }^{\prime }4\pi m_{e}c/\left(6e{B}^{\prime }{\gamma }^2\right)$ and $K_n$ are the modified Bessel functions of order n. The total power emitted at the frequency ${\nu }^{\prime }$ is obtained integrating over the electron distribution:

(42)$P_{{\nu }^{\prime }}^{\prime syn}\left({\nu }^{\prime }\right)=\int P_{{\nu }^{\prime }}^{\prime syn}({\nu }^{\prime },\gamma )\frac{dN}{d\gamma }d\gamma .$

The SSC radiation emitted by an electron with Lorentz factor $\gamma$ can be calculated as:

(43)$P_{{\nu }^{\prime }}^{\prime SSC}({\nu }^{\prime },\gamma )=\frac{3}{4}h{\sigma }_{T}c\frac{{\nu }^{\prime }}{{\gamma }^2}\int \frac{d\widetilde{{\nu }^{\prime }}}{\widetilde{{\nu }^{\prime }}}{n}_{\widetilde{{\nu }^{\prime }}}K(\gamma ,{\nu }^{\prime },\widetilde{{\nu }^{\prime }}),$

where $n_{\widetilde{{\nu }^{\prime }}}$ is the photon density of synchrotron photons and the integration is performed over the entire synchrotron spectrum. Integration over the electron distribution provides the total SSC emitted power at frequency ${\nu }^{\prime }$.

2.3.5. Pair Production

Pair production by photon-photon annihilation is particularly important for a correct estimate of the radiation spectrum in the GeV-TeV band. Indeed, some of the emitted VHE photons are lost due to their interaction with photons at lower energies (typically X-ray photons). As a result, the observed flux is attenuated and the resulting spectrum at VHE is modified. Here we follow the treatment presented in [83]. The cross section of the process ${\sigma }_{\gamma \gamma }$ as a function of ${\beta }^{\prime }$, the centre-of-mass speed of the electron and positron is given by:

(44)${\sigma }_{\gamma \gamma }\left({\beta }^{\prime }\right)=\frac{3}{16}{\sigma }_T(1-{\beta }^{\prime 2})[(3-{\beta }^{\prime 4})ln\left(\frac{1+{\beta }^{\prime }}{1-{\beta }^{\prime }}\right)-2{\beta }^{\prime }(2-{\beta }^{\prime 2})]$

where:

(45)${\beta }^{\prime }({\omega }_t,{\omega }_s,\mu )={\left[1-\frac{2}{{\omega }_t{\omega }_s(1-\mu )}\right]}^{\frac{1}{2}}$

and ${\omega }_t=h{\nu }_t^{\prime }/m_e{c}^2$ with ${\nu }_t^{\prime }$ being the target photon frequency, ${\omega }_s=h{\nu }^{\prime }/m_e{c}^2$ with ${\nu }^{\prime }$ being the source photon frequency and $\mu =cos\varphi$, where $\varphi$ is the scattering angle. Then, it is possible to derive the annihilation rate of photons into electron-positron pairs as:

(46)$R({\omega }_t,{\omega }_s)=c{\int }_{-1}^{{\mu }_{max}}\frac{d\mu }{2}(1-\mu ){\sigma }_{\gamma \gamma }({\omega }_t,{\omega }_s,\mu ),$

where ${\mu }_{max}=max(-1,1-2/{\omega }_s{\omega }_t)$ coming from the requirement ${\beta }^{\prime 2}>0$. Considering $x={\omega }_t{\omega }_s$, it is possible to derive asymptotic limits for $R({\omega }_t,{\omega }_s)\equiv R\left(x\right)$ in two regimes. For $x\to 1$ (i.e., near the threshold condition) $R\left(x\right)\to c{\sigma }_T/2{(x-1)}^{3/2}$, while for $x\gg 1$ (i.e., ultra-relativistic limit) $R\to \frac{3}{4}c{\sigma }_{T}lnx/x$. An accurate and simple approximation, which takes into account both regimes, is given by:

(47)$R\left(x\right)\approx 0.652c{\sigma }_T\frac{x^2-1}{x^3}ln\left(x\right)\mathrm{H}(x-1),$

where $H(x-1)$ is the Heaviside function [83]. The approximation accurately reproduces the behavior near the peak at $x_{peak}\sim 3.7$ and over the range $1.3<x<{10}^4$, which usually dominates during the calculations. A comparison between Equations (46) and (47) is given in Figure 3, where the goodness of the approximation adopted in the mentioned x range can be observed.

The impact of the flux attenuation due to pair production mechanism on the GRB spectra is estimated in terms of the optical depth value ${\tau }_{\gamma \gamma }$. From its definition:

(48)${\tau }_{\gamma \gamma }\left({\nu }^{\prime }\right)={\sigma }_{{\nu }^{\prime }{\nu }_t^{\prime }}{n}^{\prime }\left({\nu }_t^{\prime }\right)\Delta R^{\prime }$

where $n^{\prime }\left({\nu }_t^{\prime }\right)$ is the number density of the target photons per unit of volume, ${\sigma }_{{\nu }^{\prime }{\nu }_t^{\prime }}$ is the cross section and $\Delta R^{\prime }$ is the width of the emission region. Introducing the cross section in terms of the annihilation rate $R\left(x\right)$ in its approximated formula and integrating over all the possible target photon frequencies:

(49)${\tau }_{\gamma \gamma }\left({\nu }^{\prime }\right)=\frac{\Delta R^{\prime }}{c}\int R({\nu }^{\prime },{\nu }_t^{\prime })n_{{\nu }^{\prime }}^{\prime }\left({\nu }_t^{\prime }\right)d{\nu }_t^{\prime }$

where ${\nu }^{\prime }$ and ${\nu }_t^{\prime }$ are the frequencies of the source and of the target interacting photons. The pair production attenuation factor can be then introduced simply multiplying the flux by a factor $(1-e^{-{\tau }_{\gamma \gamma }})/{\tau }_{\gamma \gamma }$. This attenuation factor will modify the GRB spectrum, giving a non-negligible contribution in the VHE domain, in particular. An example of the modification of a GRB spectrum due to pair production can be observed in Figure 4. Here, the flux emitted in the afterglow external forward shock scenario by synchrotron and SSC radiation and the flux attenuation due to pair production have been calculated with a numerical code. For a set of quite standard afterglow parameters and assuming $\Delta R^{\prime }=R/\Gamma$, the attenuation of the observed flux due to pair production become relevant above 0.2 TeV, and it reduces the flux by ∼ 30% at 1 TeV and by ∼ 70% at 10 TeV.

Similar considerations can also be conducted for the electron/positron production. Assuming that the electron and positron arises with equal Lorentz factor $\gamma$ and that $x_{peak}\sim 3.7$, a photon with energy ${\omega }_s\gg 1$ will mostly interact with a target photon of energy ${\omega }_t\approx 3/{\omega }_s$. Then, from the energy conservation condition:

(50)$2\gamma ={\omega }_s+\frac{3}{{\omega }_s}\approx {\omega }_s=\frac{h{\nu }^{\prime }}{m_e{c}^2}$

The $e^{\pm }$ production can be observed as an additional source term for the distribution of accelerated particles. As a result, an additional injection term $Q_e^{pp}$ to be inserted in the kinetic equation (Equation (31)) is calculated as:

(51)$Q_e^{pp}(\gamma ,t^{\prime })=4\frac{m_e{c}^2}{h}{n}_{{\nu }^{\prime }}(\frac{2\gamma m_e{c}^2}{h},t^{\prime })\int d{\nu }_t^{\prime }{n}_{{\nu }^{\prime }}({\nu }_t^{\prime },t^{\prime })R(\frac{2\gamma m_e{c}^2}{h},{\nu }_t^{\prime })$

2.3.6. Comparison with Analytical Approximations

In order to compare results from the numerical method described in the previous section and analytical prescriptions available in the literature, we give an example in Figure 5. The analytical prescriptions for the synchrotron and the SSC component are calculated following [38,48]. In [38], the synchrotron spectra and light-curves are derived assuming a power-law distribution of electrons in an expanding relativistic shock, cooling only by synchrotron emission. The dynamical evolution is described following BM76 equations for an adiabatic blast-wave expanding in a constant density medium. The resulting emission spectrum (green dashed lines in Figure 5) is described with a series of sharp broken power-laws. The SSC component associated to the synchrotron emission was computed, as a function of the afterglow parameters, in [48]. In this work, calculations are performed assuming that the scatterings occur in Thomson regime. Modifications to the synchrotron spectrum caused by strong SSC electron cooling are also detailed.

From the comparison proposed in Figure 5, it can be clearly observed that analytical and numerical results are in general in good agreement. Both curves follow the same behavior except for the high-energy part of the SSC component. Here, the KN scattering regime, which is not taken into account in the analytical approximation, becomes relevant. As a result, the numerical calculations differ from the analytical ones showing a peak and a cutoff in the SSC spectrum due to the KN effects.

Nevertheless, there are a few minor discrepancies between the two methods. The numerically-derived spectrum is very smooth around the break frequencies, with the result that the theoretically expected slope (e.g., the one predicted by the analytical approximations) is reached only in regions of the spectrum that lie far from the breaks, i.e., is reached only asymptotically. This puts into questions simple methods for discriminating among different regimes and different density profiles based on closure relations, which are relations between the spectral and the temporal decay indices [10,38]. Regarding the flux normalization, there are minor discrepancies between the numerical and analytical results. This is due to the fact that in analytical prescriptions it is assumed that the radiation is entirely emitted at the characteristic synchrotron frequency. On the contrary, in the numerical derivation, the full synchrotron spectrum of a single electron is summed up over the whole electron distribution. Similar considerations apply to the SSC component when comparing with the analytical spectra. Moreover, the discrepancies observed between analytical and numerical SSC spectra are amplified by the differences observed in the target synchrotron spectra.

In general, this comparison shows that the numerical treatment is a powerful tool able to predict the multi-wavelength GRB emission in a more accurate way than the analytical prescriptions. The latter ones, however, are still giving valid approximations of the overall spectral shape. The main limitation of analytical estimates arises when TeV observations are involved. The importance of KN corrections is evident in this band and should be properly treated for a correct interpretation of the TeV spectra, as will be shown in Section 4.

3. Open Questions

As predicted by the basic standard model presented in the previous section, the afterglow emission is the result of particle acceleration and radiative cooling occurring in two different regions: the forward and the reverse shock. The temporal and spectral behavior of the two emission components can be inferred after the jet/blastwave dynamics, acceleration mechanisms and the radiation processes are modeled (Section 2). The general agreement between model predictions and observations convincingly proves that the long-lasting radio-to-GeV radiation is indeed produced in interactions between the ejecta and the external medium. Moreover, the radiative mechanisms involved and the nature of emitting particles are well established, with synchrotron (and possibly SSC) from the accelerated electrons (either at the forward or reverse shock) being the source of the detected radiation.

Despite the general success of the external shock scenario, there are several, longstanding open issues which represent a serious challenge for our present understanding of the afterglow emission and the GRB phenomenon in general. Moreover, even when observations seem to be in qualitative agreement with predictions, the extraction of the model parameters (which would give important feedbacks on our understanding of particle acceleration and GRB environments) is limited by the large degeneracy among parameters and lack of solid inputs from theoretical considerations.

Afterglow emission studies have not experienced relevant progresses in the last years, with observations and techniques that are the same since the launch of the Swift satellite. The recent discovery of TeV radiation from GRBs is opening the possibility to renovate and boost afterglow studies, with major impacts on the general understanding of GRB sources.

In this section, we list and comment on those aspects still lacking a clear explanation, and in particular we selected topics which might largely benefit from observations and detections in the VHE regime.

3.1. X-ray Flares

Observations of the afterglow emission in the X-ray and optical often display behaviors that are not predicted by the standard scenario, and require the inclusion of additional emission components contributing to the detected radiation. In the standard external forward shock scenario, the afterglow light-curves in the X-ray and optical band are expected to decay following a power-law or a broken power-law behavior, where the breaks are interpreted as the cooling or injection frequency crossing the observed band [38,47,54]. The advent of Swift-XRT and the increasing number of optical follow-up observations performed by ground-based robotic telescopes have highlighted the presence, in a good fraction of cases, of unexpected features in the early time afterglow, such as flares and plateaus [49,57].

Flares are episodes of sudden rebrightenings characterized by a very fast rise of the flux, followed by an exponential decay profile. Comprehensive studies of X-ray afterglows demonstrate that an X-ray flare is observed in ∼33% of the GRBs [84,85]. The times at which they are observed span a very wide range, from around ∼30 s up to ∼${10}^6$ s after the trigger time. The time where the flare peaks is shown in Figure 6 ($T_{pk}$, x-axis) for a large sample of 468 X-ray flares in long GRBs. Most of the flares occur within ${10}^3$ s, even though there are many cases of flares occurring several hours after the burst. The width of the flare $\omega$ is found to evolve linearly with time to larger values following the trend $\omega \sim 0.2T_{pk}$ [84]. The average and peak luminosities L of the flare also display a dependence from $T_{pk}$, with $L\propto T_{pk}^{-2.7}$ at least for early time ($T_{pk}<{10}^3$ s) flares [84,85]. When also including late time flares [86,87], a shallower index is obtained, around ∼$-1.2$. The energy emitted during flare episodes is quite large and, for early time flares, is around ∼10% of the prompt emission or sometimes even comparable [88].

Flares have also been detected in the optical, although the sample of optical flares is far smaller than the X-ray one [89]. A statistical study of optical flares detected by Swift/UVOT demonstrates that most of them correlate with and share similar temporal properties to simultaneous X-ray flares. Nevertheless, there are a few dozen of GRBs for which no X-ray flaring activity is observed simultaneously with optical flares [89].

Flares are believed to have an inner origin and to be associated with a prolonged activity of the GRB central engine [57,90,91,92,93,94]. However, the relatively long timescales on which they are detected represent a challenge for the model. Many questions are still open, such as the location of the emitting region, what is powering the flares, and whether late time flares have a different origin than flares detected at early times.

Speculations about possible signatures of X-ray flares in the GeV-TeV range are present in the literature [95,96,97,98]. Assuming that flares have an internal origin and are produced at $R<R_{dec}$, forward shock electrons will be exposed to the flare radiation, producing an IC emission component by up-scattering the flare photons. Following these estimates, the IC component peaks at $\sim 100$ GeV and has a flux comparable to the X-ray flux. Alternatively, GeV-TeV radiation associated to flares can be produced by the SSC mechanism, where electrons responsible for X-ray synchrotron flares also upscatter these photons to higher energies. The process is considered less interesting for TeV radiation because the peak of this SSC component is expected to be around 1 GeV [95], due to a relatively low minimum Lorentz factor ${\gamma }_{min}\sim 60$. Such value is estimated from theoretical considerations where ${\gamma }_{min}\simeq 60{ϵ}_{e,-1}({\Gamma }_{IS}-1)$ for $p=2.5$, ${ϵ}_e=0.1$ and a relative shock Lorentz factor ${\Gamma }_{IS}$ of the order of unity. We notice that the recent estimates of the minimum electron Lorentz factor in the late prompt emission of GRBs [29] may modify these predictions, and place the expected SSC around 100 GeV. The luminosity of this component will strongly depend on the size of the emitting region. As a result, the detection of flares in GeV-TeV band can provide relevant information to identify the properties at the emitting region and the production site of the flaring activity.

To understand what the chances of current and future VHE ground-based instruments are in contributing to the study of flares, we perform some simplified estimates. The MAGIC telescopes observed 138 GRBs in almost ∼16.5 years, from 2005 up to June 2021 [99]. More than half of them (74 events) have been observed with delays from shorter than ${10}^3$ s, which means ∼4.5 GRBs yr${}^{-1}$, and 37 events observed with delays shorter than 100 s (i.e., 2.2 GRB yr${}^{-1}$). Considering that ∼33% of the long GRBs have an X-ray flare and considering the distribution of their peak times (see Figure 6), we estimate that ∼1 GRBs/yr is the rate of GRBs with an X-ray flare occurring during MAGIC observations.

Let us go a bit further and estimate the detectability of a putative ∼${10}^2$ GeV counterpart of X-ray flares. For the flux of the GeV-TeV flare, we consider as reference value the X-ray flux, and discuss what happens if a similar or ten times smaller flux is emitted at ∼${10}^2$ GeV.

We collect the X-ray flux of a large sample of flares from the catalog of X-ray flares presented in [87]. The results are shown in Figure 6. The average flux of the flare and the flare peak time correlate, and the orange line represents the best fit. To perform the estimates, we consider two different flare peak times, $T_{pk}={10}^2$ s and $T_{pk}={10}^3$ s. The typical average fluxes at those times are $F=1\times {10}^{-9}$ erg cm${}^{-2}$ s${}^{-1}$ and $F=1\times {10}^{-10}$ erg cm${}^{-2}$ s${}^{-1}$, respectively. Assuming that a similar amount of flux is emitted around 100 GeV we can compare these values with the differential sensitivity as a function of the observing time of IACT instruments. Figure 7 [100] shows the sensitivity for several telescopes to the detection of a point-like source at five standard deviations significance as a function of the exposure time and for four selected energies. Considering that the width of the flare is related to the peak time following the relation $\omega \sim 0.2T_{pk}$, we can compare the flare fluxes estimated at $T_{pk}={10}^2$ s and $T_{pk}={10}^3$ s with the differential sensitivity for observing the time of $t_{obs}=20$ s and $t_{obs}=200$ s. The flare fluxes lie close to the differential sensitivity of the MAGIC telescopes (for 100 GeV at $t_{obs}=20$ s is ∼$1.0\times {10}^{-9}$ erg cm${}^{-2}$ s${}^{-1}$ and at $t_{obs}=200$ s is ∼$5.0\times {10}^{-10}$ erg cm${}^{-2}$ s${}^{-1}$). This indicates that MAGIC telescopes can barely detect such a flare. Moreover, Extragalactic Background Light (EBL) attenuation reduces the flux, which is why we are making the estimates at 100 GeV, where the attenuation is still small. We conclude that MAGIC would be able to detect (or place relevant constraints on) only the brightest X-ray flares (as can be observed in Figure 6, the correlation has a large spread, and flares at ${10}^2$ or ${10}^3$ s can easily have fluxes one order of magnitude larger than what is assumed here).

Concerning future instruments, the Cherenkov Telescope Array (CTA³) will have a sensitivity which is almost one order of magnitude lower than the MAGIC one and similar slewing capabilities. The same estimates performed for MAGIC can be applied to CTA, with the advantage that CTA will have a northern and southern sites, approximately doubling the possibility to follow GRBs within short time-scales. This is a promising indication that the CTA array will be potentially able to detect possible counterparts at $E\sim 30$ the GeV of X-ray flares, provided that this counterpart has a flux that is no less than ten times smaller than what is detected in X-rays. As a result, it can play a major role in exploring and improving our knowledge of flares and their connection with prompt emission and with the prolonged activity of the central engine.

3.2. Density Profile of the External Medium

Following the established connection between long GRBs and the core-collapse of massive stars, the jet is expected to produce the afterglow while propagating in the wind of the star in its free-streaming phase. Afterglow radiation of long GRBs should then be produced in the interaction with a medium with a radial density profile $n\propto R^{-2}$. However, several investigations have demonstrated that about half of the long GRBs are better explained if the blast-wave is assumed to run into a medium with constant density. We revise the evidence in support of the constant density medium and discuss the difficulties in reconciling these observations with expectations on the environment surrounding long GRB progenitors.

Long GRBs originate from the core-collapse of massive stars, most likely rapidly rotating, and with a possible evidence of a preference for low-metallicity. The most convincing evidence in support of this paradigm is the association with type Ic supernovae and the proximity of GRBs to young star-forming regions. While the connection of long GRBs (or at least with the bulk of the population) with the core-collapse of massive stars is solid, the role of metallicity and rotation in the launch of the GRB jets, and the identification of the progenitor star, is still uncertain. The progenitor is usually identified with Wolf-Rayet stars, massive stars ($M>20M_{\odot }$) in the final stages of their evolution, characterized by powerful winds and a high mass loss rate [101]. The wind from the star is expected to interact and deeply modify the environment where the GRB explodes and leaves imprints on its afterglow emission.

More in detail, the interaction between the stellar wind and the ISM four concentric regions with different properties is expected to form. In the inner part (i.e., close to the star), the circumburst medium is permeated by the free-streaming wind, producing a density with radial profile $n\propto R^{-2}$. The density is related to the mass loss rate $\dot{M}$ and to the velocity $v_w$ of the free-streaming stellar wind by:

(52)$n\left(R\right)=\frac{\dot{M}}{4\pi R^2{m}_p{v}_w}.$

A termination shock separates the unshocked from the shocked wind: the latter forms a hot bubble of thermalized wind material, with a nearly constant density profile, as the formation of pressure and density gradients is prevented by the high sound speed inside the bubble. The hot bubble, in its outer part, is enclosed by a shell of shocked ISM, surrounded by the unshocked ISM. The GRB jet is supposed to trill its way in this stratified medium [10].

To understand where most of the afterglow evolution occurs, we have to estimate the deceleration radius $R_d$ and the non-relativistic radius $R_{NR}$ (i.e., the radius where the blast-wave has decelerated to non-relativistic velocity) and compare them to the termination shock radius. For typical parameters ($\dot{M}={10}^{-5}{M}_{\odot }$ yr${}^{-1}$ and $v_w={10}^3$ km s${}^{-1}$), the fit to numerical models of Wolf-Rayet stars [102] give the following relation between the termination shock radius and the density of the unshocked ISM: $R_T=10n_{ISM}^{-1/2}$ pc, where $n_{ISM}$ is the density of the unshocked ISM. From the blast-wave dynamics, the deceleration and the non-relativistic radius are $R_d=6\times {10}^{-5}{E}_{52}{v}_{w,3}/\left({\dot{M}}_{-5}{\Gamma }_{0,2}\right)$ pc and $R_d=0.6E_{52}{v}_{w,3}/{\dot{M}}_{-5}$ pc, respectively. It is evident how the complete evolution of the afterglow radiation occurs well inside the free-streaming region.

In afterglow modeling of long GRBs it is then customary to assume a density profile described by Equation (52), where $\dot{M}$ and $v_w$ are treated as unknown parameters (normalized to the typical values of a Wolf-Rayet star) combined in one single free model parameter $A_{☆}$: $n\left(R\right)=3\times {10}^{35}{A}_{☆}{R}^{-2}$. Despite this robust prediction, the modeling of afterglow observations shows that in a relevant fraction of cases, observations are better explained by adopting a circumburst medium with a constant density $n=n_0$.

The fraction of this case varies depending on the method and on the selected sample, and is on average about 50% [103,104,105,106].

To place the termination shock at least inside the non-relativistic radius, one should invoke a very large density of the ISM, $n\gtrsim {10}^5$ cm${}^{-3}$, typical of dense cores of molecular clouds: $R_T=0.03(n_{ISM}/{10}^5$ cm${}^{-3}{)}^{-1/2}$ pc. Density profiles for different ISM densities are shown in Figure 8, upper panel. Alternatively, one can try to variate the wind parameters. How the termination shock radius changes for different values of $\dot{M}$ and $v_w$ is shown in the bottom panel of Figure 8. A very low mass loss rate $\dot{M}={10}^{-7}{M}_{\odot }$yr${}^{-1}$ (which may find a justification in case of a low-metallicity star) is needed to bring the termination shock radius below 1 pc (for $n_{ISM}=10$ cm${}^{-3}$). With this low mass loss rate, the deceleration and non-relativistic radius increase ($R_d\sim 6\times {10}^{-3}$ pc and $R_{NT}\sim 60$ pc), placing the termination shock still after the deceleration radius but well within the non-relativistic radius, allowing for part of the observed emission to develop into a constant density environment. By increasing the blast-wave energy, the deceleration radius can further approach $R_T$. This suggests that it is more likely for a very energetic GRB to cross the termination shock at early times and then expand in a ISM-like medium, as compared to a faint GRB. An indication of an average larger $E_{\gamma }$ in GRBs with a wind-like medium as compared to GRBs with a ISM-like medium has been found in [107], but is in contrast with results from the study performed by [106] on a larger sample.

The parameter space for which part of the afterglow emission can indeed be produced in the ISM-like density profile of the shocked wind is very limited, as it corresponds to the most energetic GRBs, low-metallicity progenitors and high-density ISM, or a combination of these factors [107]. These considerations on the diversity of $E_k,\dot{M},v_w$ and ISM density may not be sufficient to explain the results of the modeling (i.e., the preference for a ISM-like environment). The fraction of GRBs, which might have these peculiar parameters can hardly account for the large fraction of GRBs for which a wind-like profile is excluded by observations. The required conditions are too extreme to be verified in half of the population. However, it is not clear if this percentage has been overestimated by present studies. To quantify the inconsistency, the first step would be to perform a dedicated study of afterglow emission to assess the percentage of long GRB afterglows that are not consistent with a wind-like environment.

Methods based on closure relations may not be valid if the spectrum is modified by Compton scattering in the Klein–Nishina regime (see also [108]). Moreover, these are based on a simple approximation of the synchrotron spectrum into power-law segments, while the wide curvature of the real synchrotron spectra might lead to incorrect estimates of the value of p if the observed frequency is in the vicinity of a synchrotron break frequency. A full modeling is then necessary to really assess the fraction of long GRBs for which an $R^{-2}$ density profile is excluded, and ultimately understand if the paradigm for the environment of GRBs should be drastically modified. Radio observations may be of great help, since the flux temporal behavior does not depend on p and is quite different in the case of constant or wind-like density profile. Similarly, the detection of SSC radiation can help solve this ambiguity.

3.3. Small Values of ${ϵ}_B$

For a long time, the typical value of ${ϵ}_B$ has been considered to variate between 0.01 and 0.1, both on the basis of theoretical considerations on particle acceleration and findings by numerical simulations. Indeed, the present understanding of the micro-physics at weakly magnetized shocks invoke the existence of self-generated micro-turbulence both behind and in front of the shock, at a level corresponding to ${ϵ}_B\sim 0.01–0.1$. This layer of intense micro-turbulence is expected on theoretical grounds and recently corroborated by numerical PIC simulations. Inferences of the value of ${ϵ}_B$ from the early modeling of afterglow radiation were broadly consistent with these numbers, confirming the presence of large self-generated fields in ultra-relativistic weakly magnetized shocks. More recently, several independent methods have provided evidence for significantly lower values.

In particular, several studies on GRBs with GeV temporally extended emission detected by LAT arrive to the same conclusions that in order to explain GeV radiation as part of the synchrotron emission, multi-wavelength observations require ${ϵ}_B={10}^{-6}$–${10}^{-3}$ [24,109,110,111,112]. Similar values have been inferred from studies that are based on radio, optical and/or X-ray emission and do not make use of high energy emission, such as [67,68,69,113]. A smaller magnetic field in the region where most of the particle cooling occurs might increase the expected relevance of the SSC component, as supported by recent detections of TeV radiation by IACTs.

Such small values of ${ϵ}_B$ may appear to be problematic [114], because strongly self-generated micro-turbulence must be present to ensure the scattering and acceleration of the particles, which otherwise would be simply advected away.

It was later pointed out that the inferred low values of magnetization might be indicative of a turbulence that is decaying on time-scales comparable with the electron cooling time [65,66]. From a theoretical perspective, indeed, the micro-turbulence is expected to decay beyond some hundreds of skin depths. This picture has been validated by PIC simulations, which however are still far from probing time-scales comparable with the dynamical time-scale of the system. Dedicated simulations the magnetic field does decay behind the shock, on a time-scale much longer than $c{\omega }_{pi}$. Immediately behind the shock, the magnetic field carries a magnetization ${ϵ}_B\sim 0.01$, which decays in time after ${10}^2$–${10}^3$ plasma times. Eventually, the magnetic field will settle to the shock-compressed value $4\Gamma B_u$, where $B_u$ is the magnetization of the upstream unperturbed medium. In this scenario, high-energy particles, which produce MeV-GeV photons, feel only the region close to the shock, where the magnetization is large, due to their short cooling time. Particles that cool on longer time-scales (and produce radio, optical and X-ray photons) cool on longer time-scales, and then in a region where the magnetic field has decayed.

The application of cooling in a decaying magnetic turbulence to four GRBs detected by LAT has proved to be very successful and even able to give indications on how fast the turbulence decays, being consistent with a power-law decay ${ϵ}_B\propto t^{-{\alpha }_t}$ with ${\alpha }_t\sim 0.5$ [66].

To understand and constrain the value of the magnetic field relevant for the particle cooling is of great importance, since an incorrect assumption or prior affects the estimates of all the other afterglow parameters, and in particular the density of the external medium [67,106].

A low value of ${ϵ}_B$ tends to increase the level of SSC luminosity for a given synchrotron luminosity. The recent detection of bright TeV emission from the afterglow of GRBs is an indication that this might indeed be the case. Existing and future TeV observations will shed a light on this issue, fostering a revision of our prejudice on the value of the magnetic field in the region where particles cool.

3.4. Variation of Micro-Physical Parameters with Time

Thanks to the increasing number of available observations on a wide range of frequencies (from radio to TeV) and times (from seconds to weeks), the basic assumption that micro-physical parameters (such as ${ϵ}_e,{ϵ}_e,p$ and ${\xi }_e$) are constant over the whole afterglow evolution can be tested. We comment on the hints (inferred from afterglow modeling) for temporal evolution of these parameters.

In case of well-sampled multi-wavelength light-curves, the modeling with synchrotron spectra is able not only to identify the location of the spectral breaks at a certain time but also evaluate their evolution in time. As a result, hints that micro-physical parameters ${ϵ}_e$ and ${ϵ}_B$ may vary with time have been found in some events with well-detailed multi-wavelength follow-up campaigns.

In [115], broad-band (from near infrared up to X-ray) afterglow data from GRB 091127 were interpreted in the standard external forward shock scenario assuming a constant-density medium. The good quality of the data allows one to identify the breaks in the light-curves and associate them with the synchrotron spectral breaks. As a result, the time evolution of the synchrotron breaks was estimated. In particular, it was calculated that the cooling break frequency ${\nu }_{cool}$ evolves as ${\nu }_{cool}\propto t^{-1.2}$, which is in contrast with synchrotron predictions for which a less steeper decay ${\nu }_{cool}\propto t^{-0.5}$ is expected. As a result, some assumptions of the standard model must be relaxed to remove the tension between observations and theoretical predictions. The most likely option able to explain the cooling break observational behavior without affecting the general interpretation of the data is to let the ${ϵ}_B$ parameter variate with time. Assuming that ${ϵ}_B\propto t^{0.49}$ the time evolution of ${\nu }_{cool}$ can be explained successfully.

In [116] for GRB 130427A modeling, in order to explain the observed fast evolution of the injection frequency ${\nu }_m\propto t^{-1.9}$, a temporal evolution of ${ϵ}_e$ is claimed. Considering that ${\nu }_m\propto {ϵ}_e^2$, a modest evolution of ${ϵ}_e$ following the trend ${ϵ}_e\propto t^{-0.2}$ is able to satisfactorily describe the observed light-curves.

A time-dependent evolution of the micro-physical parameters has also been proposed in order to explain the features observed in the early afterglow phase, which are not predicted by the external forward shock scenario such as X-ray afterglow plateaus, chromatic breaks, and afterglow rebrightenings [117,118,119,120].

Information from TeV observations can be certainly exploited in order to reduce the uncertainty on the values of the micro-physical parameters. The expansion of the broad band afterglow observations to a new spectral window will be a further test and a challenge for the multi-wavelength modeling based on the standard external forward shock scenario. In particular, the time evolution of the different energetic components, also including TeV emission, will give new insights useful to investigating the evolution of the micro-physical parameters. A first proof is provided by the well-sampled multi-wavelength emission observed for GRB 190114C, one of the few GRBs detected so far at TeV energies. The broadband emission can be explained only by invoking the evolution of the micro-physical parameters with time [121], as will be discussed in the next section.

3.5. Maximum Synchrotron Photon Energy

One of the expectations from Fermi-LAT observations of GRB afterglows was the identification of a spectral cutoff in the afterglow synchrotron spectrum marking the maximum energy of synchrotron photons [122,123]. Such a cutoff has not been firmly identified. Its location is directly connected with the shock micro-physics conditions and the maximum energy at which electrons can be accelerated. This maximum energy is typically estimated equating the timescale for synchrotron cooling and the acceleration timescale, where acceleration is assumed to proceed at the Bohm level, considered as the fastest rate. This estimate returns, hence the maximum energy of the accelerated particles. Assuming that the accelerated particles are electrons:

(53)$t_L^{\prime }=\frac{r_L}{c}=\frac{\gamma m_{e}c}{e{B}^{\prime }}$

where $r_L$ is the Larmor radius. For each crossing, the electrons gain energy by a factor $\sim 2$. On the other hand, the energy losses for synchrotron radiation on this timescale are:

(54)$\delta E^{\prime }=t_L^{\prime }{P}^{\prime }=\frac{\gamma m_{e}c}{e{B}^{\prime }}\frac{{\sigma }_{T}c{\gamma }^2{B}^{\prime 2}}{6\pi }=\frac{1}{6\pi }\frac{{\sigma }_T{m}_e{c}^2{B}^{\prime }{\gamma }^3}{e}$

The particle stops to gain energy when:

(55)$\delta E^{\prime }=\gamma m_e{c}^2$

Therefore, the maximum Lorentz factor for electrons can be derived:

(56)${\gamma }_{max}=\sqrt{\frac{3\pi e}{{\sigma }_T{B}^{\prime }}}$

The corresponding maximum synchrotron photon energy is:

(57)$h{\nu }_{max}^{\prime }=\frac{e{B}^{\prime }{\gamma }_{max}^{2}h}{2\pi m_{e}c}$

which for electrons is ∼50 MeV in their rest frame.

Similar considerations can also be conducted for protons. Following the same arguments presented below, one obtains:

(58)${\gamma }_{cool,p}=\frac{6\pi m_p^{3}c}{{\sigma }_T{m}_e^2{B}^{\prime 2}{t}^{\prime }}$

for the cooling Lorentz factor and:

(59)${\gamma }_{max,p}=\sqrt{\frac{3\pi e{m}_p^2}{{\sigma }_T{m}_e^2{B}^{\prime }}}$

for the maximum Lorentz factor, which sets a maximum photon energy of ∼100 GeV. Synchrotron emission is less efficient for protons so they are less affected by cooling and they can reach higher maximum Lorentz factors than the electrons.

Within this framework, it is expected that observations in the GeV band can be exploited to identify the presence of a cut-off in the HE tail domain. At the current stage, Fermi-LAT observations indicate that the afterglow component of the HE energy emission is usually modeled with a single power-law component with index ∼$-2$ and with no evidence of spectral evolution in time [124] and HE cut-offs.

The absence of the cut-off in the observational data may be explained in several ways. The most likely interpretations are the limited sensitivity of the LAT instrument in the GeV range and the possible contamination due to the rising of the SSC spectral component. As a result, the synchrotron cut-off is hidden behind the VHE spectral component, which can be detected in the GeV-TeV domain. This implies that TeV observations are fundamental in order to disentangle between the two spectral components and infer the cutoff of the synchrotron spectrum.

Another possible interpretation is that the lack of a cutoff in the observational data is genuine. In this case, the synchrotron emission can exceed the limit assumed for the maximum photon energy and extend in the GeV-TeV domain. This interpretation can be tested with VHE observations. The extension of the HE power-law derived by LAT up to the TeV domain should be consistent with VHE data, and no spectral hardening in the GeV-TeV band should be observed. Such a scenario is in contrast with the standard particle acceleration model presented in Section 2. A deep revision of our understanding of acceleration mechanisms is required in order to make TeV emission from synchrotron radiation possible. In particular, a mechanism, able to accelerate electrons up to PeV energies, is needed.

The calculation performed so far assumes the presence of a uniform magnetic field $B^{\prime }$ throughout the shock-heated plasma. If this assumption is rejected, it is possible to consider a non-uniform magnetic field, stronger close to the shock front and decaying downstream. Following the calculation of [125]. the magnetic field $B^{\prime }$ can be expressed in terms of the distance from the shock front x as:

(60)$B^{\prime }\left(x\right)=B_s^{\prime }{\left(\frac{x}{L_p}\right)}^{-\eta }+B_w^{\prime },$

where $B_s^{\prime }$ and $B_w^{\prime }$ are, respectively, the strongest and the weakest magnetic field strengths, $\eta$ is the power-law decaying index, and $L_p$ is the field decay length scale, which is estimated as [126]:

(61)$L_p={\left[\frac{m_p{\Gamma }_s{c}^2}{4\pi n{e}^2}\right]}^{1/2}=2.2\times {10}^7\frac{{\Gamma }_s}{n}\mathrm{cm}$

where ${\Gamma }_s$ is the shock front Lorentz factor and n is the number density of the accelerated particles in the shocked fluid comoving frame. As a consequence, the Larmor radius $r_L$ increases with the distance from the shock front, since $B^{\prime }\left(x\right)$ becomes weaker and an electron travelling downstream will be likely sent back upstream when $r_L\le x$. When considering the case $B_s^{\prime }\gg B_w^{\prime }$ and $x\gg L_p$ the particles will lose most of their energy in the region of low magnetic field. Therefore, from the condition that losses in the low magnetic field region should be greater than losses in the high magnetic field region, after some algebra the following condition is obtained:

(62)${\left(\frac{B_s^{\prime }}{B_w^{\prime }}\right)}^2\lesssim \frac{r_l}{L_p}$

valid for $\eta >1/2$ and $x_0/L_p\gg 1$ where $x_0$ is the width of the high magnetic field region. Considering that $x_0/L_p\equiv {(B_s^{\prime }/B_w^{\prime })}^{1/\eta }$, Equation (62) states that the Larmor radius in the high magnetic field region is larger than the actual width of the region and electrons will be barely deflected in such portion of the shocked plasma. As a result, it is possible to calculate the maximum Lorentz factor for electrons that lose most of their energy in the weak magnetic field region following the same conditions presented for the uniform magnetic field case:

(63)${\gamma }_{max}=\sqrt{\frac{3\pi e}{{\sigma }_T{B}_w^{\prime }}}$

As a result, the maximum synchrotron photon energy is given by:

(64)$h{\nu }_{max}^{\prime }=\frac{e{\gamma }_{max}^{2}h}{2\pi m_{e}c}\left(\frac{B_s^{\prime }}{B_w^{\prime }}\right)$

which is greater than the one calculated in Equation (57) by a factor $B_s^{\prime }/B_w^{\prime }$. Numerical calculations [41] show that this ratio can be larger than ∼${10}^2$. As a result, photons of energies ≳ 100 GeV can be produced via synchrotron process when assuming a non-uniform magnetic field, which decays downstream of the shock front and with particles losing most of their energies in the weakest field region.

In both interpretations presented here, TeV observations are fundamental in order to investigate with unprecedented details the possible presence or absence of the synchrotron cutoff spectrum. This have also a direct impact on the study of the possible radiation mechanisms responsible for the VHE component in GRBs.

3.6. Prompt Emission Efficiency

The overall efficiency ${\eta }_{\gamma }$ of the prompt emission mechanism is the result of three processes: the efficiency ${\eta }_{diss}$ of the (still unidentified) mechanism responsible for the dissipation of the jet energy, the efficiency ${ϵ}_e$ of the acceleration mechanism in converting the dissipated energy into random energy of the electrons, and the radiative efficiency ${ϵ}_{rad}$ of the electrons: ${\eta }_{\gamma }={\eta }_{diss}{ϵ}_e{ϵ}_{rad}$. Provided that it is reasonable to assume a fast cooling regime for the prompt emission (${ϵ}_{rad}=1$), the overall prompt efficiency is limited by the capability of the dissipation mechanism in extracting the kinetic or magnetic energy of the jet and the capability of the particle acceleration process to convey a fraction of this energy into the non-thermal electron population. The value of the efficiency then provides a fundamental clue to placing constraints on the origin of energy dissipation in GRBs, which is still quite uncertain, discriminating between matter and magnetic dominated jets.

From the definition of ${\eta }_{\gamma }=E_{\gamma }/E_0$ (where $E_0=E_{\gamma }+E_k$ is the initial explosion energy), we can write the relation $E_k=(1-{\eta }_{\gamma })/{\eta }_{\gamma }{E}_{\gamma }$. The parameter ${\eta }_{\gamma }$ can then be estimated from the comparison of the energy $E_{\gamma }$ emitted in the prompt phase and the energy $E_k$ left in the jet after the end of the prompt emission (i.e., at the beginning of the afterglow phase). While the former is directly estimated from observations, the latter can be inferred only indirectly, from the modeling of afterglow radiation.

One of the most adopted methods to infer $E_k$ for large samples of GRBs is to rely on the X-ray luminosity and use it as a proxy for the energy content of the blast-wave [127,128,129,130,131]. This method is solid as long as the X-ray band lies above $max({\nu }_m,n{u}_c)$ and is not affected by inverse Compton cooling. If these two conditions are verified, then the electrons emitting X-ray photons are in a fast cooling regime and their cooling is dominated by synchrotron losses. The luminosity produced is then proportional to the energy content of the accelerated electrons $E_k{ϵ}_e$. Assuming a value (typically 0.1) for ${ϵ}_e$, then $E_k$ can be estimated. Investigations based on the X-ray emission have inferred very large values of ${\eta }_{\gamma }$, between 0.5 and 0.9 [49,129,132,133].

The very same approach can also be applied to 100 MeV-GeV photons detected by the LAT, under the assumption that these are synchrotron photons. A strong correlation between the GeV luminosity and $E_{\gamma ,iso}$ has been indeed found, supporting the possibility that GeV photons lie in the high-energy part of the synchrotron spectrum, where the afterglow luminosity is proportional to $E_k{ϵ}_e$ and can be used, similarly to the X-ray luminosity, to estimate $E_k$ [134]. A study by [70] revealed that the energetics $E_k$ inferred independently from X-ray and GeV luminosities on a sample of 10 GRBs are inconsistent with each other. The authors show that the inconsistency is solved if ${\nu }_c>{\nu }_X$ (where ${\nu }_X$ is the X-ray frequency), or if Compton losses are important in the X-ray band. Full modeling of the GeV, X-ray and optical data support this scenario. In both cases, ${ϵ}_B$ is required to be much smaller than usually assumed, with values in the range ${10}^{-7}$–${10}^{-3}$. This analysis shows that the GeV band is a much better proxy for $E_k$, since it is above ${\nu }_c$ and is not affected by inverse Compton cooling, due to Klein–Nishina suppression. Adopting GeV luminosities as a proxy for $E_k$, the estimated values of $E_k$ increase, also affecting the estimates of ${\eta }_{\gamma }$, which are around 5–10% [24].

A correct estimate of ${\eta }_{\gamma }$ is extremely important, since its value is related to the mechanism dissipating energy in the jet. Since internal shocks can hardly reach values of ${\eta }_{\gamma }$ larger than 10%, values around 50–90% have been used to argue that internal shocks are not a viable mechanism to explain prompt emission in GRBs, and more efficient mechanisms should be considered (e.g., magnetic reconnection). If the efficiency is, however, smaller than initially estimated, internal shocks may still be a viable solution. Moreover, different estimates of ${\eta }_{\gamma }$ lead to different estimates on the total initial jet energy $E_0=E_{\gamma ,iso}+E_k$, with repercussions on the energy budget of GRBs and finally on their progenitors and mechanisms for jet launching. Small values of ${ϵ}_B$ may then relax the problem with very large prompt efficiency, which is definitely unreasonable for internal shocks, but also difficult to attain for magnetic reconnection models (for a discussion, see e.g., [19,20,135]).

A scenario where the magnetic field strength is relatively low in the emitting region implies a stronger SSC emission. Recent TeV detection of GRBs support this scenario, and provide additional observations to constrain the magnetic field. Moreover, as shown by the first detections by IACTs (Section 4), the energy in the TeV component is comparable to the energy in X-rays, providing better estimates for the energy budget in the afterglow phase. Future detections from a larger sample of GRBs can help in assessing the energy budget of the jet more precisely during the afterglow emission, add important information to constrain the efficiency of the prompt emission and favor or exclude some dissipation scenarios.

3.7. Fraction of Accelerated Particles

As described in Section 2, the representation of the relativistic shock acceleration in GRB afterglow relies on some free parameters. These values describe the energy equipartition between particles and magnetic field and the non-thermal accelerated particle distribution.

In particular, the parameter ${\xi }_e$ is responsible to account for the fraction of particles (here we consider electrons, but the same considerations are also valid for protons) accelerated into a non-thermal distribution. This means that from relativistic shock theory it is expected that a fraction $1-{\xi }_e$ of electrons is heated into a thermal distribution rather than a non-thermal one (see Figure 9). For simplicity, it is usually assumed in afterglow studies that ${\xi }_e=1$, i.e., all the particles are accelerated into a non-thermal distribution. Such an assumption is used in order to avoid the large degeneracy which affects the GRB parameters when including this additional free value. In particular, afterglow modeling predictions obtained assuming ${\xi }_e=1$ for parameters $E_k$, n, ${ϵ}_e$ and ${ϵ}_B$ cannot be distinguished from those estimated for any ${\xi }_e$ in the range $m_e/m_p<{\xi }_e<1$ and afterglow parameters $E_k^{\prime }=E_k/{\xi }_e$, $n^{\prime }=n/{\xi }_e$, ${ϵ}_e^{\prime }={\xi }_e{ϵ}_e$ and ${ϵ}_B^{\prime }={\xi }_e{ϵ}_B$ [73]. This can be proven considering the dependencies of the jet dynamics and shock energy equipartition on the model parameters. As shown in BM76 and by previous calculations on the evolution of a relativistic blastwave, in the self-similar regime the bulk Lorentz factor evolve as $\Gamma \propto {(E_k/n)}^{1/2}$. As a result, the same flow evolution can be obtained with different values of $E_k$ and n as long as their ratio is preserved. The fraction of energy given to the magnetic field is reduced by a factor ${\xi }_e$ but at the same time the energy density given to the shock is increased by a factor $1/{\xi }_e$ so that the magnetic field energy density is the same in the scenario when including ${\xi }_e$ or not. Analogous considerations can be also conducted for the number and the energy density of the electrons so that their values are preserved. As a result, in principle it is not possible to distinguish between the two parameter sets obtained for ${\xi }_e$ or for any value $m_e/m_p<{\xi }_e<1$. This implies that afterglow model parameters, when considering ${\xi }_e\ne 1$, are estimated with an uncertainty of factor of $m_e/m_p$, and the afterglow observations do not directly constrain their values (e.g., $E_k$ or ${ϵ}_B$) but rather are a fraction of their value multiplied or divided by ${\xi }_e$ (e.g., $E_k/{\xi }_e$ or ${ϵ}_B{\xi }_e$).

It is possible to obtain information on the value of ${\xi }_e$ through PIC simulations or indirect features of the thermal component on the synchrotron afterglow spectra. As mentioned in the previous section, PIC simulations of weakly magnetized shocks have found that the downstream population include around ∼3% of the electrons, which are accelerated into a non-thermal distribution. In case the efficiency is low (around ∼10% or less) the presence of a large population of thermal electrons may affect the afterglow radiation spectra. The thermal electrons are distributed at lower energies than the non-thermal ones since $\eta \gamma m_e{c}^2$≪$\gamma m_p{c}^2$, where $\eta \ll m_p/m_e$ is a factor describing the ignorance on the plasma physics governing electron heating beyond $\gamma m_e{c}^2$. As a result, the synchrotron radio component emission may be affected through the production of a new emission component from thermal electrons (for $\eta \gg 1$ and moderate $1/{\xi }_e$) or a large self-absorption optical depth (for ${\xi }_e\ll 1$), which may be visible in a time scale of few hours. Possible effects of the thermal component are discussed in [136,137,138].

Information from the TeV component cannot completely solve the degeneracy between afterglow parameters and cannot provide additional clues on the non-thermal electron distribution. However, it can provide unique fundamental data useful to constrain the other afterglow parameters less constrained such as ${ϵ}_B$ and the density. This will also impact the ${\xi }_e$ calculation since it can help to reduce the degeneracy between the sets of solutions available in the parameter space. Indeed, the multi-wavelength modeling of GRB 190829A (detected at TeV energies by H.E.S.S.) showed that the only way to explain all the detected radiation, from radio to TeV, is to introduce the parameter ${\xi }_e$ in the modeling, which is constrained by data to be ${\xi }_e<0.13$ [72].

4. Discovery of a TeV Emission Component in GRBs

The robust theoretical framework developed throughout the years to explain the afterglow radiation predicts that, to some extent, GRBs should be TeV emitters (Section 2). Observations in the HE band, and in particular the presence of ≳GeV photons with energies up to ∼100 GeV, support this possibility. On the other hand, from the observational side, the search for such emission is hampered by several drawbacks. Space-born telescopes, such as Fermi-LAT, sensitive up to few hundred GeV, have an hard time with GRBs due to their low $\gamma$-ray photon flux at the highest energies (∼${10}^2$ GeV), caused by their cosmological distance and strong EBL absorption. These difficulties can be overcome by the much larger effective area of IACTs in the common energy range of sensitivity (50–300 GeV). As a downside, IACTs have a small field of view (a few degrees wide), higher low-energy threshold (≳50–100 GeV), and reduced duty cycle (less than 10%).

In the last decades, IACTs have performed a huge effort to become instruments suitable for GRB observations. In particular, the efforts have been focused in two directions: (i) the development of fast repointing systems to promptly react to GRB alerts and start observations with delays of a few tens of seconds after the trigger time, and (ii) the extension of the energy threshold below 100 GeV, important to reduce the impact of the EBL attenuation on the detection probability of cosmological GRBs.

After a decade of VHE observations resulting in non-detections, the first announcement of GRBs detected by IACTs arrived in 2019, thanks to the MAGIC and H.E.S.S. telescopes [139]. These detections have firmly established that GRBs can be bright sources of TeV radiation. Somewhat unexpectedly, VHE emission was also detected several hours/a few days after the GRB onset, and up to energies of ≳3 TeV. The timescales of the detections place the origin of the emission in the afterglow phase. The TeV emission has been studied and interpreted in a multi-wavelength context, in order to evaluate the properties and the nature of the responsible radiation mechanisms. In particular, investigations have focused on SSC, external inverse Compton (EIC), and synchrotron radiation.

In this section, all GRBs for which a detection (significance $>5\sigma$) or a hint of detection (significance between 3 and 5$\sigma$) has been claimed by Cherenkov telescopes, are presented. These are, in total, six events (one short and five long): GRB 160821B (Section 4.1), GRB 180720B (Section 4.2), GRB 190114C (Section 4.3), GRB 190829A (Section 4.4), GRB 201015A (Section 4.5) and GRB 201216C (Section 4.6). For each event, we start with a brief description of the prompt and afterglow multi-wavelength observations. Then, we describe VHE observations and summarize the main results. These detections being a novelty, and some of them laying close to the sensitivity detection threshold of the instrument, we describe in detail the VHE data analysis, the calculation of the significance excess at the GRB position (following the usual prescription used for VHE sources presented in [140]), and the methods adopted for the derivation of the spectral energy distribution (SED) and of the light-curves. For each GRB, we also present the interpretations that have been put forward in the literature. A discussion on the main common properties and differences among this initial population of VHE GRBs and with respect to the whole GRB population is presented in Section 5, where we also address the question of what we have been learning from these few detections.

In this section, all quoted times refer to the time elapsed from the trigger time $T_0$ of the Swift-BAT or Fermi-GBM instrument, as will be specified. Photon indices are given in the notation $N_{\nu }\propto {\nu }^{\alpha }$, while temporal indices are defined by $F\left(t\right)\propto t_T^{\beta }$.

4.1. GRB 160821B

GRB 160821B is a short GRB at $z=0.162$ detected by the Swift-BAT [141] on 21 August 2016 at $T_0=22:29:13$ UT and by the Fermi-GBM [142]. The analysis of MAGIC observations shows a $\sim 3\sigma$ excess at the GRB position.

4.1.1. General Properties and Multi-Wavelength Observations

The BAT prompt spectrum ($T_{90}=0.48$ s) is well described by a power-law with index $\alpha =-0.11\pm 0.88$ and an exponential high-energy cutoff, corresponding to a peak energy with $E_p=46.3\pm 6.4$ keV [143]. The GBM prompt spectrum ($T_{90}=1.088\pm 0.977$ s) is fitted with a cutoff power-law function as well, with $E_p=92\pm 28$ keV. Being located at redshift $z=0.162$, this is one of the nearest short GRBs detected up to date. Its isotropic-equivalent energy $E_{\gamma ,iso}\sim 1.2\times {10}^{49}$ erg is toward the low energy edge of the known distribution for short GRBs [144]).

Afterglow observations are available in the radio, optical, X-ray and (V)HE band. Fermi-LAT observations were performed from the trigger time up to 2315 s and from 5285 s to 8050 s, and did not reveal any significant excess in the 0.3–3 GeV band [145]. Swift-XRT observations [146] started 57 s after the trigger time and revealed a complex behavior of the X-ray afterglow light curve. An initial plateau is followed by a steep decay at around ${10}^2$ s. Then, a power-law decay with index ∼$-0.8$ is observed after ∼${10}^3$ s [147,148]. Optical observations were performed by several instruments [149,150,151,152], revealing the presence of a fading source with a magnitude $r=22.6\pm 0.1$ mag $0.95$ h after $T_0$. The identification of the host galaxy allowed the measurement of the spectroscopic redshift $z=0.162$. The GRB is located in the outskirts of the host spiral galaxy, at ∼15 kpc projected distance from its center [148,153]. A fading radio afterglow counterpart was observed at 6 GHz by VLA starting from 3.6 h after the burst trigger [154].

4.1.2. VHE Observations and Results

The MAGIC telescopes started the follow-up of GRB 160821B with a very short delay of 24 s after $T_0$ and continued observations for about 4 h [145]. The observations were performed with a relatively high Night Sky Background (NSB) (2–8 times brighter than in dark nights) due to the presence of the Moon, and in mid-high zenithal angle conditions (from 34${}^{\circ }$ to ${55}^{\circ }$). Unfortunately, the first ∼1.7 h of the data were strongly affected by clouds. As a result, dedicated and optimized software configurations were used. A more stringent image cleaning with respect to dark conditions was applied to take into account the spurious contribution coming from the high NSB. The analysis required cuts optimized on the Crab Nebula and on Mrk421 observed in similar conditions, and correction factors for the low atmospheric transmission, calculated thanks to the LIDAR facility [155].

The pre-trial analysis showed the presence of a $3.1\sigma$ ($2.9\sigma$ post-trial⁴) significant excess at the GRB position provided by Swift-XRT (see Figure 10). The flux has been estimated for energies above $0.5$ TeV assuming a power-law spectrum with photon index $\alpha =-2$. In the first 1.7 h, where data taking was affected by bad atmospheric transmission, only flux upper limits could be derived. This time window has been divided into two time intervals (24–1216 s and 1258–6098 s) and the derived upper limits are, respectively, $1.1\times {10}^{-11}$ cm${}^{-2}$s${}^{-1}$ and $5.4\times {10}^{-12}$ cm${}^{-2}$ s${}^{-1}$. In the subsequent 2.2 h (6134–14,130 s) assuming that the signal is real, a flux value could be calculated and is $9.9\pm 4.8\times {10}^{-13}$ cm${}^{-2}$ s${}^{-1}$. For the same time interval, a flux upper limit has also been estimated and gives $3.0\times {10}^{-12}$ cm${}^{-2}$ s${}^{-1}$. All the mentioned fluxes are shown in Figure 11 (red symbols) and refer to the observed values, i.e., without correcting for EBL absorption. Upper limits have been calculated at a 95% confidence level following the prescriptions of [156]. The low (∼3$\sigma$) significance estimated did not allow us to obtain an unfolded spectrum. As a result, in the SED (Figure 12) the reconstructed flux in the third bin (6134–14,130 s) over the energy range 0.5–5 TeV is represented as an error box. The statistical error on the photon flux has been taken into account, while, for simplicity, the systematic error for the assumed spectral index was neglected. The flux inferred by MAGIC observations in the 0.5–5 TeV energy range would imply a TeV luminosity at least 5 times larger (when de-absorbed by EBL) than the luminosity emitted in the X-ray band.

4.1.3. Interpretation

A modeling of the multi-wavelength observations, including MAGIC data, has been presented by the MAGIC Collaboration in [145]. The emission is interpreted as the sum of several components, dominating at different times and in different energy bands:

• synchrotron and SSC emission from electrons accelerated at the forward shock; this is in general the dominant emission component;

• synchrotron emission from electrons accelerated by the reverse shock, which is found to dominate the radio band until $t\sim 4.8$ h;

• kilonova emission, powered by freshly synthesized r-process elements released in neutron star mergers; this component is found to dominate the optical/nIR from around 1 to 4 days [148,153];

• an X-ray extended emission component, widely attributed to long-lasting activity of the central engine, here dominating the X-ray band for $t<{10}^3$ s.

In performing this multi-component modeling, the synchrotron and SSC forward shock emission have been calculated with a one-zone numerical code (see [71] for details), while the reverse shock and kilonova emission contributions have been taken from [148]. Only X-ray data at $t>{10}^3$ s have been included in the modeling, to exclude the extended emission component. The broad-band modeling is shown over-plotted to the light-curves in Figure 11 (solid lines) and to the SED between 1.7 and 4 h in Figure 12.

A very good agreement between data and modeling is found in radio (green lines and points), optical (yellow and pink lines and points) and X-rays (blue lines and points). A large degeneracy is present in the parameters, and the data modeling only allows us to identify the ranges for the permitted values of each parameter. These are reported in Table 1 and we note that they are very similar to those estimated in [148] and in a later work by [157]. In the allowed parameter space defined by radio, optical and X-ray observations, different combinations of the parameters predicts different SSC fluxes at 1 TeV are found, reaching at most $F_{SSC}^{\left(1TeV\right)}\sim 2\times {10}^{-13}$ erg cm${}^{-2}$ s${}^{-1}$. This value, when attenuated by EBL, is at least one order of magnitude fainter than the one inferred from the data analysis of the MAGIC observations. In other words, the parameter space constrained by the observations at lower frequencies is unable to account for such energetic TeV emission, and the SSC forward shock scenario fails to reproduce the observations, provided that the hint of excess found by MAGIC is a real signal from the source.

An alternative scenario that has been explored is the external inverse Compton (EIC) scenario, investigated by [157]. These authors first consider a one-zone SSC model, and reach similar conclusions to those presented by the MAGIC Collaboration [145]: the SSC mechanism predicts a TeV flux around 1–2 orders of magnitude lower than the MAGIC observations (see Figure 13, orange curves). The alternative EIC scenario is then considered by the authors, where the seed photons are provided by the extended X-ray emission and the X-ray plateau. The extended emission and the plateau are fitted using two phenomenological functions. The energy spectrum of the late-prompt emission is described by a broken power-law (see Figure 13, top and bottom panels). For the EIC model, the VHE spectrum is inferred for three different observed times ($t=1.1,1.8,2$ h) and compared to the MAGIC flux averaged between 1.7 and 4 h. As can be observed in Figure 13 (bottom panel), the model flux at 2 h under-predicts the MAGIC flux (the MAGIC observed flux, green shaded area, should be compared with the EBL-absorbed model flux). We conclude that the EIC model is also unable to explain the large TeV flux suggested by MAGIC observations.

4.2. GRB 180720B

GRB 180720B is a long GRB at $z=0.654$ triggered on 20 July 2018 by the Fermi-GBM [158] at $T_0=14:21:39$ UT and 5 s later by the Swift-BAT instrument [159]. The H.E.S.S. telescopes observed and detected GRB 180720B about 11 h after the prompt emission, with a ∼5$\sigma$ statistical significance.

4.2.1. General Properties and Multi-Wavelength Observations

The extremely bright prompt emission of this event, which is the seventh in brightness among the GRBs detected by the Fermi-GBM until then, lasted for $T_{90}=48.9\pm 0.4$ s and released an isotropic-equivalent energy $E_{\gamma ,iso}=6.0\pm 0.1\times {10}^{53}$ erg in the 50–300 keV range.

Multi-wavelength afterglow observations covered the entire electromagnetic spectrum (see Figure 11). A significant signal was detected by Fermi-LAT from the trigger time up to 700 s, with the highest photon energy of 5 GeV detected 137 s after the burst trigger [160]. The Swift-XRT telescope observed and identified a bright afterglow starting from 90 s. This was still visible almost 30 days after the trigger time. The late-time light curve (from $2\times {10}^3$ s to $4\times {10}^6$ s) can be modeled with an initial power-law decay with an index $-1.{19}_{-0.02}^{+0.01}$ followed by a break at $t_{break}=8\times {10}^4$ s to an index of $-1.{55}_{-0.05}^{+0.04}$⁵. Several optical observations [161,162,163,164,165,166,167,168,169,170] revealed the presence of a counterpart and allowed us to estimate the redshift value of $z=0.654$. The optical afterglow was observed to be slowly fading at an almost constant rate from around 10–11 h after the trigger time [171,172] as discussed in [173]. Radio observations (not shown in the figure) were also performed starting from ∼1.7 days after the burst showing a steep power-law decaying emission [173].

4.2.2. VHE Observations and Results

The H.E.S.S. telescopes followed-up the event for ∼2 h starting from 10.1 h, revealing the presence of a source with a $5.3\sigma$ pre-trial significance ($5.0\sigma$ post-trial⁶). The observation was performed in standard dark and good weather conditions with a mean zenith angle of 31.5${}^{\circ }$. Another observation was performed under similar conditions 18 days after the previous one with results consistent with background events. The inferred flux at ∼11 h and the flux upper limit at 18 d are shown in Figure 14 (red symbols).

The observed spectrum in the 0.1–0.44 TeV energy band has been fitted both with a power-law (panel on the left in Figure 15) and with a power-law with a cutoff $F_{int}=F_{0,int}{\left(\frac{E}{E_0}\right)}^{-{\gamma }_{int}}$, to describe an intrinsic power-law spectrum affected by the EBL attenuation $e^{-\tau (E,z)}$ (see Figure 15, panel on the right). With reference to the second fit, the analysis returns a photon index ${\gamma }_{int}=1.6\pm 1.2$ (statistical) $\pm 0.4$ (systematic) and a flux normalization $F_{0,int}=(7.52\pm 2.03$ (statistical) ${}_{-3.84}^{+4.53})$ (systematic) $\times {10}^{-10}$ TeV${}^{-1}$ cm${}^{-2}$ s${}^{-1}$, evaluated at an energy $E_{0,int}=0.154$ TeV.

4.2.3. Interpretation

The H.E.S.S. Collaboration explored two possible radiation mechanisms to explain the VHE emission from GRB 180720B [4]: synchrotron emission and SSC radiation. A full modeling is not performed, and the discussion and comparison among the two different scenarios is based on estimates of the typical and maximal electron energy necessary in the two cases and on the comparison between spectral and temporal indices in different energy ranges. A synchrotron spectrum with a flat ($\alpha \sim -2$) slope extending from X-ray to VHE could model the emission with one single broad component and explain the similarity between the H.E.S.S., Fermi-LAT, and Swift-XRT luminosities, and the consistency among their photon index values. The large error on the VHE photon index, however, is not placing strong constraints, leaving open both the possibility of a consistency with the extrapolation of the synchrotron spectrum but also the possibility of a spectral hardening, indicative of a second component. A synchrotron origin of ${10}^2$ GeV photons would require to find a process able to accelerate electrons up to PeV energies, which is in excess of the maximum electron energy achievable in external shocks (for a discussion, see Section 2.3.1). Adopting the standard Bhom limit, >100 GeV emission 10 h after the burst would require a huge bulk Lorentz factor $\Gamma \sim 1000$, which at these late times is really unlikely. As a result, these strong requirements disfavor the synchrotron emission as responsible of the VHE component in GRB 180720B.

The SSC scattering, on the contrary, arises as a natural candidate. A full broad-band modeling of GRB 180720B data in this scenario is presented in [174]. A numerical code reproducing the synchrotron and SSC emission in the afterglow shocks has been used (see [175]). The resulting light-curves and SED in the H.E.S.S. observational time window are shown in Figure 16. The full emission is explained as afterglow forward shock radiation (except for the initial peak in the optical and X-ray curves at $t\sim {10}^2$ s, which is attributed to reverse shock emission). In the case of a constant-density ISM environment, the parameters that best reproduce the data are: $E_k={10}^{54}$ erg, $n=0.1$ cm${}^{-3}$, ${ϵ}_e=0.1$, ${ϵ}_B={10}^{-4}$, ${\Gamma }_0=300$ and $p=2.4$. As it can be noticed, the equipartition factor ${ϵ}_B$ needs to assume quite a low value in order to explain the observations. A stellar wind-like environment is discarded by the authors on the basis of the comparison between the expected flux at ∼1–10 GeV, following the prescriptions of [48] ($\gtrsim 2\times {10}^{-6}$ erg cm${}^{-2}$ s${}^{-1}$ at $t\approx 100$ s) and the one observed by Fermi-LAT (∼${10}^{-8}$ erg cm${}^{-2}$ s${}^{-1}$ at $t\approx 100$ s). A low magnetic field equipartition factor is derived from the condition $E_{KN}\gtrsim$ 0.44 TeV at $t\approx 10$ h, where $E_{KN}$ is the energy at which the KN scattering becomes relevant. A transition energy between the synchrotron and SSC component of ∼1 GeV is derived. Such value falls into the Fermi-LAT energy range and is compatible with a hardening of the spectrum in the VHE band. However, since the Fermi-LAT sensitivity is above the predicted flux of GRB 180720B at 10 h in the GeV band, the data cannot firmly confirm the presence of this transition.

4.3. GRB 190114C

GRB 190114C is a long GRB at redshift $z=0.42$ triggered by the Swift-BAT [176] on 14 January 2019 at $T_0=20:57:03$ UT, and by the Fermi-GBM [177]. The event was detected also by several other space $\gamma$-ray instruments such as AGILE, INTEGRAL/SPI-ACS, Insight/HXMT, and Konus-Wind [71]. MAGIC detected GRB 190114C starting ∼60 s after $T_0$ with a significance above $50\sigma$.

4.3.1. General Properties and Multi-Wavelength Observations

The duration of the prompt emission is $T_{90}\approx 116$ s as measured by Fermi-GBM and $T_{90}\approx 362$ s by Swift-BAT. However, the prompt light curve showed a multi-peak structure only for about 25 s, suggesting that the remaining activity, which is characterized by a smooth power-law decay and recorded by these instruments, may already be the afterglow emission. Support to such an interpretation is also obtained from a joint spectral and temporal analysis of the Fermi-GBM and Fermi-LAT data [178]. The total radiated prompt energy is $E_{\gamma ,iso}=(2.5\pm 0.1)\times {10}^{53}$ erg in the energy range 1–10${}^4$ keV [179].

Extensive follow-up observations from several different instruments from GeV to radio are available. Light-curves are shown in Figure 17. Fermi-LAT observations started since the beginning of the prompt phase. A GeV counterpart was detected from $T_0$ to 150 s, when the burst left the LAT field of view and remained outside it until 8600 s. When LAT resumed observations, a significant signal was still detected at a flux level ∼2 ×${10}^{-10}$ erg cm${}^{-2}$ s${}^{-1}$ (0.1–1 GeV). After ∼60 s from the burst trigger, Swift-XRT started follow-up observations, which covered in total ∼${10}^6$ s. The light-curve in the 1–10 keV energy band is consistent with a power-law decay $F\propto t^{\alpha }$ with $\alpha =-1.36\pm 0.02$ [71]. NuSTAR and XMM-Newton observations are also available around 1–2 days. The NIR, optical and UV data were taken from around ∼100 s. The early emission is particularly bright and is interpreted as dominated by the reverse shock component [180]. Afterwards, the decay rate flattens and then steepens again after ∼3 ×${10}^4$ s (see Figure 17). The Nordic Optical Telescope measured a redshift of $z=0.4245\pm 0.0005$ [181], which was then confirmed by Gran Telescopio Canarias [182]. Radio and sub-mm data were taken from ∼${10}^4$ s and exhibit an achromatic behavior, possibly dominated by the reverse shock in the sub-mm range, followed by emission at late times with nearly constant flux.

4.3.2. VHE Observations and Results

After receiving (at 22 s after the BAT trigger time) and validating (at 50 s) the GRB alert, the MAGIC telescopes started observing GRB 190114C at 57 s and operated stably from 62 s, starting from a zenith angle of 55.8${}^{\circ }$. Observations lasted until 15,912 s, when a zenith angle of 81.14${}^{\circ }$ was reached. The observation was performed in good weather conditions but in presence of the moon, resulting in a night sky background approximately six times higher than the standard dark night conditions. The results of the offline analysis demonstrate a clear detection above the $50\sigma$ level in the first 20 min of observation [3].

The light-curve (see Figure 18, upper panel) for the intrinsic flux (i.e., corrected for the EBL absorption) in the 0.3–1 TeV range was derived starting from 62 s and up to 2454 s. The TeV light curve is well described by a power-law with temporal decay index ${\beta }_T=-1.60\pm 0.07$, steeper than the one exhibited by the X-ray flux. The temporal evolution of the intrinsic spectral photon index ${\alpha }_{int}$ of the TeV differential photon spectrum is shown in the bottom panel. A constant value of ${\alpha }_{int}$$\approx -2$ is consistent with the data, considering the statistical and systematic errors, but there is evidence for a softening of the spectrum with time. The spectral fit in the 0.2–1 TeV energy range for the time-integrated emission (62–2454 s) returns ${\alpha }_{obs}=-5.34\pm 0.22$ and ${\alpha }_{i}nt=-2.{22}_{-0.25}^{+0.23}$ for the observed and EBL-corrected spectrum, respectively.

4.3.3. Interpretation

The properties of the VHE light curve and spectrum of GRB 190114C were studied by the MAGIC Collaboration in [3]. The PL behavior, the absence of variability, and the relatively long timescale of the emission support the evidence that the VHE component belongs to the afterglow phase. An estimate of the amount of radiated power in the TeV range can be derived, assuming that the afterglow onset is at ∼6 s [178]. In this case, the energy radiated in the TeV band is ∼10% of the isotropic-equivalent energy of the prompt emission $E_{\gamma ,iso}$ considering the temporal evolution estimated from the MAGIC light-curve.

The energy of the photons observed by the MAGIC telescopes was compared with the maximum energy of synchrotron photons assuming two possible scenarios for the radial profile of the external density, namely constant and wind-like (Figure 19). These estimates of the maximum energy are based on the widely adopted limit on the maximum electron Lorentz factor set by equating the acceleration at Bhom rate with the synchrotron cooling rate (see Equation (57)). Adopting this limit, synchrotron emission can not account for the TeV photons detected by MAGIC, and a different radiation mechanism must be invoked.

An additional, model independent indication for the presence of a spectral component other than synchrotron is evident after multi-wavelength simultaneous SEDs are built. In [71], the VHE data were rebinned into five time intervals and XRT, BAT, GBM and LAT data were added when available, i.e., in the first two time intervals (Figure 20). The spectrum demonstrates a double-peaked behavior with a first peak in the X-ray band and the second one in the VHE band. The Fermi-LAT data play a particularly important role in revealing the shape of the SED, as they show a dip in the flux, strongly supporting an interpretation of the whole SED as a superposition of two distinct components.

Following these considerations, in [71], the SSC mechanism is explored. The broad-band emission is modeled with a numerical code reproducing the synchrotron and SSC radiation in the external forward shock scenario, including the proper KN cross section and the effects of $\gamma -\gamma$ annihilation.

The predicted spectra and lightcurves are compared with the data in Figure 21 and Figure 22. Acceptable modeling of the multi-wavelength afterglow spectra have been found for a constant medium with $E_k\gtrsim 3\times {10}^{53}$ erg, ${\epsilon }_e\approx 0.05-0.15$, ${\epsilon }_B\approx (0.05-1)\times {10}^{-3}$, $n\approx 0.5-5$ cm${}^{-3}$ and $p\approx 2.4-2.6$. It is found that, the peak of the SSC component being below 200 GeV, the KN suppression and the $\gamma -\gamma$ internal absorption have a non-negligible role in shaping the peak of the VHE spectrum. The modeling reproduces the XRT, LAT and TeV emission very well (solid blue curve in Figure 21 and solid blue, green and red curves in Figure 22), while it overproduces both the optical and radio flux at late times (solid violet, yellow and cyan curves in Figure 22). According to [71], a similar fit is found, also assuming a wind-like profile for the external density. In this case, the parameters are $E_k=4\times {10}^{53}$ erg, ${\epsilon }_e=0.6$, ${\epsilon }_B=1\times {10}^{-4}$, $A_{*}=0.1$ and $p=2.4$. Very interestingly, the modeling shows that the late LAT observation (around ${10}^4$ s) is completely dominated by SSC emission (red-dashed curve in Figure 22). A different type of modelization is also investigated by [71], under the requirement to model optical data. In this case (dotted curves in Figure 22), the fit is very good for optical, X-ray and LAT observations, but fails in reproducing the MAGIC light-curve.

The values inferred for the GRB afterglow parameters are similar to those used for past GRB afterglow studies at lower frequencies. This is an indication that the SSC component can be a relatively common process for GRB afterglows, since it does not require peculiar values of the parameters to be explained.

Several other successful modelings of GRB 190114C data within the synchrotron and SSC external forward shock scenario have been published in the literature [174,184,185,186]. A summary of the parameters inferred by different works can be found in Table 2.

In [174], the X-ray, optical and LAT data before 100 s are attributed to reverse shock emission or prompt contribution. A constant-density environment for the circumburst material is assumed. A time-averaged SED (50–150 s) is estimated, demonstrating that at GeV energies a transition between the synchrotron and SSC component can be identified. From the re-analysis of LAT data, a hard photon index ($1.76\pm 0.21$) is derived, in agreement with the hardening of the spectrum caused by the rising of the SSC component. Differently from what is observed by [71,184] $\gamma -\gamma$, the absorption does not contribute significantly in shaping the VHE spectrum.

A similar interpretation is given in [185], although the inferred value of ${ϵ}_B$ is larger (${ϵ}_B\sim {10}^{-3}$). A consistent modeling of the multi-wavelength observations as synchrotron and SSC radiation is found in both the ISM and wind-like scenarios. The SED at 80 s (see Figure 2 and 5 in [185]) and the broad-band light curves (Figure 3 and 6 in [185]) are reproduced, despite at 10${}^3$ s, the model predictions in the 0.3–1 TeV band and X-rays are slightly brighter than the observed data.

In [186], analytical approximations are adopted for the description of the synchrotron and SSC components. In addition, the KN cut-off energy and the $\gamma -\gamma$ absorption contribution are calculated and compared with the data. A wind-like environment was used for the circumburst medium. The SEDs in the time intervals 68–110 s and 110–180 s are modeled, and the two values of the KN cut-off energies calculated at these times are ∼3.7 TeV and $2.1$ TeV. This implies that the KN effect is relevant only at TeV levels and the VHE data can be modeled assuming that the SSC scattering is in the Thompson regime. The $\gamma -\gamma$ absorption is also considered negligible since the estimated attenuation factor is way lower than the one due to EBL attenuation, and it reaches values around unity only for energies ≳ 1 TeV.

In [184], the multi-wavelength data were fitted with a single-zone numerical code with an exact calculation of KN cross-sections as well as the attenuation due to the pair production mechanism. Smoothed analytical approximations for the electron injection function were used. A systematic scan over a 4-dimensional parameter space was performed to search for the best-fit solution at early and later times. The SED calculated at 90 s and 150 s (Figure 3 and Figure 9, respectively in [184]) are found to be well described by a fast cooling regime. The KN effect and the pair production mechanism shape the VHE spectrum significantly. It is estimated that ≃10% of the total emitted power, i.e., ≃25% of initially produced IC power, is absorbed.

4.4. GRB 190829A

GRB 190829A is a nearby ($z=0.078$) long GRB triggered by Swift-BAT [187] and Fermi-GBM [188]. The Fermi-GBM trigger time is $T_0=19:55:53.13$ UTC. H.E.S.S. detected this GRB over three consecutive nights, with a significance of $21.7\sigma$ during the first night (∼4 h after the GRB trigger).

4.4.1. General Properties and Multi-Wavelength Observations

The prompt emission detected by the two instruments consists of two episodes, with the first one observed in the time interval from the trigger time to 4 s and the second brighter episode from 47 s to 61 s. The two episodes have very different spectral properties: the first one is described by a power-law with index $-1.41\pm 0.08$ and an exponential high-energy cutoff function with $E_p=130\pm 20$ keV, the second one can be described with a Band function with $E_p=11\pm 1$ keV, $\alpha =-0.92\pm 0.62$ and ${\beta }_T=-2.51\pm 0.01$[189]. The (isotropic equivalent) prompt emission energy inferred from the spectral analysis of Fermi-GBM data is $E_{\gamma ,iso}\sim 2\times {10}^{50}$ erg.

A multi-wavelength observational campaign of the event was performed covering the entire electromagnetic spectrum. The event was not detected in the HE range by Fermi-LAT. Nevertheless, ULs have been reported in the MeV-GeV band up to $3\times {10}^4$ s [190]. The Swift-XRT started observations at 97.3 s and detected a bright X-ray afterglow, which was monitored until ∼7.8 × ${10}^6$ s [191]. The X-ray light curve in the 0.3–10 keV energy range (observer frame) shows a peculiar behavior with an initial steep decay phase followed by a plateau and a strong flare episode (Figure 23, upper panel, blue points). After the flare, the standard afterglow phase starts with a decay following a power-law with a possible steepening around 10 days. In the UV/optical/NIR band, the event was followed by several instruments. The redshift was estimated to be $z=0.0785\pm 0.005$ [192], which makes this event one of the closest GRBs ever detected. Starting from 4.5–5.5 days after the GRB trigger, an associated supernova was reported [193]. Moreover, in the optical data, a flare is observed simultaneously with the one in the X-rays. In the radio band, the detection was reported by several instruments, starting from ∼1 day after the trigger [194,195,196,197]. The radio flux initially slowly increases and then starts to decay after 20–30 days.

4.4.2. VHE Observations and Results

The afterglow emission of GRB 190829A was followed-up by the H.E.S.S. telescopes starting at 4.3 h and continued for three consecutive nights (up to 55.9 h). The observations were performed using the four medium-size telescopes of H.E.S.S. for a total amount of 13 h divided respectively in 3.6 h (starting at 4.3 h), 4.7 h (starting at 27.2 h) and 4.7 h (starting at 51.2 h). The statistical significance at the GRB position found in the three nights are, respectively, $21.7\sigma$, $5.5\sigma$ and $2.4\sigma$.

Spectral analysis was performed for the first two nights fitting the observed photon spectrum with a power-law model. The following values are found: ${\alpha }_{obs}=-2.59\pm 0.09$ (stat.) $\pm 0.23$ (syst.) in the 0.18–3.3 TeV (first night) and ${\alpha }_{obs}=-2.46\pm 0.22$ (stat.) $\pm 0.14$ (syst.) in the 0.18–1.4 TeV energy range (second night). Fitting a power-law attenuated by EBL, the photon indices inferred for the intrinsic spectrum are: ${\alpha }_{int}=-2.06\pm 0.10$ (stat.) $\pm 0.26$ (syst.) in the 0.18–3.3 TeV energy range (first night) and ${\alpha }_{int}=-1.86\pm 0.26$ (stat.) $\pm 0.17$ (syst.) in the 0.18–1.4 TeV energy range (second night). The photon indices in each night are consistent and within the systematical uncertainties with those of the simultaneous X-ray emission. Combining all three nights, the photon index is ${\alpha }_{int}=-2.07\pm 0.09$ (stat.) $\pm 0.23$ (syst.) in the 0.18–3.3 TeV energy range.

The light-curve in the 0.2–4.0 TeV energy range derived up to 56 h is compared in Figure 24 with the XRT light-curve and the LAT upper limits. The time-evolving flux was satisfactorily modeled with a power-law decay $F\left(t\right)\propto t^{\alpha }$ with $\alpha =-1.09\pm 0.05$. Such a decay index is similar to the X-ray one derived in the same time interval (${\alpha }_X=-1.07\pm 0.09$).

4.4.3. Interpretation

The interpretation of the VHE emission from GRB 190829A is debated and different radiation mechanisms including synchrotron, SSC or EIC emission have been proposed so far to explain the origin of the TeV emission.

The H.E.S.S. Collaboration [5] investigated both the synchrotron and the SSC emission in the external forward shock as a responsible radiation mechanism of the observed TeV component. Multi-wavelength data collected simultaneously with H.E.S.S. observations in the first two nights were modeled separately with a time-independent numerical code using the Markov-chain Monte Carlo (MCMC) approach to explore the parameter space. The results of the fitting show that the SSC mechanism fails to explain the VHE emission. The low Lorentz bulk factor predicted by the observations ($\Gamma \lesssim 10$) implies that the SSC emission occurs in KN cross-scattering regime. As a result, a steep spectrum, inconsistent with the observational VHE data, is obtained (see Figure 25, light blue shaded area). Possible improvements between the data and the model foresee a higher $\Gamma$, which is in contrast with the observations, or the presence of an additional hard component in the distribution of the accelerated electrons. However, this latter solution implies extreme assumptions on the density of the circumburst medium ($n_0={10}^{-5}$ cm${}^{-3}$ in the case of the strong magnetic field or $n_0={10}^5$ cm${}^{-3}$ for weak magnetic field) and a SED strongly dominated by the SSC component, which is inconsistent with the data. A better fitting of the observational data can be obtained when considering an alternative model where the maximum electron energy set by the radiative losses is ignored. In such scenario, the synchrotron emission is able to extend up to TeV energies and the observational broad-band data are described by a single synchrotron component (see Figure 25, orange shaded area). The SSC contribution is negligible while the $\gamma -\gamma$ absorption shapes the VHE spectrum. The single synchrotron component scenario provides a better fit (>$5\sigma$) to the multi-wavelength data. On the other hand, this interpretation requires unknown acceleration processes or non-uniform magnetic field strength in the emission region, as described for GRB 180720B (see Section 4.2).

A complete multi-wavelength modeling of the GRB 190829A data, including contribution of the synchrotron and SSC emission for both the forward and reverse shocks, and considering a constant-density environment, is presented in [72]. The predicted broad-band light curves and the SED at the time of the H.E.S.S. detection are shown in Figure 23. A MCMC approach was adopted in order to estimate the best-fit parameters for the multi-wavelength modeling. The resulting values of the parameters related with the forward shock scenario are shown in Table 3. In contrast with the H.E.S.S. Collaboration results, the VHE emission is well reproduced with the SSC external forward shock scenario. The usual simplified assumption that ${\xi }_e=1$ is excluded by the fit, which provides acceptable solutions only for ${\xi }_e\lesssim 6.5\times {10}^{-2}$. Moreover, an isotropic-equivalent kinetic energy at the afterglow onset $E_k=2.5_{-1.3}^{+1.9}\times {10}^{53}$ erg is estimated. Considering the observed GBM prompt energy, such a value implies that the prompt efficiency is $\eta =1.2_{-0.5}^{+1.0}\times {10}^{-3}$, which is much lower than the typical values derived from the previous GRB studies. The other parameters ($n_0$, ${ϵ}_e$ and ${ϵ}_B$) are found to be similar to the ones estimated for GRB 190114C.

A two-component off-axis jet model has also been investigated [198]. Such a model proposes that the GRB jet is observed off-axis (${\theta }_{view}=1.{78}^{\circ }$) and it consists of a narrow (${\theta }_{jet}=0.{86}^{\circ }$) fast ($\Gamma =350$) jet and a slow ($\Gamma =20$) co-axial jet. The former jet component is responsible for the emission of SSC photons in the VHE band. The calculation of the SSC flux at the time of the H.E.S.S. detection is conducted following the prescriptions of [48], considering only the Thompson scattering regime.

An EIC plus SSC scenario has also been proposed for the production of the VHE component [199]. The seed photons belong to the long-lasting X-ray flare observed for GRB 190829A, which can be up-scattered to TeV energies. A numerical calculation of the afterglow dynamics and radiative processes have been used to model the observational data. For $t\sim {10}^3$–${10}^4$ s, the EIC component dominates the VHE emission, while for later times ($t\gtrsim 3\times {10}^4$ s) the EIC gradually decays and the SSC component becomes relevant. The initial afterglow kinetic energy used for the modeling ($E_k={10}^{52}$ erg) suggests that GRB 190829A is not a typical low-luminosity GRB but it may have much higher kinetic energy.

4.5. GRB 201015A

GRB 201015A is a long GRB at $z=0.426$ detected by the Swift-BAT [200] on 15 October 2020, at $T_0=22:50:13$ UT. The Fermi-GBM instrument did not trigger the event, but the targeted search revealed a transient source consistent with the Swift-BAT location [201]. MAGIC observations show a possible detection with a significance of ∼3.5$\sigma$.

4.5.1. General Properties and Multi-Wavelength Observations

The (isotropic equivalent) prompt emission energy inferred from spectral analysis of Fermi-GBM data is $E_{\gamma ,iso}=(1.1\pm 0.2)\times {10}^{50}$ erg [202]. The prompt duration is $T_{90}=9.78\pm 3.47$ s (15–350 keV band). The BAT time-average spectrum in the time interval 0–10 s is well fitted by a power-law model with a photon index ${\beta }_T=-3.03\pm 0.68$, suggesting a low peak energy $E_p<10$ keV [203].

Swift-XRT [204] follow-up the event starting only 3214 s after $T_0$ due to observational constraints. The light curve up to almost 1 day is well described by a power-law with decay index $\alpha =-1.{49}_{-0.21}^{+0.24}$. Late-time observations performed by the Chandra X-ray Observatory [205] and Swift-XRT [206] from ∼8 days up to ∼21 days showed a flattening of the X-ray light curve, i.e., a flux level inconsistent (higher) with the extrapolation of the power-law decay rate at early times. Optical observations confirmed the presence of an afterglow counterpart from around 168 s [207]. The optical light-curves showed a clear initial rise with a peak around 200 s followed by a decay [208]. A bright radio counterpart (flux density ∼1.3 $\times {10}^{-4}$ Jy at 6 GHz 1.4 days after the burst) was also detected by several instruments [209,210,211]. Late-time optical observations identified an associated supernova rising from 5 days after the burst reaching its maximum flux around 12–20 days after $T_0$ [212,213]. The measurement of the redshift was reported by the GTC ($z=0.426$) [214] and then confirmed by the NOT ($z=0.423$) [215] instrument.

4.5.2. VHE Observations and Results

The final results from VHE data analysis of GRB 201015A have not been published yet. Preliminary information reported here have been released in [216,217]. Observations of GRB 201015A were performed by the MAGIC telescopes starting 33 s after the trigger time, under dark conditions, with a zenith angle ranging from 24${}^{\circ }$ up to 48${}^{\circ }$, and lasted for about 4 h. In the second half of the data, taking the presence of passing clouds affected the observation for ∼0.45 h. These data were removed and the remaining ones were analyzed with the standard MAGIC analysis software. Offline analysis showed a possible excess with a $3.5\sigma$ significance at the GRB position (see Figure 26) and a significant spot in the sky map. The energy threshold of the analysis is calculated to be 140 GeV from Monte Carlo simulated $\gamma$-ray data.

4.6. GRB 201216C

GRB 201216C is a long GRB at $z=1.1$ triggered by the Swift-BAT at $T_0=23:07:31$ UT on 16 December 2020 [218]. Fermi-GBM also detected the event with a slightly different trigger time (6 s before the Swift-BAT) [219]. MAGIC detected GRB 201216C with a significance of ∼6$\sigma$.

4.6.1. General Properties and Milti-Wavelength Observations

The duration is estimated as $T_{90}=48\pm 16$ s in the 15–350 keV band by Swift-BAT [220] and around $29.9$ s in the 50–300 keV band by Fermi-GBM. The light curve shows a multiple peak structure with a main peak around 20 s after the trigger time. The time-averaged GBM spectrum in the first 50 s is best fit by a Band function with $E_p=326\pm 7$ keV, $\alpha =-1.06\pm 0.01$ and ${\beta }_T=-2.25\pm 0.03$. The isotropic equivalent energy $E_{\gamma ,iso}$ in the 10–1000 keV band is $(4.71\pm 0.16)\times {10}^{53}$ erg, as calculated from the fluence measured by Fermi-GBM.

Fermi-LAT observed the GRB starting from around 3500 s and up to 5500 s. No significant emission was reported [221]. Swift-XRT began the observation at $t=2966.8$ s due to an observational constrain. A fading source was detected following a broken power-law behavior with decay indices of $1.{97}_{-0.09}^{+0.10}$ and $1.{07}_{-0.10}^{+0.15}$ and a break at 9078 s [222]. Optical observations were also performed by several instruments. The r-band light curve made along with VLT data point [223] and inferred data from FRAM-ORM [224] show a power-law decay in flux with the index equal to 1. The Liverpool Telescope observations, performed around 177 s after the trigger time, seems to be around the peak of the optical afterglow [225]. The HAWC observatory followed-up the event but no significant detection was identified in the TeV band [226]. Redshift estimation of $z=1.1$ was performed by the ESO VLT [227].

4.6.2. VHE Observations and Results

Final results from VHE data analysis of GRB 201216C have not been published yet. Preliminary information reported here have been released in [228,229]. MAGIC observations and data taking of GRB 201216C started with a delay of 56 s after the Swift-BAT trigger time. The observation lasted for 2.2 h and was performed in optimal atmospheric condition and in absence of the moon. The zenith angle ranged from 37${}^{\circ }$ to 68${}^{\circ }$. The low level of night sky background was also allowed to retain the low energy events and therefore obtain a low energy threshold compared to the other GRBs observed. To keep as many low-energy events as possible, an image cleaning method to extract dimmer Cherenkov showers initiated by gamma rays than the standard method was adopted.

The signal significance was calculated to be $6.0\sigma$ pre-trial ($5.9\sigma$ post-trial⁷) for the first 20 min of observation (see Figure 27). A preliminary time-integrated spectrum for the first 20 min of observation was produced. Due to the strong absorption effect by EBL, a very steep power-law decay was found for the observed spectrum, especially for the events with energies higher than a few hundred of GeV. The intrinsic spectrum, corrected for the EBL absorption, was found to be consistent with a flat single power-law until 200 GeV above, in which no significant spectral points have been derived. A preliminary light curve in the time interval 56 s–2.2 h was also calculated. After 50 min, only upper limits on the emitted flux have been derived, as no significant emission was found after this time. The preliminary results are consistent with a monotonically decaying light curve fitted with a power-law.

5. The New TeV Spectral Window: Discussion

After decades of searches, MAGIC and H.E.S.S. observations have unequivocally proven that (long) GRBs can be accompanied by a significant amount of TeV emission during the afterglow phase. Table 4 summarizes the main properties of the GRBs detected by IACTs, and presented in detail in the previous section. The list also includes two events (namely GRB 160821B and GRB 201015A), where only a hint of excess (i.e., with a significance at ∼3–4 $\sigma$) was found. For the other four events, namely GRB 180720B, GRB 190114C, GRB 190829A and GRB 201216C, the detections are robust (>5$\sigma$). The table lists several properties, such as duration $T_{90}$ and total emitted energy $E_{\gamma ,iso}$ of the prompt emission, redshift, and information on the IACT detection (the starting time $T_{delay}$ of observations elapsed since the trigger time $T_0$, the energy range where photons have been detected, the name of the telescope and the significance of the excess). GRB 160821B is the only one belonging to the short class; the other five being long GRBs.

In this section, we address the question why these GRBs have been detected, whether they have peculiar properties and whether they show some common behaviors that may be at the basis of the production of TeV radiation. To do that, one should be careful, since these GRBs have been followed-up under very different observational conditions and with very different time delays after the trigger time, and they span quite a large range of redshifts (from 0.078 to 1.1). Keeping in mind these differences, which have a strong impact on the detection capabilities of IACTs, we compare the observed and intrinsic properties of the population of GRBs at VHE, highlighting their similarities and differences, and discuss how they compare to the whole population.

5.1. Observing Conditions

Low zenith angles, fast repointing, dark nights, low redshift, and highly energetic events have always been considered as optimal, if not necessary, conditions to have some chances for GRB detections with IACTs. On the other hand, these first VHE GRBs have demonstrated that GRBs can have a level of TeV emission large enough to be detected by the current generation of IACTs, even under non-optimal conditions. GRB 190114C was observed with a zenith angle > ${55}^{\circ }$ and in the presence of the moon. Both conditions imply a higher energy threshold (typically ≳ 0.2 TeV) and require a dedicated data analysis. Another example is GRB 160821B, that was observed with a NSB 2–8 times higher than the standard dark night conditions. Moreover, significant VHE excess was found not only in case of short delays (less than hundreds of seconds) from the burst trigger but, somewhat surprisingly, also at quite late times, i.e., with delays of several hours or even days, as in the case of GRB 180720B and GRB 190829A, respectively. This showed the importance of pointing a GRB also at relatively late times, in cases fast follow-up observations are not feasible.

Optimal observing conditions and short delays remain, however, crucial to detecting GRBs at a higher redshift, for which the impact of EBL is already large at a few hundreds GeV. This explains how the detection of a GRB has been possible at $z=1.1$ (GRB 201216C): in this case, optimal observing conditions allowed us to reach a low energy threshold of the sensitivity window (∼70 GeV). The excess of signal was indeed found only below 200 GeV (more precisely, between 70 and 200 GeV) where the attenuation by the EBL is still limited.

5.2. Redshift and the Impact of EBL

The redshift of the detected GRBs covers a broad range, from $z=0.079$ (GRB 190829A) to $z=1.1$ (GRB 201216C). The impact of the EBL attenuation on the spectrum is severely changing, depending on the redshift value and on the photon energy. For redshift z ∼ 0.4, the impact becomes relevant for energies ≳ 0.2 TeV with a flux attenuation of ∼50% for $0.2$ TeV and almost ∼99.5% for 1 TeV [183]. For nearby events ($z\lesssim$ 0.1–0.2), the effect of EBL is less severe and becomes relevant only for energies ≳ 0.3 TeV, reaching an attenuation factor of an order of magnitude only for energies ≳ 2 TeV. As a result, the GRB observed photon indices and the energy range of detected TeV photons differs significantly between the events. GRBs with redshift $z>0.4$ such as GRB 190114C or GRB 180720B have very steep photon indices and they are detected in the lower energy range up to 0.44 TeV for GRB 180720B and 1.0 TeV for GRB 190114C. Spectral analysis from GRB 201216C are not yet public but preliminary results indicate that the emission is concentrated in the lower energy band between 0.1–0.2 TeV. Nearby GRBs with redshift $z\lesssim$ 0.1–0.2 such as GRB 160821B or GRB 190829A show a less steep photon spectrum (around $-2.5$) and the TeV detection range extends above 1 TeV. The detection of several GRBs with significant value of redshift ($z>0.4$) is robust proof that IACTs can overcome the limitations due to the EBL absorption and can expand the VHE detection horizon at the current stage up to $z=1.1$. On the other hand, it is evident that detection of nearby GRBs is fundamental in order to more robustly explore the spectral shape, unbiased by the EBL effect, which is a non-negligible source of uncertainty for higher redshifts.

5.3. Energetics

In terms of $E_{\gamma ,iso}$, the VHE GRB sample spans more than three orders of magnitude from ∼${10}^{49}$ erg up to ∼6 $\times {10}^{53}$ erg. The five long GRBs detected follows the Amati correlation, as shown in Figure 28. GRB 160821B, the only short GRB of the sample, is consistent with the existence of a possible Amati-like correlation for short GRBs, with this event falling in the weak-soft part of the correlation. The detections of GRB 190829A and GRB 201015A show that an event does not need to be extremely energetic in terms of isotropic-equivalent prompt energy in order to produce a TeV emission with (intrinsic) luminosity comparable to the X-ray luminosity. As a result, sources with $E_{\gamma ,iso}\sim {10}^{50-51}$ erg are not excluded as possible TeV emitters, even though their detection is possible only for relatively low redshift. This reduces the available volume, and hence the detection rate of similar events. In any case, this is relevant also for short GRBs which are less energetic than long ones, with typical isotropic energies falling within the ∼10${}^{49-52}$ erg range [144].

5.4. X-ray Lightcurves

The comparison between TeV and X-ray light-curves suggests an intimate connection between the emission in these two bands, both in terms of emitted energy and luminosity decay rate. In Figure 29, the XRT afterglow light-curves (luminosity versus rest-frame time) in the 0.3–10 keV energy range are compared with the VHE light-curves (integrated over different energy ranges, depending on the detection window, see Table 4). Different colors refer to the six different GRBs. The VHE luminosity is shown with empty circles.

Considering the X-ray luminosity, the GRB sample can be divided into two groups: GRB 190114C, GRB 180720B and GRB 201216C display large and clustered X-ray luminosity (at $t\sim {10}^4$ s their luminosity is around 1–5 $\sim {10}^{47}$ erg s${}^{-1}$) and their light curves almost overlap for the entire afterglow phase. The other three GRBs (GRB 190829A, GRB 201015A and GRB 160821B) are much fainter in terms of X-ray luminosity (at least two orders of magnitude at $t\sim {10}^4$ s). This is consistent with the fact that they also have a smaller $E_{\gamma ,iso}$. The correlation between X-ray afterglow luminosity and prompt $E_{\gamma ,iso}$ is found in the bulk of the long GRB population, and these GRBs make no exception.

Observations in the VHE band (empty circles in Figure 29) reveal that the VHE luminosities observed in the afterglow phase are in general smaller but comparable to the simultaneous X-ray luminosity, implying that almost an equal amount of energy is emitted in the two energy bands. Any theory aimed at explaining the origin of the TeV radiation should explain the origin of these similarities. Concerning the decay rate, observations are still not conclusive. The decay rate of the TeV emission is available only for two events. For GRB 190829A, the temporal indices in X-ray and VHE are very similar, while for GRB 190114C, the VHE clearly decays faster than the X-ray emission.

For GRB 190114C at $t\sim 380$ s, the VHE luminosity $L_{VHE}$ is ∼1.5–2.5 $\times {10}^{48}$ erg s${}^{-1}$ and the X-ray one $L_X$ is ∼0.6–1.0 $\times {10}^{49}$ erg s${}^{-1}$. As a result, the power radiated in the VHE band is about ∼25% of the X-ray one. Similarly for GRB 190829A at $t\sim 4.5$ h the VHE luminosity $L_{VHE}$ is ∼4.0–8.5 $\times {10}^{44}$ erg s${}^{-1}$, which is around ∼15–20% of the corresponding X-ray one ($L_X$∼ 2.0–5.0 $\times {10}^{45}$ erg s${}^{-1}$). For GRB 180720B at $t\sim 2\times {10}^4$ s the VHE luminosity $L_{VHE}$ is ∼9 $\times {10}^{47}$ erg s${}^{-1}$ and the X-ray one $L_X$ is ∼1.5–2.5 $\times {10}^{49}$ erg s${}^{-1}$. In this case, the power radiated in the VHE band is around ∼35–60% of the X-ray one.

5.5. The TeV Contribution to the Multi-Wavelength Modeling

Modeling of multi-wavelength afterglow data provide important insights concerning the GRB afterglow physics. In particular, the VHE data were crucial to investigate (i) the radiation mechanisms responsible for the production of photons between 10–100 GeV already detected by LAT; (ii) the environmental conditions at the GRB site; (iii) the free parameters which describe the shock micro-physics, and in particular, the self-generated magnetic field.

The modelings proposed so far in literature to explain the VHE component have considered two different radiation processes at the origin of the TeV emission: SSC and synchrotron. In the first case, the VHE emission is interpreted as a distinct spectral component from the synchrotron radiation dominating from radio to ∼ GeV energies, which provides the seed photons that are upscattered at higher energies by the same electron population. In the second scenario, the VHE emission is observed as the extension of the synchrotron spectrum up to TeV energies.

In principle, a simultaneous SED covering the X-ray, HE and VHE range should be sufficient to discriminate among these two different possibilities. A hardening of the spectrum from GeV to TeV energies should be the smoking gun for the presence of a distinct component. In reality, the uncertainties in the spectral slope at VHE (caused by the uncertainty on the EBL and on the narrow energy range of TeV detection) can make the distinction hard to perform. In this case, LAT observations are of paramount importance to reveal the presence of a dip in the SED, which would also prove the need to invoke a different origin for the VHE emission. This seems to be the case for GRB 190114C, for which, at least in one SED, the LAT flux strongly suggests a dip in the GeV flux and hence the presence of the characteristic double bump observed for a synchrotron-SSC emission (Figure 20). For GRB 190829A, LAT provides only an upper limit, which is not constraining for modeling the shape of the SED (Figure 25). In this GRB, an interpretation of the whole SED in terms of synchrotron radiation cannot be excluded, although a modeling as SSC radiation has been proven to be successful [72] (Figure 23). For the other events detected at VHE, either the data do not allow for building a proper SED with simultaneous multi-wavelength observations, or they are not yet public. Despite this, the SSC emission seems to be the most viable mechanism able to explain the TeV data. A firm conclusion on the responsible radiation mechanism has not been reached yet and future detections will be crucial for deeper investigations.

Assuming one of the two scenarios, TeV data coupled with broad band observations at lower energies can be exploited to give additional information on the details of the afterglow external forward shock scenario.

Concerning the shock micro-physics, several modelings have suggested the possibility that the fraction of electrons accelerated in non-thermal distribution, ${\xi }_e$, is different from the standard value of one, which is usually assumed. In a few GRBs’ modelings, namely for GRB 190114C [185] and GRB 190829A [72,198], the introduction of a ${\xi }_e<1$ was essential in order to consistently fit the observational data. In particular in [72], the requirement for a low value of ${\xi }_e\lesssim 6.5\times {10}^{-2}$ was required to provide an acceptable fit of the data. The other modelings assume a greater value of ${\xi }_e$, around ∼0.3. Further detections will be exploited in order to verify if such an indication could be present also in other events.

Some considerations can also be drawn for the equipartition parameters ${ϵ}_e$ and ${ϵ}_B$. These values, especially the former one, are usually well unconstrained and can span several orders of magnitudes. The TeV modeling described so far suggest that around ∼10% of the energy is given to the electrons, while a lower value (from ${10}^{-5}$ to ${10}^{-3}$) is given to the magnetic field. larger values of ${ϵ}_B$ such as 0.1–0.01, which are considered in an external shock scenario, are excluded. Moreover, some results can also be interpreted as an indication of an evolution in time of these parameters. In Figure 22, the modeling of the broad band light curves of GRB 190114C is shown. Two different modeling are presented: one optimized for the early time X-ray, HE and VHE observations (solid line) and one optimized for the late time lower energy bands (dotted line). This is due to the fact that the model that reproduces the early time data over-predicts the late time optical and radio observations. This result points towards the possibility that some of the fixed parameters of the afterglow theory (e.g., the electron and magnetic field equipartition parameters) may evolve in time. A further clue of the presence of time-dependent shock micro-physics parameters is derived from the low frequency multi wavelength modeling of GRB 190114C presented in [121]. In order to model the optical and the radio data, it is required that the micro-physical parameters evolve with time as ${ϵ}_e\propto t^{-0.4}$ and ${ϵ}_B\propto t^{0.1}$ in the ISM case and ${ϵ}_B\propto t^{0.76}$ for the stellar wind scenario.

An issue that still is not solved by TeV observations is the discrimination between constant and wind-like profile for the ambient density. It is expected that long GRBs occur in wind-like environments. Nevertheless, at the current stage there seems to be no preference between such an environment and a constant ISM one, which is able to well reproduce the observational data. Therefore conclusive answers on the topic cannot be drawn yet.

In conclusion, the current population of GRBs at VHE already demonstrate quite broad properties, spanning more than three orders of magnitude in $E_{\gamma ,iso}$ and more than two orders of magnitude in terms of afterglow luminosity and ranging in redshift between 0.079–1.1. The afterglow X-ray and VHE emission have comparable fluxes and decay slopes. The afterglow emitted power in the VHE band seems to constitute from $15\%$ up to $60\%$ of the X-ray one. Data modeling suggest that the responsible VHE radiation mechanism is the SSC emission, although different mechanisms (e.g., synchrotron radiation, EIC) cannot be completely excluded and a conclusive answer cannot be given yet. Multi-wavelength modeling show no preferences concerning the GRB environments between an ISM or wind-like scenario and indicate that shock micro-physics parameters, which seem to be able to reproduce VHE emission, are ${\epsilon }_e\sim 0.1$ and ${\epsilon }_B\sim {10}^{-5}$–${10}^{-3}$. Such features can be an indication of the universality of TeV emission in GRBs. It is then expected that a larger sample of GRBs than the current one will be detected in the VHE band, including also short GRBs for which, at the current stage, there are no confirmed detections except for the hint of excesses observed for GRB 160821B.

6. Conclusions and Future Prospects

The recent discoveries performed by the current generation of Cherenkov telescopes in the VHE band have opened a new observational spectral window on GRBs. The presence of a TeV afterglow component has been unequivocally proven and the studies on the currently available samples have demonstrated the potential that such detections have in probing several long-standing open questions in the GRB field. These first studies have focused on the identification of the responsible radiation mechanism, which is the first issue to address, and the comparison of the energetics, luminosity, and temporal behavior of the VHE component with respect to emission at lower frequencies. A modeling of multi-wavelength data covering from radio up to TeV energies was performed, giving interesting insights on the shock micro-physics conditions.

Limitations to the robust use of VHE data for afterglow modeling are imposed by the severe modification of the intrinsic spectrum cased by the energy-dependent flux-attenuation induced by EBL. GRBs with redshift $z>0.4$, four out of six in the current VHE sample, are strongly affected by EBL absorption starting from hundreds of GeV. This implies large uncertainties on the shape and photon index of the intrinsic VHE spectrum. As a result, firm conclusions on the origin and spectral regime of the TeV component cannot be drawn yet. The low-energy extension of the range of sensitivity of IACTs is then fundamental for reaching a larger rate of detections and a more robust determination of the spectral index of the TeV component.

A full comprehension and exploitation of TeV data is expected to be reached thanks to the next generation of Cherenkov telescopes. The Cherenkov Telescope Array (CTA) will be a huge step forward for the detection of GRBs in the VHE band. The major upgrades with respect to the current generation of Cherenkov telescopes that will impact GRB observations are: (i) a lower energy threshold (≲30 GeV), (ii) a larger effective area at multi-GeV energies ($\sim {10}^4$ times larger than Fermi-LAT at 30 GeV) and (iii) a rapid slewing capability (180 degrees azimuthal rotation in 20 s). Moreover, its planned mixed-size array of large, medium and small size telescopes (called LST, MST and SST, respectively) situated at two sites in the northern and southern hemispheres will provide a full sky coverage from few tens of GeV up to hundreds of TeV. CTA will have a much better sensitivity and a broader energy range with respect to current ground-based facilities. A comparison is shown in Figure 30. At the present stage, the first prototype of the LSTs has been built and is operative under a commissioning phase⁹ in the northern site at the Roque de los Muchachos Observatory in La Palma. Despite these performance improvements, the expected CTA detection rates of GRBs will be influenced anyway by the relatively low duty cycle affecting IACTs and by the synergies with other instruments. Indeed, Cherenkov telescopes’ repointing relies on external triggers coming from space satellites. Assuming that currently operating space telescopes will be still operative, GRB alerts will be mostly provided by Swift-BAT and partially by the Fermi-GBM, and in the future by the French-Chinese mission Space-based multi-band astronomical Variable Objects Monitor (SVOM [232]).

BAT observes around 92 GRBs per year with a typical localization error of a few arcmin [233]. The good localization (later refined by XRT to a few arcsec) is fundamental for Cherenkov telescopes, given their limited field of view (e.g., about $4^{\circ }$ for the LSTs and $7^{\circ }$ for the MSTs). The GBM provides a much higher number of alerts, around 250 per year but with a larger localization error, from 1–3${}^{\circ }$ up to ${10}^{\circ }$, which makes follow-up with IACTs very challenging. In the case of such large localization errors, CTA can exploit the so-called divergent mode for observations, which is currently under study [234]. In this pointing strategy, each telescope points to a position in the sky that is slightly offset to extend the field of view. Concerning future instruments, SVOM is expected to provide Swift-like alerts at a rate of ∼60–80 GRBs/yr with a localization error <$1^{\circ }$, including also 10 GRBs/yr with redshift $z<1$.

Available estimates of the CTA detection rate of GRBs are reported in [235]. These studies were performed before the discovery of TeV emission. They are based on Swift-like alerts (triggered by Swift-BAT or SVOM) and Fermi-GBM alerts. The predicted detection rate is around a few GRBs per year, depending on the energy threshold of the observation and on the observation delay [235]. An updated study that considers current knowledge of TeV emission in the afterglow of GRBs is in progress [236].

Despite for decades GRBs’ hunting by Cherenkov telescopes has been primarily focused on reaching low energy thresholds in order to explore the multi-GeV band, these first detections have shown that photons above TeV energies can be produced in GRBs and can be detected. This is mostly valid only for nearby events of redshift below 0.1–0.2, where EBL attenuation is not too severe. The exploration of the GRB emission component above 1 TeV can be of potential interest for SSTs and for the ASTRI Mini-Array. The ASTRI Mini-Array, currently under construction, will be an array of nine imaging atmospheric dual-mirror Cherenkov telescopes at the Teide Observatory site, expected to deliver the first scientific results in 2023. After the detection of GRB 190114C, the capabilities of the ASTRI Mini-Array in detecting and performing spectral studies of an event similar to the MAGIC GRB have been explored [237]. GRB 190114C has been taken as a template to simulate possible GRB emission from a few seconds to hours, and has been extrapolated to 10 TeV on the bases of model predictions.

The results demonstrate that the instrument will be able to detect afterglow TeV emission from an event such as GRB 190114C up to $\sim 200$ s (see the comparison between the GRB observed flux at 1 TeV and the differential ASTRI Mini-Array sensitivity in Figure 31). By moving the GRB at a smaller redshift (down to $z=0.078$, the redshift of the TeV GRB 190829A), the time for which the GRB is detectable increases up to ∼${10}^5$ (however, the light-curve in this case should be re-scaled by the lower energetics of nearby events). Nearby GRBs are then potential target of interests for the ASTRI Mini-Array. These are certainly rare events, but their detection will provide a wealth of information, with spectra that can be characterized up to several TeV [237].

In conclusion, after decades of huge efforts, current ground-based VHE facilities have started a new era in the comprehension and study of GRB physics. Their breakthrough detections allow unprecedented studies. As discussed in this review, many open questions in afterglow physics can largely benefit from the inclusion of TeV data. The first detections are providing glimpses of such a huge potential. Luckily, we are at the dawn of the VHE era thanks to the upcoming CTA observatory, which will assure major upgrades in sensitivity, energy range, temporal resolution, and sky coverage. Future observations, if complemented by simultaneous observations in X-rays and at ∼ GeV energies, will play a paramount role to improve our knowledge on the physics of GRB during the afterglow phase and hopefully also in the prompt phase. In particular, the afterglow SSC one-zone model will be tested to understand whether it can grasp the main properties of the VHE emission or if a revision of our comprehension on the particle acceleration processes, shock micro-physics and radiation mechanisms is needed.

Author Contributions

All authors, D.M. and L.N., have contributed to writing. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester (https://www.swift.ac.uk/xrt_curves/).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Notes

1. https://magic.mpp.mpg.de/ (accessed on 27 February 2022).

2. https://www.mpi-hd.mpg.de/hfm/HESS/ (accessed on 27 February 2022).

3. https://www.cta-observatory.org (accessed on 27 February 2022).

4. A trial factor of 2 is considered due to the two sets of analysis cuts used for MAGIC data analysis.

5. https://www.swift.ac.uk/xrt_live_cat/00848890/ (accessed on 27 February 2022).

6. post-trial significance is estimated by accounting for the previously well-localized GRBs observed by H.E.S.S. in the same array configuration as of GRB 180720B.

7. See note 4 above.

8. https://www.swift.ac.uk/burst_analyser/ (accessed on 27 February 2022).

9. https://www.lst1.iac.es/ (accessed on 27 February 2022).

10. https://www.cta-observatory.org/science/ctao-performance/ (accessed on 27 February 2022).

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Spectral index ([Forumla omitted. See PDF.]) of the electrons accelerated in shocks as a function of the shock velocity. Curves refer to the equation derived by [61] under the hypothesis of isotropic, small-angle scattering and is a generalization of the non-relativistic formula. Symbols show the comparison with numerical studies. Different curves refer to different assumptions on the type of shock (see [61] for details): in all cases, the value of the spectral index approaches the same value [Forumla omitted. See PDF.] (corresponding to [Forumla omitted. See PDF.]) in the ultra-relativistic limit.

Figure 2. PIC simulations: temporal evolution of the downstream spectrum of ions (upper panel) and electrons (bottom) for mass ratio [Forumla omitted. See PDF.], shock Lorentz factor [Forumla omitted. See PDF.] and magnetization [Forumla omitted. See PDF.]. The evolution is followed until [Forumla omitted. See PDF.]. Inset (a): mean post-shock ion (red) and electron (blue) energy (in units of the bulk energy of the upstream flow). The dashed blue line shows the electron energy at injection. Inset (b): temporal evolution of the maximum Lorentz factor of ions (red) and electrons (blue). For comparison, the black dashed line shows the scaling [Forumla omitted. See PDF.]. From [42]. ©AAS. Reproduced with permission.

Figure 3. Comparison between the exact annihilation rate (Equation (46)) and the approximated formula (Equation (47)). The ratio between the two curves in the range ([Forumla omitted. See PDF.]) is shown in the bottom panel. In this range, the ratio is always [Forumla omitted. See PDF.].

Figure 4. Spectral energy distribution in the GRB afterglow external forward shock scenario estimated with the numerical code presented in this section for [Forumla omitted. See PDF.] s. The effect of the pair production attenuation is clearly evident in the VHE tail. The set of afterglow parameters used are the following: [Forumla omitted. See PDF.] erg, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] cm[Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 5. Comparison between spectra estimated with analytical approximations and numerical calculations. For the analytical method, the synchrotron emission (green dashed line) is estimated starting from [38] (SPN98, see legend) while the SSC (black dashed line) is taken from [48] (SE01). Vertical lines mark the break frequencies. The results from the numerical code described in Section 2.3 are shown with solid blue and red lines. This example shows the spectrum calculated at [Forumla omitted. See PDF.] s for [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] erg, [Forumla omitted. See PDF.] cm[Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 6. Average fluxes of Swift-XRT flares in the 0.3–10 keV energy range versus peak times [Forumla omitted. See PDF.]. The blue points are the 468 flares observed by Swift from April 2005 to March 2015 (collected from [87]). The orange line is the resulting linear regression, which gives the expression: [Forumla omitted. See PDF.].

Figure 7. Differential sensitivity as a function of the observation time for several HE and VHE instruments (Fermi-LAT, MAGIC, VERITAS and CTA) at four selected energies (75, 100, 150, 250 GeV). From [100].

Figure 8. Density profile produced by a star wind as a function of the distance from the central object. The [Forumla omitted. See PDF.] profile characterizes the region where the wind freely streams, which is separated from the shocked wind (with nearly constant density) by the termination shock. Top panel: the impact of different values of the ISM density on the termination shock is shown. For all the curves, it is assumed [Forumla omitted. See PDF.] yr[Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] km s[Forumla omitted. See PDF.], while the ISM density varies from 10 to [Forumla omitted. See PDF.] cm[Forumla omitted. See PDF.], with the termination shock moving to smaller distances when the density increases. Bottom panel: for a fixed ISM density [Forumla omitted. See PDF.] cm[Forumla omitted. See PDF.], the impact of different [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] is shown. Black-solid: [Forumla omitted. See PDF.] yr[Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] km s[Forumla omitted. See PDF.]; red-solid: [Forumla omitted. See PDF.] yr[Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] km s[Forumla omitted. See PDF.]; black-dotted: [Forumla omitted. See PDF.] yr[Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] km s[Forumla omitted. See PDF.]; red-dotted: [Forumla omitted. See PDF.] yr[Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] km s[Forumla omitted. See PDF.]. From [102]. ©AAS. Reproduced with permission.

Figure 9. Distribution of accelerated electrons from a relativistic shock. The two bumps correspond, respectively, to the distribution of the [Forumla omitted. See PDF.] fraction of particles heated into a thermal component with post shock energy [Forumla omitted. See PDF.] and the fraction [Forumla omitted. See PDF.] of particles accelerated into a non-thermal component with post shock energy [Forumla omitted. See PDF.]. Adapted from [73]. ©AAS. Reproduced with permission.

Figure 10. GRB 160821B: angular distance distribution [Forumla omitted. See PDF.] between the nominal source position and the reconstructed event arrival directions. The gray histogram represents the background events, while the black points with blue crosses are the [Forumla omitted. See PDF.]-like events. The vertical dashed line marks the [Forumla omitted. See PDF.] cut value and defines the region in which excess events and signal significance are calculated. From [145]. ©AAS. Reproduced with permission.

Figure 11. GRB 160821B: multi-wavelength light curves (from radio to TeV) and their modeling according to [145]. ©AAS. Reproduced with permission. The different contributions from the forward shock (FS), reverse shock (RS) and kilonova are shown (see legend).

Figure 12. GRB 160821B: modeling of the simultaneous multi-wavelength SED at approximately ∼3 h according to [145] (©AAS. Reproduced with permission.), for the same parameters used to model the lightcurves in Figure 11. The shaded areas show the sensitivity energy range of the different instruments. The MAGIC error box on the reconstructed flux is also shown. Synchrotron (solid black line), intrinsic SSC (before EBL absorption, dashed black line) and SSC emission with EBL attenuation (solid red line) estimated from the numerical modeling are shown.

Figure 13. GRB 160821B: modeling of multi-wavelength light curves (top panel) and simultaneous SED (bottom panel) in radio, optical, X-ray and VHE band presented by [157]. ©AAS. Reproduced with permission. Two scenarios are considered: SSC (dashed orange lines) and EIC (dash-dotted green lines). In the top panel, the flux contribution in the different energy bands due to the forward shock (FS) both in the SSC and in the EIC scenario, the extended X-ray emission (EE) and the X-ray plateau emission (PL) are shown. In the bottom panel for the VHE band, both the intrinsic (thin lines) and the absorbed (thick lines) VHE curves are shown for the SSC and EIC scenarios for three different observed times. For the lower energies, the broken power-law of the late prompt emission (dashed violet line) and the expected synchrotron SEDs (solid red line) are displayed. The MAGIC observations (green symbols and green shaded area) are absorbed by EBL.

Figure 14. GRB 180720B: broad-band light curves and photon index evolution including optical, X-ray and VHE data. For the late time H.E.S.S. observation (18 days after the trigger) an UL has been derived. A flux temporal decay with index [Forumla omitted. See PDF.] (similar to what observed in optical and X-ray) is also shown in black dashed line. From [4].

Figure 15. GRB 180720B: spectral fit to the observed emission in the 0.1–0.44 TeV energy range. In the left panel (a), the observed spectrum is assumed to be described by a power-law model, while in the panel on the right (b), the spectrum is an power-law attenuated by EBL. Statistical and systematic uncertainties at a [Forumla omitted. See PDF.] confidence level are shown with shaded areas. From [4].

Figure 16. GRB 180720B: modeling of the broad-band light curves (left panel) and SED at the time of the H.E.S.S. detection (right panel) proposed by [174]. Both the synchrotron and the SSC contribution to the total flux are shown (see legend). In the SED, X-ray and H.E.S.S. data are shown, respectively, with the green and the blue boxes.

Figure 17. GRB 190114C: light-curves at different frequencies. From [71].

Figure 18. GRB 190114C. Upper panel: MAGIC light-curve (red circles), compared with the XRT (green band) and LAT (red band) emission. Bottom panel: temporal evolution of the intrinsic spectral photon index in the MAGIC data analysis time bins. From [3].

Figure 19. GRB 190114C: comparison between the distribution of the observed [Forumla omitted. See PDF.]-ray events binned in time and energy (blue shaded areas) and the limiting curves for synchrotron maximum energy. Two different density radial profiles are considered for the derivation of the theoretical curves: constant (red-dotted curve) and wind-like (red-dashed curve). From [3].

Figure 20. GRB 190114C SEDs from soft X-rays to 1 TeV in five different time intervals. The contour regions for flux uncertainties are also shown. For MAGIC and LAT the 1[Forumla omitted. See PDF.] error of their best-fit power-law functions are displayed. For Swift data, the 90% confidence contours for the XRT-BAT joint fit is shown. MAGIC spectra are EBL-corrected assuming the model by [183]. From [71].

Figure 21. GRB 190114C: modeling of two SEDs in consecutive time intervals. The different curves refer to: the observed spectrum (thin solid line), the EBL deabsorbed spectrum (thick blue line) and the SSC component neglecting the effects of internal [Forumla omitted. See PDF.]-[Forumla omitted. See PDF.] opacity (dashed line). From [71].

Figure 22. Modeling of broad band GRB 190114C light curves. Description of the different modelings used is given in Section 4.3.3. From [71].

Figure 23. GRB 190829A: modeling of multi-wavelength light curves (upper panel) and SED (bottom panel) proposed by [72]. The 90% and 50% credible intervals from the fit are shown in lighter shades.

Figure 24. GRB 190829A: multi-wavelength light curves (A) (upper panel) and photon index evolution (bottom panel) in the X-ray, HE and VHE band (B). The BAT prompt light curve is shown in the inset (C). From [5].

Figure 25. GRB 190829A: modeling of X-ray, LAT and H.E.S.S. data proposed by the H.E.S.S. collaboration for the two time intervals with VHE and X-ray detections. Two scenarios are investigated for the TeV emission: synchrotron and SSC. H.E.S.S. flux contours are displayed considering the statistical uncertainty. The synchrotron and SSC component are shown in dashed and dash-dotted lines, respectively. The shaded areas represent the 68% confidence intervals determined from the posterior probability distribution of the MCMC parameter fitting for the standard SSC model (light blue) and for the model without maximum energy for synchrotron emission (orange). From [5].

Figure 26. GRB 201015A: distribution of the angular distance [Forumla omitted. See PDF.] between the reconstructed event arrival directions and the nominal source position. The gray histogram represents background events while the black point with blue crosses are the [Forumla omitted. See PDF.]-like events. The vertical dashed line describes the [Forumla omitted. See PDF.] cut value and defines the region in which excess events and signal significance are calculated. From [217].

Figure 27. [Forumla omitted. See PDF.] angular distance distribution between reconstructed event arrival directions and nominal source position for GRB 201216C. The gray histogram represents background events while the black point with blue crosses are the [Forumla omitted. See PDF.]-like events. The vertical dashed line describes the [Forumla omitted. See PDF.] cut value and defines the region in which excess events and signal significance are calculated. From [229].

Figure 28. The Amati ([Forumla omitted. See PDF.]–[Forumla omitted. See PDF.]) correlation for a sample of 136 long GRBs (grey dots, from [230]) and a sample of 11 short GRBs (empty blue squares) detected by Swift. The corresponding power-law fit for the sample of long GRBs and the [Forumla omitted. See PDF.] scatter of the distribution of points around the best fits are shown. The six GRBs detected in VHE are also added in the plot. Adapted from [231].

Figure 29. X-ray and VHE luminosity versus rest-frame time for the six GRBs detected in the TeV domain. The X-ray light-curve in the 0.3–10 keV energy range is taken from the Swift Burst Analyzer webpage 8.

Figure 30. CTAO differential sensitivity10 (defined as the minimum flux needed to obtain a 5-standard-deviation detection of a point-like source) for 50 h of observations with the Northern and Southern array compared to the sensitivity of several other Cherenkov telescopes and with Fermi-LAT (1 year).

Figure 31. Light-curve of GRB 190114C at 1 TeV (dotted purple curve) compared to the sensitivity of the ASTRI Mini-Array. The yellow dashed and green dot-dashed curves show how GRB 190114C rescaled at redshift [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], respectively, would appear. From [237].

Table 1

GRB 160821B. List of the best fit parameters inferred from multi-wavelength modeling of the afterglow radiation by different authors.

	${\mathit{E}}_{\mathit{k}}$	log(${\mathit{ϵ}}_{\mathit{e}}$)	log(${\mathit{ϵ}}_{\mathit{B}}$)	log(n)	p	${\mathit{\xi }}_{\mathit{e}}$	${\mathit{\theta }}_{\mathit{j}}$
	erg			cm${}^{-3}$			rad
MAGIC Coll.	${10}^{51}$–${10}^{52}$	[−1; −0.1]	[−5.5; −0.8]	[−4.85; −0.24]	2.2–2.35	1	/
Troja + 2019	${10}^{50}$–${10}^{51}$	[−0.39; −0.05]	[−3.1; −1.1]	[−4.2; −1.7]	2.26–2.39	1	0.08–0.50
Zhang + 2021 (SSC)	$3\times {10}^{51}$	−0.52	−5	−1.3	2.3	0.5	0.15
Zhang + 2021 (EIC)	$2\times {10}^{51}$	−0.3	−6	−1	2.5	0.1	0.1

Table 2

GRB 190114C: parameters inferred by different authors from the modeling of observations with a synchrotron-SSC scenario.

	${\mathit{E}}_{\mathit{k}}$	${\mathit{ϵ}}_{\mathit{e}}$	${\mathit{ϵ}}_{\mathit{B}}$	n	p	${\mathit{\xi }}_{\mathit{e}}$
	erg			cm${}^{-3}$
MAGIC Coll.	≳3 ×${10}^{53}$	0.05–0.15	0.05–1 × ${10}^{-3}$	0.5–5	2.4–2.6	1
Wang + 2019	$6\times {10}^{53}$	0.07	$4\times {10}^{-5}$	0.3	2.5	1
Asano + 2020	${10}^{54}$	0.06	$9\times {10}^{-4}$	1	2.3	0.3
Asano + 2020	${10}^{54}$	0.08	$1.2\times {10}^{-3}$	0.1 (wind)	2.35	0.3
Joshi + 2021	$4\times {10}^{54}$	0.03	0.012	$2\times {10}^{-2}$ (wind)	2.2	1
Derishev + 2021	$3\times {10}^{53}$	0.1	2–6 $\times {10}^{-3}$	2	2.5	1

Table 3

Parameters for modeling of GRB 190829A.

	${\mathit{E}}_{\mathit{k}}$	${\mathit{ϵ}}_{\mathit{e}}$	${\mathit{ϵ}}_{\mathit{B}}$	n	p	${\mathit{\xi }}_{\mathit{e}}$	${\mathit{\theta }}_{\mathit{j}}$
	erg			cm${}^{-3}$			rad
Hess Coll. (SSC)	$2.0\times {10}^{50}$	0.91	5.9–7.7 $\times {10}^{-2}$	1.	2.06–2.15	1.	/
Hess Coll. (Sync)	$2.0\times {10}^{50}$	0.03–0.08	≈1	1.	2.1	1.	/
Salafia + 2021	1.2–4.4 $\times {10}^{53}$	0.01–0.06	1.2–6.0 $\times {10}^{-5}$	0.12–0.58	2.01	<6.5 × ${10}^{-2}$	0.25–0.29
Zhang + 2021	$9.8\times {10}^{51}$	0.39	$8.7\times {10}^{-5}$	0.09	2.1	0.34	0.1

Table 4

List of the GRBs observed by IACTs with a firm detection (significance $>5\sigma$) or a hint of detection (3–4$\sigma$) above 100 GeV. The $T_{90}$ and $E_{\gamma ,iso}$ refer to the duration and total emitted energy of the prompt emission; the redshift is listed in column 3; $T_{delay}$ is the time delay between the trigger time $T_0$ and the time when IACT observations started; $E_{range}$ defines the energy range of the detected photons. The name of the telescope which made the observation and the significance of the detection are listed in the last column.

	${\mathit{T}}_{90}$	${\mathit{E}}_{\mathit{\gamma },\mathbf{iso}}$	z	${\mathit{T}}_{\mathbf{delay}}$	${\mathit{E}}_{\mathbf{range}}$	IACT (Sign.)
	s	erg		s	TeV
160821B	0.48	$1.2\times {10}^{49}$	0.162	24	0.5–5	MAGIC ($3.1\sigma$)
180720B	48.9	$6.0\times {10}^{53}$	0.654	3.64 × ${10}^4$	0.1–0.44	H.E.S.S. ($5.3\sigma$)
190114C	362	$2.5\times {10}^{53}$	0.424	57	0.3–1	MAGIC (>50$\sigma$)
190829A	58.2	$2.0\times {10}^{50}$	0.079	1.55 × ${10}^4$	0.18–3.3	H.E.S.S. ($21.7\sigma$)
201015A	9.78	$1.1\times {10}^{50}$	0.42	33	0.14	MAGIC ($3.5\sigma$)
201216C	48	$4.7\times {10}^{53}$	1.1	56	0.1	MAGIC ($6.0\sigma$)

References

1. Kouveliotou, C.; Meegan, C.A.; Fishman, G.J.; Bhat, N.P.; Briggs, M.S.; Koshut, T.M.; Paciesas, W.S.; Pendleton, G.N. Identification of Two Classes of Gamma-Ray Bursts. Astrophys. J.; 1993; 413, L101. [DOI: https://dx.doi.org/10.1086/186969]

2. Mazets, E.P.; Golenetskii, S.V.; Ilinskii, V.N.; Panov, V.N.; Aptekar, R.L.; Gurian, I.A.; Proskura, M.P.; Sokolov, I.A.; Sokolova, Z.I.; Kharitonova, T.V. Catalog of cosmic gamma-ray bursts from the KONUS experiment data. Astrophys. Space Sci.; 1981; 80, pp. 3-83. [DOI: https://dx.doi.org/10.1007/BF00649140]

3. MAGIC Collaboration Acciari, V.A.; Ansoldi, S.; Antonelli, L.A.; Arbet Engels, A.; Baack, D.; Babić, A.; Banerjee, B.; Barres de Almeida, U.; Barrio, J.A. et al. Teraelectronvolt emission from the γ-ray burst GRB 190114C. Nature; 2019; 575, pp. 455-458. [DOI: https://dx.doi.org/10.1038/s41586-019-1750-x]

4. Abdalla, H.; Adam, R.; Aharonian, F.; Ait Benkhali, F.; Angüner, E.O.; Arakawa, M.; Arcaro, C.; Armand, C.; Ashkar, H.; Backes, M. et al. A very-high-energy component deep in the γ-ray burst afterglow. Nature; 2019; 575, pp. 464-467. [DOI: https://dx.doi.org/10.1038/s41586-019-1743-9]

5. HESS Collaboration Abdalla, H.; Aharonian, F.; Ait Benkhali, F.; Angüner, E.O.; Arcaro, C.; Armand, C.; Armstrong, T.; Ashkar, H.; Backes, M. et al. Revealing x-ray and gamma ray temporal and spectral similarities in the GRB 190829A afterglow. Science; 2021; 372, pp. 1081-1085. [DOI: https://dx.doi.org/10.1126/science.abe8560]

6. Rees, M.J.; Meszaros, P. Unsteady Outflow Models for Cosmological Gamma-Ray Bursts. Astrophys. J.; 1994; 430, L93. [DOI: https://dx.doi.org/10.1086/187446]

7. Thompson, C. A model of gamma-ray bursts. Mon. Not. RAS; 1994; 270, pp. 480-498. [DOI: https://dx.doi.org/10.1093/mnras/270.3.480]

8. Spruit, H.C.; Daigne, F.; Drenkhahn, G. Large scale magnetic fields and their dissipation in GRB fireballs. Astron. Astrophys.; 2001; 369, pp. 694-705. [DOI: https://dx.doi.org/10.1051/0004-6361:20010131]

9. Rees, M.J.; Meszaros, P. Relativistic fireballs—Energy conversion and time-scales. Mon. Not. RAS; 1992; 258, 41. [DOI: https://dx.doi.org/10.1093/mnras/258.1.41P]

10. Chevalier, R.A.; Li, Z.Y. Gamma-Ray Burst Environments and Progenitors. Astrophys. J.; 1999; 520, pp. L29-L32. [DOI: https://dx.doi.org/10.1086/312147]

11. Woosley, S.E. Gamma-Ray Bursts from Stellar Mass Accretion Disks around Black Holes. Astrophys. J.; 1993; 405, 273. [DOI: https://dx.doi.org/10.1086/172359]

12. Paczyński, B. Are Gamma-Ray Bursts in Star-Forming Regions?. Astrophys. J.; 1998; 494, pp. L45-L48. [DOI: https://dx.doi.org/10.1086/311148]

13. Abbott, B.P.; Abbott, R.; Abbott, T.D.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.E.A. GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. Phys. Rev. Lett.; 2017; 119, 161101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.119.161101]

14. Usov, V.V. Millisecond pulsars with extremely strong magnetic fields as a cosmological source of γ-ray bursts. Nature; 1992; 357, pp. 472-474. [DOI: https://dx.doi.org/10.1038/357472a0]

15. Thompson, T.A.; Metzger, B.D.; Bucciantini, N. Proto-Magnetars as GRB Central Engines: Uncertainties, Limitations, & Particulars. Deciphering the Ancient Universe with Gamma-ray Bursts; American Institute of Physics Conference Series Kawai, N.; Nagataki, S. American Institute of Physics: College Park, MD, USA, 2010; Volume 1279, pp. 81-88. [DOI: https://dx.doi.org/10.1063/1.3509357]

16. Metzger, B.D.; Giannios, D.; Thompson, T.A.; Bucciantini, N.; Quataert, E. The protomagnetar model for gamma-ray bursts. Mon. Not. RAS; 2011; 413, pp. 2031-2056. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2011.18280.x]

17. Dall’Osso, S.; Stratta, G.; Guetta, D.; Covino, S.; De Cesare, G.; Stella, L. Gamma-ray bursts afterglows with energy injection from a spinning down neutron star. Astron. Astrophys.; 2011; 526, A121. [DOI: https://dx.doi.org/10.1051/0004-6361/201014168]

18. Giannios, D.; Spruit, H.C. Spectra of Poynting-flux powered GRB outflows. Astron. Astrophys.; 2005; 430, pp. 1-7. [DOI: https://dx.doi.org/10.1051/0004-6361:20047033]

19. Zhang, B.; Yan, H. The Internal-collision-induced Magnetic Reconnection and Turbulence (ICMART) Model of Gamma-ray Bursts. Astrophys. J.; 2011; 726, 90. [DOI: https://dx.doi.org/10.1088/0004-637X/726/2/90]

20. Lazarian, A.; Zhang, B.; Xu, S. Gamma-Ray Bursts Induced by Turbulent Reconnection. Astrophys. J.; 2019; 882, 184. [DOI: https://dx.doi.org/10.3847/1538-4357/ab2b38]

21. Lazzati, D.; Ghisellini, G.; Celotti, A. Constraints on the bulk Lorentz factor in the internal shock scenario for gamma-ray bursts. Mon. Not. RAS; 1999; 309, pp. L13-L17. [DOI: https://dx.doi.org/10.1046/j.1365-8711.1999.02970.x]

22. Sironi, L.; Petropoulou, M.; Giannios, D. Relativistic jets shine through shocks or magnetic reconnection?. Mon. Not. RAS; 2015; 450, pp. 183-191. [DOI: https://dx.doi.org/10.1093/mnras/stv641]

23. Fan, Y.; Piran, T. Gamma-ray burst efficiency and possible physical processes shaping the early afterglow. Mon. Not. RAS; 2006; 369, pp. 197-206. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2006.10280.x]

24. Beniamini, P.; Nava, L.; Piran, T. A revised analysis of gamma-ray bursts’ prompt efficiencies. Mon. Not. RAS; 2016; 461, pp. 51-59. [DOI: https://dx.doi.org/10.1093/mnras/stw1331]

25. Ghisellini, G.; Celotti, A.; Lazzati, D. Constraints on the emission mechanisms of gamma-ray bursts. Mon. Not. RAS; 2000; 313, pp. L1-L5. [DOI: https://dx.doi.org/10.1046/j.1365-8711.2000.03354.x]

26. Preece, R.D.; Briggs, M.S.; Mallozzi, R.S.; Pendleton, G.N.; Paciesas, W.S.; Band, D.L. The Synchrotron Shock Model Confronts a “Line of Death” in the BATSE Gamma-Ray Burst Data. Astrophys. J.; 1998; 506, pp. L23-L26. [DOI: https://dx.doi.org/10.1086/311644]

27. Oganesyan, G.; Nava, L.; Ghirlanda, G.; Celotti, A. Detection of Low-energy Breaks in Gamma-Ray Burst Prompt Emission Spectra. Astrophys. J.; 2017; 846, 137. [DOI: https://dx.doi.org/10.3847/1538-4357/aa831e]

28. Oganesyan, G.; Nava, L.; Ghirlanda, G.; Celotti, A. Characterization of gamma-ray burst prompt emission spectra down to soft X-rays. Astron. Astrophys.; 2018; 616, A138. [DOI: https://dx.doi.org/10.1051/0004-6361/201732172]

29. Oganesyan, G.; Nava, L.; Ghirlanda, G.; Melandri, A.; Celotti, A. Prompt optical emission as a signature of synchrotron radiation in gamma-ray bursts. Astron. Astrophys.; 2019; 628, A59. [DOI: https://dx.doi.org/10.1051/0004-6361/201935766]

30. Ravasio, M.E.; Oganesyan, G.; Ghirlanda, G.; Nava, L.; Ghisellini, G.; Pescalli, A.; Celotti, A. Consistency with synchrotron emission in the bright GRB 160625B observed by Fermi. Astron. Astrophys.; 2018; 613, A16. [DOI: https://dx.doi.org/10.1051/0004-6361/201732245]

31. Ravasio, M.E.; Ghirlanda, G.; Nava, L.; Ghisellini, G. Evidence of two spectral breaks in the prompt emission of gamma-ray bursts. Astron. Astrophys.; 2019; 625, A60. [DOI: https://dx.doi.org/10.1051/0004-6361/201834987]

32. Meszaros, P.; Rees, M.J. Relativistic Fireballs and Their Impact on External Matter: Models for Cosmological Gamma-Ray Bursts. Astrophys. J.; 1993; 405, 278. [DOI: https://dx.doi.org/10.1086/172360]

33. Mészáros, P.; Rees, M.J. Optical and Long-Wavelength Afterglow from Gamma-Ray Bursts. Astrophys. J.; 1997; 476, pp. 232-237. [DOI: https://dx.doi.org/10.1086/303625]

34. Sari, R.; Piran, T. Predictions for the Very Early Afterglow and the Optical Flash. Astrophys. J.; 1999; 520, pp. 641-649. [DOI: https://dx.doi.org/10.1086/307508]

35. Blandford, R.D.; McKee, C.F. Fluid dynamics of relativistic blast waves. Phys. Fluids; 1976; 19, pp. 1130-1138. [DOI: https://dx.doi.org/10.1063/1.861619]

36. Sedov, L.I. Propagation of strong shock waves. J. Appl. Math. Mech.; 1946; 10, pp. 241-250.

37. Taylor, G. The Formation of a Blast Wave by a Very Intense Explosion. II. The Atomic Explosion of 1945. Proc. R. Soc. Lond. Ser. A; 1950; 201, pp. 175-186. [DOI: https://dx.doi.org/10.1098/rspa.1950.0050]

38. Sari, R.; Piran, T.; Narayan, R. Spectra and Light Curves of Gamma-Ray Burst Afterglows. Astrophys. J.; 1998; 497, pp. L17-L20. [DOI: https://dx.doi.org/10.1086/311269]

39. Kobayashi, S.; Sari, R. Optical Flashes and Radio Flares in Gamma-Ray Burst Afterglow: Numerical Study. Astrophys. J.; 2000; 542, pp. 819-828. [DOI: https://dx.doi.org/10.1086/317021]

40. Vanthieghem, A.; Lemoine, M.; Plotnikov, I.; Grassi, A.; Grech, M.; Gremillet, L.; Pelletier, G. Physics and Phenomenology of Weakly Magnetized, Relativistic Astrophysical Shock Waves. Galaxies; 2020; 8, 33. [DOI: https://dx.doi.org/10.3390/galaxies8020033]

41. Sironi, L.; Spitkovsky, A. Particle Acceleration in Relativistic Magnetized Collisionless Electron-Ion Shocks. Astrophys. J.; 2011; 726, 75. [DOI: https://dx.doi.org/10.1088/0004-637X/726/2/75]

42. Sironi, L.; Spitkovsky, A.; Arons, J. The Maximum Energy of Accelerated Particles in Relativistic Collisionless Shocks. Astrophys. J.; 2013; 771, 54. [DOI: https://dx.doi.org/10.1088/0004-637X/771/1/54]

43. Keshet, U.; Katz, B.; Spitkovsky, A.; Waxman, E. Magnetic Field Evolution in Relativistic Unmagnetized Collisionless Shocks. Astrophys. J.; 2009; 693, pp. L127-L130. [DOI: https://dx.doi.org/10.1088/0004-637X/693/2/L127]

44. Kumar, P.; Zhang, B. The physics of gamma-ray bursts & relativistic jets. Phys. Rep.; 2015; 561, pp. 1-109. [DOI: https://dx.doi.org/10.1016/j.physrep.2014.09.008]

45. Zhang, B. The Physics of Gamma-Ray Bursts; Cambridge University Press: Cambridge, UK, 2018; [DOI: https://dx.doi.org/10.1017/9781139226530]

46. Nava, L. High-energy emission from gamma-ray bursts. Int. J. Mod. Phys. D; 2018; 27, 1842003. [DOI: https://dx.doi.org/10.1142/S0218271818420038]

47. Panaitescu, A.; Kumar, P. Analytic Light Curves of Gamma-Ray Burst Afterglows: Homogeneous versus Wind External Media. Astrophys. J.; 2000; 543, pp. 66-76. [DOI: https://dx.doi.org/10.1086/317090]

48. Sari, R.; Esin, A.A. On the Synchrotron Self-Compton Emission from Relativistic Shocks and Its Implications for Gamma-Ray Burst Afterglows. Astrophys. J.; 2001; 548, pp. 787-799. [DOI: https://dx.doi.org/10.1086/319003]

49. Nousek, J.A.; Kouveliotou, C.; Grupe, D.; Page, K.L.; Granot, J.; Ramirez-Ruiz, E.; Patel, S.K.; Burrows, D.N.; Mangano, V.; Barthelmy, S. et al. Evidence for a Canonical Gamma-Ray Burst Afterglow Light Curve in the Swift XRT Data. Astrophys. J.; 2006; 642, pp. 389-400. [DOI: https://dx.doi.org/10.1086/500724]

50. Nardini, M.; Ghisellini, G.; Ghirlanda, G.; Tavecchio, F.; Firmani, C.; Lazzati, D. Clustering of the optical-afterglow luminosities of long gamma-ray bursts. Astron. Astrophys.; 2006; 451, pp. 821-833. [DOI: https://dx.doi.org/10.1051/0004-6361:20054085]

51. Zaninoni, E.; Bernardini, M.G.; Margutti, R.; Oates, S.; Chincarini, G. Gamma-ray burst optical light-curve zoo: Comparison with X-ray observations. Astron. Astrophys.; 2013; 557, A12. [DOI: https://dx.doi.org/10.1051/0004-6361/201321221]

52. Bernardini, M.G.; Margutti, R.; Mao, J.; Zaninoni, E.; Chincarini, G. The X-ray light curve of gamma-ray bursts: Clues to the central engine. Astron. Astrophys.; 2012; 539, A3. [DOI: https://dx.doi.org/10.1051/0004-6361/201117895]

53. Sari, R.; Piran, T. Hydrodynamic Timescales and Temporal Structure of Gamma-Ray Bursts. Astrophys. J.; 1995; 455, L143. [DOI: https://dx.doi.org/10.1086/309835]

54. Granot, J.; Sari, R. The Shape of Spectral Breaks in Gamma-Ray Burst Afterglows. Astrophys. J.; 2002; 568, pp. 820-829. [DOI: https://dx.doi.org/10.1086/338966]

55. Nava, L.; Sironi, L.; Ghisellini, G.; Celotti, A.; Ghirlanda, G. Afterglow emission in gamma-ray bursts - I. Pair-enriched ambient medium and radiative blast waves. Mon. Not. RAS; 2013; 433, pp. 2107-2121. [DOI: https://dx.doi.org/10.1093/mnras/stt872]

56. Zhang, B.; Mészáros, P. Gamma-Ray Burst Afterglow with Continuous Energy Injection: Signature of a Highly Magnetized Millisecond Pulsar. Astrophys. J.; 2001; 552, pp. L35-L38. [DOI: https://dx.doi.org/10.1086/320255]

57. Zhang, B.; Fan, Y.Z.; Dyks, J.; Kobayashi, S.; Mészáros, P.; Burrows, D.N.; Nousek, J.A.; Gehrels, N. Physical Processes Shaping Gamma-Ray Burst X-Ray Afterglow Light Curves: Theoretical Implications from the Swift X-Ray Telescope Observations. Astrophys. J.; 2006; 642, pp. 354-370. [DOI: https://dx.doi.org/10.1086/500723]

58. Laskar, T.; Berger, E.; Margutti, R.; Perley, D.; Zauderer, B.A.; Sari, R.; Fong, W.F. Energy Injection in Gamma-Ray Burst Afterglows. Astrophys. J.; 2015; 814, 1. [DOI: https://dx.doi.org/10.1088/0004-637X/814/1/1]

59. Shapiro, P.R. Relativistic blast waves that accelerate. Astrophys. J.; 1980; 236, pp. 958-969. [DOI: https://dx.doi.org/10.1086/157823]

60. Sironi, L.; Keshet, U.; Lemoine, M. Relativistic Shocks: Particle Acceleration and Magnetization. Space Sci. Rev.; 2015; 191, pp. 519-544. [DOI: https://dx.doi.org/10.1007/s11214-015-0181-8]

61. Keshet, U.; Waxman, E. Energy Spectrum of Particles Accelerated in Relativistic Collisionless Shocks. Phys. Rev. Lett.; 2005; 94, 111102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.94.111102]

62. Weibel, E.S. Spontaneously Growing Transverse Waves in a Plasma Due to an Anisotropic Velocity Distribution. Phys. Rev. Lett.; 1959; 2, pp. 83-84. [DOI: https://dx.doi.org/10.1103/PhysRevLett.2.83]

63. Spitkovsky, A. Particle Acceleration in Relativistic Collisionless Shocks: Fermi Process at Last?. Astrophys. J.; 2008; 682, L5. [DOI: https://dx.doi.org/10.1086/590248]

64. Böhm, D. The Characteristics of Electrical Discharges in Magnetic Fields; Guthrie, A.; Wakerling, R.K. McGraw-Hill: New York, NY, USA, 1949.

65. Lemoine, M. Synchrotron signature of a relativistic blast wave with decaying microturbulence. Mon. Not. RAS; 2013; 428, pp. 845-866. [DOI: https://dx.doi.org/10.1093/mnras/sts081]

66. Lemoine, M.; Li, Z.; Wang, X.Y. On the magnetization of gamma-ray burst blast waves. Mon. Not. RAS; 2013; 435, pp. 3009-3016. [DOI: https://dx.doi.org/10.1093/mnras/stt1494]

67. Santana, R.; Barniol Duran, R.; Kumar, P. Magnetic Fields in Relativistic Collisionless Shocks. Astrophys. J.; 2014; 785, 29. [DOI: https://dx.doi.org/10.1088/0004-637X/785/1/29]

68. Wang, X.G.; Zhang, B.; Liang, E.W.; Gao, H.; Li, L.; Deng, C.M.; Qin, S.M.; Tang, Q.W.; Kann, D.A.; Ryde, F. et al. How Bad or Good Are the External Forward Shock Afterglow Models of Gamma-Ray Bursts?. Astrophys. J. Lett.; 2015; 219, 9. [DOI: https://dx.doi.org/10.1088/0067-0049/219/1/9]

69. Zhang, B.B.; van Eerten, H.; Burrows, D.N.; Ryan, G.S.; Evans, P.A.; Racusin, J.L.; Troja, E.; MacFadyen, A. An Analysis of Chandra Deep Follow-up Gamma-Ray Bursts: Implications for Off-axis Jets. Astrophys. J.; 2015; 806, 15. [DOI: https://dx.doi.org/10.1088/0004-637X/806/1/15]

70. Beniamini, P.; Nava, L.; Duran, R.B.; Piran, T. Energies of GRB blast waves and prompt efficiencies as implied by modelling of X-ray and GeV afterglows. Mon. Not. RAS; 2015; 454, pp. 1073-1085. [DOI: https://dx.doi.org/10.1093/mnras/stv2033]

71. MAGIC Collaboration Veres, P.; Bhat., P.N.; Briggs, M.S.; Cleveland, W.H.; Hamburg, R.; Hui, C.M.; Mailyan, B.; Preece, R.D.; Roberts, O.J. et al. Observation of inverse Compton emission from a long γ-ray burst. Nature; 2019; 575, pp. 459-463. [DOI: https://dx.doi.org/10.1038/s41586-019-1754-6]

72. Salafia, O.S.; Ravasio, M.E.; Yang, J.; An, T.; Orienti, M.; Ghirlanda, G.; Nava, L.; Giroletti, M.; Mohan, P.; Spinelli, R. et al. Multi-wavelength view of the close-by GRB~190829A sheds light on gamma-ray burst physics. arXiv; 2021; arXiv: 2106.07169

73. Eichler, D.; Waxman, E. The Efficiency of Electron Acceleration in Collisionless Shocks and Gamma-Ray Burst Energetics. Astrophys. J.; 2005; 627, pp. 861-867. [DOI: https://dx.doi.org/10.1086/430596]

74. Ackermann, M.; Ajello, M.; Asano, K.; Atwood, W.B.; Axelsson, M.; Baldini, L.; Ballet, J.; Barbiellini, G.; Baring, M.G.; Bastieri, D. et al. Fermi-LAT Observations of the Gamma-Ray Burst GRB 130427A. Science; 2014; 343, pp. 42-47. [DOI: https://dx.doi.org/10.1126/science.1242353] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24263133]

75. Nakar, E.; Ando, S.; Sari, R. Klein-Nishina Effects on Optically Thin Synchrotron and Synchrotron Self-Compton Spectrum. Astrophys. J.; 2009; 703, pp. 675-691. [DOI: https://dx.doi.org/10.1088/0004-637X/703/1/675]

76. Petropoulou, M.; Mastichiadis, A. On the multiwavelength emission from gamma ray burst afterglows. Astron. Astrophys.; 2009; 507, pp. 599-610. [DOI: https://dx.doi.org/10.1051/0004-6361/200912970]

77. Pennanen, T.; Vurm, I.; Poutanen, J. Simulations of gamma-ray burst afterglows with a relativistic kinetic code. Astron. Astrophys.; 2014; 564, A77. [DOI: https://dx.doi.org/10.1051/0004-6361/201322520]

78. Bošnjak, Ž.; Daigne, F.; Dubus, G. Prompt high-energy emission from gamma-ray bursts in the internal shock model. Astron. Astrophys.; 2009; 498, pp. 677-703. [DOI: https://dx.doi.org/10.1051/0004-6361/200811375]

79. Chang, J.S.; Cooper, G. A Practical Difference Scheme for Fokker-Planck Equations. J. Comput. Phys.; 1970; 6, pp. 1-16. [DOI: https://dx.doi.org/10.1016/0021-9991(70)90001-X]

80. Rybicki, G.B.; Lightman, A.P. Radiative Processes in Astrophysics; Wiley: Hoboken, NJ, USA, 1979.

81. Jones, F.C. Calculated Spectrum of Inverse-Compton-Scattered Photons. Phys. Rev.; 1968; 167, pp. 1159-1169. [DOI: https://dx.doi.org/10.1103/PhysRev.167.1159]

82. Ghisellini, G.; Guilbert, P.W.; Svensson, R. The Synchrotron Boiler. Astrophys. J.; 1988; 334, L5. [DOI: https://dx.doi.org/10.1086/185300]

83. Coppi, P.S.; Blandford, R.D. Reaction rates and energy distributions for elementary processes in relativistic pair plasmas. Mon. Not. RAS; 1990; 245, pp. 453-507. [DOI: https://dx.doi.org/10.1093/mnras/245.3.453]

84. Chincarini, G.; Mao, J.; Margutti, R.; Bernardini, M.G.; Guidorzi, C.; Pasotti, F.; Giannios, D.; Della Valle, M.; Moretti, A.; Romano, P. et al. Unveiling the origin of X-ray flares in gamma-ray bursts. Mon. Not. RAS; 2010; 406, pp. 2113-2148. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2010.17037.x]

85. Margutti, R.; Bernardini, G.; Barniol Duran, R.; Guidorzi, C.; Shen, R.F.; Chincarini, G. On the average gamma-ray burst X-ray flaring activity. Mon. Not. RAS; 2011; 410, pp. 1064-1075. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2010.17504.x]

86. Bernardini, M.G.; Margutti, R.; Chincarini, G.; Guidorzi, C.; Mao, J. Gamma-ray burst long lasting X-ray flaring activity. Astron. Astrophys.; 2011; 526, A27. [DOI: https://dx.doi.org/10.1051/0004-6361/201015703]

87. Yi, S.X.; Xi, S.Q.; Yu, H.; Wang, F.Y.; Mu, H.J.; Lü, L.Z.; Liang, E.W. Comprehensive Study of the X-Ray Flares from Gamma-ray Bursts Observed by Swift. Astrophys. J. Lett.; 2016; 224, 20. [DOI: https://dx.doi.org/10.3847/0067-0049/224/2/20]

88. Beniamini, P.; Kumar, P. X-ray flares in GRBs: General considerations and photospheric origin. Mon. Not. RAS; 2016; 457, pp. L108-L112. [DOI: https://dx.doi.org/10.1093/mnrasl/slw003]

89. Yi, S.X.; Yu, H.; Wang, F.Y.; Dai, Z.G. Statistical Distributions of Optical Flares from Gamma-Ray Bursts. Astrophys. J.; 2017; 844, 79. [DOI: https://dx.doi.org/10.3847/1538-4357/aa7b7b]

90. Burrows, D.N.; Romano, P.; Falcone, A.; Kobayashi, S.; Zhang, B.; Moretti, A.; O’Brien, P.T.; Goad, M.R.; Campana, S.; Page, K.L. et al. Bright X-ray Flares in Gamma-Ray Burst Afterglows. Science; 2005; 309, pp. 1833-1835. [DOI: https://dx.doi.org/10.1126/science.1116168]

91. Willingale, R.; Genet, F.; Granot, J.; O’Brien, P.T. The spectral-temporal properties of the prompt pulses and rapid decay phase of gamma-ray bursts. Mon. Not. RAS; 2010; 403, pp. 1296-1316. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2009.16187.x]

92. Fan, Y.Z.; Wei, D.M. Late internal-shock model for bright X-ray flares in gamma-ray burst afterglows and GRB 011121. Mon. Not. RAS; 2005; 364, pp. L42-L46. [DOI: https://dx.doi.org/10.1111/j.1745-3933.2005.00102.x]

93. Lazzati, D.; Perna, R. X-ray flares and the duration of engine activity in gamma-ray bursts. Mon. Not. RAS; 2007; 375, pp. L46-L50. [DOI: https://dx.doi.org/10.1111/j.1745-3933.2006.00273.x]

94. Maxham, A.; Zhang, B. Modeling Gamma-Ray Burst X-Ray Flares Within the Internal Shock Model. Astrophys. J.; 2009; 707, pp. 1623-1633. [DOI: https://dx.doi.org/10.1088/0004-637X/707/2/1623]

95. Wang, X.Y.; Li, Z.; Mészáros, P. GeV-TeV and X-Ray Flares from Gamma-Ray Bursts. Astrophys. J.; 2006; 641, pp. L89-L92. [DOI: https://dx.doi.org/10.1086/504151]

96. Ajello, M.; Arimoto, M.; Asano, K.; Axelsson, M.; Baldini, L.; Barbiellini, G.; Bastieri, D.; Bellazzini, R.; Berretta, A.; Bissaldi, E. et al. Bright Gamma-Ray Flares Observed in GRB 131108A. Astrophys. J.; 2019; 886, L33. [DOI: https://dx.doi.org/10.3847/2041-8213/ab564f]

97. Wang, K.; Dai, Z.G. GeV Emission during X-Ray Flares from Late Internal Shocks: Application to GRB 100728A. Astrophys. J.; 2013; 772, 152. [DOI: https://dx.doi.org/10.1088/0004-637X/772/2/152]

98. He, H.N.; Zhang, B.B.; Wang, X.Y.; Li, Z.; Mészáros, P. Origin of the GeV Emission during the X-Ray Flaring Activity in GRB 100728A. Astrophys. J.; 2012; 753, 178. [DOI: https://dx.doi.org/10.1088/0004-637X/753/2/178]

99. Berti, A. MAGIC View of Gamma-ray Bursts at Very High Energies, Talk at 16th Marcel Grossmann meeting, online, July 2021. Available online: https://indico.icranet.org/event/1/contributions/450/attachments/112/126/berti_MAGIC_GRBs_MG16.pdf (accessed on 27 February 2022).

100. Fioretti, V.; Ribeiro, D.; Humensky, T.B.; Bulgarelli, A.; Maier, G.; Moralejo, A.; Nigro, C. The Cherenkov Telescope Array sensitivity to the transient sky. Proceedings of the 36th International Cosmic Ray Conference (ICRC2019); Madison, WI, USA, 24 July–1 August 2019; Volume 36, 673.

101. Woosley, S.E. Models for Gamma-Ray Burst Progenitors and Central Engines. arXiv; 2011; arXiv: 1105.4193

102. Fryer, C.L.; Rockefeller, G.; Young, P.A. The Environments around Long-Duration Gamma-Ray Burst Progenitors. Astrophys. J.; 2006; 647, pp. 1269-1285. [DOI: https://dx.doi.org/10.1086/505590]

103. Panaitescu, A.; Kumar, P. Properties of Relativistic Jets in Gamma-Ray Burst Afterglows. Astrophys. J.; 2002; 571, pp. 779-789. [DOI: https://dx.doi.org/10.1086/340094]

104. Schulze, S.; Klose, S.; Björnsson, G.; Jakobsson, P.; Kann, D.A.; Rossi, A.; Krühler, T.; Greiner, J.; Ferrero, P. The circumburst density profile around GRB progenitors: A statistical study. Astron. Astrophys.; 2011; 526, A23. [DOI: https://dx.doi.org/10.1051/0004-6361/201015581]

105. Li, L.; Wu, X.F.; Huang, Y.F.; Wang, X.G.; Tang, Q.W.; Liang, Y.F.; Zhang, B.B.; Wang, Y.; Geng, J.J.; Liang, E.W. et al. A Correlated Study of Optical and X-Ray Afterglows of GRBs. Astrophys. J.; 2015; 805, 13. [DOI: https://dx.doi.org/10.1088/0004-637X/805/1/13]

106. Gompertz, B.P.; Fruchter, A.S.; Pe’er, A. The Environments of the Most Energetic Gamma-Ray Bursts. Astrophys. J.; 2018; 866, 162. [DOI: https://dx.doi.org/10.3847/1538-4357/aadba8]

107. van Marle, A.J.; Langer, N.; Achterberg, A.; García-Segura, G. Forming a constant density medium close to long gamma-ray bursts. Astron. Astrophys.; 2006; 460, pp. 105-116. [DOI: https://dx.doi.org/10.1051/0004-6361:20065709]

108. Srinivasaragavan, G.P.; Dainotti, M.G.; Fraija, N.; Hernandez, X.; Nagataki, S.; Lenart, A.; Bowden, L.; Wagner, R. On the Investigation of the Closure Relations for Gamma-Ray Bursts Observed by Swift in the Post-plateau Phase and the GRB Fundamental Plane. Astrophys. J.; 2020; 903, 18. [DOI: https://dx.doi.org/10.3847/1538-4357/abb702]

109. Kumar, P.; Barniol Duran, R. On the generation of high-energy photons detected by the Fermi Satellite from gamma-ray bursts. Mon. Not. RAS; 2009; 400, pp. L75-L79. [DOI: https://dx.doi.org/10.1111/j.1745-3933.2009.00766.x]

110. Kumar, P.; Barniol Duran, R. External forward shock origin of high-energy emission for three gamma-ray bursts detected by Fermi. Mon. Not. RAS; 2010; 409, pp. 226-236. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2010.17274.x]

111. Barniol Duran, R.; Kumar, P. Evidence for mild deviation from power-law distribution of electrons in relativistic shocks: GRB 090902B. Mon. Not. RAS; 2011; 417, pp. 1584-1600. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2011.19369.x]

112. He, H.N.; Wu, X.F.; Toma, K.; Wang, X.Y.; Mészáros, P. On the High-energy Emission of the Short GRB 090510. Astrophys. J.; 2011; 733, 22. [DOI: https://dx.doi.org/10.1088/0004-637X/733/1/22]

113. Barniol Duran, R. Constraining the magnetic field in GRB relativistic collisionless shocks using radio data. Mon. Not. RAS; 2014; 442, pp. 3147-3154. [DOI: https://dx.doi.org/10.1093/mnras/stu1070]

114. Piran, T.; Nakar, E. On the External Shock Synchrotron Model for Gamma-ray Bursts’ GeV Emission. Astrophys. J.; 2010; 718, pp. L63-L67. [DOI: https://dx.doi.org/10.1088/2041-8205/718/2/L63]

115. Filgas, R.; Greiner, J.; Schady, P.; Krühler, T.; Updike, A.C.; Klose, S.; Nardini, M.; Kann, D.A.; Rossi, A.; Sudilovsky, V. et al. GRB 091127: The cooling break race on magnetic fuel. Astron. Astrophys.; 2011; 535, A57. [DOI: https://dx.doi.org/10.1051/0004-6361/201117695]

116. van der Horst, A.J.; Paragi, Z.; de Bruyn, A.G.; Granot, J.; Kouveliotou, C.; Wiersema, K.; Starling, R.L.C.; Curran, P.A.; Wijers, R.A.M.J.; Rowlinson, A. et al. A comprehensive radio view of the extremely bright gamma-ray burst 130427A. Mon. Not. RAS; 2014; 444, pp. 3151-3163. [DOI: https://dx.doi.org/10.1093/mnras/stu1664]

117. Ioka, K.; Toma, K.; Yamazaki, R.; Nakamura, T. Efficiency crisis of swift gamma-ray bursts with shallow X-ray afterglows: Prior activity or time-dependent microphysics?. Astron. Astrophys.; 2006; 458, pp. 7-12. [DOI: https://dx.doi.org/10.1051/0004-6361:20064939]

118. Kong, S.W.; Wong, A.Y.L.; Huang, Y.F.; Cheng, K.S. Variation of microphysics in wind bubbles: An alternative mechanism for explaining the rebrightenings in Gamma-ray burst afterglows. Mon. Not. RAS; 2010; 402, pp. 409-416. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2009.15886.x]

119. Huang, L.Y.; Wang, X.G.; Zheng, W.; Liang, E.W.; Lin, D.b.; Zhong, S.Q.; Zhang, H.M.; Huang, X.L.; Filippenko, A.V.; Zhang, B. GRB 120729A: External Shock Origin for Both the Prompt Gamma-Ray Emission and Afterglow. Astrophys. J.; 2018; 859, 163. [DOI: https://dx.doi.org/10.3847/1538-4357/aaba6e]

120. Panaitescu, A. Swift gamma-ray burst afterglows and the forward-shock model. Mon. Not. RAS; 2007; 379, pp. 331-342. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2007.11956.x]

121. Misra, K.; Resmi, L.; Kann, D.A.; Marongiu, M.; Moin, A.; Klose, S.; Bernardi, G.; de Ugarte Postigo, A.; Jaiswal, V.K.; Schulze, S. et al. Low frequency view of GRB 190114C reveals time varying shock micro-physics. Mon. Not. RAS; 2021; 504, pp. 5685-5701. [DOI: https://dx.doi.org/10.1093/mnras/stab1050]

122. Guilbert, P.W.; Fabian, A.C.; Rees, M.J. Spectral and variability constraints on compact sources. Mon. Not. RAS; 1983; 205, pp. 593-603. [DOI: https://dx.doi.org/10.1093/mnras/205.3.593]

123. de Jager, O.C.; Harding, A.K. The Expected High-Energy to Ultra–High-Energy Gamma-Ray Spectrum of the Crab Nebula. Astrophys. J.; 1992; 396, 161. [DOI: https://dx.doi.org/10.1086/171706]

124. Ajello, M.E.A. A Decade of Gamma-Ray Bursts Observed by Fermi-LAT: The Second GRB Catalog. Astrophys. J.; 2019; 878, 52. [DOI: https://dx.doi.org/10.3847/1538-4357/ab1d4e]

125. Kumar, P.; Hernández, R.A.; Bošnjak, Ž.; Barniol Duran, R. Maximum synchrotron frequency for shock-accelerated particles. Mon. Not. RAS; 2012; 427, pp. L40-L44. [DOI: https://dx.doi.org/10.1111/j.1745-3933.2012.01341.x]

126. Medvedev, M.V.; Loeb, A. Generation of Magnetic Fields in the Relativistic Shock of Gamma-Ray Burst Sources. Astrophys. J.; 1999; 526, pp. 697-706. [DOI: https://dx.doi.org/10.1086/308038]

127. Kumar, P. The Distribution of Burst Energy and Shock Parameters for Gamma-Ray Bursts. Astrophys. J.; 2000; 538, pp. L125-L128. [DOI: https://dx.doi.org/10.1086/312821]

128. Berger, E.; Kulkarni, S.R.; Frail, D.A. A Standard Kinetic Energy Reservoir in Gamma-Ray Burst Afterglows. Astrophys. J.; 2003; 590, pp. 379-385. [DOI: https://dx.doi.org/10.1086/374892]

129. Lloyd-Ronning, N.M.; Zhang, B. On the Kinetic Energy and Radiative Efficiency of Gamma-Ray Bursts. Astrophys. J.; 2004; 613, pp. 477-483. [DOI: https://dx.doi.org/10.1086/423026]

130. Berger, E. The Prompt Gamma-Ray and Afterglow Energies of Short-Duration Gamma-Ray Bursts. Astrophys. J.; 2007; 670, pp. 1254-1259. [DOI: https://dx.doi.org/10.1086/522195]

131. D’Avanzo, P.; Salvaterra, R.; Sbarufatti, B.; Nava, L.; Melandri, A.; Bernardini, M.G.; Campana, S.; Covino, S.; Fugazza, D.; Ghirlanda, G. et al. A complete sample of bright Swift Gamma-ray bursts: X-ray afterglow luminosity and its correlation with the prompt emission. Mon. Not. RAS; 2012; 425, pp. 506-513. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2012.21489.x]

132. Granot, J.; Königl, A.; Piran, T. Implications of the early X-ray afterglow light curves of Swift gamma-ray bursts. Mon. Not. RAS; 2006; 370, pp. 1946-1960. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2006.10621.x]

133. Zhang, B.; Liang, E.; Page, K.L.; Grupe, D.; Zhang, B.B.; Barthelmy, S.D.; Burrows, D.N.; Campana, S.; Chincarini, G.; Gehrels, N. et al. GRB Radiative Efficiencies Derived from the Swift Data: GRBs versus XRFs, Long versus Short. Astrophys. J.; 2007; 655, pp. 989-1001. [DOI: https://dx.doi.org/10.1086/510110]

134. Nava, L.; Vianello, G.; Omodei, N.; Ghisellini, G.; Ghirlanda, G.; Celotti, A.; Longo, F.; Desiante, R.; Barniol Duran, R. Clustering of LAT light curves: A clue to the origin of high-energy emission in gamma-ray bursts. Mon. Not. RAS; 2014; 443, pp. 3578-3585. [DOI: https://dx.doi.org/10.1093/mnras/stu1451]

135. Beniamini, P.; Granot, J. Properties of GRB light curves from magnetic reconnection. Mon. Not. RAS; 2016; 459, pp. 3635-3658. [DOI: https://dx.doi.org/10.1093/mnras/stw895]

136. Ressler, S.M.; Laskar, T. Thermal Electrons in Gamma-Ray Burst Afterglows. Astrophys. J.; 2017; 845, 150. [DOI: https://dx.doi.org/10.3847/1538-4357/aa8268]

137. Warren, D.C.; Dainotti, M.; Barkov, M.V.; Ahlgren, B.; Ito, H.; Nagataki, S. A Semianalytic Afterglow with Thermal Electrons and Synchrotron Self-Compton Emission. Astrophys. J.; 2022; 924, 40. [DOI: https://dx.doi.org/10.3847/1538-4357/ac2f43]

138. Giannios, D.; Spitkovsky, A. Signatures of a Maxwellian component in shock-accelerated electrons in GRBs. Mon. Not. RAS; 2009; 400, pp. 330-336. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2009.15454.x]

139. Mirzoyan, R.; Noda, K.; Moretti, E.; Berti, A.; Nigro, C.; Hoang, J.; Micanovic, S.; Takahashi, M.; Chai, Y.; Moralejo, A. et al. MAGIC detects the GRB 190114C in the TeV energy domain. GRB Coord. Netw.; 2019; 23701, 1.

140. Li, T.P.; Ma, Y.Q. Analysis methods for results in gamma-ray astronomy. Astrophys. J.; 1983; 272, pp. 317-324. [DOI: https://dx.doi.org/10.1086/161295]

141. Siegel, M.H.; Barthelmy, S.D.; Burrows, D.N.; Lien, A.Y.; Marshall, F.E.; Palmer, D.M.; Sbarufatti, B. GRB 160821B: Swift detection of a short burst. GRB Coord. Netw.; 2016; 19833, 1.

142. Stanbro, M.; Meegan, C. GRB 160821B: Fermi GBM Detection. GRB Coord. Netw.; 2016; 19843, 1.

143. Palmer, D.M.; Barthelmy, S.D.; Cummings, J.R.; Gehrels, N.; Krimm, H.A.; Lien, A.Y.; Markwardt, C.B.; Norris, J.P.; Sakamoto, T.; Siegel, M.H. et al. GRB 160821B: Swift-BAT refined analysis. GRB Coord. Netw.; 2016; 19844, 1.

144. Berger, E. Short-Duration Gamma-Ray Bursts. Annu. Rev. Astron Astrophys.; 2014; 52, pp. 43-105. [DOI: https://dx.doi.org/10.1146/annurev-astro-081913-035926]

145. Acciari, V.A.; Ansoldi, S.; Antonelli, L.A.; Arbet Engels, A.; Asano, K.; Baack, D.; Babić, A.; Baquero, A.; Barres de Almeida, U.; Barrio, J.A. et al. MAGIC Observations of the Nearby Short Gamma-Ray Burst GRB 160821B. Astrophys. J.; 2021; 908, 90. [DOI: https://dx.doi.org/10.3847/1538-4357/abd249]

146. Evans, P.A.; Goad, M.R.; Osborne, J.P.; Beardmore, A.P. GRB 160821B: Enhanced Swift-XRT position. GRB Coord. Netw.; 2016; 19837, 1.

147. Sbarufatti, B.; Burrows, D.N.; Osborne, J.P.; Page, K.L.; Mingo, B.; Maselli, A.; Melandri, A.; D’Avanzo, P.; Roegiers, T.G.R.; Siegel, M.H. GRB 160821B: Swift-XRT refined Analysis. GRB Coord. Netw.; 2016; 19841, 1.

148. Troja, E.; Castro-Tirado, A.J.; Becerra González, J.; Hu, Y.; Ryan, G.S.; Cenko, S.B.; Ricci, R.; Novara, G.; Sánchez-Rámirez, R.; Acosta-Pulido, J.A. et al. The afterglow and kilonova of the short GRB 160821B. Mon. Not. RAS; 2019; 489, pp. 2104-2116. [DOI: https://dx.doi.org/10.1093/mnras/stz2255]

149. Jeong, S.; Park, I.H.; Hu, Y.; Scarpa, R.; Castro-Tirado, A.J. GRB 160821B: GTC follow-up observation. GRB Coord. Netw.; 2016; 19847, 1.

150. Troja, E.; Tanvir, N.; Cenko, S.B.; Levan, A.; Barnes, J.; Castro-Tirado, A.; Fruchter, A.S.; Gehrels, N.; Greiner, J.; Kawai, N. et al. GRB 160821B: HST detection of the optical and IR counterpart. GRB Coord. Netw.; 2016; 20222, 1.

151. Xu, D.; Malesani, D.; de Ugarte Postigo, A.; Gafton, E.; Rivero Losada, I. GRB 160821B: NOT optical afterglow candidate. GRB Coord. Netw.; 2016; 19834, 1.

152. Levan, A.J.; Wiersema, K.; Tanvir, N.R.; Malesani, D.; Xu, D.; de Ugarte Postigo, A. GRB 160821B: WHT afterglow confirmation and redshift of candidate host. GRB Coord. Netw.; 2016; 19846, 1.

153. Lamb, G.P.; Tanvir, N.R.; Levan, A.J.; de Ugarte Postigo, A.; Kawaguchi, K.; Corsi, A.; Evans, P.A.; Gompertz, B.; Malesani, D.B.; Page, K.L. et al. Short GRB 160821B: A Reverse Shock, a Refreshed Shock, and a Well-sampled Kilonova. Astrophys. J.; 2019; 883, 48. [DOI: https://dx.doi.org/10.3847/1538-4357/ab38bb]

154. Fong, W.; Alexander, K.D.; Laskar, T. GRB 160821B: VLA 5.0 GHz observations. GRB Coord. Netw.; 2016; 19854, 1.

155. Fruck, C.; Gaug, M.; Zanin, R.; Dorner, D.; Garrido, D.; Mirzoyan, R.; Font, L. A novel LIDAR-based Atmospheric Calibration Method for Improving the Data Analysis of MAGIC. arXiv; 2014; arXiv: 1403.3591

156. Rolke, W.A.; López, A.M.; Conrad, J. Limits and confidence intervals in the presence of nuisance parameters. Nucl. Instrum. Methods Phys. Res. A; 2005; 551, pp. 493-503. [DOI: https://dx.doi.org/10.1016/j.nima.2005.05.068]

157. Zhang, B.T.; Murase, K.; Yuan, C.; Kimura, S.S.; Mészáros, P. External Inverse-Compton Emission Associated with Extended and Plateau Emission of Short Gamma-Ray Bursts: Application to GRB 160821B. Astrophys. J.; 2021; 908, L36. [DOI: https://dx.doi.org/10.3847/2041-8213/abe0b0]

158. Roberts, O.J.; Meegan, C. GRB 180720B: Fermi GBM observation. GRB Coord. Netw.; 2018; 22981, 1.

159. Siegel, M.H.; Burrows, D.N.; Deich, A.; Gropp, J.D.; Kennea, J.A.; Laporte, S.J.; Lien, A.Y.; Moss, M.J.; Page, K.L.; Palmer, D.M. et al. Correction: Swift trigger #848890 is GRB 180720B (not GRB 180720A). GRB Coord. Netw.; 2018; 22975, 1.

160. Bissaldi, E.; Racusin, J.L. GRB 180720B: Fermi-LAT detection. GRB Coord. Netw.; 2018; 22980, 1.

161. Sasada, M.; Nakaoka, T.; Kawabata, M.; Uchida, N.; Yamazaki, Y.; Kawabata, K.S. GRB 180720B: Kanata 1.5m optical/NIR observation. GRB Coord. Netw.; 2018; 22977, 1.

162. Martone, R.; Guidorzi, C.; Kobayashi, S.; Mundell, C.G.; Gomboc, A.; Steele, I.A.; Cucchiara, A.; Morris, D. GRB 180720B: LCO Haleakala possible bright optical candidate. GRB Coord. Netw.; 2018; 22976, 1.

163. Reva, I.; Pozanenko, A.; Volnova, A.; Mazaeva, E.; Kusakin, A.; Krugov, M. GRB 180720B: TSHAO optical observations. GRB Coord. Netw.; 2018; 22979, 1.

164. Watson, A.M.; Butler, N.; Becerra, R.L.; Lee, W.H.; Roman-Zuniga, C.; Kutyrev, A.; Troja, E. GRB 180720B: COATLI Optical Detection. GRB Coord. Netw.; 2018; 23017, 1.

165. Schmalz, S.; Graziani, F.; Pozanenko, A.; Volnova, A.; Mazaeva, E.; Molotov, I. GRB 180720B: ISON-Castelgrande optical observations. GRB Coord. Netw.; 2018; 23020, 1.

166. Zheng, W.; Filippenko, A.V. GRB 180720B: KAIT Optical Observations. GRB Coord. Netw.; 2018; 23033, 1.

167. Lipunov, V.; Gorbovskoy, E.; Tiurina, N.; Vlasenko, D.; Kornilov, V.; Kuznetsov, A.; Chazov, V.; Gorbunov, I.; Zimnukhov, D.; Kuvshinov, D. et al. GRB 180720B: MASTER Global Net OT observations. GRB Coord. Netw.; 2018; 23023, 1.

168. Itoh, R.; Murata, K.L.; Tachibana, Y.; Harita, S.; Morita, K.; Shiraishi, K.; Iida, K.; Oeda, M.; Adachi, R.; Niwano, S. et al. GRB 180720B: MITSuME Akeno optical observations. GRB Coord. Netw.; 2018; 22983, 1.

169. Izzo, L.; Kann, D.A.; de Ugarte Postigo, A.; Thoene, C.C.; Bensch, K.; Blazek, M.; Diaz-Martin, M.C.; Rodriguez-Llano, S. GRB 180720B: OAJ optical observations. GRB Coord. Netw.; 2018; 23040, 1.

170. Covino, S.; Fugazza, D. GRB 180720B: REM photometry. GRB Coord. Netw.; 2018; 23021, 1.

171. Jelinek, M.; Strobl, J.; Hudec, R.; Polasek, C. GRB 180720B: D50 optical observations. GRB Coord. Netw.; 2018; 23024, 1.

172. Kann, D.A.; Izzo, L.; Casanova, V. GRB 180720B: OSN detection, fading. GRB Coord. Netw.; 2018; 22985, 1.

173. Rhodes, L.; van der Horst, A.J.; Fender, R.; Monageng, I.M.; Anderson, G.E.; Antoniadis, J.; Bietenholz, M.F.; Böttcher, M.; Bright, J.S.; Green, D.A. et al. Radio afterglows of very high-energy gamma-ray bursts 190829A and 180720B. Mon. Not. RAS; 2020; 496, pp. 3326-3335. [DOI: https://dx.doi.org/10.1093/mnras/staa1715]

174. Wang, X.Y.; Liu, R.Y.; Zhang, H.M.; Xi, S.Q.; Zhang, B. Synchrotron Self-Compton Emission from External Shocks as the Origin of the Sub-TeV Emission in GRB 180720B and GRB 190114C. Astrophys. J.; 2019; 884, 117. [DOI: https://dx.doi.org/10.3847/1538-4357/ab426c]

175. Liu, R.Y.; Wang, X.Y.; Wu, X.F. Interpretation of the Unprecedentedly Long-lived High-energy Emission of GRB 130427A. Astrophys. J.; 2013; 773, L20. [DOI: https://dx.doi.org/10.1088/2041-8205/773/2/L20]

176. Gropp, J.D.; Kennea, J.A.; Klingler, N.J.; Krimm, H.A.; Laporte, S.J.; Lien, A.Y.; Moss, M.J.; Palmer, D.M.; Sbarufatti, B.; Siegel, M.H. GRB 190114C: Swift detection of a very bright burst with a bright optical counterpart. GRB Coord. Netw.; 2019; 23688, 1.

177. Hamburg, R.; Veres, P.; Meegan, C.; Burns, E.; Connaughton, V.; Goldstein, A.; Kocevski, D.; Roberts, O.J. GRB 190114C: Fermi GBM detection. GRB Coord. Netw.; 2019; 23707, 1.

178. Ravasio, M.E.; Oganesyan, G.; Salafia, O.S.; Ghirland a, G.; Ghisellini, G.; Branchesi, M.; Campana, S.; Covino, S.; Salvaterra, R. GRB 190114C: From prompt to afterglow?. Astron. Astrophys.; 2019; 626, A12. [DOI: https://dx.doi.org/10.1051/0004-6361/201935214]

179. Ajello, M.e.a. Fermi and Swift Observations of GRB 190114C: Tracing the Evolution of High-energy Emission from Prompt to Afterglow. Astrophys. J.; 2020; 890, 9. [DOI: https://dx.doi.org/10.3847/1538-4357/ab5b05]

180. Laskar, T.; Alexander, K.D.; Gill, R.; Granot, J.; Berger, E.; Mundell, C.G.; Barniol Duran, R.; Bolmer, J.; Duffell, P.; van Eerten, H. et al. ALMA Detection of a Linearly Polarized Reverse Shock in GRB 190114C. Astrophys. J.; 2019; 878, L26. [DOI: https://dx.doi.org/10.3847/2041-8213/ab2247]

181. Selsing, J.; Fynbo, J.P.U.; Heintz, K.E.; Watson, D. GRB 190114C: NOT optical counterpart and redshift. GRB Coord. Netw.; 2019; 23695, 1.

182. Castro-Tirado, A.J.; Hu, Y.; Fernandez-Garcia, E.; Valeev, A.; Sokolov, V.; Guziy, S.; Oates, S.; Jeong, S.; Pandey, S.B.; Carrasco, I. et al. GRB 190114C: Refined redshift by the 10.4m GTC. GRB Coord. Netw.; 2019; 23708, 1.

183. Domínguez, A.; Primack, J.R.; Rosario, D.J.; Prada, F.; Gilmore, R.C.; Faber, S.M.; Koo, D.C.; Somerville, R.S.; Pérez-Torres, M.A.; Pérez-González, P. et al. Extragalactic background light inferred from AEGIS galaxy-SED-type fractions. Mon. Not. RAS; 2011; 410, pp. 2556-2578. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2010.17631.x]

184. Derishev, E.; Piran, T. GRB afterglow parameters in the era of TeV observations: The case of GRB 190114C. arXiv; 2021; arXiv: 2106.12035[DOI: https://dx.doi.org/10.3847/1538-4357/ac2dec]

185. Asano, K.; Murase, K.; Toma, K. Probing Particle Acceleration through Broadband Early Afterglow Emission of MAGIC Gamma-Ray Burst GRB 190114C. Astrophys. J.; 2020; 905, 105. [DOI: https://dx.doi.org/10.3847/1538-4357/abc82c]

186. Joshi, J.C.; Razzaque, S. Modelling synchrotron and synchrotron self-Compton emission of gamma-ray burst afterglows from radio to very-high energies. Mon. Not. RAS; 2021; 505, pp. 1718-1729. [DOI: https://dx.doi.org/10.1093/mnras/stab1329]

187. Dichiara, S.; Bernardini, M.G.; Burrows, D.N.; D’Avanzo, P.; Gronwall, C.; Gropp, J.D.; Kennea, J.A.; Klingler, N.J.; Krimm, H.A.; Kuin, N.P.M. et al. GRB 190829A: Swift detection of a burst consistent with a galaxy at z = 0.08. GRB Coord. Netw.; 2019; 25552, 1.

188. Fermi GBM Team. GRB 190829A: Fermi GBM Final Real-time Localization. GRB Coord. Netw.; 2019; 25551, 1.

189. Chand, V.; Banerjee, A.; Gupta, R.; Dimple,; Pal, P.S.; Joshi, J.C.; Zhang, B.B.; Basak, R.; Tam, P.H.T.; Sharma, V. et al. Peculiar Prompt Emission and Afterglow in the H.E.S.S.-detected GRB 190829A. Astrophys. J.; 2020; 898, 42. [DOI: https://dx.doi.org/10.3847/1538-4357/ab9606]

190. Piron, F.; Longo, F.; Axelsson, M.; Arimoto, M.; Racusin, J.L.; Bissaldi, E. Fermi-LAT Collaboration. GRB 190829A: Fermi-LAT Upper Limits. GRB Coord. Netw.; 2019; 25574, 1.

191. Osborne, J.P.; Page, K.L.; Capalbi, M.; Perri, M.; D’Elia, V.; Burrows, D.N.; Gropp, J.D.; Tohuvavohu, A.; Beardmore, A.P.; Dichiara, S. et al. GRB 190829A: Swift-XRT refined Analysis. GRB Coord. Netw.; 2019; 25568, 1.

192. Valeev, A.F.; Castro-Tirado, A.J.; Hu, Y.D.; Fernand ez-Garcia, E.; Sokolov, V.V.; Carrasco, I.; Castellon, A.; Garcia Alvarez, D.; Rivero, M. GRB 190829A: 10.4m GTC spectroscopy. GRB Coord. Netw.; 2019; 25565, 1.

193. Lipunov, V.; Balakin, F.; Gorbovskoy, E.; Kornilov, V.; Tyurina, N.; Balanutsa, P.; Kuznetsov, A.; Vladimirov, V.; Vlasenko, D.; Gorbunov, I. et al. GRB 190829A: MASTER confirmation of GROND SN. GRB Coord. Netw.; 2019; 25652, 1.

194. Laskar, T.; Bhandari, S.; Schroeder, G.; Fong, W.; Alexander, K.D.; Berger, E.; Chornock, R.; Coppejans, D.; Margutti, R.; Ayache, E. et al. GRB 190829A: ATCA cm-band detection. GRB Coord. Netw.; 2019; 25676, 1.

195. Monageng, I.; van der Horst, A.J.; Woudt, P.A.; Bottcher, M. GRB 190829A: MeerKAT radio observation. GRB Coord. Netw.; 2019; 25635, 1.

196. de Ugarte Postigo, A.; Bremer, M.; Kann, D.A.; Schulze, S.; Thoene, C.C.; Blazek, M.; Bensch, K.; Agui, J.F.; Izzo, L.; Perley, D.A. et al. GRB 190829A: NOEMA detection of the mm afterglow. GRB Coord. Netw.; 2019; 25589, 1.

197. Chandra, P. GRB 190829A: Detection of radio afterglow with the uGMRT. GRB Coord. Netw.; 2019; 25627, 1.

198. Sato, Y.; Obayashi, K.; Yamazaki, R.; Murase, K.; Ohira, Y. Off-axis jet scenario for early afterglow emission of low-luminosity gamma-ray burst GRB 190829A. Mon. Not. RAS; 2021; 504, pp. 5647-5655. [DOI: https://dx.doi.org/10.1093/mnras/stab1273]

199. Zhang, B.T.; Murase, K.; Veres, P.; Mészáros, P. External Inverse-Compton Emission from Low-luminosity Gamma-Ray Bursts: Application to GRB 190829A. Astrophys. J.; 2021; 920, 55. [DOI: https://dx.doi.org/10.3847/1538-4357/ac0cfc]

200. D’Elia, V.; Ambrosi, E.; Barthelmy, S.D.; D’Ai, A.; Gropp, J.D.; Klingler, N.J.; Lien, A.Y.; Palmer, D.M.; Sbarufatti, B.; Siegel, M.H. et al. GRB 201015A: Swift detection of a burst. GRB Coord. Netw.; 2020; 28632, 1.

201. Fletcher, C.; Veres, P. Fermi-GBM Team. Fermi GBM Sub-Threshold Detection of GRB 201015A. GRB Coord. Netw.; 2020; 28663, 1.

202. Minaev, P.; Pozanenko, A. GRB 201015A: Classification as long GRB. GRB Coord. Netw.; 2020; 28668, 1.

203. Markwardt, C.B.; Barthelmy, S.D.; Cummings, J.R.; D’Elia, V.; Krimm, H.A.; Laha, S.; Lien, A.Y.; Palmer, D.M.; Sakamoto, T.; Stamatikos, M. et al. GRB 201015A: Swift-BAT refined analysis (a soft short pulse with a tail emission). GRB Coord. Netw.; 2020; 28658, 1.

204. D’Ai, A.; Gropp, J.D.; Kennea, J.A.; Tohuvavohu, A.; Page, K.L.; Beardmore, A.P.; Evans, P.A.; Melandri, A.; D’Elia, V. Swift-XRT Team. GRB 201015A: Swift-XRT refined Analysis. GRB Coord. Netw.; 2020; 28660, 1.

205. Gompertz, B.; Levan, A.; Tanvir, N.; Fruchter, A.; Cucchiara, A.; Greiner, J.; Hjorth, J.; Kangas, T.; Lamb, G.; Lyman, J. et al. GRB 201015A: Late X-ray Detections with Chandra. GRB Coord. Netw.; 2020; 28822, 1.

206. D’Elia, V. Swift Team. GRB 201015A: Swift-XRT late time observations. GRB Coord. Netw.; 2020; 28857, 1.

207. Lipunov, V.; Gorbovskoy, E.; Kornilov, V.; Tyurina, N.; Balanutsa, P.; Kuznetsov, A.; Balakin, F.; Vladimirov, V.; Vlasenko, D.; Gorbunov, I. et al. Swift GRB 201015A: Global MASTER-Net observations report. GRB Coord. Netw.; 2020; 28634, 1.

208. Grossan, B.; Maksut, Z.; Kim, A.; Krugov, M.; Smoot, G.F.; Linder, E. GRB 201015A NUTTelA-TAO/BSTI Early Measurements (Preliminary). GRB Coord. Netw.; 2020; 28674, 1.

209. Fong, W.; Schroeder, G.; Rastinejad, J.; Hajela, A. GRB 201015A: 6 GHz VLA radio afterglow candidate detection. GRB Coord. Netw.; 2020; 28688, 1.

210. Giarratana, S.; Giroletti, M.; Marcote, B.; Ghirlanda, G.; Ribó, M.; Paredes, J.M. GRB 201015A radio afterglow with e-MERLIN. GRB Coord. Netw.; 2020; 28939, 1.

211. Marcote, B.; Ribó, M.; Paredes, J.M.; Giroletti, M.; Giarratana, S.; Ghirlanda, G. GRB 201015A: 5-GHz radio afterglow detection with the EVN. GRB Coord. Netw.; 2020; 29028, 1.

212. Pozanenko, A.; Belkin, S.; Volnova, A.; Moskvitin, A.; Burhonov, O.; Kim, V.; Krugov, M.; Rumyantsev, V.; Klunko, E.; Inasaridze, R.Y. et al. GRB 201015A: Optical observations and supernova identification. GRB Coord. Netw.; 2020; 29033, 1.

213. Rossi, A.; Benetti, S.; Palazzi, E.; D’Avanzo, P.; D’Elia, V.; De Pasquale, M. CIBO Collaboration. GRB 201015A: Evidence of supernova in LBT spectra. GRB Coord. Netw.; 2021; 29306, 1.

214. de Ugarte Postigo, A.; Kann, D.A.; Blazek, M.; Agui Fernandez, J.F.; Thoene, C.; Gomez Velarde, G. GRB 201015A: Redshift from GTC/OSIRIS. GRB Coord. Netw.; 2020; 28649, 1.

215. Izzo, L.; Malesani, D.B.; Zhu, Z.P.; Xu, D.; de Ugarte Postigo, A.; Pursimo, T. GRB 201015A: Redshift confirmation. GRB Coord. Netw.; 2020; 28661, 1.

216. Blanch, O.; Gaug, M.; Noda, K.; Berti, A.; Moretti, E.; Miceli, D.; Gliwny, P.; Ubach, S.; Schleicher, B.; Cerruti, M. et al. MAGIC observations of GRB 201015A: Hint of very high energy gamma-ray signal. GRB Coord. Netw.; 2020; 28659, 1.

217. Suda, Y.; Artero, M.; Asano, K.; Berti, A.; Nava, L.; Noda, K.; Terauchi, K. MAGIC Collaboration. Observation of a relatively low luminosity long duration GRB 201015A by the MAGIC telescopes. Proceedings of the 37th International Cosmic Ray Conference—PoS (ICRC2021); Berlin, Germany, 12 July 2021; [DOI: https://dx.doi.org/10.22323/1.395.0797]

218. Beardmore, A.P.; Gropp, J.D.; Kennea, J.A.; Klingler, N.J.; Laha, S.; Lien, A.Y.; Moss, M.J.; Page, K.L.; Palmer, D.M.; Tohuvavohu, A. et al. GRB 201216C: Swift detection of a burst. GRB Coord. Netw.; 2020; 29061, 1.

219. Malacaria, C.; Veres, P.; Meegan, C.; Bissaldi, E. Fermi GBM Team. GRB 201216C: Fermi GBM detection. GRB Coord. Netw.; 2020; 29073, 1.

220. Ukwatta, T.N.; Barthelmy, S.D.; Beardmore, A.P.; Krimm, H.A.; Laha, S.; Lien, A.Y.; Markwardt, C.B.; Palmer, D.M.; Sakamoto, T.; Stamatikos, M. GRB 201216C: Swift-BAT refined analysis. GRB Coord. Netw.; 2020; 29080, 1.

221. Bissaldi, E.; Omodei, N.; Kocevski, D.; Axelsson, M.; Longo, F.; Moretti, E. Fermi-LAT Collaboration. GRB 201216C: Fermi-LAT upper limit. GRB Coord. Netw.; 2020; 29076, 1.

222. Evans, P.A.; Osborne, J.P.; Bernardini, M.G.; Ambrosi, E.; Capalbi, M.; Tohuvavohu, A.; Sbarufatti, B.; Burrows, D.N.; Beardmore, A.P. Swift-XRT Team. GRB 201216C: Swift-XRT refined Analysis. GRB Coord. Netw.; 2021; 29278, 1.

223. Izzo, L.; Malesani, D.B.; Kann, D.A. GRB 201216C: VLT afterglow candidate discovery. GRB Coord. Netw.; 2020; 29066, 1.

224. Jelinek, M.; Strobl, J.; Karpov, S.; Masek, M.; Janecek, P.; Jurysek, J.; Ebr, J.; Cunniffe, R.; Travnicek, P.; Prouza, M. GRB 201216C: FRAM-ORM afterglow confirmation. GRB Coord. Netw.; 2020; 29070, 1.

225. Shrestha, M.; Melandri, A.; Smith, R.; Steele, I.A.; Kobayashi, S.; Mundell, C.; Gomboc, A.; Guidorzi, C. GRB 201216C: Liverpool Telescope First Hour Observations. GRB Coord. Netw.; 2020; 29085, 1.

226. Ayala, H. GRB201216C: No significant detection in HAWC. GRB Coord. Netw.; 2020; 29086, 1.

227. Vielfaure, J.B.; Izzo, L.; Xu, D.; Vergani, S.D.; Malesani, D.B.; de Ugarte Postigo, A.; D’Elia, V.; Fynbo, J.P.U.; Kann, D.A.; Levan, A.J. et al. GRB 201216C: VLT X-shooter spectroscopy and potential high redshift of a VHE-emitting GRB. GRB Coord. Netw.; 2020; 29077, 1.

228. Blanch, O.; Longo, F.; Berti, A.; Fukami, S.; Suda, Y.; Loporchio, S.; Micanovic, S.; Green, J.G.; Pinter, V.; Takahashi, M. et al. GRB 201216C: MAGIC detection in very high energy gamma rays. GRB Coord. Netw.; 2020; 29075, 1.

229. Fukami, S.; Berti, A.; Loporchio, S.; Suda, Y.; Nava, L.; Noda, K.; Bosnjak, Z.; Asano, K.; Longo, F. MAGIC Collaboration. Very-high-energy gamma-ray emission from GRB 201216C detected by MAGIC. Proceedings of the 37th International Cosmic Ray Conference—PoS (ICRC2021); Berlin, Germany, 12 July 2021; Volume 395, 788. [DOI: https://dx.doi.org/10.22323/1.395.0788]

230. Nava, L.; Salvaterra, R.; Ghirlanda, G.; Ghisellini, G.; Campana, S.; Covino, S.; Cusumano, G.; D’Avanzo, P.; D’Elia, V.; Fugazza, D. et al. A complete sample of bright Swift long gamma-ray bursts: Testing the spectral-energy correlations. Mon. Not. RAS; 2012; 421, pp. 1256-1264. [DOI: https://dx.doi.org/10.1111/j.1365-2966.2011.20394.x]

231. D’Avanzo, P.; Salvaterra, R.; Bernardini, M.G.; Nava, L.; Campana, S.; Covino, S.; D’Elia, V.; Ghirlanda, G.; Ghisellini, G.; Melandri, A. et al. A complete sample of bright Swift short gamma-ray bursts. Mon. Not. RAS; 2014; 442, pp. 2342-2356. [DOI: https://dx.doi.org/10.1093/mnras/stu994]

232. Wei, J.; Cordier, B.; Antier, S.; Antilogus, P.; Atteia, J.L.; Bajat, A.; Basa, S.; Beckmann, V.; Bernardini, M.G.; Boissier, S. et al. The Deep and Transient Universe in the SVOM Era: New Challenges and Opportunities—Scientific prospects of the SVOM mission. arXiv; 2016; arXiv: 1610.06892

233. Lien, A.; Sakamoto, T.; Barthelmy, S.D.; Baumgartner, W.H.; Cannizzo, J.K.; Chen, K.; Collins, N.R.; Cummings, J.R.; Gehrels, N.; Krimm, H.A. et al. The Third Swift Burst Alert Telescope Gamma-Ray Burst Catalog. Astrophys. J.; 2016; 829, 7. [DOI: https://dx.doi.org/10.3847/0004-637X/829/1/7]

234. Donini, A.; Gasparetto, T.; Bregeon, J.; Pierro, F.D.; Longo, F.; Maier, G.; Moralejo, A.; Vuillaume, T. The Cherenkov Telescope Array Performance in Divergent Mode. Proceedings of the 36th International Cosmic Ray Conference (ICRC2019); Madison, WI, USA, 24 July–1 August 2019; Volume 36, 664.

235. Inoue, S.; Granot, J.; O’Brien, P.T.; Asano, K.; Bouvier, A.; Carosi, A.; Connaughton, V.; Garczarczyk, M.; Gilmore, R.; Hinton, J. et al. Gamma-ray burst science in the era of the Cherenkov Telescope Array. Astropart. Phys.; 2013; 43, pp. 252-275. [DOI: https://dx.doi.org/10.1016/j.astropartphys.2013.01.004]

236. Grazia Bernardini, M.; Bissaldi, E.; Bosnjak, Z.; Carosi, A.; D’Avanzo, P.; Di Girolamo, T.; Inoue, S.; Gasparetto, T.; Ghirlanda, G.; Longo, F. et al. POSyTIVE—A GRB population study for the Cherenkov Telescope Array (ICRC-2019). arXiv; 2019; arXiv: 1908.01544

237. Stamerra, A.; Saturni, F.G.; Green, J.G.; Nava, L.; Lucarelli, F.; Antonelli, L.A. TeV Transients with the ASTRI Mini-Array: A case study with GRB 190114C. Proceedings of the 37th International Cosmic Ray Conference—PoS (ICRC2021); Berlin, Germany, 12 July 2021; Volume 395, 890. [DOI: https://dx.doi.org/10.22323/1.395.0890]