© ESO, 2014

1. Introduction

Gamma-ray burst (GRB) afterglows are produced by relativistic electrons radiating mainly via the synchrotron and inverse Compton mechanisms. According to the standard afterglow model, the electrons are accelerated to highly relativistic energies at two shock fronts, the forward shock and the reverse shock, which are the result of the interaction between the relativistic jet from the GRB central engine and the surrounding medium (for reviews, see, e.g., Piran 2004; Mészáros 2006).

The earliest afterglow models invoke pure synchrotron radiation from the forward shock in a constant-density interstellar medium (ISM) or a wind-type environment and yield analytic time-evolving synchrotron spectra of the decelerating blast wave (Sari et al. 1998; Chevalier & Li 2000; Granot & Sari 2002). The role of inverse Compton scattering of the synchrotron photons has also been investigated, typically relying on an approximate treatment of the scattering process because no analytic solution for the inverse Compton spectrum is available (Panaitescu & Meszaros 1998; Chiang & Dermer 1999; Panaitescu & Kumar 2000; Sari & Esin 2001).

The GRB observations by the Swift and Fermi satellites have revealed some surprising features in the afterglow and prompt light curves, resulting in a need to improve the models for GRB emission. For example, the Fermi/LAT telescope has observed >100 MeV emission from several GRBs. In the literature, the high-energy radiation has been attributed to the prompt emission (e.g., Abdo et al. 2009), to pure synchrotron radiation from the external shock (Gao et al. 2009; Ghisellini et al. 2010; Kumar & Barniol Duran 2010), to a combination of external synchrotron photons and synchrotron self-Compton emission (Tam et al. 2013; Wang et al. 2013; Liu et al. 2013; Fan et al. 2013), and to a superposition of the prompt and afterglow emission (Maxham et al. 2011). Another possibility is that some of the high-energy emission stems from prompt photons being Compton scattered to higher energies by the afterglow-emitting electrons (Beloborodov 2005b; He et al. 2012; Beloborodov et al. 2013; Fan et al. 2013).

Owing to the Swift observations, it has been discovered that a typical X-ray light curve begins with a phase of steeply decaying flux, which is often followed by a shallow decay segment. The late-time afterglow, on the other hand, can often be explained by the standard synchrotron model (Granot & Kumar 2006; Fan & Piran 2006). Energy injection to the blast wave is currently the most popular explanation for the shallow decay phase observed both in X-ray and optical afterglows (Granot & Kumar 2006; Fan & Piran 2006; Nousek et al. 2006; Zhang et al. 2006; Panaitescu & Vestrand 2011; Li et al. 2012). Other models introduced to explain the shallow decay phase include the evolution of microphysical parameters (Panaitescu et al. 2006; Granot et al. 2006), emission due to an outflow ejected before the prompt GRB (Yamazaki 2009; Birnbaum et al. 2012), an off-axis viewing angle of the jet (Eichler & Granot 2006), dust scattering of X-rays (Shao & Dai 2007), late prompt emission (Ghisellini et al. 2007; Murase et al. 2011), and an adiabatic evolution of the shock following a radiative phase (Dermer 2007). It has also been suggested that the main contribution to the afterglow could come from a long-lived reverse shock, which may also explain the shallow decay phase in the X-ray afterglows (Uhm & Beloborodov 2007; Genet et al. 2007).

Models aiming to explain all the different slopes seen in the light curves have also been presented, including accretion of different layers of the progenitor star (Kumar et al. 2008) and the curvature effect that is usually only invoked to explain the early steep X-ray decay (Qin 2008).

The evolution of the GRB blast wave is described well by the self-similar solution by Blandford & McKee (1976), which is valid in the deceleration phase while the blast is still highly relativistic. The evolution in the late non-relativistic phase is given by the Sedov-Taylor solution (Sedov 1959; Taylor 1950). A mechanical model of the relativistic blast ensuring mass, energy, and momentum conservation has been presented by Beloborodov & Uhm (2006), and it is nearly identical to the Blandford-McKee solution at late times after the shock has started to decelerate. However, the mechanical model gives an accurate description of the blast also before the deceleration time, while earlier models unphysically assume an equal pressure at the forward and reverse shock. The evolution of the shell in the mildly relativistic phase can be found by means of hydrodynamic simulations, which can then be coupled to a radiation code to find the radiation spectrum from the shock. Such simulations can also be applied to calculate the afterglow emission for an observer with an off-axis viewing angle (van Eerten et al. 2010a).

Results of one- and two-dimensional hydrodynamic simulations of the blast wave have been presented by Kobayashi et al. (1999), Meliani et al. (2007), Mimica et al. (2009), and Ramirez-Ruiz & MacFadyen (2010) but without discussing the radiation mechanism of the afterglow. A calculation of the synchrotron radiation from the blast has been coupled to the hydrodynamic simulations of Downes et al. (2002), Zhang & MacFadyen (2009), van Eerten et al. (2010b, 2011), and Wygoda et al. (2011). However, none of these works calculate the afterglow component due to Compton scattering, which is expected to appear at high energies.

Simulations including an accurate treatment of both synchrotron and Compton processes, as well as pair production, have been presented by Petropoulou & Mastichiadis (2009) (PM09), who use the solution of Blandford & McKee (1976) to evaluate the evolution of the emitting shell. The code developed by PM09 is similar to the one presented in this paper, but it does not account for the electron heating due to synchrotron self-absorption.

For the first time, we present simulations of afterglow emission from the forward shock with a relativistic kinetic code that treats synchrotron emission and absorption, Compton scattering, and electron-positron pair production/annihilation in a self-consistent way. The kinetic equations determining the time evolution of the electron and photon distributions are solved simultaneously at each timestep. We also consider electron heating due to synchrotron self-absorption, which shapes the electron distribution at low energies.

Our treatment accounts for the fact that electrons injected at different times also have different cooling histories. For example, the magnetic field that determines the synchrotron cooling rate evolves while the electrons are cooling. It follows that there are no sharp cooling breaks in the electron distribution (Uhm & Zhang 2014).

The current version of the code applies a one-zone model of the emission region. It does not account for the different locations of the particles behind the shock and assumes a constant magnetic field throughout the shell. A more accurate treatment of synchrotron emission would require a model of the magnetic field structure behind the shock. Also, knowledge of the spatial photon and electron distributions is required for an exact calculation of Compton scattering.

As an example of the applications of the code, we report the results of a simulation where the afterglow-emitting electrons interact with an external source of photons roughly corresponding to prompt GRB emission. The shocked electrons are expected to upscatter a small fraction of the prompt photons to GeV−TeV energies as long as the prompt emission overlaps with the shocked electrons (Beloborodov 2005b; Fan et al. 2005). Some of the high-energy photons then produce pairs with the prompt MeV photons, which in turn are able to scatter radiation to higher energies.

In addition, we compare the forward shock emission in a wind environment with the emission in a constant-density ISM. The results indicate that a power-law electron distribution can thermalize at low energies thanks to synchrotron self-absorption heating in a wind medium with a typical density structure expected from the surroundings of a Wolf-Rayet star. Along with the ambient density, the importance of thermalization is mainly determined by the fraction of shock-generated energy given to the magnetic field. Our simulations imply that the thermalized electrons are unlikely to produce an observable signature in the afterglow spectrum.

In our final example, we study the difference between the forward shock radiation due to Maxwellian and power-law electron injection. The standard afterglow model assumes that the injection function is a pure power law, even though a large fraction of the shock-generated energy goes to a thermal population of electrons. We find that pure Maxwellian injection can lead to a flattening in the X-ray light curve. The flux during this phase is found to be very low compared to Swift/XRT detections for a constant-density ISM, but detectable flux levels during the shallow decay may be achieved in a wind-type environment.

2. Physical model of the afterglow

2.1. Hydrodynamic evolution

The relativistic shell ejected by the GRB central engine initially propagates into the surrounding medium with a constant Lorentz factor Γ₀. After sweeping an external mass M₀/Γ₀, where M₀ is the initial mass of the ejecta, the shell starts to decelerate according to a self-similar solution found by Blandford & McKee (1976). The self-similar solution is no longer valid if the reverse shock is long-lived (Uhm et al. 2012) or if there is significant lateral spreading of the shell after the jet break time (Rhoads 1999). In the rest of the paper, except for Sects. 3.1 and 3.3, the quantities measured in the fluid comoving frame are indicated by primes and the unprimed quantities are given in the observer frame. However, the electron Lorentz factor γ and dimensionless momentum $\mathit{p}\mathrm{\equiv }\sqrt{{\mathit{\gamma }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\mathit{,}$(1)along with the electron number density per unit energy (N(γ)) or momentum (N(p)) interval, are always expressed in the fluid frame but left unprimed to simplify the notation. In the following, we assume an adiabatic hydrodynamic evolution of the blast wave, which means negligible radiation losses. This assumption is valid even if the electrons cool rapidly as long as the main fraction of the shock-generated energy is given to the protons, which do not radiate efficiently.

The bulk Lorentz factor of the shocked fluid evolves with radius r approximately as $\mathrm{\Gamma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{(}\mathit{r}\mathit{/}{\mathit{R}}_{\mathit{B}}{\mathrm{)}}^{\mathrm{1}\mathrm{-}\mathit{g}}\mathit{,}$(2)where the value of the index g depends on the structure of the ambient medium. If the medium has a constant density, $\begin{array}{ccc}\mathit{g}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(3)where R_B ≡ R_dec/2^2/3 is the approximate radius at which the Lorentz factor starts to decrease from its initial value Γ₀, and R_dec is the so-called deceleration radius, i.e., the radius where the shell has lost half of its initial kinetic energy and Γ = Γ₀/2 (Rees & Meszaros 1992). According to this definition, ${\mathit{R}}_{\mathrm{dec}}\mathrm{=}{\left(\frac{\mathrm{3}{\mathit{E}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }{\mathit{n}}_{\mathrm{0}}{\mathit{m}}_{\mathrm{p}}{\mathit{c}}^{\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{0}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathit{,}$(4)where E₀ is the isotropic equivalent energy of the blast wave after the prompt GRB emission, and n₀ the electron number density of the external medium.

In the case of a wind-type medium with a density profile n(r) = A_wr^-2, where A_w is a constant giving the number of particles per unit length, $\begin{array}{ccc}\mathit{g}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(5)where R_B ≡ R_dec/4, and ${\mathit{R}}_{\mathrm{dec}}\mathrm{=}\frac{{\mathit{E}}_{\mathrm{0}}}{\mathrm{4}\mathit{\pi }{\mathit{A}}_{\mathrm{w}}{\mathit{m}}_{\mathrm{p}}{\mathit{c}}^{\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{·}$(6)The shock radius and the Lorentz factor of the shocked ejecta are related to the fluid comoving time t′ according to $\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{d}\mathit{r}}{\mathit{c}\mathrm{\Gamma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(7)and the observed time can be obtained from $\mathrm{d}\mathit{t}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{d}\mathit{r}}{\mathrm{2}\mathit{c}{\mathrm{\Gamma }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(8)where z is the redshift of the GRB. Before the deceleration time, i.e., r<R_dec, $\mathit{r}\mathrm{=}\frac{\mathrm{2}\mathit{c}{\mathrm{\Gamma }}_{\mathrm{0}}^{\mathrm{2}}\mathit{t}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathrm{·}$(9)For r ≫ R_dec, the radius evolves in time as $\begin{array}{ccc}\mathit{r}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(10)and the bulk Lorentz factor evolves as $\begin{array}{ccc}\mathrm{\Gamma }\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(11)One of the main uncertainties of the afterglow model is the strength of the magnetic field in the emission region. We adopt the typical approach to the problem and assume that a fraction ϵ_B of the shock-generated energy goes into the energy of the shock-compressed interstellar magnetic field. The comoving magnetic field B′ then evolves in the emission region with Γ and n as ${\mathit{B}}^{\mathrm{\prime }}\mathrm{=}\mathit{c}\mathrm{\Gamma }\sqrt{\mathrm{32}\mathit{\pi }{\mathit{m}}_{\mathrm{p}}{\mathit{ϵ}}_{\mathit{B}}\mathit{n}}\mathit{,}$(12)determined from the shock jump conditions (Blandford & McKee 1976).

2.2. Distribution of the accelerated electrons

The forward and reverse shocks are thought to be capable of accelerating electrons to relativistic energies according to a power-law distribution (e.g., Achterberg et al. 2001) of the form $\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{\gamma }}\mathrm{\equiv }\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{\propto }{\mathit{\gamma }}^{\mathrm{-}\mathit{s}}$(13)for γ_min ≤ γ ≤ γ_max, where s is the power-law index, and γ the random Lorentz factor of an electron measured in the fluid frame. However, simulations by Spitkovsky (2008) and Martins et al. (2009) show that shock acceleration actually leads to a hybrid distribution with low-energy Maxwellian electrons connected to a high-energy power-law tail. An afterglow model that assumes such a distribution is discussed by Giannios & Spitkovsky (2009).

The number of electrons injected at a shock per unit time is obtained by multiplying the number density of electrons in the downstream frame, Γn, by the volume swept by the shock per unit time, 4πr²βc. Dividing this quantity by the volume of the shell, 4πr²ΔR′, one finds that the number of injected electrons per unit time and volume is $\mathit{ṅ}\begin{array}{c}\mathrm{\prime }\\ \mathrm{el}\end{array}\mathrm{=}\frac{\mathrm{\Gamma }\mathit{n\beta c}}{\mathrm{\Delta }{\mathit{R}}^{\mathrm{\prime }}}\mathrm{·}$(14)Assuming a pure power-law distribution, the minimum Lorentz factor determined from the shock jump conditions (Sari et al. 1996) is ${\mathit{\gamma }}_{\min }\mathrm{=}{\mathit{ϵ}}_{\mathrm{e}}\frac{\mathit{s}\mathrm{-}\mathrm{2}}{\mathit{s}\mathrm{-}\mathrm{1}}\frac{{\mathit{m}}_{\mathrm{p}}}{{\mathit{m}}_{\mathrm{e}}}\mathrm{\Gamma }\mathit{,}$(15)where ϵ_e is the fraction of the shock-generated energy given to the electrons. This relation holds for γ_min ≪ γ_max.

The maximum Lorentz factor γ_max is determined either by the radiation losses of the electrons, the age of the flow, or the saturation limit discussed by Sironi et al. (2013), among others. Typically it has been assumed that the value of γ_max and the corresponding synchrotron frequency are too high to affect the observed afterglow properties. After the detection of high-energy (>100 MeV) emission from several GRBs, which may be part of the synchrotron or inverse Compton component of the afterglow, it has become more important to determine γ_max accurately to study whether there are electrons at high enough energies to produce >100 MeV synchrotron radiation.

2.3. Radiation processes

2.3.1. Synchrotron emission and self-absorption

The cooling of the shocked electrons is mainly determined by synchrotron losses, with a significant contribution from adiabatic cooling to the distribution of the low-energy electrons. It is expected that the cooled electrons are distributed according to N(γ) ∝ γ^−s−1 above the injection energy γ_min (Kardashev 1962). If the electrons are cooling to energies below the injection energy (the fast cooling case), N(γ) ∝ γ^-2 for γ<γ_min as long as the electrons do not thermalize due to self-absorption heating.

During the late stages of a typical afterglow, the magnetic field is low, and only the most energetic electrons are cooling (slow cooling). The uncooled electrons below a critical energy γ_c are distributed as N(γ) ∝ γ^−s similarly to the injection function. The different parts of the electron distribution correspond to the power-law segments in the emergent radiation spectrum as described by Granot & Sari (2002). The characteristic observed synchrotron photon frequency of a relativistic electron is $\mathit{\nu }\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Gamma }{\mathit{\gamma }}^{\mathrm{2}}\mathit{e}{\mathit{B}}^{\mathrm{\prime }}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{2}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}\mathit{c}}\mathrm{·}$(16)The spectral slope of the electron distribution gradually changes around the cooling energy γ_c, which is defined by Sari et al. (1998) according to the relation ${\mathit{\gamma }}_{\mathrm{c}}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\mathrm{=}{\mathit{P}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\gamma }}_{\mathrm{c}}\mathrm{\right)}{\mathit{t}}^{\mathrm{\prime }}\mathit{,}$(17)where ${\mathit{P}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\gamma }}_{\mathrm{c}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}{\mathit{\sigma }}_{\mathrm{T}}\mathit{c}{\mathit{\gamma }}^{\mathrm{2}}\frac{\mathrm{(}{\mathit{B}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{8}\mathit{\pi }}$(18)is the synchrotron power of an electron with γ ≫ 1 in the comoving frame. The relation is based on the definition that a cooling electron loses an energy comparable to its own initial energy in the lifetime of the flow. Correspondingly, the comoving cooling time of an electron with a Lorentz factor γ is ${\mathit{t}}_{\mathrm{cool}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{6}\mathit{\pi }{\mathit{m}}_{\mathrm{e}}\mathit{c}}{{\mathit{\sigma }}_{\mathrm{T}}\mathrm{(}{\mathit{B}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}\mathit{\gamma }}\mathrm{·}$(19)In reality, a sharp break at the energy γ_c does not appear because the actual electron distribution consists of different electron populations with their own cooling histories. The radiation power P′(γ_c) evolves in time as the magnetic field B′ decreases, and an integration over P′(γ_c)dt′ is required to find the exact energy lost by each electron population. Both our code and the one by Petropoulou & Mastichiadis (2009) calculate the shape of the electron distribution according to the kinetic equation, and the results show that the transition between the cooled and uncooled segments in the distribution is very gradual. This effect is taken into account by, say, Uhm & Zhang (2014), who also point out that the different locations of the electron populations behind the shock shape the outgoing radiation spectrum. In our current one-zone treatment, the structure of the electron distribution and the magnetic field behind the shock are not resolved. A multi-zone approach is expected to contribute further to the curvature of the particle distributions, and the shape of the resulting spectrum will be a result of the emissions from varying distances behind the shock by electron populations with different cooling histories. Additional smoothing of the cooling break would also be provided by emission from large angles, but the code assumes at this stage that all the emission comes from the line of sight to the observer.

In addition to being cooled by synchrotron emission, low-energy electrons can be heated due to self-absorption of the synchrotron photons (e.g., Ghisellini et al. 1988), leading to some thermalization of the electron distribution. Full thermalization of the low-energy electrons takes roughly one synchrotron cooling time. The cooling time for a given electron energy depends on the magnetic field as ∝(B′)^-2 ∝ n^-1Γ^-2 (see Eq. (12)). Because the density in a wind-type medium at small radii is considerably larger than in a typical constant-density ISM, the magnetic field is also higher and correspondingly the cooling time is shorter. This implies that electron thermalization due to self-absorption is more likely to occur in a wind environment.

2.3.2. Inverse Compton scattering

Some of the synchrotron afterglow photons are Compton scattered to higher energies by the same electrons that emit the synchrotron radiation. If the prompt GRB photons overlap with the afterglow-emitting region, some of these photons are also upscattered by the relativistic electrons up to >TeV energies. Compton scattering depends on the photon density and consequently the geometry of the emission region, which is discussed in Sect. 3.2.

In the Thomson regime, the location of the inverse Compton peak is determined by the scattering of the peak photons of the νF_ν Band spectrum against the peak electrons of the γ²N(γ) distribution. However, if the scattering of the peak photons against the peak electrons is suppressed due to the Klein-Nishina (K-N) effect, the peak of the Compton scattered photons is located at a lower energy. Noting that the location of the peak should correspond to the highest possible energy that can be transferred to a photon from a peak electron with a Lorentz factor γ = γ_pk, the observed peak energy becomes ${\mathit{E}}_{\mathrm{pk}\mathit{,}\mathrm{IC}}\mathrm{~}\frac{\mathrm{2}\mathrm{\Gamma }{\mathit{\gamma }}_{\mathrm{pk}}{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathrm{·}$(20)

2.3.3. Pair production

Electron-positron pair production and annihilation may have some impact on the observed afterglow spectra, especially if a fraction of the prompt photons is Compton scattered to high energies, after which they are able to produce pairs with the unscattered GRB photons. Defining the dimensionless photon energy as $\mathit{x}\mathrm{\equiv }\frac{\mathit{h\nu }}{{\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}}\mathit{,}$(21)high-energy photons of energy x_HE produce pairs with target photons exceeding the threshold energy ${\mathit{x}}_{\mathrm{thr}}\mathrm{~}\frac{\mathrm{4}{\mathrm{\Gamma }}^{\mathrm{2}}}{{\mathit{x}}_{\mathrm{HE}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(22)This value of x_thr is obtained from the threshold condition x_thrx_HE(1 + z)²(1 − cosθ) > 2, where θ is the half-angle inside which the upscattered photons propagate. The condition states that the invariant four-product of the photon momenta should be greater than 2 to be able to produce two leptons with zero kinetic energy in the center of momentum frame. With a typical angle θ ~ 1/Γ, one obtains Eq. (22). For target photons just above x_thr, the optical depth for pair production reaches its maximum value (e.g., Zdziarski 1988) ${\mathit{\tau }}_{\mathit{\gamma \gamma }}\mathrm{\approx }\frac{{\mathit{\sigma }}_{\mathrm{T}}}{\mathrm{5}}{\mathit{n}}_{\mathrm{ph}}\mathit{r}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\theta }\mathrm{)}\mathrm{\approx }\frac{{\mathit{\sigma }}_{\mathrm{T}}}{\mathrm{5}}{\mathit{n}}_{\mathrm{ph}}\frac{\mathit{r}}{\mathrm{2}{\mathrm{\Gamma }}^{\mathrm{2}}}\mathit{,}$(23)where n_ph is the number density of the target photons above the threshold energy.

For hard bursts, the pair production opacity is maximal for high-energy photons for which the threshold energy is x_thr ~ x_pk, where x_pk is the peak energy of the νF_ν prompt GRB spectrum. The distribution of prompt photons is typically described by the so-called Band function (Band et al. 1993) $\begin{array}{ccc}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{x}}\mathrm{\propto }\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(24)where dn/dx is the photon number density per dimensionless energy interval. Because the number density of the photons of energy x is n_ph ~ xdn/dx ∝ x^α+1 or n_ph ∝ x^β+1 and the opacity is proportional to n_ph according to Eq. (23), one finds that ${\mathit{\tau }}_{\mathit{\gamma \gamma }}\mathrm{\propto }{\mathit{x}}_{\mathrm{thr}}^{\mathit{\alpha }\mathrm{+}\mathrm{1}}\mathrm{\propto }{\mathit{x}}_{\mathrm{HE}}^{\mathrm{-}\mathit{\alpha }\mathrm{-}\mathrm{1}}$ for ${\mathit{x}}_{\mathrm{HE}}\mathit{>}\mathrm{4}{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{x}}_{\mathrm{pk}}^{-1}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-2}$ and ${\mathit{\tau }}_{\mathit{\gamma \gamma }}\mathrm{\propto }{\mathit{x}}_{\mathrm{thr}}^{\mathit{\beta }\mathrm{+}\mathrm{1}}\mathrm{\propto }{\mathit{x}}_{\mathrm{HE}}^{\mathrm{-}\mathit{\beta }\mathrm{-}\mathrm{1}}$ for ${\mathit{x}}_{\mathrm{HE}}\mathrm{\le }\mathrm{4}{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{x}}_{\mathrm{pk}}^{-1}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-2}$. For a typical GRB with α ~ −1, the opacity remains relatively constant at ${\mathit{x}}_{\mathrm{HE}}\mathit{>}\mathrm{4}{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{x}}_{\mathrm{pk}}^{-1}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-2}$. One also notes that the opacity increases for softer and decreases for harder bursts.

The number density n_ph is dominated by the prompt photons while they are going through the shell of ejecta and can be estimated as ${\mathit{n}}_{\mathrm{ph}}\mathrm{~}\frac{\mathit{L}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{c}{\mathit{E}}_{\mathrm{pk}}}\mathit{,}$(25)which is obtained by dividing the prompt GRB luminosity L by the volume covered by the photons per unit time, 4πr²c, and by the average energy of a photon, which we now assume to be the prompt peak energy E_pk.

Using the expression for n_ph together with Eq. (23), one finds that τ_γγ ∝ r^-1Γ^-2 for high-energy photons for which the threshold energy is E_pk. The prompt photons mainly overlap with the ejected shell during the coasting phase when Γ = Γ₀ and r ∝ t. In this case, the optical depth evolves as τ_γγ ∝ r^-1 ∝ t^-1. After the deceleration time, r ∝ t^1/4 and Γ ∝ r^−3/2 in the ISM case, from which it follows that τ_γγ ∝ r² ∝ t^1/2. The probability for pair production is thus expected to decrease during the coasting phase but increase after R_B if the prompt photons are still passing through the emission region. The rising opacity is due to the increasing angular spread of the emitted high-energy photons, which offsets the decrease due to the declining target photon density. Taking into account the beaming of the prompt radiation, the suppression of pair production is not exponential even at τ_γγ > 1. In this case there always exists an escape cone at sufficiently small angles, within which high-energy photons can escape. In a wind medium, where Γ ∝ r^−1/2, τ_γγ ∝ r⁰ after the deceleration time, that is to say, the optical depth remains constant.

The above discussion on the pair production opacity applies to the high-energy photons for which the threshold energy of target photons is E_pk, and the energy of the former is generally changing with time. However, if one assumes the canonical photon index α ~ −1, the photon number density n_ph ∝ x^α+1 = x⁰ below the peak energy. In this case, the density of target photons below E_pk is still given by Eq. (25) and the pair production opacities obtained above apply to any high-energy photon for which the threshold energy is below E_pk.

Once the prompt photons have crossed the emission region, the synchrotron emission provides the target photons for pair production with upscattered high-energy photons. The number of synchrotron photons providing the pair production opacity is proportional to the number of electrons emitting at x_thr. The energy of these electrons can be estimated from ${\mathit{x}}_{\mathrm{thr}}\mathrm{=}\mathrm{4}{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{x}}_{\mathrm{HE}}^{-1}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-2}\mathrm{=}{\mathit{x}}_{\mathrm{syn}}\mathrm{\propto }{\mathrm{\Gamma }}^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}$ (Eq. (16)), where x_syn is the dimensionless synchrotron photon energy. This yields that the electron energy corresponding to the threshold photon energy is γ_thr = constant for a fixed x_HE. In the slow cooling regime the number of these electrons is N_e,thr ∝ r³(γ_thr/γ_min)^−s+1 ∝ r³Γ^s−1 (Eq. (15)). For example, in the ISM case with s = 2 one gets an optical depth τ_γγ ∝ r²Γ ∝ t^1/8 (Eqs. (10) and (11)); i.e., the pair production opacity is increasing in time, although very slowly.

3. Numerical treatment

3.1. Kinetic equations

The numerical code we use to simulate GRB afterglows is based on the code developed by Vurm & Poutanen (2009), which has successfully been applied to, say, modeling prompt GRB (Vurm et al. 2011) and black hole emission (Veledina et al. 2011, 2013). The code calculates the time-evolving particle distributions in an astrophysical plasma by solving the full relativistic kinetic equations for each particle species without any energy limitations. This means that the equations include both differential and integral terms, depending on the nature of each radiation process. If the energy losses of a particle are approximately continuous, differential terms are used to describe the process. An integral term is necessary if the particle loses a considerable fraction of its energy due to a single interaction. In this section, all quantities except for r and Γ are expressed in the fluid comoving frame and are left unprimed for ease of notation.

The equations for photons and electrons in the fluid comoving frame both have the same general form (for a derivation, see, e.g., Blumenthal & Gould 1970) $\frac{\mathit{\partial N}}{\mathit{\partial t}}\mathrm{=}\mathit{Ṅ}\mathrm{syn}\mathrm{+}\mathit{Ṅ}\mathrm{cs}\mathrm{+}\mathit{Ṅ}\mathrm{pp}\mathrm{+}\mathit{Ṅ}\mathrm{ad}\mathrm{+}{\mathit{Q}}_{\mathrm{inj}}\mathrm{-}\frac{\mathit{N}}{{\mathit{t}}_{\mathrm{esc}}}\mathrm{-}\frac{\mathit{N}}{{\mathit{t}}_{\exp }}\mathit{,}$(26)where N ≡ N(x) (for photons) or N ≡ N(γ) (for electrons) is the particle number density distribution in the shocked region per dimensionless photon or electron energy interval. The subscripts syn, cs, pp and ad correspond to synchrotron emission and absorption, Compton scattering, pair production and adiabatic cooling, respectively. Each process produces a term on the right-hand side of the equation (for a detailed description of the numerical treatment of the processes, see Vurm & Poutanen 2009). The source term Q_inj gives the contribution of newly injected particles. The electron injection function is assumed to be a pure power law (Eq. (13)) in most examples presented in this paper. However, the injected electrons have a Maxwellian distribution in some of the examples in Sect. 4.3. The term N/t_esc accounts for the escape of particles from the emission region and N/t_exp gives the dilution of the particle densities due to the expansion of the emission region, t_esc and t_exp being the characteristic time scales for these processes. The adiabatic cooling and density dilution terms are discussed in Sect. 3.3. The photon escape time t_esc is evaluated by solving the plane-parallel radiative diffusion equation and is equal to ${\mathit{t}}_{\mathrm{esc}}\mathrm{=}\frac{\mathrm{3}\mathrm{\Delta }\mathit{R}}{\mathrm{4}\mathit{c}}\left[\mathrm{1}\mathrm{+}\frac{\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\ast }}}}{\sqrt{\mathit{ϵ}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{\ast }}}\mathrm{)}}\right]\mathit{,}$(27)where ΔR is the comoving radial dimension of the emission region, ${\mathit{\tau }}_{\mathrm{\ast }}\mathrm{\equiv }\mathrm{2}\mathrm{\Delta }\mathit{R}\sqrt{\mathrm{(}{\mathit{\alpha }}_{\mathrm{syn}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{pp}}\mathrm{)}\mathrm{(}{\mathit{\alpha }}_{\mathrm{syn}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{pp}}\mathrm{+}{\mathit{\alpha }}_{\mathrm{cs}}\mathrm{)}}$ gives the effective optical depth and ϵ ≡ (α_syn + α_pp)/(α_syn + α_pp + α_cs) is the probability for photon absorption, α_syn, α_pp and α_cs being the extinction coefficients due to synchrotron (syn) and pair production (pp) absorption and Compton scattering (cs).

The code recalculates the radius r and bulk Lorentz factor Γ at each timestep according to Eqs. (2) and (7) and evaluates the observed time from Eq. (8). The other time-dependent quantities such as B, ṅ_el and γ_min (Eqs. (12), (14) and (15)) can then be evaluated to obtain, for instance, the updated synchrotron emissivities.

The main difference between our code and the similar one described in PM09 is that we include the second-order differential terms that correspond to electron heating and diffusion due to synchrotron and Compton processes. The differential terms in the kinetic equation for electrons take the form $\mathit{Ṅ}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial \gamma }}\left[{\mathit{A}}_{\mathrm{e}}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{-}{\mathit{B}}_{\mathrm{e}}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\frac{\mathit{\partial N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}}{\mathit{\partial \gamma }}\right]$(28)for all continuous processes, such as synchrotron processes, Compton scattering in the Thomson regime and adiabatic cooling, and the coefficients A_e and B_e depend on the process of interest. Synchrotron self-absorption also contributes a first- and second-order differential term in the electron equation, where the coefficients A_e and B_e depend on the number density of photons and the synchrotron emissivity of an electron.

The kinetic equations are discretized and solved on finite grids of photon and electron/positron energies. Both the photon and electron energy grids consist of 200 points, with the photon grid ranging from x = 10^-11 to x = 10⁸ (E = 5 × 10^-6eV to E = 50 TeV) and the electron grid ranging from p = 10^-4 to p = 10⁸.

Because of the high-energy boundary of the electron grid, the maximum energy of the electrons is evaluated as ${\mathit{\gamma }}_{\max }\mathrm{=}\begin{array}{c}\min \end{array}\left({\mathrm{10}}^{\mathrm{8}}\mathit{,}\sqrt{\frac{\mathrm{3}\mathit{\pi e}}{{\mathit{\sigma }}_{\mathrm{T}}{\mathit{B}}^{\mathrm{\prime }}}}\right)\mathit{.}$(29)The maximum Lorentz factor ${\mathit{\gamma }}_{\max }\mathrm{=}\sqrt{\mathrm{3}\mathit{\pi e}\mathit{/}\mathrm{\left(}{\mathit{\sigma }}_{\mathrm{T}}{\mathit{B}}^{\mathrm{\prime }}\mathrm{\right)}}$ is obtained by comparing the acceleration time to the synchrotron cooling time. This is the value used by our code as long as it does not exceed γ = 10⁸, the highest electron energy of the grid.

3.2. Size of the emission region

In the simulations presented here, we consider only the radiation from the forward shock. The current code applies a one-zone approximation by assuming that the particles behind the shock front are homogeneously distributed and that the magnetic field has a uniform value in the region of interest.

The emission region is the spherical shell between the forward shock and the contact discontinuity that separates the shocked external medium and the shocked GRB ejecta. The Lorentz factor of the forward shock is ${\mathrm{\Gamma }}_{\mathrm{sh}}\mathrm{=}\sqrt{\mathrm{2}}\mathrm{\Gamma }$ (Blandford & McKee 1976), which means that the contact discontinuity is moving away from the shock at a velocity c/3. This leads us to define the radial extent of the shell as $\mathrm{\Delta }{\mathit{R}}^{\mathrm{\prime }}\mathrm{=}\mathit{c}{\mathit{t}}^{\mathrm{\prime }}\mathit{/}\mathrm{3}$(30)and the comoving volume of the emission region becomes ${\mathit{V}}^{\mathrm{\prime }}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathit{c}{\mathit{t}}^{\mathrm{\prime }}\mathit{/}\mathrm{3.}$(31)Here it is not taken into account that the shell is most likely not spherical but a fraction of a narrow jet. However, the luminosity per unit solid angle is the same in both of these cases before the so-called jet break time, when the beaming angle of the radiation exceeds the opening angle of the jet.

3.3. Adiabatic cooling and density dilution

In addition to the radiation processes discussed in this paper, the code accounts for adiabatic particle cooling due to the spreading emission region. In this section, all quantities except for r are given in the fluid frame and are left unprimed. The adiabatic cooling term for electrons is $\mathit{Ṅ}\mathrm{ad}\mathrm{\equiv }\frac{\mathit{\partial N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial \gamma }}\left[\mathit{\gamma ̇}\mathrm{ad}\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\right]\mathit{,}$(32)where is the cooling rate due to adiabatic expansion. The cooling rate is obtained from the pressure P and internal energy E of the electrons. For a monoenergetic population of N_e electrons of energy γ in a volume V, $\mathit{P}\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}}\frac{{\mathit{N}}_{\mathrm{e}}}{\mathit{V}}\mathit{\gamma }{\mathit{\beta }}_{\mathrm{e}}^{\mathrm{2}}\mathit{,}$(33)where β_e is the random velocity of an electron in units of c, and $\mathit{E}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{N}}_{\mathrm{e}}\mathit{.}$(34)From the first law of thermodynamics, dE = −PdV, one then obtains the cooling rate $\mathit{\gamma ̇}\mathrm{ad}\mathrm{=}\mathrm{-}\frac{\mathit{\gamma }{\mathit{\beta }}_{\mathrm{e}}^{\mathrm{2}}}{\mathrm{3}}\frac{\mathrm{d}\ln \mathit{V}}{\mathrm{d}\mathit{t}}\mathrm{·}$(35)Evaluating the time derivative of the comoving volume (Eq. (31)), the cooling rate for a constant value of g (see Eqs. (3)−(6)) becomes $\mathit{\gamma ̇}\mathrm{ad}\mathrm{=}\mathrm{-}\frac{\mathit{\gamma }{\mathit{\beta }}_{\mathrm{e}}^{\mathrm{2}}}{\mathrm{3}\mathit{t}}\left(\frac{\mathrm{2}}{\mathit{g}}\mathrm{+}\mathrm{1}\right)\mathit{,}$(36)where t is the comoving lifetime of the shock. Because the simulation extends from the coasting phase of the relativistic shell to the deceleration phase, the value changes from g = 1 to g = 5/2 (constant density) or g = 3 /2 (wind) according to Eqs. (3) and (5) during the deceleration. This change is gradual in reality, but our approximation of the Blandford-McKee solution has the side effect that the time derivative of the volume and the cooling term of Eq. (36) are discontinuous at r = R_B. In order to avoid the discontinuity, the cooling rate used in the simulations is $\mathit{\gamma ̇}\mathrm{ad}\mathrm{=}\mathrm{-}\frac{\mathit{\gamma }{\mathit{\beta }}_{\mathrm{e}}^{\mathrm{2}}}{\mathrm{3}}\left(\frac{\mathrm{2}}{{\mathit{t}}_{\mathrm{ad}}}\mathrm{+}\frac{\mathrm{1}}{\mathit{t}}\right)\mathit{,}$(37)where $\begin{array}{ccc}{\mathit{t}}_{\mathrm{ad}}\mathrm{\equiv }\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(38)t_B being the comoving time corresponding to the transition radius R_B. This approach guarantees that the adiabatic cooling term is continuous and that the solution is accurate both at early and late times.

The term giving the dilution of the electron density (the last term on the right-hand side of Eq. (26)) is obtained by considering a constant number of particles with density n = ∫N(γ)dγ in an expanding volume V: $\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}\mathrm{(}\mathit{n}{\mathit{V}}^{\mathrm{)}}\mathrm{=}\mathrm{0.}$(39)This relation is equivalent to $\frac{\mathrm{d}\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\frac{\mathrm{d}\ln \mathit{V}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}}{{\mathit{t}}_{\exp }}\mathit{,}$(40)where ${\mathit{t}}_{\exp }\mathrm{=}{\left(\frac{\mathrm{2}}{{\mathit{t}}_{\mathrm{ad}}}\mathrm{+}\frac{\mathrm{1}}{\mathit{t}}\right)}^{-1}$(41)and t_ad is given by Eq. (38).

3.4. Test simulations

Fig. 1 Simulated afterglow spectrum in the observer frame (top panel) and the electron distribution in the fluid frame (bottom panel) at radius r = 3.4 × 1017 cm. The particle distributions have been obtained with four different combinations of radiative and adiabatic processes. Synchrotron emission/absorption, Compton scattering, pair production/annihilation and adiabatic energy losses are indicated in the figure by the abbreviations syn, IC, pp and ad, respectively. In the simulation where only synchrotron processes are included, self-absorption heating of the electrons is not accounted for. The analytic synchrotron solution by Sari et al. (1998) is also shown in the figures. The simulation parameters are E0 = 1053erg, n0 = 1 cm-3, ϵe = 0.1, ϵB = 10-3, Γ0 = 400 and s = 2.3. The high-energy cutoff of the electron distribution has a constant value γmax = 4 × 107, and the GRB is located at a redshift z = 1. In the bottom panel, N(γ) ≡ dn/dγ is the number density of electrons per dimensionless energy interval. The distribution function has been multiplied by σTΔR′ to find the Thomson optical depth and by γs in order to visualize which electrons are cooling: a flat segment in the distribution means that the slope is the same as that of the injection function, demonstrating that the electrons are uncooled. The figure may be compared with Figs. 4 and 5 of PM09.

In order to test the validity of our code, we have compared our simulation results with those obtained by PM09 who have developed a similar numerical code. The photon spectra and electron distributions at r = 3.4 × 10¹⁷ cm (Fig. 1) calculated by our code should be compared with Figs. 4 and 5 of PM09, which present the particle distributions obtained with the same set of parameters. In these simulations, a constant value of the maximum electron energy, γ_max = 4 × 10⁷, is assumed. Our Fig. 1 shows the results of several simulations with different combinations of radiation processes. All the photon spectra in this section and the rest of the paper are presented in the observer frame, with energies boosted by a factor of 2Γ/(1 + z) from the flow frame.

The two codes produce electron distributions with very similar shapes: the cooled electrons populate the high-energy end of the distribution going as N(γ) ∝ γ^−s−1, with the slope of the distribution gradually approaching the slope of the injection function at lower energies. The electrons below γ_min ~ 600 have a nearly flat N(γ) distribution in PM09, whereas we obtain a much steeper decline of the electron density. To study whether this difference could be due to the inclusion of synchrotron self-absorption heating, we turned off this process for a test simulation, but this had no visible impact on the electron distribution. Because the shell is in the deceleration phase and γ_min is decreasing, such a steep cutoff may appear if γ_min declines faster than the electrons are cooled adiabatically. However, the electrons below γ_min do not carry a large fraction of the total electron energy and have no observable impact on the afterglow emission.

Fig. 2 Left panels: time-evolving observer frame photon spectrum without absorption on extragalactic background light (top panel) and electron distribution (bottom panel) from the forward shock together with an additional photon injection term that roughly represents prompt GRB emission. The electron distributions are plotted as a function of the dimensionless electron momentum p (Eq. (1)) instead of γ (note that p → γ for γ ≫ 1). The number of electrons per unit dimensionless momentum is N(p) ≡ dn/dp. The parameters of the forward shock emission are E0 = 1.4 × 1055 erg, n0 = 3 × 10-2 cm-3, ϵe = 0.2, ϵB = 10-6, Γ0 = 800 and s = 2.4. The injected photons are distributed according to a Band function, with the parameters α = −0.61, β = −3.8 and Epk = 730 keV (as defined in Eq. (24)). The Band function cuts off at E = 1 GeV and the injection lasts for t = 22 s, after which the photon luminosity decreases exponentially with time. The GRB takes place at z = 1.8. The black long-dashed lines show the photon and electron distributions at t = 0.1 s without pair production, and the solid magenta line in the bottom left panel represents the positron distribution at t = 0.1 s. Right panels: the evolution of the observer frame photon spectrum (top panel) and electron distribution (bottom panel) without photon injection. The forward shock parameters are same as in the simulation presented in the left panels.

The normalization of our electron distribution is lower by about half an order of magnitude than that obtained by PM09, but the origin of this difference is unclear so far. In our radiation spectrum, the spectral slopes and the positions of the peak and break frequencies appear identical to those in PM09, but the normalization of the flux is again slightly different. It is notable that the relative magnitudes of the synchrotron and inverse Compton components are clearly different in the two simulations. A possible cause for this discrepancy is the slight difference in the geometries assumed in the simulations.

The electron distribution resulting from synchrotron cooling without adiabatic cooling or self-absorption heating is also presented in Fig. 1, showing that especially the low-energy electrons are strongly affected by adiabatic cooling. The distribution around the cooling break is highly curved without this process, and there is a clear cutoff below the injection energy γ_min since the electrons no longer have a way to cool to lower energies.

The radiation spectrum from this simulation is also presented in Fig. 1. The curved part of the electron distribution produces a corresponding curved segment in the radiation spectrum. The low-energy end of the spectrum, which basically consists of radiation from monoenergetic electrons at γ_min, has the same shape in all cases. All in all, adiabatic cooling has a clear observable effect on the emergent spectrum but the contribution of synchrotron self-absorption heating seems negligible.

The analytic synchrotron solution by Sari et al. (1998) is included in Fig. 1 for comparison with the numerical solution. It is clear that a spectrum consisting of pure power-law segments with sharp breaks deviates from the actual curved synchrotron radiation component, and fitting the simple analytic model to observed afterglow spectra may lead to an inaccurate determination of the forward shock parameters.

4. Examples

4.1. High-energy emission due to Compton scattering of MeV photons

As an example of an important application of the code, we present the results of two simulations where the injected electrons interact with an external source of photons distributed according to a Band function (Eq. (24)), which represents a typical prompt GRB spectrum. Some of the GRB photons are scattered to higher energies by the shock-accelerated electrons, producing an additional spectral component extending up to TeV energies in the observer frame. The high-energy end of the spectrum can then be further modified as some of the high-energy photons produce electron-positron pairs with the MeV photons. The ≳TeV component may be observable in the case of a low-redshift GRB; otherwise the very high-energy gamma-rays are absorbed by the extragalactic background light. The simulations presented here demonstrate that our code can be applied to study whether the >100 MeV GRB emission is due to Compton scattering of the prompt photons at the external shocks. However, further modifications of the code are needed to improve especially the treatment of the early high-energy emission.

Fig. 3 Time-evolving observer frame photon spectrum without absorption on extragalactic background light (top panel) and electron distribution (bottom panel) from the forward shock with the same photon injection term as in Fig. 2. The other simulation parameters are E0 = 1.4 × 1055erg, n0 = 3 cm-3, ϵe = 0.8, ϵB = 10-8, Γ0 = 450 and s = 2.4. The redshift is z = 1.8. The black and red long-dashed lines show the photon and electrons distributions at t = 0.1 s and t = 10 s without pair production, and the positron distribution at t = 0.1 s is presented by a solid magenta line in the bottom panel.

The current version of the code calculates the emission from the forward shock only, although the reverse shock emission can be significant in the case of a long burst. In addition, electron-positron pair production in the external medium due to the interaction with the prompt gamma-rays is likely to affect the afterglow emission (Beloborodov 2005a). This effect is not accounted for in the code at the moment. However, both the reverse shock emission and the pair loading of the external medium are expected to shape the afterglow at relatively early times, and thus the late-time afterglow emission is likely to be fairly consistent with the results of the forward shock simulation.

Figure 2 shows the simulated time-evolving forward shock spectrum and the corresponding electron distribution both with and without an additional Band injection term, and Fig. 3 presents the results of a simulation with the same Band injection function but different external shock parameters. The positron distribution at t = 0.1 s is also presented in the cases where pair production is important, meaning in the simulations with a photon injection term. At the moment, the positron distribution has a lower numerical accuracy than that of the electrons.

The simulation parameters have been chosen to reproduce the Swift/XRT afterglow light curve of GRB 090902B, which was a bright burst with a long-lasting high-energy component (Abdo et al. 2009). The Band function used in the simulations is kept constant in time for simplicity and corresponds to the time-integrated spectral fit presented in Table 1 of Abdo et al. (2009). The Band photon density n_ph in the observer frame is evaluated from Eq. (25), where the luminosity is obtained from dividing the isotropic equivalent energy of the prompt GRB, E_iso = 3.6 × 10⁵⁴ erg, by the burst duration, t = 22 s. The number of injected photons per unit time and area is then ~cn_ph in the observer frame, and the injection rate per unit volume in the fluid frame is $\mathit{ṅ}\begin{array}{c}\mathrm{\prime }\\ \mathrm{ph}\end{array}\mathrm{~}\frac{\mathit{c}{\mathit{n}}_{\mathrm{ph}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}{\mathrm{2}\mathrm{\Gamma \Delta }{\mathit{R}}^{\mathrm{\prime }}}\mathit{,}$(42)where the shell thickness ΔR′ is defined in Eq. (30). We do not inject the additional power-law component that is observed in the prompt emission below ~50 keV and above 100 MeV.

The allowed forward shock parameter space based on the late-time afterglow data has been calculated by Kumar & Barniol Duran (2010), who claim that the >100 MeV emission can be explained as pure synchrotron radiation. The simulations presented in this section correspond to two points in the allowed n − ϵ_B space given in Figure 3 of their paper, with E₀ and ϵ_e being evaluated from their Eqs. (12) and (13). Our Fig. 2 shows the results in the case of a lower density and a higher magnetic energy fraction than in Fig. 3. We assume that the ejecta begin to decelerate at the end of the prompt burst. This assumption is used to evaluate the initial bulk Lorentz factor of the emitting shell, Γ₀. The extragalactic background light absorption of the ≳TeV photons is not accounted for in the example presented here, and these photons would in fact be unobservable from the redshift of GRB 090902B, z = 1.8. However, a very high-energy component could be observed in the case of a low-redshift burst with similar parameters.

The left-hand panels in Fig. 2 show that the Band emission at energies E ≲ 1 GeV dominates over the underlying forward shock synchrotron radiation at t = 0.1 s, while a double-peaked inverse Compton component appears at higher energies. The break at ~100 GeV is due to pair production with the peak Band photons, which are located at the threshold energy for pair production with the photons at ~100 GeV according to Eq. (22). For target photons with a Band distribution, the optical depth peaks at ~10x_thr, where the lowest point between the two VHE spectral peaks is located. The maximum optical depth can be approximated using Eqs. (23) and (25) to obtain τ_γγ ~ a few. The importance of pair production can be confirmed by looking at the electron and positron distributions, which are nearly identical below γ_min ~ 10⁵. The VHE peak in the radiation spectrum is due to the scattering of Band photons against the peak electrons with γ_pk = γ_min. According to Eq. (20), the peak energy is E ~ 30 TeV, which agrees with the figure.

The right-hand panels, showing the results of the simulation without the Band injection term, confirm that the synchrotron spectral component at t = 0.1 s is clearly below the Band emission, and the contribution from the Compton scattered synchrotron photons is also negligible. In the simulation with the additional photon injection term, a low-energy tail forms in the electron distribution. This is mainly caused by pair production with a contribution from scatterings in the K-N regime, which can take a large fraction of the energy of an electron in one scattering.

At t = 10 s, close to the deceleration time, the IC component has only one peak left in both simulations. The optical depth for pair production is expected to decrease until the deceleration time according to τ_γγ ∝ t^-1 and thus should drop by a factor of 100 between t = 0.1 s (where τ_γγ ~ 2) and t = 10 s, which is in agreement with the simulation results. In the left panels, the prompt emission dominates as the source of seed photons for inverse Compton, as the resulting high-energy emission is ~2 orders of magnitude more luminous than with synchrotron self-Compton (SSC) emission alone. The high-energy part of the synchrotron bump has become visible at t = 10 s and dominates the emission at E ~ 100 MeV−10 GeV. The inverse Compton peak is still located at E ~ 30 TeV, and it seems that the prediction for the location of the peak is very robust during the prompt emission due to the hard spectrum of the soft target photons. In the slow cooling regime and in the absence of significant pair opacity, the spectral slope below the ~30 TeV peak mimics the low-energy slope of the prompt spectrum. Because the peak energy is proportional to Γγ_min and γ_min ∝ ϵ_eΓ, a measurement of the peak energy would provide us with the combination ϵ_eΓ².

After the end of the photon injection, the spectrum quickly becomes identical to the pure forward shock emission, as can be seen by comparing the spectra at t = 10³ s and t = 10⁵ s in the left and right panels. The electrons are unable to cool by emitting synchrotron radiation because such a small fraction of the shock energy is given to the magnetic field. The SSC cooling of the electrons on synchrotron radiation is not effective either, which can be seen from the electrons being distributed as N(γ) ∝ γ^−s between γ_min and γ_max. The synchrotron spectral component produced by these electrons now corresponds to a slow cooling spectrum as described by Sari et al. (1998).

The simulation results presented in Fig. 3 correspond to a higher density and smaller magnetic energy density than those in Fig. 2. The VHE spectrum at t = 0.1 s is similar to the one in Fig. 2. In this case, the spectrum at t = 10 s is also shaped by the inclusion of pair production, which leads to an increased energy content in electrons at γ<γ_min and correspondingly more energy in the upscattered photons between ~3 GeV and ~1 TeV. After the prompt emission is gone, the SSC component is very prominent compared to the synchrotron bump because of the high value of ϵ_e/ϵ_B used in the simulation. Even though a large fraction of the shock-generated energy, ϵ_e = 0.8, is given to the electrons now, it is still valid to assume an adiabatic evolution of the shock because the radiative cooling is inefficient.

It is notable that a power-law component similar to the one observed during the prompt phase of GRB 090902B (Abdo et al. 2009) is not seen in our simulations. According to the simple model presented here, the early power-law component cannot be of external origin. The low-energy part of the component might be detectable above the Band emission in the case of a low γ_min (Eq. (15)). Because we assume that the deceleration time is approximately equal to the duration of the prompt GRB, we can use Eqs. (4) and (9) in our paper, together with Eqs. (12) and (13) in Kumar & Barniol Duran (2010) to find that ${\mathit{\gamma }}_{\min }\cong \mathrm{6}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\mathrm{(}\mathit{n}{\mathit{ϵ}}_{\mathit{B}}^{\mathrm{3}}{\mathrm{)}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{16}}$; i.e., γ_min depends very weakly on n and ϵ_B. From the allowed n − ϵ_B space presented in Fig. 3 of Kumar & Barniol Duran (2010), one finds that γ_min cannot be much lower than ~10⁵.

Fig. 4 Observer frame light curves in the Swift/XRT band (0.3−10 keV, red curves at the top of the figure), the Fermi/LAT band (0.1−300 GeV, green curves in the middle), and a very high energy (VHE) 0.1−100 TeV band (blue curves at the bottom). The light curves correspond to the simulations presented in Figs. 2 and 3.

Figure 4 shows the light curves corresponding to the simulations presented in Figs. 2 and 3 (Band injection is included in Examples 1 and 3, but not in Example 2) in three different energy bands: the XRT and LAT bands and a very high-energy (VHE) band ranging from 0.1 TeV to 100 TeV, where the only contribution is from the inverse Compton emission due to the upscattering of either synchrotron or external photons. Because the Band function is kept constant in time, the flux in both the XRT and LAT bands is also constant in the beginning of the simulations with the photon injection term. The inverse Compton flux in the VHE band rises at early times as the energy of the shocked electrons increases and the suppression of the flux due to pair production becomes less important.

The LAT light curves are dominated by synchrotron emission between t ~ 10 s and 10³ s. At 10³ s the synchrotron component falls below the LAT band and there is a break in the light curves, followed by a more slowly decaying SSC emission. The X-ray emission from the forward shock is buried under the prompt emission and emerges only after the latter has ended. In the TeV band the emission is dominated by the upscattered prompt photons as long as they are available, resulting in substantially higher luminosity at early times than would be obtained by SSC alone.

The simulated XRT band light curves are consistent with the observed XRT data of GRB 090902B (Pandey et al. 2010) as intended, but the observed LAT light curve (Abdo et al. 2009) is not reproduced in the current simulations. This is simply because we assume a lower value of γ_max (Eq. (29)) than Kumar & Barniol Duran (2010), who claim that the high-energy emission consists of synchrotron photons. The maximum photon energy produced by the synchrotron process is ~2Γ × 50 MeV/(1 + z) in the observer frame (Guilbert et al. 1983), determined from the balance between the cooling and acceleration times for electrons radiating in the background magnetic field B′ (for synchrotron losses taking place in a magnetic field generated by the Weibel instability, see Sironi et al. 2013). The most energetic photon observed from GRB 090902B, with E = 33 GeV in the observer frame, requires a bulk Lorentz factor Γ ~ 900 at t = 82 s to be consistent with the synchrotron model. The Lorentz factor at this time can be estimated from Eq. (11). For E₀ ~ 10⁵⁵erg and n₀ = 10^-3 cm^-3, the lowest possible external density according to Kumar & Barniol Duran (2010), we find that Γ ~ 400 at t = 82 s, not supporting a synchrotron origin for this photon.

The examples in this section show that a luminous ≳TeV component can arise due to inverse Compton scattering in the forward shock, although this radiation is observable only in low-redshift cases. The Cherenkov telescope MAGIC is the most likely candidate to catch the upscattered prompt radiation in the VHE range due to its fast slewing capability, but MAGIC has so far been unable to detect any emission from GRBs (Albert et al. 2007; Aleksić et al. 2010, 2014). However, the lack of detections does not rule out the existence of a ~10 TeV spectral component because of the high redshifts of the observed bursts.

Fig. 5 Left panels: time-evolving observer frame photon spectrum (top panel) and electron distribution (bottom panel) resulting from synchrotron emission and absorption together with adiabatic cooling at the forward shock in a wind-type medium. The simulation parameters are E0 = 1054erg, Aw,35 ≡ Aw/(1035 cm-1) = 3, ϵe = 0.1, ϵB = 0.1, Γ0 = 500 and s = 2.5. The redshift is z = 1. Right panels: observed photon spectra (top panel) and electron distributions (bottom panel) in the case of a constant-density medium with a typical density n0 = 1 cm-3. The other simulation parameters are the same as in the simulation presented in the left panels.

If the external density were higher than in our examples, such as in the case of a wind-type medium, the inverse Compton flux would also increase and possibly be able to reproduce the magnitude of the GeV photon flux in the LAT band. This idea is supported by the recent simulations of Beloborodov et al. (2013), who show that the GeV emission of GRB 080916C is consistent with inverse Compton scattering of the prompt emission. The work of Beloborodov et al. (2013) is mainly concerned with the early stage of the outflow and does not include SSC emission, which is important once the prompt radiation is gone.

Our current model cannot explain the extended duration of the observed LAT emission from GRB 090902B, which is however reproduced by Beloborodov et al. (2013). The delay of the arrival times of the high-energy photons compared to the MeV emission may result from the different propagation angle of the Compton scattered prompt photons compared to that of the unscattered photons, which is not accounted for in the code at the moment. It must be noted that a full treatment of the afterglow requires a more complete model of the prompt emission, the inclusion of the reverse shock emission and the effect of the pair loading of the external medium.

4.2. Constant-density ISM vs. wind medium

The density structure of the ambient medium has an impact on the evolution of the bulk Lorentz factor Γ, the injection rate of the electrons and the strength of the shock-compressed magnetic field. In a wind-type medium with n(r) = A_wr^-2 all the synchrotron-emitting electrons are likely to be fast-cooling at early times because of the high magnetic field, provided that the magnetization parameter ϵ_B is not too low. Synchrotron cooling is counteracted by self-absorption at the low-energy end of the distribution. As long as the synchrotron cooling time of the electrons is shorter than the lifetime of the flow for all electron energies, the low-energy electrons are able to thermalize due to self-absorption.

The difference between the pure synchrotron emission from the forward shock in a wind-type and a constant-density medium is illustrated in Fig. 5, which shows the simulated photon and electron distributions for a typical set of parameters. Compton scattering and pair production are left out of this simulation in order to study the effect of synchrotron self-absorption heating on the particle distributions. Again, we assume that the distribution of the shock-accelerated electrons is a pure power law.

The deceleration radius (Eqs. (4) and (6)) is reached earlier in a wind medium than a constant-density ISM. In the simulations presented in Fig. 5, R_dec = 7.1 × 10¹⁴ cm in the wind case and R_dec = 8.6 × 10¹⁶ cm in the ISM case. Because most of the blast energy is dissipated at this radius, the peak luminosity is higher in the wind case because the same amount of energy is emitted in a shorter time compared to the ISM environment.

At the observed time t = 0.1 s, the ambient particle density in the wind case is n = 1.1 × 10⁶ cm^-3, explaining the large difference between the number of shocked electrons in the two simulations. At this moment, all the electrons are fast-cooling in both simulations due to the large fraction of energy given to the magnetic field. However, the magnetic field in the wind case is much higher because of the high density, even though the shock has already entered the deceleration phase before t = 0.1 s. The electrons have now been able to cool below γ_min and form a clear Maxwellian component centered at γ ~ 3. The Maxwellian electrons have had plenty of time to thermalize because their synchrotron cooling time in the fluid comoving frame, ${\mathit{t}}_{\mathrm{cool}}^{\mathrm{\prime }}\mathrm{~}\mathrm{0.2}\mathrm{s}$ (Eqs. (12) and (19)) is clearly shorter than the comoving lifetime of the flow, t′ = 40 s. The electrons at higher energies are distributed according to N(γ) ∝ γ^-2 below γ_min and N(γ) ∝ γ^−s−1 between γ_min and γ_max, as expected. It is notable that the electrons would be able to cool to much lower energies in a simulation without the self-absorption heating term. Thus, the inclusion of this term is necessary for the modeling of a wind-type medium with a high magnetic compactness. The thermalized electrons, however, only produce a very small bump in the corresponding synchrotron spectrum.

At t = 10 s, the Maxwellian component in the electron distribution in a wind medium is already smaller because of the decreasing magnetic field. The bump gradually disappears, and the slope of the electron distribution at t = 10⁵ s between γ_c ~ 600 and γ_min ~ 3 × 10³ is becoming slightly shallower because these electrons are no longer able to cool. At this moment, the electron and photon distributions in the ISM case look very similar to those in the wind case.

Fig. 6 Comparison of the radiation spectrum (top panel) and electron distribution (bottom panel) at t = 0.1 s from Fig. 5 (black solid line) to four cases where the value of one parameter is changed at a time. In two simulations, the wind density is decreased by a factor of 10 (red dashed line) and 100 (blue dashed line) from the fiducial value. In the other two cases, the magnetic energy fraction is decreased by a factor of 102 (red dash-dot line) and 104 (blue dash-dot line).

Fig. 7 Time evolution of the observer frame photon spectrum (top panel) and electron distribution (bottom panel) in a wind medium when all the relevant radiative processes (synchrotron, Compton and pair production) are accounted for in the simulation. The positron distribution at t = 0.1 s is also shown as a solid magenta line in the bottom panel. The simulation parameters are the same as in Fig. 5.

The simulations presented here show that the synchrotron self-absorption heating of electrons can produce a prominent Maxwellian component in the electron distribution for certain parameter sets. Figure 6 shows how the size of the Maxwellian bump at t = 0.1 s decreases if the parameters ϵ_B and A_w are changed from the fiducial values used to obtain Fig. 5 and how this affects the corresponding synchrotron spectrum. When the magnetic energy fraction ϵ_B is decreased, the synchrotron cooling time is longer and the electrons are unable to cool to very low energies or thermalize efficiently.

The magnetic field and the size of the thermalized bump also decrease if the density of the wind goes down: when the coefficient A_w,35 ≡ A_w/(10³⁵ cm^-1) is decreased from A_w,35 = 3 by a factor of 100, the low-energy bump in both the electron and photon distributions has nearly disappeared. It should be noted that the bulk Lorentz factor Γ is lower in the simulation with A_w,35 = 3 compared to the two cases with different A_w,35, because the shock is already in its deceleration phase at t = 0.1 s with Γ ~ 300 in the case with the highest density. As a result, the shock radius is slightly smaller in the case of a decreasing Γ, but only by a factor of ~ a few.

In general, the magnetic field (see Eq. (12)) decreases at radii r<R_dec (Eq. (6)) due to the radial dependence of the wind density: ${\mathit{B}}^{\mathrm{\prime }}\mathrm{\propto }\mathrm{\Gamma }{\mathit{n}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\propto }{\mathrm{\Gamma }}_{\mathrm{0}}{\mathit{r}}^{-1}\mathrm{\propto }{\mathrm{\Gamma }}_{\mathrm{0}}^{-1}{\mathit{t}}^{-1}$, so a larger Γ₀ indicates a lower magnetic field B′ at a given observer time t as long as the blast is its coasting phase. The synchrotron cooling time of the electrons thus becomes longer, working against the formation of a Maxwellian component.

The deceleration radius R_dec depends on the parameters E₀,A_w and Γ₀ but the hydrodynamic evolution at r > R_dec is only affected by E₀ and A_w according to Eqs. (10) and (11). The magnetic field then depends on these parameters and time as ${\mathit{B}}^{\mathrm{\prime }}\mathrm{\propto }\mathrm{\Gamma }{\mathit{n}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{\propto }\mathrm{(}{\mathit{A}}_{\mathrm{w}}^{\mathrm{3}}{\mathit{E}}_{\mathrm{0}}^{-1}{\mathit{t}}^{-3}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}$ (naturally, B′ also depends on ϵ_B). Together with the width of the emission region, ΔR′ ∝ Γt (Eq. (30)), the magnetic field determines the magnetic compactness l_B ≡ ΔR′U_Bσ_T/(m_ec²), where U_B = (B′)²/(8π) is the magnetic energy density. An electron with a random Lorentz factor γ is able to cool due to synchrotron emission if the ratio t/t_cool ~ γl_B (see Eq. (19)) is large. Because ${\mathit{l}}_{\mathrm{B}}\mathrm{\propto }\mathrm{\Delta }{\mathit{R}}^{\mathrm{\prime }}\mathrm{(}{\mathit{B}}^{\mathrm{\prime }}{\mathrm{)}}^{\mathrm{2}}\mathrm{\propto }{\mathit{A}}_{\mathrm{w}}^{\mathrm{5}\mathit{/}\mathrm{4}}{\mathit{E}}_{\mathrm{0}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{4}}{\mathit{t}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{4}}$, one can obtain lower values of t/t_cool at a given time by increasing E₀ or decreasing A_w, which may prevent the thermalization at low energies.

In reality, pair production has a significant effect on the external shock emission as long as the wind density is high. Figure 7 presents the time-evolving radiation spectrum and electron distribution from a simulation where synchrotron, Compton and pair production processes are all accounted for. At t = 0.1 s, the interplay of these processes results in a complex shape of the electron and positron distributions for γ< 100. Numerous leptons appear at these energies because of pair production, and are able to thermalize at low energies due to the high magnetic compactness. However, the Maxwellian leptons are now non-relativistic and do not result in an observable signature in the corresponding photon spectrum. The spectrum at t = 0.1 s is clearly different from the pure synchrotron case in Fig. 5. In Fig. 7, the electrons at γ< 100 produce a distinct spectral segment between ~100 eV and 100 keV, and pair production shapes the synchrotron and Compton components at high energies. The pair production opacity decreases at later times and the spectrum becomes similar to the pure synchrotron case, with the obvious exception that it extends to higher energies because of Compton scattering.

4.3. Emission from a Maxwellian distribution of injected electrons

Because a Maxwellian component can contain as much as ~90% of the energy of the shocked electrons (Spitkovsky 2008), we have also simulated the forward shock emission resulting from a Maxwellian electron injection function. The relativistic Maxwellian has the general form $\mathit{N}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{\propto }\mathit{\gamma }\sqrt{{\mathit{\gamma }}^{\mathrm{2}}\mathrm{-}\mathrm{1}}\exp \left[\mathrm{-}\frac{\mathit{\gamma }}{\mathit{k}{\mathit{T}}^{\mathrm{\prime }}\mathit{/}\left({\mathit{m}}_{\mathrm{e}}{\mathit{c}}^{\mathrm{2}}\right)}\right]\mathit{,}$(43)where T′ is the temperature of the electrons. The temperature can be determined from $\mathrm{3}\mathit{k}{\mathit{T}}^{\mathrm{\prime }}\mathrm{~}{\mathit{ϵ}}_{\mathrm{e}}\mathrm{\Gamma }{\mathit{m}}_{\mathrm{p}}{\mathit{c}}^{\mathrm{2}}\mathit{,}$(44)where ϵ_eΓm_pc² is the average energy transferred to a shocked electron.

Fig. 8 Left panels: time evolution of the observer frame photon spectrum (top panel) and electron distribution (bottom panel) resulting from relativistic Maxwellian electron injection in a constant-density ISM. The radiation processes are synchrotron emission/absorption, Compton scattering and pair production/annihilation. The simulation parameters are E0 = 1052erg, n0 = 1 cm-3, ϵe = 10-2, ϵB = 10-4 and Γ0 = 200, and the redshift is z = 1. The Maxwellian injection function has the temperature T′ ~ ϵeΓmpc2/k. Right panels: evolving particle distributions in the case of a power-law electron injection with s = 2.5. The other parameters are the same as in the simulation with Maxwellian injection.

A comparison between the simulated photon and electron distributions resulting from Maxwellian and power-law electron injection in a constant-density ISM is presented in Fig. 8. Synchrotron, Compton and pair production processes are all included in these examples. The main difference between the two cases is the width of the synchrotron component, which extends to higher energies when a large fraction of the electron energy is in the power-law electrons above the peak Lorentz factor.

In both simulations, the low-energy part of the synchrotron spectrum at t = 0.1 s and t = 10 s consists of the tail emission from the peak electrons, having the spectral shape νF_ν ∝ ν^4/3. The synchrotron spectral component is naturally broader in the case of power-law electron injection, since the radiating electrons in the case of Maxwellian injection are nearly monoenergetic. The ~10 TeV inverse Compton peak is due to the scattering events by the electrons near the K-N limit. At t = 10⁵ s, a second Compton scattering order appears as the first scattering order moves to ~100 eV. At all times, a Maxwellian electron injection results in a spectrum consisting of distinctly separate components, which is very different from the case of power-law electron injection.

Fig. 9 Observer frame light curves in the Swift/XRT band (0.3−10 keV) resulting from different parameter sets. The black solid line in both panels shows the fiducial light curve resulting from Maxwellian electron injection obtained with the same parameters as the spectra in Fig. 8. The green dash-dot-dot-dot line in the top panel represents the emission due to a power-law electron injection function. The other curves represent simulations with a Maxwellian injection term where one parameter at a time is varied from the fiducial value. In the top panel, the dotted (dashed) lines show the light curves when the blast energy E0 (external density n0) is increased/decreased by a factor of 100 (10). In the bottom panel, the long-dashed (dash-dot) lines correspond to the emission when the electron energy fraction ϵe (magnetic energy fraction ϵB) is increased/decreased by a factor of 10 (100).

Fig. 10 Evolving photon (top panel) and electron (bottom panel) distributions from a simulation with Maxwellian electron injection in a wind-type medium. The magenta line in the bottom panel shows the positron distribution at t = 0.1 s. The simulation parameters are the same as in Fig. 8, except for ϵe = 0.1 and Aw,35 = 3.

Fig. 11 Observer frame Swift/XRT band light curves in the case of Maxwellian electron injection in a wind medium. The solid line is the fiducial light curve corresponding to Fig. 10. The dotted (dash-dot) curve is obtained by increasing E0 (ϵB) by a factor of 100 from the fiducial value.

Figure 9 shows the light curves in the Swift/XRT band resulting from both the simulations in Fig. 8 compared to several cases where one parameter at a time is changed. It is notable that a flattening or a slight rebrightening appears in all the cases with a Maxwellian electron injection. The flat segment is produced when the synchrotron component moves out of the observed energy range and is gradually replaced by the inverse Compton emission, similarly as described by Petropoulou et al. (2011), who simulated the external shock emission from power-law electrons with a narrow energy distribution. A similar effect is also seen in the simulations of Stern & Poutanen (2004). They considered prompt GRB emission due to a narrow electron distribution that has developed due to the balance between heating and synchrotron/SSC cooling. In this case, a flat section in the BATSE band light curve develops when the first Compton scattering order moves out of the observed band and is replaced by the second scattering order. A spectrum due to power-law electrons does not contain such dramatic changes of the spectral slope, and consequently there is no shallow decay segment in the light curve.

In all the cases shown in Fig. 9, the flux during the flattening or rebrightening is very low. The flux is higher when E₀, n₀, ϵ_e or ϵ_B is increased, but higher values of E₀,ϵ_e and ϵ_B also lead to a delayed start of the flat segment. One can estimate that the flattening begins when the synchrotron frequency of the peak Maxwellian electrons decreases below the observed waveband. The peak electrons have γ ~ Γϵ_em_p/m_e (Eq. (44)) and the characteristic synchrotron energy of the emitted photons is ${\mathit{x}}_{\mathrm{syn}}\mathrm{\propto }\mathrm{\Gamma }{\mathit{\gamma }}^{\mathrm{2}}{\mathit{B}}^{\mathrm{\prime }}\mathrm{\propto }{\mathit{E}}_{\mathrm{0}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{ϵ}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{ϵ}}_{\mathit{B}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{t}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}$ (see Eqs. (16), (11) and (12)). Increased values of E₀,ϵ_e and ϵ_B thus correspond to a higher synchrotron energy at a given time t, delaying the moment when the synchrotron bump crosses the observing window. However, a higher external density n₀ does not significantly delay the appearance of the shallow decay.

A more luminous shallow decay phase may occur in a wind-type medium. The particle distributions resulting from a Maxwellian electron injection function in a wind-type medium are presented in Fig. 10. Similarly as in the case of power-law injection in a wind environment (Fig. 7), pair production shapes the electron distribution for γ< 10³ at t = 0.1 s, where it is practically identical to the positron distribution. However, we are now interested in the spectral evolution at later times when the shallow decay may appear as the observing window becomes dominated by inverse Compton photons. Figure 11 shows that the flux during the shallow decay can be high enough to be detected in the wind case. In the simulation presented here, the flux decay is steeper during this phase than in the ISM case, but still consistent with many of the observed XRT afterglows. The flux increases with higher values of, say, E₀ and ϵ_B, but simultaneously the start of the shallow decay is delayed, as described above.

The simulations with different parameter sets show that a flattening in the X-ray light curve is likely to take place if the electron injection function is Maxwellian. In the ISM case, the flux during the simulated flattening appears to be much lower than a typical shallow decay flux observed by XRT. A wind environment, on the other hand, can lead to a detectable photon flux during the shallow decay.

5. Summary and conclusions

In this paper, we have studied the modeling of GRB forward shock emission for the first time with a code that treats all the radiation processes self-consistently. The main advantages of our code are the inclusion of synchrotron self-absorption heating and Compton scattering as mechanisms for electron thermalization, and the exact treatment of Compton scattering and pair production. We have shown that our results are in very good agreement with both analytic estimations and numerical simulations that have been presented in the literature. Inclusion of an accurate calculation of radiative transport, a multi-zone treatment of the emitting region, and the emission from the reverse shock will be the topics of future research.

We presented the results of test simulations where external MeV photons, distributed similarly to prompt GRB emission, interact with the shock-accelerated electrons and are Compton scattered to TeV energies. Such VHE emission is likely to dominate the synchrotron self-Compton emission at times up to t ≳ 10 s in GRBs with similar parameter sets and might be an interesting target for observations of low-redshift bursts with future instruments such as the Cherenkov Telescope Array (Inoue et al. 2013). We also showed that the inverse Compton flux is not always quenched by electron-positron pair production, owing to the decrease in the pair production opacity before the deceleration time. A very prominent inverse Compton peak appears during the prompt emission, the location of which is determined by the maximum energy given to the photons by the peak electrons.

Even though the test simulations show that the inverse Compton scattering of photons distributed according to a Band function is likely to produce a high-energy spectral component, the current model does not explain the typical long duration of the observed LAT emission in the case of GRB 090902B. A possible explanation for the LAT data is that the unscattered prompt photons are traveling within a smaller solid angle than the scattered photons, and the arrival time of the upscattered radiation is delayed as a result (Beloborodov 2005b; Beloborodov et al. 2013). This effect is not accounted for in the current version of the code. However, the large-angle effect cannot explain the LAT emission if it lasts much longer than the prompt GRB.

As another example, we compared the synchrotron emission in a wind-type and a constant-density ISM medium for a set of typical parameters and studied the effect of different simulation parameters on the thermalization of low-energy electrons. As long as the synchrotron cooling time is shorter than the lifetime of the flow, the low-energy end of a power-law electron distribution is strongly thermalized due to synchrotron self-absorption heating, which is usually neglected in similar numerical codes. The shape of the emerging radiation spectrum is still nearly identical to the case with pure power-law electrons, although a small bump is seen in the synchrotron spectrum thanks to the Maxwellian shape of the low-energy electrons. However, the simulation shows that it is necessary to include the self-absorption heating term in the electron kinetic equation to prevent the electrons from cooling down to the low-energy edge of the grid as long as thermalization is important.

Additionally, we have demonstrated that a flattening or rebrightening segment in the forward shock light curve is expected in the case of purely Maxwellian electron injection in a waveband where the observed flux gradually becomes dominated by the inverse Compton component, as described by Petropoulou et al. (2011). Such flattenings are observed in a large fraction of X-ray light curves, but their physical origin is still under debate. In the ISM simulations, the flux during the flat segments is much lower than in the case of a typical shallow decay phase observed by Swift/XRT. However, the shallow decay flux in the wind case could be detectable.

Acknowledgments

T.P. received the funding for this work from the Finnish Doctoral Program in Astronomy and Space Physics.

References

All Figures

Fig. 1 Simulated afterglow spectrum in the observer frame (top panel) and the electron distribution in the fluid frame (bottom panel) at radius r = 3.4 × 1017 cm. The particle distributions have been obtained with four different combinations of radiative and adiabatic processes. Synchrotron emission/absorption, Compton scattering, pair production/annihilation and adiabatic energy losses are indicated in the figure by the abbreviations syn, IC, pp and ad, respectively. In the simulation where only synchrotron processes are included, self-absorption heating of the electrons is not accounted for. The analytic synchrotron solution by Sari et al. (1998) is also shown in the figures. The simulation parameters are E0 = 1053erg, n0 = 1 cm-3, ϵe = 0.1, ϵB = 10-3, Γ0 = 400 and s = 2.3. The high-energy cutoff of the electron distribution has a constant value γmax = 4 × 107, and the GRB is located at a redshift z = 1. In the bottom panel, N(γ) ≡ dn/dγ is the number density of electrons per dimensionless energy interval. The distribution function has been multiplied by σTΔR′ to find the Thomson optical depth and by γs in order to visualize which electrons are cooling: a flat segment in the distribution means that the slope is the same as that of the injection function, demonstrating that the electrons are uncooled. The figure may be compared with Figs. 4 and 5 of PM09.

In the text

Fig. 2 Left panels: time-evolving observer frame photon spectrum without absorption on extragalactic background light (top panel) and electron distribution (bottom panel) from the forward shock together with an additional photon injection term that roughly represents prompt GRB emission. The electron distributions are plotted as a function of the dimensionless electron momentum p (Eq. (1)) instead of γ (note that p → γ for γ ≫ 1). The number of electrons per unit dimensionless momentum is N(p) ≡ dn/dp. The parameters of the forward shock emission are E0 = 1.4 × 1055 erg, n0 = 3 × 10-2 cm-3, ϵe = 0.2, ϵB = 10-6, Γ0 = 800 and s = 2.4. The injected photons are distributed according to a Band function, with the parameters α = −0.61, β = −3.8 and Epk = 730 keV (as defined in Eq. (24)). The Band function cuts off at E = 1 GeV and the injection lasts for t = 22 s, after which the photon luminosity decreases exponentially with time. The GRB takes place at z = 1.8. The black long-dashed lines show the photon and electron distributions at t = 0.1 s without pair production, and the solid magenta line in the bottom left panel represents the positron distribution at t = 0.1 s. Right panels: the evolution of the observer frame photon spectrum (top panel) and electron distribution (bottom panel) without photon injection. The forward shock parameters are same as in the simulation presented in the left panels.

In the text

Fig. 3 Time-evolving observer frame photon spectrum without absorption on extragalactic background light (top panel) and electron distribution (bottom panel) from the forward shock with the same photon injection term as in Fig. 2. The other simulation parameters are E0 = 1.4 × 1055erg, n0 = 3 cm-3, ϵe = 0.8, ϵB = 10-8, Γ0 = 450 and s = 2.4. The redshift is z = 1.8. The black and red long-dashed lines show the photon and electrons distributions at t = 0.1 s and t = 10 s without pair production, and the positron distribution at t = 0.1 s is presented by a solid magenta line in the bottom panel.

In the text

	Fig. 4 Observer frame light curves in the Swift/XRT band (0.3−10 keV, red curves at the top of the figure), the Fermi/LAT band (0.1−300 GeV, green curves in the middle), and a very high energy (VHE) 0.1−100 TeV band (blue curves at the bottom). The light curves correspond to the simulations presented in Figs. 2 and 3.
In the text

Fig. 5 Left panels: time-evolving observer frame photon spectrum (top panel) and electron distribution (bottom panel) resulting from synchrotron emission and absorption together with adiabatic cooling at the forward shock in a wind-type medium. The simulation parameters are E0 = 1054erg, Aw,35 ≡ Aw/(1035 cm-1) = 3, ϵe = 0.1, ϵB = 0.1, Γ0 = 500 and s = 2.5. The redshift is z = 1. Right panels: observed photon spectra (top panel) and electron distributions (bottom panel) in the case of a constant-density medium with a typical density n0 = 1 cm-3. The other simulation parameters are the same as in the simulation presented in the left panels.

In the text

Fig. 6 Comparison of the radiation spectrum (top panel) and electron distribution (bottom panel) at t = 0.1 s from Fig. 5 (black solid line) to four cases where the value of one parameter is changed at a time. In two simulations, the wind density is decreased by a factor of 10 (red dashed line) and 100 (blue dashed line) from the fiducial value. In the other two cases, the magnetic energy fraction is decreased by a factor of 102 (red dash-dot line) and 104 (blue dash-dot line).

In the text

Fig. 7 Time evolution of the observer frame photon spectrum (top panel) and electron distribution (bottom panel) in a wind medium when all the relevant radiative processes (synchrotron, Compton and pair production) are accounted for in the simulation. The positron distribution at t = 0.1 s is also shown as a solid magenta line in the bottom panel. The simulation parameters are the same as in Fig. 5.

In the text

Fig. 8 Left panels: time evolution of the observer frame photon spectrum (top panel) and electron distribution (bottom panel) resulting from relativistic Maxwellian electron injection in a constant-density ISM. The radiation processes are synchrotron emission/absorption, Compton scattering and pair production/annihilation. The simulation parameters are E0 = 1052erg, n0 = 1 cm-3, ϵe = 10-2, ϵB = 10-4 and Γ0 = 200, and the redshift is z = 1. The Maxwellian injection function has the temperature T′ ~ ϵeΓmpc2/k. Right panels: evolving particle distributions in the case of a power-law electron injection with s = 2.5. The other parameters are the same as in the simulation with Maxwellian injection.

In the text

Fig. 9 Observer frame light curves in the Swift/XRT band (0.3−10 keV) resulting from different parameter sets. The black solid line in both panels shows the fiducial light curve resulting from Maxwellian electron injection obtained with the same parameters as the spectra in Fig. 8. The green dash-dot-dot-dot line in the top panel represents the emission due to a power-law electron injection function. The other curves represent simulations with a Maxwellian injection term where one parameter at a time is varied from the fiducial value. In the top panel, the dotted (dashed) lines show the light curves when the blast energy E0 (external density n0) is increased/decreased by a factor of 100 (10). In the bottom panel, the long-dashed (dash-dot) lines correspond to the emission when the electron energy fraction ϵe (magnetic energy fraction ϵB) is increased/decreased by a factor of 10 (100).

In the text

	Fig. 10 Evolving photon (top panel) and electron (bottom panel) distributions from a simulation with Maxwellian electron injection in a wind-type medium. The magenta line in the bottom panel shows the positron distribution at t = 0.1 s. The simulation parameters are the same as in Fig. 8, except for ϵe = 0.1 and Aw,35 = 3.
In the text

	Fig. 11 Observer frame Swift/XRT band light curves in the case of Maxwellian electron injection in a wind medium. The solid line is the fiducial light curve corresponding to Fig. 10. The dotted (dash-dot) curve is obtained by increasing E0 (ϵB) by a factor of 100 from the fiducial value.
In the text