© The Authors 2024

1. Introduction

The first detection of gravitational waves (GWs) from a binary neutron star (BNS) merger, GW 170817 (Abbott et al. 2017b), was followed by the detection of several electromagnetic counterparts (Abbott et al. 2017a and references therein). The short gamma-ray burst GRB 170817A, was detected ∼1.7 s after the GW signal (Goldstein et al. 2017; Savchenko et al. 2017); a fast-decaying thermal transient in the visible and infrared range, a kilonova, was observed for ∼10 days after the merger (see e.g. Villar et al. 2017; Tanvir et al. 2017; Cowperthwaite et al. 2017); and a non-thermal afterglow was observed from radio to X-rays for more than three years (see e.g. Balasubramanian et al. 2021; Troja et al. 2020; Hajela et al. 2019). Among the many advances made possible by this exceptional multi-messenger event, GW 170817/GRB 170817A is in particular the first direct association between a BNS merger and a short gamma-ray burst (GRB) and it has helped us to better understand the relativistic ejection associated with events like this. Indeed, the study of the relativistic ejecta benefited not only from an exceptional multi-wavelength follow-up, but also from observing conditions that were very different from those for other short GRBs: a much smaller distance (GW 170817 was hosted in NGC 4993 at ∼40 Mpc; Palmese et al. 2017; Cantiello et al. 2018), and a significantly off-axis observation (${32}_{-13}^{+10}\pm 1.7$ deg derived from the GW signal using the accurate localisation and distance of NGC 4993, Finstad et al. 2018).

The prompt short GRB is puzzling: it is extremely weak despite the short distance, but its peak energy is above 150 keV (Goldstein et al. 2017). It is very unlikely that it is produced by internal dissipation in the ultra-relativistic core jet like in other GRBs at cosmological distance (Matsumoto et al. 2019). GRB 170817A was rather emitted in mildly relativistic and mildly energetic material on the line of sight: The shock-breakout emission when the relativistic core jet emerges from the kilonova ejecta is a promising mechanism (Bromberg et al. 2018; Gottlieb et al. 2018).

The interpretation of the afterglow is better understood. The slow rise of the light curve until its peak after ∼120 − 160 days indicates that this was an off-axis observation of a decelerating jet surrounded by a lateral structure (see e.g. D’Avanzo et al. 2018) that may have been inherited from the early interaction with the kilonova ejecta. This was also invoked for the weak prompt emission (Gottlieb et al. 2018). The relativistic ejecta seen off-axis is also confirmed by the compactness of the source and its apparent superluminal motion as measured by very long baseline interferometry (VLBI) imagery (Mooley et al. 2018; Ghirlanda et al. 2019; Mooley et al. 2022). The light curve rise is dominated by the lateral structure of the relativistic ejecta, while its peak and decay are dominated by the deceleration of the ultra-relativistic jet core, which requires a kinetic energy comparable to values commonly found in short GRBs (see e.g. Ghirlanda et al. 2019).

A lateral structure like this may be a common feature in GRBs due to the early propagation of the relativistic ejecta, either through the infalling envelope of the stellar progenitor (long GRBs), or through the post-merger ejecta that produces the kilonova emission (short GRBs) (Bromberg et al. 2011). Possible signatures of this lateral structure in GRBs at cosmological distance viewed slightly off-axis were recently discussed, for instance to explain the complex phenomenology observed in the early afterglow (Beniamini et al. 2020a; Oganesyan et al. 2020; Ascenzi et al. 2020; Duque et al. 2022) or the non-standard decay of the afterglow of the extremely bright GRB 221009A (O’Connor et al. 2023; Gill & Granot 2023). Accounting for the lateral structure of the jet in afterglow models has thus become necessary not only for GW 170817, but for other cosmic GRBs as well.

Most models of the off-axis afterglow of GW 170817 only considered synchrotron emission and are therefore limited to the spectral range of observations, from radio to X-rays. However, upper limits on the flux at very high energy (VHE) were obtained by the high energy stereoscopic system telescopes (H.E.S.S.) at two epochs, a few days after the merger (Abdalla et al. 2017) and around its peak (Abdalla et al. 2020). Less constraining upper limits were also obtained by the high-altitude water Cherenkov gamma-ray observatory (HAWC) at early times (Galván et al. 2019) and by the major atmospheric gamma-ray imaging Cherenkov telescope (MAGIC) around the afterglow peak (Salafia et al. 2022a). A discussion of the constraints associated with these upper limits requires that the synchrotron self-Compton (SSC) emission is included in afterglow models, that is, the inverse Compton (IC) scattering of the synchrotron photons by the relativistic electrons that emit them. Recent detections of the VHE afterglow of several long GRBs¹ by H.E.S.S. (GRB 180720B; Abdalla et al. 2019; GRB 190829A; H.E.S.S. Collaboration 2021), MAGIC (GRB 190114C; MAGIC Collaboration 2019a,b; GRB 201216C; Blanch et al. 2020), and the large high altitude air shower observatory (LHAASO) (GRB 221009A; LHAASO Collaboration 2023) challenge purely synchrotron afterglow models (see however H.E.S.S. Collaboration 2021) and suggest the dominant contribution of a new emission process at VHE (see the discussion of possible processes at VHE by Gill & Granot 2022), most probably, SSC emission, as suggested for instance by the modelling of the afterglows of GRB 190114C (MAGIC Collaboration 2019b; Wang et al. 2019; Derishev & Piran 2021), GRB 190829A (Salafia et al. 2022b), or GRB 221009A (Sato et al. 2023).

While SSC emission can be efficiently computed analytically in the Thomson regime (Panaitescu & Kumar 2000; Sari & Esin 2001), several studies have pointed out that the Klein-Nishina (KN) attenuation at high energy also needs to be accounted for (see e.g. Nakar et al. 2009; Murase et al. 2011; Beniamini et al. 2015; Jacovich et al. 2021; Yamasaki & Piran 2022). In the KN regime, the IC power of an electron strongly depends on its Lorentz factor. This affects the cooling of the electron distribution, and it therefore also impacts the synchrotron spectrum (Derishev et al. 2001; Nakar et al. 2009; Bošnjak et al. 2009), which can significantly differ from the standard prediction given by Sari et al. (1998) for the purely synchrotron case.

Motivated by these recent advances in GRB afterglow studies, we propose in this paper a consistent model of GRB afterglow emission from a decelerating laterally structured jet where electrons accelerated at the external forward shock radiate at all wavelengths by synchrotron emission and SSC scattering. The treatment of SSC in the KN regime follows the approach proposed by Nakar et al. (2009), extended to account for the maximum Lorentz factor of the accelerated electrons and the attenuation at high energy due to pair creation. The numerical implementation of the model was optimised to be computationally efficient. This allowed to use Bayesian statistics to infer the parameters. We then applied this model to the multi-wavelength observations of the afterglow of GW 170817 to understand the physical constraints brought by the H.E.S.S. upper limits and to discuss whether the detection of post-merger VHE afterglows may become possible in the future, especially in the coming era of the Cherenkov telescope array (CTA; see CTA Consortium et al. 2019).

The model and its assumptions are detailed in Sect. 2 and are then tested and compared to other afterglow models in Sect. 3. In Sect. 4 we show the results obtained when fitting the afterglow of GW 170817. We discuss the predicted VHE emission in the light of the H.E.S.S. upper limit. We study the conditions for future post-merger detections of VHE afterglows in Sect. 5 and highlight that VHE emission is favoured by a higher external density compared to GW 170817, which may help to probe the population of fast-merging binaries, if it exists. Our conclusions are summarised in Sect. 6.

2. Modelling the VHE afterglow of structured jets

Our aim is to model the GRB afterglow at all wavelengths from radio bands to the TeV range. We limit ourselves to the contribution of the forward external shock propagating in the external medium and leave the extension of the model to a future work that includes the reverse-shock contribution at early times. Motivated by the observations of GW 170817, we rather focus here on two other aspects: the lateral structure of the jet, and the detailed calculation of the IC spectral component, including KN effects.

In this section, we describe our model: The assumptions for the structure of the initial relativistic outflow (Sect. 2.1), the dynamics of its deceleration (Sect. 2.2), the assumptions for the acceleration of electrons and the amplification of the magnetic field at the shock front (Sect. 2.3), the detailed calculation of the emission in the comoving frame of the shocked external medium (including synchrotron and IC components; Sect. 2.4), and finally, the integration over equal-arrival time surfaces to compute the observed flux in the observer frame (light curves and spectra; Sect. 2.5). Our numerical implementation of this model optimises the computation time (typically a few seconds for computing light curves and/or spectra at different frequencies on a laptop) so that we can explore the parameter space with a Bayesian approach for the data fitting, as described in Sect. 4. In case of ambiguity, a physical quantity q is written q without a prime in the fixed source frame, and q′ with a prime in the comoving frame of the emitting material.

2.1. Relativistic outflow: Geometry and structure

We consider a laterally structured jet. The initial energy per solid angle ϵ₀(θ) and the initial Lorentz factor Γ₀(θ) decrease with θ, which is the angle from the jet axis. We define θ_c as the opening angle of the core, that is, the ultra-relativistic/ultra-energetic central part of the jet. We then write

$$\begin{array}{c}\hfill {ϵ}_0(\theta )={ϵ}_0^{\mathrm{c}}\times f_a\left(\frac{\theta }{{\theta }_{\mathrm{c}}}\right)\end{array}$$(1)

and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_0(\theta )=1+({\mathrm{\Gamma }}_0^{\mathrm{c}}-1)\times g_b\left(\frac{\theta }{{\theta }_{\mathrm{c}}}\right),\end{array}$$(2)

where ϵ₀^c is the initial energy per solid angle, and Γ₀^c is the initial Lorentz factor, both on the jet axis (θ = 0). f_a and g_b are the normalised profiles for the energy and Lorentz factor in the lateral structure.

Different lateral structures have been suggested in the literature, following the photometric and VLBI observations of the GRB afterglow of GW 170817. In this paper, we model this afterglow assuming the power-law structure used by Duque et al. (2019),

$$\begin{array}{cc}\hfill f_a\left(x\right)& =\{\begin{array}{ccc}1& \mathrm{if}& x<1\\ x^{-a}& \mathrm{if}& x\ge 1,\end{array}\hfill \\ \hfill g_b\left(x\right)& =\{\begin{array}{ccc}1& \mathrm{if}& x<1\\ x^{-b}& \mathrm{if}& x\ge 1.\end{array}\hfill \end{array}$$(3)

This prescription allows a direct comparison with a top-hat jet (f_a(x) = g_b(x) = 1 for x ≤ 1 and 0 otherwise). For comparison with other studies (see Sect. 3), we also considered the following possible structures:

• A power-law jet from Gill & Granot (2018), f_a(x) = X^−a and g_b(x) = X^−b, with $X=\sqrt{1+x^2}$.

• A power-law jet as defined in afterglowpy, from Ryan et al. (2020), $f_a\left(x\right)={\left(\sqrt{1+x^2/a}\right)}^{-a}$. The model used in afterglowpy is limited to the self-similar evolution of the jet, which is independent of the initial value of the Lorentz factor. Therefore, g_b is not specified in this case.

• A Gaussian jet from Gill & Granot (2018), f_a(x) = g_b(x) = max(e^−x2/2;e^−x_max2/2), where x_max = θ_max/θ_c and θ_max is defined as the angle where β_0, min(θ_max) = 0.01 (by default) or any other specific value. We note that Ryan et al. (2020) used a similar parametrisation of the lateral structure for Gaussian jets in afterglowpy, with θ_max = 90 deg.

The initial total energy of the jet and of its core are given by

$$\begin{array}{cc}\hfill E_0& =2{\displaystyle {\int }_0^{\pi /2}}{ϵ}_0(\theta )2\pi sin\theta \mathrm{d}\theta ,\hfill \\ \hfill E_0^{\mathrm{c}}& =2{\displaystyle {\int }_0^{{\theta }_{\mathrm{c}}}}{ϵ}_0(\theta )2\pi sin\theta \mathrm{d}\theta .\hfill \end{array}$$(4)

The factor of 2 accounts for the counter-jet. The usual isotropic-equivalent energy of the core jet equals

$$\begin{array}{c}\hfill E_{0,\mathrm{iso}}^{\mathrm{c}}=\frac{E_0^{\mathrm{c}}}{1-cos{\theta }_{\mathrm{c}}}.\end{array}$$(5)

For the power-law structure considered here (Eq. (3)), we have $E_0^{\mathrm{c}}=4\pi {ϵ}_0^{\mathrm{c}}(1-cos{\theta }_{\mathrm{c}})$ and $E_{0,\mathrm{iso}}^{\mathrm{c}}=4\pi {ϵ}_0^{\mathrm{c}}$.

We did not consider a possible radial structure of the outflow that may also very well be present due to the variability of the central engine. This radial structure may be at least partially smoothed out during the early propagation, for example via internal shocks (Kobayashi et al. 1997; Daigne & Mochkovitch 1998, 2000), and is expected to have more impact on the reverse shock (see e.g. Uhm & Beloborodov 2007; Genet et al. 2007), which is not included in the current version of the model. Therefore, $E_{0,\mathrm{iso}}^{\mathrm{c}}$, Γ₀^c, θ_c, a, and b (or β_0,min for the Gaussian jet) are the only four free parameters needed to fully describe the initial structure of the jet before the deceleration starts.

In our numerical implementation of the model, we suppressed the lateral structure above the maximum angle θ_max. This angle is taken as the maximum between θ_ϵ, defined as the angle up to which a fraction (1 − ϵ) of the jet energy is contained, and the viewing angle θ_v (see Sect. 2.5). We therefore solved for θ_ϵ the equation

$$\begin{array}{c}\hfill (1-ϵ)E_0=2{\displaystyle {\int }_0^{{\theta }_{ϵ}}}{ϵ}_0(\theta )2\pi sin\theta \mathrm{d}\theta .\end{array}$$(6)

We use ϵ = 0.01 in the following. Ryan et al. (2020) introduced a similar parameter θ_W to minimise computation time. Taking θ_max = max(θ_ϵ;θ_v) allowed us to keep a precise calculation at early times even for very large viewing angles.

In practice, the lateral structure is discretised into N + 1 components, i = 0 being the core jet and i = 1 → N being rings at increasing angles in the lateral structure. Each component is defined by θ_min, i ≤ θ ≤ θ_max, i with θ_min, 0 = 0 and θ_max, 0 = θ_c for the core jet; θ_min, i = θ_{max, i − 1} for i = 1 → N, and θ_max, N = θ_max. Each component is treated independently with fixed limits in latitude (no lateral spreading). The dynamics is computed such that R_i(t) = R(θ_i; t) and Γ_i(t) = Γ(θ_i; t), where R(θ, t) and Γ(θ, t) are given by the solution for the dynamics of the deceleration discussed below, and where θ_i = (θ_min, i + θ_max, i)/2 for i ≥ 1 and θ₀ = 0. We define the successive θ_min, i and θ_max, i such that each component carries an equal amount of energy (except for the core, i = 0), although a linear increase can also be chosen. Finally, we typically choose a discretisation comprised of N = 15 components, which we find to be a good compromise between model accuracy and computational efficiency.

2.2. Dynamics of the deceleration

The structured outflow decelerates in an external medium, with an assumed density profile as a function of radius R,

$$\begin{array}{c}\hfill {\rho }_{\mathrm{ext}}(R)=\frac{A}{R^s},\end{array}$$(7)

with s = 0 for a constant-density external medium as considered in the following to model the afterglow of GW 170817. In this case, we write A = n_extm_p, with n_ext the external medium density. We focus on the dynamics of the shocked external medium at the forward shock. It is computed assuming that (i) the dynamics of each ring of material at angle θ is independent of other angles, and that (ii) the lateral expansion of the outflow is negligible. These two assumptions are questionable at late times, close to the transition to the Newtonian regime (Rhoads 1997).

The Lorentz factor Γ(θ; R) of the shocked material at angle θ and radius R is computed in a simplified way, using energy conservation (see e.g. Panaitescu & Kumar 2000; Gill & Granot 2018),

$$\begin{array}{ccc}\hfill \mathrm{\Gamma }(\theta ;R)& =\hfill & \hfill \frac{{\mathrm{\Gamma }}_0(\theta )}{2m(\theta ;R)}\times \\ \hfill & [-1+\sqrt{1+4m(\theta ;R)+4{\left(\frac{m(\theta ;R)}{{\mathrm{\Gamma }}_0(\theta )}\right)}^2}],\hfill \end{array}$$(8)

where

$$\begin{array}{c}\hfill m(\theta ;R)=\frac{M_{\mathrm{ext}}(R)}{M_{\mathrm{ej}}(\theta )/{\mathrm{\Gamma }}_0(\theta )}={\left(\frac{R}{R_{\mathrm{dec}}(\theta )}\right)}^{3-s},\end{array}$$(9)

with M_ext(R) = ∫₀^Rρ_ext(R) R²dR the swept-up mass per unit solid angle, and M_ej(θ) = ϵ₀(θ)/Γ₀(θ)c² the ejected mass per unit solid angle. The deceleration radius of the material at angle θ is defined by

$$\begin{array}{c}\hfill R_{\mathrm{dec}}(\theta )={\left(\frac{(3-s){ϵ}_0(\theta )}{A{\mathrm{\Gamma }}_0^2(\theta )c^2}\right)}^{1/(3-s)}.\end{array}$$(10)

For a power-law structure as defined by Eqs. (3) and (4), the deceleration radius scales as $R_{\mathrm{dec}}(\theta )\propto {\theta }^{\frac{2b-a}{3-s}}$. Best-fit models of the afterglow of GW 170817 usually have 2b > a so that the deceleration radius is larger at high latitude (see Sect. 4.2). More generally, Beniamini et al. (2020b) derived the conditions for the observed afterglow light curve to have a single peak, as in the case of GW 170817. For an off-axis observation with a large viewing angle, these conditions are Γ₀^cθ_c > 1; 2b > a/(4 − s) (i.e. 8b > a for a uniform external medium), b > −lnΓ₀^c/lnθ_c, and ${\theta }_{\mathrm{v}}>{\theta }_{\mathrm{c}}{\left({\mathrm{\Gamma }}_0^{\mathrm{c}}{\theta }_{\mathrm{c}}\right)}^{\frac{1}{b-1}}$.

This description allows us to characterise the early-time dynamics in the coasting phase R ≪ R_dec(θ) where the Lorentz factor remains constant. It continuously branches out to the relativistic self-similar evolution (Blandford & McKee 1976) for $R_{\mathrm{dec}}(\theta )\ll R\ll R_{\mathrm{N}}(\theta )={\mathrm{\Gamma }}_0^{\frac{2}{3-s}}(\theta )R_{\mathrm{dec}}(\theta )$, and finally, to the Sedov-Taylor phase (Sedov 1946; Taylor 1950) in the non-relativistic regime for R ≫ R_N(θ). For R ≫ R_dec(θ), the self-similar evolution becomes independent of the value of the initial Lorentz factor Γ₀(θ).

When the Lorentz factor Γ(θ; R) and the velocity $\beta (\theta ;R)=\frac{\sqrt{{\mathrm{\Gamma }}^2(\theta ;R)-1}}{\mathrm{\Gamma }(\theta ;R)}$ are known, the corresponding time t in the source frame is given by

$$\begin{array}{c}\hfill t(\theta ;R)={\displaystyle {\int }_0^R}\frac{\mathrm{d}r}{\beta (\theta ;r)c}.\end{array}$$(11)

This leads to the solution R(θ, t) and Γ(θ, t) that we use for the dynamics of each component of the structured jet.

The physical conditions in the shocked medium are easily deduced from the shock-jump conditions in the strong-shock regime (Blandford & McKee 1976). In the comoving frame, the mass density equals

$$\begin{array}{c}\hfill {\rho }_{\ast }(\theta ;R)=(4\mathrm{\Gamma }(\theta ;R)+3){\rho }_{\mathrm{ext}}(R)\end{array}$$(12)

and the internal energy per unit mass is given by

$$\begin{array}{c}\hfill {ϵ}_{\ast }(\theta ;R)=(\mathrm{\Gamma }(\theta ;R)-1)c^2.\end{array}$$(13)

We assumed an adiabatic index of 4/3, which is valid as long as ϵ_* ≫ c². Finally, the timescale of the adiabatic cooling of the shocked region due to the spherical expansion is given in the comoving frame by

$$\begin{array}{c}\hfill t_{\mathrm{dyn}}^{\prime }(\theta ;R)=\frac{R}{\mathrm{\Gamma }(\theta ;R)\beta (\theta ;R)c}=\frac{R}{c\sqrt{{\mathrm{\Gamma }}^2(\theta ;R)-1}}.\end{array}$$(14)

2.3. Accelerated electrons and amplified magnetic field

We consider a shocked region where the physical conditions in the comoving frame are given by the mass density ρ_*, the internal energy per unit mass ϵ_*, and the dynamical timescale t′_dyn (i.e. the characteristic timescale of the adiabatic cooling due to the spherical expansion), and we assume that the emission is produced by non-thermal shock-accelerated electrons that radiate in a local turbulent magnetic field amplified at the shock. In practice, we only consider the forward external shock, so that in this study, ρ_*, ϵ_*, and t′_dyn are given by Eqs. (14)–(16).

We use the following standard parametrisation of the microphysics at the shock: (i) a fraction ϵ_B of the internal energy is injected in the magnetic field,

$$\begin{array}{c}\hfill u_{\mathrm{B}}=\frac{{B^{\prime }}^2}{8\pi }={ϵ}_{\mathrm{B}}{\rho }_{\ast }{ϵ}_{\ast },\end{array}$$(15)

leading to

$$\begin{array}{c}\hfill B^{\prime }=\sqrt{8\pi {ϵ}_{\mathrm{B}}{\rho }_{\ast }{ϵ}_{\ast }};\end{array}$$(16)

(ii) a fraction ϵ_e of the internal energy is injected into non-thermal electrons that represent a fraction ζ of all available electrons. Their number density (cm⁻³) and energy density (erg ⋅ cm⁻³) in the comoving frame are therefore given by

$$\begin{array}{c}\hfill n_{\mathrm{e}}^{\mathrm{acc}}=\zeta \frac{{\rho }_{\ast }}{m_{\mathrm{p}}}\end{array}$$(17)

and

$$\begin{array}{c}\hfill u_{\mathrm{e}}={ϵ}_{\mathrm{e}}{\rho }_{\ast }{ϵ}_{\ast }.\end{array}$$(18)

We assume a power-law distribution at injection with an index −p and 2 < p < 3. This leads to the following distribution of accelerated electrons (cm⁻³):

$$\begin{array}{c}\hfill n(\gamma )=(p-1)\frac{n_{\mathrm{e}}^{\mathrm{acc}}}{{\gamma }_{\mathrm{m}}}{\left(\frac{\gamma }{{\gamma }_{\mathrm{m}}}\right)}^{-p}\mathrm{for}{\gamma }_{\mathrm{m}}\le \gamma \le {\gamma }_{\max },\end{array}$$(19)

where the minimum Lorentz factor at injection equals

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{m}}=\frac{p-2}{p-1}\frac{{ϵ}_{\mathrm{e}}}{\zeta }\frac{m_{\mathrm{p}}}{m_{\mathrm{e}}}\frac{{ϵ}_{\ast }}{c^2}.\end{array}$$(20)

As discussed below, it is important to take a realistic estimate of the maximum Lorentz factor γ_max up to which electrons can be accelerated at the shock into account for a discussion of the GeV-TeV afterglow emission. This was not included in the study by Nakar et al. (2009), on which we base our calculation of the emission. Eqs. (21) and (22) above were obtained by assuming that the maximum Lorentz factor γ_max is much higher than γ_m. This is fully justified as we have typically γ_max/γ_m > 10⁶ in the best-fit models of GW 170817 at the peak (Sect. 4). The maximum electron Lorentz factor γ_max is evaluated by imposing that the acceleration timescale always remains shorter than the radiative and dynamical timescales,

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{\prime }(\gamma )\le min(t_{\mathrm{rad}}^{\prime }(\gamma );t_{\mathrm{dyn}}^{\prime }).\end{array}$$(21)

The acceleration timescale is written as a function of the Larmor time R′_L/c following Bohm’s scaling,

$$\begin{array}{c}\hfill t_{\mathrm{acc}}^{\prime }(\gamma )=K_{\mathrm{acc}}\frac{R_{\mathrm{L}}^{\prime }(\gamma )}{c}=K_{\mathrm{acc}}\frac{\gamma m_{\mathrm{e}}c}{e{B}^{\prime }},\end{array}$$(22)

where K_acc ≥ 1 is a dimensionless factor. In practice, electrons at γ_max usually cool fast, and the maximum electron Lorentz factor is determined by the radiative timescale, t′_rad(γ). Its evaluation is non-trivial when IC scattering is taken into account (see Sect. 2.4). The resulting detailed calculation of γ_max is explained in Appendix A. In the following, we assume K_acc = 1. Our value of γ_max should therefore be considered as an upper limit for the true value of the maximum electron Lorentz factor.

2.4. Emissivity in the comoving frame

Our calculation of the emission from non-thermal electrons is mostly taken from Nakar et al. (2009). The population of electrons cools down by adiabatic cooling, synchrotron radiation, and via SSC scattering. As we are interested in the VHE afterglow emission, we include the attenuation due to pair production in Sect. 2.4.8, but do not include the emission of secondary leptons. We focus on the afterglow of GW 170817, for which the radio observations did not show any evidence for absorption, and we therefore do not include the effect of synchrotron self-absorption in the current version of the model. In this subsection, all quantities are written in the comoving frame of the considered shocked region. We therefore omit the prime to simplify the notations.

2.4.1. Synchrotron and inverse-Compton radiation for a single electron

2.4.1.1. Synchrotron radiation.

The synchrotron power of a single electron (erg ⋅ s⁻¹) is given by

$$\begin{array}{c}\hfill P^{\mathrm{syn}}(\gamma )=K_1{B}^2{\gamma }^2,\end{array}$$(23)

with $K_1=K_P\frac{{\sigma }_{\mathrm{T}}c}{6\pi }$ and K_P a dimensionless parameter, and its synchrotron frequency equals

$$\begin{array}{c}\hfill {\nu }_{\mathrm{syn}}(\gamma )=K_{2}B{\gamma }^2,\end{array}$$(24)

with $K_2=K_{\nu }\frac{e}{2\pi m_{\mathrm{e}}c}$ and K_ν a dimensionless parameter. The dimensionless parameters K_P and K_ν depend on the assumptions made on the pitch angle of electrons relative to the magnetic field lines and different values are found in afterglow models, as listed in Table 1. By default, we use the same values as Gill & Granot (2018), taken from Granot et al. (1999), who obtained them by fitting the broken power-law approximation of the synchrotron power of a power-law distribution of electrons to the exact calculation.

The corresponding synchrotron power at frequency ν (erg ⋅ s⁻¹ ⋅ Hz⁻¹) is given by

$$\begin{array}{c}\hfill P_{\nu }^{\mathrm{syn}}(\gamma )=P_{\max }\mathrm{\Phi }\left(\frac{\nu }{{\nu }_{\mathrm{syn}}(\gamma )}\right),\end{array}$$(25)

with

$$\begin{array}{c}\hfill P_{\max }=\frac{P^{\mathrm{syn}}(\gamma )}{{\nu }_{\mathrm{syn}}(\gamma )}=\frac{K_1}{K_2}B=K_{P_{\max }}\frac{{\sigma }_{\mathrm{T}}{m}_{\mathrm{e}}{c}^2}{3e}B,\end{array}$$(26)

K_Pmax = K_P/K_ν and ∫₀^∞Φ(x)dx = 1. In practice, we adopt a simplified shape $\mathrm{\Phi }(x)=\frac{4}{3}{x}^{1/3}$ for x ≤ 1 and 0 otherwise.

2.4.1.2. Inverse-Compton scattering.

The spectral number density of seed photons is n_ν^syn (cm⁻³ ⋅ Hz⁻¹), which is assumed to result only from the synchrotron radiation of the electrons on which they also scatter (SSC). The IC power of an electron (erg ⋅ s⁻¹) is then given by

$$\begin{array}{c}\hfill P^{\mathrm{SSC}}(\gamma )={\displaystyle {\int }_0^{\infty }}\mathrm{d}\nu n_{\nu }^{\mathrm{syn}}{\sigma }_{\mathrm{T}}{f}_{\mathrm{KN}}(w)c{K}_{\mathrm{IC}}{\gamma }^2{g}_{\mathrm{KN}}(w)h\nu ,\end{array}$$(27)

where w = γhν/m_ec²; and f_KN and g_KN are the KN corrections to the cross-section and to the mean energy of the scattered photon. In the Thomson regime (w ≪ w_KN), the cross-section is σ_T and the mean energy of scattered photons is K_IC γ²hν for a seed synchrotron photon with a frequency ν, so that f_KN = g_KN = 1 for w ≪ w_KN. By default, we use the standard values K_IC = 4/3 and w_KN = 1, as in Nakar et al. (2009). We note, however, that Yamasaki & Piran (2022) recently suggested that a better agreement with a detailed numerical spectral calculation of the SSC component was achieved with w_KN = 0.2, especially at the transition between the Thomson and the KN regime. For simplicity, we do not write in Eq. (29) the dependence on the angle between the upscattered photon and the scattering electron explicitly. This was discussed in Nakar et al. (2009). We assume an isotropic seed radiation field. We only account for single IC scattering as the effects of multiple scatterings are expected to be negligible because the afterglow is produced in the optically thin regime. The Compton parameter is defined by

$$\begin{array}{c}\hfill Y(\gamma )=\frac{P^{\mathrm{SSC}}(\gamma )}{P^{\mathrm{syn}}(\gamma )}.\end{array}$$(28)

This Compton parameter is constant only when SSC scattering occurs in the Thomson regime. We take the KN regime into account, which explains the dependence on the electron Lorentz factor. This affects not only the SSC emission, but also the synchrotron component. Following Nakar et al. (2009) and keeping the same notations, we simplify the IC cross-section by assuming that scattering is entirely suppressed in the KN regime: f_KN(w) = 0 for w ≥ w_KN. Then, a synchrotron photon produced by an electron with Lorentz factor γ can only be upscattered by electrons with Lorentz factors below a certain limit $\widehat{\gamma }$ to remain in the Thomson regime,

$$\begin{array}{c}\hfill \widehat{\gamma }=\frac{w_{\mathrm{KN}}{m}_{\mathrm{e}}{c}^2}{h{\nu }_{\mathrm{syn}}(\gamma )}\propto {\gamma }^{-2}.\end{array}$$(29)

We note that $\widehat{\gamma }$ is a decreasing function of γ: If the energy of the synchrotron photon is higher, the maximum energy of any electron on which it could be upscattered in the Thomson regime must be lower. In an equivalent way, electrons with Lorentz factor γ only scatter photons below the energy

$$\begin{array}{c}\hfill h\tilde{\nu }=w_{\mathrm{KN}}{m}_{\mathrm{e}}{c}^2/\gamma .\end{array}$$(30)

Then, it is convenient to define $\stackrel{\sim }{\gamma }$ such that ${\nu }_{\mathrm{syn}}(\stackrel{\sim }{\gamma })=\stackrel{\sim }{\nu }$,

$$\begin{array}{c}\hfill \tilde{\gamma }=\sqrt{w_{\mathrm{KN}}\gamma m_{\mathrm{e}}{c}^2/h{\nu }_{\mathrm{syn}}(\gamma )}\propto {\gamma }^{-1/2}.\end{array}$$(31)

Electrons with a Lorentz factor γ can only scatter synchrotron photons produced by electrons with Lorentz factors below $\stackrel{\sim }{\gamma }$. The Lorentz factor γ_self is defined by

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{self}}={\widehat{\gamma }}_{\mathrm{self}}={\tilde{\gamma }}_{\mathrm{self}}={\left(\frac{w_{\mathrm{KN}}{m}_{\mathrm{e}}{c}^2}{h{K}_{2}B}\right)}^{1/3}.\end{array}$$(32)

The definition of $\stackrel{\sim }{\nu }$ allows us to simplify the expression of the SSC power given by Eq. (29):

$$\begin{array}{c}\hfill P^{\mathrm{SSC}}(\gamma )=K_{\mathrm{IC}}{\sigma }_{\mathrm{T}}c{\gamma }^2{\displaystyle {\int }_0^{\tilde{\nu }}}\mathrm{d}\nu u_{\nu }^{\mathrm{syn}},\end{array}$$(33)

with u_ν^syn = n_ν^syn × hν the energy density of the seed synchrotron photons (erg ⋅ cm⁻³ ⋅ Hz⁻¹). Then,

$$\begin{array}{c}\hfill Y(\gamma )=\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{{\displaystyle {\int }_0^{\tilde{\nu }}}\mathrm{d}\nu u_{\nu }^{\mathrm{syn}}}{u_{\mathrm{B}}}.\end{array}$$(34)

Assuming a strict suppression in the KN regime leads to a simple expression for the SSC power at frequency ν (erg ⋅ s⁻¹ ⋅ Hz⁻¹) for a single electron with a Lorentz factor γ,

$$\begin{array}{c}\hfill P_{\nu }^{\mathrm{SSC}}(\gamma )={\sigma }_{\mathrm{T}}c{u}_{{\nu }_{\mathrm{seed}}=\nu /\left(K_{\mathrm{IC}}{\gamma }^2\right)}^{\mathrm{syn}}\end{array}$$(35)

if hν ≤ K_ICγm_ec² and 0 otherwise.

2.4.2. Radiative regime of a single electron

Following Sari et al. (1998), the first expected break in the electron distribution occurs at the critical electron Lorentz factor γ_c, defined by the condition t_rad(γ) = t_dyn, so that high-energy electrons with γ ≫ γ_c are radiatively efficient (fast cooling), and low-energy electrons with γ ≪ γ_c mainly cool via the adiabatic expansion and are radiatively inefficient (slow cooling). The radiative timescale is given by

$$\begin{array}{c}\hfill t_{\mathrm{rad}}(\gamma )=\frac{\gamma m_{\mathrm{e}}{c}^2}{P_{\mathrm{syn}}(\gamma )+P_{\mathrm{SSC}}(\gamma )}=\frac{m_{\mathrm{e}}{c}^2}{[1+Y(\gamma )]K_1{B}^2\gamma }.\end{array}$$(36)

Then, self-consistently computing γ_c requires us to solve the equation

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{c}}[1+Y({\gamma }_{\mathrm{c}})]={\gamma }_{\mathrm{c}}^{\mathrm{syn}}=\frac{m_{\mathrm{e}}{c}^2}{K_1{B}^2{t}_{\mathrm{dyn}}},\end{array}$$(37)

where ${\gamma }_{\mathrm{c}}^{\mathrm{syn}}$ is the value of the critical Lorentz factor obtained when only the synchrotron radiation is taken into account. Our procedure to compute γ_c and Y(γ_c) in the most general case where the SSC cooling impacts the electron distribution is discussed below in Sect. 2.4.6.

2.4.3. Electron distribution and associated emission

2.4.3.1. Electron distribution.

Following Nakar et al. (2009), we compute the time-averaged distribution of electrons $\overline{n}(\gamma )$ over the dynamical timescale t_dyn by keeping only the dominant term for the electron cooling: The instantaneous power of an electron is either dominated by the synchrotron and SSC radiation (1+Y(γ))P^syn(γ) for γ > γ_c or by the adiabatic cooling γ m_ec²/t_dyn if γ < γ_c. This leads to

$$\begin{array}{c}\hfill \overline{n}(\gamma )=\frac{n_{\mathrm{e}}^{\mathrm{acc}}}{1+Y(\gamma )}\frac{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}}{{\gamma }_{\mathrm{m}}^2}\times \{\begin{array}{cc}{(\gamma /{\gamma }_{\mathrm{m}})}^{-2}& \mathrm{if}{\gamma }_{\mathrm{c}}<\gamma <{\gamma }_{\mathrm{m}}\hfill \\ {(\gamma /{\gamma }_{\mathrm{m}})}^{-(p+1)}& \mathrm{if}{\gamma }_{\mathrm{m}}<\gamma <{\gamma }_{\max }\hfill \end{array}\end{array}$$(38)

in fast-cooling regime (γ_m > γ_c) and to

$$\begin{array}{ccc}\hfill \overline{n}(\gamma )& =\hfill & \hfill n_{\mathrm{e}}^{\mathrm{acc}}\frac{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}}{{\gamma }_{\mathrm{m}}^2}{\left(\frac{{\gamma }_{\mathrm{c}}}{{\gamma }_{\mathrm{m}}}\right)}^{-(p+1)}\\ \hfill & \times \{\begin{array}{cc}{(\gamma /{\gamma }_{\mathrm{c}})}^{-p}/(1+Y({\gamma }_{\mathrm{c}}))& \mathrm{if}{\gamma }_{\mathrm{m}}<\gamma <{\gamma }_{\mathrm{c}}\hfill \\ {(\gamma /{\gamma }_{\mathrm{c}})}^{-(p+1)}/(1+Y(\gamma ))& \mathrm{if}{\gamma }_{\mathrm{c}}<\gamma <{\gamma }_{\max }\hfill \end{array}\hfill \\ \hfill \end{array}$$(39)

in slow-cooling regime (γ_m < γ_c). We define γ_m, c = max(γ_m;γ_c) and γ_c, m = min(γ_m;γ_c) so that γ_m, c = γ_m and γ_c, m = γ_c in fast cooling and γ_m, c = γ_c and γ_c, m = γ_m in slow cooling. The minimum electron Lorentz factor therefore is γ_min = γ_c, m. We define the corresponding synchrotron frequencies ν_m, c = ν_syn(γ_m, c) and ν_c, m = ν_syn(γ_c, m). We also always only keep the dominant term in the radiated power, that is, 1 + Y(γ)≃1 if Y(γ) < 1 and Y(γ) otherwise. As described in detail in Nakar et al. (2009), the resulting distribution shows several breaks in addition to γ_m and γ_c: A break is expected at γ₀ such that Y(γ₀) = 1, and several additional breaks are possible and must be identified in an iterative way, as described in Sect. 2.4.6.

Following Sari et al. (1998) and Nakar et al. (2009), we approximate $\overline{n}(\gamma )$ and the associated synchrotron spectrum by broken power laws to allow for a semi-analytical calculation. Electrons Lorentz factors are normalised by γ_m, c, that is, x = γ/γ_m, c, and photon frequencies by ν_m, c, that is, y = ν/ν_m, c. Hence, the normalised synchrotron frequency of an electron with normalised Lorentz factor x is y = x². In the following, the notations introduced for characteristic electron Lorentz factors or photon frequencies are implicitly conserved for the normalised quantities. For instance, x_self = γ_self/γ_m, c, $\widehat{x}=\widehat{\gamma }/{\gamma }_{\mathrm{m},\mathrm{c}}=x_{\mathrm{self}}^3/x^2$, and $\stackrel{\sim }{y}=\stackrel{\sim }{\nu }/{\nu }_{\mathrm{m},\mathrm{c}}={\stackrel{\sim }{x}}^2$.

2.4.3.2. Normalised electron distribution.

The electron distribution is given by

$$\begin{array}{c}\hfill \overline{n}(\gamma )=\frac{1}{I_0}\frac{n_{\mathrm{e}}^{\mathrm{acc}}}{{\gamma }_{\mathrm{m},\mathrm{c}}}f(x),\end{array}$$(40)

where f(x) is the normalised broken power-law electron distribution with N breaks and power-law segments,

$$\begin{array}{c}\hfill f(x)=\{\begin{array}{cc}\hfill x^{-p_1}& \mathrm{if}{x}_1\le x\le x_2\hfill \\ \hfill x_2^{p_2-p_1}{x}^{-p_2}& \mathrm{if}{x}_2\le x\le x_3\hfill \\ \hfill x_2^{p_2-p_1}{x}_3^{p_3-p_2}{x}^{-p_3}& \mathrm{if}{x}_2\le x\le x_3\hfill \\ \hfill ⋮& ⋮\hfill \\ \hfill x_2^{p_2-p_1}{x}_3^{p_3-p_2}\cdots x_N^{p_N-p_{N-1}}{x}^{-p_N}& \mathrm{if}{x}_N\le x\le x_{N+1}\hfill \end{array}\end{array}$$(41)

with x₁ = γ_min/γ_m, c and x_N + 1 = γ_max/γ_m, c. The number of breaks in the purely synchrotron case or when all IC scatterings are in the Thomson regime is N = 2, as discussed in Sect. 2.4.4 and 2.4.5. In the general case, where KN corrections cannot be neglected, more breaks appear in the electron distribution, as discussed in Sect. 2.4.6. All possible cases are provided in Appendix B. Typically, N = 2 to 5, except in some very specific regions of the parameter space (case III in Nakar et al. 2009). This is also discussed in Appendix B.

The dimensionless integral I₀ is given by

$$\begin{array}{c}\hfill I_0={\displaystyle {\int }_{x_1}^{x_{N+1}}}f(x)\mathrm{d}x\end{array}$$(42)

to conserve the number of electrons, ${\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}}\overline{n}(\gamma )\mathrm{d}\gamma =n_{\mathrm{e}}^{\mathrm{acc}}$.

2.4.3.3. Normalised synchrotron spectrum.

The synchrotron power per electron and per unit frequency p_ν^syn (erg ⋅ s⁻¹ ⋅ Hz⁻¹ ⋅ electron⁻¹) averaged over the timescale t_dyn is deduced from Eq. (27) using the simplified shape for Φ. This leads to a broken power-law synchrotron spectrum,

$$\begin{array}{cc}\hfill p_{\nu }^{\mathrm{syn}}& =\frac{1}{n_{\mathrm{e}}^{\mathrm{acc}}}{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}}\mathrm{d}\gamma \overline{n}(\gamma )P_{\nu }^{\mathrm{syn}}(\gamma )\hfill \end{array}$$(43)

$$\begin{array}{cc}& =\frac{{\gamma }_{\mathrm{m},\mathrm{c}}^2{m}_{\mathrm{e}}{c}^2}{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}{\nu }_{\mathrm{m},\mathrm{c}}{t}_{\mathrm{dyn}}}\times \frac{4}{3}{y}^{1/3}{\displaystyle {\int }_{max(x_1;y^{1/2})}^{x_{N+1}}}\mathrm{d}x\frac{f(x)}{x^{2/3}}\hfill \\ \hfill & \simeq \frac{1}{I_0}\frac{I_2}{J_0}\frac{{\gamma }_{\mathrm{m},\mathrm{c}}^2{m}_{\mathrm{e}}{c}^2}{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}{\nu }_{\mathrm{m},\mathrm{c}}{t}_{\mathrm{dyn}}}\times g(y),\hfill \end{array}$$(44)

where the normalised broken power-law spectral shape is obtained by keeping the dominant term in the integral of f(x)/x^2/3,

$$\begin{array}{c}\hfill g(y)=\{\begin{array}{cc}\hfill x_1^{\frac{1}{3}-p_1}{y}^{\frac{1}{3}}& \mathrm{if}y<x_1^2\hfill \\ \hfill y^{-\frac{p_1-1}{2}}& \mathrm{if}{x}_1^2<y<x_2^2\hfill \\ \hfill x_2^{p_2-p_1}{y}^{-\frac{p_2-1}{2}}& \mathrm{if}{x}_2^2<y<x_3^2\hfill \\ \hfill x_2^{p_2-p_1}{x}_3^{p_3-p_2}{y}^{-\frac{p_3-1}{2}}& \mathrm{if}{x}_3^2<y<x_4^2\hfill \\ \hfill ⋮& ⋮\hfill \\ \hfill x_2^{p_2-p_1}{x}_3^{p_3-p_2}\cdots x_N^{p_N-p_{N-1}}{y}^{-\frac{p_N-1}{2}}& \mathrm{if}{x}_N^2<y<x_{N+1}^2\hfill \end{array}.\end{array}$$(45)

The dimensionless factors I₂ and J₀ are defined by

$$\begin{array}{cc}\hfill I_2& ={\displaystyle {\int }_{x_1}^{x_{N+1}}}\mathrm{d}x{x}^{2}f(x)\hfill \end{array}$$(46)

$$\begin{array}{cc}\hfill J_0& ={\displaystyle {\int }_0^{x_{N+1}^2}}\mathrm{d}yg(y)\hfill \end{array}$$(47)

to ensure that the total synchrotron power per electron (erg ⋅ s⁻¹ ⋅ electron⁻¹) is conserved,

$$\begin{array}{cc}\hfill p^{\mathrm{syn}}& ={\displaystyle {\int }_0^{\infty }}\mathrm{d}\nu p_{\nu }^{\mathrm{syn}}=\frac{I_2}{I_0}\frac{{\gamma }_{\mathrm{m},\mathrm{c}}^2{m}_{\mathrm{e}}{c}^2}{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}{t}_{\mathrm{dyn}}}\hfill \\ \hfill & =\frac{1}{n_{\mathrm{e}}^{\mathrm{acc}}}{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}}\mathrm{d}\gamma \overline{n}(\gamma )P^{\mathrm{syn}}(\gamma ).\hfill \end{array}$$(48)

Finally, we define ν_p as the peak frequency of the synchrotron spectrum, that is, the frequency ν where ν²p_ν^syn is maximum. This peak frequency corresponds to synchrotron photons emitted by electrons at a Lorentz factor γ_p, that is, ν_p = ν_syn(γ_p). We note i_p the corresponding index of the break in the normalised distributions: x_p = x_ip and y_p = y_ip = x_p².

2.4.3.4. Normalised SSC spectrum.

The SSC power per electron and per unit frequency p_ν^SSC (erg ⋅ s⁻¹ ⋅ Hz⁻¹ ⋅ electron⁻¹) emitted over the timescale t_dyn is deduced from Eq. (37) with

$$\begin{array}{c}\hfill u_{\nu }^{\mathrm{syn}}=n_{\mathrm{e}}^{\mathrm{acc}}{p}_{\nu }^{\mathrm{syn}}{t}_{\mathrm{dyn}}.\end{array}$$(49)

This leads to

$$\begin{array}{cc}\hfill p_{\nu }^{\mathrm{SSC}}& =\frac{1}{n_{\mathrm{e}}^{\mathrm{acc}}}{\displaystyle {\int }_{{\gamma }_{\min }}^{{\gamma }_{\max }}}\mathrm{d}\gamma \overline{n}(\gamma )P_{\nu }^{\mathrm{SSC}}(\gamma )\hfill \\ \hfill & \simeq \hfill & \hfill \frac{1}{I_0}\frac{I_2}{J_0}\frac{{\tau }_{\mathrm{T}}}{I_0}\frac{{\gamma }_{\mathrm{m},\mathrm{c}}^2{m}_{\mathrm{e}}{c}^2}{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}{\nu }_{\mathrm{m},\mathrm{c}}{t}_{\mathrm{dyn}}}\times G(y),\end{array}$$(50)

where

$$\begin{array}{c}\hfill {\tau }_{\mathrm{T}}=n_{\mathrm{e}}^{\mathrm{acc}}{\sigma }_{\mathrm{T}}c{t}_{\mathrm{dyn}}\end{array}$$(51)

is the Thomson optical depth, and where the normalised spectral shape is given by

$$\begin{array}{c}\hfill G(y)={\displaystyle {\int }_{max(x_1;\frac{1}{K_{\mathrm{IC}}}\frac{h{\nu }_{\mathrm{m},\mathrm{c}}}{{\gamma }_{\mathrm{m},\mathrm{c}}{m}_{\mathrm{e}}{c}^2}y)}^{x_{N+1}}}\mathrm{d}xf(x)g(y_{\mathrm{seed}}=\frac{1}{K_{\mathrm{IC}}}\frac{y}{{\gamma }_{\mathrm{m},\mathrm{c}}^2{x}^2}).\end{array}$$(52)

The cutoff at x_KN = γ_KN/γ_m, c with ${\gamma }_{\mathrm{KN}}=\frac{1}{K_{\mathrm{IC}}}\frac{h\nu }{m_{\mathrm{e}}{c}^2}$ is a direct consequence of the strict suppression assumed in the KN regime: The VHE flux due to the scattering in the KN regime is neglected. Yamasaki & Piran (2022) introduced a factor f_KN in their calculations to compensate for this approximation, which we do not include here. The KN regime affects the high-energy part of the emitted spectrum: Above a photon energy hν = K_ICγ_minm_ec², we have γ_KN > γ_min and the SSC emission is reduced. The SSC emission is even entirely suppressed (G(y) = 0) at very high photon energies hν > K_ICγ_maxm_ec², corresponding to γ_KN > γ_max.

In practice, G(y) is computed exactly in our numerical implementation. We also note that J₀ ≃ 2I₂ when only the leading term is kept in the integrals. However, all dimensionless factors I₀, I₂, and J₀ are also computed exactly in our numerical implementation, which ensures the continuity of the electron distribution and emitted spectrum at the transition from the fast- to the slow-cooling regime, as well as between the different possible cases in either regime (see Sect. 2.4.6).

2.4.4. Purely synchrotron case

The standard purely synchrotron case in which the SSC emission is neglected (Y(γ) = 0 for all electrons) was fully described in Sari et al. (1998). We have γ_c = γ_c^syn in this case. In the fast-cooling regime (γ_m > γ_c), the normalised distributions are given by

$$\begin{array}{cc}\hfill f(x)& =\{\begin{array}{cc}\hfill x^{-2}& \mathrm{if}{x}_1<x<x_2\hfill \\ \hfill x_2^{p-1}{x}^{-(p+1)}& \mathrm{if}{x}_2<x<x_3\hfill \end{array},\hfill \end{array}$$(53)

$$\begin{array}{cc}\hfill g(y)& =\{\begin{array}{cc}\hfill x_1^{-5/3}{y}^{1/3}& \mathrm{if}y<x_1^2\hfill \\ \hfill y^{-1/2}& \mathrm{if}{x}_1^2<y<x_2^2\hfill \\ \hfill x_2^{p-1}{y}^{-p/2}& \mathrm{if}{x}_2^2<y<x_3^2\hfill \end{array},\hfill \end{array}$$(54)

with x₁ = γ_c/γ_m, x₂ = 1 and x₃ = γ_max/γ_m. In the slow-cooling regime (γ_m < γ_c), they are given by

$$\begin{array}{cc}\hfill f(x)& =\{\begin{array}{cc}\hfill x^{-p}& \mathrm{if}{x}_1<x<x_2\hfill \\ \hfill x_2{x}^{-(p+1)}& \mathrm{if}{x}_2<x<x_3\hfill \end{array},\hfill \end{array}$$(55)

$$\begin{array}{cc}\hfill g(y)& =\{\begin{array}{cc}\hfill x_1^{1/3-p}{y}^{1/3}& \mathrm{if}y<x_1^2\hfill \\ \hfill y^{-(p-1)/2}& \mathrm{if}{x}_1^2<y<x_2^2\hfill \\ \hfill x_2{y}^{-p/2}& \mathrm{if}{x}_2^2<y<x_3^2\hfill \end{array},\hfill \end{array}$$(56)

with x₁ = γ_m/γ_c, x₂ = 1 and x₃ = γ_max/γ_c. Calculations using this prescription are labelled “no SSC” in the following.

2.4.5. Synchrotron self-Compton case in the Thomson regime

If the KN regime is neglected, all IC scatterings occur in the Thomson regime, and the Compton parameter is the same for all electrons: Y(γ) = Y^no KN = cst. The corresponding solution was given by Sari & Esin (2001): The normalised distributions f(x) and g(y) are the same as in the purely synchrotron case above, but the value of the critical Lorentz factor γ_c is decreased due to the IC cooling, γ_c = γ_c^syn/(1+Y^no KN). Using Eqs. (36), (46), and (51), we obtain

$$\begin{array}{c}\hfill Y^{\mathrm{no}\mathrm{KN}}(1+Y^{\mathrm{no}\mathrm{KN}})=\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{p-2}{p-1}\frac{{ϵ}_{\mathrm{e}}}{{ϵ}_{\mathrm{B}}}\frac{{\gamma }_{\mathrm{m},\mathrm{c}}}{{\gamma }_{\mathrm{c},\mathrm{m}}}\frac{I_2}{I_0}.\end{array}$$(57)

As I₂/I₀ is a function of γ_c/γ_m = γ_c^syn/γ_m/(1 + Y^no KN), this is an implicit equation to be solved numerically to obtain Y^no KN and γ_c. The limits for (1+Y^no KN)Y^no KN are $\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{{ϵ}_{\mathrm{e}}}{{ϵ}_{\mathrm{B}}}$ for γ_m ≫ γ_c and $\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{{ϵ}_{\mathrm{e}}}{{ϵ}_{\mathrm{B}}}\frac{1}{3-p}{\left(\frac{{\gamma }_{\mathrm{m}}}{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}}\right)}^{p-2}$ for γ_m ≪ γ_c. Calculations using this prescription are labelled “SSC (Thomson)” in the following.

2.4.6. Self-consistent calculation of the electron distribution and the Compton parameter in the general case

2.4.6.1. Normalised Compton parameter.

In the general case where the KN suppression at high energy is included, the Compton parameter Y(γ) is a decreasing function of the electron Lorentz factor. From Eqs. (36), (46), and (51), we obtain

$$\begin{array}{c}\hfill Y(\gamma )(1+Y({\gamma }_{\mathrm{c}}))=\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{p-2}{p-1}\frac{{ϵ}_{\mathrm{e}}}{{ϵ}_{\mathrm{B}}}\frac{{\gamma }_{\mathrm{m},\mathrm{c}}}{{\gamma }_{\mathrm{c},\mathrm{m}}}\frac{I_2}{I_0}\frac{1}{J_0}{\displaystyle {\int }_0^{\tilde{y}}}\mathrm{d}yg(y).\end{array}$$(58)

For low Lorentz factors γ, $\stackrel{\sim }{\nu }$ becomes very large, so that for $\stackrel{\sim }{\nu }>{\nu }_{\mathrm{syn}}({\gamma }_{\max })$, we recover the Thomson regime, where the right-hand side of Eq. (61) is formally the same as in Eq. (60), even if the integrals I₀ and I₂ can be different when the normalised electron distribution f(x) is different. As described in Nakar et al. (2009), when only the leading term in the integral of g(y) is kept, the corresponding scaling law for the Compton parameter is Y(γ) = cst if $\stackrel{\sim }{\nu }$ is above the peak frequency ν_p and $Y(\gamma )\propto {\gamma }^{-\frac{3-\stackrel{\sim }{p}}{2}}$ otherwise, where $\stackrel{\sim }{p}$ is the slope of electron distribution in the power-law segment of the electron distribution including $\stackrel{\sim }{\gamma }$. We therefore introduce a normalised Compton parameter defined by

$$\begin{array}{c}\hfill Y(\gamma )(1+Y({\gamma }_{\mathrm{c}}))=\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{p-2}{p-1}\frac{{ϵ}_{\mathrm{e}}}{{ϵ}_{\mathrm{B}}}\frac{{\gamma }_{\mathrm{m},\mathrm{c}}}{{\gamma }_{\mathrm{c},\mathrm{m}}}\frac{I_2}{I_0}h(x),\end{array}$$(59)

where h(x) follows this scaling and is normalised so that Eq. (62) has the exact limit in the Thomson regime. This leads to

$$\begin{array}{c}\hfill h(x)=\{\begin{array}{cc}\hfill 1& \mathrm{if}x<{\widehat{x}}_{\mathrm{p}}\hfill \\ \hfill {\widehat{x}}_{\mathrm{p}}^{\frac{3-p_{i_{\mathrm{p}}-1}}{2}}{x}^{-\frac{3-p_{i_{\mathrm{p}}-1}}{2}}& \mathrm{if}{\widehat{x}}_{\mathrm{p}}<x<{\widehat{x}}_{i_{\mathrm{p}}-1}\hfill \\ \hfill {\widehat{x}}_{\mathrm{p}}^{\frac{3-p_{i_{\mathrm{p}}-1}}{2}}{\widehat{x}}_{i_{\mathrm{p}}-1}^{\frac{p_{i_{\mathrm{p}}-1}-p_{i_{\mathrm{p}}-2}}{2}}{x}^{-\frac{3-p_{i_{\mathrm{p}}-2}}{2}}& \mathrm{if}{\widehat{x}}_{i_{\mathrm{p}}-1}<x<{\widehat{x}}_{i_{\mathrm{p}}-2}\hfill \\ \hfill ⋮& ⋮\hfill \\ \hfill {\widehat{x}}_{\mathrm{p}}^{\frac{3-p_{i_{\mathrm{p}}-1}}{2}}{\widehat{x}}_{i_{\mathrm{p}}-1}^{\frac{p_{i_{\mathrm{p}}-1}-p_{i_{\mathrm{p}}-2}}{2}}\cdots & \\ \hfill \cdots {\widehat{x}}_2^{\frac{p_2-p_1}{2}}{x}^{-\frac{3-p_1}{2}}& \mathrm{if}{\widehat{x}}_2<x<{\widehat{x}}_1\hfill \\ \hfill {\widehat{x}}_{\mathrm{p}}^{\frac{3-p_{i_{\mathrm{p}}-1}}{2}}{\widehat{x}}_{i_{\mathrm{p}}-1}^{\frac{p_{i_{\mathrm{p}}-1}-p_{i_{\mathrm{p}}-2}}{2}}\cdots & \\ \hfill \cdots {\widehat{x}}_2^{\frac{p_2-p_1}{2}}{\widehat{x}}_1^{\frac{p_1}{2}-\frac{1}{6}}{x}^{-4/3}& \mathrm{if}x>{\widehat{x}}_1.\hfill \end{array}\end{array}$$(60)

We recall that i_p is the index of the break corresponding to the peak frequency of the synchrotron spectrum. In most cases (see Appendix B), we have i_p = 2, leading to

$$\begin{array}{c}\hfill h(x)=\{\begin{array}{cc}\hfill 1& \mathrm{if}x<{\widehat{x}}_2\hfill \\ \hfill {\widehat{x}}_2^{\frac{3-p_1}{2}}{x}^{-\frac{3-p_1}{2}}& \mathrm{if}{\widehat{x}}_2<x<{\widehat{x}}_1\hfill \\ \hfill {\widehat{x}}_2^{\frac{3-p_2}{2}}{\widehat{x}}_1^{\frac{p_1}{2}-\frac{1}{6}}{x}^{-4/3}& \mathrm{if}x>{\widehat{x}}_1.\hfill \end{array}\end{array}$$(61)

2.4.6.2. Self-consistent solution for the normalised electron distribution.

Following Nakar et al. (2009), we define²γ₀ by Y(γ₀) = 1 and assume 1 + Y(γ) = Y(γ) for γ < γ₀ and 1 otherwise. From Eqs. (40) and (41), this shows that the IC cooling only affects the electron distribution in the interval

$$\begin{array}{c}\hfill max({\gamma }_{\mathrm{c}};{\widehat{\gamma }}_{\mathrm{p}})<\gamma <{\gamma }_0.\end{array}$$(62)

Therefore, the solution is entirely determined by the orderings of γ_m, γ_c, γ₀, ${\widehat{\gamma }}_{\mathrm{m}}$, and $\widehat{{\gamma }_{\mathrm{c}}}$. All possible cases and the corresponding solutions are listed in Appendix B. In some cases, subcases are introduced as breaks can appear at ${\widehat{\gamma }}_0$, ${\widehat{\widehat{\gamma }}}_{\mathrm{m}}$, ${\widehat{\widehat{\gamma }}}_{\mathrm{c}}$, etc. All the most relevant cases for GRB afterglows were described in Nakar et al. (2009), along with the detailed method for obtaining the corresponding solution for the electron distribution. For completeness, we list in Appendix B the additional cases that allow us to fully describe the parameter space, even in regions that are unlikely to be explored in GRBs.

2.4.6.3. Calculation of the critical Lorentz factor γ_c.

From Eq. (62), we obtain the following equation for Y(γ_c):

$$\begin{array}{c}\hfill Y({\gamma }_{\mathrm{c}})=\{\begin{array}{cc}A\frac{{\gamma }_{\mathrm{m}}}{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}}\frac{I_2}{I_0}h(x_{\mathrm{c}})& \mathrm{if}{\gamma }_{\mathrm{m}}>{\gamma }_{\mathrm{c}}\hfill \\ A\frac{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}}{{\gamma }_{\mathrm{m}}}\frac{I_2}{I_0}h(x_{\mathrm{c}})& \mathrm{if}{\gamma }_{\mathrm{m}}<{\gamma }_{\mathrm{c}}\mathrm{and}Y({\gamma }_{\mathrm{c}})<1\hfill \\ {\left(A\frac{{\gamma }_{\mathrm{c}}^{\mathrm{syn}}}{{\gamma }_{\mathrm{m}}}\frac{I_2}{I_0}h(x_{\mathrm{c}})\right)}^{1/3}& \mathrm{if}{\gamma }_{\mathrm{m}}<{\gamma }_{\mathrm{c}}\mathrm{and}Y({\gamma }_{\mathrm{c}})>1,\hfill \end{array}\end{array}$$(63)

where $A=\frac{3}{4}\frac{K_{\mathrm{IC}}}{K_P}\frac{p-2}{p-1}\frac{{ϵ}_{\mathrm{e}}}{{ϵ}_{\mathrm{B}}}$ is a constant as the microphysics parameters ϵ_B, ϵ_e, and p are assumed to be constant during the forward-shock propagation. We note that x_c = 1 in the slow-cooling regime γ_m < γ_c. A second useful relation is obtained from Eq. (62) and the definition of γ₀,

$$\begin{array}{c}\hfill Y({\gamma }_{\mathrm{c}})=\frac{Y({\gamma }_{\mathrm{c}})}{Y({\gamma }_0)}=\frac{h(x_{\mathrm{c}})}{h(x_0)}.\end{array}$$(64)

Then, γ_c is computed from the following implicit equation, which is solved iteratively,

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{c}}=\{\begin{array}{cc}\hfill {\gamma }_{\mathrm{c}}^{\mathrm{syn}}& \mathrm{if}Y({\gamma }_{\mathrm{c}})<1\hfill \\ \hfill {\gamma }_{\mathrm{c}}^{\mathrm{syn}}/Y({\gamma }_{\mathrm{c}})& \mathrm{if}Y({\gamma }_{\mathrm{c}})>1.\hfill \end{array}\end{array}$$(65)

We start by assuming Y(γ_c) = Y^no KN as defined in Sect. 2.4.5 and then iterate as follows:

1. Compute γ_c from Eq. (68).

2. Knowing γ_m, γ_c, and γ_self, scan the different possible cases for the ordering of γ_m, γ_c, γ₀, ${\widehat{\gamma }}_{\mathrm{m}}$, and $\widehat{{\gamma }_{\mathrm{c}}}$ (and additional characteristic Lorentz factors for some cases) as listed in Appendix B and compute γ₀ from Eq. (67), which can be analytically inverted as provided for each case in Appendix B. We stop the scan of the possible cases when the ordering of all characteristic Lorentz factors including the obtained value of γ₀ is the correct one.

3. Having identified the correct case and therefore knowing the expression of f(x), compute the integrals I₀ and I₂.

4. Compute an updated value of Y(γ_c) from Eq. (66).

5. Start a new iteration at step 1 until convergence.

We stop the procedure when the relative variation of γ_c in an iteration falls below ϵ_tol. We use ϵ_tol = 10⁻⁴, which is reached in most cases in fewer than 20 iterations. In the context of a full afterglow light-curve calculation, as we expect Y to vary smoothly during the jet propagation, we use the value of Y(γ_c) obtained at the previous step of the dynamics instead of Y^no KN to start the iterative procedure. Then, we usually reach convergence in only a few iterations.

2.4.7. Final calculation of the emission in the comoving frame

The self-consistent spectrum in the comoving frame, that is, p_ν^syn and p_ν^SSC, the synchrotron and SSC powers per electron and per unit frequency averaged over the timescale t_dyn, can be computed by the following procedure based on the previous paragraphs:

1. The input parameters are $n_{\mathrm{e}}^{\mathrm{acc}}$, t_dyn, B, u_e/u_B = ϵ_e/ϵ_B, γ_m, and p.

2. The following quantities can immediately be deduced: γ_c^syn from Eq. (39), γ_self from Eq. (34), τ_T from Eq. (54), and Y^no KN from Eq. (60).

3. Then, the iterative procedure described above allows us to identify the spectral regime among all the possibilities listed in Appendix B, and we can compute γ_c. At the end of this step, all breaks and slopes in the functions f(x) and h(x) are known. The corresponding integrals I₀ and I₂ are computed by a numerical integration of Eqs. (44) and (48).

4. The function g(y) can immediately be deduced from f(x), using Eq. (47). The corresponding integral J₀ is computed by a numerical integration of Eq. (49), and the synchrotron spectrum per electron, p_ν^syn, is deduced from Eq. (46).

5. Finally, the function G(y) is computed by an numerical integration of Eq. (55) and the SSC spectrum per electron, p_ν^SSC, is deduced from Eq. (53).

2.4.8. Pair production

At VHE, pair production $\gamma \gamma {\stackrel{e}{\to }}^{+}{e}^{-}$ can start to contribute to the SSC flux depletion. We account for this mechanism by a simplified treatment, assuming an isotropic distribution of the low-energy seed photons, as for the SSC calculation, and approximating the cross-section for pair production by a Dirac function at twice the threshold of the interaction. Then, VHE photons of frequency ν can produce pairs by interacting with low-energy photons at a frequency ν_seed = 2(m_ec²)²/(h²ν), and the characteristic timescale of this interaction is given by

$$\begin{array}{c}\hfill t_{\gamma \gamma }(\nu )=\frac{h}{{\sigma }_{\mathrm{T}}c}{\left[u_{{\nu }_{\mathrm{seed}}=2\frac{{\left(m_{\mathrm{e}}{c}^2\right)}^2}{h^2\nu }}\right]}^{-1},\end{array}$$(66)

with u_ν = u_ν^syn + u_ν^SSC = n_e^acc(p_ν^syn+p_ν^SSC)t_dyn, leading to

$$\begin{array}{c}\hfill t_{\gamma \gamma }(\nu )=\frac{h}{{\tau }_{\mathrm{T}}}{\left[{(p_{{\nu }_{\mathrm{seed}}}^{\mathrm{syn}}+p_{{\nu }_{\mathrm{seed}}}^{\mathrm{SSC}})}_{{\nu }_{\mathrm{seed}}=2\frac{{\left(m_{\mathrm{e}}{c}^2\right)}^2}{h^2\nu }}\right]}^{-1}.\end{array}$$(67)

For very high frequencies ν, the seed photons for pair production have a low frequency and p_νseed^syn + p_νseed^SSC ≃ p_νseed^syn. The total emitted power per electron p_ν = p_ν^syn + p_ν^SSC is then corrected for by a factor $\frac{t_{\gamma \gamma }(\nu )}{t_{\mathrm{dyn}}}(1-e^{-t_{\mathrm{dyn}}/t_{\gamma \gamma }(\nu )})$, which in practice only attenuates the SSC component at high frequency, where t_γγ(ν)≪t_dyn.

The implementation of the model allows us to select the level of approximation for the SSC emission (synchrotron only, Thomson, or full calculation) and to include or exclude the attenuation due to the pair production.

2.5. Observed flux

When the emissivity in the comoving frame is known (Sect. 2.4), the flux measured by a distant observer with a viewing angle θ_v is computed by an integration over equal-arrival time surfaces, taking into account relativistic Doppler boosting and relativistic beaming, as well as the effect of cosmological redshift. This leads to the following expression of the observer frame flux density (erg ⋅ s⁻¹ ⋅ cm⁻² ⋅ Hz⁻¹) measured at time t_obs^z and frequency ν_obs^z(Woods & Loeb 1999):

$$\begin{array}{cc}\hfill F_{{\nu }_{\mathrm{obs}}^z}\left(t_{\mathrm{obs}}^z\right)& =\frac{1+z}{4\pi D_{\mathrm{L}}^2}{\displaystyle {\int }_0^{\infty }}\mathrm{d}r{\displaystyle {\int }_0^{\pi }}\mathrm{d}\psi {\displaystyle {\int }_0^{2\pi }}\mathrm{d}\varphi r^{2}sin\psi \hfill \\ \hfill & \times {\left[{\mathcal{D}}^2(r,\psi ,t)4\pi j_{{\nu }^{\prime }}^{\prime }(r,\psi ,\varphi ,t)\right]}_{t=\frac{t_{\mathrm{obs}}^z}{1+z}+\frac{r}{c}cos\psi },\hfill \end{array}$$(68)

where we use spherical coordinates (r;ψ;ϕ) with the polar axis equal to the line of sight, ψ and ϕ the colatitude and longitude, and the jet axis in the direction (ψ=θ_v,ϕ = 0) (see Fig. 1), and where z and D_L are the redshift and the luminosity distance of the source. The angle ψ differs from θ used in Sect. 2.2, which is measured from the jet axis. The time t (source frame) and frequency ν′ (comoving frame) are given by

$$\begin{array}{cc}\hfill t& =\frac{t_{\mathrm{obs}}^z}{1+z}+\frac{r}{c}cos\psi \hfill \end{array}$$(69)

$$\begin{array}{cc}\hfill {\nu }^{\prime }& =(1+z)\frac{{\nu }_{\mathrm{obs}}^z}{\mathcal{D}(r,\psi ,t)}.\hfill \end{array}$$(70)