Afterglow Light Curves from Off-Axis GRB Jets in Stratiﬁed Circumburst Medium

: We study the gamma-ray burst (GRB) afterglow light curves produced by an off-axis jet in a stratiﬁed circumburst medium and summarize the temporal indices of the coasting phase, the deceleration phase, the Newtonian phase, and the deep Newtonian phase for various viewing angles and power-law indices of medium density. Generally, the afterglow light curves of off-axis GRBs in the homogeneous interstellar medium have a steep rise arising due to jet deceleration. In the stratiﬁed medium, the ﬂux rises is more shallow but peaks earlier for the same viewing angle due to faster deceleration of the jet running into the denser stratiﬁed medium, compared with the case of the interstellar medium (ISM). Observations of off-axis bursts will possibly increase over the coming years due to the arrival of the multi-messenger era and the forthcoming surveys in multiple bands. The temporal indices of off-axis afterglows derived in the paper will provide a reference for comparison with the observations and can diagnose the circumburst environment. The numerical code calculating the afterglow light curve also can be used to ﬁt the multi-wavelength light curves.


Introduction
In the gamma-ray burst (GRB) standard model, afterglows are produced by the interaction between a fireball shell and the homogeneous interstellar medium (ISM) [1,2]. A pair of shocks, forward shock and reverse shock, will be developed when the shell begins to decelerate considerably after sweeping up enough circumburst medium (CBM). The forward shock starts to experience a self-similar expansion subsequently [3]. Electrons in the shocked medium are accelerated. Assuming a power-law distribution of the electrons and considering the dynamic evolution of the shock, one can find the muti-segment broken power-law radiation spectrum and light curve [4], which is called the standard afterglow model. The power-law decay of light curves predicted by the standard afterglow model has been widely tested and confirmed by abundant observations. However, the afterglow light curves in some bursts displayed a marked achromatic steepening break lately, suggesting that at least some GRB fireballs are possibly not spherical but collimated jets [5,6]. The extremely high isotropic gamma-ray energy in some bursts also can be explained in the jet model. During the ultra-relativistic phase, the light curves of the afterglow from a jet are not different from those from a spherical fireball due to the relativistic beaming effect if the line of sight (LOS) is within the jet. Thus the standard afterglow model still applies in this phase until the jet edge is seen or the lateral expansion becomes significant. After that, the light curve will become steeper [7,8].
A natural consequence of the jet model is the presence of off-axis bursts. Some GRBs have spectral peaks of tens of keV or less, much lower than those of typical bursts which were believed to arise from off-axis jets with viewing angles larger than jet half opening angles e.g., [9,10] or from off-axis structured jets (e.g., [11]). It is also possible that if the observer is off-axis, we could miss GRBs but observe their late afterglows when the beaming effect becomes weak due to jet deceleration. This is the so-called "orphan afterglow", whose detectability has been extensively investigated [7,[12][13][14][15][16][17][18][19][20]. The detection of an orphan afterglow will be essential evidence of the GRB jet model. However, until now, no definite orphan afterglow detection has been reported (however, see [21] for some candidates). Off-axis afterglow light curves have been studied for various jet models [13,20,[22][23][24][25]. Generally, the light curves have a rapid rise initially and then peak when the jet Lorentz factor Γ decreases down to γ ∼ 1/θ 0 , where θ 0 is the jet half opening angle. The rising amplitude depends on the viewing angle if the CBM is homogeneous.
Another factor that can considerably affect the afterglow light curve is the circumburst environment. Observations have supported the model that long GRBs are produced in the collapse of a massive star [26][27][28]. Thus, the CBM density is expected to be inhomogeneous and decrease in radius due to the wind from the progenitor star. For a free wind with constant mass-loss rate and wind speed, the number density of CBM scales as n ∝ R −2 with the radius R of the afterglow shock [29][30][31][32]. The power-law slope of density can deviate from two due to variable mass-loss rate and wind speed, which are poorly known. Short bursts usually happen outside of their host galaxy due to a long migration with a large kick speed before the merger of a compact binary, where the CBM density transits from the ISM to the intergalactic medium (IGM). Their CBM is usually homogeneous and is more rarefied than that of long bursts [33]. However, it is not excluded that some short bursts happen in a density fluctuation region, which may have an ununiform density profile, approximately described by a power-law profile (stratified medium). Indeed, there are individual bursts, whose light curves can be explained by the wind medium model e.g., [34]. Some temporal properties of the afterglow's forward and reverse shocks in the stratified medium have been investigated [35,36]. Compared with the ISM case, in the stratified medium, the density is large in small radii. Thus the jet deceleration is fast at the beginning and will gradually become slower due to the decreasing medium density with the increase in the radius. This will lead to different light curves from those in the standard model.
Recently, a well-known short GRB 170817A was found, which was the first GRB associated with a gravitational-wave event [37] and marked the dawn of multi-messenger astronomy. Compared with other GRBs, the afterglow of the burst has a multi-wavelength rise with a power-law slope of ∼0.8, lasting for as long as ∼160 days, which has never been found in other bursts. The isotropic energy of the burst is 3-4 orders of magnitude lower than the conventional short bursts. These unconventional properties strongly suggest that this burst arises from a structured jet seen off-axis [38][39][40][41][42][43][44][45]. Thus this detection of the burst opened the era of the observation of off-axis bursts.
Stimulated by the long rise of GRB 170817A, here we systematically study the off-axis afterglow light curves from a uniform GRB jet in the stratified medium (although this burst may not be the case), paying more attention to the analytical scalings and numerical calculation of the light curves. The scalings and light curves could serve as a reference to the observations of off-axis bursts in the future. The paper is organized as follows. In Section 2, we analytically derive the off-axis light curves in the stratified medium in various regimes. In Section 3, we numerically calculate the light curves using a dynamic model to test the analytical results. The summary and discussion are provided in the last section. Table 1 lists the acronyms used in the paper. Table 1. Alphabetically ordered list of the acronyms used in this paper.

CBM
circumburst medium GRB gamma-ray burst IGM intergalactic medium

The Analytical Light Curves from an Off-Axis Jet in the Stratified Medium
Consider a relativistic outflow at a radius R from the burst source collimated into an initial opening angle θ 0 running into the CBM. The CBM can be the ISM or the wind-like medium whose density can be described by a simple power-law profile n ∝ R −2 e.g., [30,31]. More generally, regardless of the exact mechanism, one can consider a power-law stratified number density as where n 0 is the density at an initial radius of R 0 . A forward shock and a reverse shock will form. Here we only consider the forward shock. When the shock sweeps up enough CBM whose mass is equal to γ −1 0 (γ 0 is the initial Lorentz factor of the ejecta) of the rest mass of the ejecta, it begins to decelerate significantly and enters a self-similar expansion, corresponding to the onset of the GRB afterglow [46,47]. Before the deceleration, the Lorentz factor of the shock is nearly constant, corresponding to the so-called coasting phase.

Coasting Phase
With the evolution of the shock, its radius will increase with the burst source time t = T + R cos θ/c. Here T is the observer time, θ is the angle between the speed of a fluid element within the jet and the viewing line, and c is the light speed. Note that in the analytical consideration, we neglect the redshift z. In the coasting phase, the shock speed β = γ 2 − 1 is a constant. The observed time is thus T = R(1 − β cos θ)/c. For the off-axis case, the emission of the jet mainly comes from the nearest angular region of the jet from the viewing line unless the viewing line is very close to the jet edge. Thus, as an approximation, we only consider the nearest point of the jet from the viewing line. Given the viewing angle θ v , for θ v − θ 0 1 and γ 1, one can find In the following analytical consideration, we focus our attention on the off-axis case where γ 2 (θ v − θ 0 ) 2 1 is easily satisfied before significant lateral expansion (we neglect the lateral expansion in the rising phase of the light curves in the analytical estimation) unless the viewing line is very close to the edge of the jet.
We can find R ∝ T and n ∝ R −s ∝ T −s . The total electron number in the shocked ejecta is N e ∝ R 3−s . The observed synchrotron typical frequency and the cooling frequency are with m e being the electron mass and e being the electron energy fraction. δ D ≈ 2γ/[1 + γ 2 (θ v − θ 0 ) 2 ] ∝ γ −1 is defined as the "effective" Doppler factor, which is the Doppler effect of the fluid element at θ v − θ 0 . The magnetic field behind the shock is B = [8π B nm p c 2 4γ(γ − 1)] 1/2 ≈ (32πγ 2 m p c 2 B n) 1/2 ∝ γR −s/2 [48], B is the magnetic field energy fraction of the internal energy and m p is the proton mass. Note that primed quantities are measured in the shock co-moving frame and that following the often-used notation, the magnetic field, energy equipartition parameters, electron energies, and electron distribution are unprimed, although they are measured in the co-moving frame. The observed flux peak is F ν,max ∝ N e δ 3 D P ν,max , where P ν,max ∝ B is the synchrotron peak spectral power. In the coasting phase, we can find that In low-frequency bands (e.g., the radio band), synchrotron self-absorption (SSA) has to be considered. Two regimes, including fast cooling (ν m > ν c ) and slow cooling (ν c > ν m ), will be involved successively due to the evolution of ν c and ν m . For typical parameters, the afterglow light curves are generally in the two regimes of ν a < ν c < ν m and ν a < ν m < ν c . The SSA frequency scales as [49,50] Hence, following the synchrotron spectra [4], the light curves can be described by and If the observed band in early time falls in the spectral segment of ν m < ν < ν c or ν > ν c > ν m , for the ISM case, we have F ν ∝ T 3 or F ν ∝ T 2 , respectively, while for the wind case (s = 2), the light curves scale as F ν ∝ T 0.5 or F ν ∝ T −0.7 (p = 2.3) for ν c < ν < ν m or ν > ν m > ν c , respectively. These are consistent with previous results (e.g., [51,52]). For a general density profile 0 < s < 2, the light curve slopes fall in between the two cases. See Table 2 for the slopes with s = 0, 0.5, 1, 1.5, 2 and p = 2.3. Table 2. Temporal indexes from an off-axis jet in the coasting phase for p = 2.3.

Deceleration Phase
In the decelerating phase, considering the adiabatic case, the isotropic equivalent The Lorentz factor of the shocked fluid scales as Thus, we obtain The SSA frequency in the regimes of ν a < ν c < ν m and ν a < ν m < ν c scales as [49,50] Thus the light curves in the deceleration phase are and The analytical temporal indexes for s = 0, 0.5, 1, 1.5, 2 with a typical value of p ∼ 2.3 for the fast and slow cooling cases in the deceleration phase are shown in Tables 3 and 4, respectively. Note that in Tables 3 and 4, we include the regimes of ν c < ν a < ν m and ν m < ν a < ν c . In the dense CBM, it is also possible that the SSA frequency is larger than max(ν m , ν c ). We derive the light curve indexes for these cases in the Appendix A. We can find if the observed frequency is in ν a < ν m < ν < ν c and s = 0, the light curve index is ∼4.1. If ν > ν c , the rising index is ∼0.7 for s = 2. These are roughly consistent with [52]. Given the rising index of ∼0.8 for light curves of GRB 170817A, the density index would be s ∼ 1.5 if the burst happened in the stratified medium. Table 3. Temporal indexes from an off-axis jet for the fast-cooling regime in the deceleration phase

Late Afterglow
This preceding derivation of light curves does not include the sideways expansion, which is not significant at early times. However, it will become significant when the jet slows down to γ ∼ 1/θ 0 [7,8]. The jet will rapidly start an approximately spherical expansion, and R ∼ const. The light curves for various viewing angles will have a common decline as F ν ∝ T −p [8] until the expansion enters the Newtonian phase when the rest of the swept-up mass-energy is equal to the shock energy. For the ISM and the wind-like medium, the light curves in this phase have been investigated [53][54][55][56][57]. For more general CBM, using the so-called Sedov-Taylor solution e.g., [58], one can find the scalings: . Thus the light curves scale as Considering the density profiles of s = 0 [2], one can find the flux scales as F ν ∝ T 3(7−5p)/10 [(7p−5)/6] and F ν ∝ T (4−3p)/2 [(8−7p)/6] for ν m < ν < ν c and ν > ν c , respectively, which are consistent with previous analytical results (e.g., [54,56]). The scalings within square brackets are for the wind medium. For later times, the shock will enter the deep Newtonian phase with γ m < 2, below which the synchrotron approximation becomes invalid. One can only consider the electrons with energy γ e > 2 and neglect the lower-energy electrons. Using the total electron number emitting synchrotron photons N e ∝ β 2 R 3−s and γ m = 2 ( [59], or see the similar treatment in [60,61]), we can find The indexes for the sideways expanding phase (when the sideways expansion becomes significant), the Newtonian phase, and the deep Newtonian phase are displayed in Table 5. For s = 0 [2] and ν m < ν < ν c , F ν ∝ T −3(p+1)/10[−(p+1)/2] , which is in agreement with previous results (e.g., [60]

Dynamic Evolution of the Afterglow Jet
The analytical estimates give the light curve slopes but neglect some details, including the geometry shape of the jet, sideways expansion of the jet, etc. To derive more precise light curves, we should appeal to numerical calculations of the dynamic evolution and radiation. The dynamic evolution of a jet has been extensively investigated e.g., [7,8,55,57,62,63]. In these studies, the sideways expansion is a crucial factor but remains an open question. Recent numerical and analytical studies have found that sideways expansion is not as fast as predicted in the analytical model (e.g., [7,57]). However, Ref. [64] reached the opposite conclusion. It appears that if the initial opening angle is quite small, sideway expansion approaches the prediction of the analytical model, while if it is large, the sideway expansion is slow [63]. Regardless of the debate, we use the recent dynamic evolution model of [44], in which the evolution equation of the jet opening angle is where u = γβ is the dimensionless velocities, and β ⊥ is the dimensionless lateral expansion speed in the burst source frame. In the dynamic equations in [44], the finding that the "conical" spreading is more accurate than the "trumpet" model [62] is also taken into account (see Ryan et al. [44] for more details on the dynamic equations). The dynamic equations appear to obtain similar light curves to those of the multi-dimensional hydrodynamic simulation.

Flux Calculation
We will numerically calculate the afterglow light curve in this section. Consider a jet with the half opening angle of θ j . The observed flux for a given observed time T can be given by [20,65] where ∆R = R 12γ 2 1 1−β cos θ is the shocked shell width within which the radiation is effectively produced, and D L is the luminosity distance. The observed time T is the function of γ, R, and θ and can be given by For a homogeneous jet, the flux can be reduced as where δ is the Doppler factor for a jet element. j ν = j ν(1+z)/δ is the co-moving emissivity, i.e., where γ max = (6πq e /σ T B) 1/2 is the maximum Lorentz factor of the injected electrons and q e is the electron charge.P ν is the pitch angle (α) averaged synchrotron spectral power, written as [66] y = ν 3γ 2 e ν L /2 and ν L is the Larmor frequency. K 5/3 (t) is the modified Bessel function [67]. n γ e is the electron distribution (the electron number per unit energy Lorentz factor per unit volume). For the fast cooling case (γ m > γ c ), the distribution is where n e = 4γn is the co-moving electron density. While for the slow cooling (γ m < γ c ) case, it is We have taken SSA into account in the calculation. τ = k ν ∆R = k ν ∆R/δ is the SSA optical depth and k ν is the SSA coefficient and it is given by [67]

Code Verification
We have developed a Fortran code to solve the dynamic equations derived by [44] and calculate the light curve. The fourth-order Runge-Kutta algorithm is used to solve R(t), u(t), and θ j (t) on a fixed logarithmically spaced grid of t [68]. For a given observed time T, the observed flux is calculated by numerically integrating over the equal arrival time surface defined by Equation (19). The Steffensen iteration method is used to find the integral interval faster than the usual iteration method. The Gaussian integration is used [68]. The integral interval is split into multiple segments to improve the integral precision. To test the code, we compare the calculated light curves with those calculated by afterglowpy (see Ryan et al. [44]), shown in Figure 1. On the whole, the agreement of the two sets of light curves is good, especially in the early and late afterglow phases. In the intermediate phase, our results are somewhat higher than that derived by afterglowpy. In detail, our results are somewhat higher than those derived by afterglowpy, especially in the optical and radio bands (top panel of Figure 1). This is due to the fact that we consider the pitch angle as an averaged synchrotron spectrum, which will lead to a slight change in the synchrotron peak frequency and spectral power. If we use the analytical synchrotron spectrum as in [44], the agreement of the two codes is fairly good (see the bottom panel in Figure 1). There is also a slight difference in the intermediate phase, which can be attributed to the different numerical methods or numerical precision. It appears that our results approach those of Boxfit (see Figure 2 in Ryan et al. [44]).  In the top and bottom panels, the light curves are derived with the numerical and analytical synchrotron spectra, respectively. The latter spectrum is also adopted in Ryan et al. [44].

Comparison with the Analytical Results
In order to compare with the analytical results, we calculate the following light curve adopting the following typical parameters: 3, E iso = 10 53 erg, z = 0.1, and d L = 1.4 × 10 27 cm. The initial Lorentz factor γ 0 = 300 is used. The CBM number density for various s is set to be the typical value of the ISM at the on-axis jet break radius R j , where the lateral expansion becomes significant, i.e., n(R j ) = n 0 (R j /R 0 ) −s = n i R −s j = 1 [22,69]. The jet break radius can be given from For s = 0, we can find at an early time, that the X-ray bands (ν X = 1 keV) are in the regime of ν X > ν m > ν c or ν X > ν c > ν m in the coasting phase and thus the slope is ∼2 (see Table 2). Then the X-ray gradually goes into the regime of ν X > ν c > ν m of the deceleration phase, so the X-ray light curve transits from 2 to 4.6 (see Table 3). The optical (ν opt = 6 × 10 14 Hz) flux rises with the slope of 3.0 in the regime of ν c > ν m > ν opt of the coasting phase and then transits from ∼ 3.0 to 2.0 with ν opt > max(ν m , ν c ) and then to 4.6, like the X-ray. The radio (ν R = 1 × 10 9 ) band is in the same regime as the optical band. Thus the radio flux rises with a slope of 3.0 and then 3.7 with ν R < ν c < ν m . Later, the jet edge is visible, whose light curve index is not easy to derive. After the whole jet is visible, the jet roughly declines with a slope of ∼p = 2.3. At a later time, the jet goes into the Newtonian phase. The flux in the three bands all decline with a slope of ∼−1.4 (see Table 5). We can find in the early and late times the analytical slopes are roughly consistent in the numerical results. At the intermediate time, the beaming effect becomes weaker with the deceleration of the jet, and the visible region increases. The jet's sideways expansion also becomes significant. The total effect on the light curves is not easy to describe with the analytical method. The case of s = 2.0 is similar (see Figure 3).

Typical Off-Axis Light Curves in a Stratified Medium
We numerically calculated multi-band light curves, including radio (1 GHz), optical (6.0 × 10 14 Hz), and X-ray (1 keV) bands with different density slopes s = 0, 1, 1.5, 2 and various viewing angles of θ v = 0, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, as shown in Figures 4-6. We find that the light curves exhibit a variety of shapes depending on different viewing angles and medium-density profiles. The resulting light curves are similar to those of [22] except in the radio band. We find the SSA is quite significant for the medium profile with large s, which is not taken into consideration in [22]. We also note that, compared with the ISM case, the off-axis light curves in the stratified medium depict a more shallow rise whose rising slopes depend on the medium density slope, but the light curves peak earlier than in the ISM case. This is because the density of the stratified medium is denser than the ISM, especially at small radii, leading to the faster deceleration of the jet and the jet edge being visible earlier. The rising duration depends on the viewing angle. The large viewing angle arises as a long rising phase but with a lower peak flux, which is harder to detect.

Summary and Discussion
In this paper, we derive the analytical off-axis light curves of GRB afterglows in the stratified external medium with a power-law profile and summarize the scalings of offaxis light curves from the afterglow onset to the deep Newtonian phase (see Tables 2-5), except for the intermediate phase. These light curve scalings can serve as a helpful tool for afterglow observers to quickly identify the spectral regime and the medium density profile from their data. Then we present the numerical afterglow light curves in radio, optical, and X-ray bands for a wide range of viewing angles in the stratified external medium. These results verify the analytical scalings for early and late afterglow. The numerical calculation code of the light curves can also be used to fit the multi-wavelength light curves fitting.
Previous works have derived the light curve scalings of off-axis afterglows for ISM (s = 0) and steady wind (s = 2) (e.g., [51][52][53][54]56,60]). However, some fittings to afterglow data have found that the medium slopes are not 0 or 2 but have other values (e.g., Yi et al. [35], Liang et al. [70], Kouveliotou et al. [71]), suggesting the CBM is not ISM and the preburst wind speed and/or mass-loss rate are not constant, otherwise leading to a medium slope of s = 2. In this paper, we generalize the scalings for arbitrary power-law circumburst density profiles, including ISM and wind. This can qualitatively diagnose the circumburst density profile from afterglow observations. Further, the observed afterglow light curves can be fitted with theoretical light curves by numerical methods, which will be used to determine the circumburst environment quantitatively.
A prominent feature of the off-axis afterglow light curves in the stratified medium is their shallow rising. If the viewing angle is large, we will observe a long and shallow rising. It is widely known that the off-axis structured jet also can give rise to a long and shallow rising in the light curve, which is used to explain the light curve of GRB 170817A (e.g., Lazzati et al. [42], Ghirlanda et al. [72], Kasliwal et al. [73]). This rising mechanism is due to the combined effect of the jet deceleration and the increase in the energy density and Lorentz factor at larger angles from the LOS (but closer to the jet core) where the emission is observable at later times. The two risings can be similar, but the underlying models may be distinguished from the late light curves. In the Newtonian or deep Newtonian phases, the decay indices of the light curves depend on the power-law slopes (s) of CBM and the power-law indexes of the electron distribution (see Table 5). The power-law indexes (p) of the electron distribution can be obtained by the observed spectrum. Thus, we can derive the power-law slopes of CBM, which, in turn, can be used to test whether the rising phase is entirely due to the stratified medium or not since the rising phase can also depend on the s and p in most cases. If a burst from a structured jet happens in the stratified CBM, both rising mechanisms are responsible for the light curves. The late afterglow light curve is crucial to diagnosing the jet structure and the density profile of CBM. A challenge is the late afterglow could be too faint to be well detected unless it happens close to us. The upcoming X-ray, optical, or radio survey (e.g., Einstein Probe, Yuan et al. [74]; the Large Synoptic Survey Telescope, LSST Science Collaboration et al. [75]) will possibly detect more off-axis bursts and provide more information on the GRB circumburst environment for us.

Conflicts of Interest:
The authors declare no conflict of interest.