Gamma-Ray Bursts at TeV Energies: Theoretical Considerations

Gamma-ray bursts (GRBs) are the most luminous explosions in the Universe and are powered by ultra-relativistic jets. Their prompt $\gamma$-ray emission briefly outshines the rest of the $\gamma$-ray sky making them detectable from cosmological distances. It is followed by, and sometimes partially overlaps with, a similarly energetic but very broadband and longer-lasting afterglow emission. While most GRBs are detected below a few MeV, over a hundred were detected at high ($\gtrsim0.1\;$GeV) energies and several have now been observed up to tens of GeV with the \textit{Fermi} Large Area Telescope (LAT). A new electromagnetic window in the very high energy (VHE) domain ($\gtrsim0.1\;$TeV) was recently opened with the detection of afterglow emission in the $(0.1$\textendash$1)\,$TeV energy band by ground-based imaging atmospheric Cherenkov telescopes. The emission mechanism for the VHE spectral component is not fully understood, and its detection offers important constraints for GRB physics. This review provides a brief overview of the different leptonic and hadronic mechanisms capable of producing VHE emission in GRBs. The same mechanisms possibly give rise to the high-energy spectral component seen during the prompt emission of many \textit{Fermi}-LAT GRBs. Possible origins of its delayed onset and long duration, well into the afterglow phase, with implications for the emission region and relativistic collisionless shock physics are discussed. Key results for using GRBs as ideal probes for constraining models of extra-galactic background light and intergalactic magnetic fields, as well as for testing Lorentz invariance violation, are presented.


Introduction
Gamma-ray bursts (GRBs) are cataclysmic events that occur at cosmological distances. (See, e.g., [1][2][3][4] for a comprehensive review.) They are the most electromagnetically luminous transient phenomena in the Universe. GRBs involve the explosive release of energy over a short timescale, producing a burst of γ-rays with isotropic-equivalent luminosity of L γ,iso ∼10 51 -10 54 erg s −1 . Their emission is powered by ultrarelativistic (with bulk Lorentz factors Γ 100) bipolar collimated outflows driven by a compact object central engine. The identity of the central engine, which could be either a black hole (BH) or a millisecond magnetar, is not entirely clear as the highly variable emission is produced far away from it at a radial distance of R ∼ 10 12 -10 16 cm. The most luminous phase of the burst, referred to as the "prompt" phase, is short-lived with a bimodal duration distribution separated at t ∼ 2 s, where the short and long GRBs have typical observed durations of t GRB ∼ 10 −0.5 s and t GRB ∼ 10 1.5 s, respectively [5]. These two classes of GRBs are also distinct spectrally, with the short GRBs being spectrally harder compared to the long GRBs that produce softer γ-rays. The long-soft GRBs are associated with the core collapse of massive ( (20-30)M ) Wolf-Rayet stars [6,7], whereas (at least some) short-hard GRBs originate in compact object mergers of two neutron stars (NSs) or a NS-BH binary [8,9]. The first-ever detection of a short GRB coincident with gravitational waves (GWs) from the merger of two NSs came from GW 170817/GRB 170817A [10,11].
Many details of the prompt GRB emission, in particular, the energy dissipation process, the exact radiation mechanism, and the transfer of radiation in the highly dynamical flow remain poorly understood. All of these different processes combine to produce a non-thermal spectrum that is often well-described by the Band-function [12], an empirical fit to the spectrum featuring a smoothly broken power law. This break manifests as a peak in the νF ν spectrum, at a mean photon energy E br 250 keV, with the asymptotic power-law photon indices below and above the peak energy having mean values of α Band −1 and β Band −2.3, respectively [13,14].
While most of the energy in the prompt GRB comes out at E E br , the featureless power-law spectrum above this energy extends beyond 100 MeV in most GRBs detected by the Fermi-Large Area Telescope (LAT) [15], with a high-energy spectral cutoff seen in only about 20% of the cases, e.g., [16], which are most likely caused by intrinsic opacity to pair production [16,17]. In rare cases, the prompt GRB spectrum shows an additional hard spectral component that extends beyond ∼10 GeV, as seen by the Fermi-LAT, and well into very high energies ( 0.1 TeV), as seen by ground-based atmospheric Cherenkov telescopes, e.g., MAGIC and H.E.S.S. (See, e.g., [18,19] for a review.) This high-energy (HE; 100 MeV) emission overlaps with the sub-MeV prompt GRB, and both the HE and very-high-energy (VHE) emissions persist throughout the afterglow phase-the much longer-lasting and broadband (X-rays/Optical/Radio) emission that follows the short-lived prompt phase. The spectral and temporal properties of the HE emission provide a glimpse into the global energetics of the bursts as well as yield important constraints on GRB physics that cannot be obtained from the sub-MeV emission alone.
The main objective of this review is to provide a concise summary of the widely discussed radiation mechanisms that may explain the spectral and temporal properties of the VHE and/or HE emission in GRBs. We first discuss several HE/VHE radiation mechanisms in Section 2 and provide some of the fundamental quantities that can be calculated and compared to observations. This is followed by a discussion of the delayed HE emission, additional prompt GRB spectral component at high energies, and long-lived HE emission seen by Fermi-LAT as well as popular theoretical explanations offered for it, along with implications for the bulk Lorentz factor Γ, in Section 3. Next, Section 4 presents an overview of the HE afterglow seen in the exceptionally bright GRB 130427A, along with several important implications for the radiation mechanism and relativistic shock acceleration physics. The recent detection of a ∼TeV afterglow emission by MAGIC and H.E.S.S in only a few GRBs and key implications of such a detection for GRB physics are discussed in Section 5. The use of HE photons from distant GRBs as a probe of extra-galactic background light (EBL), inter-galactic magnetic field, and Lorentz invariance violation are the topics of discussion in Section 6. Finally, in Section 7, we end this review with important outstanding questions in GRB physics and present closing remarks in Section 8.

Relevant High-Energy or Very-High-Energy Emission Mechanisms
There are several HE/VHE γ-ray emission mechanisms that operate wherever particles (leptons and hadrons) are accelerated to or generated with high Lorentz factors (LFs). In GRBs, the emission regions can be either internal to the relativistic outflow, e.g., at internal shocks or magnetic reconnection sites, or external to it, e.g., in the shocked external medium behind the external forward (afterglow) shock, or even at larger distances from the outflow. Below, we review some of the widely discussed processes that are capable of producing HE to VHE γ-ray photons. Other more detailed reviews on this topic are [4,20].

Electron Synchrotron Emission
Relativistic electrons with LFs γ e 1 cool by emitting synchrotron photons when gyrating around magnetic field lines with comoving magnetic field strength B (all primed quantities are in the comoving/fluid rest frame). At collisionless shocks (internal or external), a fraction ξ e of the electrons are accelerated into a non-thermal power-law energy distribution, dN/dγ e ∝ γ −p e for γ m ≤ γ e ≤ γ M and 2 p 3, that holds a fraction e of the post-shock internal energy density, and arises due to Fermi acceleration [21][22][23]. The minimal LF of this distribution is where Γ ud is the relative LF between the regions upstream and downstream of the shock front. The resulting (observed) optically thin synchrotron spectrum in this case comprises multiple power-law segments joined smoothly at characteristic break energies [24,25] (shown here for fiducial parameters relevant for prompt emission for which Γ 1 and β 1), whereh = h/2π with h being the Planck's constant, σ T is the Thomson cross section, e is the elementary charge, m e is the electron rest mass, and c is the speed of light. The energy E m corresponds to the characteristic synchrotron frequency (ν m ) of minimal energy electrons with LF γ m , and the cooling break energy E c corresponds to the cooling frequency (ν c ) of electrons with LF γ c = (6πm e c 2 /σ T )(Γβ/B 2 R) ≈ 2.2R 13 Γ 3 2.5 f −1 σ,−2 L −1 iso,52 that are cooling at the dynamical time, such that their synchrotron cooling time, t syn = 6πm e c/σ T B 2 γ e , equals the dynamical time, t cool = t dyn = R/Γβc. For some model parameters, γ c < 1, which is obviously unphysical, but instead represents very rapid cooling of particles to non-relativistic velocities in less than the dynamical time [26]. As a result, relativistically hot particles only occupy a thin layer behind the shock which is a fraction γ c of the comoving width ∆ of the ejecta shell, where the electrons are cold in the remaining majority of the shell. In the above equations, we have expressed the comoving magnetic field in terms of the more useful quantities, using the fact that the total isotropic-equivalent power of the outflow can be written in terms of L k,iso and L B,iso , the kinetic energy and magnetic field powers, respectively. As a result, is the fraction of total power carried by the magnetic field with σ = L B,iso /L k,iso being the outflow magnetization, and L γ,iso is the isotropic-equivalent γ-ray luminosity which is a fraction γ of the total power. This yields the comoving B-field strength B ≈ 1.8 × 10 4 f 1/2 σ,−2 L 1/2 iso,52 R −1 13 Γ −1 2.5 G with β 1 for an ultra-relativistic flow. The ordering of the break energies depends on whether the electrons are in the fast cooling regime, for which E c < E m , or the slow cooling regime, with E m < E c . This relative ordering also decides the values of the spectral indices of the flux density F E for the different power law segments, The emission in the power-law segment above the spectral peak energy (max(E c , E m )) can only extend up to the maximum synchrotron energy E syn,max . This energy depends on the efficiency of the acceleration process while the charged particles (electrons or protons) lose energy to synchrotron cooling. The typical timescale t acc over which particles, say the electrons with LF γ e , are accelerated as they are scattered across the relativistic shock is at best the Larmor time t L = γ e m e c/eB , i.e., t L /t acc = κ acc ≤ 1. Their radiative cooling timescale t c is at most t syn as any additional radiative cooling besides synchrotron (e.g., inverse-Compton) would only shorten t c , i.e., t c = κ c t syn with κ c ≤ 1. Equating the acceleration and radiative cooling timescales, t acc = t c , yields the maximum LF attained by the electrons, γ M = (6πeκ/σ T B ) 1/2 where κ = κ acc κ c ≤ 1. These electrons then radiate at the characteristic synchrotron energy, e.g., [27][28][29][30][31] where α F = e 2 /hc 1/137 is the fine structure constant, and κ is a factor expected to be of order unity that depends mainly on the details of particle acceleration and diffusion in the shock downstream and upstream. It is therefore challenging to explain VHE photons as arising from synchrotron emission by electrons. In addition, depending on the compactness of the emission region, emission can be suppressed due to e ± -pair production via γγ-annihilation (γγ → e − e + ), e.g., [32][33][34][35][36][37]. This poses more of a problem for the prompt emission and less so for the afterglow. Alternatively, the VHE photons can be explained by proton synchrotron emission (see Section 2.2) or synchrotron self-compton (SSC; see Section 2.3) emission by the same electron population that produced the seed synchrotron radiation.

Proton Synchrotron Emission
High-energy protons that are accelerated at shocks (like the electrons) to LFs γ p can also cool by emitting synchrotron photons in magnetized regions [38,39]. However, the emitted power per particle (P syn ∝ σ T γ 2 is much smaller where it is suppressed by a factor (m e /m p ) 2 (1836) −2 3 × 10 −7 with respect to that for electrons when γ p = γ e (and suppressed by the square of this factor for E p = E e m p c 2 ) since the Thomson scattering cross section for protons is much smaller, σ T,p = (m e /m p ) 2 σ T , than that of electrons. To compensate for this suppression, the magnetic field in the emission region must be larger than that obtained in a leptonic synchrotron scenario, so much so that the magnetic field energy would hold a good fraction of the total energy [40]. The characteristic synchrotron energy of minimal energy protons is E m,p = (γ m,p /γ m,e ) 2 (m e /m p )E m ≈ [ξ e p /ξ p e ] 2 (m e /m p ) 3 E m (assuming the electrons and protons hold fractions e and p = 1 − e − B of the post-shock internal energy, and that fractions ξ e and ξ p of the electrons and protons, respectively, form a power-law energy distribution), and the cooling break energy is E c,p = (m p /m e ) 5 E c , with the corresponding LF γ c,p = (m p /m e ) 3 γ c . As a result, the maximum LF of protons accelerated at the same shock as electrons is γ M,p = (m p /m e )γ M , which yields E p,syn,max 13(1 + z) −1 κ p Γ 2 TeV, e.g., [39].
Recent suggestions replacing electron with proton synchrotron emisson have been made to explain the apparent low-energy (below the spectral peak) spectral breaks that are difficult to explain with electron synchrotron emission, e.g., [41]. However, knowing that protons are inefficient at radiating away their internal (or random-motion) energy as compared to electrons, the significant reduction in radiative efficiency must be compensated by having a much larger total energy budget, a requirement that may be too demanding, e.g., [42,43]. Moreover, in such a scenario it would also be very difficult to suppress the much more efficient radiation from the electrons for it to not over-power that from the protons.

Synchrotron Self-Compton (SSC)
A distribution of relativistic electrons can inverse-Compton scatter some of the same synchrotron photons that it produced, leading to a synchrotron self-Compton emission. When the energy of the incoming synchrotron photon in the rest frame of the scattering electron is much smaller than the electron's rest energy, E syn ∼ γ e E syn m e c 2 , then the scattering occurs in the Thomson regime (where the electron's recoil can be neglected) and is called elastic or coherent. The scattered photon emerges with an energy of E SSC ∼ γ e E SSC ∼ = γ e E syn ∼ γ 2 e E syn . The additional cooling of particles due to inverse-Compton scattering introduces a factor of (1 + Y) in the cooling time, such that t c = t syn /(1 + Y). Here Y(γ e ) ≡ P IC (γ e )/P syn (γ e ) is the Compton-y parameter given by the ratio of the power radiated in the IC component to that in the synchrotron component.
In the Thomson regime, where U γ is the energy density of the seed synchrotron emission that is IC scattered by the electrons. If this seed radiation arises from shock-heated electrons, U γ = ηβU e /(1 + Y), e.g., [44], where η = min[1, (ν m /ν c ) (p−2)/2 ] is the fraction of electron energy radiated away in synchrotron and IC photons, and β is the downstream velocity relative to the shock front (and is order unity for a relativistic shock). With U e = e U int and U B = B U int , where e and B are the fractions of the total internal energy behind the shock (U int ) that goes into accelerating electrons and generating the magnetic fields, the expression for Y simplifies to, e.g., [26,[44][45][46][47], When η e B , then Y 1, and Compton cooling is negligible. Otherwise, the extra cooling also means that the maximum particle LF is reduced,γ M = (1 + Y) −1/2 γ M , and likewisẽ E syn,max = (1 + Y) −1 E syn,max . The characteristic spectral break energies of the SSC spectrum corresponding to that of the synchrotron spectrum are The maximum energy of an inverse-Compton scattered photon is E IC,max = γ e m e c 2 , and since γ e ≤γ M for a power-law electron distribution, TeV .
When the energy of the incoming photon in the rest frame of the scattering electron exceeds the rest mass energy of the electron, E syn ∼ γ e E syn > m e c 2 , the recoil suffered by the electron can no longer be ignored, and quantum corrections need to be taken into account. The scattering no longer occurs in the Thomson regime, and the correct scattering cross section in this case is the Klien-Nishina cross section (σ KN ), which depends on the energy of the incoming photon [48]. For incoming photon energy x = hν /m e c 2 1, the scattering cross section is highly suppressed, with σ KN (x) ∝ x −1 . Moreover, the electron recoil implies that E SSC (γ e ) ∼ γ e m e c 2 = E KN (γ e ) in this limit. Therefore, IC scattering can efficiently cool an electron with LF γ e only for seed synchrotron photons with energies E syn < m e c 2 /γ e . Thus, accounting of Klein-Nishina effects causes the Compton-Y parameter of each electron to depend on its LF, Y = Y(γ e ) ≡ P IC /P syn ≈ U γ [E syn < m e c 2 /γ e ]/U B . This may cause interesting modifications of the spectrum (in both the synchrotron and SSC components) when η e B [47]. Notice that since Y = Y(γ e ) may vary between different electrons, it is natural to define the global Compton-Y parameter byȲ = L IC /L syn , which is the mean value of Y(γ e ) weighted by the synchrotron emissivity. Therefore, the SSC flux is suppressed above the photon energy where E = E m (E = E c ) are the energies where the synchrotron νF ν spectrum peaks in the fast (slow) cooling scenario. Likewise, the spectral peak of the SSC spectrum occurs at in the fast (slow) cooling case. The ratio of the spectral peak flux is given by , second order SSC scatterings also occur in the Thomson regime, and ifȲ is not much smaller than unity, a third spectral peak can appear, e.g., [49,50].

External Inverse-Compton
External inverse-Compton emission (EIC) arises when the softer seed photons are inverse-Compton scattered to high energies by relativistic electrons in a location physically distinct from where the seed photons were produced. This can occur in several different ways, e.g., (i) seed photons produced in internal dissipation and upscattered by forward-shock or reverseshock-heated electrons [51][52][53][54][55][56][57][58], (ii) seed photons produced in the reverse shock and upscattered by forward shock-heated-electrons [59,60], (iii) seed photons produced in the forward shock and upscattered by reverse shock-heated electrons [59,60], (iv) externally produced ambient seed photons, e.g., from the accretion disk [61] or the massive star progenitor's envelope [62], are upscattered by cold electrons in the relativistic outflow in a process also referred to as bulk Compton scattering or Compton drag, (v) photospheric seed photons in the relativistic baryon-poor jet upscattered by the shocked electrons in the shock transition layer between the baryon-poor jet and baryon-loaded envelope [63], and (vi) seed photons provided by the cocoon [64], after it breaks out of the dynamical ejecta in a NS-NS merger, or that from the AGN disk [65], if the merger occurs inside the disk of an AGN, that are IC upscattered to VHE γ-rays by electrons energized in the dissipation of prolonged jets powered by late-time central engine activity. As an illustrative example, below we summarize the important points for the simplest case in scenario (i) and provide estimates of the maximum photon energy obtained in this process when the X-ray flare emission overlaps with the external forward shock electrons [53].

IC Scattering of X-ray Flare Photons by External Forward Shock Electrons
As the relativistic ejecta plows through the circumburst medium (CBM), with density ρ = AR −k where R is the radial distance from the central engine, it is slowed down. In the process, two shocks are formed where the shocked regions are separated by a contact discontinuity that has a bulk LF Γ. The forward shock runs ahead of the contact discontinuity with bulk LF Γ fs = √ 2Γ, sweeping up the CBM and shock-heating it. The reverse shock moves backward (in the rest frame of the contact discontinuity) into the ejecta, decelerating and shock-heating it. In the following, we adopt the thin-shell case for which the reverse shock is Newtonian (or mildly relativistic). Alternatively, the reverse shock becomes relativistic before crossing the ejecta shell in the thick-shell case, which we will not discuss here (but see Sari and Piran [66]). Most of the isotropic-equivalent kinetic energy of the ejecta (E k,iso ) is transferred to the kinetic and internal energy of the shock-heated swept up CBM behind the forward shock at the deceleration radius, (11) where A = m p n = 1.67 × 10 −24 n 0 g cm −3 for k = 0 (ISM) and A =Ṁ/4πv w = 5 × 10 11 A g cm −1 for k = 2 (wind medium; A = 1 corresponds to a mass loss rate ofṀ = 10 −5 M yr −1 with a wind speed of v w = 10 8 cm s −1 ), and Γ 0 1 is the initial LF of the relativistic ejecta at which it coasts for R < R dec . For R > R dec , the blast wave dynamics become self-similar, and the bulk LF of the shocked material decays as a power law in R [67]. The transition for Γ(R) from the coasting to the self-similar power-law phase is smooth in general, but here we use the broken power-law approximation in Equation (11) for simplicity.
The LF of minimal energy power-law electrons accelerated at collisionless shocks is given by γ m in Equation (1). For electrons accelerated at the forward shock Γ ud = Γ(R) 1, in which case the minimal particle LF for R > R dec is given by (for p = 2.5) at the apparent time t = (1 + z)R/2(4 − k)cΓ 2 = 10 3 t 3 s. Here the factor ζ ≡ Γ 2 ct/R(1 + z) represents a one-zone approximation and is taken here to be 1/2(4 − k), which is appropriate along the LoS (corresponding to the radial time t r ) if Γ is taken to be that of the shock front, Γ sh . If instead, it is taken to be that of the matter just behind the shock, Γ ≈ Γ sh / √ 2 then ζ = 1/4(4 − k) along the LoS. Since there is significant contribution to the observed flux up to angles θ 1/Γ from the LoS, one should also account for the angular time t θ = R/2cΓ 2 (R) along the equal arrival time surface from the shock front. Finally, the exact value of ζ also depends on the effective thickness of the radiating shell [24,68,69] and any value is only as good as the one-zone approximation it represents.
If the spectral peak (of νF ν ) energy of the X-ray flare is E x in the observer frame, its energy in the comoving frame of the blast wave is for X-ray flare photons that are tightly beamed in the radial direction and catch up with the electrons behind the shock with (almost) radial velocity vectors. When the forwardshock electrons are in the fast cooling regime, the peak of the IC spectral component corresponds (without accounting for Klein-Nishina effects) to upscattering of ∼E x seed photons (flare photons considered monoenergetic here for simplicity) by ∼ γ f ,m electrons, e.g., [53], The spectrum of this GeV flash is expected to have power-law spectral indices d ln F ν /d ln ν of approximately −1/2 and −p/2 below and above the energy E IC,pk . Klein-Nishina effects start to become important for electrons with LF γ e ≥ γ e,KN = m e c 2 /E x = 3.3 × 10 4 (1 ,KN E x , and depending on the ratio of electron LFs, the Klein-Nishina suppression of the IC scattered spectrum can occur at energies below or above E IC,pk for ψ < 1 or ψ > 1, respectively.
VHE primary γ-ray photons with energy E from GRBs dominantly annihilate with the much softer EBL photons having energy over the mean free path length of The produced e − and e + will share the energy of the primary γ-ray photon equally and have a typical LF γ e = E(1 + z)/2m e c 2 ≈ 10 6 (1 + z)E 1 TeV . These pairs will then IC scatter the more numerous and softer CMB photons, with temperature T CMB (z) = 2.73(1 + z) K and mean This secondary HE emission is dubbed "Pair Echoes", and it arrives with a characteristic time delay with respect to the primary HE emission due to the pairs being deflected by the weak intergalactic magnetic field (IGMF) present in cosmic voids that are much less dense in comparison to filaments and clusters (much higher and highly structured magnetic fields are expected in cosmic filaments (B ∼ 10 −9 -10 −7 G [79]) and galaxy clusters (B ∼ 10 −7 -10 −6 G [80]), where the secondary pairs are expected to produce synchrotron pair echoes [81] with ). The pairs IC cool over a characteristic distance where U CMB (z) = aT 4 CMB (z) is the CMB radiation energy density, and a is the radiation constant. Assuming that the pair front expands spherically over a distance λ IC,cool with particles at a typical LF γ e , the radial delay suffered by the secondary HE emission with respect to the primary one is of the order t delay = (1 + z)λ IC,cool (1 − β e )/β e c ∼ (1 + z)λ IC,cool /2γ 2 e c for β e 1 when γ e 1. The pair echo will also be temporally smeared out but over a much larger angular time where θ 1 and λ tot = λ γγ + λ IC,cool , due to light travel time effects over the θ ∼ 1/γ e angular size of the emission region centered at the observer's line-of-sight. Another angular delay is caused by the deflections of the pairs in the intergalactic magnetic field IGMF; ref. [82]. If the coherence length scale of the IGMF is r IGMF < λ IC,cool , then the root mean square angular deflection is ( θ 2 where κ B is an order unity factor that depends on the spectrum of the magnetic field as a function of the coherence length [74,75], and r L = γ e m e c 2 /eB IGMF is the Larmor radius. The corresponding angular time over which the pair echo will be smeared is t ang,B (1 + z) θ 2 B,def λ tot /2c. For extremely energetic pairs with γ e 1, the two timescales, t ang and t ang,B , can become smaller than t VHE , the duration of the primary VHE emission (which could be either prompt and/or afterglow). In addition, for γ e 1 the mean free path for VHE γ-ray photons (λ γγ ) can become smaller than the cooling distance (λ IC,cool ) of the produced pairs, in which case λ tot ≈ λ IC,cool and t ang ≈ t ang,IC = (1 + z)λ IC,cool /2γ 2 e c. Therefore, the correct timescale over which the pair echo signal will be smeared out is t ang + t ang,B + t ang,IC + t VHE ∼ max(t ang , t ang,B , t ang,IC , t VHE ). In the top panel of Figure 1, we show the different timescales as a function of the particle Lorentz factor γ e . Only at very large γ e does the timescale t ang,IC dominate t ang due to a sharp decline in λ γγ caused by the sharp rise in the number density of target CMB photons for γ-rays with E >TeV. Other timescales are the angular time due to IC cooling (∆t IC ) and deflection of the produced pairs by the IGMF (∆t B ; shown for two different IGMF magnetic field strengths) and the angular time associated to the mean free path over which the VHE γ-ray photons produce pairs (∆t A ). The thick solid line highlights the dominant timescale for a given particle Lorentz factor γ e 1.25 × 10 6 (E echo /GeV) 1/2 ⇔ E echo 0.64γ 2 e,6 GeV. Figure from [71] (©AAS. Reproduced with permission.). (Bottom-left) The (observed) model primary and secondary (pair-echo) VHE γ-ray spectral fluences (E γ φ γ = F γ dt, for flux density F γ ) from GRBs at different redshifts (top to bottom): z = {0.1, 0.3, 0.5, 1, 3, 5}. The intrinsic primary spectrum is assumed to be a broken power-law: γ is the peak photon energy and E max γ is the intrinsic high-energy cutoff. The intrinsic spectrum in the figure assumes α = 1, β = 2.2, E b γ = 300 keV, and E max γ = 10 TeV. Figure from [73] (©AAS. Reproduced with permission.). (Bottom-right) The observed pair-echo spectrum shown for different IGMF strengths (with coherence length scale r IGMF = 100 pc) and at different times (t obs = 10 2 s, 10 4 s, 10 6 s) for a source at a fixed redshift of z = 1. The primary prompt emission spectrum is assumed to be a power-law with photon index β = 2.2 above the peak energy E b γ = 500 keV with a cut-off energy of E max γ = 10 TeV, where the prompt GRB has a duration of t GRB = 50 s and luminosity L γ,iso = 10 53 erg s −1 . Figure from [74] (©AAS. Reproduced with permission.).
In the bottom-left panel of Figure 1, we show example model fluence spectra of the primary and secondary (pair-echo) spectra that can be observed from GRBs at different redshifts. The maximum energy of the intrinsic GRB spectrum is assumed to be E max γ = 10 TeV; therefore, the maximum energy of the produced pairs is ∼5 TeV. As a result, the energies of IC scattered CMB photons can reach ∼100 GeV, but these photons may also get absorbed en route to us. Above ∼100 GeV, the contribution from IC upscatted CIB photons becomes important, producing an additional bump in the spectrum. The bottom-right panel shows the pair-echo spectrum at different apparent times and for different IGMF strengths from a source at a fixed redshift z = 1. For a given IGMF, the flux at high energies decays much more rapidly with time compared to the hard power-law at low energies. This is a result of shorter IC cooling times and shorter delay times t ang,B for pairs with larger γ e . Since weaker IGMFs have shorter t ang,B times, the flux is higher initially but decays much faster in comparison to stronger fields that have longer t ang,B times [74].
One of the main advantages (see Section 6 for pair echoes as probes of the IGMF) of detecting pair echoes is that it offers the only way to reconstruct the primary VHE emission from GRBs which would otherwise be attenuated due to VHE photons pair producing on EBL photons.

High-Energy γ-Rays From Pion Decay
Two HE photons are produced directly in the decay of a neutral pion π 0 → 2γ, in which each photon escapes with an energy E γ = m π 0 c 2 /2 67.5 MeV in the rest frame of the pion that is moving with LF γ π 0 in the fluid frame. These photons are then detected with energy E γ (1 + z) ∼ Γγ π 0 E γ 7Γ 2 GeV. Neutral pions can be produced via the following collisional processes between protons (p), neutrons (n), and photons (γ): The most important of the above hadronic energy-loss mechanisms is the photohadronic process, where a photon interacts with a proton, at a threshold photon energy of E γ,th = (m π + m 2 π /2m p )c 2 150 MeV in the rest frame of the proton, to produce a pion, e.g., [83]. When a typical γ-ray photon with energy E γ = (1 + z)E γ /Γ interacts with a proton in the flow having LF γ p , the scattering cross section for the ∆ + resonance peaks when the energy of the photon in the proton's rest frame is The above three collisional processes also produce charged pions (π + and π − ), which then decay to muons that further decay to produce electrons and positrons that can then produce HE synchrotron photons. The most important for producing HE photons is again the ∆ + resonance that also yields where µ + is the anti-muon and ν µ andν µ are its neutrino and anti-neutrino, and ν e is the electron neutrino. Approximately 20% of the proton's energy goes into π + , which is further equally distributed between the pion's decay products [85]. This produces a high-energy The photo-hadronic process, if operating in GRBs, opens up prospects for detecting high-energy (∼ 10 14 eV) neutrinos by km-scale ground-based detectors [85], e.g., IceCube [86]. The intrinsic ratio between muon and electron neutrinos at the source is expected to be 2:1 (with no τ neutrinos), but vacuum oscillations between the three neutrino flavors (ν e , ν µ , ν τ ) may yield equal distributions at Earth. The intrinsic ratio at the source can be different from 1:2:0 when the neutrinos are produced inside a star, e.g., a jet that propagates inside a blue supergiant.
In this case, resonant flavor oscillations in matter due to the Mikheyev-Smirnov-Wolfenstein effect will alter the intrinsic ratios at the source [87]; therefore, after vacuum oscillations, the ratios observed at Earth will also be different from 1:1:1. The intrinsic flavor ratios can further be modified at high energies due to electromagnetic and adiabatic energy losses of muons and pions [88] as well as due to matter oscillations [87,89] at the source. This would lead to energy-dependent, unequal flavor ratios measured at Earth.
Detection of neutrinos from GRBs can only happen for very bright GRBs with γ-ray fluences 10 −4 erg cm −2 [90][91][92]. Thus, far neutrino searches by IceCube have come out empty, even in the case of very bright GRBs, e.g., 130427A [93], and deeper upper limits have offered strong constraints on GRB physics and neutrino production therein [94,95].
Apart from the e ± -pairs produced in the photo-hadronic cascades, additional secondary e ± -pairs result from γγ → e + + e − that can have important effects on both the low and high-energy parts of the spectrum [96,97]. Such pair cascades can also be important in other hadronic scenarios, namely proton synchrotron emission as discussed in Section 2.2. The injected pair spectrum in this case has d ln n e /d ln γ e −2 which yields a relatively flat νF ν synchrotron spectrum. At low energies, i.e., below the peak of the sub-MeV Band component, the synchrotron emission from secondary pairs might dominate and make the spectrum too soft when compared with observations (see bottom-left panel of Figure 4). However, if the secondary pairs are stochastically accelerated (or heated) by MHD/plasma turbulence, then a low-energy photon index of α ∼ −1 that matches observations can be produced [98]. Above the Band-component peak energy, the spectrum is modified due to IC scattered emission by the secondary pairs.

High-Energy γ-Rays from the Bethe-Heitler Process
The Bethe-Heitler process is a photo-hadronic interaction in which the e ± -pairs are produced directly, The differential cross section for this process [99,100] strongly depends on the angle θ ± between the incoming photon and the outgoing e ± in the proton's rest frame. It peaks sharply near θ ± ∼ 1/γ ± , where γ − (γ + ) is the LF of the electron (positron) in the proton's rest frame. When the proton's LF in the jet comoving frame, γ p , is much larger than that of the produced pairs in the proton's rest frame, with γ p γ ± 1, then the pairs are produced with LF in the jet comoving frame of [101]. For typical prompt emission spectral peak energies of E pk m e c 2 /(1 + z), the Bethe-Heitler process is less efficient (by a factor of ∼ 10 2 ) in producing pairs compared to the ∆ + resonance when the LF of produced pairs is γ ± 10 6 . However, for γ ± 10 3 , it can be much more efficient, while for 10 3 γ e 10 6 its efficiency depends on the spectral index of the prompt emission [101]. The high-energy pairs produced in the process can then give rise to HE to VHE photons via synchrotron or IC emission.

GRB Prompt HE Emission-Observations vs. Theory
HE emission in the energy range of (0.1-100) GeV has been detected by the Fermi-LAT in more than 170 GRBs [15]. Prior to Fermi, emission in this energy range (but below ∼20 GeV) was also detected by the Energetic Gamma Ray Experiment Telescope (EGRET) aboard the now defunct Compton Gamma-Ray Observatory (CGRO) from a handful of GRBs [102][103][104][105].
In most of the Fermi-LAT GRBs and that detected by EGRET, the broadband prompt emission spectrum is described by a single Band-like spectral component, generally peaking in the (0.1-1) MeV range and also extending to high energies. In rare cases, it shows a clear spectral cutoff in the (20-350) MeV energy range that is interpreted as a result of the opacity of HE photons to γγ-annihilation within the source [17,106]. This may be intrinsically more common, the low observed fraction being a matter of sensitivity, as such a cutoff appears in 20% out of a sample of bright Fermi Gamma-ray Burst Monitor (GBM) bursts when performing a joint GBM-LAT spectral fit [16]. Moreover, many of the bright Fermi-LAT GRBs show a second spectral component, in addition to the softer Band-like component, that dominates the HE emission [107][108][109], typically well fit by a power law and sometimes showing a cutoff (also likely due to intrinsic γγ-opacity). A similar HE spectral component was also seen by EGRET [105]. Overall, the HE emission seen by the Fermi-LAT shows three remarkable features [15]:

1.
Extra HE spectral component: An extra power-law spectral component that extends to high energies and which is distinct from the typical sub-MeV Band component appears in several bright LAT detected GRBs.

2.
Delayed onset: The onset of this HE emission is delayed relative to the softer γ-rays near the spectral peak, with typical delays of a few to several seconds (0.1 del t del 40 s) for long-soft GRBs and a few tenths of a second (0.05 s t LAT 1 s) for short-hard GRBs.
In the following, we discuss possible origins of the HE spectral component and its delayed onset with respect to the sub-MeV emission (also see, e.g., [113] for a review).
There are two main emission regions from where the HE spectral component can be produced. The first is internal to the outflow in which the emission arises due to dissipation of kinetic energy, e.g., via internal shocks, or magnetic energy, e.g., due to magnetic reconnection, and it occurs at smaller radii before the outflow is significantly slowed down by its interaction with the circumburst medium. In this case, the emission is expected to be highly variable, with t v /t 1 where t v is the variability timescale, and correlated with the sub-MeV prompt emission, which is seen in all cases (as in this case, the two arise from the same outflow, albeit possibly at different radii). The second region is the external forward (afterglow) shock in which case the emission is produced by shock-heated swept-up circumburst medium. In contrast to the prompt emission, the lightcurve is expected to be much smoother, with t v ∼ t, and decaying after its peaks at t t GRB . Such behavior was also observed in many cases. In many LAT GRBs, there is initially a variable GeV emission followed by a smooth tail with a spectral change in the transition, suggesting a transition between prompt and afterglow GeV emission. Furthermore, upon closer inspection, in many cases the delayed onset is caused by the fact that the first spike in the prompt GRB lightcurve is missing at ∼ GeV energies, and only subsequent spikes appear in ∼ GeV and coincide with those at sub-MeV energies.

Delayed Onset of the Fermi-LAT HE Emission
In both the long-soft and short-hard GRBs detected by Fermi-LAT, the HE emission is generally delayed by t del ∼ (0.1-40) s in the former and t del ∼ (0.05-1) s in the latter. While formally t del reaches values as high as 10 4 s in rare cases of both populations [15], these are mostly cases where the GRB was outside the LAT FoV at the time of the GRB trigger and likely do not have a similar physical origin. In the majority of GRBs, the onset of LAT HE emission occurs before the softer prompt γ-ray emission recorded by Fermi-GBM is over. A number of different scenarios have been proposed to explain the delayed onset, which we briefly discuss below.

Forward External Shock Emission
The shock-heated electrons behind the forward shock radiate synchrotron photons that produce the broadband afterglow emission, whose lightcurve peaks at the apparent time (assuming a thin-shell case, for which t dec > t GRB ) In this scenario [112,[114][115][116][117], t dec is the relevant timescale to explain t del . Furthermore, for R > R dec , the proper velocity of the blast wave starts to decline as u(R) = Γ(R)β(R) ∝ R −(3−k)/2 as more mass is swept up, and the dynamical evolution of the blast wave becomes self-similar [67]. For an adiabatic (constant energy with negligible radiative losses) relativistic spherical blast wave, the flux density for ν > max(ν m , ν c ), the frequency regime relevant for HE afterglow emission, scales as F ν (t) ∝ ν −p/2 t −(3p−2)/4 for t > t dec [24,25]. If the blast wave is radiative (a short-lived early phase where its energy decreases over time due to radiative losses), the flux density has the scaling F ν (t) ∝ ν −p/2 t −(6p−2)/7 [24]. (Even for fast cooling only a fraction e of the internal energy generated at the afterglow shock is radiated away, and a similar fraction of the total energy is radiated over each dynamical time, so the blast wave may be far from being fully radiative as assumed in this scaling.) Evidently, the forward shock emission generally obeys a closure relation, whereby the temporal and spectral indices are coupled by virtue of their dependence on the electron energy distribution power-law index p. Kumar and Barniol Duran [115,116] showed that three Fermi-LAT GRBs obeyed this closure relation of an adiabatic blast wave, with p = 2.4 ± 0.06 (GRB 080916C) and p = 2.2 ± 0.2 (GRB 090510, 090902B), that yielded d log F ν /d log t = −(3p − 2)/4 = 1.15 − 1.3 consistent with the observed value to within 1-σ uncertainty. An additional argument in favor of this scenario is that the LAT emission lightcurve shows a very smooth decay, which is expected for afterglow emission. A caveat here is that this applies mainly to the long-lived LAT emission at t t GRB , whereas the variable ∼ GeV emission seen at t ≤ t GRB in bright LAT GRBs cannot be afterglow emission and is most likely prompt emission (especially when it is temporally correlated with MeV spikes in the prompt GRB lightcurve, e.g., [107][108][109]118,119]) The temporal evolution of the observed (isotropic-equivalent) luminosity is another useful probe for the origin of the HE emission. From energy conservation, E ∝ Γ 2 R (3−k) , and with Γ ∝ R (k−3)/(1+δ) where R ∝ Γ 2 t, the time evolution of the isotropic-equivalent energy can be obtained, E iso ∝ t [(δ−1)(3−k)]/(7+δ−2k) [120]. For an adiabatic (δ = 1) blast wave, E iso ∝ t 0 , as it should be, and for radiative (δ = 0) blast wave E iso ∝ t (k−3)/(7−2k) . The observed luminosity then follows with L γ,iso ∝ E iso /t ∝ t [(δ−3)(3−k)−(δ+1)]/(7+δ−2k) , which yields L γ,iso ∝ t −1 for the adiabatic case and L γ,iso ∝ t (3k−10)/(7−2k) for the radiative case. The left panel of Figure 3 shows the radiative afterglow model fit to the LAT lightcurve of GRB 080916C [112]. This agreement presents a strong argument in favor of the synchrotron afterglow origin of the late-time LAT HE emission. However, the main LAT peak is too sharp to arise from the afterglow onset corresponding to the outflow deceleration time, and instead matches the second ∼ MeV peak, so it is more likely associated with the prompt GRB emission, while the temporally smoother afterglow GeV emission likely starts dominating later, after a few tens of seconds, with a rather shallow decay slope (∼ t −1 ).
There are two major hurdles for this scenario. First, many LAT GRBs show a peak in the GeV emission while the prompt emission is still active, which is difficult to explain with synchrotron emisson from the external forward shock. In the thin-shell afterglow shock scenario [66], the peak of the HE emission will occur at t = t dec = (1 + z)R dec /2cΓ 2 0 which is always larger than the duration of the prompt GRB emission, T GRB = (1 + z)∆ 0 /c, given by the shell crossing time of the ejecta shell of initial thickness ∆ 0 . Alternatively, in the thick-shell case, t dec ∼ T GRB . Second, this model cannot explain the detection of VHE photons at late times when t > t dec where the detected photons have energies much larger than E syn,max [121]. Both of these arguments suggest that yet another mechanism might be responsible for the LAT emission.

Inverse-Compton GeV Flash
The shock-heated electrons behind the forward shock at radius R can be Compton-cooled by prompt emission ∼ MeV photons emitted at a smaller radius R prompt R as the radiation front overlaps with the blast wave [51,52,56]. When the prompt emission photons travel ahead of the blast wave, a small fraction is scattered by the yet unshocked electrons in the CBM at large angles from the radial direction. The scattered photons then produce e ± -pairs via γγ-annihilaton on the radially expanding (collimated) prompt emission radiation front. The created pairs further scatter the prompt photons, causing exponential pair-creation and the resultant high multiplicity (with M ± 10 5 ) pair-loading of the CBM ahead of the forward shock [122][123][124][125][126]. Scattering of the prompt radiation by the pair-loaded CBM also imparts momentum to the pairs and pre-accelerates them to a typical LF γ pre = (1 − β 2 pre ) −1/2 < Γ bw , where Γ bw = (1 − β 2 bw ) −1/2 is the bulk LF of the blast wave, i.e., material just behind the forward shock. As the blast wave sweeps up the pair-loaded CBM, with a relative LF Γ rel = Γ bw γ pre (1 − β bw β pre ) ≈ Γ bw /γ pre (1 + β pre ), the shock-heated pairs are thermalized with γ th ∼ Γ rel when M ± 10 3 . This model assumes that only a small number of particles are accelerated into a power-law energy distribution and most of the energy resides with the quasi-thermal pairs. The radiative efficiency of the shock-heated pairs is almost 100% during the GeV flash; therefore, the blast wave does not start to evolve adiabatically until all the prompt emission photons have overtaken it.
In this scenario, the peak of the LAT emission occurs at t pk = (1 + z)R pk /2cΓ 2 pk where R pk ∼ 10 16 cm (in the case of GRB 080916C) is the radius where the LF of the electrons behind the forward shock is γ th,pk ∼ 50 so that the IC scattered emission peaks in the GeV energy range. At R < R pk , the contrast between Γ bw and γ pre is small and therefore γ th < γ th,pk . This contrast grows over larger radii and γ th = γ th,pk at R = R pk , and for R > R pk , the contrast is much larger which yields γ th > γ th,pk and produces VHE emission at TeV energies. The right panel of Figure 3 shows the model fit to the LAT lightcurve of 080916C from Beloborodov et al. [56].

Synchrotron Emission from Protons Accelerated at the External Forward Shock
In the hadronic scenario, a proton-synchrotron emission model with a strong comoving B-field (B ) can explain the delayed onset of the LAT emission [40,127,128]. Just like electrons, protons are also accelerated at the external blast wave to energies where they can radiate ∼ GeV to TeV synchrotron radiation. This radiation is further processed into e ± -pairs via γγannihilation where the produced pairs then radiate sub-GeV synchrotron photons. The onset of HE emission is delayed due to two effects. First, protons are accelerated over the Larmor time to achieve a maximum LF γ M,p (see Section 2.2 for definition), which causes a delay of at least This should also be a lower limit on the variability time, as the local emission cannot turn on or off faster than this. Second, as shown in the model put forth by Razzaque et al. [40], it takes a finite amount of time for the peak of the proton synchrotron radiation spectrum, which peaks at higher energies at early times, to move into the LAT energy range.
This scenario requires a strong magnetization of the shocked material downstream of the blast wave to explain the delays in the LAT emission onset. A major weakness of this model is that it is radiatively inefficient and therefore requires a large amount of energy in accelerated protons that must be injected with minimum LF of γ m,p 10 6 [42,101]. Furthermore, as the proton synchrotron cooling break sweeps across the observed energy band, the spectral index should change from d ln F ν /d ln ν = (1 − p)/2 to −p/2, where p is the power-law index of the proton energy distribution, n p (γ p ) ∝ γ −p p for γ p,min ≤ γ p ≤ γ p,max . However, no such spectral change has been observed in the delayed LAT emission. Although it is possible that this spectral component has only been observed at energies above the cooling break due to the limiting sensitivity of the Fermi-LAT at high energies, it would be too much of a coincidence to have happened in all LAT bursts that show delayed emission.

SSC Emission
The delay time of the LAT emission in the SSC scenario depends on the time it takes for the IC-scattered radiation field to build up in the LAT energy band. That depends on the temporal evolution of the Compton-y parameter which must become larger than unity for IC scattering to become the dominant particle cooling mechanism. Detailed one-zone numerical simulations [128][129][130] of prompt GRB emission show that under certain conditions SSC emission in the LAT energy range can be delayed with respect to the sub-MeV synchrotron component due to the time it takes to build up the seed synchrotron photon field in the emitting region (of the order of its light crossing time). However, in many cases, the temporal delay is insufficient to explain the observed ones and remains limited to t delay < t v , where t v is the variability timescale. This effect would also lead to a systematic delay of the GeV emission w.r.t the sub-MeV emission for each spike in the prompt lightcurve. In practice, the observed delay typically reflects the first spike being absent in the GeV, with subsequent spikes coinciding in MeV and GeV.

Distinct HE Spectral Component
Many bright Fermi-LAT GRBs show a distinct HE spectral component in addition to the Band-like spectrum, where the latter represents the canonical prompt emission spectrum peaking in the ∼(0.1-1) MeV energy range. This additional component has been modeled as a powerlaw, sometimes with a high-energy cutoff, in addition to the Band component. Such a component was required by the data in GRB 090227B [131] [133]. The first such detection of an additional component, however, was made by EGRET in GRB 941017 [105]. In most cases, the additional power-law component extends to low energies (∼ few keV) and exceeds the Band component below a few tens of keV, forming a low-energy excess. At high energies, this power-law component is detected up to ∼ 10 −0.5 -10 1.5 GeV with photon index α PL ∼ −1.9 to −1.5. In some cases (e.g., GRB 090926A and GRB 190114C), however, this component shows a high-energy turnover at early times before becoming a strict power-law as the spectral break moves above the LAT energy window. Example time-integrated and time-resolved spectra for the short-hard GRB 090510 [108] and long-soft GRB 090926A [109] are shown in Figure 4.  [108], ©AAS. Reproduced with permission.) and long-soft GRB 090926A (right; [109], ©AAS. Reproduced with permission.). In both cases, the spectrum shows a low-and high-energy excess which is fit by a power-law or cutoff power-law component that has a distinct origin from the main Band component. Bottom: Theoretical modeling of the spectrum of GRB 090510 using a hadronic scenario with photo-hadronic cascades (left; [97], ©AAS. Reproduced with permission) and that of GRB 080916C using the pre-accelerated and pair-loaded ISM in the afterglow model (right; [56], ©AAS. Reproduced with permission.).
The low-energy excess presents a challenge for leptonic scenarios, e.g., SSC, that can only explain the excess emission above the peak of the sub-MeV emission and remain subdominant to the main synchrotron (Band) component at low energies. Such an excess can be produced in hadronic models featuring direct proton synchrotron emission [40] or photohadronic cascades [96]. Theoretical modeling of the prompt spectrum of GRB 090510 using the latter model is shown in the bottom-left panel of Figure 4 for two different ratios of U B /U γ = {10 −3 , 10 −1 }, where U B is the energy density of the comoving B-field, and U γ is the energy density of the Band component (shown with a black dashed curve; not produced by the same secondary pairs that IC scatter the Band component). In this case, the low-and high-energy excesses are given by synchrotron and inverse-Compton emissions from the secondary e ± -pair cascades (shown by thin red curves without absorption). The peak around ∼ GeV arises from absorption due to γγ-annihilation.
Alternatively, models that attribute the origin of the additional power-law component to afterglow emission also find it challenging to explain the low-energy excess. As mentioned earlier, in simple external forward shock models, e.g., [115,116], the dominant contribution from the afterglow occurs at t = t dec t GRB ; therefore, these models cannot explain the origin of the additional power-law component in the prompt emission spectrum. In the model put forth by Beloborodov et al. [56], the inverse-Compton emission (as shown in the bottom-right panel of Figure 4) from shock-heated electrons behind the forward shock that sweeps up pre-accelerated and pair-loaded ISM remains sub-dominant at low-energies and therefore cannot explain the low-energy excess.
The low-energy excess is not always modeled as a low-energy extension of the power-law component that dominates at high energies above the Band component. In some cases, it has been interpreted as a combination of a Band plus photospheric (quasi-thermal) components that jointly produce this excess, e.g., [134][135][136]. This degeneracy produced by different spectral models describing the same data equally well further adds to the complexity of the underlying emission mechanism.

Long-Lived HE Emission
At an early time, when the sub-MeV prompt emission is still active, the HE emission detected by the Fermi-LAT shows significant temporal variability, which in many cases, e.g., [118] is correlated with the sub-MeV emission. This can be attributed to the HE emission having originated in the same spatial region as the sub-MeV component. After the prompt emission ceases and the afterglow commences, the HE emission shows a temporally smooth and long-lasting decay with d log F ν /d log t ≈ −1 and a standard deviation of 0.8 (in some cases a broken power-law fit to the lightcurve is statistically preferred) [15]. This is often referred to as the LAT extended emission (see Figure 2 and left panel of Figure 5). Since it lacks the short timescale variability and lasts much longer, it is naturally interpreted as the HE tail of the afterglow emission from the external forward shock.
As discussed in Section 3.1, synchrotron afterglow emission from an adiabatic [115,116] (or possibly in some cases from a radiative [112]) blast wave can very well explain the temporal decay index of the LAT extended emission. The agreement with multi-wavelength observations (see, e.g., [116]) suggests that it certainly is a strong candidate for the late-time extended emission, even though this model may not be the correct description for the early time (with its delayed onset and low-energy spectral excess if produced by the same non-thermal component) LAT emission.
The one problem this scenario faces is the detection of (V)HE photons at late times. In several GRBs, HE photons with observed energy E 10 GeV arrive at t ∼ 10 2 -10 3 s, much after the cessation of the prompt GRB emission [15]. The origin of such photons using the standard leptonic synchrotron afterglow scenario is difficult to explain as they violate E syn,max . According to this limit, to produce photons with energy (1 + z)E 10 GeV in the cosmological rest-frame of the source would require bulk Γ > 10 2 at late times, which is nearly impossible. Therefore, our assumptions regarding particle acceleration at shock fronts must be revised (see further discussion in Section 4 below). The alternative is SSC emission which would manifest as an additional spectral component in the LAT energy band and would also be detected at very high energies. flects the size of the emitting region, and that the MeV and GeV emissions around the time of the 73-GeV photon at T0 + 19 s are cospatial, the requirement that the optical depth due to gg opacity be less than 1 then implies that the minimum bulk Lorentz factor is Gmin ¼ 455 þ16 −13 . Here, a SBPL fit to the GBM spectrum in the interval 11.5 to 33.0 s (table S1) and a minimum variability time scale of 0.04 T 0.01 s are used (17). The cospatial as-sumption is, however, questionable given the different time histories in the MeV and GeV emission. Moreover, values of Gmin that are smaller by a factor of 2 to 3 can be realized for models with time-dependent g-ray opacity in a thin-shell model (21).
The delayed onset of the LAT-detected emission with respect to the GBM-detected emission is an important clue to the nature of GRBs (13). For GRB 130427A, the LAT-detected emission becomes harder and more intense after the GBMdetected emission has faded (Fig. 3). This suggests that the GeV emission is produced later than the keV-MeV emission and in a different region. In particular, if the keV-MeV emission comes from interactions within the outflow itself, the GeV emission arises from the outflow's interactions with the circumburst medium.

Constraints on Bulk Γ
Since GRBs are extremely luminous sources, a typical photon near the νF ν peak with energy E ∼ E pk ∼ m e c 2 would see a large optical depth τ γγ 1 to pair production (γγ → e + e − ) [138]. For a Newtonian source, this would imply a huge compactness ≡ σ T U γ R/m e c 2 (Thomson optical depth of pairs if all photons pair produce), where U γ is the radiation field energy density, which would result in a nearly blackbody spectrum, in stark contrast with the observed prompt GRB non-thermal spectrum. The solution to this so-called compactness problem, is that the emission region must be moving towards us at ultra-relativistic speeds with bulk LF Γ 10 2 [139][140][141]. The observed energy E cut where the prompt GRB spectrum would display a cutoff due to γγ-annihilation is sensitive to Γ; therefore, an observation of such a cutoff yields a direct estimate of Γ, which is difficult to obtain otherwise. Spectral cutoffs have only been observed in a handful of GRBs, e.g., [17,106,109], and in most cases, the spectrum above the νF ν peak is a featureless power law extending to some E max , the maximum photon energy detected by the instrument. In such cases, only lower limits can be placed on the bulk LF with Γ > Γ min . The maximum possible bulk LF for a given E cut is given by Γ max = (1 + z)E cut /m e c 2 , and the true bulk LF is Γ = min[Γ min , Γ max ], e.g., [37].
In several bright GRBs detected by Fermi-LAT, estimates of Γ min have been obtained for a given E max using a simple one-zone analytic formalism, e.g., [33][34][35], with Γ min ≈ 900 for GRB 080916C [119], Γ min ≈ 1200 for GRB 090510 [108], and Γ min ≈ 1000 for GRB 090902B [107]. When a more sophisticated formalism [36,142,143] that includes temporal, spatial, and angular dependence of the radiation field, and which is verified by numerical simulations [37], is applied, it yields a Γ min estimate smaller by a factor of ∼2-3 (similar results were obtained by [144,145]). In GRBs that show a high-energy spectral cutoff, bulk LF of Γ ∼ 100-400 have been obtained using detailed numerical models [17].
In addition to its phenomenal prompt emission, the afterglow emission GRB 130427A provided extremely interesting physical insights. In particular, as described below, its temporal and spectral analyses challenge the most widely accepted model for the afterglow phase of GRBs. Its prompt emission lasted about T 90 = 276 ± 5 s (at 15-150 keV; [147]), and after it subsided, the observed emission was clearly dominated by the afterglow [137,147,154], showing a smooth power-law flux decay as well as a typical afterglowlike spectrum (see Figure 5). Moreover, while the reverse shock emission appears to dominate at early times at low frequencies, it does not dominate beyond the optical (even at early times), where the observed emission is dominated by the forward shock all along [154]. The values of the temporal and spectral indices in the power-law segment ν m < ν < ν c (PLS G from [25]) , F ν ∝ ν −0.69±0.01 t −1.30±0.05 , imply an external density profile ρ ext ∝ R −k with k = 1.4 ± 0.2 [137], suggesting that the GRB progenitor star's wind mass loss rate to velocity ratio (Ṁ w /v w ∝ R 2−k orṀ w (t) ∝ v w (t) 3−kt 2−k wheret = R/v w (t) is the wind ejection time prior to the stellar explosion leading to the GRB) slightly decreased towards the end of its life.
Most importantly, the high-energy emission from the afterglow of GRB 130427A was not only detected by the Fermi-LAT for 20 hours (see top-left panel of Figure 5), but also included multiple high-energy photons up to very late times (see bottom-left panel of Figure 5) that were clearly in excess of the maximum synchrotron photon energy, E syn,max . This upper limit on the energy of synchrotron photons is derived by equating the electron acceleration and synchrotron radiative cooling timescales, assuming a single acceleration and emission region [27,29,155,156]. While there was some evidence of E syn,max violation in previous Fermi-LAT GRBs (e.g., [29,119]), in those cases on the one hand the violation was weaker (by a smaller factor and with fewer photons of less statistical significance), and on the other hand a different emission mechanism was a viable alternative explanation. In GRB 130427A, the long-lasting (∼1 day) Fermi-LAT afterglow included a 32 GeV photon after 34 ks, and altogether five > 30 GeV photons after >200 s (with probability > 99.9% of being associated with GRB 130427A). All five significantly exceed E syn,max , by factors of at least 6.25 for k = 0 and 9.20 for k = 2 (using Equation (4) of [29]).
This has led to suggestions that the Fermi-LAT high-energy photons were not synchrotron radiation, but instead arose from a distinct high-energy spectral component [121,157]. All such options require that the high-energy part of the Fermi-LAT detected energy range (above a few to several GeV, depending on the exact time) should be dominated by a distinct spectral component, while lower energies are dominated by the usual afterglow synchrotron spectral component. Such an option was considered by [137], who fit the SED from optical to GeV at two epochs (∼1.5 an 5 days) where observations were also performed by the Nuclear Spectroscopic Telescope ARray (NuSTAR) in the 3-79 keV energy range (see right panels of Figure 5). The spectrum was fit to a detailed synchrotron afterglow model [25], which provided a good fit at both epochs. Moreover, at the first and more constraining epoch (∼1.5 days), the Fermi-LAT flux that is extrapolated by a factor of 2 in time agrees very well with this synchrotron afterglow model. Furthermore, both this, as well as the simultaneous upper limit by the Fermi-LAT [121,137], and more importantly the nearly simultaneous VHE upper limit by the Very Energetic Radiation Imaging Telescope Array System (VERITAS) [158], hardly leaves any room for a distinct high-energy spectral component.
Therefore, this strongly suggests that the late-time Fermi LAT high-energy photons in GRB 130427A are indeed afterglow synchrotron radiation. This provides the strongest direct observational support for a genuine violation of E syn,max by synchrotron photons. As the latter arise from the afterglow forward shock, this challenges our understanding of particle acceleration and magnetic field amplification in relativistic collisionless shocks. In particular, at least one of the assumptions in the derivation of E syn,max must be incorrect, requiring a modification of our understanding of afterglow shock physics.
While this potential problem was known before, these results from GRB 130427A [137] have made it much harder to circumvent (and the VHE upper limit by VERITAS played an important role). A possible solution to this problem may lie in modifying the assumption of a single uniform region where both the acceleration of electrons and radiation from them occurs. Instead, one can allow for a lower magnetic field acceleration region and a higher magnetic field synchrotron radiation region (e.g., [30,159]). Such a situation may arise for diffusive shock acceleration (Fermi Type I) if the tangled shock-amplified magnetic field decays on a short length scale behind the shock front. In this case, most of the high-energy radiation is emitted just behind the shock from where the magnetic field has not decayed significantly, while the highest-energy electrons are accelerated further downstream where the magnetic field is lower [30]. This puzzle is still far from being resolved and poses a serious challenge to our understanding of relativistic collisionless shock physics.

Very-High-Energy (TeV) Afterglow
The detection by EGRET of MeV-GeV photons over ∼90 min from GRB 940217 [103], as well as the hard additional spectral component in the prompt emission of GRB 941017 [105], led to the consideration of SSC prompt [129,[160][161][162][163] and afterglow radiation [44][45][46]59,164] and searches for VHE TeV photons by ground-based detectors. TeV emission is expected to be observed only from relatively nearby GRBs due to absorption of VHE γ-rays by γγ-annihilation on EBL photons from the more distant sources. The Universe starts to become opaque to VHE photons with E 1 TeV for redshifts z 0.08 [165]. Early efforts at detecting VHE radiation from GRBs were made by the Milagro instrument, an extended air shower detector, and hints of VHE photon (3σ) detection from GRB 970417A were found by Milagrito [166] (see [167] for prompt TeV γ-ray emission model for this detection), the smaller and less sensitive prototype detector. Over the last two decades, imaging atmospheric Cherenkov telescopes (IACTs), namely the Very Energetic Radiation Imaging Telescope Array System (VERITAS; [168]), the Major Atmospheric Gamma Imaging Cherenkov (MAGIC; [169]), and the high-energy Stereoscopic System (H.E.S.S; [170]) have been routinely monitoring for VHE radiation from GRBs.
These efforts bore fruit in January of 2019 when MAGIC announced the ( 50σ) detection of VHE (∼0.2-1 TeV) photons from GRB 190114C [171]. Figure 6 shows the multi-wavelength lightcurve and broadband afterglow spectrum of this burst. GRB 190114C had a redshift of z = 0.424, and its prompt emission was detected by several space-based γ-ray instruments [133,172] that measured an isotropic-equivalent (1 keV to 10 GeV) energy release of E iso 3 × 10 53 erg over a duration of ∼25 s (shown by the dashed vertical line in Figure 6). MAGIC detected afterglow VHE γ-ray photons from ∼60 s to ∼2400 s and measured an E iso,TeV 4 × 10 51 erg, which is only a lower limit due to the late start of observations and could be as high as ∼10% of the energy released in softer γ-rays. This was the first time that time-resolved afterglow spectra all the way up to TeV energies were obtained in any GRB observed to date. This naturally has important implications for GRB afterglow physics and overall energetics of the system. After GRB 190114C, a few other GRBs (160821B [174]with a low detection significance of ∼3σ, 180720B [175], 190829A [176], 201216C [177]) were reported to have been detected at sub-TeV to TeV energies by both H.E.S.S and MAGIC (see the reviews by Nava [18], Noda and Parsons [19] for more details).

Key Results Implications
In the following, we briefly discuss the most important implications for GRB physics from the detection and theoretical modeling of ∼ TeV afterglow emission.

IC Emission Is Needed to Explain the VHE γ-Rays
As shown in the right panel of Figure 6, the hardening of the MAGIC-detected VHE spectrum with respect to the LAT detected HE spectrum in GRB 190114C indicates the presence of an additional spectral component. It simply cannot be explained with synchrotron emission from the external forward shock alone. Several works that use analytical/semi-analytical [178][179][180][181][182] and numerical models [183,184] have now been devoted to explaining the ∼ TeV emission as SSC or a combination of EIC and SSC [185]. Significant differences between the (semi-)analytical and numerical models arise due to inclusion and more accurate handling of some of the non-linear processes, such as pair cascades due to internal γγ absorption and KN effects. In the end, the obtained shock microphysical parameters indicate that these bursts are not very different from the ones that are not detected with a VHE component, which may suggest that SSC afterglow emission is rather common. In that case, it becomes important to take into account the energy radiated in the SSC component to understand the global energetics of the bursts. For example, the energy in the SSC component was ∼40% of that radiated in the main synchrotron afterglow component for GRB 190114C [171]. Similar inferences regarding the total energy budget were also drawn before and around the first GeV detections from GRBs by the Fermi-LAT. It was later shown that on average E GeV /E MeV 0.1, and at best the two become comparable for rare individual LAT GRBs. In GRB 190114C the detected TeV emission is from the afterglow. Therefore, it does not affect the prompt GRB energy budget.
An alternative to IC emission that can explain the VHE TeV γ-rays is photohadronic emission, as demonstrated by Sahu and Fortín [186], Sahu et al. [187]. In their model, VHE γ-ray photons are produced via the pγ → ∆ + process, which then produces neutral pions that decay into γ-ray photons (see Section 2.6). The seed photons that interact with the protons can be of synchrotron or SSC origin as produced in the afterglow forward shock.

Constraints on Shock Microphysical Parameters
When fitting afterglow observations, the parameter space is usually degenerate, and unique values of the shock microphysical parameters cannot be obtained. The parameter space consists of E k,iso , the isotropic-equivalent total kinetic energy of the flow, n, the number density of the ISM or its normalization (n ∝ AR −2 ) for a wind circumburst medium, and the shock microphysical parameters e , B , and ξ e . The power-law index p of the particle energy distribution is uniquely determined from the broadband spectrum. It is generally assumed that ξ e = 1, in which case the remaining four afterglow model parameters (E * k,iso , n * , * e , * B ) can be uniquely determined using (ν a , ν m , ν c , F ν,max ), leaving the degeneracy due to ξ e [188] where all values of (E k,iso , n, e , B ) = (ξ −1 e E * k,iso , ξ −1 e n * , ξ e * e , ξ e * B ) fit the data equally well for any m e /m p < ξ e ≤ 1. This degeneracy may possibly be broken when accounting for the emission, absorption, or propagation effects of the thermal electrons [189][190][191][192][193]. In addition to these parameters, multi-wavelength modeling of TeV afterglows can potentially be used to constrain γ M and in turn the acceleration efficiency (κ acc ) of relativistic collisionless shocks. Such constraints can then be used for comparison with first principles PIC simulations of Weibel-mediated collisonless shocks that do predict the value of γ M .
In Table 1, we list the afterglow fit parameters of GRBs that were detected at very high energies. In all cases, the energy deposited in power-law electrons is much larger than that in the shock-generated B-field, i.e., B e . Consequently, this yields the Compton-y parameter larger than unity which results in producing a bright SSC component detected by MAGIC and H.E.S.S. In most works, the afterglow shock microphysical parameters are taken to be constant throughout the afterglow evolution; however, Misra et al. [194] report the possibility of evolving microphysical parameters to explain the long-term radio/mm afterglow of GRB 190114C.

VHE γ-Rays as Electromagnetic Counterparts of Binary NS Mergers
The detection of afterglow TeV γ-rays in long-soft GRBs has opened up the prospect of also detecting VHE emission in short-hard GRBs. Electromagnetic emission coincident with GWs was first detected from the binary NS merger in GW 170817/GRB 170817A. An impressive multiwavelength follow-up by a number of ground-and space-based observatories tracked its peculiar afterglow emission. No ∼TeV γ-rays were detected [196,198,199] for this relatively nearby (∼40 Mpc) event as the relativistic jet was observed off-axis, and none have been detected from other short-hard GRBs. Short-hard GRBs are detected more numerously at redshifts z < 1 with a mean redshift of z ≈ 0.5 in comparison to z ≈ 2 for the long-soft GRBs. Therefore, attenuation of VHE γ-rays by the EBL is not as extreme for the short-hard GRBs (e.g., MAGIC detected GRB 190114C had a redshift of z ≈ 0.42) as it is for the more distant long-soft GRBs. Catching the VHE emission in time from short-hard GRBs will require high sensitivity (due to lower fluences), shorter telescope slew times, as well as a large field-of-view (since catching the prompt emission would require the GRB to be in the field-of-view without slewing). Fast follow-ups of GW triggers by existing (MAGIC, H.E.S.S, VERITAS, HAWC) and upcoming observatories, namely the Cherenkov Telescope Array (CTA; [200]), will play a pivotal role in the next several years.

Studying Non-GRB Physics
Since GRBs are the most electromagnetically luminous events in the Universe, they are observed up to cosmological distances. They emit HE and VHE photons that travel over cosmological distances on the way from the source to us. This can naturally be used to probe various processes involving these photons that may occur along their way to our detectors, which provide unique and valuable information about cosmology or basic physics.

Constraining the Extragalactic Background Light Models with GRBs
The detection of HE photons from distant GRBs proves to be an excellent probe of the extragalactic background light (EBL) [165,[201][202][203][204][205][206][207], which is the cumulative star light emitted in the UV/optical to infrared energy range, i.e., ∼10 −3 eV to 10 eV (∼0.1 µm to 10 3 µm), by all the stars in the Universe. The EBL is difficult to constrain otherwise due to contamination by the zodiacal and Galactic foreground light [208]. After the cosmic microwave background, the EBL is the second dominant component that contributes to the diffuse radiation that pervades entire space. Star light with wavelength 2µm is highly absorbed by dust in the host galaxy, with only a fraction that escapes and contributes to the EBL. On the other hand, the dust re-radiates the star light but adds to the EBL in the infrared. As discussed in Section 2.5, VHE γ-ray photons with E γ 1 TeV interact with the infrared background and produce e ± -pairs, while lower energy photons typically interact with UV/optical/NIR photons emitted directly by stars. The attenuation of the HE spectra of TeV sources, including GRBs, caused by this effect can be used to constrain models of EBL, with the assumption that the intrinsic spectrum can be extrapolated to higher energies using the lower energy part of the spectrum that is not affected by such attenuation.
Before the MAGIC detection of ∼ TeV γ-rays in GRB 190114C, constraints on EBL models were placed using the observations of HE photons from only a few GRBs. For example, the fast evolution and baseline models of Stecker et al. [165] were disfavored at the >3σ level by the detection of a 33.4 GeV photon from GRB 090902B which was at a redshift of z = 1.822 [107]. Observations of higher redshift GRBs at high energies offer a better chance of constraining EBL models. A 13.2 GeV photon was detected from GRB 080916C having a redshift of z = 4.35. The opacity of the Universe to such a photon is shown in Figure 7 that compares different EBL models. The suppression of the ∼ TeV spectrum of a GRB due to the EBL was first clearly seen in GRB 190114C [171,173], as shown in Figure 6.

Probing the Intergalactic Medium B-Field
The space between galaxies is expected to be permeated by a very weak (B IGMF 10 −20 G) magnetic field [209][210][211][212], the origin, strength, and coherence length of which are poorly understood. This intergalactic magnetic field (IGMF) possibly acted as the seed magnetic field in galaxies and galaxy clusters, which was amplified to typical strengths of ∼ µG by a dynamo mechanism as well as flux conserving collapse during their formation. Therefore, its origin predates structure formation in the Universe. Since such field amplification processes are absent in the voids between galaxies, which would have otherwise erased the initial magnetic field properties, study of the IGMF can provide important insights into the origin of the seed field in galaxies, and it can be used to constrain physical processes in the early Universe that may have generated it. Contamination of the primordial IGMF is possible via magnetized outflows from active galaxies and galactic winds driven by star formation.
One of the ways to study the IGMF, albeit indirectly, is by detecting ∼ GeV pair echos (see Section 2.5) created by the IC scattering of CMB photons by e ± pairs that in turn are produced by ∼TeV γ-rays (from an astrophysical source) annihilating with EBL photons [82,213]. This causes several observable effects, e.g., time delays between the TeV and GeV signal [70,71,214] and extended γ-ray haloes [215], that can be used to constrain the properties of the IGMF. Non-observation of GeV γ-rays from persistent sources, e.g., TeV blazars, were used to derive a lower limit of B IGMF 10 −16 G for a coherence length of 10 kpc to 1 Mpc [216,217]. The potential problem with persistent sources, when GeV radiation is observed, is that the pair echo photons can overlap with the intrinsic emission causing contamination. In that regard, GRBs serve as better probes since there could be a clear temporal separation between the short-lived prompt emission and the detection of the longer-lived GeV pair echo.
Before the detection of ∼ TeV γ-rays from GRB 190114C, constraints on the IGMF were obtained using GRB 130427A from (i) the VERITAS upper limits at 100 GeV at 0.82 days, (ii) Fermi-LAT detection of a 32 GeV photon at 34.4 ks post-trigger, that could not be explained as synchrotron afterglow radiation, and (iii) Fermi-LAT upper limits at GeV energies at late times [218]. From the non-detection of GeV emission in Fermi-LAT observations of the location of GRB 190114C over a period up to 3 months, Wang et al. [219] derive a lower limit of B IGMF ≈ 2 × 10 −20 (λ B /0.1 Mpc) −1/2 G for the coherence length of λ B < 0.1 Mpc. On the other hand, Dzhatdoev et al. [220] argued that the Fermi-LAT flux upper limits were insufficient in constraining the IGMF and that the results of Wang et al. [219] are in error due to overestimation of the pair echo intensity.

Lorentz Invariance Constraints
One of the tenets of the special theory of relativity is that the speed of light is a Lorentz invariant, i.e., it is the same in two Lorentz frames regardless of their relative motion. In particular, it is independent of the photon energy (or wavelength), which is frame-dependent. The underlying assumption is that this remains true at all length scales (or wavelengths) in nature, no matter how small. However, quantum effects are expected to strongly affect the nature of spacetime at the Planck scale, corresponding to a length scale of l Planck = √h G/c 3 ≈ 1.62 × 10 −33 cm, or equivalently at an energy scale of E Planck = M Planck c 2 = √h c 5 /G ≈ 1.22 × 10 19 GeV. Such effects can possibly lead to Lorentz invariance violation (LIV) where the speed of light changes with photon energy [221][222][223]. In some theories of quantum gravity, this effect can lead to dispersion as the photon propagates in the vacuum of space, such that its speed varies as where E γ is the photon energy, E QG is the quantum gravity energy scale (expected to be ∼ E Planck ), and the sign ambiguity depends on the dynamical framework. Since E QG is much larger than observed photon energies, the change in velocity is rather minute. However, this effect can accumulate over cosmological distance scales D, which makes GRBs ideal probes of LIV.
This LIV effect would manifest as an arrival time difference between photons having different energies, with ∆t LIV ≈ ±(∆E γ /E QG )D/c and ∆E γ = E γ,high − E γ,low . Since GRBs show temporal variability as short as ∆t ∼ ms, their large distances D ∼ 10 28 L 28 cm can probe quantum gravity energy scales approaching the Planck scale using ∼ GeV photons, E QG ≈ 3 × 10 19 (∆E γ /GeV)D 28 ∆t −1 −2 GeV ∼ E Planck . This technique was employed in the case of short-hard GRB 090510 that emitted a 31 GeV photon 0.829 s after the burst onset and which coincided in time with the last of the seven pulses comprising the prompt emission [118]. By using an unbinned analysis, in both energy and time, testing different dispersion coefficients that would yield velocity differences of ∆v ≈ E γ /E QG , and then maximizing ∆v so that it yields the sharpest lightcurve, Abdo et al. [118] obtained a lower limit of |∆t/∆E γ | < 30 ms GeV −1 (at the 99% confidence level) or equivalently E QG > 1.2E Planck . Using the same data for the short-hard GRB 090510, this limit was somewhat improved using a more refined analysis [224], and a Planck-scale limit was also derived on space-time fuzziness and stochastic LIV [225], which are motivated by the notion of space-time foam.

Outstanding Questions
Even with only a handful of detections in the VHE domain, new questions have emerged. We briefly highlight some of the fundamental questions that may be resolved with future TeV detections as well as improved modeling of radiation processes.
(a) What makes GRBs TeV bright? : All TeV bright GRBs are also very bright in prompt γ-rays as well as in their X-ray afterglow emission. In fact, apart from GRB 190829A, the rest of the TeV bright GRBs have high prompt γ-ray fluences that put them among the top 1%, see Figure 1 of [19], as is also evident from their high E γ,iso 10 53 erg from Table 1.
Although not all MeV-bright GRBs were observed at TeV energies, it begs the question why no TeV emission was detected from, e.g., 130427A (one of the most energetic GRBs with E γ,iso ≈ 1.4 × 10 54 erg) by VERITAS and HAWC, and whether we would have seen TeV γ-rays from all such GRBs. The majority of the highly energetic GRBs are also more distant with z > 1.0 see, e.g., Figure 3 of [19], which makes it challenging to detect their TeV emission due to suppression by γγ-annihilation on EBL photons. Internal absorption due to γγ-annihilation of IC photons (that produce the TeV component) on the seed synchrotron photons (that produce the X-ray afterglow) can also become important [184] and may perhaps be enough in some bursts to significantly suppress the VHE emission. Detailed semi-analytic numerical models including the effects of pair cascades and Klein-Nishina suppression that can explain the multi-wavelength spectral and temporal evolution may shed more light on the properties of the emission region. (b) What causes the delayed onset of the Fermi-LAT emission?: The delayed onset of the HE emission w.r.t the ∼ MeV prompt γ-rays as seen by the Fermi-LAT has been interpreted as the peak of the standard afterglow emission [115,116], IC GeV flash from the preaccelerated pair-rich circumburst medium swept up by the external forward shock [52,56], acceleration time of protons in hadronic emission scenarios [40], and the timescale over which the SSC radiation field builds up [128]. The latter two scenarios have some difficulties with, respectively, the global energetics and limitation on the delay duration, and the former two struggle with producing the observed variability at early times (while the sub-MeV prompt emission is ongoing) as both invoke emission from the external forward shock. Future and more sensitive observations of such delayed emission will be important in distinguishing between the different models. (c) What mechanism produces the Fermi-LAT extended emission?: The smooth temporal decay of the GeV Fermi-LAT extended emission naturally favors its origin in emission arising from the external forward shock. The main question here is whether the emission is entirely synchrotron radiation from non-thermal shock-heated electrons [112,115,116], the standard scenario, or IC radiation from mostly quasi-thermal electrons as the blast wave encounters pair-rich and pre-accelerated circumburst medium [52,56]. The latter scenario can only operate as long as there are softer seed photons that can Compton cool the thermal electrons. At early times, they are the prompt sub-MeV photons that overlap the afterglow shock, and at later times softer photons can be of synchrotron origin [226]. Detailed numerical models of blast waves propagating into pair-enriched media and the comparison of afterglow lightcurves with observations over the entire duration of the LAT extended emission can shed more light on this issue. (d) What mechanism produces multi-GeV photons at late times?: The detection of 10 GeV photons in several GRBs at late times (t ∼ 10 2 -10 3 s) is puzzling. When their origin is interpreted as the standard afterglow synchrotron emission from shock-heated electrons, for which strong evidence came from the broad band (optical to GeV) SED fits of the afterglow of GRB 130427A [137], it challenges our understanding of particle acceleration at relativistic collisionless shocks since the photon energy clearly violates E syn,max . The alternative is IC (either SSC or EIC) afterglow emission, which can produce HE photons at late times. A prime example is GRB 190114C from which HE and VHE photons were detected by the Fermi-LAT (t 150 s) and MAGIC (t 2400 s), respectively, at late times. Future such events with multi-wavelength constraints, especially at VHEs, along with numerical simulations of particle acceleration at shock fronts will be able to shed more light on this issue.
To answer the above questions, both leptonic and hadronic models have been discussed. Important constraints on the latter scenario are offered by multi-messenger observations that include follow-up and monitoring of GRBs by neutrino detectors [227]. Even non-detections offer very useful information about the underlying radiation mechanism. However, the prospect of turning these non-detections into detections, or at least providing more stringent upper limits, is looking better with the installation of km 3 -scaled neutrino detectors, namely Baikal-GVD [228] and KM3NeT [229], in the next few to several years.

Closing Remarks
The detection of afterglow TeV γ-rays in a few GRBs, first reported for GRB 190114C, has opened up a new window for understanding the properties of relativistic collisionless shocks and radiation processes that operate near the shock fronts. VHE emission was anticipated in GRBs for some time, but it remained undetected for decades, garnering only upper limits from ground-based imaging atmospheric Cherenkov detectors. Detailed spectral modeling of the afterglow TeV emission is now shedding new light on the global energetics of the system leading to better constraints on the prompt γ-ray emission efficiency. Moreover, detection of TeV emission during the prompt-GRB phase would help pin down its illusive emission mechanism(s). In most cases, one-zone SSC emission is the most favored radiation mechanism for producing afterglow TeV photons, however, with only a few sources the details of SSC emission from shock-heated relativistic electrons (or e ± -pairs) are not entirely clear. Future, multi-messenger and perhaps more sensitive, observations from low redshift GRBs will offer better opportunities to constrain microphysical processes at shock fronts.
Author Contributions: Writing-original draft preparation-review and editing, R.G. and J.G. All authors have read and agreed to the published version of the manuscript.