Modeling Gamma-ray burst Afterglow observations with an Off-axis Jet emission

Gamma-ray bursts (GRBs) are fascinating extragalactic objects. They represent a fantastic opportunity to investigate unique properties not exhibited in other sources. Multi-wavelength afterglow observations from some short- and long-duration GRBs reveal an atypical long-lasting emission that evolves differently from the canonical afterglow light curves favoring the off-axis emission. We present an analytical synchrotron afterglow scenario, and the hydrodynamical evolution of an off-axis top-hat jet decelerated in a stratified surrounding environment. The analytical synchrotron afterglow model is shown during the coasting, deceleration (off- and on-axis emission), and the post-jet-break decay phases, and the hydrodynamical evolution is computed by numerical simulations showing the time evolution of the Doppler factor, the half-opening angle, the bulk Lorentz factor, and the deceleration radius. We show that numerical simulations are in good agreement with those derived with our analytical approach. We apply the current synchrotron model and describe successfully the delayed non-thermal emission observed in a sample of long and short GRBs with evidence of off-axis emission. Furthermore, we provide constraints on the possible afterglow emission by requiring the multi-wavelength upper limits derived for the closest Swift-detected GRBs and promising gravitational-wave events.

GRB 140903A (Troja et al. 2016) and GRB 160821B (Lü et al. 2017;Stanbro & Meegan 2016). On the other hand, several searches for afterglow emission around the closest bursts ( 200 Mpc) reported by the Burst Alert Telescope (BAT) instrument onboard the Neil Gehrels Swift Observatory (Dichiara et al. 2020) and the GW events by Advanced LIGO and Advanced VIRGO detectors (Abbott et al. 2021;The LIGO Scientific Collaboration et al. 2021) have been performed without successful, but setting multi-wavelength upper limits.
The density profile of the medium surrounding a burst has been addressed previously in different contexts. For instance, for the case of SNe, Chevalier (1982) studied the interaction of an adiabatic flow in a circumstellar density profile for Type II SNe of the form ∝ r −2 . Ever since, subsequent authors have adopted this proposal for modeling the circumstellar medium, and it has even been applied in the research of different SN types. Examples of such studies are those by Blondin et al. (1996), Soderberg et al. (2006), Kotak et al. (2004), and Chevalier (1984), among others. Nevertheless, a generalization of this power law has also been considered. A later study, Moriya & Tominaga (2012) showed that the diversity might explain the spectral diversity of Type II luminous SNe in the density slope of the surrounding dense wind. To this effect, they proposed a wind density structure in the form ∝ r −k . They noticed that the ratio of the diffusion timescale in the optically thick region of the wind and the shock propagation timescale after the shock breakout strongly depends on the stratification parameter k, which led to differences in the spectral SN evolution. On the other hand, the requirement of a stratified environment condition has been proposed in some cases for modeling the multi-wavelength afterglows, such as in work by Yi et al. (2013). The authors analyzed more than one dozen GRBs and concluded that the circumburst environment could be neither a homogeneous nor a stellar-wind medium but something in between, with a general density distribution with a stratification parameter in the range 0.4 ≤ k ≤ 1.4. A more recent example of an analysis that links SN and GRB emission with a stratified environment is the one by Izzo et al. (2020), in which the authors studied SN 2020bvc. They found an excellent agreement with the GRB-associated, broad-lined Ic SN 1998bw; thus, it was categorized as a young broad-lined Ic SN. They also noted that its X-ray light curve was consistent with simulations of an off-axis GRB afterglow in a stratified medium with k = 1.5; thus, this event represented the first case of an off-axis GRB discovered via its associated SN. It was later argued by Ho et al. (2020b), however, that such a model would predict an 8.5 GHz radio light curve several orders of magnitude more luminous than their measurements. Nevertheless, they stated that an off-axis jet could not be ruled out, and future radio observations would be needed.
In this work, we extend the analytical synchrotron afterglow scenario of the off-axis top-hat jet used to describe the multiwavelength afterglow observations in GRB 170817A (see Fraija et al. 2019b) adding several ingredients. Here we consider i) the circumburst external medium as stratified with a profile density ∝ r −k with k in the range of 0 ≤ k < 3, ii) the synchrotron radiation in self-absorption regime, iii) the afterglow emission during the transition from off-axis to on-axis before the lateral expansion phase (relativistic phase), iv) the hydrodynamical evolution computed by numerical simulations and v) a fraction of electrons accelerated by the shock front. We apply the current model to describe the delayed non-thermal emission observed in GRB 080503, GRB 140903A, GRB 150101B, GRB 160821B, and SN 2020bvc, and to provide constraints on the possible afterglow emission using multi-wavelength upper limits associated with the closest Swift-detected sGRBs and the promising GW events. This paper is arranged as follows: Section 2 presents the analytical synchrotron scenario and an hydrodynamical evolution of an off-axis top-hat jet decelerated in a stratified surrounding environment. In Section 3, we apply the proposed analytical model to describe the multi-wavelength observations of a sample of bursts and provide constraints to other ones. In Section 4, we present our conclusions.
2. OFF-AXIS TOP HAT MODEL 2.1. Radiative model Accelerated electrons in forward-shock models are distributed in accordance with their Lorentz factors (γ e ) and are described by the electron power index p as N (γ e ) dγ e ∝ γ −p e dγ e for γ m ≤ γ e , where γ m = m p /m e g(p)ε e (Γ − 1)ζ −1 e is the minimum electron Lorentz factor with Γ is the bulk Lorentz factor, m p and m e the proton and electron mass, ε e is the fraction of energy given to accelerate electrons, ζ e denotes the fraction of electrons that were accelerated by the shock front (Fan & Piran 2006) and g(p) = p−2 p−1 . The comoving of the magnetic field strength in the blastwave B 2 /(8π) = ε B e is derived from the energy density e = [(γΓ + 1)/(γ − 1)](Γ − 1)n(r)m p c 2 withγ the adiabatic index (Huang et al. 1999) and its respective fraction given to magnetic field (ε B ). Hereafter, we adopt the unprimed and prime terms to refer them in the observer and comoving frames, respectively. The termγ is the adiabatic index and n(r) = A k r −k =Ṁ W 4πvW r −k , where v W is the wind velocity andṀ W is the mass-loss rate. The sub-index k lies in the range 0 ≤ k ≤ 3, with k = 0 the constant-density medium (A 0 = n), and k = 2 the stellar wind ejected by its progenitor (A 2 A W 3 × 10 35 cm −1 ) where A W is the density parameter. The cooling electron Lorentz factor is γ c = (6πm e c/σ T )(1 + Y ) −1 Γ −1 B −2 t −1 , where σ T is the cross-section in the Thomson regime and Y is the Compton parameter (Sari & Esin 2001;Wang et al. 2010). The synchrotron spectral breaks and the synchrotron radiation power per electron in the comoving frame are given by ν i = q e /(2πm e c)γ 2 i B and P ν m √ 3q 3 e /(m e c 2 )B , respectively, with hereafter the subindex i = m and c for the characteristic and cooling break, and the constants q e and c the elementary charge and the speed of light, respectively (e.g., see Sari et al. 1998;Fraija 2015). The synchrotron spectral breaks in the self-absorption regime are derived from ν a,1 = ν c τ 3 5 0,m , ν a,2 = ν m τ 2 p+4 0,m and ν a,3 = ν m τ 3 5 0,c with the optical depth given by τ 0,i 5 3−k qen(r)r B γ 5 i with r the shock radius (Panaitescu & Mészáros 1998). Considering the total number of emitting electrons N e = (Ω/4π) n(r) 4π 3−k r 3 and also taking into account the transformation laws for the solid angle (Ω = Ω /δ 2 D ), the radiation power (P νm = δ D /(1 + z)P ν m ) and the spectral breaks (ν i = δ D /(1 + z)ν i ), the maximum flux given by synchrotron radiation is Weinberg 1972) is the luminosiy distance, r = δ D /(1 + z)Γβct is the shock radius, and δ D = 1 Γ(1−µβ) is the Doppler factor with µ = cos ∆θ, β = v/c with v the velocity of the material, and ∆θ = θ obs − θ j is given by the viewing angle (θ obs ) and the half-opening angle of the jet (θ j ). We assume for the cosmological constants a spatially flat universe ΛCDM model with H 0 = 69.6 km s −1 Mpc −1 , Ω M = 0.286 and Ω Λ = 0.714 (Planck Collaboration et al. 2016 We consider the dynamical equations proposed by Huang et al. (1999Huang et al. ( , 2000. The dynamical evolution of the relativistic outflow into the circumburst medium can be described by where c s = γ(γ − 1)(Γ − 1)c 2 / (1 +γ(Γ − 1)) withγ ≈ (4Γ + 1)/3Γ, M ej is the initial value of the ejected mass and is the radiative efficiency with = 0 in the adiabatic regime and = 1 in the fully radiative regime. The previous equations are consistent with the self-similar solution during the ultra-relativistic (Blandford-McKee solution;Blandford & McKee 1976) and the Newtonian phase (Sedov-Taylor solution), respectively, and consider the beaming effect (Rhoads 1999).
The observed quantities are integrated over the equal arrival time surface (EATS) determined by (Waxman 1997)

Analytical approach
During the coasting phase (before the deceleration phase), the relativistic outflow is not affected by the circumburst medium, so the bulk Lorentz factor is constant Γ cp = Γ 0 and the radius evolves as r = cβ 0 t/[(1 + z)(1 − β 0 µ)] with β 0 = Γ 2 0 − 1/Γ 0 . During the deceleration phase, the relativistic outflow transfers a large amount of its energy to the circumstellar medium driving a forward shock. Considering the adiabatic evolution of the forward shock with an isotropic equivalent-kinetic energy E = 4π 3 m p c 2 A k r 3 Γ 2 (Blandford-McKee solution;Blandford & McKee 1976) and a radial distance r = cβ of t/[(1 + z)(1 − β of µ)], the bulk Lorentz factor evolves as with β of = Γ 2 of − 1/Γ of . The deceleration time scale t dec can be defined using eq. 4. As the bulk Lorentz factor becomes Γ 1/∆θ, the observed flux becomes in our field of view. During the on-axis emission, the bulk Lorentz factor in the adiabatic regime evolves as with the radius r 2β on Γ 2 on ct/(1 + z) and β on = Γ 2 on − 1/Γ on , before the outflow enters to the post-jet-break decay phase (Γ 1/θ j ). During the post-jet-break decay phase, the bulk Lorentz factor evolves as with the shock's radius and velocity given by r 2β le Γ 2 le ct/(1 + z) and β le = Γ 2 le − 1/Γ le , respectively. We summarize the evolution of the bulk Lorentz factor as The minimum and cooling electron Lorentz factors, the spectral breaks, the maximum flux and the synchrotron light curves during the coasting, deceleration (off-and on-axis) and post-jet-break decay phases are shown in Appendix A and Tables 1, 2 and 3. Figure 1 presents examples of the time evolution of the Doppler factor, the half-opening angle, the bulk Lorentz factor, and the deceleration radius, all in a constant circumburst medium (k = 0). The solid lines characterize the numerical simulations, while the dashed lines stand-in for the theoretical approximations, which are detailed in Appendix A. Lines in black color corresponds to quantities observed off-axis and in gray color the ones observed on-axis.

Comparison of our model with previous simulations
For the case of δ D , the upper left panel shows an excellent agreement between the simulations and the model at early times up to t obs ≈ 10 2 s. After this point, there are very slight variations between both solutions, specifically in the steepness of the rise and fall of this quantity. The numerical simulation predicts a sharper peak, while the theoretical approximation anticipates a wider profile. Despite this slight difference, both curves drop in parallel, showing that their delinquent behavior follows the same power law. In the case of the half-opening angle θ j (lower left panel), as in the case of δ D , the early time evolution is the same between both curves. The difference in their behavior is first made apparent at t obs ≈ 10 2 s. At this time, the theoretical solution presents a faster rise than the simulation. The initial variation between both curves keeps increasing and becomes substantial at times ≥ 10 6 s.
The upper right panel compares the bulk Lorentz factors when the emission is off-and on-axis. For the off-axis case, the same features mentioned in the last two panels are repeated, namely perfect agreement at early times and variations at late times. However, the on-axis emission is different as there are slight differences initially, but both solutions tend to be the same as time progresses. Finally, the right panel on the bottom presents the shock's radius for the on-and off-axis cases. For the on-axis curves, the simulation and the theoretical approximation differ early, but they unite late. There is also a difference in their shapes at the beginning, as the theoretical curve is a straight line, which corresponds to a single power law. At the same time, the simulation shows that there is a transition from a flatter curve to a steeper one. The case of the off-axis emission presents similar differences. In general, the simulation and the theoretical variables exhibited in Figure 1 are in good agreement.

Analysis and description of the synchrotron light curves
The analytical synchrotron afterglow model during the coasting, deceleration (off-and on-axis emission), and the post-jetbreak decay phases are shown in Appendix A. Table 1 shows the evolution of the synchrotron light curves, and Table 2 shows closure relations of synchrotron radiation as a function of k during the coasting, deceleration (off-and on-axis emission), and the post-jet-break decay phases in the stratified environment. We can note in both tables that a break is expected around the transition time between fast-and slow-cooling regimes during the deceleration phase when the afterglow emission is seemed off-axis, but not during another phase. Figures 2 and 3 show the predicted synchrotron light curves produced by the deceleration of the off-axis top-hat jet in the circumburst medium described by a density profile with k = 0, 1.0, 1.5 and 2.0, respectively. Panels from top to bottom correspond to radio wavelength at 6 GHz, optical at R-band and X-rays at 1 keV for typical values of GRB afterglow. 1 These figures display that regardless of the viewing angle, the X-ray, optical, and radio fluxes increase gradually, reach the maximum value, and finally, they begin to decrease. It can be observed that as the viewing angle increases, the maximum value moves to later times. For the chosen parameters, in most cases the maximum value lies in the range (10 −1 − 10 2 ) days for 0 ≤ θ obs ≤ 60 deg. This is different, however, when the medium is like a stellar wind, as the right column of Figure 3 shows that the maximum happens later, namely in the range (10 −1 − 10 4 ) days. Figures 2 and 3 show that the bump is less evident in the light curves with 1.0 ≤ k ≤ 2.0. Then, a clear rebrightening in a timescale from days to weeks with GW detection could be associated with the deceleration of the off-axis jet launched by a BNS or BH-NS merger. The synchrotron fluxes in all these panels lie in the slow-cooling regime, although for a different set of parameters with values (for instance, the equivalent kinetic energy E = 10 54 erg, the uniform-density medium n ≈ 1 cm −3 , the equipartition parameters ε e = 0.1 and ε B = 10 −4 ), these would lie in the fast-cooling regime. Table 4 shows the evolution of the density parameter in each cooling condition of the synchrotron afterglow model. For instance, the synchrotron light curve as a function of the density parameter in the cooling condition ν a,1 < ν m < ν < ν c is given by F ν : ∝ A α k k , with α k = p+5 4 , 11−p 4 , 2 4−k and 3−p 4(3−k) , for the coasting, deceleration (off-and on-axis) and post-jet-break decay phases, respectively. Any transition between a stratified environment and density-constant medium could be detected during the post jet-break phase and the deceleration phase when the afterglow emission (for this cooling condition) is seemed on-axis but not off-axis. Table 4 shows that in general during the coasting phase and the deceleration phase when the afterglow emission is seemed off-axis cannot be detected a transition between different environments. Table 4 displays that synchrotron fluxes do not depend on the density parameter when these are observed at the lowest-energy (ν < min{ν a,1 , ν a,2 , ν a,3 }) and the highest-energy (max{ν m , ν c } < ν) frequencies when the afterglow emission is seemed off-and on-axis, respectively. During the post jet-break phase, any transition between a stratified environment and density-constant medium will be better observed in low-energy frequencies, such as radio wavelengths.
The uniform medium with constant density is expected for sGRBs associated with binary compact objects, and the stratified medium with non-constant density is expected for lGRBs related to massive stars with different evolution at the end of their lives (Ramirez-Ruiz et al. 2005;van Marle et al. 2006). For instance, an external stratified medium with 0.4 ≤ k ≤ 1.4 was found by Yi et al. (2013) after modeling the afterglow emission in a GRB example, and a density profile with k > 2 was reported by Kumar et al. (2008) after studying the accretion of the stellar envelope by a BH as the possible origin of the plateau phase in X-ray light curves. Giblin et al. (1999) analyzed the prompt gamma-ray emission in the BATSE 2 detected burst GRB 980923. The light curve exhibited a main prompt episode lasting ∼ 32 s followed by a smooth emission tail that lasted ∼ 390 s. The authors found that the spectrum in the smooth tail evolved as the synchrotron cooling break t −0.52±0.12 , concluding that the gamma-ray emission was associated with the afterglow evolution and also it had begun during the prompt gamma-ray episode. Afterward, spectra analyses of GRB tails were done to identify early afterglows (Barthelmy et al. 2005;Yamazaki et al. 2006). We could identify the off-axis synchrotron emission from an off-axis outflow analyzing the spectral break evolution. In this case, the spectral breaks of the synchrotron radiation generated from the deceleration of the off-axis jet in the relativistic phase evolve as a ν m ∝ t − 6−k 3 and ν c ∝ t 2+k 2 , respectively. For instance, with k = 1 the characteristic and cooling breaks evolve as ν m ∝ t − 5 3 and ν c ∝ t 3 2 , respectively, which fully different to those breaks that evolve when the afterglow emission is on-axis (e.g., ν m ∝ t − 3 2 and ν c ∝ t − 1 6 ).

A SAMPLE OF SOME BURSTS WITH EVIDENCE OF OFF-AXIS AFTERGLOW EMISSION
3.1. Multi-wavelength Observations GRB 080503 -On 2008 May 3 at 12:26:13 UTC the Swift-BAT instrument detected and triggered on the short burst GRB 080503 Ukwatta et al. 2008). The prompt episode was evaluated in the (15 − 150) keV energy range and reported with a duration of (0.32 ± 0.07) s, while its observed flux was measured to be (1.2 ± 0.2) × 10 −7 erg cm −2 s −1 . The Swift X-ray Telescope (XRT) instrument and Chandra ACIS-S satellite also performed subsequent observations in the X-ray band Perley et al. 2009). For the case of Swift-XRT, it took data on the burst in the timeframe from ∼ 82 s to 1 day after the initial Swift-BAT trigger. On the other hand, Chandra measured the burst throughout two observational campaigns, the first from 2008 May 07 (19:18:23 UTC) to 08 (04:09:59 UTC), during which an X-ray flux was detected, which coincided with the location of the optical afterglow. The second campaign took place from 2008 May 25 (18:11:36) to 26 (03:04:28), during which the X-ray source was monitored. There was a lack of detection, but constraining limits were provided. Several efforts were also taken to monitor this burst in the optical energy range. Observations and upper limits were obtained with the Swift Ultra-Violet Optical Telescope (UVOT) instrument (UVOT; Brown & Mao 2008), the Hubble Space Telescope (HST) using the Wide-Field Planetary Camera (WFPC2) in the F606W, F450W and F814W bands Perley et al. 2008f,b,c), and with the Keck-I telescope equipped with LRIS (Perley et al. 2008e). The Gemini-N observatory also observed the burst using the Multi-Object Spectrograph (GMNOS) through the g, r, i and z optical bands and NIRI through the Ks band (Perley et al. 2008a;Bloom & Perley 2008). Regarding the radio energy band, the Karl G. Jansky Very Large Array (VLA) was used to observe GRB 080503 at a frequency of 8.46 GHz without any detection but providing a 3-sigma upper limit (Frail & Chandra 2008).
GRB 140903A -The Swift-BAT detected GRB 140903A at 15:00:30 UT on 2014 September 14 ). This burst was located with coordinates R.A.= 238.036 • and Dec = +27.578 • (J2000). The BAT light curve exhibited a peak with a duration of 0.45 s (Cummings 2014). The XRT instrument began observing the position of this burst 59 s after the BAT trigger and monitored the X-ray afterglow until this emission faded below the detector sensitivity threshold. The MAXI-GSC observed the position reported by BAT 12 s after the trigger. Although no object was detected, upper limits were derived at the (4−10) keV energy range. The Chandra satellite began observing the X-ray afterglow ∼ 2.7 s after the BAT trigger (Sakamoto et al. 2014). This burst was observed in the optical r-band (20.4 ± 0.5) mag by the 2-m telescope Faulkes Telescope North ∼ 15.6 hours after the BAT trigger (Dichiara et al. 2014). The Nordic Optical Telescope (NOT) observed the field of this burst, reporting an optical emission of (20.1 ± 0.5) and ( NaID in absorption and H-beta in emission at the usual redshift of z = 0.351. Later, the detection of optical variability, together with a coincident radio detection (Fong 2014), confirmed the host association of this redshift ).
GRB 150101B -The Swift-BAT and Fermi-GBM detected GRB 150101B at 15:23:35 and 15:24:34.468 UT on 2015 January 01, respectively (Cummings 2015;Stanbro 2015;Burns et al. 2018). Data analysis of Swift-BAT revealed a bright source, constraining the location at R.A.=188.044 • and Dec: -10.956 • (J2000). The γ-ray pulse in the (15 − 150) keV band consisted of a single pulse with duration and fluence of (0.012 ± 0.001) s and F γ = (6.1 ± 2.2) × 10 −8 erg cm −2 , respectively (Lien et al. 2016). Recently, Burns et al. (2018) presented a new analysis of fine timescales revealing a two-component structure; a short hard spike followed by a longer soft tail. The total duration of the prompt episode shown as a two-component structure was 0.08 ± 0.93 s and the fluence was (1.2 ± 0.1) × 10 −7 erg cm −2 for the main peak and (2.0 ± 0.2) × 10 −8 erg cm −2 for the soft tail. The Chandra X-ray Observatory ACIS-S reported two observations, 7.94 and 39.68 days after the BAT trigger, with durations of ∼ 4.1 hours each one (Fong et al. 2016). Optical/near IR observations and upper limits were collected with Magellan/Baade using IMACS with r, g, i and z optical bands, Very Large Telescope (VLT) equipped with the High Acuity Wide field K-band Imager I (HAWK-I) with the J, H, K and Y optical bands (Fong 2015;Levan et al. 2015), TNG using NICS with the J optical band (D'Avanzo et al. 2015), Gemini-S instrumented with GMOS (r band) (Fong et al. 2015), UKIRT with the instrument WFCAM (J and K bands) and HST using the WFC3 with the F160W and F606 W bands (Fong et al. 2016). Spectroscopic observations in the wavelength range 530 − 850 nm of 2MASX J12320498-1056010, revealed several prominent absorption features that could associate GRB 150101B to an early-type host galaxy located at z = 0.1343 ± 0.0030 (Levan et al. 2015;Fong et al. 2016).
SN 2020bvc -On 2020 February 04 14:52:48 UTC SN 2020bvc was first detected by the ASAS-SN Brutus instrument using the g-Sloan filter with a reported location of R.A.=14 h : 33 m : 57.024 s and Dec=+40 • : 14 : 36.85 (J2000). It was associated to the host galaxy UGC 09379, with a redshift of z = 0.025235 (Stanek 2020). A later report confirmed this association and redshift and based on the blue featureless continuum and the absolute magnitude at the discovery of -18.1 classified this event as a young core-collapse supernova (Hiramatsu et al. 2020). A later analysis (Perley et al. 2020) of the spectrum obtained with the SPRAT spectrograph on the 2 m robotic Liverpool Telescope showed excellent agreement with the GRB-associated, broad-lined SN Ic 1998bw, thus it was categorized as a young broad-lined Ic SN. It was noticed to have an extremely fast rise by a steep decay in the two days following the first detection. It was later shown to rise towards a second peak. The decay temporal index reported by Izzo et al. (2020) was α dec = 1.35±0.9. Twelve days later, on February 16, the Very Large Array (VLA) observed the position of SN 2020bvc and detected a point source with a flux density of 66 µJy in the X-band and a luminosity of 1.3 × 10 27 erg s −1 Hz −1 (Ho 2020). On the same day, a 10ks Chandra observation was obtained. The data was reduced, and the spectrum fitted with a power-law source model with a flux of approximately 10 −14 erg cm −2 s −1 (Ho et al. 2020a).

Analysis and Discussion
GRB 080503 -We apply the Bayesian statistical approach of Markov-Chain Monte Carlo (MCMC) simulations to determine the best-fit values of the parameters that characterize the multi-wavelength afterglow observations with their upper limits (e.g., see Fraija et al. 2019c). A set of eight parameters, {E, n, p, ε B , Γ, ε e , θ j and θ obs }, is required by our synchrotron afterglow model evolving in a constant-density medium to describe the multi-wavelength observations. A total of 17300 samples and 4100 tuning steps is used to describe the entire dataset. Normal distributions are used to characterize all the parameters. Corner plots illustrate the best-fit values and the median of the posterior distributions of the parameters, as shown in Figure 5, respectively. In this figure, the best-fit values are highlighted in green, and the median of the posterior distributions are presented in Table 5.
The multi-wavelength afterglow observations, together with the fits computed using the synchrotron off-axis model evolving in a constant-density medium are shown in Figure 4. The synchrotron light curves obtained with the same electron population and displayed in the optical (R and g) and X-ray bands support the scenario of one-emitting zone inside an off-axis outflow. The fact that the beaming cone of synchrotron radiation reaches our line of sight is compatible with the best-fit value of the viewing angle θ obs = 15.412 +0.268 −0.269 deg and the re-brightening in all bands at ∼ one day. The best-fit value of the electron spectral index p = 2.319 +0.049 −0.049 matches the typical value observed when forward-shocked electrons radiate by synchrotron emission (see, e.g. Kumar & Zhang 2015). It confirms that these multi-wavelength observations originate during the afterglow.
The best-fit value of the constant-density medium n = 4.221 +0.102 −0.103 × 10 −2 cm −3 indicates that GRB 080503 took place in a medium with very low density comparable to an intergalactic density environment with ∼ 10 −3 cm −3 . The best-fit values of bulk Lorentz factor (Γ = 2.939 +0.101 −0.078 × 10 2 ) and the equivalent-kinetic energy (E = 2.156 +0.294 −0.295 × 10 52 erg) indicate that synchrotron radiation is emitted from a narrowly collimated outflow. Perley et al. (2009) analyzed the multi-wavelength observations at ∼ one day. The authors dismissed the kilonova-like emission proposed by Li & Paczyński (1998) and gave an afterglow explanation, pointing out that the X-ray and optical data had similar evolutions. They hypothesized that the late optical and X-ray bumps might be explained in a slightly off-axis jet or a refreshed shock. The faint afterglow compared to the intense prompt emission could be described in the very low circumburst medium. Gao et al. (2015) claimed that under certain requirements on the bulk Lorentz factor and the beaming angle of the relativistic jet, refreshed shocks in the synchrotron forward-and reverse-shock scenario could adequately describe the late re-brightening in GRB 080503. Finally, Gao et al. (2015) proposed that the late optical peak was due to the emission from a magnetar-powered "mergernova", and the X-ray hump from magnetic dissipation of the magnetar dipole spin-down luminosity. According to our analysis, these observations at ∼ one day are consistent with the afterglow emission found off-axis. It is worth noting that measurements of the linear polarization of the optical emission could discriminate between an on-and off-axis scenario and gravitational waves would provide information about the progenitor (a merger of NS -NS, BH -NS or a stable NS. GRB 140903A -We once again conducted MCMC simulations with a set of eight parameters, {E, n, p, ε B , Γ, ε e , θ j and θ obs }, to find the best-fit values that describe the multi-wavelength afterglow observations with their upper limits. A set of eight parameters, {E, n, p, ε B , Γ, ε e , θ j and θ obs }, is required. To represent the entire observations in this scenario, a total of 16200 samples and 4150 tuning steps were used. Figure 7 displays the best-fit values and the median of the posterior distributions of the parameters. In Table 5, the best-fit values are shown in green, and the median of the posterior distributions is presented. The multi-wavelength observations of GRB 140903A are shown in Figure 6, together with the fits derived using the synchrotron off-axis model evolving in a homogenous density. The best-fit values of the viewing angle θ obs = 5.162 +0.271 −0.267 deg and the half-opening angle θ j = 3.210 +0.080 −0.081 deg indicate that the relativistic outflow was slightly off-axis. The best-fit value of the constant density medium n = 4.219 +0.102 −0.101 × 10 −2 cm −3 indicates that GRB 140903A took place in a medium with very low density comparable to an intergalactic density environment. The best-fit values of bulk Lorentz factor (Γ = 3.627 +0.100 −0.100 × 10 2 ) and the equivalent-kinetic energies (E = 3.163 +0.290 −0.296 × 10 52 erg) indicate that synchrotron emission is radiated from a narrowly collimated jet. The best-fit values of the microphysical parameters ε e = 5.104 +0.298 −0.302 × 10 −2 and ε B = 4.050 +0.909 −0.921 × 10 −3 , and the best-fit value of the spectral index of the shocked electrons 2.073 +0.048 −0.050 , are similar to those reported in Troja et al. (2016). This spectral index is in the typical range to those accelerated in forward shocks (see, e.g. Kumar & Zhang 2015), thus reaffirming the afterglow as its origin. Troja et al. (2016) reported and analyzed the afterglow observations of GRB 140903A for the first two weeks. Applying the fireball scenario, the authors demonstrated that this burst was caused by a collimated jet seen off-axis and was also connected with a compact binary object. The X-ray "plateau" seen in GRB 140903A was attributed to the energy injection into the decelerating blast wave by Zhang et al. (2017). The authors then modelled the late afterglow emission, which required a half-opening angle of ≈ 3 deg, similar to the value found with our model. GRB 140903A was formed in a collimated outflow observed off-axis that decelerates in a uniform density, according to our findings. GRB 150101B -We use MCMC simulations with the eight parameters used for GRB 080503 and GRB 140903A to find the best-fit values that model the X-ray afterglow observations with the optical upper limits. To represent the entire data in this case, a total of 15900 samples and 4400 tuning steps are used. Figure 9 displays the best-fit values and the median of the posterior distributions of the parameters. In Table 5, the best-fit values are shown in green, and the median of the posterior distributions is presented.
The X-ray, optical and radio observations and upper limits, as well as the fit obtained using the synchrotron off-axis model evolving in homogeneous density are shown in Figure 8. The left-hand panel shows the light curves at 1 keV (gray), R-band (blue), F606W filter (orange) and F160W filter (dark green), and the right-hand panel displays the broadband SEDs at 2 (red) and 9 (green) days. The red area corresponds to the spectrum of AT2017gfo, which is adapted by Troja et al. (2018). The best-fit values of the viewing angle θ obs = 14.114 +2.327 −2.179 deg and the half-opening angle θ j = 6.887 +0.662 −0.682 deg can explain the lack of X-ray emission during the first day. The best-fit value of the constant density medium n = 0.164 +0.021 −0.021 × 10 −2 cm −3 indicates that GRB 150101B, like other short bursts, happened in an environment with very low density. The best-fit values of bulk Lorentz factor (Γ = 4.251 +0.468 −0.453 × 10 2 ) and the equivalent-kinetic energy (E = 1.046 +0.120 −0.124 × 10 52 erg) suggest that synchrotron afterglow emission is emitted from a narrowly collimated jet. The values of the spectral index of the electron population, the circumburst density, the microphysical parameters, and the viewing angle disfavor the isotropic cocoon model reported in Troja et al. (2018) and are consistent with the values of an outflow when homogeneous density is taken into account. The best-fit values of the spectral index of the shocked electrons 2.150 +0.217 −0.215 is similar to those reported in synchrotron afterglow models (see, e.g. Kumar & Zhang 2015).
GRB 160821B -We require MCMC simulations with eight parameters used in the previous bursts to find the best-fit values that describe the multi-wavelength afterglow observations with their upper limits. To characterize the complete data in this scenario, a total of 17300 samples and 7400 tuning steps are used. Figure 11 exhibits the best-fit values and the median of the posterior distributions of the parameters. In Table 5, the best-fit values and the medians of the posterior distributions are presented.
The multi-wavelength observations since 0.2 days after the GBM trigger are shown in Figure 10, together with the fits found requiring the synchrotron off-axis model evolving in a homogeneous density. The left-hand panel exhibits the light curves at 1 keV (gray), z-band (purple), F606W filter (dark red), R-band (salmon), F110W filter (cyan), F160W filter (blue sky), K s -band (blue), X-channel (olive) and C-channel (emerald green), and the right-hand panel shows the broadband SEDs of the X-ray, optical and radio afterglow observations at 2 h (red), 2 days (blue) and 4 days (green). The shaded areas in blue and green correspond to blackbody spectra with decreasing temperatures firstly suggested in Jin et al. (2018) and then confirmed by Lamb et al. (2019) and Troja et al. (2019a). The kilonova emission is the most natural explanation for the "new" radiation component. The Swift-XRT-UVOT data were received from the public database from the official Swift website. 3 The C-band displays radio data, the white, v, b, u, UVW1, UVW2, and UVM2 bands display Swift-UVOT data, and the 1 keV band displays XRT data. Using the conversion factor proposed in Evans et al. (2010), the flux density of XRT data is extrapolated from 10 keV to 1 keV. The best-fit values of circumburst density n = 0.869 +0.093 SN 2020bvc -We use MCMC simulations to find the best-fit values of the parameters that describe the X-ray afterglow observations. A set of eight parameters, {E, A w , p, ε B , Γ, ε e , θ j and θ obs }, is used to describe the X-ray observations. To characterize the complete data in this scenario, a total of 17100 samples and 7200 tuning steps are used. Figure 13 exhibits the best-fit values and the median of the posterior distributions of the parameters. In Table 5, the best-fit values are shown in green, and the medians of the posterior distributions are presented. We assumed a stratified medium with a parameter k = 1.5, consistent with the proposal by Izzo et al. (2020). Our best fit values, however, are slightly different. We propose that the emission is due to an off-axis jet that is ≈ 5 times more energetic and with half of the off-axis angle compared to the values of Izzo et al. (2020). This discrepancy is due to a different choice of the electron velocity distribution index p, as our MCMC simulation suggested p = 2.313 +0.037 −0.035 , in contrast with the value of p = 2 used by the cited authors. Figure 12 shows the X-ray observations of SN 2020bvc with the best-fit synchrotron light curve generated by the deceleration of an off-axis jet in a medium with stratification parameter k = 1.5.
Our results are consistent with the X-ray observations before and after ∼ 4 days since the trigger time; when the observed flux increases and decreases, respectively. Initially, the flux increases with a minimum rise index of α m,ris > 1.65 and later the observed flux decreases with α dec = −1.35 ± 0.09 (Izzo et al. 2020). The allowed value of the minimum rise index is estimated considering a simple power-law function, the X-ray upper limit and the maximum flux. For instance, given the best-fit value of p = 2.313 +0.037 −0.035 and for 0 < k < 1.5, the temporal rise and decay indexes are 1.65 ± 0.03 ≤ α ris ≤ 4.53 ± 0.03 and α dec = −1.24 ± 0.03, respectively, for ν c < ν (see Table 1). For k > 1.5, the expected rise index would be α ris < 1.65, which cannot reproduce the X-ray observations. The synchrotron scenario from on-axis outflow in a very-low density environment is disfavored for p ∼ 2. The closure relations of synchrotron on-axis model from an outflow decelerating in a stratified environment for k in general can be estimated. During the coasting and the deceleration phases the bulk Lorentz factor evolves as Γ ∝ t 0 and Γ ∝ t k−3 8−2k , respectively. Therefore, the synchrotron flux during the slow-cooling regime evolves as F ν ∝ t 12−k(p+5) 4 ν − p−1 2 for ν m < ν < ν c and ∝ t 8−k(p+2) 4 ν − p 2 for ν c < ν during the coasting phase, and F ν ∝ t − 12(p−1)+k(5−3p) 4(4−k) ν − p−1 2 for ν m < ν < ν c and ∝ t − 3p−2 4 ν − p 2 for ν c < ν, during the deceleration phase. It is worth noting that for the cooling condition ν c < ν and with a value of p = 2.5, the temporal evolution is only consistent for k 0.3 The best-fit values of bulk Lorentz factor (Γ = 2.291 +0.100 −0.100 × 10 2 ), the equivalent-kinetic energy (E = 2.38 +0.01 −0.01 × 10 51 erg) and the half-opening angle θ j = 2.121 +0.078 −0.079 deg indicate that synchrotron emission is produced from a narrowly collimated outflow decelerating in an external medium. The best-fit values of the viewing angle θ obs = 12.498 +0.268 −0.281 deg and the half-opening angle θ j = 2.121 +0.078 −0.079 deg are consistent with the lack of early multi-wavelength observations. Lü et al. (2012) discovered a correlation between the bulk Lorentz factors and the isotropic gamma-ray luminosities in a sample of GRBs. Fan et al. (2012) showed that the correlation of these parameters were consistent with the parameters predicted in the photospheric emission model. Figure 14 shows the diagram of the bulk Lorentz factors and the isotropic gamma-ray luminosities of sGRBs described in this work (red) with those (gray) reported in Lü et al. (2012). For off-axis sGRBs, we found an empirical correlation Γ = a(L/10 52 erg) b with a = (3.27 ± 0.39) × 10 2 and b = −(4.9 ± 2.0) × 10 −2 .  100) keV energy range. The duration and measured fluence in the energy range of 15 -150 keV were 128 ± 16 ms and (5.9 ± 3.2) × 10 −8 erg cm −2 , respectively (Parsons et al. 2005). Levan et al. (2008) provided the specifics of Swift's deep optical and infrared observations. According to the authors, no X-ray nor optical/IR afterglow was detected to deep limits, and no residual optical or IR emission was observed.
GRB 070810B -Swift-BAT was triggered by GRB 070810B at 15:19:17 UTC on August 10, 2007, with a reported location of R.A. = 00 h 35 m 46 s and Dec = +08 • 50 07 (J2000) with a positional accuracy of 3 (Marshall et al. 2007). The KANATA 1.5-m telescope, the Xinlong TNT 80 cm telescope, the 2-m Faulkes Telescope South, the Shajn 2.6 m telescope, and the Keck I telescope (HST) conducted follow-up observations after the initial detection, which are summarized in Bartos et al. (2019). From the whole observational campaign, only the Shajn telescope detected a source inside the error box of GRB 070810B (Rumyantsev et al. 2007).
GRB 080121 -Swift-BAT was triggered by GRB 080121 at 21:29:55 UTC on January 21, 2008, with a reported location of R.A. = 09 h 08 m 56 s , Dec = +41 • 50 29 (J2000) and a positional accuracy of 3 . In the (15 − 150) keV energy range, the duration and measured fluence were (0.7 ± 0.2) s and (3 ± 2) × 10 −8 erg cm −2 s −1 , respectively, according to Cummings & Palmer (2008). Follow-up observations were carried out 2.3 days following the burst using the Swift/UVOT and the Swift/XRT. Within the Swift-BAT error circle, however, no X-ray afterglow candidate or sources were discovered (Cucchiara & Schady 2008;Troja & Burrows 2008). Within the Swift-BAT error circle, two galaxies were found, indicating a redshift of z ∼ 0.046 for GRB 080121, however the isotropic energy released would be several orders of magnitude lower than usual short-hard bursts (Perley et al. 2008d).   Table 6. We report two values for each parameter; the upper values correspond to synchrotron light curves on the left-hand panels and lower ones on the right panels. For typical values of GRB afterglow reported in Table 6, the synchrotron emission is ruled out for a density of n = 1 cm −3 , but not for n = 10 −2 cm −3 . The value of the uniform-density medium ruled out in our model is consistent with the mean value reported for sGRBs (e.g., see Berger 2014). It is worth noting that for values of ε B < 10 −4.3 and ε e < 10 −1.2 , the value of the density n = 1 cm −3 would be allowed.

Promising GW events in the third observing run (O3) that could generate electromagnetic emission
The Advanced LIGO and Advanced Virgo produced 56 non-retracted alerts of gravitational waves candidates during the O3 run, covering almost one year of operations (from 2019 April 01 to 2020 March 27). Nevertheless, three of them have a probability of being terrestrial larger than 50%. The O3 observing run was divided into two epochs associated to "O3a" (from April 01 to September 30) and "O3b" (from November 01 to March 27, 2020). The GW events in the O3a and O3b runs are listed in the Gravitational Wave Transient (GWTC-2) Catalog 2 (Abbott et al. 2021) and (GWTC-3) Catalog 3 (The LIGO Scientific Collaboration et al. 2021), respectively, where from GCNs there were 31 and 22 candidate events discovered during O3a and O3b respectively. The promising candidates that are consistent with a source with m 2 < 3M and that could generate electromagnetic emission are GW190425, GW190426 152155, GW190814 in GWTC-2 (Abbott et al. 2021) and GW191219 163120, GW200105 162426, GW200115 042309, GW200210 092254 in GWTC-3 (The LIGO Scientific Collaboration et al. 2021). Table 7 enumerates the main characteristics of these candidates. Figure 16 shows the multi-wavelength upper limits of GW events in GWTC-2 and GWTC-3 consistent with a source with m 2 < 3M and that could generate electromagnetic emission and the synchrotron light curves from an off-axis outflow decelerating in a constant-density medium with n = 1 cm −3 (left panels) and 10 −2 cm −3 (right panels). The synchrotron light curves are presented at 1 keV (green), R-band (brown) (from Becerra et al. (2021a)) and 3 GHz (red). The parameter values used are E = 5 × 10 50 erg, θ j = 3 deg, θ obs = 6 deg, Γ = 100, ε e = 0.1, p = 2.5, ζ e = 1.0 and ε B = 10 −2 . The left-hand panels associated to the S190425z, S190426c, S190814bv and S200115j events show that a uniform-density is ruled out for n = 1 cm −3 , but not for n = 10 −2 cm −3 . For instance, we can note that in the panel related to S190426c the synchrotron emission at 1 keV, at the R-band and at 3 GHz is above the upper limits around ∼ one day for n = 1 cm −3 , and in the panel associated to S190814bv the synchrotron flux at the R-band and 3 GHz is above the upper limits, but not at 1 keV. The value of the constant-density medium ruled out in our model is consistent with the densities derived by Dobie et al. (2019); Ackley et al. (2020); Gomez et al. (2019) using distinct off-axis jet models. We need further observations on different timescales and energy bands to derive tighter constraints.

CONCLUSIONS
We have extended the synchrotron off-axis model presented in Fraija et al. (2019b) initially proposed to describe the multiwavelength afterglow observations in GRB 170817A. In the current model, we have considered i) the circumburst external medium as stratified with a profile density ∝ r −k with k in the range of 0 ≤ k ≤ 3, ii) the synchrotron radiation in selfabsorption regime, iii) the afterglow emission during the transition from off-axis to on-axis before the lateral expansion phase (relativistic phase), iv) the hydrodynamical evolution computed by numerical simulations and v) a fraction of electrons accelerated by the shock front. The time evolution of the Doppler factor, the half-opening angle, the bulk Lorentz factor, and the deceleration radius computed by numerical simulations are in good agreement with those derived with an analytical approach. The advantage of this general approach (with a density profile A k ∝ r −k ) is that this model allows us to take into account not only both a homogeneous medium (k = 0) and a wind-like medium (k = 2) but regions with non-standard stratification parameters, such as k = 1.0, 1.5 or 2.5.
We have calculated the synchrotron light curves and presented the closure relations in a stratified environment, including the self-absorption regime for all cooling conditions during the coasting, deceleration (off-and on-axis emission), and the post-jetbreak decay phases. We have noted that a break is expected around the transition between fast-and slow-cooling regimes during the deceleration phase when the afterglow emission seems off-axis, but not during other stages. We have analyzed the behavior of the flux for different parameters of the density distribution. We have noticed that the behavior during the relativistic phase approaches flatness as the stratification parameter is raised. On the other hand, we have shown that the time evolution of the light curves after the jet break is independent of k, so this model gives freedom to explain the early-time evolution of the radiation while keeping the long-time results invariant. Furthermore, we have derived the change of the density parameter in the entire phase. We have shown that: In general, during the coasting and the deceleration phases, when the afterglow emission seems off-axis, a transition between different environments cannot be detected, contrary to the post-jet-break phase and the deceleration phase with on-axis emission. The synchrotron fluxes do not depend on the density parameter when these are observed at the lowest and the highest frequencies when the afterglow emission seems off-and on-axis, respectively. During the post-jet-break phase, any transition between a stratified environment and density-constant medium will be better observed in low frequencies, such as radio wavelengths.
In particular, we have applied our model to describe the delayed non-thermal emission observed in a sample of bursts with evidence of off-axis emission. In accordance with the best-fit values for sGRBs, we found that i) the constant-density medium required to model the multi-wavelength observations is low (10 −2 ≤ n ≤ 0.4 cm −3 ), indicating that the central engine are located in a low density circumstellar medium, ii) the spectral indexes of the shocked electrons (2.1 ≤ p ≤ 2.3) are in the range of those reported after the description of the afterglow observations (see, e.g. Kumar & Zhang 2015;Fraija et al. 2017;Becerra et al. 2019aBecerra et al. ,b,c, 2021b, and iii) the half-opening angles (1.5 ≤ θ j ≤ 6.6 deg), the bulk Lorentz factors (130 ≤ Γ ≤ 450) and the equivalent-kinetic energies (0.2 ≤ E ≤ 5.7 × 10 52 erg) provide evidence of narrowly collimated outflow expanding into a constant-density environment. The previous results confirm that the multi-wavelength observations are emitted from the GRB afterglow and indicate a merger of compact objects (two NSs or NS-BH) as possible progenitors of these bursts. The low-density medium agrees with the larger offsets of sGRBs compared with lGRBs. Regarding SN 2020bvc, we found that an atypical stratification parameter of k = 1.5 is required, supporting the CC-SN scenario. The best-fit values of the half-opening angle 2.121 +0.078 −0.079 deg, the viewing angle 12.498 +0.268 −0.281 deg, the equivalent-kinetic energy 2.38 +0.11 −0.10 × 10 51 erg and the bulk-Lorentz factor 2.291 +0.100 −0.100 × 10 2 provide evidence of the scenario of off-axis GRB afterglow.
We have applied the current model to provide constraints on the possible afterglow emission using multi-wavelength upper limits associated with the closest Swift-detected sGRBs (< 200 Mpc) and the promising GW events. We have shown that the value of the constant-density medium is ruled out, which is consistent with the mean value of densities reported in for sGRBs (e.g., see Berger 2014) and those derived by S190814bv event (e.g., see Dobie et al. 2019;Ackley et al. 2020;Gomez et al. 2019) using different models. To derive tighter constraints, further observations on timescales from months to years post-merger phase are required.