GRB Afterglow of the Sub-relativistic Materials with Energy Injection

Sub-relativistic materials launched during the merger of binary compact objects and the core-collapse of massive stars acquire velocity structures when expanding in a stratified environment. The remnant (either a spinning magnetized neutron star (NS) or a central black hole) from the compact-object or core-collapse could additionally inject energy into the afterglow via spin-down luminosity or/and by accreting fall-back material, producing a refreshed shock, modifying the dynamics, and leading to rich radiation signatures at distinct timescales and energy bands with contrasting intensities. We derive the synchrotron light curves evolving in a stratified environment when a power-law velocity distribution parametrizes the energy of the shock, and the remnant continuously injects energy into the blastwave. As the most relevant case, we describe the latest multi-wavelength afterglow observations ($\gtrsim 900$ days) of the GW170817/GRB 170817A event via a synchrotron afterglow model with energy injection of a sub-relativistic material. The features of the remnant and the synchrotron emission of the sub-relativistic material are consistent with a spinning magnetized NS and the faster"blue"kilonova afterglow, respectively. Using the multi-band observations of some short-bursts with evidence of kilonova, we provide constraints on the expected afterglow emission.

In the sub-relativistic regime, the interaction of the decelerated material with the surrounding circumburst medium has been explored to interpret the multi-wavelength observations in timescales from days to several years as synchrotron afterglow models (e.g., see Wijers et al. 1997;Dai & Lu 1999;Huang et al. 1999;Livio & Waxman 2000;Huang & Cheng 2003;Sironi & Giannios 2013;Barniol Duran & Giannios 2015). In most of the cases, the isotropic-equivalent kinetic energy of the materials launched during the coalescence of BCOs and the CC-SNe has been described by a power-law (PL) velocity distribution (e.g., see Tan et al. 2001, and the references therein).
The canonical X-ray light curve exhibits a typical shape that consists of four distinct PL segments ∝ t −α with a great flare (e.g., see Zhang et al. 2006;Nousek et al. 2006). The initial steep decay with a temporal decay index 3 α 5, the normal decay phase with 0.9 α 1.5 and the late abrupt decay with 1.9 α 2.4 have been explained in terms of the end of main episode, the standard synchrotron forward-shock model (Sari et al. 1998) and the post-jetbreak decay phase (Vaughan et al. 2006), respectively. There is, however, one segment that occurs between the end of the prompt phase and the normal decay, a shallower than usual decay with −0.1 α 0.7. This so-called "plateau" phase has been explained in several scenarios such as continuous energy injection from the central engine (either a spinning magnetized NS or a central BH) into the blastwave (Barthelmy et al. 2005;King et al. 2005;Dai et al. 2006;Perna et al. 2006;Proga & Zhang 2006;Burrows et al. 2005;Chincarini et al. 2007;Dall'Osso et al. 2017;Becerra et al. 2019a,b), stratified ejecta Jin et al. 2007;Hascoët et al. 2017), ejecta with a wide range of Lorentz factors (Rees & Mészáros 1998;Kumar & Piran 2000;Sari & Mészáros 2000;Zhang & Mészáros 2002;Fraija et al. 2019c) and variation on microphysical parameters (Fan & Piran 2006;Ioka et al. 2006;Fraija et al. 2020). Several modeling efforts of multi-wavelength afterglows evoking energy injection by central engine have been widely explored (e.g., see Laskar et al. 2015;Zhao et al. 2021;Becerra et al. 2019a;Pereyra et al. 2022;Fraija et al. 2021b). For instance, Laskar et al. (2015) described GRB 100418A, GRB 100901A, GRB 120326A and GRB 120404A and found that the majority of the kinetic energy of the relativistic jet in each burst was carried by slow-moving ejecta, thus indicating a correlation between the injection rates and the Lorentz factor distribution.
On August 17, 2017, a gravitational wave (GW) event (GW170817; Abbott et al. 2017a,b) was linked with a faint gamma-ray prompt emission of GRB 170817A (Goldstein et al. 2017;Savchenko et al. 2017). Immediately, GRB 170817A was followed by an extensive observational campaign covering radio, optical, and X-ray bands (e.g., see Abbott et al. 2017b;Kisaka et al. 2017;D'Avanzo et al. 2018a, and references therein). The observations of the non-thermal spectrum of GRB 170817A gathered during the first ≈ 900 days after the initial merger were analyzed by several authors. It was shown that they were consistent with synchrotron forward-shock emission generated by the deceleration of an off-axis structured jet with an opening angle θ j ≈ 5 • that was observed from a viewing angle in the range of 15 • ≤ θ obs ≤ 25 • Kasliwal et al. 2017b;Lamb & Kobayashi 2017Resmi et al. 2018;Margutti et al. 2017a;Lazzati et al. 2017;Gottlieb et al. 2018b;Fraija et al. 2019b;Gottlieb et al. 2018b;Hotokezaka et al. 2018;Fraija et al. 2019c). In some proposed models, the off-axis structured jet is formed with an off-axis jet with a cocoon (Lazzati et al. 2017;Gottlieb et al. 2018b;Fraija et al. 2019b) and a shock breakout Hotokezaka et al. 2018;Fraija et al. 2019c;Urrutia et al. 2021). Recently, Hajela et al. (2021) analyzed the latest X-ray and radio observations of GRB 170817A collected with the Chandra X-ray Observatory, the Very Large Array (VLA), and the MeerKAT radio interferometer about 3.3 years after the initial merger. These new observations did not agree with the best-fit synchrotron curves from the off-axis jet model, thus reporting evidence of a new X-ray emission component. Given these contrasting observations, the authors offered the solution to explain this phenomenon in the context of either radiation from accretion processes on the compact-object remnant or a KN afterglow.
The study of properties of KNe has a great impact for the present-day field of study, especially given the link between GWs, short-duration gamma-rays, and KN emission. The merger of two NS associated with GW170817, GRB 170817A, and AT 2017gfo have provided the needed tools to predict the KN emission and its characteristics. While the prompt episode and the early afterglow are produced in internal and external shocks by an ultra-relativistic and extremely collimated jet, the KN transient is associated with a quasiisotropic emission easier to detect at angles far away from those emitted from a collimated jet (Metzger & Berger 2012). Despite the advantageous prospects for detection, only four transient events with different brightness to AT 2017gf have been classified as KNe. They are associated to GRB 050709 (Jin et al. 2016), GRB 060614 (Yang et al. 2015), GRB 130603B (Tanvir et al. 2013;Berger et al. 2013) and GRB 160821B (Kasliwal et al. 2017a;Troja et al. 2019).
Recently, Fraija et al. (2021a) presented the afterglow light curves generated by the deceleration of sub-relativistic masses ejected from the merger of BCOs and the death of massive stars. The authors assumed that a PL velocity distribution describes the isotropic-equivalent kinetic energy of these masses and that the sub-relativistic ejected masses were decelerated, in turn, by a stratified-density environment. As a particular case, to explain the multi-wavelength observations of the gravitational event GW170817/GRB 170817A at ∼ 900 days, they constrained the parameter space of the synchrotron light curves of a sub-relativistic mass ejected during the merger of two NSs and decelerated in a constantdensity environment. The synchrotron radiation of the subrelativistic material was consistent with the faster "blue" KN afterglow. Inspired by the new observations of this GW event at 3.3 years after the initial merger (Hajela et al. 2021), in this paper, we extend the synchrotron model presented in Fraija et al. (2021a) including the continuous energy injection from the central engine (either a spinning magnetized NS or BH remnant) into the blastwave through a numerical approach and analytic arguments. In addition, we apply the current model to potential candidates of sGRB events with evidence of a KN. The paper is organized as follows: Section 2 presents the dynamical evolution of the afterglow when the central engine continuously injects energy into the blastwave. We show an analytical solution and numerical approach. In Section 3, we show a synchrotron model with energy injection from a spinning magnetized NS and BH remnants. Section 4 shows the analysis of the multi-wavelength light curves using typical values of the GRB afterglow. In Section 5, we apply our model to several potential candidates including GW170817/GRB 170817A, and finally, in Section 6, we summarize. We consider the convention Q x = Q 10 x in c.g.s. units and 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).

Synchrotron radiation
We consider electrons accelerated in the forward shock which evolves in a stratified external environment with a density profile described by ρ(r) = A k r −k =Ṁ W 4πvW r −k , whereṀ W is the mass-loss rate and v W is the wind velocity. There are two common choices of stratification, the value of k = 0 corresponds to the constant-density medium (A 0 = n), and k = 2 is associated to the density of the stellar wind ejected by its progenitor (A 2 A W 3 × 10 35 cm −1 ) with A W the density parameter. We assume that the shockedaccelerated electrons in the forward shocks can be described by a single PL energy distribution dn dγe ∝ γ −p e for γ e ≥ γ m with p the spectral index and γ m the Lorentz factor of the lowest-energy electrons.
The post-shock magnetic field evolves as B ∝ A plicitly written in Fraija et al. 2021a).
Given the synchrotron spectral breaks and the maximum flux density, the synchrotron light curves in each cooling condition evolve as and We want to emphasize that the synchrotron spectrum is always in the slow-cooling regime (νm < νc) and the spectrum in the fastcooling regime (νc < νm) is derived for completeness since it is not relevant for the timescales investigated here. It is worth noting that the peak flux density is always at the peak of the spectrum.
The velocity β in the case of the coasting and deceleration phase without energy injection was derived in Fraija et al. (2021a). In the following, we derive the evolution of the velocity when the central engine continuously injects energy into the blastwave.

Dynamical evolution
Energy injection by the central engine on the GRB afterglow can produce refreshed shocks. The luminosity injected from the central engine into the blastwave can be described by (e.g., Zhang et al. 2006 where q is the energy injection index, Linj is the initial luminosity and tc is the characteristic timescale. The isotropic-equivalent kinetic energy can be estimated as Given q = 1, the energy does not evolve with time and the standard synchrotron light curves are recovered (Sironi & Giannios 2013;Barniol Duran & Giannios 2015;Fraija et al. 2021a), and for q > 1, the decreasing value of isotropic-equivalent kinetic energy is not considered. The energy injection could be due to magnetic spin-down from a spinning magnetized NS (q = 0; Ruffert et al. 1997;Zhang & Mészáros 2001) and a fallback material onto a central BH (q ≤ 1) (Proga & Zhang 2006;Barthelmy et al. 2005;King et al. 2005;Dai et al. 2006;Perna et al. 2006;Proga & Zhang 2006;Burrows et al. 2005;Chincarini et al. 2007;Lei et al. 2013;Wu et al. 2013;Dall'Osso et al. 2017). It is relevant to mention that the spinning magnetized NS could also accrete material (e.g., see Metzger et al. 2018). The luminosity due to magnetic spin-down (L sd ) or accreting BH (LBH) scenario can be converted into flux through the efficiency in converting its spindown/accreting energy to radiation (η) and the beaming factor of the wind (f b = 1 − cos θj) with θj the half-opening angle. In both scenarios, the initial luminosity can be written as Linj = η f b Lj, and for typical and similar values of η and f b , Linj ≈ Lj with j = sd or BH.
During the deceleration phase, the ejected mass acquires a velocity structure, the velocity of matter in the front of the ejected mass is faster than the one that moves in the back (Sari & Mészáros 2000). Tan et al. (2001) studied the acceleration of the ejected mass with relativistic and sub-relativistic velocities. They found that the isotropic-equivalent kinetic energy in the sub-and ultrarelativistic limit can be expressed as a PL velocity distribution given by E k (≥ β) ∝ β −5.2 for β 1 and E k (≥ βΓ) ∝ (βΓ) −1.1 for βΓ 1 (with Γ = 1/1 − β 2 ), respectively. 7 Here, we consider the sub-relativistic regime, so the isotropic-equivalent kinetic energy distribution is given by whereẼ is the fiducial energy. We consider the power-law velocity distribution in the sub-relativistic regime with the values of α in the range 3 ≤ α ≤ 5.2 presented in Tan et al. (2001). It is worth noting that these values of α in the sub-relativistic regime have been widely used (e.g., see Hotokezaka & Piran 2015;Metzger 2017;Hajela et al. 2019).
The total isotropic-equivalent kinetic energy is given by the superposition of the energy injection (Eq. 5) and the energy distribution (Eq. 6). In the sub-relativistic regime, the ejected material is described by the Sedov-Taylor solution as where E β =Ẽβ −α and Et =Ê t tc 1−q withÊ = 1 1−q tc Linj for q < 1 and mp the proton mass.

Analytical solution
The Sedov-Taylor solution can be solved analytically in the asymptotic cases; Et E β and E β Et. Each limiting case leads to a different velocity; they are given by Both cases may be written with just one expression: 7 The polytropic index np = 3 is used.
where the case Et E β is obtained by setting q = 1 with E =Ẽ, while the case E β Et is obtained by setting α = 0 with E =Ê. It is worth mentioning that the deceleration time can be obtained from Eq. 9, and is presented in Section A.1. The blastwave radius (r = βt/(1 + z)) can be written as The standard equations in constant-density medium are recovered when q = 1 and α = 0 (i.e., β ∝ t − 3 5 and r ∝ t 2 5 ; Sironi & Giannios (2013)). The dynamics, spectral breaks and synchrotron light curves derived in Fraija et al. (2021a) are recovered for Et E β . Using Eqs. (9) and (10), we report in the Appendix the equations of the dynamics, the synchrotron spectral breaks, the flux density and the light curves in the fast-and slow-cooling regime. It is worth noting that the synchrotron spectrum can lie in the slow-or fast-cooling regime, depending on the parameter values.

Numerical approach: Comparison with analytic solution
We solve Eq. (7) numerically using the bisection method 8 and plot the evolution of the shock's velocity for different parameter values, as shown in Figure 1. The rows represent different choices of luminosity, namely the top one corresponds to a value of L inj = 10 43 erg s −1 and the lower one to L inj = 10 45 erg s −1 . The column on the left presents the velocity's development for different choices of the circumburst density profile ∝ r −k with k = 0, 1, 1.5, 2 and 2.5. The center column displays its transformation according to different values of the isotropic-equivalent kinetic energy of the outermost matter's PL distribution index with α = 3.0, 4.0 and 5.0. On the rightmost column, we show the evolution of the velocity for varying possibilities of the energy injection index, namely q = 0, 0.5 and 1. For the middle and right-hand panels, the solid lines represent a choice of constant-density medium (k = 0), while the dashed lines correspond to stellar wind (k = 2).
All panels show two distinct behaviours. At early times (before approximately 10 years), the forward shock expands into the circumburst medium uninhibited with a constant velocity β 0 , the so-called "coasting phase" which is represented in the panels from Figure 1, as horizontal lines. Once the shock has interacted with enough material, the coasting phase comes to an end and the deceleration phase commences. It is governed by the solution of equation (7) and can be seen in Figure 1 as the decrease in the velocity after its constant segment. 8 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.bisect.html Upon comparison of the middle panels between both rows, it can be seen that when the dominant contribution to the energy of the shock is the fiducial energyẼ (represented by the top row) the behaviour of the shock's velocity depends more strongly on the velocity distribution's PL index α. This is made apparent from the separation of the curves for different values of this parameter in the deceleration phase, which is more pronounced in the top row for both constant-density and wind-like environment.
On the other hand, the opposite behaviour is noticed when the same comparison is performed between the rightmost panels. In this case, the energy injection parameter q's variation is more easily observed whenÊ is the dominant component, as is made clear by the lower row in both types of medium considered. It is also apparent that when the energy that is injected into the shock is substantial, it may reach a quasi-constant value, where the energy injected equals the energy lost by interaction with the surrounding medium. This behaviour is best exemplified by the dashed q = 0 curve in the lower panel, where it can be seen that the shock's velocity becomes quasi-constant even at late times.

ENERGY INJECTION FROM A SPINNING MAGNETIZED NS AND AN ACCRETING BH
A spinning magnetized NS and BH remnants are created from the merger of a BCO system (Price & Rosswog 2006;Siegel et al. 2013;Kiuchi et al. 2014) and the death of a massive star (Usov 1992;Wheeler et al. 2000;Thompson et al. 2004;Bucciantini et al. 2007;Metzger et al. 2011). These remnants could accrete material and inject energy into the blastwave. The spinning magnetized NS or the central BH inject energy due to either a magnetic spin-down (Ruffert et al. 1997;Zhang & Mészáros 2001) or accretion (Proga & Zhang 2006;Barthelmy et al. 2005;King et al. 2005;Dai et al. 2006;Perna et al. 2006;Proga & Zhang 2006;Burrows et al. 2005;Chincarini et al. 2007;Lei et al. 2013;Wu et al. 2013;Dall'Osso et al. 2017), respectively. For instance,  summarized a variety of final remnants from a merger of BCOs, which might inject energy into the blastwave; a hyper-massive NS (HMNS) that collapses into a BH in short timescales, a non-spinning BH, and a spinning magnetized NS and BH. In the case of a spinning magnetized NS or BH, the authors reported that the late-time activity via spin-down power or accreting material could be expected up to timescales of years. In this case, the synchrotron light curves from ejected materials would be modified by energy injection into the blastwave.

A spinning magnetized NS with fall-back accretion
Rapid spinning magnetized NSs called "millisecond magnetars" are potential candidates for long and short GRBs. The energy reservoir of a millisecond magnetar is the total rotation energy which is given by where P is the spin period associated to an angular frequency Ω = 2π/P and I 1.3 × 10 45 M 3 2 ns,1.4 g cm 2 (Lattimer & Schutz 2005) is the NS moment of inertia with M ns = 1.4 M the NS mass. The merger of two NSs or CC-SN usually leave a fraction of the stellar progenitor bound to the NS. This fraction of material will begin to rotate into an accretion disk and to fall-back over a long period (Chevalier 1989;Rosswog 2007;Woosley & Heger 2012;Quataert & Kasen 2012). The fall-back accretion rate can be written as (Metzger et al. 2018) where M fb is the accreting mass over a characteristic fallback time t fb . Once the millisecond magnetar is formed, the NS might be subject to fall-back accretion. This accretion depends on the dipole magnetic moment (µ), and the Alfén (r m ), the co-rotation (r c ) and cylinder (r lc ) radii. The spin evolution of an accreting magnetar is given by (Piro & Ott 2011) where the spin-down terms due to the torque and the accretion are and respectively, with G N the gravitational constant, µ = BR 3 ns where R ns = 1.2 × 10 6 cm is the NS radius and B is the strength of the dipole magnetic field. For details, see Metzger et al. (2018). In this scenario, the isotropic-equivalent kinetic energy due to injection would be E t = L inj dt = η f b L sd dt. The left-hand panels of Figure 2 show the spin-down luminosity and the synchrotron forward-shock light curves from the millisecond magnetar with an initial spin period P 0 = 10 −3 s (Metzger et al. 2018), an accreting mass M fb = 0.8M (Metzger et al. 2018), a half-openig angle θ j = 30 • (Dhawan et al. 2020) and an efficiency η = 0.1 (Xiao & Dai 2019). The spin-down luminosities exhibit a "plateau" phase for an energy injection index of q = 0 (as indicated), and different time scales are observed for four different sets of parameters. The black solid curve corresponds to a magnetic field strength of B = 10 14 G and a characteristic fallback time of t fb = 10 3 s (Metzger et al. 2018), the gray solid line takes the values B = 10 16 G and t fb = 10 3 s, the black dashed one takes B = 10 14 G and t fb = 10 7 s and the gray dashed curve represents B = 10 16 G and t fb = 10 7 s. It is worth nothing that the gray dashed curve displays a spindown luminosity ≈ 10 42 erg s −1 during the "plateau" phase. Light curves show one or two "plateau" phases depending on the values of the magnetic field, spin period and the characteristic timescale of fallback. The first "plateau" related with a precursor is ∼ two order of magnitude less than the second one, which is associated with the prompt emission. For small values of the characteristic timescale of fallback and larger values of magnetic field, the light curves exhibit a 'plateau" in a timescale of seconds. Larger values of the characteristic timescale and small magnetic fields lead to a small luminosity during the "plateau" phase. The first "plateau" is explained when the Alvén radius is larger than the co-rotation one (r c r m ). In this case, the spin-down luminosity becomes . For a large value of t sd the spin-down luminosity decays very slowly. The second "plateau" can be explained once the equilibrium is reached and an analytic solution of Eq. 13 can be derived (see Metzger et al. 2018;Fraija et al. 2021b). In this case, the spindown luminosity becomes L sd ≈ 10 40.7 erg s −1 B where the typical values of the accreting mass M fb = 0.8M and the characteristic fall-back time t fb = 10 8 s are used. The profile of the spin-down luminosity (Eq. 17) agrees with those profiles shown in the upper left-hand panel in Figure 2. The black solid curve displays an uninterrupted drop during 10 3 days. The gray solid line presents a very similar behaviour to its black counterpart as it decays with no interruption with the same power law. The main difference between both curves being that the luminosity achieved by this solution is smaller by about two orders of magnitude. The black dashed curve shows an initial "plateau" phase for approximately 10 −1 days. Afterwards it decreases by about two orders of magnitude until it reaches its second plateau phase after a couple of days. This phase happens for roughly one month after which the luminosity drops once again with a power law slightly less steep than the one present in the previously mentioned solid curves. Finally, the gray dashed profile presents the same type of behaviour as its black analogue, however this solution starts off with larger values than the black line at early times and then drops below them during the black curve's first "plateau" phase. Afterwards, the gray profile remains below the black one.
The lower left-hand panel shows the light curves in this same scenario in three energy bands: radio (3 GHz), optical (R-band) and X-ray (1 keV). All curves assume evolution in a constant density medium k = 0 with n = 5.88×10 −3 cm −3 , E = 10 48 erg, z = 0.023, p = 2.15, B = 10 −3 , e = 10 −1 , α = 3.0 and β = 0.3 (Metzger 2017;Hajela et al. 2019;Fraija et al. 2021a;Hajela et al. 2021). The values of the magnetic field strength and characteristic time scale were chosen to be the same as the ones of the gray dashed line from the upper panel.
All three light curves grow with the same power law during the first 10 4 days. After this point in time only the X-ray profile reaches its maximum and starts to drop, while the other two continue to rise, albeit with a smaller slope. The panel shows that the flux density of a profile increases with the energy of the band, that is, that the X-ray curve has the largest flux density, while the radio one has the smallest value.

Fall-back material onto a BH
The Eddington luminosity for pure ionized heavy elements is where M BH = 2.3 M is the mass of the remnant BH, σ T is the Thompson cross section and Y e ≈ 0.4.
The fall-back material onto the remmant BH could create an accretion disk, powering a new material via Blandford-Znajek (BZ;Blandford & Znajek 1977), or neutrinoannihilation mechanism (Popham et al. 1999). The BZ jet power from a BH with mass M BH and angular momentum J BH can be described as (Lee et al. 2000) is the dimensionless spin parameter (e.g., see Tchekhovskoy et al. 2008 Wu et al. 2013). The termṀ BH is the accretion rate onto the BH given by (Kumar et al. 2008a,b) where τ vis is the viscous timescale, t 0 is the starting time of the accretion andṀ fb is the accretion rate described by (Chevalier 1989;MacFadyen & Woosley 1999;MacFadyen et al. 2001;Zhang et al. 2008) withṀ p and t p the fallback rate and the time at the peak, respectively. For details, see Wu et al. (2013). In this scenario, the isotropic-equivalent kinetic energy due to injection The right-hand panels of Figure 2 are analogous to their left-hand counterparts, which were explained in the previous subsection, but now in the context of the fall-back material onto a BH scenario. The upper panel presents the BZ luminosity with q = 2/3 (as indicated in the dashed gray line) with different time scales. This behavior corresponds to an accretion rate ofṀ fb ∝ t 1 2 . The change of the slope in the BZ luminosity is associated with the variation of the accretion rate fromṀ fb ∝ t 1 2 to ∝ t − 5 3 (see Eq. 21). It is worth noting that the BH torus system has no effect on the afterglow evolution at later times, which corresponds to an accretion rate of ∝ t − 5 3 (e.g., MacFadyen et al. 2001;Janiuk et al. 2004;Zhang et al. 2006). We consider typical values of starting time scale of accretion t 0 = 1 s, time scale at the peak t p = 10 3 s , BH mass M BH = 2.3 M and dimensionless spin parameter a = 0.7 Zhao et al. 2021). Once again, four different sets of parameters were considered. The black solid curve corresponds to a viscous timescale of τ vis = 10 9 s and a fallback rate oḟ M p = 10 −6 M s −1 (Zhao et al. 2021), the gray solid line takes the values 10 6 s and 10 −5 M s −1 , the black dashed one takes 10 9 s and 10 −4 M s −1 and the gray dashed curve represents 10 12 s and 10 −5 M s −1 . In a similar fashion to the left-hand panel, both of the solid curves have a similar behaviour. They both decrease following the same power law and the black profile is greater than the gray one by about two orders of magnitude. In the case of the dashed lines, the black solution decreases uninterrupted, but it changes its power law at around 10 4 days. The gray curve also decreases with no interruption and follows the same slope as the black one at early times, but it changes its slope at approximately 10 2 days. It is also interesting to note, that similarly to the left panel, the gray dashed profile starts off with a greater luminosity than its black counterpart, but it falls below it at later times, namely in the vicinity of 10 3 days.
The lower right-hand panel shows the synchrotron forward-shock light curves in this scenario in the same energy bands as in the panel to its left. The deceleration was once again assumed to be in a wind-like medium k = 2 with A W = 10 −2 cm −3 and the values of the viscous timescale and fallback rate were chosen to be the same as the ones of the gray dashed line from the upper panel. The rest of the parameters were selected asẼ = 10 48 erg, z = 0.023, p = 2.15, B = 10 −3 , e = 10 −1 , α = 3.0 and β = 0.3 (Metzger 2017;Hajela et al. 2019;Fraija et al. 2021a;Hajela et al. 2021). In all cases the BZ luminosity is much higher than the Eddington luminosity.
The optical and X-ray bands have the same exact behaviour, namely that their light curve reaches its peak very early, it remains constant until ∼ 10 days and from this moment onward its flux density decreases. It is worth to note that both of these curves follow the same power laws; the difference between them being that the optical band lies two orders of magnitude above the X-ray band. The case of the radio light curve is different, as it can be observed that it reaches its maximum later than the previously mentioned profiles. Its constant phase also lasts less time and its behaviour at late times is according to a power law with a steeper slope than the one of the other two energy bands.
Upon comparison of the two upper panels of Figure 2, it can be observed that, in general, the luminosities in the context of a BH decrease in time less steeply than the ones from a millisecond magnetar. It can also be noted that for a BH, there is no appearance of the so-called "plateau" phase, which is a characteristic that could discriminate between both scenarios.
In the case of the lower panels of the aforementioned Figure, it can be concluded that the flux density increases in smaller timescales in the BH scenario, as in this case the maximum is reached at ∼ 10 −1 days, while in the case of the magnetar it is obtained at ∼ 10 3 days. The same observation can be made for the drop of the flux, as it can be noticed from the Figure that the light curves in the right panel begin to decrease at ∼ 10 1 days, while on the left panel a drop is not apparent in none of the shown energy bands even at ∼ 10 5 days.
An analytic solution of Eq. 19 can be derived for t p < t In this case, the term e ± t τ vis ≈ 1 and therefore the BZ luminosity becomes withṀ p,−6 =Ṁ BH 10 −6 M s −1 . The profile of the BZ luminosity (Eq. 22) agrees with those profiles shown in the upper right-hand panel in Figure 2.
It is worth mentioning that the neutrino-annihilation luminosity could not describe late-time activities (e.g., see Zhao et al. 2021) and therefore, it cannot be considered in this work.  Figure, in ascending order, presents the radiative behaviour produced by the interaction between the sub-relativistic ejecta and its surrounding medium, which is described by a density profile A k r −k with k = 0, 1, 1.5, 2 and 2.5, respectively. Panels from top to bottom correspond to the electromagnetic bands in radio at 1.6 GHz, optical at the R-band and X-rays at 1 keV for E = 10 49 erg, B = 10 −3 , e = 10 −1 and z = 0.023 (Metzger 2017;Hajela et al. 2019;Fraija et al. 2021a;Hajela et al. 2021). The left-hand panels show the light curves for p = 2.6 with α = 3.0, 4.0 and 5.0, and the right-hand panels show the light curves for α = 3.0 with p = 2.2, 2.8 and 3.4.
We present the predicted synchrotron light curves in the previously mentioned energy bands for typical values of GRB afterglows. Most light curves peak on timescales from several months to a few years, which is in agreement with the observations of some SNe, such as SN 2014C (Margutti et al. 2017b) and SN2016aps (Nicholl et al. 2020). There is an outlier, however, present in the flux in the radio band for p = 3.4 and α = 3.0, which reaches its maximum in a timescale of a couple of days for stratified media when k > 1.5. The synchrotron light curves are shown using a stratified medium with density profile ∝ r −k with k = 0, 1, 1.5, 2 and 2.5, which covers both long and short GRB progenitors. The constant-density medium (k = 0) is usually related to short GRBs which stem from the merger of two NSs, while the stratified medium (1 ≤ k ≤ 2.5) is only associated to long GRBs from dying massive stars with different mass-loss evolution. For instance, in their article, Yi et al. (2013) investigated the evolution of the emission of forward-reverse shocks propagating in a medium described by the previously mentioned PL distribution. They applied their model to 19 GRBs and found that the density profile index took values of 0.4 ≤ k ≤ 1.4, with a typical value of k ∼ 1. This value was also obtained by Liang et al. (2013), who analyzed a bigger sample of 146 GRBs.
All Figures show that, during the early stages of the evolution, there is an increase in the flux. During this epoch, when the sub-relativistic material decelerates in a constant-density medium, the light curve grows steeply. For more stratified media, this growth is not as evident, as the example presented in Figure 7, where the early-time behaviour is more gradual. That is to say, the amount of time during which the flux grows depends on the stratification of the surrounding medium. This is exemplified by Figures 3 and 4, where the rise in the flux is evident, while in subsequent Figures this growth in not observed, which means that this phase happens much quicker as the stratification is increased. This result implies that if such a flattening or rebrightening at timescales from months to years in the light curve was observed together with GW detection, then this would be associated with the deceleration of a sub-relativistic material launched during the merger of two NSs. Otherwise, we show that an observed flux that gradually decreases on timescales from months to years could be associated with the deceleration of a sub-relativistic material launched during the death of a massive star with different mass-loss evolution at the end of its life. It is worth noting that all these results are for onaxis observers, and for an off-axis observer the flux would decrease as the viewing angle between the material and the observer increases. For a relativistic off-axis component in the outflow, the spectral breaks and the maximum flux are corrected by the Doppler factor (δ D ) as where Pν m = δD/(1 + z)P ν m is the radiation power and Ne = (Ω/4π) n(r) 4π 3−k r 3 is the total number of emitting electrons with r = δ D /(1 + z)Γβct the shock radius and the transformation law for the solid angle as Ω = Ω /δ 2 D . The Doppler factor is defined as δ D = 1 Γ(1−µβ) with µ = cos ∆θ, Γ the bulk Lorentz factor and ∆θ = θ obs − θ j given by the viewing angle (θ obs ) and the half-opening angle of the jet (θ j ). Table 2 shows the evolution of the density parameter in each cooling condition of the synchrotron afterglow model. For instance, the synchrotron light curve in the slow-cooling regime as a function of the density parameter is given by k for ν c < ν. Any variation of the density will be better observed in low-energy frequencies, such as radio. Additionally, it shows that variations of the density profile index are more apparent in the radio light curve when compared to the other fluxes in the other energy bands. This change in the density profile is also more easily appreciated for large values of α, namely 4.0 and 5.0. Therefore, a transition between density profiles will be more easily observed in the radio band with high values of the velocity distribution parameter. Figure 8 shows several synchrotron light curves in order to compare the effects of the continuous injection of energy into the blastwave. It is divided into two columns, the one on the left considers a constant-density medium (k = 0) with the millisecond magnetar remnant, while the one on the right takes into account a stellar wind (k = 2) with the scenario of fall-back material onto a BH. The panels from top to bottom correspond to radio (1.6 GHz), optical (R-band) and X-ray (1 keV), respectively. Every panel shows two curves, the dashed line represents the evolution of the flux density with energy injection, and the solid line stands for an evolution with no injection of energy.

Comparison: with and without energy injection
It can be seen in all panels on the left that regardless of injection of energy or not, both solutions are the same at early times. However, all solid curves remain several orders of magnitude below their dashed counterparts after the evolution has reached times of approximately 20 days. There also seems to be a difference in timescales when the maximum is reached, as the dashed curves reach it within the limits of the plot, while the solid lines continue their upward trend. This is noted by observing that for late times (> 10 4 days) the rise in the light curve is steeper when there is energy injection, while the solution represented by the dashed lines reaches its peak with a smaller slope and then begins to drop.
On the other hand, the panels on the right present contrasting behaviour, as both the case with energy injection and the case without reach the same peak flux, the only difference between them being the decay of the light curves. For a stellar-wind environment, the late-time behaviour is the opposite of the constant medium case, namely that the flux density drops and this decay is less sharp when there is energy injection.
Upon comparison between the right and left columns, it can also be observed that the early-time behaviour is different in the optical and X-ray bands. Regardless of whether there is energy injection or not, the flux density in the constant medium increases until its peak. For a stellar-wind environment, however, the light curve remains constant and begins to decay very slowly. This behavior is due to the change in the energy injection. At early times it evolves as q = 2/3 (up to 10 3 days), and at later times withṀ fb ∝ t − 5 3 the BH torus system has no effect on the afterglow evolution (e.g., Mac-Fadyen et al. 2001;Janiuk et al. 2004;Zhang et al. 2006). It is worth noting that after 10 3 days both curves decrease with the similar slope, as expected. We emphasize that for q = 1 the standard synchrotron light curves are recovered (e.g., Zhang et al. 2006).

SYNCHROTRON EMISSION FROM DIFFERENT EJECTED MATERIALS AND APPLICATIONS
It is believed that sub-energetic GRBs are quasi-spherical explosions whose dominating components are the subrelativistic materials, which contribute approximately 99.9% of the explosion's energy. The mildly relativistic materials, on the other hand, correspond to only ≈ 0.1% (e.g., see Margutti et al. 2014;Modjaz et al. 2020). There is wide agreement in the community that the origin of sGRBs and lGRBs is closely related to the merger of BCOs and the death of massive stars leading to KNe and SNe, respectively. In addition to KN and SN materials, other types of materials are launched into the circumstellar medium with different velocities and, as such, will contribute at distinct timescales in distinct energy bands with contrasting intensities. In the following we will give a brief introduction about the values of masses, the isotropic-equivalent kinetic energies and velocities of each decelerated material that is ejected during the merger of two NSs, namely the dynamical ejecta, the shock breakout material, the disk wind and the cocoon material.
Shock breakout material -A shock at the interface between the two NSs is formed the moment directly after their coalescence. This shock manages to break out from the NS core to the crust at sub-relativistic velocities (β in 0.25 e.g., see Kyutoku et al. 2014;Metzger et al. 2015). When the shocked material reaches half of the escape velocity it converts a fraction of the shock-heated internal energy into kinetic energy and it escapes the merger into a nearly vacuum region (for details see Kyutoku et al. 2014;Fraija et al. 2019c). The shock breakout material's properties depend on the mass, radius and velocity of the merger remnant. Numerical simulations indicate that the material mass, the kinetic energy and the velocities lie in the ranges of 10 −6 M ej 10 −4 M , 10 47 Ẽ 10 50.5 erg and βΓ 0.8, respectively (e.g., see Kyutoku et al. 2014;Metzger et al. 2015).
Disk wind -The coalescence of the NS binary will end in a tidal disruption and some of the material of the stars will be shed, forming an accretion disk around the central remnant. This component of the sub-relativistic material represents a significant portion of the total material mass and might dominate over other constituents (Siegel & Metzger 2017). The mass of the accretion disk will depend on the initial NS spins and is located in the range of 10 −3 M ej 0.3 M (Shibata & Taniguchi 2006;Hotokezaka et al. 2013). The disk's kinetic energy and the velocities lie in the ranges of 10 47 Ẽ 10 50 erg and 0.03 βΓ 0.1, respectively (e.g., see Dessart et al. 2009;Metzger & Fernández 2014;Fernández et al. 2015).
Cocoon material -The GRB jet will deposit energy as it travels through the neutrino-driven or magnetically driven wind (previously expelled during the merger of two NSs). The energy deposited laterally will produce a cocoon with an energy similar to that of the jet's electromagnetic emission. Murguia-Berthier et al. (2014) looked into the conditions required for cocoon formation as a function of the jet's luminosity. A weak cocoon emission was predicted, regardless of the magnitude of the jet's luminosity. In particular, when Nagakura et al. (2014) numerically examined a low-luminosity jet, the authors concluded that a hot cocoon enclosing the jet would form. The cocoon would break free and spread along the axis of the relativistic jet as soon as it would reach the shock-breakout material. The external pressure would then drop dramatically beyond the breakout material, allowing the cocoon to accelerate and expand relativistically until it became transparent. The material mass liberated in the cocoon, the kinetic energy and the velocities lie in the ranges of 10 −6 M ej 10 −4 M , 10 47 Ẽ 10 50.5 erg and 0.2 βΓ 10, respectively (e.g., see Nagakura et al. 2014;Murguia-Berthier et al. 2014;Lazzati et al. 2017Lazzati et al. , 2018Nakar & Piran 2017;Gottlieb et al. 2018a). Figure 9 shows the synchrotron light curves with energy injection by a spinning magnetized NS remnant and generated by materials ejected from the merger of two NSs such as the dynamical ejecta, the cocoon material, the shock breakout material and the wind ejecta. The synchrotron light curves correspond to the (from top to bottom) radio (1.6 GHz) and X-ray (1 keV) bands, respectively for k = 0. Taking into account the velocities of the wind and the shock breakout material, we also consider the trans-relativistic (TR; β ∼ 0.8) and Deep-Newtonian (DN; β ∼ 0.08) regimes, respectively. The TR and DN timescales during the deceleration phase are given in Appendix. The light curves are shown for n = 10 −2 cm −3 , α = 3, P 0 = 10 −3 s, B = 10 16 G, t fb = 5 × 10 9 s, B = 10 −2 , e = 10 −1 and p = 2.2 and the pair of values ofẼ = 10 50 erg and β = 0.2 for the dynamical ejecta, 10 48 erg and 0.3 for the cocoon material, 10 48.5 erg and 0.8 for the shock breakout material and 10 50 erg and 0.07 for the wind. Besides, we show for completeness the synchrotron afterglow radiation from an on-axis and off-axis outflow with a viewing angle of 30 • , which are specified in the top panel. The synchrotron light curves from the off-axis outflow are considered as detailed in Fraija et al. (2019b).
The disk wind, the dynamical ejecta, the cocoon and the shock breakout peak at 1.2 × 10 5 days, 7150 days, 522.5 days and 56.1 days, respectively. The total contri-bution of synchrotron light curves exhibits a brightening once the off-axis emission decreases at a few years, so that, depending on the parameter values and conditions, the synchrotron emission from the ejected materials could be detected or not (see, Murguia-Berthier et al. 2014;Nagakura et al. 2014, for the cocoon material).
The top panel shows the evolution of the integrated light curve in the radio band. The early-time behaviour of the radiation due to the jet's emission is emphasized by showing different solutions for distinct viewing angles. This panel shows that, as the viewing angle is increased, the light curve at early times decreases and takes longer to enter our line of sight. Once the jet enters on-axis, the behaviour is independent of the initial opening angle, however, the light curve becomes dominated by the emission from the other emitted materials.
This last point is highlighted in the subsequent lower panel, where each component's contribution is explicitly plotted. It can be seen that at timescales of ≈ 10 3 −10 4 days, the flux density is dominated by the emission from the dynamical ejecta and the shock breakout, while at late times ( 10 5 days) both the dynamical and the disk wind ejecta have the upper hand as the most influential constituents of the synchrotron light curve. Contrastingly, the cocoon emission lies a couple of orders of magnitude below the other contributions and is not observed. Figure 10 shows the synchrotron light curves with energy injection produced by a fall-back material onto a BH and generated when ejected materials and a relativistic outflow decelerate in a stellar wind environment. The synchrotron light curves correspond to the (from top to bottom) radio (1.6 GHz) and X-ray (1 keV) bands, respectively for k = 2, A 2 = 3 × 10 36 cm −1 (A W = 10), t 0 = 1 s, t p = 10 3 s, M BH = 2.3 M , a = 0.7, τ vis = 10 9 s, M p = 10 −6 M s −1 , α = 3.0, B = 5 × 10 −3 , e = 10 −1 and p = 3.2, and the pair of values 10 48 erg and 0.3 for the cocoon material, and 10 48.5 erg and 0.8 for the shock breakout material.
For a stellar-wind surrounding medium, the jet's contribution on the radio light curve is much shorter than for the case of a constant-density medium, as the top panel shows that it lasts less than a day, while the left-hand panels show that this phase has a duration 10 3 days. On the other hand, the light curve at 1 keV shows that the flux density is completely dominated by the jet's emission, as the cocoon's and the shock breakout's contributions lie several orders of magnitude below the jet's. Therefore, a high variation in the behaviour of the light curves while comparing the different energy bands could be a hint for an outflow evolving in a medium that is stratified. The synchrotron light curves at the radio bands show that, depending on the parameter values, the afterglow emission from the cocoon and shock breakout materials could be detected on timescales of days. However, if the outflow is chocked, radio fluxes could be detected on timescales of hours.
GRBs could be successful or chocked, this is determined by the range of values in the observables such as luminosity, duration and bulk Lorentz factor (e.g., see MacFadyen et al. 2001;Mészáros & Waxman 2001;Murguia-Berthier et al. 2014;Fraija 2014;Nagakura et al. 2014;Sobacchi et al. 2017;Bromberg et al. 2011). As an indication of this behavior, successful GRBs might be less frequent than choked ones, only limited by the ratio of SNe (types Ib/c and II) to lGRB rates (Totani 2003;Ando & Beacom 2005). Some SNe of type Ic-BL not connected with GRBs have been suggested to arise from events such as off-axis GRBs or failed jets (e.g., see Izzo et al. 2019Izzo et al. , 2020Beniamini et al. 2020). One exponent of such a possibility is the failed burst GRB 171205A, which besides being associated to SN 2017iuk, exhibited material with high expansion velocities β ∼ 0.4 interpreted as mildly relativistic cocoon material (Izzo et al. 2019). As another example, in the context of off-axis GRBs that enter on-axis after some time, Izzo et al. (2020) found that the X-ray observations from the nearby SN 2020bvc were consistent with the afterglow emission generated by an off-axis jet with viewing angle of 23 • when it decelerated in a circumburst medium with a density profile with k = 1.5. However, as no prompt emission was detected, the authors implied that this was a hint for the first orphan GRB detected through its associated SN emission.
As follows we apply the current model to describe the latest multi-wavelength afterglow observations ( 900 days) of the GW170817/GRB 170817A event, and using multiwavelength upper limits associated with i) promising GW events in GWTC-2 and GWTC-3 that could generate electromagnetic emission, ii) short-bursts with the lowest-redshifts (100 ≤ z ≤ 200 Mpc), and iii) evidence of KNe, we provide constraints on the possible afterglow emission.

GW170817/GRB 170817A event
GW radiation from the merger of two NSs (Metzger 2017) is expected together with a short gamma-ray prompt and an UV-optical-IR KN emission in timescales of ∼ 1 s and a few days, respectively (Li & Paczyński 1998;Rosswog 2005;Metzger et al. 2010;Kasen et al. 2013;Metzger 2017). A KN classified as "blue" and "red" is a transient powered by radioactive decay of unstable heavy nuclei via the rapid neutron capture (r-process) synthesized in merger ejecta. The "blue" KN situated at the polar regions has low opacity and fast velocity β 0.3 and the "red" KN positioned at the equatorial plane has high opacity due to the Lanthanidebearing matter and slower velocity β 0.1 (Metzger & Fernández 2014;Perego et al. 2014;Wanajo et al. 2014;Miller et al. 2019).
As follows, we present the GW170817/GRB 170817A observations, focusing on and describing the latest ones with synchrotron afterglow radiation from a sub-relativistic material, i.e., a KN afterglow emission, and when the spinning magnetized NS remnant is accreting and injecting energy into afterglow. We consider the characteristics of the "blue" KN as used in , and references therein). The observations of the nonthermal spectrum of GRB 170817A gathered during the first ≈ 900 days after the initial merger were analyzed by several authors and it was shown that they were consistent with synchrotron radiation from an off-axis structured jet decelerated in a constant-density medium. This relativistic jet observed from a viewing angle of 15 • ≤ θ obs ≤ 25 • was described with an opening angle θ j ≈ 5  Hajela et al. (2021) analyzed the latest X-ray and radio observations of GRB 170817A collected with the Chandra X-ray Observatory, the Very Large Array (VLA), and the MeerKAT radio interferometer about 3.3 years after the initial merger, and reported evidence of a new X-ray emission component. This new measurement was not in agreement with the synchrotron off-axis afterglow model in constantdensity. Given these contrasting properties, the authors offered the solution to explain this phenomena in the framework of either radiation from accretion processes on the compact-object remnant or a KN afterglow.

Analysis, Description and Discussion
The observations of the non-thermal spectrum of GRB 170817A gathered during the first ≈ 900 days after the initial merger have been modelled with synchrotron forwardshock emission generated by the deceleration of a relativistic off-axis jet Kasliwal et al. 2017b;Lamb & Kobayashi 2017Resmi et al. 2018;Margutti et al. 2017a), a cocoon (Lazzati et al. 2017;Gottlieb et al. 2018b;Fraija et al. 2019b) and a shock breakout Hotokezaka et al. 2018;Fraija et al. 2019c) in a constant-density medium. While the syncrotron radiation from an off-axis jet peaked at ≈ 110 − 130 days, the synchrotron radiation from the relativistic cocoon material with a bulk Lorentz factor (Γ c 4) increased gradually during the first weeks, after reached a maximum flux at ∼ 15 − 45 days and decreased afterwards (Lazzati et al. 2017;Fraija et al. 2019b). A similar description was performed considering the relativistic shock breakout material (e.g., see Fraija et al. 2019c). On the other hand, Figure 9 shows that disk wind ejecta peaks at time scales as longer as ≈ 10 5 s (see eq. A3). Therefore, we only consider the sub-relativistic decelerated material with the typical parameters of the dynamical ejecta which peaks at timescales of years (see Eq. A2), as shown in Figure 9.
We use the X-ray, optical and radio observations of GRB 170817A displayed in Fraija et al. (2021a) and Hajela et al. (2021), together with the best-fit curve found by the offaxis jet with cocoon model, which is introduced in Fraija et al. (2019b). To describe the latest multi-band observations through the sub-relativistic decelerated material in a constant-density medium, we constrain the parameter space that reproduces their synchrotron light curves. Figures 11  and 12 show the parameter space allowed for β = 0.3 and 0.4, respectively. We consider a spinning magnetized NS as the remnant of the merger of two NSs which is continuously injecting energy into the blastwave due to magnetic spin-down. For a typical value of efficiency η = 0.1 (Xiao & Dai 2019) and the half-opening angle of the KN AT2017gfo associated with GW170817 (θ j = 30 • ; Dhawan et al. 2020), the luminosity injected to the afterglow becomes similar to the spin-down luminosity. Both figures are displayed as a function of the magnetic field of the NS remnant, the isotropic-equivalent kinetic energyẼ and the constant-density medium n = A 0 that describe the latest multi-band observations for values of the synchrotron afterglow model e = 10 −1 , B = 10 −3 , α = {3.0, 4.0, 5.0} and p = {2.05, 2.15} (Bauswein et al. 2013;Metzger 2017;Fong et al. 2019;Mooley et al. 2018b;Fraija et al. 2019a;Metzger 2019;Kathirgamaraju et al. 2019;Troja et al. 2020;Fraija et al. 2021a). We can see that these parameter spaces are strongly degenerate.
The three columns in Figures 11 and 12 correspond to the values α = 3.0, 4.0 and 5.0 from left to right, respectively. In a similar manner, the two rows correspond to the values p = 2.05 and 2.15 from top to bottom, respectively. Figure  11 shows that approximately the same parameter space of the magnetic field and isotropic-equivalent kinetic energy is allowed for both values of the electron energy distribution index p. Nevertheless, there are differences in the allowed values of the circumburst density, as Figure 11 shows that for p = 2.05, larger densities (n ∼ 10 −2 cm −3 ) are preferred, while for p = 2.15, densities of the order n ∼ 10 −3 cm −3 are favored. On the other hand, upon increase of α, it is shown that the 3D parameter space shrinks, which means that smaller values of α are able to give stronger constraints on the rest of the parameters. Figure 12 is the same as 11, but it considers β = 0.4. Overall, the regions of allowed energy and magnetic field parameter space are very similar to the ones from Figure 11 and the same behaviour as described in the previous paragraph can be observed. The main difference between both figures is that for β = 0.4 lower number densities are preferred, as shown by the deeper red color of the plots. For both values of velocity (β = 0.3 and 0.4), t fb = 5 × 10 9 s and 10 13.7 B 10 16.3 G, the values of p = 2.15, n < 5 × 10 −3 cm −3 and α = 5 are preferred, although the value of α = 3.0 is not discarded. Figure 13 shows the multi-band afterglow observations of GRB 170817A, the best-fit curve from a relativistic structured off-axis jet (solid lines; Fraija et al. 2019b) and several allowed curves (dotted, dashed and dash-dotted lines) from the sub-relativistic material shown in this work. The afterglow observations are shown at X-rays, optical bands and radio wavelengths, and the synchrotron light curves are obtained at 1 keV (blue), 2.1 eV (red), 6 GHz (black) and 3 GHz (gray). In each panel we consider the PL indexes α = 3.0 (dotted line), 4.0 (dashed line) and 5.0 (dash-dotted line) for β = 0.3 (panels above) and β = 0.4 (panels below) with p = 2.05 (left) p = 2.15 (right). We can see that for different sets of parameters; we can obtain similar results about the description of the latest observations. It indicates, as expected, that our results are not unique but are possible solutions because the synchrotron equations are degenerate in these parameters.
It is relevant to mention that the value of the mean opacity for which the KN ejecta is transparent k s ≈ β 2 −1 M −1 ej,−1.3 t 2 th,3 ≈ 10 3.5 g −1 cm 2 agrees with the radiation transfer simulations, which is 10 g −1 cm −2 for lanthaniderich ejecta (Barnes & Kasen 2013;Shibata & Hotokezaka 2019). The luminosity would be estimated as L bol ≈ 10 38 erg s −1 t 1.3 th,3 M ej,−1.3 (Shibata & Hotokezaka 2019), where t th is the timescale when the KN eject enters the thin regime. It is worth noting that the features of the synchrotron emission of the sub-relativistic material is consistent with the faster "blue" kilonova afterglow.
We have considered the spinning magnetized NS scenario and discarded the fall-back material on BH scenario proposed in subsection 3.2 because the rate of fall-back accretion estimated at early times with the parameters used to describe the latest observations are fully different from the rate of fall-back accretion used in hydrodynamical simulations (e.g., see MacFadyen et al. 2001;Zhang et al. 2008). It can be demonstrated as follows. At t = t p = 1000 days, the expectation accretion rate and the BZ luminosity (Eq. 20) areṀ BH ≈Ṁ p (t/t p ) −5/3 with τ vis ≈ t p and L BZ ≈ 10 39 erg s −1 G(a)Ṁ p,−2 (t/1000 days) −5/3 , respectively. The previous derivation is similar to that resulted found by Hajela et al. (2021) after modelling the latest observations at 1000 days. However, extrapolating the rate of fall-back accretion at early times ∼ 1 s, it would bė M fb ≈ 10 −17 M s −1 , which is different from the rate of fall-back accretion used in hydrodynamical simulations (e.g., see MacFadyen et al. 2001;Zhang et al. 2008). It is worth noting that in the BH scenario the rate of fall-back accretion at early times is ∝ (t/t p ) 1/2 instead of (t/t p ) −5/3 , as considered by Hajela et al. (2021).

Short GRBs with evidence of a KN
Candidates discussed in the literature with evidence of a KN emission are GRB 050709 (Jin et al. 2016), GRB 060614 (Yang et al. 2015), GRB 130603B (Tanvir et al. 2013;Berger et al. 2013) andGRB 160821B (Kasliwal et al. 2017a;Troja et al. 2019). As follows, we present the four claimed KN observations, and then we show the synchrotron light curves with a set of allowed and ruled out parameters, assuming the characteristics of the "blue" KN.

Multi-band observations
GRB 050709 -GRB 050709 was detected on 2005 July 9 at 22:36:37 UT by the Soft X-Ray Camera (SXC), the Wide-Field X-Ray Monitor (WXM) and the French Gamma Telescope (FREGATE) instruments on board the High Energy Transient Explorer 2 satellite (HETE) with a reported location of R.A. = +23 h 01 m 30 s , Dec = −38 • 58 33 (J2000) (Villasenor et al. 2005). The prompt emission had an approximate duration of 0.5 seconds in the form of a hard spike in the 3-400 keV energy band, which was followed by an extended X-ray emission lasting ∼ 130 seconds (Jin et al. 2016). The accurate location of the burst led to the first-ever identification of the optical afterglow of a short-hard burst in ground-based experiments and HST (Hjorth et al. 2005;Fox et al. 2005). This, in turn, led to the determination of its host galaxy, which lied at redshift z = 0.16.  (Gehrels et al. 2006). However, subsequent observations showed that the event lacked an associated supernova, which is expected for lGRBs and that its temporal lag and peak luminosity were in line with those of shortduration GRBs (Jin et al. 2016). This contrasting behaviour led to denoting GRB 060614 as a 'hybrid GRB'.  (Troja et al. 2019), which places it in the sGRB class. Upon analysis of the X-ray light curves, Troja et al. (2019) found that there is evidence for continued energy injection from a long-lived central engine. On the other hand, the optical and near IR observations showed behaviour consistent with a KN. Figure 14 presents four columns, where each one corresponds to a different short GRB with evidence of KN emission. Each panel shows the multi-band afterglow observations of bursts with evidence of a KN emission and the synchrotron light curves from the cocoon (upper) with β = 0.3 and the shock breakout (lower) with β = 0.8 decelerating in a constant-density medium with n = 1 cm −3 (dashed lines) and 10 −2 cm −3 (dotted lines). The synchrotron light curves are presented at 1 keV (blue), R-band (gold) and 5 GHz (green). The disk wind and dynamical ejecta are not displayed because these decelerated materials peak at timescales longer than ≈ 10 3 days, and the multi-band afterglow observations with the respective upper limits are reported at timescales from days to ∼ one month. We consider a spinning magnetized NS with accretion as remnant of merger of two NSs. The parameter values used are P 0 = 10 −3 s, B = 7 × 10 14 G, t fb = 5 × 10 5 s,Ẽ = 2.1 × 10 49 erg, e = 0.3, B = 0.1, p = 2.05 and α = 3.0. For GRB 050709, the synchrotron emission at the F814W band from the shock breakout material with a velocity β = 0.8 is ruled out for both n = 1 cm −3 and n = 10 −2 cm −3 , but not for the cocoon material. For GRB 060614, the synchrotron emission at the R-band from the shock breakout material is ruled out for both n = 1 cm −3 and n = 10 −2 cm −3 , but not for the cocoon material. For GRB 130603B, the synchrotron curves from the shock breakout and cocoon are allowed at all bands. For GRB 160821B, the synchrotron curves at 5 GHz and at the R-band from the cocoon material are ruled out for a density of n = 1 cm −3 , but not for n = 10 −2 cm −3 . The synchrotron curves at 1 keV from the shock breakout and cocoon are allowed are allowed. Similarly, all synchrotron curves from the cocoon are allowed. The value of the uniform-density medium with n = 1 cm −3 is ruled out in our model for GRB 050709, GRB 060614 and GRB 160821B, but nor for GRB 130603B. This result is consistent with the mean value reported for sGRBs (e.g., see Berger 2014).

Analysis and Description
Continuous energy injection by the progenitor on the afterglow could generate a long-lived reverse shock at very long timescales, and shocked-accelerated electrons in this region would radiate modifying the forward-shock light curves, as found in some GRB afterglows (e.g., see Chevalier & Fransson 2006;van der Horst et al. 2007van der Horst et al. , 2014Laskar et al. 2018). In a forthcoming manuscript, we will present a detailed analysis of this long-lived reverse shock scenario.

SUMMARY
We have extended the synchrotron model presented in Fraija et al. (2021a) and presented the dynamics of deceleration of a sub-relativistic material when the central engine (a remnant of either a spinning magnetized NS or a fall-back material onto BH) injects energy into the blastwave and the external medium is stratified with a density profile A k r −k with 0 ≤ k ≤ 2.5. We have considered different profiles of the energy injection and also GRB progenitors. The energy injection index q = 0 is connected with a spinning magnetized NS and q, in general, with fall-back material onto a central black hole. The total isotropic-equivalent kinetic energy was introduced as the superposition of the energy distribution E β and the energy injection E t . When the condition E t E β is satisfied, the synchrotron light curve mimics the ones without energy injection. Otherwise, for E t E β the continuous energy injection into the afterglow dominates and many differences are observed. The constant-density medium (k = 0) is associated with the death of massive stars and the merger of two NSs, and the stratified medium (1 ≤ k ≤ 2.5) is only expected with the death of massive stars. We have presented the synchrotron light curves in radio at 1.6 GHz, optical at the R-band and X-rays at 1 keV with typical values of GRB afterglows and the scenarios of spinning magnetized NS and fall-back material onto BH with different characteristic timescales for a generic remnant located at d z = 100 Mpc. The synchrotron light curves exhibit the maximum flux on timescales from days to years, although if the remnant injects large amounts of energy, the maximum flux could be expected on timescales of hours. These light curves exhibit that, during the early stages, there is an increase in the flux. During this epoch, when the subrelativistic ejecta decelerates in a constant-density medium, the light curve grows steeply. For more stratified media, this growth is not as evident. For instance, in the light curves with a density profile with k = 2.5, the early-time behaviour is more gradual. This result implies that if such a flattening or rebrightening at timescales from months to years in the light curve was observed together with GW detection, then this would be associated with the deceleration of a sub-relativistic material launched during the merger of two NSs. Otherwise, we showed that an observed flux that gradually decreases on timescales from months to years could be associated with the deceleration of a sub-relativistic material launched during the death of a massive star with different mass-loss evolution at the end of its life. We have shown that variations of the stratification parameter are more apparent in the radio light curve when compared to the other fluxes in the other energy bands. This change in the density profile is also more easily appreciated for large values of α, (i.e., α = 4.0 and 5.0). Therefore, a transition between density profiles will be more easily observed in the radio band with high values of α.
The two NS merger ejects sub-relativistic material with distinct velocities which are decelerated by the external medium, generating the synchrotron light curves at different frequencies peaking at timescales from days to years. It is important to mention that before the ejected materials are in the sub-relativistic regime, these might begin momentarily in the TR regime (i.e., with β ∼ 0.8). Similarly, they could finally with the passage of time (∼ 10 3 years) lie in the DN regime with a velocity (i.e., β ∼ 0.08). Therefore, we have introduced timescales in both regime. For instance, with the values of parameters used we have shown that the dynamical ejecta, the cocoon, the shock breakout and the wind peak at 7150 days, 522.5 days, 56.1 days and 321.8 years, respectively. The total contribution of synchrotron light curves exhibits a brightening once the off-axis emission decreases at a few years, so that, depending on the parameter values and conditions, the synchrotron emission from the sub-relativistic materials could be detected or not. The synchrotron light curves could be associated with GW detections. We showed that, in the case of a failed or an offaxis GRB, the non-thermal emission generated by the deceleration of sub-relativistic materials could be detected at early times. In the case of an on-axis GRB, the afterglow emission originated from deceleration of the relativistic jet would have to decrease substantially so that the afterglow emission from the sub-relativistic ejecta could be observed. In addition, we gave an important tool to distinguish the afterglow emission among the sub-relativistic materials from the relativistic jet through the evolution of the synchrotron flux.
We have applied our model to describe the latest multi-wavelength afterglow observations (> 900 days) of GW170817/GRB 170817A, and constraints on the afterglow emission of some short-bursts with evidence of KNe. Regarding to the GW170817/GRB 170817A event, we have constrained the parameter space that reproduces the synchrotron light curves evolving in a constant-density medium from the faster "blue" KN afterglow. We have considered a spinning magnetized NS as remnant of the merger of two NSs which is continuously injecting energy due to spin-down luminosity into the blastwave. We plot the parameter space that describes the X-ray observations and is below the radio upper limits as a function of e = 10 −1 , B = 10 −3 , α = {3.0, 4.0, 5.0} and p = {2.05, 2.15}. We have shown that for both values of velocity (β = 0.3 and 0.4), t fb = 5×10 9 s and 10 13.7 B 10 16.3 G, the values of p = 2.15, n < 5 × 10 −3 cm −3 and α = 5 are preferred, although the value of α = 3.0 cannot be discarded. The value allowed of α = 3.0 in our theoretical model agrees with the values found in numerical simulations (e.g., see Bauswein et al. 2013) and used for describing the KN emission (Metzger 2017(Metzger , 2019, and α = 5.0 is more consistent with those reported in Kathirgamaraju et al. (2019); Hajela et al. (2021). It is worth noting that the features of the synchrotron emission of the sub-relativistic material is consistent with the faster "blue" KN afterglow.

ACKNOWLEDGMENTS
We would like to mention our deeply gratitude to the anonymous referee for his or her careful reading of the paper and useful recommendations that helped improve the clarity of the manuscript. We thank Tanmoy Laskar, Paz Beniamini, Bin-bin Zhang and Bing Zhang for useful discussions. NF acknowledges financial support from UNAM-DGAPA-PAPIIT through grant IN106521. RLB acknowledges support from CONACyT postdoctoral fellowships and the support from the DGAPA/UNAM IG100820 and IN105921. Due to the fact that materials ejected from the merger of two NSs and the gravitational collapse have an ample range of velocities, we consider in addition for this subsection the trans-relativistic (TR) and the Deep Newtonian (DN) regimes.
Trans-relativistic regime -In this regime, the kinetic energy of the shock in a stratified medium is given by E = 4π 3 σm p β 2 Γ 2 A k r 3−k (Blandford & McKee 1976) where the parameter σ is a function of velocity (σ = 0.73 − 0.38β; Huang et al. 1998). The deceleration time becomes t dec,TR 23.1 days 1 + z 1.022 When the bulk Lorentz factor Γ → 1, the TR timescale approaches the Newtonian timescale.

A.2. Synchrotron emission
During the deceleration phase, the post-shock magnetic field evolves as, B ∝ t − 2(2+q)+k(1−q+α) 2(α+5−k) . The Lorentz factors of the lowest-energy electrons and of the higher energy electrons, which are efficiently cooled by synchrotron emission are The corresponding synchrotron break frequencies are given by In the self-absorption regime, the synchrotron break frequencies are νa,1 = ν 0 a,1 The spectral peak flux density becomes The quantities γ 0 m , γ 0 c , ν 0 m , ν 0 c and F 0 ν,max given in eqs. A4, A5, A6 and A7 are reported in Table 1 for k=0, 1, 1.5, 2 and 2.5.

Flux density (mJy)
Time since burst (day) Jet (θ obs = 30º) Figure 9. Synchrotron light curves generated by the deceleration of the relativistic jet, cocoon, the shock breakout, the dynamical, and the wind ejecta in the constant-density medium (k = 0). The left-hand panels shows the light curves from an on-axis jet and the right-hand ones from an off-axis jet with a viewing angle of 30 • . The spinning magnetized NS scenario with the parameter values P0 = 10 −3 s, B = 10 16 G and 10 -1 10 0 10 1 10 2 10 3