The afterglow emission from a stratified jet in GRB 170817A

GRB 170817A is known as the electromagnetic counterpart of the first binary neutron stars merger event. It was triggered by the Fermi Gamma-Ray Burst Monitor ∼ 1.7 s after the triggered of gravitational waves (GW 170817) (e.g., Abbott et al. 2017a,b). Later, the X-ray, radio, and optical afterglow of GRB 170817A were detected in sequence on ∼ 9, ∼16, and ∼ 110 days post-merger (e.g., Alexander et al. 2017; Haggard et al. 2017; Hallinan et al. 2017; Lyman et al. 2018; Margutti et al. 2017,2018; Troja et al. 2017). Actually, the related optical emission was detected earlier at ∼ 15 days post-merger, but it was thought to be the kilonova emission (e.g., Abbott et al. 2017b; Arcavi et al. 2017; Kasliwal et al. 2017; Pian et al. 2017). GRB 170817A is an unconventional short gamma-ray burst (SGRB) which is reflected in two aspects. First, the prompt emission is low-luminosity and the isotropic equivalent energy is E_{γ, iso} ∼ 6 × 10⁴⁶ erg (e.g., Margutti et al. 2018). Such an isotropic equivalent energy is ∼ four orders of magnitude lower than the classical SGRBs for which the typical energy is E_{γ, iso}∼10⁵⁰ – 10⁵² erg (Berger 2014; Fong et al. 2015). Second, the afterglow light curves of the burst have a super-long and slowly rising phase (with F_ν ∝ T^0.8) with a duration of ∼ 160 days (Mooley et al. 2018a).

Given the peculiar properties of the burst, there are mainly four potential scenarios. First, a classical top-hat jet (THJ) is viewed on-axis. This scenario needs a low-energy jet to produce the low-luminosity gamma-ray emission in the prompt phase, but such a low-energy jet is unlikely to break out the ejecta whose mass is suggested to be ≈ 0.05M_⊙ by the observed ultraviolet-optical-infrared counterpart (Kasliwal et al. 2017). In addition, this scenario cannot produce the slowly rising phase of the afterglow. Second, a classical THJ is viewed off-axis. The slowly rising phase of afterglow has ruled out this scenario, since an off-axis observed THJ will lead to a rise with the slope ≳ 3, much steeper than the observed F_ν ∝ T^0.8 (e.g., Lazzati et al. 2018; Margutti et al. 2017; Mooley et al. 2018a).

The other two scenarios involve the radial and angular structures of ejecta and, third, a quasi-spherical shell with an energy injection model (hereafter 'QSSEI'). The QSSEI model has a radial profile of energy, which is distributed as E(> u) ∝ u^−k (e.g., Gill & Granot 2018; Mooley et al. 2018a; Huang & Li 2018a), where u = Γ β is the proper velocity of the ejecta and k is the energy injection index. Such energy distribution indicates that there is more energy residing in the slower ejecta. As the fastest ejecta will gradually be decelerated by the external medium, the slower and more energetic ejecta will catch up with it and energize it (e.g., Sari & Mészáros 2000; Nakamura & Shigeyama 2006; Gill & Granot 2018). This scenario can explain the slowly rising phase of the afterglow but it is hard to produce a sharply decaying observed flux density as F_ν ∝ T^−2.2 (Mooley et al. 2018b,c). Moreover, by taking a global network of 32 radio telescopes, Ghirlanda et al. (2019) has constrained the apparent source size to be smaller than 2.5 milliarcseconds at 90% confidence level and the isotropic outflow scenario has been ruled out.

The fourth scenario is a wide-angle mildly relativistic cocoon with a successful off-axis jet (e.g., Kasliwal et al. 2017; Mooley et al. 2018c). The scenario consists of an energetic narrow core (with Lorentz factor Γ ≳ 100) and a less energetic wide-angle cocoon (with Lorentz factor Γ ≲ 10) with an angular profile of energy. The line of sight is off-axis of the jet core but within the cocoon. With the decelerating of the outflow, the visible area of the outflow is increasing and including more and more of the energetic region. Hence, the observed flux will gradually rise until the emission of the jet core comes into our view entirely, and then the observed flux drops sharply. Moreover, the angular displacement of the radio flux centroid is consistent with the observed for this scenario (Mooley et al. 2018c).

GRB 170817A appears to be the case of the fourth scenario. However, by studying the observed distribution of E_X /E_γ (the ratio of the isotropic equivalent early X-ray afterglow to prompt γ-ray energy), Beniamini & Nakar (2019) found that in order to be consistent with the observations, most of observed (long) GRB should be from a narrow region around the core if the jet has angular structure. Thus GRB 170817A is a unique burst, which could be from an off-core region. Given the uniqueness of the burst, it is worth exploring other possibilities besides the angular structure model. In this paper, we consider a stratified jet (with radial speed structure) model to explain the afterglow data of GRB 170817A. The stratified jet model can also interpret the early shallow decay on timescales of 10² – 10⁵ s in X-ray afterglows. Such a radial speed structure may originate from the central engine activities (Rees & Mészáros 1998; Sari & Mészáros 2000; Nakamura & Shigeyama 2006; Gill & Granot 2018). Both of the lateral spreading and non-spreading jets are considered. The paper is organized as follows. The model and methods are described in Section 2, which are including the model description, the dynamic evolution of the afterglow shock, the flux calculation, the analytical light curves from a stratified jet, and the estimation of the source size and the angular displacement of flux centroid. Fitting results and analysis are given in Section 3. Finally, we make a conclusion and discussion of our results in Section 4.

2.1. Model

2.2. Dynamic Evolution of the Afterglow Shock

The dynamic evolution of the jet directly determine the afterglow light curve. In the ultrarelativistic phase, the jet opening angle nearly remains a constant. But the lateral expansion becomes significant when the jet is decelerated down to Γ∼ 1/θ_j (Rhoads 1999), or later (e.g., van Eerten et al. 2010), where θ_j is the half opening angle of the jet. The jet opening angle evolution can be given by (Granot & Piran 2012)

$$\begin{eqnarray}&&\frac{d{\theta }_{j}}{dR}=\frac{1}{{\Gamma }^{2}{\theta }_{j}R},\end{eqnarray} \tag{ 1 }$$

where R is the radius of the jet. For a non-spreading jet, we have $\frac{d{\theta }_{j}}{dR}=0.$

The evolution of the afterglow shock's swept-up mass m can be described by

$$\begin{eqnarray}&&\frac{dm}{dR}=2\pi (1-\cos {\theta }_{j}){R}^{2}n{m}_{p},\end{eqnarray} \tag{ 2 }$$

where n is the circumburst medium number density and m_p is the proton mass. We adopt a thin shell approximation in this work. In the radially structured ejecta model, ∼ 90% of the injecting energy is shared by the forward-shocked medium and ∼ 10% is shared by the reverse-shocked material, and the reverse shock is Newtonian or trans-relativistic for the thin shell approximation case (Yu et al. 2007). Considering the energy fraction shared by the reverse-shocked material is very small, we take an approximation that the whole injecting energy is converted into the forward-shocked medium. Thus we obtain (Huang & Li 2018a)

$$\begin{eqnarray}&&m{c}^{2}[g(\Gamma)-1]\approx E(\gt \Gamma \beta)={E}_{0}{(\Gamma \beta)}^{-k},\end{eqnarray} \tag{ 3 }$$

where $g(\Gamma)=\frac{{\Gamma }^{3}+(\hat{\gamma }-1)({\Gamma }^{2}-1)(\Gamma -1)}{\Gamma }$, $\hat{\gamma }=(4\Gamma +1)/(3\Gamma)$ is the adiabatic index of the shocked medium, c is the light speed, and E(>Γ β) = E₀(Γ β)^−k is the kinetic energy of the ejecta. Note that Equation (3) is available for both relativistic and Newtonian cases. Combining Equations (1), (2) and (3), we can obtain the evolution of θ_j , m and Γ over R.

2.3. Flux Calculation

The observed flux for a THJ at a given observed time T can be given by (Granot et al. 1999; Woods & Loeb 1999)

$$\begin{eqnarray}\begin{array}{ll}{F}_{\nu }(T)= & \frac{1+z}{{D}_{L}^{2}}\displaystyle {\int }_{0}^{{\theta }_{v}+{\theta }_{j}}2\Delta \phi ({\theta }_{v},{\theta }_{j},\theta)\\ & \times {R}^{2}{\delta }^{3}{{I}^{\prime}}_{{\nu }^{\prime}}\cos \theta \sin \theta d\theta,\end{array}\end{eqnarray} \tag{ 4 }$$

where D_L is the luminosity distance of the burst, z is the redshift, δ = 1/Γ(1 – β cos θ) is the Doppler factor, θ_v is the viewing angle, and θ is the angle between the speed of a volume element within the jet and the viewing line. The observed time T is a function of Γ, R, and θ, which is

$$\begin{eqnarray}&&\frac{T}{1+z}=\displaystyle {\int }_{0}^{R}\frac{dr}{\beta c}-\frac{R{\rm{\cos }}\theta }{c}.\end{eqnarray} \tag{ 5 }$$

Δ ϕ is given by

$$\begin{eqnarray}\begin{array}{ll} & \Delta \phi =\\ & \left\{\begin{array}{ll}0 & ({\theta }_{v}\gt {\theta }_{j},\theta \lt {\theta }_{v}-{\theta }_{j})\\ \pi & ({\theta }_{v}\lt {\theta }_{j},\theta \lt {\theta }_{j}-{\theta }_{v})\\ {\cos }^{-1}(\frac{\cos {\theta }_{j}-\cos \theta \cos {\theta }_{v}}{\sin {\theta }_{v}\sin \theta }) & ({\rm{otherwise}})\end{array}.\right.\end{array}\end{eqnarray} \tag{ 6 }$$

${{I}^{\prime}}_{{\nu }^{\prime}}={{j}^{\prime}}_{{\nu }^{\prime}}{\Delta }^{\prime}$ is the radiation intensity in the jet comoving frame, and ${{j}^{\prime}}_{{\nu }^{\prime}}=\frac{{N}_{e}}{4\pi {R}^{2}{\Delta }^{\prime}}\frac{{{P}^{\prime}}_{{\nu }^{\prime},max}}{4\pi }f({\nu }^{\prime})$ is the comoving emissivity. Δ' is the comoving width of the shocked material and N_e = 4 π R³ n/3 is the total number of swept-up electrons in the shocked fluid. ${{P}^{\prime}}_{{\nu }^{\prime},max}={\phi }_{p}\frac{\sqrt{3}{q}_{e}^{3}B}{{m}_{e}{c}^{2}}$ (Wijers & Galama 1999) and f(ν') are respectively the comoving peak spectral power and dimensionless spectrum in the comoving frame, where ϕ_p is the dimensionless peak flux and q_e is the electron charge. $B=\sqrt{8\pi n{m}_{p}{c}^{2}{\epsilon }_{B}(\Gamma -1)(4\Gamma +3)}$ (Sari & Piran 1995; Sari et al. 1998) is is the magnetic field strength of the comoving frame, where _B is the energy fraction of the magnetic field. For the fast cooling and slow cooling cases, f(ν') are respectively given by (Sari et al. 1998)

$$\begin{eqnarray}\begin{array}{ll} & f({\nu }^{\prime})=\\ & \left\{\begin{array}{ll}{({\nu }^{\prime}/{{\nu }^{\prime}}_{c})}^{1/3} & ({\nu }^{\prime}\lt {{\nu }^{\prime}}_{c})\\ {({\nu }^{\prime}/{{\nu }^{\prime}}_{c})}^{-1/2} & ({{\nu }^{\prime}}_{c}\lt {\nu }^{\prime}\lt {{\nu }^{\prime}}_{m})\\ {({{\nu }^{\prime}}_{m}/{{\nu }^{\prime}}_{c})}^{-1/2}{({\nu }^{\prime}/{{\nu }^{\prime}}_{m})}^{-p/2} & ({\nu }^{\prime}\gt {{\nu }^{\prime}}_{m})\end{array},\right.\end{array}\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}\begin{array}{ll} & f({\nu }^{\prime})=\\ & \left\{\begin{array}{ll}{({\nu }^{\prime}/{{\nu }^{\prime}}_{m})}^{1/3}\,\,\, & ({\nu }^{\prime}\lt {{\nu }^{\prime}}_{m})\\ {({\nu }^{\prime}/{{\nu }^{\prime}}_{m})}^{-(p-1)/2} & ({{\nu }^{\prime}}_{m}\lt {\nu }^{\prime}\lt {{\nu }^{\prime}}_{c})\\ {({{\nu }^{\prime}}_{c}/{{\nu }^{\prime}}_{m})}^{-(p-1)/2}{({\nu }^{\prime}/{{\nu }^{\prime}}_{c})}^{-p/2} & ({\nu }^{\prime}\gt {{\nu }^{\prime}}_{c})\end{array},\right.\end{array}\end{eqnarray} \tag{ 8 }$$

where ${{\nu }^{\prime}}_{m}=\frac{3{x}_{p}{\gamma }_{m}^{2}{q}_{e}B}{4\pi {m}_{e}c}$ and ${{\nu }^{\prime}}_{c}=0.286\frac{3{\gamma }_{c}^{2}{q}_{e}B}{4\pi {m}_{e}c}$ are respectively the comoving synchrotron typical frequency and the cooling frequency (Wijers & Galama 1999). x_p is the dimensionless peak frequency of the spectrum and x_p ≈ ϕ_p ≈ 0.6 are adopted (Wijers & Galama 1999). p is the electron spectrum index. ${\gamma }_{m}=\frac{{m}_{p}}{{m}_{e}}\frac{p-2}{p-1}{\epsilon }_{e}(\Gamma -1)$ and ${\gamma }_{c}\simeq \frac{6\pi {m}_{e}c}{{\sigma }_{T}{B}^{2}\Gamma T}$ are respectively the minimum LF of the injected electrons and the cooling LF of the electrons, where _e is the energy fraction of the electrons, m_e is the electron mass, and σ_T is Thomson scattering cross section. Note that the synchrotron self-absorption (SSA) is not considered in our calculation, since the rising phase of the observed radio, optical, and X-ray light curves of GRB 170817A are in the same spectral segment. The observed flux thus can be written as

$$\begin{eqnarray}\begin{array}{ll}{F}_{\nu }(T) & =\frac{1+z}{6\pi {D}_{L}^{2}}\displaystyle {\int }_{max({\theta }_{v}-{\theta }_{j},0)}^{{\theta }_{v}+{\theta }_{j}}\Delta \phi ({\theta }_{v},{\theta }_{j},\theta)\\ & \times {R}^{3}(T,\Gamma)n{\delta }^{3}f({\nu }^{\prime}){{P}^{\prime}}_{{\nu }^{\prime},max}\cos \theta \sin \theta d\theta.\end{array}\end{eqnarray} \tag{ 9 }$$

2.4. The Analytical Light Curve Scalings from a Stratified Jet

We attempt to give the analytical light curve from an SJ with the radial profile of energy. The circumburst medium (CBM) number density n is considered as a constant in this paper. We first consider a spherical outflow. The swept up mass of the external shock is $m=\frac{4}{3}\pi {R}^{3}n{m}_{p}$. For the ultrarelativistic case, we have Γ ≫ 1 and β ∼ 1. Thus we get $g(\Gamma)-1={u}^{2}(1+\frac{{\beta }^{2}}{3})\approx \frac{4{u}^{2}}{3}\approx \frac{4{\Gamma }^{2}}{3}$, and combining Equation (3), we obtain

$$\begin{eqnarray}&&\Gamma ={[\frac{9{E}_{0}}{128\pi {T}^{3}n{m}_{p}{c}^{5}}]}^{\frac{1}{(8+k)}}\propto {T}^{\frac{-3}{(8+k)}},\end{eqnarray} \tag{ 10 }$$

where the observed time T ≈ R/2 Γ² c is used. We can find the scalings: $R\propto {T}^{\frac{2+k}{8+k}}$, $B\simeq \Gamma \sqrt{32\pi n{m}_{p}{c}^{2}{\epsilon }_{B}}\propto \Gamma$, γ_m ∝ Γ, γ_c ∝ Γ⁻³ T⁻¹, N_e ∝ R³, and ${P}_{\nu,max}\simeq \Gamma {{P}^{\prime}}_{{\nu }^{\prime},max}\propto \Gamma B\propto {\Gamma }^{2}$, where P_{ν, max} is the observed peak spectral power. Thus we obtain (Sari et al. 1998; Wijers & Galama 1999)

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\nu }_{m}\simeq \Gamma {{\nu }^{\prime}}_{m}\propto \Gamma {\gamma }_{m}^{2}B\propto {\Gamma }^{4}\propto {T}^{\frac{-12}{(8+k)}}\\ {\nu }_{c}\simeq \Gamma {{\nu }^{\prime}}_{c}\propto \Gamma {\gamma }_{c}^{2}B\propto {\Gamma }^{-4}{T}^{-2}\propto {T}^{\frac{-2(2+k)}{(8+k)}}\\ {F}_{\nu,max}\simeq {N}_{e}{P}_{\nu,max}/4\pi {D}_{L}^{2}\propto {N}_{e}{P}_{\nu,max}\propto {R}^{3}{\Gamma }^{2}\propto {T}^{\frac{3k}{(8+k)}}\end{array}\,,\right.\end{eqnarray} \tag{ 11 }$$

where ν_m , ν_c , and F_{ν, max} are respectively the synchrotron typical frequency, the cooling frequency, the peak flux in the observer frame. There are two types of spectra, depending on the order between ν_m and ν_c (Sari et al. 1998). For the fast cooling (ν_m > ν_c ) regime, the observed flux is given by

$$\begin{eqnarray}{F}_{\nu }\equiv {F}_{\nu,max}f(\nu)={F}_{\nu,max}\left\{\begin{array}{ll}{(\nu /{\nu }_{c})}^{1/3}\propto {T}^{\frac{(11k+4)}{3(8+k)}} & (\nu \lt {\nu }_{c})\\ {(\nu /{\nu }_{c})}^{-1/2}\propto {T}^{\frac{2(k-1)}{(8+k)}} & ({\nu }_{c}\lt \nu \lt {\nu }_{m})\\ {({\nu }_{m}/{\nu }_{c})}^{-1/2}{(\nu /{\nu }_{m})}^{-p/2}\propto {T}^{\frac{-(k+6p-4)}{(8+k)}} & (\nu \gt {\nu }_{m})\end{array},\right.\end{eqnarray} \tag{ 12 }$$

while for the slow cooling (ν_c > ν_m ) regime, the observed flux is

$$\begin{eqnarray}{F}_{\nu }={F}_{\nu,max}\left\{\begin{array}{ll}{(\nu /{\nu }_{m})}^{1/3}\propto {T}^{\frac{(4+3k)}{(8+k)}} & (\nu \lt {\nu }_{m})\\ {(\nu /{\nu }_{m})}^{-(p-1)/2}\propto {T}^{\frac{3(2-2p+k)}{(8+k)}} & ({\nu }_{m}\lt \nu \lt {\nu }_{c})\\ {({\nu }_{c}/{\nu }_{m})}^{-(p-1)/2}{(\nu /{\nu }_{c})}^{-p/2}\propto {T}^{\frac{-(k+6p-4)}{(8+k)}} & (\nu \gt {\nu }_{c})\end{array}.\right.\end{eqnarray} \tag{ 13 }$$

where f(ν) is the dimensionless spectrum in the observer frame. The above analytical results are based on the assumption of a spherical outflow, but we consider a collimated jet in this paper. Actually, if the following two conditions are met, the above analytical results are still applicable to the jet case. The first is that the lateral expansion is not important. The second is that the jet edge is not seen (1/Γ < θ_v + θ_j ). If the lateral expansion is not important but the jet edge is seen (1/Γ > θ_v + θ_j ), then the flux in Equations (11), (12) and (13) should be multiplied by a reduction factor of ${\theta }_{j}^{2}/{(1/\Gamma)}^{2}={\Gamma }^{2}{\theta }_{j}^{2}$ (e.g., Panaitescu et al. 1998; Gao et al. 2013). Thus we get ${F}_{\nu,max}\propto {R}^{3}{\Gamma }^{4}\propto {T}^{\frac{(3k-6)}{(8+k)}}$, and the observed flux for the fast and slow cooling regimes can be respectively given by

$$\begin{eqnarray}{F}_{\nu }={F}_{\nu,max}\left\{\begin{array}{ll}{(\nu /{\nu }_{c})}^{1/3}\propto {T}^{\frac{(11k-14)}{3(8+k)}} & (\nu \lt {\nu }_{c})\\ {(\nu /{\nu }_{c})}^{-1/2}\propto {T}^{\frac{2(k-4)}{(8+k)}} & ({\nu }_{c}\lt \nu \lt {\nu }_{m})\\ {({\nu }_{m}/{\nu }_{c})}^{-1/2}{(\nu /{\nu }_{m})}^{-p/2}\propto {T}^{\frac{-(k+6p+2)}{(8+k)}} & (\nu \gt {\nu }_{m})\end{array},\right.\end{eqnarray} \tag{ 14 }$$

and

$$\begin{eqnarray}{F}_{\nu }={F}_{\nu,max}\left\{\begin{array}{ll}{(\nu /{\nu }_{m})}^{1/3}\propto {T}^{\frac{(3k-2)}{(8+k)}} & (\nu \lt {\nu }_{m})\\ {(\nu /{\nu }_{m})}^{-(p-1)/2}\propto {T}^{\frac{3(k-2p)}{(8+k)}} & ({\nu }_{m}\lt \nu \lt {\nu }_{c})\\ {({\nu }_{c}/{\nu }_{m})}^{-(p-1)/2}{(\nu /{\nu }_{c})}^{-p/2}\propto {T}^{\frac{-(k+6p+2)}{(8+k)}} & (\nu \gt {\nu }_{c})\end{array}.\right.\end{eqnarray} \tag{ 15 }$$

Note that the pure edge effect has no effect on the Newtonian phase since Γ ∼ 1 in this phase. Given Γ ∼ 1 and ε ≪ 1, we can get $g(\Gamma)-1={u}^{2}(1+\frac{{\beta }^{2}}{3})\approx {u}^{2}\approx {\beta }^{2}$, and combining Equation (3), we obtain

$$\begin{eqnarray}&&\beta \approx {[\frac{3{E}_{0}}{4\pi {R}^{3}n{m}_{p}{c}^{2}}]}^{\frac{1}{(2+k)}}\propto {R}^{\frac{-3}{(2+k)}}\,.\end{eqnarray} \tag{ 16 }$$

Considering R ≈ ε c T, one can find the scalings: $\beta \propto {T}^{\frac{-3}{(5+k)}}$, $R\propto {T}^{\frac{2+k}{5+k}}$, $B\propto \beta \propto {T}^{\frac{-3}{(5+k)}}$, ${\gamma }_{m}\propto {\beta }^{2}\propto {T}^{\frac{-6}{(5+k)}}$, ${\gamma }_{c}\propto {B}^{-2}{T}^{-1}\propto {T}^{\frac{(1-k)}{(5+k)}}$, ${\nu }_{m}\propto {\gamma }_{m}^{2}B\propto {T}^{\frac{-15}{(5+k)}}$, ${\nu }_{c}\propto {\gamma }_{c}^{2}B\propto {T}^{\frac{-(1+2k)}{(5+k)}}$, and ${F}_{\nu,max}\propto {N}_{e}{P}_{\nu,max}\propto {R}^{3}B\propto {T}^{\frac{3(1+k)}{(5+k)}}$. Thus the observed flux in the Newtonian phase can be given by

$$\begin{eqnarray}{F}_{\nu }=\left\{\begin{array}{ll}{(\nu /{\nu }_{m})}^{-(p-1)/2}{F}_{\nu,max}\propto {T}^{\frac{3(2k-5p+7)}{2(5+k)}} & ({\nu }_{m}\lt \nu \lt {\nu }_{c})\\ {({\nu }_{c}/{\nu }_{m})}^{-(p-1)/2}{(\nu /{\nu }_{c})}^{-p/2}{F}_{\nu,max}\propto {T}^{\frac{(4k-15p+20)}{2(5+k)}} & (\nu \gt {\nu }_{c})\end{array}.\right.\end{eqnarray} \tag{ 17 }$$

These results, if there is no energy injection or energy injection has ceased, are consistent with those derived by some authors (Dai & Lu 1999; Frail et al. 2000; Livio & Waxman 2000; Huang & Cheng 2003; Zhang & MacFadyen 2009).

For more later times, γ_m < 2 is reached, below which the synchrotron approximation becomes invalid. We can neglect the emission from electrons with γ_e < 2. Using the treatment of Granot et al. (2006), the total electron number emitting synchrotron photons N_e ∝ β² R³ and γ_m = 2 (or see the similar treatment in Sironi & Giannios 2013 and Huang & Li 2018b), we can find

$$\begin{eqnarray}{F}_{\nu }=\left\{\begin{array}{ll}{(\nu /{\nu }_{m})}^{-(p-1)/2}{F}_{\nu,max}\propto {T}^{\frac{3(2k-p-1)}{2(5+k)}} & {\nu }_{m}\lt \nu \lt {\nu }_{c}\\ {(\nu /{\nu }_{c})}^{-p/2}{({\nu }_{c}/{\nu }_{m})}^{-(p-1)/2}{F}_{\nu,max}\propto {T}^{\frac{(4k-3p-4)}{2(5+k)}}\,\, & \nu \gt {\nu }_{c}\end{array}.\right.\end{eqnarray} \tag{ 18 }$$

When _e is small, it is possible that γ_m < 2 arrives earlier than the Newtonian phase of the bulk LF (see the fitting results of GRB 170817A in the next section). In this situation, the velocity of the jet is relativistic, while γ_m < 2 is reached. If the jet edge is seen, the peak flux should be corrected as ${F}_{\nu,max}\propto {\Gamma }^{2}(\Gamma -1){N}_{e}{P}_{\nu,max}\approx {\Gamma }^{3}{N}_{e}{P}_{\nu,max}\propto {T}^{\frac{3(k-1)}{(8+k)}}$, where the factor Γ – 1 is taking into account the fraction of the relativistic electrons (γ_e > 2) in the total electrons. Thus the observed flux can be given by

$$\begin{eqnarray}{F}_{\nu }=\left\{\begin{array}{ll}{(\nu /{\nu }_{m})}^{-(p-1)/2}{F}_{\nu,max},\propto {T}^{\frac{3(-2+k-p)}{(8+k)}} & {\nu }_{m}\lt \nu \lt {\nu }_{c}\\ {(\nu /{\nu }_{c})}^{-p/2}{({\nu }_{c}/{\nu }_{m})}^{-(p-1)/2}{F}_{\nu,max}\propto {T}^{\frac{(2k-3p-2)}{(8+k)}}\,\, & \nu \gt {\nu }_{c}\end{array}.\right.\end{eqnarray} \tag{ 19 }$$

Note that all the scalings in this section are applicable to the case of on-axis (θ_v ≲ θ_j ) observation and the case of ${\theta }_{v}-{\theta }_{j}\lesssim \frac{1}{\Gamma }$.

2.5. Estimation of the Source Size and Angular Displacement of Flux Centroid

We estimate the observed angular size and the angular displacement of flux centroid of GRB afterglow image by establishing the following coordinate system (e.g., Gill & Granot 2018). The jet symmetry axis is selected as the z-axis, while the line of sight is the $\tilde{z}$-axis and is in the x – z plane. The y and $\tilde{y}$ axes coincide. The observed image lies in the $\tilde{x}-\tilde{y}$ plane. The observed image has two mutually perpendicular scales, where ${R}_{\perp \tilde{x}}$ is the size perpendicular to the line of sight and is parallel to the $\tilde{x}$-axis, while ${R}_{\perp \tilde{y}}$ is the size perpendicular to the line of sight and is parallel to the $\tilde{y}$-axis. For a given observed time T, ${R}_{\perp \tilde{x}}$ can be estimated as

$$\begin{eqnarray}{R}_{\perp \tilde{x}}\approx \left\{\begin{array}{ll}{R}_{g}\sin (\frac{1}{\Gamma }) & ({\theta }_{j}-{\theta }_{v}\gtrsim \frac{1}{\Gamma })\\ [{R}_{g}\sin (\frac{1}{\Gamma })-{R}_{e}\sin ({\theta }_{v}-{\theta }_{j})]/2 & (0\lt {\theta }_{v}-{\theta }_{j}\lesssim \frac{1}{\Gamma })\end{array},\right.\end{eqnarray} \tag{ 20 }$$

where R_g and R_e are the shock radii on θ = 1/Γ and on the jet edge of θ = θ_v – θ_j at the observed time T, respectively. The radii can be derived by using the dynamic evolution of Γ – R relation and Equation (5). ${R}_{\perp \tilde{y}}$ can be estimated by

$$\begin{eqnarray}{R}_{\perp \tilde{y}}\approx \left\{\begin{array}{ll}{R}_{g}\sin (\frac{1}{\Gamma }) & ({\theta }_{j}-{\theta }_{v}\gtrsim \frac{1}{\Gamma })\\ {R}_{g}\sin (\frac{1}{\Gamma })\sin \Delta {\phi }_{g} & ({\theta }_{v}\gt {\theta }_{j},{\theta }_{v}-{\theta }_{j}\lt \frac{1}{\Gamma }\lt {\theta }_{m})\\ {R}_{m}\sin {\theta }_{m}\sin \Delta {\phi }_{m} & ({\theta }_{v}\gt {\theta }_{j},{\theta }_{m}\le \frac{1}{\Gamma })\end{array},\right.\end{eqnarray} \tag{ 21 }$$

where $\Delta {\phi }_{g}={\cos }^{-1}(\frac{\cos {\theta }_{j}-\cos (\frac{1}{\Gamma })\cos {\theta }_{v}}{\sin {\theta }_{v}\sin (\frac{1}{\Gamma })})$, $\Delta {\phi }_{m}={\cos }^{-1}(\frac{\cos {\theta }_{j}-\cos {\theta }_{m}\cos {\theta }_{v}}{\sin {\theta }_{v}\sin {\theta }_{m}})$, θ_m = cos⁻¹(cosθ_j cosθ_v ), and R_m is the shock radius on θ = θ_m at the observed time T. Hence, the observed angular sizes at the two mutually orthogonal directions are given by

$$\begin{eqnarray}&&\left\{\begin{array}{l}{\theta }_{\perp \tilde{x}}\approx \frac{{R}_{\perp \tilde{x}}}{{D}_{A}}\approx {(1+z)}^{2}{R}_{\perp \tilde{x}}/{D}_{L}\\ {\theta }_{\perp \tilde{y}}\approx \frac{{R}_{\perp \tilde{y}}}{{D}_{A}}\approx {(1+z)}^{2}{R}_{\perp \tilde{y}}/{D}_{L}\end{array},\right.\end{eqnarray} \tag{ 22 }$$

where D_A = D_L /(1 + z)² is the angular distance of the burst. The flux centroid is moving along the $\tilde{x}$-axis, and by taking an approximation, we obtain

$$\begin{eqnarray}{\tilde{x}}_{fc}(T)\approx \left\{\begin{array}{ll}0 & ({\theta }_{j}-{\theta }_{v}\gtrsim \frac{1}{\Gamma })\\ {R}_{v}\sin {\theta }_{v} & ({\theta }_{v}\gt {\theta }_{j},{\theta }_{v}+{\theta }_{j}\lesssim \frac{1}{\Gamma })\\ {R}_{e}\sin ({\theta }_{v}-{\theta }_{j}) & ({\theta }_{v}\gt {\theta }_{j},{\theta }_{v}-{\theta }_{j}\gt \frac{1}{\Gamma })\\ [{R}_{g}\sin (\frac{1}{\Gamma })+{R}_{e}\sin ({\theta }_{v}-{\theta }_{j})]/2 & ({\theta }_{v}\gt {\theta }_{j},{\theta }_{v}-{\theta }_{j}\lt \frac{1}{\Gamma },{\theta }_{v}+{\theta }_{j}\gg \frac{1}{\Gamma })\end{array},\right.\end{eqnarray} \tag{ 23 }$$

where R_v is the shock radius on the jet symmetry axis of θ = θ_v at the observed time T. Thus the angular displacement and the apparent velocity of the flux centroid are respectively given by

$$\begin{eqnarray}&&\left\{\begin{array}{l}\Delta {\theta }_{fc}\approx \frac{|{\tilde{x}}_{fc}({T}_{2})-{\tilde{x}}_{fc}({T}_{1})|}{{D}_{A}}\approx {(1+z)}^{2}|{\tilde{x}}_{fc}({T}_{2})-{\tilde{x}}_{fc}({T}_{1})|/{D}_{L}\\ {\beta }_{app}\approx \frac{|{\tilde{x}}_{fc}({T}_{2})-{\tilde{x}}_{fc}({T}_{1})|}{c({T}_{2}-{T}_{1})}\end{array}.\right.\end{eqnarray} \tag{ 24 }$$

We use a package MceasyFit in the light curve fitting. The package is based on the Markov-Chain Monte Carlo (MCMC) method, which is described in Zhang et al. (2016). We fit 3 GHz and 5.5 GHz radio, optical 5.1 × 10¹⁴ Hz (F606W), and X-ray (1keV) afterglow light curves. The afterglow data in the four bands are taken from the following papers: Hallinan et al. (2017); Alexander et al. (2018); Dobie et al. (2018); Margutti et al. (2018); Mooley et al. (2018a,b); Piro et al. (2019); Lyman et al. (2018); Lamb et al. (2019); Nynka et al. (2018); Troja et al. (2018,2019,2020), and Hajela et al. (2019). There are nine free parameters in our fitting, i.e.,

$$\begin{eqnarray}&&\Phi =[{\theta }_{0},{\theta }_{v},{\epsilon }_{B},{\epsilon }_{e},{E}_{1},n,{\Gamma }_{{\rm{\min }}},k,p],\end{eqnarray} \tag{ 25 }$$

where θ₀ is the initial half opening angle of the jet, and E₁ is a constant, which is related to E₀ by E₀ = E₁ × (1 – cos θ₀). In the fitting, we confined the free parameters as θ₀ ∈ [0.01,0.45], θ_v ∈ [0,0.6], log₁₀(_B ) ∈ [–7,–0.1], log₁₀(E₁) ∈ [53,66], log₁₀(n) ∈ [–7.5,1], Γ_min ∈ [2,10], k ∈ [2,10], p ∈ [2.05,2.6], and log₁₀(_e ) ∈ [–6,–0.1]. Note that the viewing angle is constrained to θ_v > θ₀ in the fitting since the superluminal motion of the centre-of-brightness on the sky in GRB 170817A was observed (Mooley et al. 2018c).

We use the SJ models to fit the afterglow datasets of GRB 170817A, and the best fit parameters for the spreading and the non-spreading cases are displayed in Table 1. The 3 GHz and 5.5 GHz radio, optical (F606W), and X-ray (1keV) light curves derived with the best fit parameters are shown in Figure 1 and Figure 2 for the spreading jet and non-spreading jet, respectively. The corner plots for the spreading jet and non-spreading jet are respectively shown in Figure 3 and Figure 4. The dynamic evolutions of the spreading and non-spreading jets are shown in Figure 5. As shown in Figure 1 and Figure 2, both the models can give good fits to the afterglow datasets of GRB 170817A with χ²/DOF ≈ 1.47 and ≈ 1.38 for the non-spreading and spreading jets, respectively. As shown in Table 1, θ₀ ≈ 8.3° and θ_v ≈ 9.9° for non-spreading jet, and θ₀ ≈ 5.7° and θ_v ≈ 8.2° for spreading jet. These results suggest that both of these two scenarios are slightly off-axis observed at the beginning of the afterglow. The energy injection indexes take the values of k ≃ 7.1 and k ≃ 8.4 for non-spreading and spreading jets, respectively. The true energy of the non-spreading and spreading jets are respectively given by E_jet ≃ E_max = E₀ (Γ_min β_min)^−k = E₁ × (1 – cosθ₀) × (Γ_min β_min)^−k ≃ 4.8 × 10⁵¹ erg and ≃ 4.6 × 10⁵² erg.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 5 shows the evolution of the LF of the fastest ejecta (upper panel), evolution of the half opening angle of the jet (middle panel), and evolution of the synchrotron typical frequency ν_m and synchrotron cooling frequency ν_c (bottom panel). As shown in Figure 5, both of the spreading and non-spreading jets' LF evolution curves show a change in the decline slope at about 160 days and it indicates the end of the energy injection. At the same time, the observed flux reaches its peak and begins to decrease. The evolution of ν_m and ν_c indicates the radio, optical, and X-ray emission are in the same spectral regime, i.e., slow cooling in the ν_m < ν < ν_c regime. Thus the light curves of the four bands in our fitting have the same slope. In the rising phase of the afterglow, the slopes of the LF evolution curves for the non-spreading and spreading jets respectively are $f\simeq -\frac{3}{8+k}\simeq -0.20$ and ≃ –0.18, while the analytical temporal indexes of the light curves are given by $\alpha \simeq \frac{3(2-2p+k)}{8+k}\simeq 0.9$ (see Eq. (13)) and ≃ 1.1 in our fittings. The analytical temporal index α ≃ 1.1 for the spreading jet is a little larger than α ≃ 0.8 given by Mooley et al. (2018a). This is due to the lateral expanding effect in the spreading jet. The jet energy will partly spread laterally and be decelerated faster (see Fig. 5, even in the rising phase with small lateral expanding effect). Thus the rising slope of the light curve would be more shallow than the analytical one and the true temporal index of the rising phase can be decreased as α ≃ 0.8. For the non-spreading jet, with the decelerates of the jet, the visible area $\theta \approx \frac{1}{\Gamma }$ has exceeded the angle of θ_v – θ_j after ∼ 35 days. For the spreading jet, the opening angle θ_j has already exceeded the viewing angle θ_v at ∼ 1 day, thus the scenario has become on-axis from slightly off-axis at early rising phase. Hence, the scalings in Section 2.4 are roughly applicable for both of the non-spreading and spreading jet models in our fittings.

Table 1. Best Fit Parameters for the SJ Model for GRB 170817A

Parameter	Spreading jet	Non-spreading jet
k	${8.42}_{-0.58}^{+0.58}$	${7.07}_{-0.44}^{+0.78}$
p	${2.18}_{-0.01}^{+0.00}$	${2.20}_{-0.01}^{+0.01}$
Γmin	${5.42}_{-1.99}^{+0.18}$	${3.18}_{-1.11}^{+0.91}$
${\theta }_{0}{(}^{^\circ })$	${5.71}_{-4.18}^{+2.72}$	${8.30}_{-1.49}^{+3.53}$
${\theta }_{v}{(}^{^\circ })$	${8.18}_{-3.70}^{+3.24}$	${9.92}_{-3.35}^{+3.27}$
log10(B )	$-{5.62}_{-3.52}^{+0.30}$	$-{5.79}_{-2.55}^{+0.48}$
log10(e )	$-{2.79}_{-0.09}^{+0.09}$	$-{2.48}_{-0.72}^{+0.14}$
log10(E1(erg))	${61.08}_{-1.11}^{+1.11}$	${57.05}_{-1.34}^{+1.34}$
log10(n(cm−3))	$-{3.65}_{-0.92}^{+0.92}$	$-{2.01}_{-0.55}^{+2.50}$
BIC	104.40	108.58
χ2	67.86	72.04
χ2/DOF	1.38	1.47
βapp	0	3.6
${\theta }_{\perp \tilde{x}}$ (mas)	6.3	1.5
${\theta }_{\perp \tilde{y}}$ (mas)	6.3	2.2

E₁ is not the true energy of the jet, and the true energy of the jet is given by E_jet ≃ E_max = E₀ (Γ_min β_min)^−k = E₁ × (1–cosθ₀) × (Γ_min β_min)^−k . Thus the energy of the non-spreading and spreading jets respectively are E_jet ≃ 4.8 × 10⁵¹erg and E_jet ≃ 4.6 × 10⁵²erg. Note that BIC = χ² + n_p × ln(n_d ), where n_p = 9 is the number of free parameters, n_d = 58 is the number of data points, and DOF = n_d – n_p = 49 is the degree of freedom.

In the decay phase, the light curves for both spreading and non-spreading jets first experience a rapid decline from ∼ 160 days, and then gradually transit to a shallow decline after ∼ 600 days (see Fig. 1) and ∼ 400 days (see Fig. 2). By numerical calculation, we found that γ_m < 2 occurs at around the peak times (∼ 160 days) for both spreading and non-spreading jets, earlier than the arrival times of the Newtonian phase Γ ∼ 1 due to the small values of _e in the two cases. The origin of the rapid decline for the spreading jet is that the lateral expansion is significant (see Fig. 5) after the light curve peaks. Although the jet edge is not seen at that time, the change of dynamic evolution due to the lateral expansion will lead to an asymptotic decay with the slope of α ∼ –p ∼ –2.2 in the light curves (Sari et al. 1999). For the non-spreading jet, the jet edge is seen (1/Γ > θ_v + θ_j ) and the corresponding analytical temporal index is $\alpha {\sim}\frac{3(k-p-2)}{(8+k)}{\sim}-\frac{-3(p+2)}{8}{\sim}-1.6$ (see Eq. (19)). Note that Equation (19) is obtained by assuming an ultrarelativistic (Γ ≫ 1) jet. However, the LF of the non-spreading jet is Γ ∼ 3 at the peak time, at which the ultrarelativistic approximation Γ ≫ 1 is not very appropriate, so Equation (19) is only marginally consistent with the numerical results. The reasons for the shallow decline after ∼ 600 days and ∼ 400 days for spreading and non-spreading jets are that the jets begin to transit to the Newtonian phase (Γ < 2) with $\alpha {\sim}\frac{3(2k-p-1)}{2(5+k)}{\sim}\frac{-3(p+1)}{10}{\sim}-1.0$ (see Eq. (18)).

We have shown that both of the non-spreading and spreading jets can provide good fits to the afterglow light curves of GRB 170817A with plausible parameters, but the observation limits of the source angular size and the apparent velocity of the flux centroid can further constrain the two jet models. The apparent velocity of the flux centroid has been obtained as β_app = 4.1 ± 0.5 between 75 and 230 days from VLBI by Mooley et al. (2018c), and the observed angular size of the radio image at 207 days is constrained as θ_⊥ < 2.5 mas at 90% confidence level by Ghirlanda et al. (2019). By taking the estimation methods in Section 2.5, we get β_app ≈ 3.6 for the non-spreading jet, while β_app ≈ 0 due to ${\theta }_{j}-{\theta }_{v}\gtrsim \frac{1}{\Gamma }$ (see Eq. (23)) at 75 and 230 days for the spreading jet. Moreover, the observed angular sizes of the source are ${\theta }_{\perp \tilde{x}}\approx 1.5$ mas and ${\theta }_{\perp \tilde{y}}\approx 2.2$ mas for the non-spreading jet, while ${\theta }_{\perp \tilde{x}}\approx {\theta }_{\perp \tilde{y}}\approx 6.3$ mas for the spreading jet. Hence, for the non-spreading jet model, both the observed angular size of the source and the apparent velocity of flux centroid meets the observation limits, while those for the spreading jet model are beyond the observation limits.

In this paper, we consider an SJ model to explain the peculiar afterglow light curves of GRB 170817A and its observed image data. We investigate two kinds of jet evolutions, including a spreading jet and a non-spreading jet. By fitting the afterglow light curves of the burst in multiple bands, we find both the spreading and non-spreading jet models can give good explanation to the light curves. However, by comparing the theoretical image properties with the observations, we find the apparent angular size of the jet and the apparent velocity of the flux centroid for the non-spreading jet model can satisfy the observation limits, while those for the spreading jet model violate the observation limits. This suggests that this burst may arise from an SJ and that the spreading of the jet is not significant, at least at ∼ 230 days post-burst. In addition, the CBM number density for the non-spreading jet model is n ≃ 10^−2.01cm⁻³, consistent with the fact that the short burst happened in the outskirt of its host galaxy where the density should fall in between the intergalactic medium density (∼ 10⁻⁶ cm⁻³) and the ISM number density (∼ 1 cm⁻³).

In the SJ model, the rising phase of the afterglow is due to the energy injection of the stratified ejecta, whose energy profile is given as E ≃ 1.2 × 10⁵⁵ × (Γ β)^−7.1 for the non-spreading jet. The total kinetic energy is E ≃ 4.8 × 10⁵¹ erg and the energy injection index is k ≃ 7.1 for the non-spreading jet. The energy injection index is a little larger than the QSSEI model (see e.g., Gill & Granot 2018; Mooley et al. 2018a; Huang & Li 2018a). Based on our analytical results in Section 2.4 (see Eq. (5), also see Gill & Granot 2018), the energy injection index is given by k ≃ (8α + 6p – 6)/(3 – α) in the slow cooling phase in ν_m < ν < ν_c spectral regime (see Eq. (5) and Fig. 5). The temporal index of the rising phase of GRB 170817A is α ≃ 0.8 (Mooley et al. 2018a). For α ≃ 0.8 and p ≃ 2.2, the energy injection index is thus k ≃ 6.2 (Gill &; Granot 2018). For the non-spreading jet, with the decelerates of the jet, the visible area $\theta \approx \frac{1}{\Gamma }$ has exceeded the angle of θ_v – θ_j after ∼ 35 days. Hence, the scalings in Section 2.4 are applicable for the non-spreading jet after ∼ 35 days. Thus the energy injection index k ≃ 7.1 for the non-spreading jet is close to k ≃ 6.2 in the QSSEI model, and the small difference comes from the jet geometry and off-axis observed scenario in the non-spreading jet model.

The steep decay of the light curves after the peak in GRB 170817A can arise from two effects. The first is the significant lateral expansion effect (e.g., Rhoads 1999; Sari et al. 1999) leads to a rapid deceleration of the jet, and thus gives rise to a steep light curve. The second is the jet edge is seen ($\frac{1}{\Gamma }\gt {\theta }_{v}+{\theta }_{j}$) and the observer will feel the deficit of flux outside the jet cone compared with the spherical shock case, which leads to a steeper decay (e.g., Panaitescu et al. 1998). For the models in this paper, the origin of the steep decaying of the light curves in GRB 170817A after the peak for the non-spreading jet is that the jet edge is seen at ∼ 160 days, while for the spreading jet it is that the lateral expansion is significant. The reasons for the shallow decay after ∼ 400 days and ∼ 600 days for non-spreading and spreading jets are that the jets begin to transit to the Newtonian phase with Γ < 2. The shallow decay is consistent with the observations. It is possible that the kilonova afterglow or the energy injection by a pulsar (Troja et al. 2020) will contribute to the late afterglow, which will lead to a slight rise at a later time. The counter-jet can also contributes at several thousand days after the GRB trigger and will somewhat change the light curves (e.g., Granot et al. 2018; Li et al. 2019). Further observations in the future might reveal the origin of the shallow decay.

We thank the anonymous referee for useful suggestions and comments. This work was supported by the National Natural Science Foundation of China (Nos. U1831135 and 11833003), Yunnan Natural Science Foundation (2014FB188). B.B.Z acknowledges support by the National Key Research and Development Programs of China (2018YFA0404204), the National Natural Science Foundation of China (grant Nos. 11833003 and U2038105), and the Program for Innovative Talents, Entrepreneur in Jiangsu. We also acknowledge the use of public data from the Fermi Science Support Center (FSSC).