#  Multiwavelength Modeling for the Shallow Decay Phase of Gamma-Ray Burst
Afterglows

Katsuaki Asano  Institute for Cosmic Ray Research, The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan

###### Abstract

We simulate the emission in the shallow decay phase of gamma-ray burst
afterglows using a time-dependent code. We test four models: the energy
injection model, evolving the injection efficiency of non-thermal electrons,
evolving the amplification of the magnetic field, and the wind model with a
relatively low bulk Lorentz factor. All of the four models can reproduce the
typical X-ray afterglow lightcurve. The spectral shape depends on not only the
parameter values at the time corresponding to the observer time but also the
past evolution of the parameters. The model differences appear in the
evolution of the broadband spectrum, especially in the inverse Compton
component. Future gamma-ray observations with imaging atmospheric Cherenkov
telescopes such as CTA will reveal the mechanism of the shallow decay phase.

Gamma-ray bursts (629); Non-thermal radiation sources (1119); Shocks (2086)

††journal: ApJ

##  1 Introduction

A large fraction of gamma-ray bursts (GRBs) show X-ray afterglows with a decay
index of $<0.5$ during initial $\sim 10^{3}$ s (Nousek et al., 2006; Liang et
al., 2007). This shallow decay phase (or plateau phase) is one of the unsolved
problems in the GRB afterglow physics. One of the most frequently discussed
models is the energy injection model (Zhang & Mészáros, 2001; Zhang et al.,
2006; Granot & Kumar, 2006), in which the energy in the forward shock region
increases with time, so that the decay of the bulk Lorentz factor becomes
slower than the standard Blandford–McKee solution (Blandford & McKee, 1976).
However, other type of models can be possible; for example, the two-component
jets or anisotropic jet structure (Eichler & Granot, 2006; Toma et al., 2006;
Granot et al., 2006), off-axis jet (Beniamini et al., 2020), evolving
microphysical parameters (Fan & Piran, 2006; Panaitescu et al., 2006; Ioka et
al., 2006), thick shell model (Leventis et al., 2014), and the decreasing
circumstellar density (wind model, Shen & Matzner, 2012; Dereli-Bégué et al.,
2022).

Optical afterglows frequently show a single power-law decay (Panaitescu &
Vestrand, 2008; Li et al., 2015), while the simultaneous X-ray afterglow is
shallower. The early optical afterglow may be dominated by the reverse shock
emission, thus optical observations may be not an ideal method to investigate
the mechanism of the shallow decay phase.

GeV–TeV gamma-ray observations can provide clues to the shallow decay phase.
GeV afterglows in initial $10^{3}$–$10^{4}$ s are detected with Fermi-LAT for
a fraction of GRBs (Ajello et al., 2019). However, Yamazaki et al. (2020)
claimed that GRBs detected with Fermi-LAT tend not to show the X-ray shallow
decay behavior. This may suggest that much deeper observations in the GeV band
are needed to detect GeV afterglows accompanying a shallow-decaying X-ray
afterglow.

Recently, imaging atmospheric Cherenkov telescopes have detected TeV
afterglows from several GRBs (MAGIC Collaboration et al., 2019; Abdalla et
al., 2019; H. E. S. S. Collaboration et al., 2021; Abe et al., 2024). The
detection of more than 64,000 photons above 0.2 TeV from GRB 221009A with
LHAASO (LHAASO Collaboration et al., 2023) is encouraging for future
observations of TeV afterglows.

In this paper, we simulate broadband emission from the shallow decay phase
with several models using our time-dependent code. The model differences
appear in the evolution of the photon spectra, especially in the inverse
Compton (IC) component. Future gamma-ray observations of GRBs will reveal the
mechanism of the shallow decay phase.

This paper is organized as follows. In §2, we explain our method to simulate
the afterglow emission. We attempt four models for the shallow decay phase,
which are presented in §3. The results for the four models are shown in §4. §5
is devoted to summary.

##  2 Method

In this paper, to simulate the evolution of GRB afterglow spectra, we adopt
the numerical code in Fukushima et al. (2017), which was also used in Asano et
al. (2020). The code is one-zone time-dependent to simulate the evolution of
the electron and photon energy distributions in the emission region, which is
the shocked circumstellar medium (CSM). The physical processes included in
this code are particle injection, synchrotron emission, synchrotron self-
absorption, inverse Compton emission, $\gamma\gamma$-absorption,
electron–positron pair creation, adiabatic cooling, and photon escape. With
the energy conservation law, the bulk Lorentz factor $\Gamma$ of the emission
region is developed with the radius $R$ from the central engine taking into
account mass loading and photon emission. The emission region is treated as a
uniform shell. The shock jump condition provides the density and the volume
(or equivalently width $W$) of the emission region under the one-zone
approximation.

Electrons are injected with a power-law index $p$ ($\dot{n}_{\rm
e}\propto\varepsilon_{\rm e}^{-p}$) with the energy fraction $\epsilon_{\rm
e}$ to the injected energy into the shell. Those two parameters regulate the
minimum Lorentz factor $\gamma_{\rm min}$ of electrons at injection. The
magnetic field evolves with the energy fraction $\epsilon_{B}$. See Fukushima
et al. (2017) for other details. The photon spectrum for an observer is
calculated consistently with the amount of escaped photons. Assuming a
spherically symmetric geometry, the code integrates photons escaped from the
emission surface taking into account the Doppler beaming and the surface
curvature to address the photon arrival time.

Most of the analytical formulae for GRB afterglow spectra (Sari et al., 1998;
Panaitescu & Kumar, 2000; Zhang & Mészáros, 2001) are based on the
approximation that the electron energy distribution is instantly determined by
the local quantities such as $\Gamma$ and the CSM density at that radius $R$,
irrespectively of the past evolution of the emission region. Our time-
dependent code is free from this instant control approximation and one-to-one
correspondence approximation between $R$ and the observer time $t_{\rm obs}$
(see e.g. Derishev & Piran, 2021).

In this paper, we change the photon escape process in the numerical code,
which affects the photon energy density. Isotropic photons in the shell escape
from both the front and rear surfaces. In Fukushima et al. (2017), the
boundary of the shell is assumed to be static in the shell rest frame, so that
the photon escape rate is $cN_{\gamma}/(2W)$, where $N_{\gamma}$ is the total
photon number in the shell. However, as the shock is relativistic in early
afterglows, the shock front in the shocked fluid rest frame propagates with a
speed of $c/3$. Only photons with the angle $\theta$ of $\cos\theta>1/3$ can
escape with the effective velocity of $c(\cos\theta-1/3)$. The rear boundary
is the contact discontinuity or rarefaction wave in actual cases, so the
boundary velocity may be non-relativistic. However, such an asymmetric
structure makes it difficult to estimate the photon density consistently in
the one-zone and isotropic approximation. For simplicity, we adopt the same
escape rate even for the rear boundary, then the escape rate becomes

$\displaystyle\dot{N}_{\gamma}=\frac{2N_{\gamma}}{9W}.$ (1)

Figure 1: Comparison of de-absorbed photon spectra (neglecting EBL absorption)
between the previous (dashed) and present (solid) methods for the photon
escape process. From top to bottom, the observed time is 3,000 s, 10,000 s,
and 30,000 s, respectively. The parameters are $z=0.5$, $E_{0}=10^{53}$ erg,
$\Gamma_{0}=300$, $n_{0}=1~{}\mbox{cm}^{-3}$, $\epsilon_{\rm e}=0.1$,
$\epsilon_{B}=10^{-3}$, and $p=2.2$.

We show the differences in photon spectra due to the different methods for the
photon escape process in Fig. 1. Here, we have assumed the standard ISM model
(constant CSM density) with the macroscopic parameters, the CSM density
$n_{0}$, the initial Lorentz factor $\Gamma_{0}$, and kinetic energy $E_{0}$.
Throughout this paper, we demonstrate the observed photon flux with a source
redshift $z=0.5$, neglecting absorption due to extra-galactic background light
(EBL). As discussed in Fukushima et al. (2017), our time-dependent numerical
code produces smoothly curved spectra, in which it is hard to identify the
cooling break $\nu_{\rm c}$ and the break frequency $\nu_{\rm m}$
corresponding to the typical emission frequency for electrons with
$\gamma=\gamma_{\rm min}$ in the conventional notation.

As shown in the figure, the slower photon escape slightly enhances the ratio
of the IC component to the synchrotron component. The differences in the
photon flux with the two methods are not larger than a factor of two.

##  3 Models for Shallow Decay

We consider four models in this paper to reproduce the shallow decay behavior
in GRB afterglows. The first one is the most promising and well-known model,
the energy injection model (Zhang & Mészáros, 2001; Zhang et al., 2006). We
test two other models in which the microscopic parameters evolve with time
(Panaitescu et al., 2006; Ioka et al., 2006). The last is the wind model with
a relatively low $\Gamma$ (Shen & Matzner, 2012; Dereli-Bégué et al., 2022).
In the first three models, one of the parameters, $E_{0}$, $\epsilon_{\rm e}$,
or $\epsilon_{B}$ evolves with the radius. After the shallow decay phase, the
final parameter values are common for those models as $E_{0}=10^{53}$ erg,
$\epsilon_{\rm e}=0.1$, and $\epsilon_{B}=10^{-3}$. The other parameters are
also common as $\Gamma_{0}=300$, $n_{0}=1~{}\mbox{cm}^{-3}$, and $p=2.2$ so
that the behaviors in the late stage should be the same.

Note that the parameter evolution in those models required for reproducing the
shallow decay lightcurves is not necessarily justified by a well-motivated
model. To phenomenologically reproduce the shallow decay behaviors, the
parameter evolution is adjusted.

###  3.1 Energy Injection (EI) Model

If the central engine acts longer than the prompt emission timescale or the
bulk Lorentz factor of ejecta decreases with time, the energy in the forward
shock region $E$ increases as the shock front decelerates. In most cases, the
X-ray band is above $\nu_{\rm c}$ and $\nu_{\rm m}$ in the shallow decay
phase. In this case, as discussed in Zhang & Mészáros (2001), if we
phenomenologically express $E\propto t^{s}\propto R^{s}$, the decay index is
approximated as $-(2p-2)/3$ or $-(p-2)/2$ for $s=1$ and 2, respectively. To
reproduce a shallower decay index than $-0.5$ with $p=2.2$, we adopt $s=2$,
which implies a constant luminosity $L_{\rm jet}\propto t^{0}$ or a
$\Gamma$-distribution of $dE/d\Gamma\propto\Gamma^{-5}$ in the long engine-
activity model and the decreasing-$\Gamma$ model, respectively.

Figure 2: Evolutions of $\Gamma\beta$ (black, left axis) and the total energy
injected into the CSM (blue, right axis) in the energy injection model. The
dashed lines are behaviors for the standard constant energy model with
$E_{0}=10^{53}$ erg. The parameters are $E_{\rm i}=10^{51}$ erg,
$\Gamma_{0}=300$, $n_{0}=1~{}\mbox{cm}^{-3}$, and $R_{\rm e}=3.2\times
10^{17}$ cm.

In this paper, introducing the initial and final energy as $E_{\rm i}$ and
$E_{0}$, respectively, we assume that the total energy increases with $R$ as

$\displaystyle E=E_{\rm i}+(E_{0}-E_{\rm i})\left(\frac{R}{R_{\rm
e}}\right)^{2},\quad\mbox{for}\quad R<R_{\rm e}.$ (2)

For $R\geq R_{\rm e}$, the total energy is constant $E_{0}$.

Figure 3: Evolutions of the electron energy density in the evolving
$\epsilon_{\rm e}$ model (upper panel) and the magnetic energy density in the
evolving $\epsilon_{B}$ model (lower panel). The dashed lines are behaviors
for the standard constant microparameter model ($\epsilon_{\rm e}=10^{-1}$,
$\epsilon_{B}=10^{-3}$). The parameters are $E_{0}=10^{53}$ erg,
$\Gamma_{0}=300$, $n_{0}=1~{}\mbox{cm}^{-3}$, and $R_{\rm e}=3.2\times
10^{17}$ cm.

The evolutions of the bulk 4-velocity $\Gamma\beta$ and the energy injected
into the CSM with $E_{\rm i}=10^{51}$ erg and $R_{\rm e}=3.2\times 10^{17}$ cm
are shown in Fig. 2. In our simulations, the initial energy is held by the
cold ejecta, and the energy is gradually transferred to the CSM via the shock
as seen in Fig. 2. While the simple analytical estimate with $E\propto R^{2}$
leads to $\Gamma\propto R^{-1/2}$, the numerical results show $\Gamma\propto
R^{-0.8}$ for $R<R_{\rm e}$. In Fig. 2, the two break radii for the solid line
correspond to the radius where $E=E_{\rm i}$ and $R_{\rm e}$. The solid line
roughly connects the two points on the dashed line (constant-$E$ model) at
those two radii. Thus, the exact decay index of $\Gamma$ depends on the
initial energy and $\Gamma_{0}$. For $R>R_{\rm e}$, the bulk Lorentz factor
decays following the standard Blandford–Mckee solution $\Gamma\propto
R^{-3/2}$.

###  3.2 Evolving $\epsilon_{\rm e}$ (EE)

Adopting the usual instant control approximation, the flux above $\nu_{\rm c}$
and $\nu_{\rm m}$ decays as $F\propto\epsilon_{\rm
e}^{p-1}\epsilon_{B}^{(p-2)/4}t_{\rm obs}^{(2-3p)/4}$. The one-to-one
correspondence approximation with $\Gamma\propto R^{-3/2}$ implies $t_{\rm
obs}\propto R^{4}$. To reproduce the shallow decay with those approximations,
we can assume $\epsilon_{\rm e}\propto R^{2}$, which leads to $F\propto t_{\rm
obs}^{-p/4}$. On the other hand, for $\nu_{\rm c}>\nu_{\rm obs}>\nu_{\rm m}$,
$F\propto\epsilon_{\rm e}^{p-1}\epsilon_{B}^{(p+1)/4}t_{\rm obs}^{3(1-p)/4}$.
Even in this case, $\epsilon_{\rm e}\propto R^{2}$ leads to a shallow decay as
$F\propto t_{\rm obs}^{(1-p)/4}$. Here, we attempt an evolving $\epsilon_{\rm
e}$ model with a floor value as

$\displaystyle\epsilon_{\rm e}=\epsilon_{\rm i}+(\epsilon_{\rm
f}-\epsilon_{\rm i})\left(\frac{R}{R_{\rm e}}\right)^{2},\quad\mbox{for}\quad
R<R_{\rm e}.$ (3)

For $R\geq R_{\rm e}$, $\epsilon_{\rm e}$ is constant $\epsilon_{\rm f}$. We
adopt $\epsilon_{\rm f}=0.1$ and $\epsilon_{\rm i}=10^{-4}$. As shown in
Figure 3, the electron energy density $U_{\rm e}$ is initially suppressed, but
the density for $R>R_{\rm e}$ is almost the same as that for the standard
constant microparameter model ($\epsilon_{\rm e}=10^{-1}$,
$\epsilon_{B}=10^{-3}$). Since $U_{\rm e}\propto\epsilon_{\rm e}\Gamma^{2}$ in
the instant control approximation, we can expect $U_{\rm e}\propto R^{-1}$
just before $R_{\rm e}$. Our numerical result shows $U_{\rm e}\propto
R^{-0.9}$.

###  3.3 Evolving $\epsilon_{B}$ (EB)

In the model with an increasing $\epsilon_{B}$, the initial magnetic field is
relatively low so that $\nu_{\rm c}$ likely becomes significantly higher than
the X-ray frequency. For $\nu_{\rm c}>\nu_{\rm obs}>\nu_{\rm m}$, the choice
of $\epsilon_{B}\propto R^{4}$ is preferable for the shallow decay as
$F\propto t_{\rm obs}^{(2-p)/2}$. Similarly to the evolving $\epsilon_{\rm e}$
model, we adopt

$\displaystyle\epsilon_{B}=\epsilon_{\rm i}+(\epsilon_{\rm f}-\epsilon_{\rm
i})\left(\frac{R}{R_{\rm e}}\right)^{4},\quad\mbox{for}\quad R<R_{\rm e}.$ (4)

For $R\geq R_{\rm e}$, $\epsilon_{B}$ is constant $\epsilon_{\rm f}$. We adopt
$\epsilon_{\rm f}=10^{-3}$ and $\epsilon_{\rm i}=10^{-6}$. The evolution of
the magnetic energy density $U_{B}$ is shown in Figure 3. Similarly to the
case in the evolving $\epsilon_{\rm e}$ model, the behavior of $U_{B}$ for
$R>R_{\rm e}$ is almost the same as that for the standard constant microscopic
parameter model. While the analytical estimate gives
$U_{B}\propto\epsilon_{B}\Gamma^{2}\propto R$ just before $R_{\rm e}$, our
numerical result shows $U_{B}\propto R^{1.1}$.

Figure 4: Evolutions of $\Gamma\beta$ for the wind model (solid) and the
standard ISM model (dashed, same as the dashed black line in Figure 2). The
parameters for the wind model are $E_{0}=2\times 10^{53}$ erg,
$\Gamma_{0}=30$, and $A_{*}=0.3$.

###  3.4 Wind Model

The progenitor stars of GRBs are expected to lose their envelope by stellar
winds. The density profile in the wind model is conventionally parameterized
as

$\displaystyle\rho=5.0\times
10^{11}A_{*}\left(\frac{R}{1~{}\mbox{cm}}\right)^{-2}~{}\mbox{g}~{}\mbox{cm}^{-3},$
(5)

where $A_{*}=1$ corresponds to a mass loss rate of
$10^{-5}M_{\sun}~{}\mbox{yr}^{-1}$ with a wind velocity of
$10^{3}~{}\mbox{km}~{}\mbox{s}^{-1}$. With the mass $M(R)$ within a radius
$R$, the deceleration radius $R_{\rm dec}$ is estimated by
$E_{0}=\Gamma_{0}^{2}M(R_{\rm dec})c^{2}$. In this wind model, we obtain

$\displaystyle R_{\rm dec}\simeq 2.0\times
10^{16}A_{*}^{-1}\left(\frac{E_{0}}{10^{53}~{}\mbox{erg}}\right)\left(\frac{\Gamma_{0}}{30}\right)^{-2}\mbox{cm}.$
(6)

The end of the coasting phase for an observer (see e.g. Fukushima et al.,
2017) is at

$\displaystyle t_{\rm dec}$ $\displaystyle\simeq$
$\displaystyle(1+z)\frac{R_{\rm dec}}{5c\Gamma_{0}^{2}}$ (7)
$\displaystyle\simeq$ $\displaystyle
150(1+z)A_{*}^{-1}\left(\frac{E_{0}}{10^{53}~{}\mbox{erg}}\right)\left(\frac{\Gamma_{0}}{30}\right)^{-4}\mbox{s}.$
(8)

During the coasting phase, the analytical estimate leads to $F\propto t_{\rm
obs}^{(2-p)/2}$ above the frequencies $\nu_{\rm m}$ and $\nu_{\rm c}$. Thus,
the coasting phase in the wind model can be regarded as a shallow decay phase
for an observer.

To elongate the shallow decay phase, we adopt $A_{*}=0.3$, $E_{0}=2\times
10^{53}$ erg, and $\Gamma_{0}=30$, which implies $t_{\rm dec}\simeq 1500$ s
for $z=0.5$. This parameter choice is determined to generate a similar flux
and lightcurves to those in the other models after the shallow decay phase. As
shown in Figure 4, the bulk Lorentz factor gradually decreases with radius,
and asymptotically behaves as $\Gamma\propto R^{-0.5}$ in the late phase.

##  4 Results

Figure 5 shows X-ray lightcurves for all the models. Those models well
reproduce the shallow decay phase. Considering observational errors, those
four lightcurves are almost identical. In the energy injection (EI) model, the
Lorentz factor at $R=R_{\rm e}$ is 22.7. The sharp break in the X-ray
lightcurve at $t_{\rm obs}=4000$ s implies a factor $K\simeq 8$ for the
correspondence

$\displaystyle t_{\rm obs}$ $\displaystyle\simeq$
$\displaystyle(1+z)\frac{R}{Kc\Gamma^{2}},$ (9)

while the peak time of the afterglow (end of the coasting phase) in the
standard model is consistent with $K\simeq 5$ (Fukushima et al., 2017). Note
that the definition of $\Gamma$ is not the shock front Lorentz factor
differently from Derishev & Piran (2021).

The models with the same $R_{\rm e}$, the evolving $\epsilon_{\rm e}$ (EE)
model, and the evolving $\epsilon_{B}$ (EB) model, show a more smooth
transition from the shallow to the normal decay phase. The lightcurve of the
wind model with our parameter choice is similar to that for the EE model. For
$R>R_{\rm e}$, the parameter evolutions for the constant CSM models (EI, EE,
EB) are the same so that all the lightcurves are almost the same. In the
analytical prediction, the decay index for the observation frequency above
$\nu_{\rm c}$ and $\nu_{\rm m}$ are common for all the four models. As $E_{0}$
is higher in the wind model, the flux in the wind model is slightly brighter
than the other lightcurves in the normal decay phase.

Figure 5: X-ray lightcurves (1 keV) for the energy injection model (EI, solid
black), the evolving $\epsilon_{\rm e}$ model (EE, blue), the evolving
$\epsilon_{B}$ model (EB, green), and the wind model (red). The dashed line is
for the standard model.

Figure 6: Evolutions of the afterglow spectra from 100 s to 10000 s for the
energy injection model (EI), the evolving $\epsilon_{\rm e}$ model (EE), the
evolving $\epsilon_{B}$ model (EB), and the wind model.

Although the X-ray flux in the EI model at the end of the shallow decay phase
is slightly brighter than those in the other models, considering observational
errors, the four lightcurves are almost identical.

Though the X-ray lightcurves are similar, a characteristic behavior in the
evolution of the broadband spectrum is seen for each model (Figure 6).
Especially, the differences in the evolution and dominancy of the IC component
are distinctly shown in Figure 6. Even with the same $\epsilon_{\rm e}$ and
$\epsilon_{B}$ in the EI and wind models, the gamma-ray flux for the wind
model is significantly brighter than the flux for the EI model in the early
stage. This is mainly due to the Klein-Nishina effect in the EI model, in
which $\gamma_{\rm m}\propto\Gamma$ is higher than that in the wind model. As
a result, the spectral peak of the IC component in the wind model ($\sim$ GeV)
is lower than those in the other models ($\sim$ TeV). In the EE model, the
high initial value of $\epsilon_{B}/\epsilon_{\rm e}$ suppresses the IC
component. In the EB model, the initial high value of $\epsilon_{\rm e}$ leads
to the high flux of the IC component. Differently from the other models, the
resultant high $\nu_{\rm c}$ causes a hard X-ray spectrum. Though such hard
X-ray spectra in the shallow decay phase are not the majority, several
examples have been found (Willingale et al., 2007) including the low
luminosity GRB 171205A (D’Elia et al., 2018; Izzo et al., 2019).

Figure 7: Electron energy evolution for the standard (dashed black), EE
(blue), and EB (green) models. The corresponding radii are $1.8\times 10^{15}$
cm, $1.8\times 10^{16}$ cm, $5.6\times 10^{16}$ cm, and $1.8\times 10^{17}$ cm
from thin to thick lines.

Figure 8: Electron energy distribution at $R=R_{\rm e}$ ($3.2\times 10^{17}$
cm, left) and $R=10^{18}$ cm (right) for the standard (dashed black), EI
(solid black), EE (blue), and EB (green) models.

In the EI and EB models, the spectral softening is the main reason for the
X-ray shallow decay, while the spectral peak significantly decreases with
time. Especially in the EB model, the flux in a higher energy band than X-ray
evolves rapidly. On the other hand, the spectral evolutions of the synchrotron
component in the EE and Wind models are not prominent. The X-ray shallow decay
is due to the gradual decay of the synchrotron components maintaining their
spectral shapes in those cases. This gradual evolution in the EE and wind
model is reflected in the X-ray lightcurves (Figure 5).

The evolution of the electron energy distributions for the standard, EE, and
EB models is shown in Figure 7. In the EE model, the initial low
$\epsilon_{\rm e}$ leads to a lower $\gamma_{\rm min}$. With increasing
$\epsilon_{\rm e}$, $\gamma_{\rm min}$ also increases, while the cooling break
energy decreases with time. In the EB model, the initial low $\epsilon_{B}$
suppresses the radiative cooling so that the electron spectrum evolves
maintaining the hard spectrum.

At $R=R_{\rm e}$, the microscopic parameters and bulk Lorentz factor are
common for the EI, EE, and EB models. In this case, the instant control
approximation leads to the same electron energy distribution for $R\geq R_{\rm
e}$. However, as shown in Figure 8, our numerical results show different
electron distributions. In the standard model (dashed line), the radiative
energy loss in the early stage appears as the lower electron energy density
than that in the EI model, because of the late energy injection at $R\sim
R_{\rm e}$. The low radiative efficiency in the EB model leads to a harder
spectrum than the other spectra. The low energy density of the EE model is due
to the small $\epsilon_{\rm e}$ in the early stage. At $R\simeq 3\times R_{\rm
e}$, the electron distributions in the three models are almost identical. At
this radius, the information on the initial history of the electron energy
distribution almost disappeared.

Figure 9: Photon spectra at the end of the shallow decay phase ($t_{\rm
obs}=4000$s) for the standard (dashed black), EI (solid black), EE (blue), and
EB (green) models.

The difference in the electron energy distribution is reflected in the photon
spectra (Figure 9). According to the differences in the electron energy
distribution at $R=R_{\rm e}$ (Figure 8), the fluxes are brighter in the order
of the EE, standard, and EI models. The relatively high ratio of
$\epsilon_{B}/\epsilon_{\rm e}$ in the EE model suppresses the IC component.
As shown in Figure 3, the magnetic field of the EB model at $R=R_{\rm e}$ is
still lower than that in the standard model. In our one-zone treatment, the
initial low-energy injection of the magnetic energy in the EB model causes the
lower magnetic field at $R=R_{\rm e}$. Though the electron energy density of
the EI model at $R=R_{\rm e}$ is close to that of the EB model (Figure 8), the
photon flux of the EB model is dimmer than that of the EI model.

Figure 10: Photon lightcurves in optical (left), 0.1 GeV (middle), and TeV
(right) bands for the standard (dashed black), EI (black), EE (blue), and EB
(green) models.

If we obtain an evolution of the broadband spectrum as shown in Figure 6, the
discrimination of the models would be easy. However, in most cases, we obtain
only lightcurves in limited several frequencies. In Figure 10, we plot optical
and gamma-ray lightcurves for each model. The optical flux for the EI model
characteristically increases until the end of the energy injection. Such
lightcurves are very rare, though there are several examples of increasing
optical fluxes such as GRB 080710 (Krühler et al., 2009), whose optical peak
time is $\sim 2\times 10^{3}$ s.

A fraction of shallow-decay GRBs show a shallowly decaying optical lightcurve
as well (see e.g. Li et al., 2015). Except for the EI model, the model
lightcurves seem consistent with such a shallow decay optical lightcurve,
though the fluxes are different. Even for the standard model, the optical
lightcurve shows a similar behavior to the shallow decay one. This is because
the optical frequency is below $\nu_{\rm m}$ and $\nu_{\rm c}$. The analytical
estimate predicts a growing flux $\propto t^{1/2}_{\rm obs}$ for the slow
cooling case in this frequency range (Sari et al., 1998; Fukushima et al.,
2017). The gradual switch from this growing phase to the normal decay phase
makes such a shallow-decay-like lightcurve. The optical break time is not
necessarily simultaneous to the X-ray one: some show an earlier break, while
some show a later break (Li et al., 2015). The optical lightcurves for the
standard and wind models break earlier, and those for the EE and EB models
break later. The variety of the optical break timings can be attributed to the
diversity of the shallow decay models.

However, a large fraction of optical lightcurves are consistent with a single
power-law decay (Panaitescu & Vestrand, 2008; Li et al., 2015), while the
X-ray lightcurve shows a shallow decay behavior. In general, the simultaneous
modeling for the X-ray and optical lightcurves with only the forward-shock
emission is difficult (e.g. Fan & Piran, 2006; Rykoff et al., 2009). The most
promising model for the early optical lightcurves is two-component models such
as the forward–reverse shock model (e.g. Perley et al., 2014; Japelj et al.,
2014; Fraija et al., 2019) and/or spine–sheath model (e.g. Racusin et al.,
2008; Sato et al., 2023). In such cases, the contaminated optical lightcurves
are not an ideal tool for model discrimination. Therefore, we focus on the
gamma-ray afterglow, which may be not polluted by the reverse shock emission
(however, see Arimoto et al., 2024).

As shown in Figure 6, except for the wind model, the synchrotron component
dominates the 0.1 GeV emission. Figure 10 shows that the 0.1 GeV lightcurves
should show shallow decay except for the EB model. At $t_{\rm obs}=100$ s, the
EB and wind models predict a significantly brighter 0.1 GeV flux than the EI
model.

More significant differences are seen in the TeV lightcurves. All the
lightcurves show shallow decay, but the flux levels are significantly
different during the shallow decay in those models. From only TeV lightcurves,
it may be difficult to distinguish the EI and Wind models, whose TeV fluxes
are similar. However, as shown in Figure 6, the Wind model predicts a softer
spectrum in the TeV energy range. The measurement of the spectral index, which
may be disturbed by the EBL absorption though, can distinguish the EI and wind
models.

##  5 Summary

We have simulated the emission from the forward shock with the time-dependent
code to reproduce the X-ray shallow decay phases. The energy injection model,
evolving $\epsilon_{\rm e}$ model, evolving $\epsilon_{B}$ model, and wind
model are considered. Our time-dependent treatment shows that even after the
microscopic parameters are common, different parameter evolutions in the early
stage affect the spectral shape at the transition phase from shallow decay to
normal decay.

The model difference is clearly seen in the evolutions of the broadband
spectra. Especially from the evolutions of the IC component, we can
distinguish the models. In the energy injection model, both the synchrotron
and IC fluxes do not evolve significantly during the shallow decay phase. If
the parameter $\epsilon_{\rm e}$ evolves with radius, the IC component is
suppressed in the early stage, and its flux does not largely evolve. In the
increasing $\epsilon_{B}$ model, both the synchrotron and IC components
drastically decrease with time, except for the low energy range (X-ray and
optical). In the wind model, the IC component dominates, and its spectral peak
is lower than those in the other models.

We expect that future observations with CTA will provide a hint for the
mechanism of the enigmatic shallow decay phase.

This work is supported by the joint research program of the Institute for
Cosmic Ray Research (ICRR), the University of Tokyo, and KAKENHI No. 22K03684,
23H04899, and 24H00025.

## References

  * Abdalla et al. (2019) Abdalla, H., Adam, R., Aharonian, F., et al. 2019, Nature, 575, 464, doi: 10.1038/s41586-019-1743-9
  * Abe et al. (2024) Abe, H., Abe, S., Acciari, V. A., et al. 2024, MNRAS, 527, 5856, doi: 10.1093/mnras/stad2958
  * Ajello et al. (2019) Ajello, M., Arimoto, M., Axelsson, M., et al. 2019, ApJ, 878, 52, doi: 10.3847/1538-4357/ab1d4e
  * Arimoto et al. (2024) Arimoto, M., Asano, K., Kawabata, K. S., et al. 2024, Nature Astronomy, 8, 134, doi: 10.1038/s41550-023-02119-1
  * Asano et al. (2020) Asano, K., Murase, K., & Toma, K. 2020, ApJ, 905, 105, doi: 10.3847/1538-4357/abc82c
  * Beniamini et al. (2020) Beniamini, P., Duque, R., Daigne, F., & Mochkovitch, R. 2020, MNRAS, 492, 2847, doi: 10.1093/mnras/staa070
  * Blandford & McKee (1976) Blandford, R. D., & McKee, C. F. 1976, Physics of Fluids, 19, 1130, doi: 10.1063/1.861619
  * D’Elia et al. (2018) D’Elia, V., Campana, S., D’Aì, A., et al. 2018, A&A, 619, A66, doi: 10.1051/0004-6361/201833847
  * Dereli-Bégué et al. (2022) Dereli-Bégué, H., Pe’er, A., Ryde, F., et al. 2022, Nature Communications, 13, 5611, doi: 10.1038/s41467-022-32881-1
  * Derishev & Piran (2021) Derishev, E., & Piran, T. 2021, ApJ, 923, 135, doi: 10.3847/1538-4357/ac2dec
  * Eichler & Granot (2006) Eichler, D., & Granot, J. 2006, ApJ, 641, L5, doi: 10.1086/503667
  * Fan & Piran (2006) Fan, Y., & Piran, T. 2006, MNRAS, 369, 197, doi: 10.1111/j.1365-2966.2006.10280.x
  * Fraija et al. (2019) Fraija, N., Dichiara, S., Pedreira, A. C. C. d. E. S., et al. 2019, ApJ, 879, L26, doi: 10.3847/2041-8213/ab2ae4
  * Fukushima et al. (2017) Fukushima, T., To, S., Asano, K., & Fujita, Y. 2017, ApJ, 844, 92, doi: 10.3847/1538-4357/aa7b83
  * Granot et al. (2006) Granot, J., Königl, A., & Piran, T. 2006, MNRAS, 370, 1946, doi: 10.1111/j.1365-2966.2006.10621.x
  * Granot & Kumar (2006) Granot, J., & Kumar, P. 2006, MNRAS, 366, L13, doi: 10.1111/j.1745-3933.2005.00121.x
  * H. E. S. S. Collaboration et al. (2021) H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2021, Science, 372, 1081, doi: 10.1126/science.abe8560
  * Ioka et al. (2006) Ioka, K., Toma, K., Yamazaki, R., & Nakamura, T. 2006, A&A, 458, 7, doi: 10.1051/0004-6361:20064939
  * Izzo et al. (2019) Izzo, L., de Ugarte Postigo, A., Maeda, K., et al. 2019, Nature, 565, 324, doi: 10.1038/s41586-018-0826-3
  * Japelj et al. (2014) Japelj, J., Kopač, D., Kobayashi, S., et al. 2014, ApJ, 785, 84, doi: 10.1088/0004-637X/785/2/84
  * Krühler et al. (2009) Krühler, T., Greiner, J., Afonso, P., et al. 2009, A&A, 508, 593, doi: 10.1051/0004-6361/200912649
  * Leventis et al. (2014) Leventis, K., Wijers, R. A. M. J., & van der Horst, A. J. 2014, MNRAS, 437, 2448, doi: 10.1093/mnras/stt2055
  * LHAASO Collaboration et al. (2023) LHAASO Collaboration, Cao, Z., Aharonian, F., et al. 2023, Science, 380, 1390, doi: 10.1126/science.adg9328
  * Li et al. (2015) Li, L., Wu, X.-F., Huang, Y.-F., et al. 2015, ApJ, 805, 13, doi: 10.1088/0004-637X/805/1/13
  * Liang et al. (2007) Liang, E.-W., Zhang, B.-B., & Zhang, B. 2007, ApJ, 670, 565, doi: 10.1086/521870
  * MAGIC Collaboration et al. (2019) MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2019, Nature, 575, 459, doi: 10.1038/s41586-019-1754-6
  * Nousek et al. (2006) Nousek, J. A., Kouveliotou, C., Grupe, D., et al. 2006, ApJ, 642, 389, doi: 10.1086/500724
  * Panaitescu & Kumar (2000) Panaitescu, A., & Kumar, P. 2000, ApJ, 543, 66, doi: 10.1086/317090
  * Panaitescu et al. (2006) Panaitescu, A., Mészáros, P., Burrows, D., et al. 2006, MNRAS, 369, 2059, doi: 10.1111/j.1365-2966.2006.10453.x
  * Panaitescu & Vestrand (2008) Panaitescu, A., & Vestrand, W. T. 2008, MNRAS, 387, 497, doi: 10.1111/j.1365-2966.2008.13231.x
  * Perley et al. (2014) Perley, D. A., Cenko, S. B., Corsi, A., et al. 2014, ApJ, 781, 37, doi: 10.1088/0004-637X/781/1/37
  * Racusin et al. (2008) Racusin, J. L., Karpov, S. V., Sokolowski, M., et al. 2008, Nature, 455, 183, doi: 10.1038/nature07270
  * Rykoff et al. (2009) Rykoff, E. S., Aharonian, F., Akerlof, C. W., et al. 2009, ApJ, 702, 489, doi: 10.1088/0004-637X/702/1/489
  * Sari et al. (1998) Sari, R., Piran, T., & Narayan, R. 1998, ApJ, 497, L17, doi: 10.1086/311269
  * Sato et al. (2023) Sato, Y., Murase, K., Ohira, Y., & Yamazaki, R. 2023, MNRAS, 522, L56, doi: 10.1093/mnrasl/slad038
  * Shen & Matzner (2012) Shen, R., & Matzner, C. D. 2012, ApJ, 744, 36, doi: 10.1088/0004-637X/744/1/36
  * Toma et al. (2006) Toma, K., Ioka, K., Yamazaki, R., & Nakamura, T. 2006, ApJ, 640, L139, doi: 10.1086/503384
  * Willingale et al. (2007) Willingale, R., O’Brien, P. T., Osborne, J. P., et al. 2007, ApJ, 662, 1093, doi: 10.1086/517989
  * Yamazaki et al. (2020) Yamazaki, R., Sato, Y., Sakamoto, T., & Serino, M. 2020, MNRAS, 494, 5259, doi: 10.1093/mnras/staa1095
  * Zhang et al. (2006) Zhang, B., Fan, Y. Z., Dyks, J., et al. 2006, ApJ, 642, 354, doi: 10.1086/500723
  * Zhang & Mészáros (2001) Zhang, B., & Mészáros, P. 2001, ApJ, 552, L35, doi: 10.1086/320255