###### Abstract

Afterglow radiation in gamma-ray bursts (GRB), extending from the radio band
to GeV energies, is produced as a result of the interaction between the
relativistic jet and the ambient medium. Although in general the origin of the
emission is robustly identified as synchrotron radiation from the shock-
accelerated electrons, many aspects remain poorly constrained, such as the
role of inverse Compton emission, the particle acceleration mechanism, the
properties of the environment and of the GRB jet itself. The extension of the
afterglow emission into the TeV band has been discussed and theorized for
years, but has eluded for a long time the observations. Recently the Cherenkov
telescopes MAGIC and H.E.S.S. have unequivocally proven that afterglow
radiation is 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.

###### keywords:

gamma-ray bursts; non-thermal processes;Cherenkov telescopes

100 1 1 0 https://doi.org/ Gamma-ray bursts afterglow physics and the VHE
domain GRB afterglow and the VHE domain Davide Miceli 1,2,∗, Lara Nava 3,4,∗
Davide Miceli, Lara Nava Miceli, D.; Nava, L. Correspondence:
<EMAIL_ADDRESS>(D.M<EMAIL_ADDRESS>(L.N.)

##  1 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\sim 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 Kouveliotou et al. (1993);
Mazets et al. (1981). 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 $\gamma$-rays to the radio band. The afterglow flux
decays smoothly as a power-law in time for weeks or months, and the typical
frequency of the radiation moves in time from the X-ray to the radio band.
Since 2019, the detection of a few long GRBs between 0.3 and 3 TeV on time-
scales from tens of seconds to a few days proved for the first time that GRBs
can also be sources of radiation in the TeV band, where they can convey a
sizable fraction (20-50%) of the total energy emitted during the afterglow
phase MAGIC Collaboration et al. (2019); Abdalla et al. (2019); H. E. S. S.
Collaboration et al. (2021).

All this prompt/afterglow emission is identified with radiation produced as a
result of the launch of an ultra-relativistic ($\Gamma\sim 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 Rees
and Meszaros (1994) or magnetic reconnection episodes Thompson (1994); Spruit
et al. ). In a second moment, the ejecta undergoes external dissipation Rees
and Meszaros (1992), triggered by interactions with the ambient medium (e.g.,
the interstellar medium or the wind of the progenitor’s star Chevalier and Li
(1999)). The two different dissipation processes occur at different typical
distances from the central engine ($R\sim 10^{13-14}$ cm and $R\sim
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 Woosley (1993); Paczyński (1998). In this model, the
core of a massive star collapses into a black hole and the accretion from the
surrounding disk powers the launch of two opposite, collimated
($\theta_{jet}\sim 5-10^{\circ}$) 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 Abbott et al. (2017)) or a neutron star and a black-
hole. An alternative model Usov (1992); Thompson et al. (2010); Metzger et al.
(2011); Dall’Osso et al. (2011) 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 Usov (1992); Thompson (1994)) 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 Usov (1992). 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
Giannios and Spruit (2005); Zhang and Yan (2011); Lazarian et al. (2019). In
particular, internal shocks are not an efficient mechanism Lazzati et al.
(1999); Sironi et al. (2015), 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 Fan
and Piran (2006). 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%
Beniamini et al. (2016), 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 Rees and Meszaros (1994); Ghisellini et al. (2000).
Inconsistencies between the expected synchrotron spectrum and the observed
spectral shape of the prompt emission Preece et al. (1998); Ghisellini et al.
(2000) 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 Oganesyan et al. (2017, 2018, 2019); Ravasio et al.
(2018, 2019).

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 Meszaros and Rees (1993); Mészáros and
Rees (1997); Sari and Piran (1999). These shocks are responsible for the
acceleration of particles and for the deceleration of the outflow, eventually
down to non-relativistic velocities Blandford and McKee (1976); Sedov (1946);
Taylor (1950). The observed radiation is the result of synchrotron radiation
from electrons accelerated at the forward shock Sari, Re’em et al. (1998). A
contribution from the reverse shock may also be relevant, typically in the
radio and optical band Sari and Piran (1999); Kobayashi and Sari (2000). 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 the theoretical side (see Vanthieghem et al. (2020) for a recent
review), especially thanks to numerical particle-in-cell (PIC) simulations
Sironi and Spitkovsky (2011); Sironi et al. (2013). Ultra-relativistic shocks
in a weakly magnetised medium are found to be efficient particle accelerators,
with $\xi_{e}\sim 1\%$ of the electrons being accelerated into a power-law
distribution with spectral index $p\sim 2.5$, carrying about $\epsilon_{e}\sim
10\%$ of the shock-dissipated energy. A strong magnetic turbulence is self-
generated by the accelerated particles counter-streaming in the upstream,
ahead of the shock, at a level of magnetization $\epsilon_{B}=0.01-0.1$ Sironi
et al. (2013). PIC simulations however are currently probing time-scales that
are orders of magnitude smaller than the dynamical time-scale of the blast-
wave. This implies that results from simulations can only be extrapolated to
the relevant time-scales, introducing a certain degree of uncertainty and
caution in using the results as inputs for the modeling of GRB afterglows.
What is still poorly understood, even though dedicated simulations are
starting to give important clues Keshet et al. (2009), is how the micro-
turbulence generated in the shock vicinity evolves (decays) with time. This is
particularly important for a proper interpretation of the observations, since
it is likely that the particles produce synchrotron and synchrotron-self
Compton (SSC) photons in a region of decayed micro-turbulence, and hence feel
a magnetic field with $\epsilon_{B}<<0.01$.

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

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

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

For recent and complete reviews on GRB’s phenomenology and theoretical
interpretation before the TeV era see Kumar and Zhang (2015); Zhang (2018). An
overview focused on high-energy emission (0.1-100 GeV) observations and
interpretation can be found in Nava (2018).

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

##  2 The afterglow model

Afterglow emission refers to all the broad-band radiation observed from a GRB
on longer timescales (minutes to months) as compared to the initial prompt
radiation detected in hard X-rays Panaitescu and Kumar (2000); Sari, Re’em et
al. (1998); Sari, Re’em and Esin (2001). 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 Nousek
et al. (2006); Nardini et al. (2006); Zaninoni et al. (2013); Bernardini et
al. (2012). 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 Sari and Piran
(1995). 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).

###  2.1 Jet dynamics

After the reverse shock has crossed the ejecta, the dynamics of the blast-wave
enters a self-similar regime (Blandford and McKee (1976), 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 ($\sim m(R_{dec})c^{2}\Gamma_{0}^{2}$) is comparable to
the initial energy ($E_{0}=M_{0}\Gamma_{0}\,c^{2}$) carried by the jet. This
deceleration radius is typically of the order of $R_{dec}\sim 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 $\Gamma_{0}$. Before reaching this
radius, the ejecta expands with constant velocity (coasting phase).

Most analytic estimates of the afterglow evolution with the purpose of
modeling data are developed for the deceleration phase, where the self-similar
BM76 solution for adiabatic blast-waves is adopted Sari, Re’em et al. (1998);
Panaitescu and Kumar (2000); Granot and Sari (2002). 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 Nava et al. (2013). This method allows to describe
the hydrodynamics of a relativistic blast-wave expanding into a medium with an
arbitrary density profile $\rho(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 $\Gamma$ 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
$\epsilon_{\rm 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 Nava et al. (2013).

####  2.1.1 Equation describing the evolution of the bulk Lorentz factor

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

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

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

$E=\Gamma Mc^{2}+\Gamma_{\rm eff}E^{\prime}_{\rm int}~{},$ (1)

where:

$\Gamma_{\it\it
eff}\equiv\frac{\hat{\gamma}\Gamma^{2}-\hat{\gamma}+1}{\Gamma}~{},$ (2)

which properly describes the Lorentz transformation of the internal energy.
Here, $M=M_{0}+m=\rho^{\prime}V^{\prime}$ is the sum of the ejecta mass
$M_{0}=E_{0}/\Gamma_{0}c^{2}$ and of the swept-up mass $m(R)$, and
$E^{\prime}_{int}=(u^{\prime}-\rho^{\prime}c^{2})V^{\prime}$ is the comoving
internal energy. The adiabatic index can be parameterized as
$\hat{\gamma}=(4+\Gamma^{-1})/3$ to obtain the expected limits
$\hat{\gamma}\simeq 4/3$ for $\Gamma\gg 1$ and $\hat{\gamma}\simeq 5/3$ for
$\Gamma\rightarrow 1$. The majority of analytical treatments use $\Gamma$
instead of $\Gamma_{\it eff}$, which implies an error up to a factor of $4/3$
in the ultra-relativistic limit Nava et al. (2013).

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}=\Gamma_{\it eff}\,dE^{\prime}_{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:

$d\left[\Gamma(M_{0}+m)c^{2}+\Gamma_{\it
eff}E^{\prime}_{int}\right]=dm\,c^{2}+\Gamma_{\it eff}dE^{\prime}_{rad}~{}.$
(3)

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

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

The first contribution, $dE^{\prime}_{sh}=(\Gamma-1)\,dm\,c^{2}$, is the
random kinetic energy produced at the shock as a result of the interaction
with an element $dm$ of circum-burst material: as pointed out by BM76, in the
post-shock frame, the average kinetic energy per unit mass
$dE^{\prime}_{sh}/dm$ is constant across the shock, and equal to
$(\Gamma-1)c^{2}$. The second term in eq. 4, $dE^{\prime}_{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^{\prime}_{rad}$, accounts for radiative losses.

From Eq. 3, it follows that the variation of the Lorentz factor is:

$\frac{d\Gamma}{dR}=-\frac{(\Gamma_{\it
eff}+1)(\Gamma-1)\,c^{2}\frac{dm}{dR}+\Gamma_{\it
eff}\frac{dE^{\prime}_{ad}}{dR}}{(M_{0}+m)\,c^{2}+E^{\prime}_{int}\frac{d\Gamma_{\it
eff}}{d\Gamma}}~{},$ (5)

from which the evolution $\Gamma(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 $\Gamma_{\rm eff}\,dE^{\prime}_{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 $\rho\propto R^{-3}$. To evaluate Eq. 5 it is necessary
to first specify $dE^{\prime}_{ad}$ and $E^{\prime}_{int}$.

####  2.1.2 Internal energy and adiabatic losses

In specific cases, the adiabatic losses and the internal energy content can be
expressed in an analytic form. The following treatment to estimate adiabatic
losses and the internal energy content of the blast-wave assumes that, right
behind the shock, the freshly shocked electrons instantaneously radiate a
fraction $\epsilon_{\rm rad}$ of their internal energy and then they cool only
due to adiabatic losses Nava et al. (2013). 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 $\epsilon_{\rm e}$ as the fraction of energy $dE^{\prime}_{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):

$\gamma_{acc,e}-1=(\Gamma-1)\epsilon_{e}/\mu_{e}\,.$ (6)

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

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

$\gamma_{rad,e}-1=(1-\epsilon_{rad})(\gamma_{acc,e}-1)=(1-\epsilon_{rad})(\Gamma-1)\frac{\epsilon_{e}}{\mu_{e}}~{}.$
(7)

The assumption of instantaneous radiative losses is verified in the fast
cooling regime ($\epsilon_{rad}\sim 1$), which is required (but not
sufficient) to have $\epsilon\sim 1$ (i.e., a fully radiative blast-wave). In
the opposite case $\epsilon_{rad}\ll 1$, the evolution is nearly adiabatic
($\epsilon\ll 1$), regardless of the value of $\epsilon_{\rm e}$, and the
details of the radiative cooling processes are likely to be unimportant for
the shell dynamics. The case with intermediate values of $\epsilon_{\rm rad}$
and $\epsilon$ 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
$\epsilon_{\rm p}$ of the energy dissipated by the shock (with
$\epsilon_{p}=1-\epsilon_{e}-\epsilon_{B}$), their mean post–shock Lorentz
factor will be:

$\gamma_{acc,p}-1=(\Gamma-1)\frac{\epsilon_{p}}{\mu_{p}}~{},$ (8)

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

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

$u^{\prime}_{\rm int}=\frac{N(\gamma-1)m\,c^{2}}{V^{\prime}}~{}.$ (9)

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

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

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

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

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

$\gamma_{\rm
ad}(R,r)-1=\left(\frac{r}{R}\right)^{3(\hat{\gamma}-1)}\left[\frac{\Gamma(R)}{\Gamma(r)}\right]^{(\hat{\gamma}-1)}(\gamma(r)-1)~{}.$
(12)

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

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

$E^{\prime}_{int}(R)\\!=\\!4\pi
c^{2}\\!\\!\\!\int_{0}^{R}\\!\\!\\!drr^{2}\left\\{\rho_{p}(r)[\gamma_{ad,p}(R,r)\\!-\\!1]\\!+\\!\rho_{e}(r)[\gamma_{ad,e}(R,r)\\!-\\!1]\right\\}~{}.$
(13)

With the help of Eq. 12, one can explicitly find $E^{\prime}_{int}(R)$ and
insert it in Eq. 5.

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

$\frac{dE^{\prime}_{ad}(R)}{dR}=-4\pi
c^{2}(\hat{\gamma}-1)\left(\frac{3}{R}-\frac{d\log\Gamma}{dR}\right)\int_{0}^{R}drr^{2}\left\\{\rho_{p}(r)(\gamma_{ad,p}-1)+\rho_{e}(r)(\gamma_{ad,e}-1)\right\\}~{}.$
(14)

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:

$\begin{split}E^{\prime}_{int}(R)&=4\pi
c^{2}\int_{0}^{R}\,dr\,r^{2}\frac{r}{R}\left[\frac{\Gamma(R)}{\Gamma(r)}\right]^{1/3}\Gamma(r)\left\\{\frac{\epsilon_{p}}{\mu_{p}}\rho_{p}+(1-\epsilon_{rad})\frac{\epsilon_{e}}{\mu_{e}}\rho_{e}\right\\}\\\
&=4\pi
c^{2}\int_{0}^{R}\,dr\,r^{2}\frac{r}{R}\left[\frac{\Gamma(R)}{\Gamma(r)}\right]^{1/3}\Gamma(r)\,\rho(r)\left(\epsilon_{p}+\epsilon_{e}-\epsilon\right)\\\
&\frac{dE^{\prime}_{ad}(R)}{dR}=-E^{\prime}_{int}(R)\left(\frac{1}{R}-\frac{1}{3}\frac{d\log\Gamma}{dR}\right)~{}.\end{split}$
(15)

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

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

$E^{\prime}_{int}(R)=\frac{9-3s}{9-2s}\Gamma
mc^{2}~{}~{},~{}~{}~{}~{}~{}~{}\frac{dE^{\prime}_{ad}(R)}{dR}=-\frac{(9-s)(3-s)}{2(9-2s)}\frac{\Gamma
mc^{2}}{R}~{},$ (16)

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

In the fully radiative regime $\epsilon=1$, which implies $E^{\prime}_{int}=0$
and $dE^{\prime}_{ad}=0$, Eq. 5 reduces to:

$\frac{d\Gamma}{dm}=-\frac{(\Gamma_{\it eff}+1)(\Gamma-1)}{M_{0}+m}~{}~{},$
(17)

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

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

$E_{0}\simeq\Gamma_{\it eff}E^{\prime}_{int}\simeq\frac{4}{3}\Gamma
E^{\prime}_{int}\simeq\frac{12-4s}{9-2s}\Gamma^{2}mc^{2}~{}~{},$ (18)

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:

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

No analytic model properly captures the transition between an adiabatic
relativistic blast-wave and the momentum–conserving snowplow, as $\epsilon$
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
$\Gamma$ as a function of the shock radius $R$.

###  2.2 Relativistic shock acceleration

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

In this section we 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
$\gamma_{min/max}$ and the spectral index $p$), ii) the acceleration
efficiency (i.e., the fraction of electrons $\xi_{e}$ and the fraction of
energy $\epsilon_{e}$ in the non-thermal population), and iii) the strength of
the self-generated magnetic field, usually quantified in terms of fraction
$\epsilon_{B}$ of the shock-dissipated energy conveyed in the magnetic field.
In particular, in order to compare with observations, the relevant
$\epsilon_{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 Sironi et al. (2015); Vanthieghem et al. (2020)), we
discuss how particle acceleration and magnetic field amplification are
incorporated in GRB afterglow modeling, and then we comment on the constraints
on the above-mentioned parameters as inferred from the comparison between the
model and the observations.

####  2.2.1 Inputs from theoretical investigations

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

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

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
Keshet and Waxman (2005) 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 Keshet and Waxman (2005) for
details): in all cases, the value of the spectral index approaches the same
value $s_{p}\sim 4.2$ (corresponding to $p\sim 2.2$) in the ultra-relativistic
limit.

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

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

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

Figure 2: PIC simulations: temporal evolution of the downstream spectrum of
ions (upper panel) and electrons (bottom) for mass ratio $m_{i}/m_{e}=25$,
shock Lorentz factor $\Gamma=15$ and magnetization $\sigma=10^{-5}$. The
evolution is followed until $t=2500\omega_{pi}^{-1}$. From Sironi et al.
(2013). ©AAS. Reproduced with permission.

An example of downstream particle spectra derived by PIC simulations is shown
in Figure 2 (Sironi et al. (2013)). The ion and electron spectra are shown for
a 2D simulation with $\Gamma=15$, ion-to-electron mass ratio $m_{i}/m_{e}=25$,
and $\sigma=10^{-5}$. The temporal evolution is followed up to
$t=2500\,\omega_{pi}^{-1}$. The formation of a non-thermal tail is clearly
visible.

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

The maximum energy $\gamma_{max}$ increases proportionally to $t^{1/2}$ (see
inset in Figure 2), slower than the commonly adopted Bhom rate Böhm (1949), in
which case $\gamma_{max}\propto t$. Extrapolating the $\gamma_{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 $\gamma_{max}\sim 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
$\epsilon_{B}\sim 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 close to the shock and in the region where particle cooling
takes place. Electrons would then cool in a region of weak magnetic field
Lemoine (2013); Lemoine et al. (2013). These considerations suggest that it
might not be correct to define a single magnetization $\epsilon_{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.

####  2.2.2 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
$\epsilon_{e}$ and $\epsilon_{B}$ of dissipated energy gained by the
accelerated particles and amplified magnetic field, the spectral index $p$ of
the accelerated particle spectrum, and the fraction $\xi_{e}$ of particles
which efficiently enter the Fermi mechanism and populate the non-thermal
distribution.

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

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

which is valid in both the ultra-relativistic and non-relativistic limits (see
e. g., Blandford and McKee (1976)). In the GRB afterglow scenario it is
usually assumed that pairs are unimportant and then the density of protons and
electrons is the same: $n_{p}=n_{e}=n$. This implies that the mass is
dominated by protons: $\rho=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:

$u^{\prime}_{\rm sh}=4\Gamma(\Gamma-1)nm_{p}c^{2}$ (21)

A fraction $\epsilon_{B}$ of this energy will be conveyed to the magnetic
field:

$u^{\prime}_{B}=\epsilon_{B}u^{\prime}_{\rm
sh}=\epsilon_{B}4\Gamma(\Gamma-1)nm_{p}c^{2}$ (22)

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

$B^{\prime}=\sqrt{32\pi\epsilon_{B}m_{p}c^{2}n(\Gamma-1)\Gamma}$ (23)

Similarly, for the accelerated electrons:

$u^{\prime}_{e}=\epsilon_{e}u^{\prime}_{sh}=\epsilon_{e}4\Gamma(\Gamma-1)nm_{p}c^{2}=\langle\gamma\rangle
m_{e}c^{2}4\Gamma\xi_{e}n$ (24)

where $\langle\gamma\rangle$ is the average random Lorentz of the accelerated
electrons:

$\langle\gamma\rangle=\frac{\epsilon_{e}}{\xi_{e}}\frac{m_{p}}{m_{e}}(\Gamma-1)$
(25)

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

$\langle\gamma\rangle=\frac{\int_{\gamma_{min}}^{\gamma_{max}}N(\gamma)\gamma
d\gamma}{\int_{\gamma_{min}}^{\gamma_{max}}N(\gamma)d\gamma}$ (26)

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

$\begin{array}[]{cc}\Bigg{[}\frac{\ln\big{(}{\frac{\gamma_{max}}{\gamma_{min}}}\big{)}}{\frac{1}{\gamma_{min}}-\frac{1}{\gamma_{max}}}\Bigg{]}=\frac{\epsilon_{e}}{\xi_{e}}\frac{m_{p}}{m_{e}}(\Gamma-1)\hskip
88.2037pt\text{if $p=2$}\\\
\Bigg{[}\frac{p-1}{p-2}\frac{\gamma_{min}^{-p+2}-\gamma_{max}^{-p+2}}{\gamma_{min}^{-p+1}-\gamma_{max}^{-p+1}}\Bigg{]}=\frac{\epsilon_{e}}{\xi_{e}}\frac{m_{p}}{m_{e}}(\Gamma-1)\hskip
56.9055pt\text{if $p\neq 2$}\end{array}$ (27)

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

$\gamma_{min}=\frac{\epsilon_{e}}{\xi_{e}}\frac{m_{p}}{m_{e}}\frac{p-2}{p-1}(\Gamma-1)$
(28)

Since $p$ is expected to be $2<p<3$, this condition is verified for
$\gamma_{max}\gg\gamma_{min}$. The minimum Lorentz factor is then not treated
as a free parameter of the model, as it is calculated from eq. 28 as a
function of the free parameters $\epsilon_{e},\xi_{e}$ and $p$. Concerning the
prescription for the value of $\gamma_{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
$\epsilon_{e}$ with $\epsilon_{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:

$\begin{array}[]{cc}\Bigg{[}\frac{\ln\big{(}{\frac{\gamma_{max,p}}{\gamma_{min,p}}}\big{)}}{\frac{1}{\gamma_{min,p}}-\frac{1}{\gamma_{max,p}}}\Bigg{]}=\frac{\epsilon_{p}}{\xi_{p}}(\Gamma-1)\hskip
79.6678pt\text{if $q=2$}\\\
\Bigg{[}\frac{q-1}{q-2}\frac{\gamma_{min,p}^{-q+2}-\gamma_{max_{p}}^{-q+2}}{\gamma_{min,p}^{-q+1}-\gamma_{max,p}^{-q+1}}\Bigg{]}=\frac{\epsilon_{p}}{\xi_{p}}(\Gamma-1)\hskip
56.9055pt\text{if $q\neq 2$}\end{array}$ (29)

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

$\gamma_{min,p}=\frac{\epsilon_{p}}{\xi_{p}}\frac{q-2}{q-1}(\Gamma-1)$ (30)

The equations for $\gamma_{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 $\Gamma_{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.

####  2.2.3 Constraints to the acceleration mechanism provided by observations

Assuming that accelerated particles have a power-law spectrum ($dN^{\rm
acc}/d\gamma\propto\gamma^{-p}$) and the cooling is dominated by synchrotron
radiation, the spectral slope $p$ can be inferred from observations of the
synchrotron spectrum and/or from the temporal decay of the lightcurves if
observations are performed at frequencies higher than the typical frequency
$\nu_{\rm m}$ of photons emitted by electrons with Lorentz factor
$\gamma_{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\sim 2$ to $p\sim 3$, suggesting that the spectrum of injected
particles does not seem to have a typical slope, at odds with theoretical
predictions. The determination of $p$ however, suffers from the uncertainties
on the spectral index inferred from optical and X-ray observations, where the
observed spectra are subject to unknown dust and metal absorption. A
derivation of $p$ from the decay rate of the lightcurves is also subject to
the correct identification of the spectra regime, and partially also to the
assumption on density profile of the external medium, which is often
unconstrained (see section 3.2).

The typical value of $\epsilon_{\rm 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, $\epsilon_{e}$ is perhaps the most well constrained
parameter of the model, and is in good agreement with values predicted by
numerical investigations Sironi et al. (2013). For the fraction $\epsilon_{\rm
B}$, on the contrary, the inferred values varies in a very wide range,
typically from $10^{-5}$ to $10^{-1}$ Santana et al. (2014); Wang et al.
(2015); Zhang et al. (2015). Recent studies that incorporate Fermi-LAT GeV
observations Beniamini et al. (2015, 2016) have shown that the typical values
estimated for $\epsilon_{B}$ can be even smaller, in the range $\sim
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 $\epsilon_{B}$ needed to model the radiation have
been tentatively interpreted as the sign of turbulence decay in the downstream
Lemoine (2013); Lemoine et al. (2013). As a consequence, even though the
turbulence is strong ($\epsilon_{B}\simeq 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 $\epsilon_{B}$
are confirmed by the modeling of recent TeV detections of afterglow radiation
from GRBs (et al. (2019); Salafia et al. (2021), see Section 4).

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

Observations so far have not been able to identify the location of a high-
energy cutoff in the synchrotron spectrum, that would reveal the maximal
energy of the synchrotron photons and then the maximum energy $\gamma_{max}$
of the accelerated electrons. Observations by Fermi-LAT are in general
consistent with a single power-law extending up to at least 1 GeV. Photons
with energies in excess of 1 GeV have been detected from several GRBs, the
record holder for Fermi-LAT being a $95$ GeV photon Ackermann et al. (2014).
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}\simeq r_{L}/c$) with
$r_{L}=E/eB$ being the Larmor radius, and is limited by synchrotron cooling
($t_{syn}\simeq 6\,\pi\,m_{e}c/\sigma_{T}B^{2}\gamma$) then $\gamma_{max}\sim
10^{7}-10^{8}$ can be reached. Even though this does not necessarily imply
that acceleration must proceed at the Bhom limit, the value of $\gamma_{max}$
inferred from the detection of GeV photons is quite large and barely
consistent with what found by PIC simulations. Whether or not the observations
are in tension with the present derivation of $\gamma_{max}$ from PIC
simulations and theoretical arguments, strongly depends on a clear
identification of the origin of photons in the GeV-TeV energy range. Present
and future observations with Imaging Atmospheric Cherenkov Telescopes (IACTs)
are the main candidates to shed light on this issue.

###  2.3 Derivation of the radiative output

The expected radiative output can be estimated by means of analytical
approximations, which provide prescriptions for the location of the
synchrotron self-absorption frequency $\nu_{sa}$, the characteristic frequency
$\nu_{m}$ emitted by electrons with Lorentz factor $\gamma_{min}$, the cooling
frequency $\nu_{c}$ emitted by electrons with Lorentz factor $\gamma_{c}$, and
the overall synchrotron flux Sari, Re’em et al. (1998); Panaitescu and Kumar
(2000). 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 Granot and Sari (2002). The associated SSC component in Thomson
regime Sari, Re’em and Esin (2001) and corrections to be applied to the
synchrotron and SSC spectra to account for the effects of Klein-Nishina Nakar
et al. (2009) 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 Blandford and McKee (1976) 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
($\epsilon_{e}$) and by the amplified magnetic field ($\epsilon_{B}$), and the
spectral slope of the accelerated electrons $p$. During the deceleration
phase, the initial bulk Lorentz factor $\Gamma_{0}$ does not play any role,
but its value determines the radius (or time) at which the deceleration
begins.

An alternative approach to estimate the expected spectra and their evolution
in time consists in solving numerically the differential equation describing
the evolution of the particle spectra and estimate the associated emission
Petropoulou and Mastichiadis (2009); Pennanen et al. (2014); Bošnjak et al.
(2009). 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 et
al. (2019).

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

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

where $\dot{\gamma}=\frac{\partial\gamma}{\partial t^{\prime}}$ is the rate of
change of the Lorentz factor $\gamma$ of an electron caused by adiabatic,
synchrotron and SSC losses and 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
MAGIChttps://magic.mpp.mpg.de/ and H.E.S.S.https://www.mpi-
hd.mpg.de/hfm/HESS/. 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(\gamma,t^{\prime})=Q^{acc}(\gamma,t^{\prime})+Q^{pp}(\gamma,t^{\prime})$
describes the injection of freshly accelerated particles
($Q^{acc}(\gamma,t^{\prime})=dN^{acc}/d\gamma\,dt^{\prime}$) and the injection
of pairs $Q^{pp}(\gamma,t^{\prime})$ 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 Chang and Cooper (1970) can be adopted.

####  2.3.1 Synchrotron and SSC cooling

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

$\dot{\gamma}_{syn}\equiv\frac{d\gamma}{dt^{\prime}}\bigg{|}_{syn}=-\frac{\sigma_{T}\gamma^{2}B^{\prime
2}}{6\,\pi\,m_{e}c}$ (32)

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

$K(\gamma,\nu^{\prime},\tilde{\nu^{\prime}})=\left\\{\begin{array}[]{ll}\frac{\varepsilon}{\tilde{\varepsilon}}-\frac{1}{4\gamma^{2}}\hskip
173.56198pt\frac{\tilde{\varepsilon}}{4\gamma^{2}}<\varepsilon<\tilde{\varepsilon}\\\
2q\ln{q}+(1+2q)(1-q)+\frac{1}{2}(1-q)\frac{(4\gamma\tilde{\varepsilon}q)^{2}}{(1+4\gamma\tilde{\varepsilon}q)}\hskip
5.69046pt\tilde{\varepsilon}<\varepsilon<\frac{4\gamma^{2}\tilde{\varepsilon}}{1+4\gamma\tilde{\varepsilon}}~{},\end{array}\right.$
(33)

where:

$\tilde{\varepsilon}=\frac{h\tilde{\nu^{\prime}}}{m_{e}c^{2}}\hskip
14.22636pt\varepsilon=\frac{h\nu^{\prime}}{m_{e}c^{2}}\hskip
14.22636ptq=\frac{\varepsilon}{4\gamma\tilde{\varepsilon}(\gamma-\varepsilon)}~{}.$
(34)

$\tilde{\varepsilon}$ and $\varepsilon$ 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. $\varepsilon<\tilde{\varepsilon}$) and the up-scattering
(i.e. $\varepsilon>\tilde{\varepsilon}$) process. The energy loss term for the
SSC can now be calculated with the equation:

$\dot{\gamma}_{SSC}=\frac{d\gamma}{dt^{\prime}}\bigg{|}_{SSC}=-\frac{3h\sigma_{t}}{4m_{e}c\gamma^{2}}\int
d\nu^{\prime}\nu^{\prime}\int\frac{d\tilde{\nu^{\prime}}}{\tilde{\nu^{\prime}}}n_{\tilde{\nu^{\prime}}}(t^{\prime})K(\gamma,\nu^{\prime},\tilde{\nu^{\prime}})~{}.$
(35)

####  2.3.2 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\gamma$ in a comoving time $dt^{\prime}$:

$\dot{\gamma}_{ad}=\frac{d\gamma}{dt^{\prime}}\bigg{|}_{ad}=-\frac{\gamma\beta^{2}}{3}\frac{d\ln
V^{\prime}}{dt^{\prime}}~{},$ (36)

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

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

and:

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

####  2.3.3 Synchrotron self-absorption (SSA)

Electrons can re-absorb low energy photons before they escape from the source
region. The absorption coefficient $\alpha_{\nu}$ can be expressed as Rybicki
and Lightman (1979):

$\alpha_{\nu}=-\frac{1}{8\pi\nu^{\prime 2}m_{e}}\int d\gamma
P^{\prime}(\gamma,\nu^{\prime})\gamma^{2}\frac{\partial}{\partial\gamma}\bigg{[}\frac{N(\gamma)}{\gamma^{2}}\bigg{]}$
(39)

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

$P^{\prime}_{\nu}\propto\left\\{\begin{array}[]{ll}\nu^{\prime 5/2}\hskip
76.82234pt\nu^{\prime}_{i}<\nu^{\prime}<\nu^{\prime}_{SSA}\\\ \nu^{\prime
2}\hskip
85.35826pt\nu^{\prime}<\nu^{\prime}_{SSA}<\nu^{\prime}_{i}~{},\end{array}\right.$
(40)

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

####  2.3.4 Synchrotron and Inverse Compton emission

Following Ghisellini et al. (1988), the synchrotron spectrum emitted by an
electron with Lorentz factor $\gamma$, averaged over an isotropic pitch angle
distribution is:

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

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

$P^{\prime syn}_{\nu^{\prime}}(\nu^{\prime})=\int{P^{\prime
syn}_{\nu^{\prime}}(\nu^{\prime},\gamma)\frac{dN}{d\gamma}d\gamma}~{}.$ (42)

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

$P^{\prime
SSC}_{\nu^{\prime}}(\nu^{\prime},\gamma)=\frac{3}{4}h\sigma_{T}c\frac{\nu^{\prime}}{\gamma^{2}}\int\frac{d\tilde{\nu^{\prime}}}{\tilde{\nu^{\prime}}}n_{\tilde{\nu^{\prime}}}K(\gamma,\nu^{\prime},\tilde{\nu^{\prime}})~{},$
(43)

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

####  2.3.5 Pair production

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

$\sigma_{\gamma\gamma}(\beta^{\prime})=\frac{3}{16}\sigma_{T}(1-\beta^{\prime
2})\Bigg{[}(3-\beta^{\prime
4})\ln{\bigg{(}\frac{1+\beta^{\prime}}{1-\beta^{\prime}}\bigg{)}-2\beta^{\prime}(2-\beta^{\prime
2})\Bigg{]}}$ (44)

where:

$\beta^{\prime}(\omega_{t},\omega_{s},\mu)=\Bigg{[}1-\frac{2}{\omega_{t}\omega_{s}(1-\mu)}\Bigg{]}^{\frac{1}{2}}$
(45)

and $\omega_{t}=h\nu^{\prime}_{t}/m_{e}c^{2}$ with $\nu^{\prime}_{t}$ being
the target photon frequency, $\omega_{s}=h\nu^{\prime}/m_{e}c^{2}$ with
$\nu^{\prime}$ being the source photon frequency and $\mu=\cos{\phi}$, where
$\phi$ is the scattering angle. Then, it is possible to derive the
annihilation rate of photons into electron-positron pairs as:

$R(\omega_{t},\omega_{s})=c\int_{-1}^{\mu_{max}}\frac{d\mu}{2}(1-\mu)\sigma_{\gamma\gamma}(\omega_{t},\omega_{s},\mu)~{},$
(46)

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

$R(x)\approx 0.652c\sigma_{T}\frac{x^{2}-1}{x^{3}}\ln{(x)}\text{ H}(x-1)~{},$
(47)

where $H(x-1)$ is the Heaviside function Coppi and Blandford (1990). The
approximation reproduces accurately the behaviour near the peak at
$x_{peak}\sim 3.7$ and over the range $1.3<x<10^{4}$ which usually dominates
during the calculations. A comparison between 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.

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 $\lesssim 7\%$.

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

$\tau_{\gamma\gamma}(\nu^{\prime})=\sigma_{\nu^{\prime}\nu^{\prime}_{t}}n^{\prime}(\nu^{\prime}_{t})\Delta
R^{\prime}$ (48)

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

$\tau_{\gamma\gamma}(\nu^{\prime})=\frac{\Delta R^{\prime}}{c}\int
R(\nu^{\prime},\nu^{\prime}_{t})n^{\prime}_{\nu^{\prime}}(\nu^{\prime}_{t})d\nu^{\prime}_{t}$
(49)

where $\nu^{\prime}$ and $\nu^{\prime}_{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^{-\tau_{\gamma\gamma}})/\tau_{\gamma\gamma}$. 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
$\Delta R^{\prime}=R/\Gamma$, the attenuation of the observed flux due to pair
production become relevant above 0.2 TeV and it reduces the flux by $\sim$
30$\%$ at 1 TeV and by $\sim$ 70$\%$ at 10 TeV.

Figure 4: Spectral energy distribution in the GRB afterglow external forward
shock scenario estimated with the numerical code presented in this section for
$t_{obs}=170$ s. The effect of the pair production attenuation is clearly
evident in the VHE tail. The set of afterglow parameters used are the
following: $E_{k}=1\times 10^{53}$ erg, $s=0$, $A_{0}=1$ cm-3,
$\epsilon_{e}=0.2$, $\epsilon_{B}=0.02$, $\Gamma_{0}=300$, $p=2.5$ and
$z=0.5$.

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

$2\gamma=\omega_{s}+\frac{3}{\omega_{s}}\approx\omega_{s}=\frac{h\nu^{\prime}}{m_{e}c^{2}}$
(50)

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

$Q_{e}^{pp}(\gamma,t^{\prime})=4\frac{m_{e}c^{2}}{h}n_{\nu^{\prime}}(\frac{2\gamma
m_{e}c^{2}}{h},t^{\prime})\int
d\nu^{\prime}_{t}n_{\nu^{\prime}}(\nu^{\prime}_{t},t^{\prime})R(\frac{2\gamma
m_{e}c^{2}}{h},\nu^{\prime}_{t})$ (51)

####  2.3.6 Comparison with analytical approximations

Figure 5: Comparison between spectra estimated with analytical approximations
and numerical calculations. For the analytical method, the synchrotron
emission (green dashed line) is estimated starting from Sari, Re’em et al.
(1998) (SPN98, see legend) while the SSC (black dashed line) is taken from
Sari, Re’em and Esin (2001) (SE01). Vertical lines mark the break frequencies.
The results from the numerical code described in section 2.3 are shown with
solid blue and red lines. This example shows the spectrum calculated at
$t=10^{4}$ s for $s=0$, $p=2.3$, $\epsilon_{e}=0.05$, $\epsilon_{B}=5\times
10^{-4}$, $E_{k}=10^{52}$ erg, $n_{0}=1$ cm-3, $\Gamma_{0}=400$, and $z=1$.

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 Sari, Re’em et al. (1998) and Sari,
Re’em and Esin (2001). In Sari, Re’em et al. (1998) the synchrotron spectra
and light-curves are derived assuming a power-law distribution of electrons in
an expanding relativistic shock, cooling only by synchrotron emission. The
dynamical evolution is described following BM76 equations for an adiabatic
blast-wave expanding in a constant density medium. The resulting emission
spectrum (green dashed lines in Figure 5) is described with a series of sharp
broken power-laws. The SSC component associated to the synchrotron emission
was computed, as a function of the afterglow parameters, in Sari, Re’em and
Esin (2001). 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 Sari, Re’em et al. (1998);
Chevalier and Li (1999). Regarding the flux normalization, there are minor
discrepancies between the numerical and analytical results. This is due to the
fact that in analytical prescriptions it is assumed that the radiation is
entirely emitted at the characteristic synchrotron frequency. On the contrary,
in the numerical derivation, the full synchrotron spectrum of a single
electron is summed up over the whole electron distribution. Similar
considerations apply to the SSC component when comparing with the analytical
spectra. Moreover, the discrepancies observed between analytical and numerical
SSC spectra are amplified by the differences observed in the target
synchrotron spectra.

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

##  3 Open Questions

As predicted by the basic standard model presented in the previous section,
the afterglow emission is the result of particle acceleration and radiative
cooling occurring in two different regions: the forward and the reverse shock.
The temporal and spectral 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
radio-to-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.

###  3.1 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 Sari, Re’em et
al. (1998); Panaitescu and Kumar (2000); Granot and Sari (2002). 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 Zhang et al. (2006); Nousek et al. (2006).

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 $\sim 33\%$ of the
GRBs Chincarini et al. (2010); Margutti et al. (2011). The times at which they
are observed span a very wide range, from around $\sim 30$ s up to $\sim
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 $\omega$ is found to evolve linearly with time to larger values
following the trend $\omega\sim 0.2T_{pk}$ Chincarini et al. (2010). The
average and peak luminosities $L$ of the flare also display a dependence from
$T_{pk}$, with $L\propto T_{pk}^{-2.7}$ at least for early time
($T_{pk}<10^{3}$ s) flares Chincarini et al. (2010); Margutti et al. (2011).
When including also late time flares Bernardini et al. (2011); Yi et al.
(2016) a shallower index is obtain, around $\sim-1.2$. The energy emitted
during flare episodes is quite large and, for early time flares, is around
$\sim 10\%$ of the prompt emission or sometimes even comparable Beniamini and
Kumar (2016).

Flares have been detected also in the optical, although the sample of optical
flares is far smaller than the X-ray one Yi et al. (2017). 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 Yi et al. (2017).

Flares are believed to have an inner origin and to be associated with a
prolonged activity of the GRB central engine Burrows et al. (2005); Willingale
et al. (2010); Fan and Wei (2005); Lazzati and Perna (2007); Maxham and Zhang
(2009); Zhang et al. (2006). 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 Wang et al. (2006); Ajello et al. (2019); Wang and
Dai (2013); He et al. (2012). Assuming that flares have an internal origin and
are produced at $R<R_{dec}$, forward shock electrons will be exposed to the
flare radiation, producing an IC emission component by up-scattering the flare
photons. Following these estimates, the IC component peaks at $\sim 100$ GeV
and has a flux comparable to the X-ray flux. Alternatively, GeV-TeV radiation
associated to flares can be produced by the SSC mechanism, where electrons
responsible for X-ray synchrotron flares also upscatter these photons to
higher energies. The process is considered less interesting for TeV radiation
because the peak of this SSC component is expected to be around 1 GeV Wang et
al. (2006), due to a relatively low minimum Lorentz factor $\gamma_{min}\sim
60$. Such value is estimated from theoretical considerations where
$\gamma_{min}\simeq 60\epsilon_{e,-1}(\Gamma_{IS}-1)$ for $p=2.5$,
$\epsilon_{e}=0.1$ and a relative shock Lorentz factor $\Gamma_{IS}$ of the
order of unity. We notice that the recent estimates of the minimum electron
Lorentz factor in the late prompt emission of GRBs Oganesyan et al. (2019) 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 $\sim 16.5$ years,
from 2005 up to June 2021 Berti, A. for the MAGIC Collaboration . More than
half of them (74 events) have been observed with delays from shorter than
$10^{3}$ s, which means $\sim 4.5$ GRBs yr-1, and 37 events observed with
delays shorter than 100 s (i.e. 2.2 GRB yr-1). Considering that $\sim 33\%$ of
the long GRBs have an X-ray flare and considering the distribution of their
peak times (see Figure 6), we estimate that $\sim 1$ GRBs/yr is the rate of
GRBs with an X-ray flare occurring during MAGIC observations.

Figure 6: Average fluxes of Swift-XRT flares in the 0.3-10 keV energy range
versus peak times $T_{pk}$. The blue points are the 468 flares observed by
Swift from April 2005 to March 2015 (collected from Yi et al. (2016)). The
orange line is the resulting linear regression which gives the expression:
$\log{F_{av}}=-6.76-1.08\times\log{T_{pk}}$.

Let us go a bit further and estimate the detectability of a putative $\sim
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 $\sim 10^{2}$ GeV.

We collect the X-ray flux of a large sample of flares from the catalog of
X-ray flares presented in Yi et al. (2016). Results are shown in Figure 6. The
average flux of the flare and the flare peak time correlate, and the orange
line represents the best fit. To perform the estimates, we consider two
different flare peak times, $T_{pk}=10^{2}$ s and $T_{pk}=10^{3}$ s. The
typical average fluxes at those times are $F=1\times 10^{-9}$ erg cm-2 s-1 and
$F=1\times 10^{-10}$ erg cm-2 s-1, respectively. Assuming that a similar
amount of flux is 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 (Fioretti et al. (2019)) 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 $\omega\sim 0.2T_{pk}$, we can compare the
flare fluxes estimated at $T_{pk}=10^{2}$ s and $T_{pk}=10^{3}$ s with the
differential sensitivity for observing 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 $\sim$ $1.0\times
10^{-9}$ erg cm-2 s-1 and at $t_{obs}=200$ s is $\sim$ $5.0\times 10^{-10}$
erg cm-2 s-1). This indicates that MAGIC telescopes can barely detect such a
flare. Moreover, Extragalactic Background Light (EBL) attenuation reduces the
flux, 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
(CTAhttps://www.cta-observatory.org) 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\sim 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.

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

###  3.2 Density profile of the external medium

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

Long GRBs originate from the core-collapse of massive stars, most likely
rapidly rotating, and with a possible evidence of a preference for low-
metallicity. The most convincing evidence in support of this paradigm is the
association with type Ic supernovae and the proximity of GRBs to young star-
forming regions. While the connection of long GRBs (or at least with the bulk
of the population) with the core-collapse of massive stars is solid, the role
of metallicity and rotation in the launch of the GRB jets, and the
identification of the progenitor star are still uncertain. The progenitor is
usually identified with Wolf-Rayet stars, massive stars ($M>20\,M_{\odot}$) in
the final stages of their evolution, characterised by powerful winds and a
high mass loss rate Woosley (2011). 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\propto
R^{-2}$. The density is related to the mass loss rate $\dot{M}$ and to the
velocity $v_{w}$ of the free-streaming stellar wind by:

$n(R)=\frac{\dot{M}}{4\,\pi\,R^{2}\,m_{p}\,v_{w}}~{}.$ (52)

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
Chevalier and Li (1999).

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

In afterglow modeling of long GRBs it is then customary to assume a density
profile described by eq. 52, where $\dot{M}$ 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_{\star}$: $n(R)=3\times
10^{35}A_{\star}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% Panaitescu and Kumar (2002); Schulze et
al. (2011); Li et al. (2015); Gompertz et al. (2018).

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

Figure 8: Density profile produced by a star wind as a function of the
distance from the central object. The $R^{-2}$ profile characterises the
region where the wind freely streams, which is separated from the shocked wind
(with nearly constant density) by the termination shock. Top panel: the impact
of different values of the ISM density on the termination shock is shown. For
all the curves, it is assumed $\dot{M}=10^{-5}M_{\odot}$ yr-1 and
$v_{w}=10^{3}$ km s-1, while the ISM density varies from 10 to $10^{4}$ cm-3,
with the termination shock moving to smaller distances when the density
increases. Bottom panel: for a fixed ISM density $n_{ISM}=10$ cm-3, the impact
of different $\dot{M}$ and $v_{w}$ is shown. Black-solid:
$\dot{M}=10^{-4}M_{\odot}$ yr-1 and $v_{w}=10^{3}$ km s-1; red-solid:
$\dot{M}=10^{-6}M_{\odot}$ yr-1 and $v_{w}=10^{3}$ km s-1; black-dotted:
$\dot{M}=10^{-5}M_{\odot}$ yr-1 and $v_{w}=500$ km s-1; red-dotted:
$\dot{M}=10^{-5}M_{\odot}$ yr-1 and $v_{w}=2\times 10^{3}$ km s-1. From Fryer
et al. (2006). ©AAS. Reproduced with permission.

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 van Marle et al. (2006).
These considerations on diversity of $E_{k},\dot{M},v_{w}$ and ISM density may
not be sufficient to explain the results of the modeling (i. e., the
preference for a ISM-like environment). The fraction of GRBs which might have
these peculiar parameters can hardly account for the large fraction of GRBs
for which a wind-like profile is excluded by observations. The required
conditions are too extreme to be verified in half of the population. However,
it is not clear if this percentage has been overestimated by present studies.
To quantify the inconsistency, the first step would be to perform a dedicated
study of afterglow emission to assess the percentage of long GRB afterglows
that are not consistent with a wind-like environment.

Methods based on closure relations may not be valid if the spectrum is
modified by Compton scattering in Klein-Nishina regime (see also
Srinivasaragavan et al. (2020)). 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.

###  3.3 Small values of $\epsilon_{B}$

For a long time the typical value of $\epsilon_{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 $\epsilon_{B}\sim 0.01-0.1$. This layer
of intense micro-turbulence is expected on theoretical grounds and recently
corroborated by numerical PIC simulations. Inferences of the value of
$\epsilon_{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 $\epsilon_{B}=10^{-6}-10^{-3}$ Kumar and Barniol Duran (2009, 2010);
Barniol Duran and Kumar (2011); He et al. (2011); Beniamini et al. (2016).
Similar values have been inferred from studies that are based on radio,
optical and/or X-ray emission and do not make use of high energy emission,
such as Barniol Duran (2014); Santana et al. (2014); Wang et al. (2015); Zhang
et al. (2015). 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 $\epsilon_{B}$ may appear to be problematic Piran and
Nakar (2010), because strongly self-generated micro-turbulence must be present
to ensure the scattering and acceleration of the particles, which otherwise
would be simply advected away.

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

The application of cooling in a decaying magnetic turbulence to four GRBs
detected by LAT has proved to be very successful and even able to give
indications on how fast the turbulence decays, being consistent with a power-
law decay $\epsilon_{B}\propto t^{-\alpha_{t}}$ with $\alpha_{t}\sim 0.5$
Lemoine et al. (2013).

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 Santana et al. (2014); Gompertz
et al. (2018).

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

###  3.4 Variation of micro-physical parameters with time

Thanks to the increasing number of available observations on a wide range of
frequencies (from radio to TeV) and times (from seconds to weeks), the basic
assumption that micro-physical parameters (such as
$\epsilon_{e},\epsilon_{e},p$ and $\xi_{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 $\epsilon_{e}$ and $\epsilon_{B}$
may vary with time have been found in some events with well detailed multi-
wavelength follow-up campaigns.

In Filgas et al. (2011) 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 $\nu_{cool}$ evolves as $\nu_{cool}\propto t^{-1.2}$
which is in contrast with synchrotron predictions for which a less steeper
decay $\nu_{cool}\propto t^{-0.5}$ is expected. As a result, some assumptions
of the standard model must be relaxed to remove the tension between
observations and theoretical predictions. The most likely option able to
explain the cooling break observational behaviour without affecting the
general interpretation of the data is to let the $\epsilon_{B}$ parameter
variate with time. Assuming that $\epsilon_{B}\propto t^{0.49}$ the time
evolution of $\nu_{cool}$ can be explained successfully.

In van der Horst et al. (2014) for GRB 130427A modeling, in order to explain
the observed fast evolution of the injection frequency $\nu_{m}\propto
t^{-1.9}$ a temporal evolution of $\epsilon_{e}$ is claimed. Considering that
$\nu_{m}\propto\epsilon_{e}^{2}$, a modest evolution of $\epsilon_{e}$
following the trend $\epsilon_{e}\propto t^{-0.2}$ is able to satisfactorily
describe the observed light-curves.

A time-dependent evolution of the micro-physical parameters has also been
proposed in order to explain the features observed in the early afterglow
phase which are not predicted by the external forward shock scenario such as
X-ray afterglow plateaus, chromatic breaks, and afterglow rebrightenings Ioka
et al. (2006); Kong et al. (2010); Huang et al. (2018); Panaitescu (2007).

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 Misra et al. (2021), as
will be discussed in the next section.

###  3.5 Maximum synchrotron photon energy

One of the expectations from Fermi-LAT observations of GRB afterglows was the
identification of a spectral cutoff in the afterglow synchrotron spectrum
marking the maximum energy of synchrotron photons Guilbert et al. (1983); de
Jager and Harding (1992). 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:

$t^{\prime}_{L}=\frac{r_{L}}{c}=\frac{\gamma m_{e}c}{eB^{\prime}}$ (53)

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

$\delta E^{\prime}=t^{\prime}_{L}P^{\prime}=\frac{\gamma
m_{e}c}{eB^{\prime}}\frac{\sigma_{T}c\gamma^{2}B^{\prime
2}}{6\pi}=\frac{1}{6\pi}\frac{\sigma_{T}m_{e}c^{2}B^{\prime}\gamma^{3}}{e}$
(54)

The particle stops to gain energy when:

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

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

$\gamma_{max}=\sqrt{\frac{3\pi e}{\sigma_{T}B^{\prime}}}$ (56)

The corresponding maximum synchrotron photon energy is:

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

which for electrons is $\sim$ 50 MeV in their rest frame.

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

$\gamma_{cool,p}=\frac{6\pi m^{3}_{p}c}{\sigma_{T}m^{2}_{e}B^{\prime
2}t^{\prime}}$ (58)

for the cooling Lorentz factor and:

$\gamma_{max,p}=\sqrt{\frac{3\pi em^{2}_{p}}{\sigma_{T}m^{2}_{e}B^{\prime}}}$
(59)

for the maximum Lorentz factor which sets a maximum photon energy of $\sim$
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 $\sim$ -2 and with no evidence of spectral evolution in time Ajello
(2019) 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^{\prime}$ 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 Kumar et al.
(2012) the magnetic field $B^{\prime}$ can be expressed in terms of the
distance from the shock front $x$ as:

$B^{\prime}(x)=B^{\prime}_{s}\bigg{(}\frac{x}{L_{p}}\bigg{)}^{-\eta}+B^{\prime}_{w}~{},$
(60)

where $B^{\prime}_{s}$ and $B^{\prime}_{w}$ are respectively the strongest and
the weakest magnetic field strengths, $\eta$ is the power-law decaying index,
and $L_{p}$ is the field decay length scale, which is estimated as Medvedev
and Loeb (1999):

$L_{p}=\Bigg{[}\frac{m_{p}\Gamma_{s}c^{2}}{4\pi
ne^{2}}\Bigg{]}^{1/2}=2.2\times 10^{7}\frac{\Gamma_{s}}{n}\hskip
8.5359pt\text{cm}$ (61)

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

$\bigg{(}\frac{B^{\prime}_{s}}{B^{\prime}_{w}}\bigg{)}^{2}\lesssim\frac{r_{l}}{L_{p}}$
(62)

valid for $\eta>1/2$ and $x_{0}/L_{p}\gg 1$ where $x_{0}$ is the width of the
high magnetic field region. Considering that
$x_{0}/L_{p}\equiv(B^{\prime}_{s}/B^{\prime}_{w})^{1/\eta}$, 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:

$\gamma_{max}=\sqrt{\frac{3\pi e}{\sigma_{T}B^{\prime}_{w}}}$ (63)

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

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

which is greater than the one calculated in eq. 57 by a factor
$B^{\prime}_{s}/B^{\prime}_{w}$. Numerical calculations Sironi and Spitkovsky
(2011) show that this ratio can be larger than $\sim 10^{2}$. As a result,
photons of energies $\gtrsim 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.

###  3.6 Prompt emission efficiency

The overall efficiency $\eta_{\gamma}$ of the prompt emission mechanism is the
result of three processes: the efficiency $\eta_{diss}$ of the (still
unidentified) mechanism responsible for dissipation of the jet energy, the
efficiency $\epsilon_{e}$ of the acceleration mechanism in converting the
dissipated energy into random energy of the electrons, and the radiative
efficiency $\epsilon_{rad}$ of the electrons:
$\eta_{\gamma}=\eta_{diss}\epsilon_{e}\epsilon_{rad}$. Provided that it is
reasonable to assume a fast cooling regime for the prompt emission
($\epsilon_{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 $\eta_{\gamma}=E_{\gamma}/E_{0}$ (where
$E_{0}=E_{\gamma}+E_{k}$ is the initial explosion energy), we can write the
relation $E_{k}=(1-\eta_{\gamma})/\eta_{\gamma}E_{\gamma}$. The parameter
$\eta_{\gamma}$ can then be estimated from the comparison of the energy
$E_{\gamma}$ emitted in the prompt phase and the energy $E_{k}$ left in the
jet after the end of the prompt emission (i.e., at the beginning of the
afterglow phase). While the former is directly estimated from observations,
the latter 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 Kumar (2000); Berger et al. (2003); Lloyd-Ronning and Zhang
(2004); Berger (2007); D’Avanzo et al. (2012). This method is solid as long as
the X-ray band lies above $max(\nu_{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}\epsilon_{e}$.
Assuming a value (typically 0.1) for $\epsilon_{e}$, then $E_{k}$ can be
estimated. Investigations based on the X-ray emission have inferred very large
values of $\eta_{\gamma}$, between 0.5 and 0.9 Lloyd-Ronning and Zhang (2004);
Granot et al. (2006); Nousek et al. (2006); Zhang et al. (2007).

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_{\gamma,iso}$ has been indeed
found, supporting the possibility that GeV photons lie in the high-energy part
of the synchrotron spectrum, where the afterglow luminosity is proportional to
$E_{k}\epsilon_{e}$ and can be use, similarly to the X-ray luminosity, to
estimate $E_{k}$ Nava et al. (2014). A study by Beniamini et al. (2015)
revealed that the energetics $E_{k}$ inferred independently from X-ray and GeV
luminosities on a sample of 10 GRBs are inconsistent with each other. The
authors show that the inconsistency is solved if $\nu_{c}>\nu_{X}$ (where
$\nu_{X}$ is the X-ray frequency), or if Compton losses are important in the
X-ray band. Full modeling of the GeV, X-ray and optical data support this
scenario. In both cases, $\epsilon_{B}$ is required to be much smaller than
usually assumed, with values in the range $10^{-7}-10^{-3}$. This analysis
shows that the GeV band is a much better proxy for $E_{k}$, since it is above
$\nu_{c}$ and is not affected by inverse Compton cooling, due to Klein-Nishina
suppression. Adopting GeV luminosities as a proxy for $E_{k}$, the estimated
values of $E_{k}$ increase, affecting also the estimates of $\eta_{\gamma}$,
which are around $5-10$% Beniamini et al. (2016).

A correct estimate of $\eta_{\gamma}$ is extremely important, since its value
is related to the mechanism dissipating energy in the jet. Since internal
shocks can hardly reach values of $\eta_{\gamma}$ larger than 10%, values
around 50-90% have been used to argue that internal shocks are not a viable
mechanism to explain prompt emission in GRBs, and more efficient mechanisms
should be considered (e.g., magnetic reconnection). If the efficiency is
however smaller than initially estimated, internal shocks may still be a
viable solution. Moreover different estimates of $\eta_{\gamma}$ lead to
different estimates on the total initial jet energy
$E_{0}=E_{\gamma,iso}+E_{k}$, with repercussions on the energy budget of GRBs
and finally on their progenitors and mechanisms for jet launching. Small
values of $\epsilon_{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., Zhang and Yan (2011); Lazarian et al. (2019); Beniamini and Granot
(2016)).

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.

###  3.7 Fraction of accelerated particles

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

Figure 9: Distribution of accelerated electrons from a relativistic shock. The
two bumps correspond respectively to the distribution of the $(1-\xi)$
fraction of particles heated into a thermal component with post shock energy
$\sim\gamma m_{e}c^{2}$ and the fraction $\xi$ of particles accelerated into a
non-thermal component with post shock energy $\sim\gamma m_{p}c^{2}$. Adapted
from Eichler and Waxman (2005). ©AAS. Reproduced with permission.

In particular, the parameter $\xi_{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-\xi_{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 $\xi_{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
$\xi_{e}=1$ for parameters $E_{k}$, $n$, $\epsilon_{e}$ and $\epsilon_{B}$
cannot be distinguished from those one estimated for any $\xi_{e}$ in the
range $m_{e}/m_{p}<\xi_{e}<1$ and afterglow parameters
$E^{\prime}_{k}=E_{k}/\xi_{e}$, $n^{\prime}=n/\xi_{e}$,
$\epsilon^{\prime}_{e}=\xi_{e}\epsilon_{e}$ and
$\epsilon^{\prime}_{B}=\xi_{e}\epsilon_{B}$ Eichler and Waxman (2005). This
can be proven considering the dependencies of the jet dynamics and shock
energy equipartition on the model parameters. As shown in BM76 and by previous
calculations on the evolution of a relativistic blastwave, in the self-similar
regime the bulk Lorentz factor evolve as $\Gamma\propto(E_{k}/n)^{1/2}$. As a
result, the same flow evolution can be obtained with different values of
$E_{k}$ and $n$ as long as their ratio is preserved. The fraction of energy
given to the magnetic field is reduced by a factor $\xi_{e}$ but at the same
time the energy density given to the shock is increased by a factor
$1/\xi_{e}$ so that the magnetic field energy density is the same in the
scenario when including or not $\xi_{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 $\xi_{e}$ or for any
value $m_{e}/m_{p}<\xi_{e}<1$. This imply that afterglow model parameters,
when considering $\xi_{e}\neq 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 $\epsilon_{B}$) but rather a fraction of their value
multiplied or divided by $\xi_{e}$ (e.g. $E_{k}/\xi_{e}$ or
$\epsilon_{B}\xi_{e}$).

It is possible to obtain information on the value of $\xi_{e}$ through PIC
simulations or indirect features of the thermal component on the synchrotron
afterglow spectra. As mentioned in the previous Section, PIC simulations of
weakly magnetized shocks have found that the downstream population include
around $\sim 3\%$ of the electrons which are accelerated into a non-thermal
distribution. In case the efficiency is low (around $\sim 10\%$ or less) the
presence of a large population of thermal electrons may affect the afterglow
radiation spectra. The thermal electrons are distributed at lower energies
than the non-thermal ones since $\eta\gamma m_{e}c^{2}$ $\ll$ $\gamma
m_{p}c^{2}$ where $\eta\ll m_{p}/m_{e}$ is a factor describing the ignorance
on the plasma physics governing electron heating beyond $\gamma m_{e}c^{2}$.
As a result, the synchrotron radio component emission may be affected through
the production of a new emission component from thermal electrons (for
$\eta\gg 1$ and moderate $1/\xi_{e}$) or a large self-absorption optical depth
(for $\xi_{e}\ll 1$) which may be visible in a time scale of few hours.
Possible effects of the thermal component are discussed in Ressler and Laskar
(2017); Warren et al. (2022); Giannios and Spitkovsky (2009).

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
$\epsilon_{B}$ and the density. This will impact also the $\xi_{e}$
calculation since it can help to reduce the degeneracy between the sets of
solutions available in the parameter space. Indeed, the multi-wavelength
modeling of GRB 190829A (detected at TeV energies by H.E.S.S.) showed that the
only way to explain all the detected radiation, from radio to TeV, is to
introduce the parameter $\xi_{e}$ in the modeling, which is constrained by
data to be $\xi_{e}<0.13$ Salafia et al. (2021).

##  4 Discovery of a TeV emission component in GRBs

The robust theoretical framework developed throughout the years to explain the
afterglow radiation predicts that, to some extent, GRBs should be TeV emitters
(section 2). Observations in the HE band, and in particular the presence of
$\gtrsim$ GeV photons with energies up to $\sim 100$ GeV, support this
possibility. On the other hand, from the observational side the search for
such emission is hampered by several drawbacks. Space-born telescopes, such as
Fermi-LAT, sensitive up to few hundred GeV, have an hard time with GRBs due to
their low $\gamma$-ray photon flux at the highest energies ($\sim 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
($\gtrsim 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 Mirzoyan et al. (2019). 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 $\gtrsim 3$ TeV. The timescales of
the detections place the origin of the emission in the afterglow phase. The
TeV emission has been studied and interpreted in a multi-wavelength context,
in order to evaluate the properties and the nature of the responsible
radiation mechanisms. In particular, investigations have focused on SSC,
external inverse Compton (EIC), and synchrotron radiation.

In this section, all GRBs for which a detection (significance $>5\sigma$) or a
hint of detection (significance between 3 and 5$\sigma$) has been claimed by
Cherenkov telescopes are presented. These are in total six events (one short
and five long): GRB 160821B (Section 4.1), GRB 180720B (Section 4.2), GRB
190114C (Section 4.3), GRB 190829A (Section 4.4), GRB 201015A (Section 4.5)
and GRB 201216C (Section 4.6). For each event we start with a brief
description of the prompt and afterglow multi-wavelength observations. Then,
we describe VHE observations and summarise the main results. Being these