Temporal Evolution of the Gamma-ray Burst Afterglow Spectrum for an Observer: GeV--TeV Synchrotron Self-Compton Light Curve

TL;DR

This paper presents a numerical simulation of gamma-ray burst afterglow emission, revealing earlier onset times and smoother light curves than previous models, with implications for understanding high-energy photon detections.

Contribution

The study introduces a comprehensive time-dependent simulation method for GRB afterglows, improving accuracy over analytical models and explaining observed high-energy photon features.

Findings

Earlier afterglow onset time than analytical estimates

Gamma-ray light curves can be smooth power-laws despite emission mechanism switching

Discussion on model parameter uncertainties and spectral-decay relations

Abstract

We numerically simulate the gamma-ray burst (GRB) afterglow emission with a one-zone time-dependent code. The temporal evolutions of the decelerating shocked shell and energy distributions of electrons and photons are consistently calculated. The photon spectrum and light curves for an observer are obtained taking into account the relativistic propagation of the shocked shell and the curvature of the emission surface. We find that the onset time of the afterglow is significantly earlier than the previous analytical estimate. The analytical formulae of the shock propagation and light curve for the radiative case are also different from our results. Our results show that even if the emission mechanism is switching from synchrotron to synchrotron self-Compton, the gamma-ray light curves can be a smooth power-law, which agrees with the observed light curve and the late detection of a 32 GeV photon in GRB 130427A. The uncertainty of the model parameters obtained with the analytical formula is discussed, especially in connection with the closure relation between spectral index and decay index.

Tables

Table 1. Table 1 Model parameters for the afterglow emission of GRB 130427A.

	$E_0$ (erg)	${\Gamma }_0$	${ϵ}_e$	${ϵ}_B$	$n_0$	$p$	DF111Discrepancy factor defined by Eq. (51).
Model A	$2.0\times {10}^{55}$	$350$	$0.03$	${10}^{-6}$	$1.0$	$2.35$	0.66
Model B	$4.5\times {10}^{55}$	$400$	$0.015$	${10}^{-5}$	$1.0$	$2.60$	1.6
Model C	$4.5\times {10}^{54}$	$400$	$0.23$	${10}^{-5}$	$0.1$	$2.60$	2.7

Equations

$$\Gamma =\frac{K_3(1/x)}{K_2(1/x)}-x,$$

$$\hat{\gamma }=\frac{x}{\Gamma -1}+1.$$

$${\Gamma }_{\mathrm{sh}}=\sqrt{\frac{(\Gamma +1)(\hat{\gamma }(\Gamma -1)+1{)}^2}{\hat{\gamma }(2-\hat{\gamma })(\Gamma -1)+2}}.$$

$$\frac{dM}{dt}=\frac{1}{\Gamma }\frac{dM}{d{t}^{\prime }}=4\pi R^{2}c{\beta }_{\mathrm{sh}}{n}_{\mathrm{ISM}}{m}_{\mathrm{p}},$$

$$\frac{d{E}_{\mathrm{sh}}^{\prime }}{d{t}^{\prime }}=\Gamma c^2\frac{dM}{d{t}^{\prime }}-\frac{d{E}_{\mathrm{rad}}^{\prime }}{d{t}^{\prime }}-\frac{d{E}_{\mathrm{ad}}^{\prime }}{d{t}^{\prime }},$$

$$E_{\mathrm{sh}}=\Gamma E_{\mathrm{sh}}^{\prime }=E_0+M{c}^2-E_{\mathrm{rad}}.$$

$$n^{\prime }=\frac{\hat{\gamma }\Gamma +1}{\hat{\gamma }-1}{n}_{\mathrm{ISM}},$$

$$\frac{d{E}_B^{\prime }}{d{t}^{\prime }}={ϵ}_B(\Gamma -1)c^2\frac{dM}{d{t}^{\prime }}.$$

$$B^{\prime }=\sqrt{\frac{8\pi E_B^{\prime }}{V^{\prime }}}.$$

$${\dot{N^{\prime }}}_{\mathrm{inj}}({\epsilon }_{\mathrm{e}}^{\prime })=\dot{N_0^{\prime }}({\epsilon }_{\mathrm{e}}^{\prime }/{\epsilon }_{\min }^{\prime }{)}^{-p}\exp (-{\epsilon }_{\mathrm{e}}^{\prime }/{\epsilon }_{\max }^{\prime }),$$

$$\frac{d{N}_{\mathrm{e}}^{\prime }}{d{t}^{\prime }}$$

$$\frac{d{E}_{\mathrm{e}}^{\prime }}{d{t}^{\prime }}$$

$$\displaystyle\frac{\partial N_{\rm e}(\varepsilon_{\rm e})}{\partial t}=\frac{\partial}{\partial\varepsilon_{\rm e}}\Bigl{[}\Bigl{(}\left<\dot{\varepsilon_{\rm e}}\right>_{\rm syn}+\left<\dot{\varepsilon_{\rm e}}\right>_{\rm IC}+\left<\dot{\varepsilon_{\rm e}}\right>_{\rm ad}$$

$$\displaystyle-\left<\dot{\varepsilon_{\rm e}}\right>_{\rm SSA}\Bigr{)}N_{\rm e}(\varepsilon_{\rm e})\Bigr{]}+\dot{N}_{{\rm e},\gamma\gamma}(\varepsilon_{\rm e})+\dot{N}_{\rm inj}(\varepsilon_{\rm e}),$$

$${⟨\dot{{\epsilon }_{\mathrm{e}}}⟩}_{\mathrm{ad}}=\frac{1}{3}\frac{\dot{V}}{V}\frac{{\epsilon }_{\mathrm{e}}^2+2{\epsilon }_{\mathrm{e}}{m}_{\mathrm{e}}{c}^2}{{\epsilon }_{\mathrm{e}}+m_{\mathrm{e}}{c}^2}.$$

$$\frac{\partial N_{\gamma }(\epsilon )}{\partial t}={\dot{N}}_{\gamma ,\mathrm{syn}}(\epsilon )+{\dot{N}}_{\gamma ,\mathrm{IC}}(\epsilon )-{\dot{N}}_{\gamma ,\gamma \gamma }(\epsilon )$$

$$-{\dot{N}}_{\gamma ,\mathrm{SSA}}(\epsilon )-{\dot{N}}_{\gamma ,\mathrm{esc}}(\epsilon ),$$

$${\dot{N}}_{\gamma ,\mathrm{esc}}(\epsilon )=\frac{c}{2W}{N}_{\gamma }(\epsilon ),$$

$$\frac{d{E}_{\mathrm{rad}}^{\prime }}{d{t}^{\prime }}=\int d{\epsilon }^{\prime }{\dot{N^{\prime }}}_{\gamma ,\mathrm{esc}}({\epsilon }^{\prime }){\epsilon }^{\prime },$$

$$E_{\mathrm{rad}}=\int d{t}^{\prime }\Gamma \frac{d{E}_{\mathrm{rad}}^{\prime }}{d{t}^{\prime }}.$$

$$\frac{d{\bar{\epsilon }}_{\mathrm{p}}^{\prime }}{d{t}^{\prime }}$$

$$\frac{d{E}_{\mathrm{ad}}^{\prime }}{d{t}^{\prime }}=\frac{M}{m_{\mathrm{p}}}{⟨{\dot{\bar{\epsilon }}}_{\mathrm{p}}⟩}_{\mathrm{ad}}+\int d{\epsilon }_{\mathrm{e}}^{\prime }{N}_{\mathrm{e}}^{\prime }({\epsilon }_{\mathrm{e}}^{\prime }){⟨\dot{{\epsilon }_{\mathrm{e}}}⟩}_{\mathrm{ad}}.$$

$${\epsilon }_{\mathrm{obs}}$$

$$t_{\mathrm{obs}}$$

$$\frac{d{N}_{\gamma }^{\prime }}{d{\Omega }^{\prime }d{S}^{\prime }d{t}^{\prime }d{\epsilon }^{\prime }}=c|\cos {\theta }^{\prime }|\frac{n_{\gamma }^{\prime }({\epsilon }^{\prime },t^{\prime })}{4\pi },$$

$$\Phi ({\epsilon }_{\mathrm{obs}},t_{\mathrm{obs}})=\frac{d{N}_{\gamma }}{d{S}_{\mathrm{obs}}d{t}_{\mathrm{obs}}d{\epsilon }_{\mathrm{obs}}}=$$

$$\int d\theta {\left(\frac{R(t_{\theta }^{\prime })}{D}\right)}^2\frac{\sin \theta |\cos {\theta }^{\prime }|c{n}_{\gamma }^{\prime }({\epsilon }_{\theta }^{\prime },t_{\theta }^{\prime })}{2{\Gamma }^2(1-\beta \cos \theta )(1-{\beta }_{\mathrm{sh}}\cos \theta )},$$

$$\cos {\theta }^{\prime }=\frac{\cos \theta -\beta }{1-\beta \cos \theta },$$

$$R_{\mathrm{dec}}\simeq {\left(\frac{3E_0}{4\pi n_{\mathrm{ISM}}{m}_{\mathrm{p}}{c}^2{\Gamma }_0^2}\right)}^{1/3},$$

$$t_{\mathrm{obs},\mathrm{pk}}\simeq (1+z)\frac{R_{\mathrm{dec}}}{2c{\Gamma }_0^2}\simeq 90(1+z)E_{52}^{\frac{1}{3}}{n}_0^{-\frac{1}{3}}{\Gamma }_2^{-\frac{8}{3}}\text{mbox}s,$$

$${\Gamma }_0\simeq 96E_{52}^{\frac{1}{8}}{n}_0^{-\frac{1}{8}}{t}_{\mathrm{pk},2}^{-\frac{3}{8}},$$

$$\Gamma \simeq \sqrt{\frac{3E_0}{4\pi n_{\mathrm{ISM}}{m}_{\mathrm{p}}{c}^2{R}^3}}.$$

$$R\simeq 1.6\times 10^{17}(1+z{)}^{-\frac{1}{4}}{n}_0^{-\frac{1}{4}}{E}_{52}^{\frac{1}{4}}{t}_{\mathrm{h}}^{\frac{1}{4}}\text{mbox}cm,$$

Full text

1. Introduction

The afterglow emission of gamma-ray bursts (GRBs) is robust evidence of electron acceleration at relativistic shocks. While the difficulty of the particle acceleration by magnetized relativistic shocks has been pointed out by several authors (e.g. Sironi & Spitkovsky, 2009; Lemoine & Pelletier, 2010; Sironi et al., 2013) from the theoretical point of view, the low magnetization implied from recent broadband observations of the afterglows (e.g. Kumar & Barniol Duran, 2009; Lemoine et al., 2013; Santana et al., 2014; Beniamini et al., 2015; Zhang et al., 2015) seems to be consistent with the theoretical argument.

The physical property of the shock and electron acceleration in the GRB afterglows has been investigated with the conventional microscopic parameters, the energy fractions of the accelerated electrons $\epsilon_{\rm e}$ and magnetic field $\epsilon_{B}$ in the downstream. The observations and standard analytical formulae of the external shock model by Sari, Piran & Narayan (1998) provide those microscopic parameters and jet parameters (e.g. Panaitescu & Kumar, 2002; Lloyd-Ronning & Zhang, 2004). In those analytical formulae, the electron energy distribution at a given radius is assumed to be a broken-power-law. The results of the two-dimensional hydrodynamical simulations by van Eerten et al. (2012) with the analytical broken-power-law formula have been widely used to fit the observed light curves (see, e.g. Guidorzi et al., 2014; Maselli et al., 2014; Zhang et al., 2015). Such multidimensional hydrodynamical simulations provide precise evolution of the shock propagation and angular structure of the collimated jet and are a powerful tool to constrain jet parameters, especially in the off-axis cases.

On the other hand, the actual electron energy distribution in the downstream of the propagating shock is not a simple broken-power-law. Petropoulou & Mastichiadis (2009) and Pennanen et al. (2014) calculated the evolution of the electron energy distribution in the afterglow. The resultant photon spectra are significantly curved around the cooling and injection break frequencies, and not the broken-power-law (see also Uhm & Zhang, 2014).

Some fraction of GRB afterglows are hard to explain with the standard external shock model (e.g. Willingale et al., 2007; Wang et al., 2015). Multizone models such as the spine-sheath structure (Racusin et al., 2008), the contribution of the reverse shock (Genet et al., 2007; Uhm & Beloborodov, 2007), or evolving microscopic parameters (Ioka et al., 2006) may be required to reconcile such exceptional afterglows. Before increasing the number of parameters following such complex models, however, we need to clarify the degree of the contradiction with the standard external shock model. In addition to the uncertainty of the electron energy distribution, the detections of the GeV afterglows with Fermi (Abdo et al., 2009; Kumar & Barniol Duran, 2010) require us to investigate seriously the effect of synchrotron self-Compton (SSC) emission on the spectrum and light curve. Especially, the detection of a 32 GeV photon at $3\times 10^{4}$ s in GRB 130427A (Ackermann et al., 2014) cannot be explained with the usual synchrotron emission for the standard evolution of the external shock. The SSC emission spectra numerically obtained (Petropoulou & Mastichiadis, 2009; Pennanen et al., 2014) are naturally different from a broken-power-law derived from the analytical description (Sari & Esin, 2001). In addition, if the shocked plasma is in the highly radiative regime as discussed in Ghisellini et al. (2010), the radiative cooling affects not only the electron energy distribution but also the evolution of the bulk Lorentz factor. When we treat all of the above nontrivial effects numerically without analytical approximations, the spectrum and light curve may deviate from the behaviors given by simple formulae.

In this paper, in order to discuss the uncertainty of the evolution of the emission from the external shock, we simulate the evolutions of the shocked material propagating in the interstellar medium (ISM). Our numerical code is based on the one-zone approximation, but the time-dependent treatment is completely applied for the bulk motion of the shell and the electron and photon energy distributions. Our method is similar to that in the previous studies (Petropoulou & Mastichiadis, 2009; Pennanen et al., 2014; Uhm & Zhang, 2014), but the light curves were not calculated in their studies. Our code consistently transforms the energy and arrival time of photons that escaped from the shocked shell into those for an observer. The spectrum for the observer at a certain time $t_{\rm obs}$ is not just the blue-shifted one in the shell comoving frame at the time $t^{\prime}$ given by the one-to-one correspondence between $t_{\rm obs}$ and $t^{\prime}$. Focusing on the light curve and spectral evolution for the observer, we discuss the differences in the results obtained with the analytical method and ours. We also discuss the switching signature from synchrotron to inverse Compton in the gamma-ray light curve. In §2, we present our computing method. The analytical formulae in §3 are compared with the numerical results in §4. Our model is applied to the afterglow of GRB 130427A in §5. The smooth gamma-ray light curve is reproduced in spite of the switching of the emission process in the GeV energy range. The conclusions are summarized in §6.

2. Model and Method

In this paper, we assume a spherically symmetric system, which may be an appropriate assumption before the jet break. We treat the shocked region propagating in the ISM as a uniform shell with a thickness $W=W^{\prime}/\Gamma$, where $\Gamma\equiv 1/\sqrt{1-\beta^{2}}$ is the bulk Lorentz factor of the shocked region. Hereafter, we denote values in the shell frame by primed characters. Under the one-zone approximation, we numerically solve the evolutions of $\Gamma$, magnetic field, and energy distributions of the photons and non-thermal electrons in the shell in a self-consistent manner. The model parameters are the total energy $E_{0}$ promptly released from the central engine, the initial bulk Lorentz factor of the ejecta $\Gamma_{0}$, the proton number density of the ISM $n_{\rm ISM}$, the spectral index $p$ and the number fraction $\eta$ of non-thermal electrons, and the energy fractions to the shock-dissipated energy, $\epsilon_{\rm e}$ and $\epsilon_{B}$, of non-thermal electrons and magnetic field, respectively. Below, we present the evolution of the shell regardless of whether the shell motion is relativistic or not. The average kinetic energy per proton just behind the shock front is $(\Gamma-1)m_{\rm p}c^{2}$; hence, we can obtain the temperature $T_{\rm sh}$ from

[TABLE]

where $x\equiv T_{\rm sh}/(m_{\rm p}c^{2})$ and $K_{n}(x)$ is the modified Bessel function of the second kind. Given $x$ and $\Gamma$, the heat capacity ratio is written as

[TABLE]

The shock jump condition (Blandford & Mckee, 1976) provides the bulk Lorentz factor of the shock front $\Gamma_{\rm sh}$ as

[TABLE]

When the shock front is propagating at a radius $R$ from the central engine as $dR/dt=c\beta_{\rm sh}$ ($dt=\Gamma dt^{\prime}$, $\beta_{\rm sh}\equiv\sqrt{1-1/\Gamma_{\rm sh}^{2}}$), the mass in the shell evolves as

[TABLE]

with the initial mass $M_{0}=E_{0}/(\Gamma_{0}-1)/c^{2}$. The total energy including the rest mass energy in the comoving frame evolves as

[TABLE]

where the first through third terms on the right-hand side express the energy injection, radiative cooling, and adiabatic cooling, respectively. For each time step, we numerically follow the evolutions of the shell mass and energy with Equations (4) and (5) and obtain $\Gamma$ from the energy conservation

[TABLE]

Since we assume a homogeneous shell, the density obtained by the jump condition,

[TABLE]

is adopted for the entire shell. According to the evolutions of $M$ and $\Gamma$, the shell volume is written as $V^{\prime}=M/(m_{\rm p}n^{\prime})$. Assuming that a fraction $\epsilon_{B}$ of the injected kinetic energy converts into the magnetic field, the magnetic energy $E_{B}$ is calculated by

[TABLE]

The magnetic field is estimated by

[TABLE]

The evolution of the electron and photon energy distributions in the shell frame is calculated with the same method as in Asano & Mészáros (2011). Non-thermal electrons (number fraction $\eta$) are assumed to obtain a fraction $\epsilon_{\rm e}$ of the injected kinetic energy. Assuming a cut-off power-law spectrum at injection

[TABLE]

for $\varepsilon^{\prime}\geq\varepsilon^{\prime}_{\rm min}$, the number and energy injection rates are written as

[TABLE]

The maximum electron energy $\varepsilon^{\prime}_{\rm max}$ is obtained by equating the acceleration time $\xi r^{\prime}_{\rm L}/c$ (or $20\xi r^{\prime}_{\rm L}/(3\beta^{2}c)$ for the non-relativistic case) and cooling time $\varepsilon_{\rm e}/\dot{\varepsilon_{\rm e}}$ due to synchrotron and inverse Compton emissions numerically obtained, where $r_{\rm L}$ is the Larmor radius. Hereafter, the Bohm factor $\xi$ is optimistically assumed as unity. Then, Equations (11) and (12) provide the normalization $\dot{N^{\prime}_{0}}$ and $\varepsilon^{\prime}_{\rm min}$ for given $dM/dt^{\prime}$, $\Gamma$, $\eta$, and $\epsilon_{\rm e}$.

In this paragraph, to explain the method for following the evolution of the electron/positron/photon energy distributions, we omit the prime symbol and express equations in the shell frame. Our numerical code practically solves the evolution equation of non-thermal electrons/positrons

[TABLE]

where $\left<\dot{\varepsilon_{\rm e}}\right>_{\rm syn}$ and $\left<\dot{\varepsilon_{\rm e}}\right>_{\rm IC}$ are the energy loss rates (positive values) due to synchrotron and inverse Compton (IC) emissions, respectively. The Klein–Nishina effect is numerically taken into account using the table of the emissivity prepared in advance with the Monte Carlo method (see Asano & Mészáros, 2011). The electron heating rate due to the synchrotron self-absorption (SSA) is also included as denoted with $\left<\dot{\varepsilon_{\rm e}}\right>_{\rm SSA}$. The extra term of electron–positron pair injection due to $\gamma\gamma$-absorption is $\dot{N}_{{\rm e},\gamma\gamma}(\varepsilon_{\rm e})$. The adiabatic cooling term $\left<\dot{\varepsilon_{\rm e}}\right>_{\rm ad}$ is calculated from the momentum evolution $\dot{p_{\rm e}}=-p_{\rm e}\dot{V}/(3V)$. Since the kinetic energy is $\varepsilon_{\rm e}=\sqrt{p_{\rm e}^{2}c^{2}+m^{2}_{\rm e}c^{4}}-m_{\rm e}c^{2}$, the cooling rate is written as

[TABLE]

The pair production, IC emission, and SSA depend on the photon density $n_{\gamma}(\varepsilon)=N_{\gamma}(\varepsilon)/V$. The photon energy distribution is obtained by solving

[TABLE]

where the first and second terms on the right-hand side represent synchrotron and IC photon production, respectively, and the third and fourth terms represent photon absorption due to $\gamma\gamma$ and SSA, respectively. Those terms are numerically calculated with the given electron and photon distributions and magnetic field. Photons escape from both the front and rear surfaces, so that the escape term is written as

[TABLE]

where the shell width $W=V/(4\pi R^{2})$.

Using the prime symbol again hereafter, the radiative cooling term in Equation (5) is obtained as

[TABLE]

and the radiation term in Equation (6) is calculated with

[TABLE]

Although we do not solve the proton energy distribution explicitly, the adiabatic cooling of protons is essential for the evolution of $\Gamma$. The energy injection rate into protons is $dE^{\prime}_{\rm p}/dt^{\prime}=(1-\epsilon_{\rm e}-\epsilon_{B})(\Gamma-1)c^{2}dM/dt^{\prime}$. The average kinetic energy of protons $\bar{\varepsilon}^{\prime}_{\rm p}\equiv m_{\rm p}E^{\prime}_{\rm p}/M$ evolves as

[TABLE]

where the last term is the same form as Equation (14) with ${\rm e}\to{\rm p}$. This simplified method provides the adiabatic energy loss rate

[TABLE]

With Equations (4), (17), and (20), the total shell energy is calculated from Equation (5). Then, we can obtain the Lorentz factor $\Gamma$ from Equation (6) for each time step.

In order to obtain the photon spectrum and light curve for an observer, we integrate photons over the shell surface. The method for the time and energy transformations is also the same as in Asano & Mészáros (2011). The energy and arrival time of photons escaping from the surface expanding toward an angle $\theta$ to the line of sight at radius $R$ are written as

[TABLE]

where $t=\int\Gamma dt^{\prime}$, $R=c\int\beta_{\rm sh}\Gamma dt^{\prime}$, and $R_{0}$ is the initial radius. In the comoving frame, the photon escape rate per unit surface per solid angle is written as

[TABLE]

where $dS^{\prime}=dS=2\pi R^{2}\sin{\theta}d\theta$ is the surface element. While the number of photons $dN^{\prime}_{\gamma}$ is obviously Lorentz invariant, the infinitesimal intervals are transformed as $dt_{\rm obs}=(1+z)(1-\beta_{\rm sh}\cos{\theta})\Gamma dt^{\prime}$, $d\varepsilon_{\rm obs}=d\varepsilon^{\prime}/\left\{(1+z)(1-\beta\cos{\theta})\Gamma\right\}$, and $d\Omega=\Gamma^{2}(1-\beta\cos{\theta})^{2}d\Omega^{\prime}$ for solid angle. Denoting the luminosity distance as $D_{\rm L}=(1+z)D$, the surface through which photons traveling toward $d\Omega$ pass is written as $dS_{\rm obs}=D^{2}d\Omega$. Then, we obtain the photon flux for an observer as

[TABLE]

where

[TABLE]

and $\varepsilon^{\prime}_{\theta}$ and $t^{\prime}_{\theta}$ are the comoving energy and time obtained from Equations (21) and (LABEL:ttrf), respectively, for given $\theta$, $\varepsilon_{\rm obs}$, and $t_{\rm obs}$. Notice that $\beta$ (equivalently $\Gamma$) and $\beta_{\rm sh}$ are also functions of $t^{\prime}_{\theta}$. Carrying out the integral in Equation (24) numerically over $\theta$, we can obtain the spectral evolution for an observer.

3. Analytical Behavior: Review

While we numerically follow the evolution of the photon spectrum for an observer with the method explained in the previous section, here we review the analytical description in Sari, Piran & Narayan (1998) to compare with our results. When the shock is ultra-relativistic ($\Gamma\gg 1$), $\hat{\gamma}\simeq 4/3$ and $\Gamma_{\rm sh}\simeq\sqrt{2}\Gamma$. Until the deceleration radius (Rees & Mészáros, 1992),

[TABLE]

the shell expands with a constant Lorentz factor $\Gamma_{0}$. The peak time of the afterglow for an observer corresponds to this radius as

[TABLE]

where $E_{0}=10^{52}E_{52}~{}\mbox{erg}$, $n_{\rm ISM}=n_{0}~{}\mbox{cm}^{-3}$, and $\Gamma_{0}=100\Gamma_{2}$. When the peak time is determined observationally, the initial Lorentz factor is estimated as

[TABLE]

where $t_{\rm obs,pk}/(1+z)=100t_{\rm pk,2}$ s. After the peak time, the shell starts to decelerate. Since the shell density is $n^{\prime}\simeq 4\Gamma n_{\rm ISM}$, the shell width becomes $W\simeq R/(12\Gamma^{2})$ in the one-zone approximation. The jump condition provides the energy density $U^{\prime}=4\Gamma^{2}n_{\rm ISM}m_{\rm p}c^{2}$. Neglecting the radiative cooling, the energy conservation implies that the Lorentz factor decreases as

[TABLE]

The one-zone approximation in the above equation has a slightly different factor from that in Sari, Piran & Narayan (1998), where the radial density structure behind the shock is taken into account to estimate $E_{0}$. The simple one-to-one correspondence for $R$ and $t_{\rm obs}$, $t_{\rm obs}\simeq(1+z)R/(4c\Gamma^{2})$, implies

[TABLE]

where $t_{\rm obs}=t_{\rm h}$ hr.

From Equations (11) and (12), we obtain the electron minimum Lorentz factor $\gamma^{\prime}_{\rm m}-1\equiv\varepsilon^{\prime}_{\rm min}/(m_{\rm e}c^{2})$ as

[TABLE]

A fraction $\epsilon_{B}$ of the energy density converts to the magnetic field as

[TABLE]

The typical synchrotron photon energy is obtained as

[TABLE]

where $(p-2)/(p-1)=1/6f_{1/6}$, $\epsilon_{\rm e}=0.1\epsilon_{{\rm e},-1}$, and $\epsilon_{B}=0.1\epsilon_{B,-1}$. Given the photon energy $\varepsilon_{\rm obs}=\varepsilon_{\rm eV}$ eV in observation, $\varepsilon_{\rm m}$ passes at an observer time

[TABLE]

In the electron energy distribution, the cooling break appears at a Lorentz factor

[TABLE]

Note that the formulation in this section neglects the effect of IC cooling. When IC is dominant for electron cooling, $\gamma^{\prime}_{\rm c}$ and its evolution will be modified. In our numerical simulations shown in the next section, those non-linear effects due to IC are automatically included. From Equation (36), the cooling break energy becomes

[TABLE]

The radius corresponding to $\varepsilon_{\rm c}=\varepsilon_{\rm obs}$ is written as

[TABLE]

which corresponds to the observer time

[TABLE]

According to the high/low relation between $\varepsilon_{\rm m}$ and $\varepsilon_{\rm c}$, the non-thermal electrons are judged as the fast cooling ($\varepsilon_{\rm m}>\varepsilon_{\rm c}$) or the slow cooling ($\varepsilon_{\rm m}<\varepsilon_{\rm c}$). In the early stage, the strong synchrotron cooling may lead to the fast cooling. The transition from the fast cooling to the slow cooling occurs at

[TABLE]

which corresponds to

[TABLE]

for the observer. To obtain the SSA frequency, assuming that all the injected electrons form a single power-law above $\min(\gamma^{\prime}_{\rm m},\gamma^{\prime}_{\rm c})$, we calculate the usual absorption formula with the synchrotron function. When $t_{\rm obs}<t_{\rm eq}$, the SSA frequency is

[TABLE]

where $t_{\rm obs}=100t_{2}$ s. For $t_{\rm obs}>t_{\rm eq}$, we obtain a constant value

[TABLE]

where

[TABLE]

which is $\sim 4.35$ for $p=2.2$.

As shown in Sari, Piran & Narayan (1998), the maximum flux, $F_{\rm max}=\varepsilon_{\rm c}\Phi(\varepsilon_{\rm c})$ for $t_{\rm obs}<t_{\rm eq}$ or $\varepsilon_{\rm m}\Phi(\varepsilon_{\rm m})$ for $t_{\rm obs}>t_{\rm eq}$, becomes constant for $t_{\rm obs}>t_{\rm obs,pk}$ as

[TABLE]

where $D_{\rm L}=10^{28}D_{28}$ cm. Normalizing the flux $F(\varepsilon_{\rm obs})=\varepsilon_{\rm obs}\Phi(\varepsilon_{\rm obs})$ by $F_{\rm max}$ at $\min(\varepsilon_{\rm m},\varepsilon_{\rm c})$, the broken-power-law formula yields the spectral evolution as follows:

[TABLE]

for $t_{\rm obs}<t_{\rm obs,eq}$, and

[TABLE]

for $t_{\rm obs}>t_{\rm obs,eq}$. In this section, we do not consider the cases of $\min(\varepsilon_{\rm m},\varepsilon_{\rm c})<h\nu_{\rm a}$, in which case the spectral shape and its evolution should be modified (e.g. Granot & Sari, 2002). For the parameter regions adopted in our simulations (see the next section), the self-absorption frequency is safely suppressed below $\min(\varepsilon_{\rm m},\varepsilon_{\rm c})$.

Before the peak time $t_{\rm pk}$, $\Gamma$ and $B^{\prime}$ are constant, so that the maximum flux increases as $F_{\rm max}\propto N_{\rm e}\propto R^{3}\propto t_{\rm obs}^{3}$. The characteristic photon energies behave as $\varepsilon_{\rm m}\propto t_{\rm obs}^{0}$ and $\varepsilon_{\rm c}\propto t_{\rm obs}^{-2}$. In the fast cooling case, $F(\varepsilon_{\rm obs})=F_{\rm max}(\varepsilon_{\rm m}/\varepsilon_{\rm c})^{-1/2}(\varepsilon_{\rm obs}/\varepsilon_{\rm m})^{-p/2}$ for $\varepsilon_{\rm obs}>\varepsilon_{\rm m}$, and $F(\varepsilon_{\rm obs})=F_{\rm max}(\varepsilon_{\rm obs}/\varepsilon_{\rm c})^{-1/2}$ for $\varepsilon_{\rm c}<\varepsilon_{\rm obs}<\varepsilon_{\rm m}$. Then, as long as $\varepsilon_{\rm obs}>\varepsilon_{\rm c}$, we obtain

[TABLE]

for $t_{\rm obs}<t_{\rm obs,pk}$. Similarly, for $\varepsilon_{\rm obs}<\varepsilon_{\rm c}$,

[TABLE]

4. Numerical Results: Spectrum and Light Curve

In our model, there are seven parameters. As a benchmark case, hereafter we adopt $E_{0}=10^{52}$ erg, $\Gamma_{0}=100$, $n_{\rm ISM}=1~{}\mbox{cm}^{-3}$, $p=2.2$, $\epsilon_{\rm e}=0.1$, $\epsilon_{B}=0.1$, and $\eta=1$. Figure 1 shows the evolution of $\Gamma\beta$ and $B^{\prime}$ numerically obtained for the benchmark parameter set. The deceleration starts at a slightly smaller radius than the deceleration radius expressed by Equation (26). The evolution of $\Gamma\beta$ agrees with the adiabatic approximation from Blandford–Mckee ($\Gamma\propto t^{-3/2}$ in the relativistic regime) to Sedov–Taylor ($\beta\propto t^{-3/5}$ in the non-relativistic regime) phases. The decay of the magnetic field also follows the evolution of $\Gamma\beta$ as expressed in Equation (32), though a slight deviation from the approximation $B^{\prime}\propto\Gamma$ is seen below $\Gamma\beta<10$, where the approximation $n^{\prime}=4\Gamma n_{\rm ISM}$ is not so accurate.

We also test the radiative case with $\epsilon_{\rm e}=0.9$ and $\epsilon_{B}=0.05$, where the shock-dissipated energy is efficiently released by radiation. The other parameters are the same as those in the benchmark case. As shown by the dashed line in Figure 1, the numerical result shows $\Gamma\propto t^{-2}$ in the relativistic regime, while the well-known formula for the radiative shock (Blandford & Mckee, 1976; Sari, Piran & Narayan, 1998) is $\Gamma\propto t^{-3}$. The analytic formula in the radiative case is based on the approximation $E_{\rm iso}\simeq\Gamma M_{0}c^{2}$ neglecting the increase of the inertia for $R<[3M_{0}/(4\pi m_{\rm p}n_{\rm ISM})]^{1/3}$ ($\simeq 2.5\times 10^{17}$ cm in this case). Even for $\epsilon_{\rm e}=0.9$, however, the shocked ISM of mass $\Delta M$ adds inertia $\geq(1-\epsilon_{\rm e})\Gamma\Delta M=0.1\Gamma\Delta M$, which is larger than $\Delta M$ for $\Gamma>10$. Actually, for $t>3\times 10^{6}$ s, the increase of the inertia cannot be negligible. In addition, the faster decay of the magnetic fields leads to the suppression of the radiative efficiency. The increase of the inertia and decrease of the radiative efficiency lead to $\Gamma\propto t^{-2}$ rather than $\Gamma\propto t^{-3}$ in this parameter set.

Figure 2 shows the evolution of $n^{\prime}_{\rm e}(\varepsilon^{\prime}_{\rm e})\equiv N^{\prime}_{\rm e}(\varepsilon^{\prime}_{\rm e})/V^{\prime}$ for the benchmark case. In our numerical code, electrons are injected intermittently. In the highest-energy region, the interval of the electron injection is longer than the cooling time scale, which results in the fluctuation of the electron distribution as seen in the figure. However, this does not practically affect the photon spectrum because of the longer photon escape time than the electron injection interval. The low-energy component below $10^{8}$ eV seen in the early stage is down scattered particles with photons.

Initially the system is in the fast cooling regime. The steady analytical solution for the fast cooling is $n^{\prime}_{\rm e}(\varepsilon^{\prime}_{\rm e})\propto\varepsilon^{\prime-2}_{\rm e}$. In our time-dependent treatment, the injection rate increases with time, so that the electron distribution below $\varepsilon^{\prime}_{\rm min}$ (e.g. $5.6\times 10^{8}$ eV at $t=10^{6.5}$ s) is harder than the steady solution. Equation (40) indicates that the electron distribution should be expressed with the slow cooling approximation for $t>10^{6.5}$ s. Actually, a sharp low-energy cutoff appears below $\varepsilon^{\prime}_{\rm min}$ for $t\geq 10^{7}$ s. The cooling break in the electron spectrum in the slow cooling regime (see, e.g. $\sim 10^{9}$ eV at $t=10^{7.5}$ s) is not so sharp; the simple broken-power-law function is not appropriate for our results.

The obtained photon spectrum for an observer in the benchmark case is rather simple as shown in Figure 3. According to Equation (41), the photon spectrum must be the shape described with the slow cooling approximation for $t_{\rm obs}>10^{3}$ s. However, the spectra are smoothly curved around the peak, so that it is hard to identify the spectral break at $\varepsilon_{\rm m}$ and $\varepsilon_{\rm c}$ (see Figure 4). In this parameter set, the synchrotron and SSC components almost merge; the spectrum seems to consist of a single component. As shown by the comparison of the solid and dot-dashed curves in Figure 3, the SSC component dominates above 0.1 GeV for $t_{\rm obs}=10^{4}$ s.

At $t_{\rm obs}=10^{4}$ s, we compare the photon spectrum of the benchmark case with the results for other parameter sets in Figure 4. When a parameter $\epsilon_{B}$ is reduced from 0.1 to $10^{-3}$ (see blue curve in Figure 4), $\varepsilon_{\rm c}$ grows into the X-ray range, and the SSC component is clearly seen at the TeV energy range. The analytical estimate implies $\gamma_{\rm c}/\gamma_{\rm m}\sim 300$ at this time. Following Sari & Esin (2001), the ratio of the IC peak flux to the synchrotron one is estimated as $(\gamma_{\rm c}/\gamma_{\rm m})^{(2-p)/2}\sqrt{\epsilon_{\rm e}/\epsilon_{B}}\sim 6$, while the numerical result shows a slightly dimmer IC flux than the synchrotron flux. This discrepancy may be partially due to the Klein–Nishina effect, but the time-dependent treatment apparently affects the flux ratio. For the result of $p=3$ (see red curve in Figure 4), the soft synchrotron component makes the SSC component easier to distinguish even for $\epsilon_{B}=0.1$ (see also the thin green line, which corresponds to the monoenergetic injection). Adopting Equations (34), (37), and (46), we also plot the analytic spectra of the synchrotron component in Figure 4, where the maximum photon energy is simply assumed as $0.1\Gamma/(1+z)$ GeV. The numerical results show curved spectra rather than the broken-power-law. Those curved features are similar to the time-dependent calculations in Petropoulou & Mastichiadis (2009, see also ()). The analytic broken-power-law formula roughly reproduces the overall spectral shape. As shown in Figure 4, however, the analytic fluxes are slightly overestimated for $\varepsilon_{\rm obs}\gtrsim\varepsilon_{\rm c}$ compared to our results.

The light curves for various photon energies are plotted in Figure 5. The X-ray peak time is 2.4 times earlier than Equation (27). If we adopt Equation (28) with the numerically obtained peak time to estimate $\Gamma_{0}$, this discrepancy leads to about 40% larger $\Gamma_{0}$ compared to the actual value. While Equation (28) is the same as the standard formula (Sari & Piran, 1999; Zhang et al., 2003; Molinari et al., 2007), the formula in Liang et al. (2010) is two times larger. If our $t_{\rm obs,pk}$–$\Gamma_{0}$ relation obtained numerically is adopted, the resultant $\Gamma_{0}$ becomes 2.8 times smaller than the results in Liang et al. (2010), in which a relation $\Gamma_{0}\simeq 182(E_{\gamma,\mbox{iso}}/10^{52}\mbox{erg})^{0.25}$ is obtained from the afterglow onset times of 17 GRB samples.

Before the peak time, the X-ray flux grows as $F\propto t_{\rm obs}^{1.4}$, while the simple analytical estimate leads to $F\propto t_{\rm obs}^{2}$ (see Equation (49)). The fitting of the early X-ray light curves for 11 GRB samples by Liang et al. (2010) shows a large scatter in the rising indices from 0.5 to $\sim 10$. At $\varepsilon_{\rm obs}=1$ eV, which is below $\varepsilon_{\rm c}$ in the early stage, the rising index in our calculation is about 2.2, which is also smaller than the analytical estimate $11/3$ (Equation (50)).

The standard analytic model (Sari, Piran & Narayan, 1998; Sari & Piran, 1999) predicts the evolution for the 1eV light curve as $F\propto t^{1/6}$ for $t_{\rm obs,pk}<t_{\rm obs}<t_{\rm obs,eq}$, $F\propto t^{1/2}$ for $t_{\rm obs,eq}<t_{\rm obs}<t_{\rm obs,m}$, and $F\propto t^{3(1-p)/4}=t^{-0.9}$ for $t_{\rm obs,m}<t_{\rm obs}<t_{\rm obs,c}$ (see Equations (47) and (48)). Since the photon spectrum is curved around $\varepsilon_{\rm m}$ or $\varepsilon_{\rm c}$, our result does not show such sharp breaks in the 1eV light curve. For $t>t_{\rm obs,c}$, the light curves below MeV seem consistent with the decay $F\propto t^{-1.15}$ depicted in Equation (48).

We have optimistically assumed the Bohm limit $\xi=1$ to maximize the maximum electron energy $\varepsilon^{\prime}_{\rm max}$. As a result, below $0.1\Gamma(1+z)^{-1}$ GeV, the synchrotron emission dominates. Thus, even if we neglect the SSC emission, the light curves below 0.1 GeV are almost unchanged, while the 0.1 TeV emission is greatly suppressed.

As we have mentioned, the analytic formula reviewed in section 3 tends to overestimate the X-ray flux (see Figure 4). Here, we define the “discrepancy factor” of the analytical flux based on Equations (34), (37), and (46) as the ratio of

[TABLE]

at 1 keV at $t_{\rm obs}=10^{4}$ s assuming $z=2$. We change the three model parameters as $E_{52}=0.1$, 1, and 10, $\epsilon_{B}=0.1$, $10^{-2}$, and $10^{-3}$, and $p=2.2$, $2.5$, and $3.0$. In total we test $3\times 3\times 3=27$ models, keeping the other parameters as $\Gamma=100$, $n_{\rm ISM}=1~{}\mbox{cm}^{-3}$, $\epsilon_{\rm e}=0.1$, and $\eta=1$. The discrepancy factors are shown in Figure 6. Although the normalization of the electron number $\dot{N^{\prime}_{0}}$ is proportional to $(p-1)$, which is neglected in the estimate of $F_{\rm max}$ in Equation (46), our results do not show a clear dependence on $p$ in the discrepancy factors. In the cases where the SSC emission is efficient, namely, smaller $\epsilon_{B}$ and larger $E_{0}$, the discrepancy factor tends to be large. This tendency agrees with the analytical discussion in Beniamini et al. (2016).

Note that the normalization of the flux in Equation (46) is basically the same as used in Lloyd-Ronning & Zhang (2004) to estimate the kinetic energy $E_{0}$ in the afterglow phase. Lloyd-Ronning & Zhang (2004) raised the problem that the prompt gamma-ray emission is too efficient compared to the remnant kinetic energy at the afterglow onset. The analytical formulation reviewed in section 3 shows that the synchrotron flux in the highest-energy region is proportional to $F_{\rm max}\varepsilon_{\rm c}^{1/2}\varepsilon_{\rm m}^{(p-1)/2}\propto E_{0}^{(p+2)/4}$ ($E^{1.05}$ for $p=2.2$). Therefore, $E_{0}$ obtained from the analytic formula results in an underestimate of $E_{0}$ by a factor close to the discrepancy factors shown in Figure 6.

Here we have assumed $\epsilon_{B}\geq 10^{-3}$ so that the X-ray-emitting electrons are in the fast cooling regime at $t_{\rm obs}=10^{4}$ s. Namely, the X-ray band is above both $\varepsilon_{\rm m}$ and $\varepsilon_{\rm c}$. However, the recent broadband observations (e.g. Kumar & Barniol Duran, 2009; Lemoine et al., 2013; Santana et al., 2014; Zhang et al., 2015) suggest much lower magnetization. Beniamini et al. (2015, 2016) pointed out that the high-efficiency problem in the prompt emission will be resolved by a very small $\epsilon_{B}$. In this case, the X-ray-emitting electrons are in the slow cooling regime, and the synchrotron flux in the X-ray band can be suppressed by the IC cooling. Those two effects may lead to a wrong estimate of $E_{0}$, if we adopt the usual fast cooling formula for the X-ray-emitting electrons.

We also test the discrepancy factor for $\epsilon_{B}=10^{-4}$–$10^{-6}$ as shown in Figure 7. In most of those cases, the X-ray band is between $\varepsilon_{\rm m}$ and $\varepsilon_{\rm c}$. The analytical broken-power-law formula tends to underestimates the flux for $\varepsilon_{\rm m}<\varepsilon_{\rm obs}<\varepsilon_{\rm c}$ contrary to the fast cooling regime ($\varepsilon_{\rm c}<\varepsilon_{\rm obs}$). Figure 7 shows that the analytical formula underestimate the flux by a factor of 1.2–3.

This implies that the total energy obtained from the analytical formula tends to be overestimated for smaller $\epsilon_{B}$, which will worsen the high-efficiency problem in the prompt emission. However, if we misunderstand the X-ray energy range as the fast cooling regime, the estimate of $E_{0}$ can be less than $1$% of the actual energy depending on the parameters. This misinterpretation can be a major factor that leads to an underestimate of $E_{0}$ as pointed out by Beniamini et al. (2015, 2016). In such cases, the discrepancy shown in Figure 7 seems negligible.

As is understood from Figures 6 and 7, the suppression of the X-ray flux by the IC cooling becomes maximum at $\epsilon_{B}\sim 10^{-3}$. For an extremely small $\epsilon_{B}$, though the IC cooling becomes relatively dominant compared to the synchrotron cooling, the cooling effect itself is negligible. Therefore, the synchrotron emissivity is not largely affected by the radiative cooling.

We also test the X-ray closure relation between the decay index $\alpha_{\rm X}$ ($F\propto t_{\rm obs}^{-\alpha_{\rm X}}$) and spectral index $\beta_{\rm X}$ ($F\propto\varepsilon_{\rm obs}^{-\beta_{\rm X}}$). The analytical formula of Equation (48) indicates $\alpha_{\rm X}=(3p-2)/4$ and $\beta_{\rm X}=p/2$ for $\varepsilon_{\rm obs}>\varepsilon_{\rm c}$. This implies the closure relation $(3/2)\beta_{\rm X}-\alpha_{\rm X}=1/2$. We may expect deviation from the closure relation in the numerical results. In addition to the 27 models in Figure 6, we change the initial Lorentz factor as $\Gamma=100$ and 300 and obtain $\alpha_{\rm X}$ and $\beta_{\rm X}$ at 1 keV and $t_{\rm obs}=10^{4}$ s assuming $z=2$. The numerical results for all 54 models are plotted in Figure 8. In most cases, the results slightly deviate from the analytical relation for the corresponding $p$. As $\epsilon_{B}$ decreases, the spectrum tends to be hard, and the flux decay tends to be shallow. Nevertheless, our numerical results in Figure 8 distribute along the analytic closure relation. Those points are slightly above the closure relation systematically.

The result for the radiative model ($\epsilon_{\rm e}=0.9$ and $\epsilon_{B}=0.05$) is also plotted in Figure 8. In the analytic model (Sari, Piran & Narayan, 1998), $\alpha_{\rm X}=(6p-2)/7=1.6$ and $\beta_{\rm X}=1.1$ in this case. The obtained values are $(\alpha_{\rm X},\beta_{\rm X})=(1.38,1.03)$. The shallower decay of $\Gamma$ than the analytic model as shown in Figure 1 causes the shallower flux decay. Even in this extreme model, $\alpha_{\rm X}\sim 2$ is hard to realize.

As we have mentioned, for $\epsilon_{B}\leq 10^{-4}$, the X-ray band is in the slow cooling regime in most cases. In this case, the closure relation becomes $(3/2)\beta_{\rm X}-\alpha_{\rm X}=0$. As shown in Figure 9, the numerical results show slight deviation from the closure relation (the red line), while the distribution of the decay indices $\alpha_{\rm X}$ is consistent with the analytical formula. The $\alpha_{\rm X}$–$\beta_{\rm X}$ distribution in this case does not contradict the observed distribution.

However, the large scatter in the samples in Willingale et al. (2007) is not explained by the model with constant $n_{\rm ISM}$, $\epsilon_{\rm e}$, $\epsilon_{B}$, or $\eta$. Our model does not take into account the shallow decay phase, which may require the energy injection with a longer timescale than $t_{\rm obs,pk}$. Though some additional parameters with respect to the energy injection may resolve the problem, such a complex model is beyond the scope in this paper.

5. Application to GRB 130427A

GRB 130427A (Ackermann et al., 2014; Maselli et al., 2014) is a very nearby GRB ($z=0.34$) with a very large isotropic energy release ($8.5\times 10^{53}$ erg; Perley et al., 2014) in gamma-ray. Surprisingly, the X-ray afterglow flux is well fitted by a simple power-law of $t^{-1.309}$ (hereafter we omit the subscript “obs”) as far as $8\times 10^{7}$ s without a signature of the jet break (De Pasquale et al., 2016). Thanks to the very large fluence, a long-lasting GeV emission as far as $7\times 10^{4}$ s was detected with Fermi. The most enigmatic problem in this GRB is the detection of a 32 GeV photon at $t\simeq 3\times 10^{4}$ s. The maximum synchrotron photon energy is limited by the balance between the energy loss and gain as $\sim 0.1$ GeV irrespective of the magnetic field. While the relativistic motion can boost the photon energy in the early stage, the Lorentz factor should be significantly reduced at the arrival time of the 32 GeV photon. Therefore, the 32 GeV photon may be emitted via SSC process (Fan et al., 2013; Liu et al., 2013). However, the GeV light curve is well fitted by a single power-law in the late phase and does not show the signature of the transition from synchrotron to SSC.

The decay indices of the optical and X-ray light curves are different, which implies that there are difficulties in the standard ISM model with a single emission component. Especially for the early stage of the afterglow, complicated models including a reverse shock component with a stellar-wind profile (Laskar et al., 2013; Perley et al., 2014; Vestrand et al., 2014), spine-sheath-like two-component jets (van der Horst et al., 2014), temporally evolving microscopic parameters ($\eta$, $\epsilon_{\rm e}$ and $\epsilon_{B}$; Maselli et al., 2014), or a non-standard radial density profile (Kouveliotou et al., 2013; van der Horst et al., 2014) have been attempted. Nevertheless, we adopt our single-component model focusing on the GeV–TeV emission in the early phase; our numerical method is ideal to calculate the GeV–TeV light curve with SSC emission consistently as explained in the previous section. We have tested three models, whose parameters are summarized in Table 1.

As we have mentioned, it may be difficult to reproduce both the X-ray and optical light curves by a single emission component with constant microscopic parameters. In model A, we give weight to the optical light curve as shown in Figure 10. For $t\lesssim 10^{4}$ s, the model flux in X-ray is dimmer than the observed one. Another emission component such as the reverse shock may be required to agree with the early X-ray light curve. The small value of the parameter $\epsilon_{\rm e}$ leads to a large value of $E_{0}$. If we adopt a higher $\epsilon_{\rm e}$, the high initial $\gamma_{\rm m}$ makes the peak time of the optical light curve delayed compared to the observed onset time (see Model C). As the spectrum for 240–270 s in Figure 10 shows, to generate such bright synchrotron flux as the early X-ray data indicate, $\varepsilon_{\rm m}$ should be in the X-ray energy range. However, such a high $\varepsilon_{\rm m}$ contradicts the decaying flux of the optical emission at this stage.

The 0.1 GeV light curve is well reproduced by our model. The model curve shows a smooth-power-law-like behavior. Around 0.1 GeV, both the synchrotron and SSC emissions contribute. Even if we artificially turn off the SSC emission (see the green dashed line in the left panel of Figure 10), the 0.1 GeV emission due to synchrotron yields a single power-law light curve until $10^{5}$ s. Therefore, the 0.1 GeV range is not so ideal to find the switching from synchrotron to SSC in the light curve.

The detections of the 95 and 32 GeV photons are not explained by synchrotron emission, as shown in the right panel of Figure 10 (see dashed lines for the model without SSC, and red and blue vertical lines for the energies of the detected photons). Even in the early period of $t\sim 200$–$300$ s, photons above 10 GeV are emitted via SSC in model A. The predicted 0.1 TeV light curve (purple) is also smooth and lasts a long time. Even at $t=10^{4}$ s, the flux at 0.1 TeV is about $5\times 10^{-10}~{}\mbox{erg}~{}\mbox{cm}^{-2}~{}\mbox{s}^{-1}$, which can be detected with CTA with a time resolution of a few hundred seconds (Funk & Hinton, 2013; Inoue et al., 2013).

While we have adopted the Bohm factor as unity in Model A, the thin black line in the right panel of Figure 10 shows the 0.1 GeV light curve with $\xi=100$. In this conservative model, the dominant emission process in the 0.1 GeV range is replaced from synchrotron to SSC in the later phase. However, the 0.1 GeV light curve is smooth even in this case.

Next, giving weight to the X-ray light curve rather than the optical one, let us try to find an acceptable model. The steeper slope of the X-ray light curve leads to a larger $p$, which makes the spectrum softer. Model B is an example of our results that agree with the observed X-ray and 0.1 GeV light curves (see left panel of Figure 11). As the green dashed line indicates, the emission at 0.1 GeV for $>10^{2}$ s is dominated by SSC in model B. The soft spectrum results in brighter optical flux than observed. Maselli et al. (2014) claimed that the optical extinction is negligible from the SED analysis. The relatively dim optical fluxes and the steep X-ray decay seem difficult to explain simultaneously by a single source model with constant microscopic parameters. However, interestingly, this model omitting the effect of the jet break is consistent with the X-ray light curve as far as $\sim 10^{8}$ s, though no signature of the jet break challenges the standard afterglow model and implies a very large energy release for this GRB. As De Pasquale et al. (2016) pointed out, the previous complex models (Kouveliotou et al., 2013; Laskar et al., 2013; Panaitescu et al., 2013; Maselli et al., 2014; Perley et al., 2014; van der Horst et al., 2014) have difficulties reconciling the observed long-lasting X-ray emission. Although we focus on the early afterglow rather than the late one, the physical parameters may be close to those of model B in the late phase. However, the predicted radio flux at 6.8 GHz is significantly brighter than the observed flux (Maselli et al., 2014).

Another example is model C, whose light curves are shown in the right panel of Figure 11. By increasing $\epsilon_{\rm e}$ and decreasing $n_{0}$, the required energy $E_{0}$ is drastically suppressed compared to models A and B. In model C, $\varepsilon_{\rm m}$ is kept higher than the optical range for a long time. The resultant optical light curve shows a late peak time, which seems inconsistent with the simple power-law decay of the observed light curve. The high $\gamma_{\rm m}$ in this model leads to the dominant contribution of the synchrotron emission at 0.1 GeV as far as a few times $10^{3}$ s, after which the contribution of SSC emission modulates the 0.1 GeV light curve. This deviation from single power-law in the 0.1 GeV light is still within the observational errors.

All the models in Table 1 have a very small value of $\epsilon_{B}$, which agrees with the results of recent studies (Kumar & Barniol Duran, 2009; Lemoine et al., 2013; Santana et al., 2014; Beniamini et al., 2015; Zhang et al., 2015). In spite of the small $\epsilon_{B}$, the discrepancy factors in the X-ray flux are of the order of unity (see Table 1). The initial magnetic fields of models A–C are 0.14, 0.49, and 0.16 G, respectively. A shock-compressed CSM magnetic field is only $4\Gamma_{0}B_{\rm CSM}\simeq 1.4(\Gamma_{0}/350)(B_{\rm CSM}/\mu\mbox{G})$ mG. Even for those small $\epsilon_{B}$, an amplification mechanism of the magnetic field is required (see, e.g. Barniol Duran, 2014).

6. Conclusions

In order to simulate the GRB afterglow emission, we have calculated the temporal evolutions of the energy distributions of electrons and photons in the shell relativistically propagating in the ISM. Physical processes such as the deceleration of the shell, photon escape, adiabatic cooling, and transformations of observables into the observer frame are consistently dealt with in our numerical code. Given the initial Lorentz factor $\Gamma_{0}$, the onset time of the afterglow in our results is significantly earlier than the previous analytical estimate. The uncertainty in the initial Lorentz factors obtained from the onset time may be larger than previously thought. When we mimic the radiative case by adopting an extreme value $\epsilon_{\rm e}=0.9$, the results show $\Gamma\propto R^{-2}$ and $F\propto t_{\rm obs}^{-1.4}$, which are significantly different from the conventional formulae.

In the fast cooling case, our results show that the electron spectrum for $\gamma_{\rm c}<\gamma<\gamma_{\rm m}$ is significantly curved and harder than the analytical estimate owing to the evolution of the injection rate. The spectral shape is highly curved around the typical energies $\varepsilon_{\rm m}$ and $\varepsilon_{\rm c}$. While the peak flux at $\varepsilon\sim\varepsilon_{\rm m}$ or $\varepsilon_{\rm c}$ is lower than the analytical estimate with the broken-power-law approximation, the discrepancy of the X-ray flux with the analytical synchrotron formula is not so large. The total energy obtained by fitting the observed light curves with the analytical formula may be underestimated by a factor of three or less. However, as Beniamini et al. (2015, 2016) pointed out, if we misunderstand that the X-ray-emitting electrons are in the fast cooling regime despite $\epsilon_{B}\ll 10^{-3}$, the total energy can be highly underestimated. This may resolve the high-efficiency problem in the prompt emission (Lloyd-Ronning & Zhang, 2004).

Our results show that even if the emission mechanism is switching from synchrotron to SSC, the gamma-ray light curves can be a smooth power-law, especially for the electron index of $p\simeq 2$–2.5. Note that we have not intentionally adjusted the parameters to suppress the light-curve signature of the switching from synchrotron to SSC. In most cases with fiducial parameter sets, it is difficult to find the time at which SSC starts contributing from only light curves.

Given the electron spectral index $p$, the SSC contribution makes the photon spectrum slightly harder than the expectation from the synchrotron formula. We have tested 54 models changing the parameters. The numerically obtained spectral index and decay index are scattered, but distribute along the analytical closure relation. To explain GRBs whose indices largely deviate from the closure relation, the evolution of the microscopic parameters may be required.

With our method, we have fitted the light curves of GRB 130427A, in which high-energy photons beyond the synchrotron limit were detected. Although our single-source model with constant microscopic parameters does not reproduce all the observed behaviors in multiple wavelengths, the combination of the synchrotron and SSC emissions from the external shock can consistently explain the smooth 0.1 GeV light curve and the detections of 95 and 32 GeV photons at $t=244$ s and 34,400 s, respectively. As long as $\epsilon_{B}\ll 1$, as the recent studies suggested, 10–100 GeV SSC emission will be expected to be detected with CTA (see Vurm & Beloborodov, 2016, for a conservative estimate of the detection rate).

First, we appreciate the valuable advise by the anonymous referee. This work is supported by Grants-in-Aid for Scientific Research nos. 15K05069, 16K05291 (K.A.), and 15K05080 (Y.F.) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan.