Gamma-ray bursts afterglow physics and the VHE domain

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 produced also 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, showing how future observations, especially {by} the CTA Observatory, are expected to give a key contribution in improving our comprehension of such elusive sources.


Introduction
Gamma-Ray Bursts (GRBs) are observed as transient sources of radiation displaying a distinctive pattern consisting 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 seconds [1,2]. The second emission phase, called afterglow, follows the prompt with a delay of tens of seconds, and is detected on a very wide range of frequencies, from γ-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 (Γ ∼ 100 − 1000) jet from a newly born compact object. The ejecta undergoes first 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 arXiv:2205.12146v1 [astro-ph.HE] 24 May 2022 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 (θ jet ∼ 5 − 10 • ) outflows. A similar scenario applies also 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 explaining more naturally 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 into 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 on 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 can be still consistent with internal shocks. The nature and efficiency of the dissipation mechanism in the prompt phase are still matter of intense debate. In any case, the radiation is expected to be produced by the accelerated electrons, which efficiently loose energy via synchrotron cooling [6,25]. Inconsistencies between the expected synchrotron spectrum and the observed spectral shape of the prompt emission [25,26] have called into question also 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 done in the last decade from 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 γ-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.

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 characterises 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 though to proceed via diffusive shock acceleration, but for which the details of the underlying physics remain still poorly constrained. Moreover, the overall luminosity of the afterglow radiation depends on the energy content of the blast-wave. Such 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).

Jet dynamics
After the reverse shock has crossed the ejecta, the dynamics of the blast-wave enters a selfsimilar 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(R dec )c 2 Γ 2 0 ) is comparable to the initial energy (E 0 = M 0 Γ 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 and on the ejecta mass M 0 and initial bulk Lorentz factor Γ 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 interested also 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 to describe the hydrodynamics of a relativistic blastwave expanding into a medium with an arbitrary density profile ρ(R) 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 Γ considered is the one of the fluid 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 to introduce time-varying radiative efficiency, either resulting from a change with time of 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].

Equation describing the evolution of the bulk Lorentz factor
The aim is to derive an equation describing the change dΓ of the bulk Lorentz factor of the fluid just behind the shock in response to the collision with a mass dm(R) = 4πR 2 ρ(R)dR encountered when the shock front moves from a distance R to R + dR and with ρ being the mass density. The change in Γ is determined by dissipation of the bulk kinetic energy, conversion of internal energy back into bulk motion, and injection of energy into the blastwave. 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Γ/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 = u int + ρ c 2 , where u int is the comoving internal energy and ρ is the comoving mass density. Applying Lorentz transformations, u = (u + p )Γ 2 − p , where p is the pressure and is related to the internal energy density by the equation of state p = (γ − 1) u int = (γ − 1) (u − ρ c 2 ), whereγ is the adiabatic index of the shocked plasma. The energy density is then given by: u = u int (γΓ 2 −γ + 1) + ρ c 2 Γ 2 , which shows how the internal energy and rest mass density transform. The total energy in the progenitor frame will be E = uV = uV /Γ, where V is the shell volume in the progenitor frame, and can be expressed as: where: which properly describes the Lorentz transformation of the internal energy. Here, M = M 0 + m = ρ V is the sum of the ejecta mass M 0 = E 0 /Γ 0 c 2 and of the swept-up mass m(R), and E int = (u − ρ c 2 )V is the comoving internal energy. The adiabatic index can be parameterized asγ = (4 + Γ −1 )/3 to obtain the expected limitsγ 4/3 for Γ 1 and γ 5/3 for Γ → 1. The majority of analytical treatments use Γ instead of Γ eff , which implies an error up to a factor of 4/3 in the ultra-relativistic limit [55].
The blast-wave energy E in Eq. 1 can change due to (i) the rest mass energy dm c 2 collected from the medium, (ii) radiative losses dE rad = Γ eff dE rad , 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: The overall change in the comoving internal energy dE int results from the sum of three contributions: The first contribution, dE sh = (Γ − 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 dE sh /dm is constant across the shock, and equal to (Γ − 1)c 2 . The second term in eq. 4, dE ad , 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, dE rad , accounts for radiative losses. From Eq. 3, it follows that the variation of the Lorentz factor is: from which the evolution Γ(R) 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 Γ eff dE ad /dR, accounting for adiabatic losses, allows to describe the re-acceleration of the fireball: this contribution, usually neglected, becomes important only when the density decreases faster than ρ ∝ R −3 . To evaluate Eq. 5 it is necessary to first specify dE ad and E int .

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 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 fast cooling regime or quasi-adiabatic regime, in which case the radiative losses do not affect the shell dynamics.
Defining e as the fraction of energy dE sh 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): γ acc,e − 1 = (Γ − 1) e /µ e .
Here, µ e = ρ e /ρ is the ratio between the mass density ρ e of shocked electrons and positrons (simply "electrons" from now on) and the total mass density of the shocked matter ρ. In the absence of electron-positron pairs µ e = m e /(m e + m p ) m e /m p .
Leptons then radiate a fraction rad of their internal energy, i.e., the energy lost to radiation is dE rad = − rad e dE sh = − dE sh , with ≡ rad e being the overall fraction of the shockdissipated energy that goes into radiation. After radiating a fraction rad of their internal energy, the mean random Lorentz factor of the freshly shocked electrons decreases down to: The assumption of instantaneous radiative losses is verified in the fast cooling regime ( rad ∼ 1), which is required (but not sufficient) to have ∼ 1 (i.e., a fully radiative blast-wave). In the opposite case rad 1, the evolution is nearly adiabatic ( 1), regardless of the value of e , and the details of the radiative cooling processes are likely to be unimportant for the shell dynamics. The case with intermediate values of rad and is harder to treat analytically, since the electrons shocked at radius R may continue to emit copiously also at larger distances, affecting the blast-wave dynamics.
A similar treatment can be adopted for protons: if protons gain a fraction p of the energy dissipated by the shock (with p = 1 − e − B ), their mean post-shock Lorentz factor will be: where µ p = ρ p /ρ is the ratio between the mass density of shocked protons ρ p and the total shocked mass density ρ. In the standard case, when pairs are absent, µ p 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 dE int = −p dV , where p is the pressure in comoving frame. For N particles with Lorentz factor γ the internal energy density is: The radial change of the Lorentz factor, as a result of expansion losses, is: To estimate the adiabatic losses, let us assume that the shell comoving volume scales as V ∝ R 3 /Γ, corresponding to a shell thickness in the progenitor frame ∼ R/Γ 2 . This scaling is correct for both relativistic and non-relativistic shocks, in the decelerating phase (BM76). For re-accelerating relativistic shocks, Shapiro [59] showed that the thickness of the region containing most of the blast-wave energy is still ∼ R/Γ 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 ∝ R 3 /Γ is assumed, the equation can be further developed analytically, and reads: 8 of 76 The comoving Lorentz factor at radius R, for a particle injected with γ(r) when the shock radius was r, will be where γ(r) is given by γ rad,e (r) (Eq. 7) for leptons, and by γ acc,p (r) (Eq. 8) for protons. Considering the proton and lepton energy densities separately, the comoving internal energy at radius R will be: With the help of Eq. 12, one can explicitly find E int (R) and insert it in Eq. 5.
The other term needed in Eq. 5 is dE ad /dR. First, we have derived (dγ/dR) ad for a single particle. Now integrating over the total number of particles, again considering protons and leptons separately, one obtains: In Eqs. 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: In the absence of significant magnetic field amplification, p + e 1 so that p + e − 1 − , and the radiative processes of the blast-wave are entirely captured by the single efficiency parameter . In the fast cooling regime rad ∼ 1 and e . In this case the term p + e − reduces to 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 ρ ∝ R −s , one obtains: where Γ ∝ R −(3−s)/2 as in the adiabatic BM76 solution has been used.
In the fully radiative regime = 1, which implies E int = 0 and dE ad = 0, Eq. 5 reduces to: which describes the evolution of a momentum-conserving (rather than pressure-driven) snowplow. Replacing Γ eff → Γ, 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 so that the BM76 normalization can be recovered by multiplying the density of external matter in Eqs. 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: No analytic model properly captures the transition between an adiabatic relativistic blastwave and the momentum-conserving snowplow, as increases from zero to unity. The simple interpolation in Eq. 19 joins the fully adiabatic BM76 solution with the fully radiative momentum-conserving snowplow.
In summary, Eqs. 5, 13 and 14, complemented with the correction in Eq. 19 (which should by applied to every occurrence of external density and external matter) completely determine the evolution of the shell Lorentz factor Γ as a function of the shock radius R.

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 (Γ GRB >> 1, decreasing to non-relativistic velocity only after several weeks or months), and the magnetization σ (i.e., the ratio between Poynting and kinetic flux in the pre-shocked fluid, σ = B 2 /(4π 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 summarise 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 γ min/max and the spectral index p), ii) the acceleration efficiency (i.e., the fraction of electrons ξ 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.

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 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 (γβ 1), the estimates of the spectral slopes converge to a universal value p = s − 2 ∼ 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/ω pe of the upcoming electrons (where ω pe = 4π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 ∼ 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 shown that as long as the fluid is ultra-relativistic (Γ > 5), the main parameter governing the acceleration is the magnetization σ, i.e. the efficiency of the process is insensitive to Γ, 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 Γ = 15, Figure 1. Spectral index (s p = p + 2) 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 s p ∼ 4.2 (corresponding to p ∼ 2.2) in the ultra-relativistic limit.
ion-to-electron mass ratio m i /m e = 25, and σ = 10 −5 . The temporal evolution is followed up to t = 2500 ω −1 pi . The formation of a non-thermal tail is clearly visible. The downstream non-thermal population is found to include around ξ 3% of the electrons, carrying e 10% of the energy. The spectral index is around p ∼ 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 γ max increases proportionally to t 1/2 (see inset in Figure 2), slower than the commonly adopted Bhom rate [64], in which case γ max ∝ t. Extrapolating the γ max behaviour 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 γ max ∼ 10 7 in the early phase of GRB afterglows, corresponding to synchrotron photon energies around 1 GeV, 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 computationally demanding) is how the magnetization evolves downstream. PIC simulations have found values of B ∼ 0.1 − 0.01 in the vicinity of the shock front. How this turbulence evolves on longer time-scales is still 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 The evolution is followed until t = 2500ω −1 pi . From [42]. ©AAS. Reproduced with permission.
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 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.

Description of shocks in GRB afterglow modeling
The theory of relativistic shock acceleration is applied to 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 ξ 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 dE sh = (Γ − 1)dmc 2 (see section 2.1), the corresponding energy density is u sh = (Γ − 1)ρ c 2 . From shock jump conditions, the density in the comoving frame ρ is related to the density of the unshocked medium (measured in the rest frame) by the equation: 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: ρ = nm 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: A fraction B of this energy will be conveyed to the magnetic field: from which it follows that the magnetic field strength B is: Similarly, for the accelerated electrons: u e = e u sh = e 4Γ(Γ − 1)nm p c 2 = γ m e c 2 4Γξ e n (24) where γ is the average random Lorentz of the accelerated electrons: 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(γ)dγ ∝ γ −p dγ for γ min ≤ γ ≤ γ max where γ min is the minimum Lorentz factor of the injected electrons and γ max is the maximum Lorentz factor at which electrons can be accelerated. To derive the relation between γ min , γ max and the model parameters, we consider the definition of the average Lorentz factor γ : and solve the integrals. Equating 25 and 26 leads to (for p = 1): A simplified equation for γ min can be obtained assuming that γ −p+2 max γ −p+2 min : Since p is expected to be 2 < p < 3, this condition is verified for γ max γ min . The minimum Lorentz factor is then not treated as a free parameter of the model, as it is calculated from eq. 28 as a function of the free parameters e , ξ e and p. Concerning the prescription for the value of γ max (for details see Section 3.5) usually it relies on the condition that radiative losses between acceleration episodes are equal to the energy gains, where energy gains proceeds 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: Solving the equations assuming that γ max,p γ min,p leads to: The equations for γ min and Eq. 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 Γ 0 expanding in a medium with density n(R). 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.

Constraints to the acceleration mechanism provided by observations
Assuming that accelerated particles have a power-law spectrum (dN acc /dγ ∝ γ −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 ν m of photons emitted by electrons with Lorentz factor γ min (this is correct both in case of fast and slow cooling regime). The estimated value of p from afterglow modeling are spread on a wide range, from p ∼ 2 to p ∼ 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 density profile of the external medium, which is often unconstrained (see section 3.2).
The typical value of e inferred from 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 values predicted by numerical investigations [42]. For the fraction 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 shown 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 0.1) in the vicinity of the shock, where particle are 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 ξ 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, ξ 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 γ 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 difficult to assess from spectral analysis whether they are consistent with the power-law extrapolation of the synchrotron spectrum or 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 r L /c) with r L = E/eB being the Larmor radius, and is limited by synchrotron cooling (t syn 6 π m e c/σ T B 2 γ) then γ max ∼ 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 γ max inferred from the detection of GeV photons is quite large and barely consistent with what found by PIC simulations. Whether or not the observations are in tension with the present derivation of γ 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.

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 ν sa , the characteristic frequency ν m emitted by electrons with Lorentz factor γ min , the cooling frequency ν c emitted by electrons with Lorentz factor γ 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 Thomson regime [48] and corrections to be applied to the synchrotron and SSC spectra to account for the effects of Klein-Nishina [75] cross section (see Section 2.3.1) are also available in 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(R) = 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 Γ 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 solving numerically the differential equation describing the evolution of the particle spectra and estimate the associated emission [76][77][78]. In this section we describe a radiative code that solves simultaneously 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: whereγ = ∂γ ∂t is the rate of change of the Lorentz factor γ of an electron caused by adiabatic, synchrotron and SSC losses and to 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 scenario will generate a very high energy spectral component, which is the target of searches by IACTs such as MAGIC 1 and H.E.S.S. 2 . In principle, also up-scattering of an external population of seed photons can 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(γ, t ) = Q acc (γ, t ) + Q pp (γ, t ) describes the injection of freshly accelerated particles (Q acc (γ, t ) = dN acc /dγ dt ) and the injection of pairs Q pp (γ, t ) produced by photon-photon annihilation.
In the next sections we explicit each one of the terms included in eq. 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.

Synchrotron and SSC cooling
The synchrotron power emitted by an electron with Lorentz factor γ 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 γ is given by: The cross section for the inverse Compton mechanism is constant and equal to the Thomson cross section (σ 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: where:ε . (34) ε and ε are the energies of the photons (normalized to the rest mass electron energy) before and after the scattering process, respectively. The two terms of Eq. 33 account respectively for the down-scattering (i.e. ε <ε) and the up-scattering (i.e. ε >ε) process. The energy loss term for the SSC can now be calculated with the equation:

Adiabatic cooling
As discussed in section 2.1, particles loose 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γ in a comoving time dt : with β being the random velocity of particles in unit of c. The comoving volume V 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 = dR/Γ(R) c the comoving volume is: and:γ 18 of 76

Synchrotron self-absorption (SSA)
Electrons can re-absorb low energy photons before they escape from the source region. The absorption coefficient α ν can be expressed as [80]: valid for any radiation mechanism at the emission frequency ν , with P (γ, ν ) being the specific power of electrons with Lorentz factor γ at frequency ν and assuming hν γm e c 2 . Thus, the SSA mechanism will affect mostly the low frequency range. This results in a modification of the lower frequency tail of the synchrotron spectrum as: assuming a power-law distribution of electrons, with ν i = min(ν m , ν cool ) and ν SSA the frequency below which the synchrotron flux is self-absorbed and the source becomes optically thick.

Synchrotron and Inverse Compton emission
Following [82], the synchrotron spectrum emitted by an electron with Lorentz factor γ, averaged over an isotropic pitch angle distribution is: where x ≡ ν 4π m e c/(6 e B γ 2 ), and K n are the modified Bessel functions of order n. The total power emitted at the frequency ν is obtained integrating over the electron distribution: The SSC radiation emitted by an electron with Lorentz factor γ can be calculated as: where nν 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 ν .

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 σ γγ as a function of β , the centre-of-mass speed of the electron and positron is given by: where: (45) and ω t = hν t /m e c 2 with ν t being the target photon frequency, ω s = hν /m e c 2 with ν being the source photon frequency and µ = cos φ, where φ is the scattering angle. Then, it is possible to derive the annihilation rate of photons into electron-positron pairs as: where x. An accurate and simple approximation which takes into account both regimes is given by: where H(x − 1) is the Heaviside function [83]. The approximation reproduces accurately the behaviour near the peak at x peak ∼ 3.7 and over the range 1.3 < x < 10 4 which usually dominates during the calculations. A comparison between Eq. 46 and Eq. 47 is given in Figure  3 where the goodness of the approximation adopted in the mentioned x range can be seen. The impact of the flux attenuation due to pair production mechanism on the GRB spectra is Figure 3. Comparison between the exact annihilation rate (eq. 46) and the approximated formula (eq. 47).
The ratio between the two curves in the range (1.3 < x < 10 4 ) is shown in the bottom panel. In this range the ratio is always 7%. estimated in terms of the optical depth value τ γγ . From its definition: where n (ν t ) is the number density of the target photons per unit of volume, σ ν ν t is the cross section and ∆R is the width of the emission region. Introducing the cross section in terms of the annihilation rate R(x) in its approximated formula and integrating over all the possible target photon frequencies: where ν and ν t 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 −τ γγ )/τ γγ . This attenuation factor will modify the GRB spectrum giving a non-negligible contribution especially in the VHE domain. An example of the modification of a GRB spectrum due to pair production can be seen 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 ∆R = R/Γ, 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 be done also for the electron/positron production. Assuming that the electron and positron arises with equal Lorentz factor γ and that x peak ∼ 3.7, a photon with energy ω s 1 will mostly interact with a target photon of energy ω t ≈ 3/ω s . Then, from the energy conservation condition: The e ± production can be seen as an additional source term for the distribution of accelerated particles. As a result, an additional injection term Q pp e to be inserted in the kinetic equation (equation 31) is calculated as:

Comparison with analytical approximations
In order to compare results from the numerical method described in previous section and analytical prescriptions available in literature, we give an example in Figure 5. The analytical prescriptions for the synchrotron and the SSC component are calculated following [38] and [48]. In [38] the synchrotron spectra and light-curves are derived assuming a powerlaw 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 seen that analytical and numerical results are in general in good agreement. Both curves follow the same behaviour 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 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 t = 10 4 s for s = 0, p = 2.3, e = 0.05, 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.

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 behaviour 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 radioto-GeV radiation is indeed produced in interactions between the ejecta and the external medium. Also, 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 which 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.

X-ray flares
Observations of the afterglow emission in X-ray and optical often display behaviours 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 X-ray and optical band are expected to decay following a power-law or a broken power-law behaviour, 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 characterised by a very fast rise of the flux, followed by an exponential decay profile. Comprehensive studies of X-ray afterglows show 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 seconds, even though there are many cases of flares occurring several hours after the burst. The width of the flare ω is found to evolve linearly with time to larger values following the trend ω ∼ 0.2T pk [84]. The average and peak luminosities L of the flare also display a dependence from T pk , with L ∝ T −2.7 pk at least for early time (T pk < 10 3 s) flares [84,85]. When including also late time flares [86,87] a shallower index is obtain, 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 been detected also 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 shows 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 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 ∼ 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 γ min ∼ 60. Such value is estimated from theoretical considerations where γ min 60 e,−1 (Γ IS − 1) for p = 2.5, e = 0.1 and a relative shock Lorentz factor Γ 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, 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 are the chances of current and future VHE ground based instruments to contribute 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]. 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 × 10 −9 erg cm −2 s −1 and F = 1 × 10 −10 erg cm −2 s −1 , respectively. Assuming that a similar amount of flux is emitter 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 5 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 ω ∼ 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 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 × 10 −9 erg cm −2 s −1 and at t obs = 200 s is ∼ 5.0 × 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, that 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 it can be seen from 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 assumed here).
Concerning future instruments, the Cherenkov Telescope Array (CTA 3 will have a sensitivity which is almost one order of magnitude lower than the MAGIC one and similar slewing capabilities. The same estimates done 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 counterpart at E ∼ 30 GeV of X-ray flares, provided that this counterpart has a flux is which no less than ten times smaller than what 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.

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 radial density profile n ∝ R −2 . However, several investigations have shown that about half of the long GRBs are better explained if the blastwave 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 are still uncertain. The progenitor is usually identified with Wolf-Rayet stars, massive stars (M > 20 M ) in the final stages of their evolution, characterised 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 leave imprints on its afterglow emission.
More in detail, from the interaction between the stellar wind and the ISM four concentric regions with different properties are 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 ∝ R −2 . The density is related to the mass loss rateṀ and to the velocity v w of the free-streaming stellar wind by: A termination shock separates the unshocked from the shocked wind: the latter forms a hot bubble of thermalised 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 blastwave has decelerated to non-relativistic velocity) and compare them to the termination shock radius. For typical parameters (Ṁ = 10 −5 M 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 = 10 n −1/2 ISM pc , where n ISM is the density of the unshocked ISM. From the blast-wave dynamics, the deceleration and the non-relativistic 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 eq. 52, whereṀ and v w are treated as unknown parameters (normalised to the typical values of a Wolf-Rayet star) combined in one single free model parameter A : n(R) = 3 × 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 cases 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 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Ṁ and v w is shown in the bottom panel of Figure 8. A very low mass loss rateṀ = 10 −7 M yr −1 (which may find a justification in case of 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 ∼ 6 × 10 −3 pc and R NT ∼ 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 γ 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 diversity of E k ,Ṁ, 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 Klein-Nishina regime (see also [108]). Moreover these are based on simple approximation of the synchrotron spectrum into power-law segments, while the wide curvature of the real synchrotron spectra might lead to an 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 behaviour does not depends on p and is quite different in case of constant or wind-like density profile. Similarly, the detection of SSC radiation can help solving this ambiguity.

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 magnetised shocks invoke the existence of self-generated micro-turbulence both behind and in front of the shock, at a level corresponding to B ∼ 0.01 − 0.1. This layer of intense microturbulence is expected on theoretical grounds and recently corroborated by numerical PIC simulations. Inferences of the value of B from early modeling of afterglow radiation were broadly consistent with these numbers, confirming the presence of a 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 requires 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 Xray 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 a small values of B may appear to be problematic [114], because strongly selfgenerated 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 timescale of the system. Dedicated simulations the magnetic field does decay behind the shock, on a time-scale much longer than cω pi . Immediately behind the shock, the magnetic field carries a magnetization B ∼ 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Γ 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 ∝ t −α t with α t ∼ 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.

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 microphysical parameters (such as e , e , p and ξ e ) are constant over the whole afterglow evolution can be testes. 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 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 ν cool evolves as ν cool ∝ t −1.2 which is in contrast with synchrotron predictions for which a less steeper decay ν cool ∝ 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 behaviour without affecting the general interpretation of the data is to let the B parameter variate with time. Assuming that B ∝ t 0.49 the time evolution of ν cool can be explained successfully.
In [116] for GRB 130427A modeling, in order to explain the observed fast evolution of the injection frequency ν m ∝ t −1.9 a temporal evolution of e is claimed. Considering that ν m ∝ 2 e , a modest evolution of e following the trend e ∝ 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 including also TeV emission will give new insights useful to investigate 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 evolution of the micro-physical parameters with time [121], as will be discussed in the next section.

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: where r L is the Larmor radius. For each crossing the electrons gain energy by a factor ∼ 2. On the other hand, the energy losses for synchrotron radiation on this timescale are: The particle stops to gain energy when: Therefore, the maximum Lorentz factor for electrons can be derived: The corresponding maximum synchrotron photon energy is: which for electrons is ∼ 50 MeV in their rest frame. Similar considerations can be done also for protons. Following the same arguments presented below one obtains: for the cooling Lorentz factor and: 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 modelled 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 seen. Such 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.
Calculation performed so far assumes the presence of a uniform magnetic field B throughout 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 can be expressed in terms of the distance from the shock front x as: where B s and B w are respectively the strongest and the weakest magnetic field strengths, η is the power-law decaying index, and L p is the field decay length scale, which is estimated as [126]: where Γ 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 (x) becomes weaker and an electron travelling downstream will be likely sent back upstream when r L ≤ x. When considering the case B s B w and x 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: valid for η > 1/2 and x 0 /L p 1 where x 0 is the width of the high magnetic field region. Considering that x 0 /L p ≡ (B s /B w ) 1/η , eq. 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 loose most of their energy in the weak magnetic field region following the same conditions presented for the uniform magnetic field case: As a result, the maximum synchrotron photon energy is given by: which is greater than the one calculated in eq. 57 by a factor B s /B w . 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 loosing 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.

Prompt emission efficiency
The overall efficiency η γ of the prompt emission mechanism is the result of three processes: the efficiency η diss of the (still unidentified) mechanism responsible for 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: η γ = η 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 provides then a fundamental clue to place 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 η γ = E γ /E 0 (where E 0 = E γ + E k is the initial explosion energy), we can write the relation E k = (1 − η γ )/η γ E γ . The parameter η γ can then be estimated from the comparison of the energy E γ 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 one 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(ν m , nu c ) and is not affected by inverse Compton cooling. If these two conditions are verified, then the electrons emitting X-ray photons are in 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 η γ , between 0.5 and 0.9 [49,129,132,133].
The very same approach can be applied also 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 γ,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 use, 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 ν c > ν X (where ν 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 ν 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, affecting also the estimates of η γ , which are around 5 − 10% [24].
A correct estimate of η γ is extremely important, since its value is related to the mechanism dissipating energy in the jet. Since internal shocks can hardly reach values of η γ 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 η γ lead to different estimates on the total initial jet energy E 0 = E γ,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 difficult to attain also 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 more precisely the energy budget of the jet during the afterglow emission and add important information to constrain the efficiency of the prompt emission and favour or exclude some dissipation scenarios.

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 ξ e is responsible to account for the fraction of particles (here we consider electrons but same considerations are valid also for protons) accelerated into a non-thermal distribution. This means that from relativistic shock theory it is expected that a fraction 1 − ξ e of electrons is heated into a thermal distribution rather than a non-thermal one (see Figure 9). For simplicity, usually in afterglow studies it is assumed that ξ e = 1, i.e. all the particles are accelerated into a non-thermal distribution. Such 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 ξ e = 1 for parameters E k , n, e and B cannot be distinguished from those one estimated for any ξ e in the range m e /m p < ξ e < 1 and afterglow parameters E k = E k /ξ e , n = n/ξ e , e = ξ e e and B = ξ 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 Γ ∝ (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 ξ e but at the same time the energy density given to the shock is increased by a factor 1/ξ e so that the magnetic field energy density is the same in the scenario when including or not ξ e . Analogous considerations can be done also 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 ξ e or for any value m e /m p < ξ e < 1. This imply that afterglow model parameters, when considering ξ e = 1, are estimated with an uncertainty of factor m e /m p and afterglow observations do not constrain directly their values (e.g. E k or B ) but rather a fraction of their value multiplied or divided by ξ e (e.g. E k /ξ e or B ξ e ).
It is possible to obtain information on the value of ξ 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 ηγm e c 2 γm p c 2 where η m p /m e is a factor describing the ignorance on the plasma physics governing electron heating beyond γ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 η 1 and moderate 1/ξ e ) or a large self-absorption optical depth (for ξ e 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 impact also the ξ 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 ξ e in the modeling, which is constrained by data to be ξ e < 0.13 [72].

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 γ-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 has been detected also 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σ) or a hint of detection (significance between 3 and 5σ) 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 summarise the main results. Being these detections 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 address also 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 ν ∝ ν α , while temporal indices are defined by F(t) ∝ t β T .

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 ∼ 3σ excess at the GRB position.

General properties and multi-wavelength observations
The BAT prompt spectrum (T 90 = 0.48 s) is well described by a power-law with index α = −0.11 ± 0.88 and an exponential high-energy cutoff, corresponding to a peak energy with E p = 46.3 ± 6.4 keV [143]. The GBM prompt spectrum (T 90 = 1.088 ± 0.977 s) is fitted with a cutoff power-law function as well, with E p = 92 ± 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 γ,iso ∼ 1.2 × 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 behaviour 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 ± 0.1 mag 0.95 hrs 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]. The multi-wavelength light-curves of GRB 160821B are shown in Figure 11.

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 • to 55 • ). 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σ (2.9σ post-trial 4 ) significance 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 α = −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 × 10 −11 cm −2 s −1 and 5.4 × 10 −12 cm −2 s −1 . In the subsequent 2.2 h (6134 − 14130 s) assuming that the signal is real, a flux value could be calculated and is 9.9 ± 4.8 × 10 −13 cm −2 s −1 . For the same time interval, also a flux upper limit has been estimated and gives 3.0 × 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 95% confidence level following the prescriptions of [156]. The low (∼ 3σ) significance estimated did not allow to obtain an unfolded spectrum. As a result, in the SED ( Figure 12) the reconstructed flux in the third bin (6134 − 14130 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.

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 ∼ 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 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 (1TeV) SSC ∼ 2 × 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 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 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). Flux [erg/cm2/s] MAGIC 1.7-4 h 6 GHz @ 4h optical @ 2h X-ray @ 3h 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, 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 it can be seen 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 also the EIC model is unable to explain the large TeV flux suggested by MAGIC observations.

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]. telescopes observed and detected GRB 180720B about 11 h after the prompt emission, with a ∼ 5σ statistical significance.

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 ± 0.4 s and released an isotropic-equivalent energy E γ,iso = 6.0 ± 0.1 × 10 53 erg in the 50-300 keV range. Multi-wavelength afterglow observations covered the entire electromagnetic spectrum (see Figure 11). 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 × 10 3 s to 4 × 10 6 s) can be modelled with an initial power-law decay with an index −1.19 +0.01 −0.02 followed by a break at t break = 8 × 10 4 s to an index of −1.55 +0.04 −0.05 [161][162][163][164][165][166][167][168][169][170] revealed the presence of a counterpart and allowed to estimate the redshift value of z = 0.654. The optical afterglow was seen 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].

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σ pre-trial significance (5.0σ post-trial 6 ). The observation was performed in standard dark and good weather conditions with a mean zenith angle of 31.5 • . 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 E E 0 −γ int , to describe an intrinsic power-law spectrum affected by the EBL attenuation e −τ(E,z) (see Figure 15, panel on the right). With reference to the second fit, the analysis returns a photon index γ int = 1.6 ± 1.2 (statistical) ±0.4 (systematic) and a flux normalization F 0,int = (7.52 ± 2.03 (statistical) +4.53 −3.84 ) (systematic) ×10 −10 TeV −1 cm −2 s −1 , evaluated at an energy E 0,int = 0.154 TeV.

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 (α ∼ −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 Γ ∼ 1000 which at these late times is really unlikely. As a result, these strong requirements disfavour the synchrotron emission as responsible of the VHE component in GRB 180720B.  [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.
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 ∼ 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 , Γ 0 = 300 and p = 2.4. As it can be noticed, the equipartition factor B needs to assume a quite 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] ( 2 × 10 −6 erg cm −2 s −1 at t ≈ 100 s) and the one observed by Fermi-LAT (∼ 10 −8 erg cm −2 s −1 at t ≈ 100 s). A low magnetic field equipartition factor is derived from the condition E KN 0.44 TeV at t ≈ 10 h where E KN is the energy at which the KN scattering becomes relevant. A transition energy between 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 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.

General properties and multi-wavelength observations
The duration of the prompt emission is T 90 ≈ 116 s as measured by Fermi-GBM and T 90 ≈ 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 characterised by a smooth power-law decay, recorded by these instruments may be already the afterglow emission. Support to such 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 γ,iso = (2.5 ± 0.1) × 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 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 ∝ t α with α = −1.36 ± 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 ± 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. 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 • . Observations lasted until 15912 s when a zenith angle of 81.14 • was reached. The observation was performed in good weather conditions but in presence of the Moon, resulting in a night sky background approximately 6 times higher than the standard dark night conditions. The results of the offline analysis show a clear detection above the 50σ level in the first 20 minutes of observation [3].

VHE observations and results
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 β T = −1.60 ± 0.07, steeper than the one exhibited by the X-ray flux. The temporal evolution of the intrinsic spectral photon index α int of the TeV differential photon spectrum is shown in the bottom panel. A constant value of α int ≈ −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 α obs = −5.34 ± 0.22 and α int = −2.22 +0. 23 −0.25 for the observed and EBL-corrected spectrum, respectively.

Interpretation
The properties of the VHE light curve and spectrum of GRB 190114C were studied by the MAGIC Collaboration in [3]. The PL behaviour, 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 γ,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 Eq. 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 shows a double-peaked behaviour 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 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 γ − γ 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 3 × 10 53 erg, ε e ≈ 0.05 − 0.15, ε B ≈ (0.05 − 1) × 10 −3 , n ≈ 0.5 − 5 cm −3 and p ≈ 2.4 − 2.6. It is found that, being the peak of the SSC component below 200 GeV, the KN suppression and the γ − γ internal absorption have a non-negligible role in shaping the peak of the VHE spectrum. The modeling reproduces very well the XRT, LAT and TeV emission (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 × 10 53 erg, ε e = 0.6, ε B = 1 × 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 MAGIC spectra are EBL-corrected assuming the model by [183]. From [71]. 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 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 showing that at GeV energies a transition between the synchrotron and SSC component can be identified. From re-analysis of LAT data, a hard photon index (1.76 ± 0.21) is derived, in agreement with the hardening of the spectrum caused by the rising of the SSC component. Differently from what seen by [71,184] γ − γ 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 ∼ 10 −3 ). A consistent modeling of the multi-wavelength observations as synchrotron and SSC radiation is found, both the ISM and wind-like scenarios. The SED at 80 s (see figure 2 and figure 5 in [185]) and the broad-band light curves ( figure 3 and figure 6) 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 γ − γ absorption contribution are calculated and compared with the data. A wind-like environment was used for the   TeV. This implies that the KN effect is relevant only at TeV levels and the VHE data can be modelled assuming that the SSC scattering is in Thompson regime. The γ − γ 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. A smoothed analytical approximations for the electron injection function was 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) 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.

General properties and multi-wavelength observations
The prompt emission detected by the two instruments consists of two episodes with the first one seen 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 ± 0.08 and an exponential high-energy cutoff function with E p = 130 ± 20 keV, the second one can be described with a Band function with E p = 11 ± 1 keV , α = −0.92 ± 0.62 and β T = −2.51 ± 0.01 [189]. The (isotropic equivalent) prompt emission energy inferred from spectral analysis of Fermi-GBM data is E γ,iso ∼ 2 × 10 50 erg.
A multi-wavelength observational campaign of the event was performed covering the entire electromagnetic spectrum (see Figure 23 and 25). 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 × 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 behaviour with an initial steep decay phase followed by a plateau and a strong flare episode ( Figure 25, 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 ± 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 has been reported [193]. Also in the optical data a flare is seen simultaneously with the one in 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.  Figure 23 with the XRT light-curve and the LAT upper limits. The time-evolving flux was satisfactorily modeled with a power-law decay F(t) ∝ t α with α = −1.09 ± 0.05. Such decay index is similar to the X-ray one derived in the same time interval (α X = −1.07 ± 0.09).

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 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 (Γ 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 24, light blue shaded area). Possible improvements between the data and the model foresee a higher Γ, 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 case of 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 24, orange shaded area). The SSC contribution is negligible while the γ − γ absorption shapes the VHE spectrum. The single synchrotron component scenario provides better fit (> 5σ) 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).  [72]. The 90% and 50% credible intervals from the fit are shown in lighter shades.
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 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 25. 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 ξ e = 1 is excluded by the fit which provides acceptable solutions only for ξ e 6.5 × 10 −2 . Moreover, an isotropic-equivalent kinetic energy at the afterglow onset E k = 2.5 +1.9 −1.3 × 10 53 erg is estimated. Considering the observed GBM prompt energy, such value implies that the prompt efficiency is η = 1.2 +1.0 −0.5 × 10 −3 which is much lower than the typical values derived from 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 model proposes that the GRB jet is seen off-axis (θ view = 1.78 • ) and it consists of a narrow (θ jet = 0.86 • ) fast (Γ = 350) jet and a slow (Γ = 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 done 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 seen 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 ∼ 10 3 − 10 4 s the EIC component dominates the VHE emission, while for later times (t 3 × 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.   The (isotropic equivalent) prompt emission energy inferred from spectral analysis of Fermi-GBM data is E γ,iso = (1.1 ± 0.2) × 10 50 erg [202]. The prompt duration is T 90 = 9.78 ± 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 photon index β T = −3.03 ± 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 α = −1.49 +0. 24 −0.21 . 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 × 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. Final results from VHE data analysis of GRB 201015A have not been published yet. Preliminary information reported here have been released in [217] and [216]. 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 • up to 48 • , 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σ 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 γ-ray data.

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 second before the Swift-BAT) [219]. MAGIC detected GRB 201216C with a significance of ∼ 6σ.

General properties and milti-wavelength observations
The duration is estimated as T 90 = 48 ± 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 ± 7 keV, α = −1.06 ± 0.01 and β T = −2.25 ± 0.03. The isotropic equivalent energy E γ,iso in the 10-1000 keV band is (4.71 ± 0.16) × 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 behaviour with decay indices of 1.97 +0.10 −0.09 and 1.07 +0.15 −0.10 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 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]. Final results from VHE data analysis of GRB 201216C have not been published yet. Preliminary information reported here have been released in [229] and [228]. 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 Moon. The zenith angle ranged from 37 • to 68 • . The low level of night sky background allowed to retain also the low energy events and therefore obtain a low energy threshold with 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.

VHE observations and results
The signal significance was calculated to be 6.0σ pre-trial (5.9σ post-trial 7 ) for the first 20 minutes of observation (see Figure 27). A preliminary time-integrated spectrum for the first 20 minutes 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 hundreds 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 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 since no significant emission was found after this time. The preliminary results are consistent with a monotonically decaying light curve fitted with a power-law.

The new TeV spectral window: discussion
After decades of searches, MAGIC and H.E.S.S. observations have unequivocally proved 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 includes also two events (namely GRB 160821B and GRB 201015A) where only a hint of excess (i.e., with a significance at ∼ 3 − 4σ) was found. For the other four events, namely GRB 180720B, GRB 190114C, GRB 190829A and GRB 201216C the detections are robust (> 5σ). The table lists several properties, such as duration T 90 and total emitted energy E γ,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 behaviour which 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 a quite 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.  Table 4: List of the GRBs observed by IACTs with a firm detection (significance > 5σ) or a hint of detection (3 − 4σ) above 100 GeV. The T 90 and E γ,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.

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 • and in 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 detect GRBs at higher redshift, for which the impact of EBL is large already at a few hundreds GeV. This explains how it has been possible the detection of a GRB at z = 1.1 (GRB 201216C): in this case, optimal observing conditions allowed 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.

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 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 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 a 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 explore more robustly the spectral shape, unbiased by the EBL effect which is a non-negligible source of uncertainty for higher redshifts.

Energetics
In terms of E γ,iso the VHE GRB sample spans more than three orders of magnitude from ∼ 10 49 erg up to ∼ 6 × 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 γ,iso ∼ 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].

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 ∼ 10 4 s their luminosity is around 1-5 ∼ 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 order of magnitude at t ∼ 10 4 s). This is consistent with the fact that they also have a smaller E γ,iso . The correlation between X-ray afterglow luminosity and prompt E γ,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 similarity. 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 ∼ 380 s the VHE luminosity L V HE is ∼ 1.5 − 2.5 × 10 48 erg s −1 and the X-ray one L X is ∼ 0.6 − 1.0 × 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 ∼ 4.5 hrs the VHE luminosity L V HE is ∼ 4.0 − 8.5 × 10 44 erg s −1 which is around ∼ 15 − 20% of the corresponding X-ray one (L X ∼ 2.0 − 5.0 × 10 45 erg s −1 ). For GRB 180720B at t ∼ 2 × 10 4 s the VHE luminosity L V HE is ∼ 9 × 10 47 erg s −1 and the X-ray one L X is ∼ 1.5 − 2.5 × 10 49 erg s −1 . In this case the power radiated in the VHE band is around ∼ 35 − 60% of the X-ray one.

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 seen 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. An 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 24). In this GRB an interpretation of the whole SED in terms of synchrotron radiation can not be excluded, although a modeling as SSC radiation has been proved to be successful [72] ( Figure  25). For the other events detected at VHE, either the data do not allow to build 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 modeling have suggested the possibility that the fraction of electrons accelerated in non-thermal distribution, ξ e , is different from the standard value of 1 which is usually assumed. In a few GRBs modeling, namely for GRB 190114C [185] and GRB 190829A [72,198] the introduction of a ξ e < 1 was essential in order to fit consistently the observational data. In particular in [72] the requirement for a low value of ξ e 6.5 × 10 −2 was required to provide acceptable fit of the data. The other modeling assume a greater value of ξ e , around ∼ 0.3. Further detections will be exploited in order to verify if such 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 order 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 external shock scenario, are excluded. Moreover, some results also can 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 which reproduce the early time data over-predict the late time optical and radio observations. This result point 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 ∝ t −0.4 and B ∝ t 0.1 in the ISM case and B ∝ 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 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 show quite broad properties, spanning more than three order of magnitude in E γ,iso and more than two order 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 ISM or wind-like scenario and indicate that shock micro-physics parameters which seems to be able to reproduce VHE emission are ε e ∼ 0.1 and ε B ∼ 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 detection except for the hint of excesses seen for GRB 160821B.

Conclusions and future prospects
The recent discoveries performed by 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 sample have shown the potential that such detections have in probing several longstanding 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 behaviour of the VHE component with respect to emission at lower frequencies. Modeling of multi-wavelength data covering from radio up to TeV energies were 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 (∼ 10 4 times larger than Fermi-LAT at 30 GeV) and iii) a rapid slewing capability (180 degrees azimuthal rotation in 20 seconds). 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 commissioning phase 10 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 anyway influenced 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 localisation (later refined by XRT to a few arcsec) is fundamental for Cherenkov telescopes, given their limited field of view (e.g., about 4 • for the LSTs and 7 • 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 • up to 10 • , which makes follow-up with IACTs very challenging. In case of such a 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 • , including also 10 GRBs/yr with redshift z < 1.
Available estimates of 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 few seconds to hours, and has been extrapolated to 10 TeV on the bases of model predictions. The results show that the instrument will be able to detect afterglow TeV emission from an event like GRB 190114C up to ∼ 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 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 characterised 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, 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 if 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.