The electromagnetic emission and the afterglow observations of the binary neutron star merger event GW170817A confirmed the association of the merger with a short gammaray burst GRB harboring a narrow 510 and powerful 10SUP49SUP10SUP50SUP erg jet. Using the 1 s long neutrinoradiation general relativistic MHD simulation of coalescing neutron stars of K. Kiuchi et al. and following the semianalytical estimates of M. Pais et al. we inject a narrow powerful unmagnetized jet into the postmerger phase. We explore different opening angles luminosities central engine durations and times after the merger. We explore early 0.1 s following the merger and late 1 s jet launches the latter is consistent with the time delay of 1.74 s observed between GW170817 and GRB 170817A. We demonstrate that the semianalytical estimates correctly predict the jets breakout and collimation conditions. When comparing our synthetic afterglow light curves to the observed radio data of GW170807 we find a good agreement for a 3  10SUP49SUP erg jet launched late with an opening angle in the range 57.
['Gamma-ray bursts', 'Neutron stars', 'Relativistic jets', 'Hydrodynamical simulations', '629', '1108', '1390', '767', 'Astrophysics - High Energy Astrophysical Phenomena']
Simulating Short GammaRay Burst Jets in Realistic Late Binary Neutron Star Merger Environments
{'Simulating short GRB jets in realistic late binary neutron star merger environments': "M/a.pc/t.pc/t.pc/e.pc/o.pc P/a.pc/i.pc/s.pc , 1, 2 T/s.pc/v.pc/i.pc P/i.pc/r.pc/a.pc/n.pc , 1 K/e.pc/n.pc/t.pc/a.pc K/i.pc/u.pc/c.pc/h.pc/i.pc ( 木 内 建 太 ) , 3, 4 /a.pc/n.pc/d.pc M/a.pc/s.pc/a.pc/r.pc/u.pc S/h.pc/i.pc/b.pc/a.pc/t.pc/a.pc ( 柴 田 大 ) 3, 4  \n1 Racah Institute for Physics, The Hebrew University, Jerusalem 91904, Israel  \n2 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy 3 Max Planck Institute for Gravitational Physics (Albert Einstein Institute) Am Mühlenberg 1, 14476, Potsdam-Golm, Germany 4 Center for Gravitational Physics and Quantum Information,  \nYukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan", 'ABSTRACT': "The electromagnetic emission and the afterglow observations of the binary neutron star merger event GW 170817A confirmed the association of the merger with a short gamma-ray burst (sGRB) harboring a narrow (5-10 · ) and powerful (10 49 -10 50 erg) jet. Using the 1 second-long neutrino-radiation-GR-MHD simulation of coalescing neutron stars of Kiuchi et al. (2023), and following the semi-analytical estimates of Pais et al. (2023) we inject a narrow, powerful unmagnetized jet into the post-merger phase. We explore different opening angles, luminosities, central engine durations, and times after the merger. We explore early (0 . 1 s following the merger) and late (1 s) jet launches; the latter is consistent with the time delay of ≈ 1 . 74 s observed between GW 170817 and GRB 170817A. We demonstrate that the semi-analytical estimates correctly predict the jets' breakout and collimation conditions. When comparing our synthetic afterglow light curves to the observed radio data of GW170807, we find a good agreement for a 3 × 10 49 ers jet launched late with an opening angle in the range ≃ 5 · -7 · .  \nKeywords: stars: jets - gamma-ray burst: general - stars: neutron - hydrodynamics", '1. INTRODUCTION': "The merger of a neutron star binary (BNS) or a neutron star and a black hole (BH) is the most promising candidate for the detection of gravitational waves and their electromagnetic counterparts: a short GRB (Eichler et al. 1989), a kilonova (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010) and their afterglows. The tentative observations of kilonovae following sGRB 130603B (Tanvir et al. 2013; Berger et al. 2013) and other sGRB (e.g., Yang et al. 2015; Jin et al. 2016, 2015; Gompertz et al. 2018; Ascenzi et al. 2019; Jin et al. 2020; Lamb et al. 2019a; Rossi et al. 2020; Fong et al. 2021) suggested such association but these were not a sufficient affirmation. The BNS merger GW 170817 provided the first evidence of the detection of a gravitational wave from a BNS merger (Abbott et al. 2017a). It was followed by a multiwavelength electromagnetic signal (Abbott et al. 2017b).  \nNumerous numerical simulations of mergers have been carried out over the years (e.g., Davies et al. 1994; Ruffert et al. 1995; Rosswog et al. 1998; Shibata 1999; Shibata & Ury¯u 2000; Rosswog et al. 2000; Oechslin et al. 2002; Shibata et al. 2005; Shibata & Taniguchi 2006; Anderson et al. 2008; Liu et al. 2008; Giacomazzo et al. 2009, 2011; Giacomazzo & Perna 2013; Piran et al. 2013; Bauswein et al.  \n2013; Rosswog et al. 2013; Hotokezaka et al. 2013; Kiuchi et al. 2014; Dionysopoulou et al. 2015; Giacomazzo et al. 2015; Sekiguchi et al. 2015, 2016; Foucart et al. 2016; Kawamura et al. 2016; Radice et al. 2016; Ruiz et al. 2016; Ciolfi et al. 2017; Dietrich et al. 2017; Ruiz & Shapiro 2017; Kiuchi et al. 2018; Ruiz et al. 2018; Ciolfi et al. 2019; Ruiz et al. 2019; Aguilera-Miret et al. 2020; Ciolfi 2020; Mösta et al. 2020; Ruiz et al. 2020; Hayashi et al. 2022; Ruiz et al. 2021; Aguilera-Miret et al. 2022; Palenzuela et al. 2022; Sun et al. 2022; Aguilera-Miret et al. 2023; Combi & Siegel 2023a,b; Most & Quataert 2023; Kiuchi et al. 2024). State-of-theart simulations rely on fully general relativistic numerical schemes, approximated neutrino transport, and a realistic equation of state for the neutron stars. Before coalescence, the system starts to eject mass because of the tidal forces. Typically, the resulting ejecta expands at 0 . 1 c (Hotokezaka et al. 2013), with a small fraction that moves faster, around 0 . 6 c or even faster (Kiuchi et al. 2017; Fujibayashi et al. 2023). A second source of ejecta is the bound material torn from the NS after the merger with enough angular momentum to form a thick torus around the newly created compact object from the fusion. Magneto-rotational instability (MRI) generates viscosity, leading to accretion (Zenati et al. 2023, 2024).  \nA viscous-driven wind emerges from the disk and carries a fraction of the disk material at low ( < 0 . 1 c) velocities (e.g., Fernández & Metzger 2013; Just et al. 2015; Fujibayashi et al. 2020).  \nThe afterglow observations of GW 170807 revealed that it involved a relativistic jet launched at an angle of about 18 · awayfromus, which means that while we did not see the direct emission from this jet, as observed in a typical GRB, a powerful relativistic jet was present in this system (e.g., Troja et al. 2017; Hallinan et al. 2017; Mooley et al. 2018a,b; Ghirlanda et al. 2019). Close examination of the 𝛾 -ray observations and the afterglow shows that GRB 170817A was not a regular sGRB. Its total isotropic equivalent energy ( ∼ 10 46 erg) is smaller by three orders of magnitude than the weakest sGRB measured so far (Gottlieb et al. 2018) and by four orders of magnitude than typical sGRBs (Nakar 2007). Both the radio (Hallinan et al. 2017) and X-ray (Troja et al. 2017) afterglows were delayed by several days compared to a regular sGRB. The compactness argument (Kasliwal et al. 2017; Matsumoto et al. 2019) reveals that the observed 𝛾 -rays must have been produced in a mildly or fully relativistic outflow with a Lorenz factor of Γ ≃ 2 -3.  \nAs the jet propagates through the expanding ejecta, it forms a high-pressure bubble (the cocoon) that engulfs the jet (e.g., Bromberg et al. 2011). A natural mechanism that could have produced the observed 𝛾 -rays in GRB 170817A is a shock breakout of the cocoon from the engulfing ejecta (Nakar & Piran 2017; Kasliwal et al. 2017; Gottlieb et al. 2018; Beloborodov et al. 2020), while the interaction of the jet and cocoon with the circumburst medium produced the afterglow.  \nThe propagation of the jet in a surrounding medium has been studied by numerous authors both analytically (e.g., Blandford & Rees 1974; Begelman & Cioffi 1989; Mészáros &Waxman 2001; Matzner 2003; Lazzati & Begelman 2005; Bromberg et al. 2011; Hamidani et al. 2019; Hamidani & Ioka 2021, 2022, 2023; García-García et al. 2024) and numerically (e.g., Martí et al. 1995, 1997; Aloy et al. 2000; MacFadyen et al. 2001; Zhang et al. 2004; Mizuta et al. 2006; Morsony et al. 2007; Wang et al. 2008; Lazzati et al. 2009; Mizuta & Aloy 2009; Morsony et al. 2010; Nagakura et al. 2011; LópezCámara et al. 2013; Ito et al. 2015; López-Cámara et al. 2016; Harrison et al. 2018; Gottlieb et al. 2018; Pavan et al. 2021; Gottlieb & Nakar 2022; García-García et al. 2023; Pavan et al. 2023). Most of those studies were conducted for jets expanding in a uniform or spherically symmetric surrounding medium (but see Nagakura et al. (2014)). Only recently, inhomogeneous post-merger environments directly imported from the outcome of a GR-MHD BNS merger simulations were considered (Nativi et al. 2020; Pavan et al. 2021; Nativi et al. 2022; Pavan et al. 2023). These studies describe how a jet dissipates energy as it propagates, enabling us to determine whether it breaks out of the surrounding matter or  \nis choked inside. They also show whether the surrounding pressure is sufficient to collimate the jet (Komissarov et al. 2009; Lyubarsky 2009).  \nWhile uniform or spherically symmetric external matter is a good approximation for stellar envelopes (relevant for long GRBs), the ejecta of the BNS merger are far from that. Pais et al. (2023) explored semi-analytically the conditions for a jet to break out from the ejecta structure that arose in a BNS-post merger environment based on the latest simulation described in Kiuchi et al. (2023). A key ingredient here is the interplay between the external pressure that acts to collimate the jet and the existence of an empty funnel along the system's rotation axis, along which the jet propagates almost freely. Pais et al. (2023) found that a jet with energy between 10 48 -10 50 erg and an opening angle of 3 · -7 · can break out from this ejecta. At the same time, they found that the breakout conditions are almost independent of the time at which they are calculated. In this work, we numerically test the validity of these results by injecting an unmagnetized relativistic jet in the same BNS merger environment, following its evolution numerically.  \nThe paper is structured as follows. We begin, in Section 2, with a brief description of the environment setup resulting from the neutrino-radiation-GRMHD simulation of a BNS merger, referring the reader to Kiuchi et al. (2023) for details. The following paragraphs briefly describe the data import procedure and the jet injection into the post-BNS environment. More details on the data structures are reported in Appendix A. In Section 3, we describe the results of several realizations of the jet with different total energy, exploring cases of jet breakout and jet choking. We calculate the afterglow light curves using the energy distribution from our simulations. We discuss the implications of our findings for GW 170817 and other sGRBs and summarise our results in Section 4. Throughout this paper, 𝑐 denotes the speed of light.", '2.1. The BNS environment data': 'We import our initial conditions describing the ejecta surrounding the remnant BH into /p.pc/l.pc/u.pc/t.pc/o.pc from the fully 3D general relativistic neutrino-radiation MHD simulation of a BNS merger described in greater detail in Kiuchi et al. (2022) and Kiuchi et al. (2023). In these simulations, the neutron star is modeled with the SFHo equation of state (Steiner et al. 2013), and the binary is composed of a 1 . 2 M ⊙ and a 1 . 5 M ⊙ neutron star, with a chirp mass consistent with the one observed in GW 170817. We employ the /l.pc/o.pc/r.pc/e.pc/n.pc/e.pc library (Gourgoulhon et al. 2016) to construct a quasi-equilibrium configuration of the irrotational BNS with the initial orbital frequency of 𝐺𝑀 0 Ω 0 / 𝑐 3 = 0 . 025, where 𝑀 0 = 2 . 7 M ⊙ is the total mass, and with an orbital eccentricity of O( 10 -3 ) . The entire run lasts one full second and is the longest simulation of a BNS  \n])  \n])  \nFigure 1. 2D slice across the 𝑥 -𝑧 plane passing through the center of the simulation box of the initial conditions taken at 𝑡 = 1 . 06 s after the merger. Left : density slice in logarithmic units superposed to the streamlines of the velocity field. Right : pressure slice in logarithmic units. In both plots, we can notice the presence of the low-pressure/low-density polar funnel.  \n<!-- image -->  \n<!-- image -->  \n]  \nm  \nk  \n[  \nz  \nmerger in the literature. The merger and the subsequent BH formation occur at ≈ 0 . 015 s and ≈ 0 . 032 s, respectively, where 𝑡 = 0 is the start time of the simulation.  \nA fraction of neutron star matter is dynamically ejected by the tidal interaction and shock heating in the merge phase and shows a fast tail with a terminal velocity up to ≃ 0 . 96 c (Fujibayashi et al. 2023). A high-mass neutron star (HMNS) is formed and survives for ∼ 20 ms after the merger before collapsing to a BH. It has a non-axisymmetric density structure and exerts a gravitational torque on the plasma, transporting angular momentum outwards. The torque leads to the formation of a torus of mass ≃ 0 . 05 M ⊙ , which surrounds a newly formed spinning BH of mass ≃ 2 . 55 M ⊙ and a dimensionless spin of ≃ 0 . 65.  \nThe poloidal magnetic field is initialized with a magnetarlike strength of 10 15 G, and the Kelvin-Helmholtz instability at the onset of the merger amplifies it further (Price & Rosswog 2006; Kiuchi et al. 2014, 2015, 2018; Aguilera-Miret et al. 2020, 2022; Palenzuela et al. 2022; Aguilera-Miret et al. 2023, 2024; Kiuchi et al. 2024). Non-axis-symmetric dynamo-driven MRI contributes to the enhancement of the magnetic field. Still, even at this level, the magnetic field is irrelevant to the overall dynamical evolution. The magnetic winding and radial motion of the fluid enhance the magnetic field in the torus until the electromagnetic energy saturates at ∼ 1% of the internal energy of the torus. At 𝑡 ≈ 0 . 04 s after the merger, the magnetic field energy is maximal and then decreases. At all stages, the magnetic field energy is subdominant to the total energy budget ( ∼ 10 49 erg vs ∼ 10 51 erg, see Kiuchi et al. (2023)). The subdominant magnetic field energy allows us to neglect the magnetic fields in our work, which focuses on jet propagation.  \nAt 𝑡 = 0 . 1 s, a thick torus of roughly 300 km in radius surrounds the compact object, and a funnel structure has already formed with its distinctive features of low density and low pressure. The torus expands outward due to MRI-driven turbulent viscosity transporting angular momentum, and as it expands, it cools adiabatically. A part of the torus matter is also ejected due to this MRI-driven turbulence as the post-merger ejecta after the neutrino cooling in the torus becomes inefficient. This occurs at 𝑡 ≈ 0 . 3 s after the merger. At around 𝑡 ≈ 1 s, the post-merger ejecta has a typical velocity of around 0.15 c, slower than the average value of 0.2 c reported in literature (see Sekiguchi et al. 2016). This is largely due to the post-merger ejecta that is launched at the outskirt of the massive torus formed around the black hole. This turbulent viscosity-driven post-merger ejection still lasts ≃ 1 s.  \nAt 𝑡 = 1 s the thick torus surrounding the BH has a size of ≃ 750 km of radius and a comparable thickness, with an average density of 𝜌 torus > 10 6 g / cm 3 . A funnel with a half-opening angle of roughly 10 · and an average density below 3 × 10 3 g / cm 3 exists along the polar axis. As the simulation progresses in time, the entire system expands in size. The total mass of the ejecta at 𝑡 = 1 s is approximately 1 . 5 × 10 -2 M ⊙ . Fig. 1 shows the initial conditions of the simulation at 𝑡 = 1 . 06safter the merger, in which the presence of two polar funnels with low density and low pressure is noticeable.  \nThe post-merger secular ejecta expand in a \'numerical vacuum" of circumbinary medium with a uniform pressure of 1 . 35 × 10 12 erg / cm 3 (in order to avoid spurious random motion) and a density of 𝜌 ≃ 0 . 17 g / cm 3 . This results in a mass surrounding the merger of approximately 1 . 2 × 10 -3 M ⊙ enclosed in a radius of 𝑅 max = 1 . 5 × 10 5 km, which is an order  \nof magnitude less than the mass of the ejecta. The fast tail of the dynamical ejecta stretching to 0.96 c gradually slows down due to interaction with the external density that occupies the "numerical vacuum\' (see Kiuchi et al. (2023) and Fujibayashi et al. (2023) for more details) leaving behind only slower, mildly-relativistic material. Apart from influencing the low-mass, the very high-velocity tail of the unbound material, this artificial vacuum doesn\'t significantly impact the ejecta\'s expansion within this simulation box. Hence, it does not influence our results.', '2.2. Computational setup': 'We construct our 2.5D initial conditions from the BNS simulation by cylindrically averaging the original 3D data (details of importing the BNS simulation data into our calculations, using /p.pc/l.pc/u.pc/t.pc/o.pc, are described in Appendix A). The original cylindrical average of the output is imported and regridded into /p.pc/l.pc/u.pc/t.pc/o.pc (details of the computational grid can be found in Appendix A.2). After the data import, we evolve the matter distribution within the code and, with a user-defined boundary prescription, we insert a jet into the environment through a nozzle of a radius 𝑟 j at the excision height 𝑧 0, imposing a reflective boundary elsewhere along 𝑧 = 𝑧 0 (and an axis-symmetric boundary along 𝑟 = 0). Our choice of the boundary condition at 𝑧 0 does not influence the outcome of our simulations since the velocity of our injected outflow is much larger than the velocity of the ejecta. In addition, due to the strong power of the jet, the newly formed cocoon quickly dominates the ejecta structure in the innermost part of the ejecta.  \nWe perform our simulations using the massively parallel multidimensional relativistic magneto-hydrodynamic code /p.pc/l.pc/u.pc/t.pc/o.pc (v4.4.2) /one.sup (Mignone et al. 2007). The code uses a finite-volume, shock-capturing scheme that integrates a system of conservation laws. The flow quantities are discretized on a rectangular computational grid enclosed by a boundary. We use the special relativistic hydrodynamics module in 2.5D cylindrical coordinates (i.e., a 2D grid in cylindrical geometry in which 𝑣 𝜙 is also integrated and evolved with periodic boundary conditions). We performed our simulations using a linear reconstruction scheme combined with a second-order Runge-Kutta time stepping. We also force the code to reconstruct the 4-velocity vectors at each time step. We chose a Taub-Matthews equation of state (Mignone & McKinney 2007), which accounts for both relativistic and non-relativistic material, i.e., 𝛾 = 5 / 3 for non-relativistic material and 𝛾 = 4 / 3 for relativistically moving matter.  \nThe choice of a 2.5D simulation over 3D has implications for both the jet propagation and the development of hydrodynamic turbulence in the flow. A direct effect of the dimen-  \nsionality is the direction of the turbulent cascade. While in 3D, the cascade goes from large to small scales, in 2D, the opposite is true, as shown by Kraichnan (1967), resulting in a flow organizing itself in large eddies and generating different turbulent structures that affect the mixing of the jet material with the ejecta. Since a 2.5D simulation lacks a degree of freedom, the jet is forced to propagate, staying on-axis. This results in a numerical effect called \'plug", which consists of matter stalling over the jet and affecting its expansion. Conversely, in 3D simulations, this feature is absent, and the jet, while carving its way out, wobbles due to flow instabilities and interaction with the surrounding ejecta. When 2D and 3D runs are compared, the difference is not too stark (Harrison et al. 2018).', '2.3. Jet injection': "We inject a jet with constant luminosity 𝐿 j, operating for an engine time 𝑡 e with an initial bulk Lorentz factor Γ 0 , j, a density 𝜌 j and a dimensionless specific enthalpy ℎ j. The hot jet is injected through a nozzle parallel to the 𝑧 -axis with a radius of 𝑟 j and with an initial opening angle, 𝜃 j, which also determines the initial Lorentz factor via Γ 0 , j ≃ 1 /( 1 . 4 𝜃 j ) (see details on this injection method at Mizuta & Ioka 2013 and Harrison et al. 2018). We place the nozzle at an initial excision height of 𝑧 0 = 𝑟 j / 𝜃 j. In all our simulations, we chose 𝑟 j = 10 km, a size comparable to the Schwarzschild radius of the central BH, allowing sufficient mesh coverage. The head cross section is given by Σ j = 𝜋𝑟 2 j . This determines 𝜌 j as:  \n𝜌 j = 𝐿 j Σ j ( ℎ j Γ 2 0 , j -ℎ eff ) 𝛽 0 , jc 3 , (1)  \nwith 𝛽 0 , j = √︃ 1 -1 / Γ 2 0 , j and ℎ eff = GLYPH<16> 5 ℎ j -√︃ 9 ℎ 2 j + 16 GLYPH<17> / 8, which is derived from Taub-Matthews EOS. We choose ℎ j = 100, which gives ℎ eff ≃ 25. This choice of the specific dimensionless enthalpy determines the terminal Lorentz factor of the system, i.e. Γ ∞ ≃ ℎ j Γ 0 , j. The jet's pressure is given by 𝑃 j = ℎ eff 𝜌 jc 2 . For details regarding the computational grid used in /p.pc/l.pc/u.pc/t.pc/o.pc, we refer to Appendix A.2. At the same time, a more detailed description of the jet injection is found in Appendix B.1.", '3.1. Collimation criterion': 'In Pais et al. (2023), we showed how the ambient pressure counterbalances the Poynting flux of a heavily magnetized jet, assuming transverse equilibrium for the magnetic field. In our case, the jet is not magnetized, and an equation of state with 𝛾 = 4 / 3 describes the relativistic pressure injected. So, we cannot strictly apply the same criteria for collimation as described in our previous work. In particular, if 𝑃 j ≃  \nFigure 2. Jet breakout conditions for injected jets with 𝑡 d = 0 . 1 s and 𝑡 d = 1 s after the merger. Top : 𝐿 j = 10 50 erg / s; bottom : 𝐿 j = 10 51 erg / s. The colored lines represent the energy required to collimate the jet to a certain opening angle according to Eq. (3). In contrast, the black solid lines represent the energy required for a jet to escape as a function of the opening angle 𝜃 j derived from Eq. (4). The dots represent the tested cases for 𝜃 j = 5 · and 10 · with triangles (stars) representing successful (choked) jets. The black dashed lines in each panel, superposed to the colored lines, represent the fit for the collimation curve expressed by Eq. (3).  \n<!-- image -->  \nFigure 3. Minimal engine time, 𝑡 e, required for a jet to escape as a function of the injection angle 𝜃 j , the luminosity, 𝐿 j , and the time delay, 𝑡 d , between the merger and jet launch. The lines represent the estimate from Eq. (5). We can see that in the case of 𝑡 d = 0 . 1 s we get 𝑡 e ≳ 0 . 05 s, while for 𝑡 d = 1 s, this value increases to 𝑡 e ≃ 0 . 1 -0 . 3 s for jets with 𝜃 j = 5-10 · . All these values for 𝑡 e are compatible with a sGRB.  \n<!-- image -->  \n( ℎ j -1 ) 𝜌 jc 2 / 4 for hot jets and if 𝜌 j is given by Eq. (1), then the equilibrium condition with ambient pressure reads:  \n𝑃 ej ≃ 𝐿 j 4 𝜋𝑟 2 Γ 2 c , (2)  \nwhere 𝑃 ej ( 𝑅, 𝑧 ) is the local ejecta pressure, and Γ = Γ ( 𝑅 ) is the Lorentz factor of the flow.  \nFor a relativistic jet confined by a non-static medium along the cross-sectional radius 𝑅 of the jet, the scaling Γ ∝ 𝑅 holds (Bromberg & Levinson 2007), which is the same scal-  \nwhen transverse equilibrium for the magnetic field is reached for magnetically dominated jets. In Pais et al. (2023), we derived the collimation curves for the simulation data of Kiuchi et al. (2023), varying the luminosity of the jet smoothly, determining the corresponding equilibrium contours between the jet pressure and the ambient pressure, and calculating the opening angle of these contours with the assumption of a quasi-conical shape for the jet. For the increase in the luminosity, the increase in the jet opening angle', 't sim = 1.360 s , t j = 0.300 s': 'Figure 4. An example of evolved jets into the BNS-post merger environment for the model 𝑡 d = 1 s after 300 ms of continuous injection for three different cases: top -𝐿 = 10 50 erg/s and 𝜃 = 5 · , middle -𝐿 = 10 51 erg/s and 𝜃 = 5 · , bottom -𝐿 = 10 51 erg/s and 𝜃 j = 10 · . We started the simulation at 𝑡 = 1 . 06 s, which is roughly at 𝑡 d = 1 s after the collapse of the HMNS. The time 𝑡 j = 0 . 3 s indicates the time since the jet injection has started. The white line represents the opening angle of the jet. From left to right, the panels represent the Γ 𝛽 factor, the density 𝜌 , the pressure 𝑃 , and the jet tracer (flagged cells with jet material advected with the fluid).  \n<!-- image -->  \nwas obtained, giving a functional relation between these two quantities. Since we use the same background data from the BNS merger as the initial condition for our jet injection, we can fit the collimation curve 𝐿 j , coll as a function of the angle 𝜃 j with a power law ≃ 𝐿 0 × 𝜃 4 j . Thus, under the assumption that the luminosity is constant over time and given an angle 𝜃 j, the corresponding collimation energy is:  \n𝐸 coll ( 𝜃 j ) = ∫ 𝑡 e 0 𝐿 j , coll ( 𝜃 j ) d 𝑡 ≈ 3 . 6 × 10 49 GLYPH<18> 𝜃 j 0 . 1 rad GLYPH<19> 4 GLYPH<16> 𝑡 e 1 s GLYPH<17> erg . (3)  \nWeremark that this formula is an approximation, and its prefactor is tailored to the specific energy distribution of the BNS simulation we are using. Since the luminosity is assumed to be a constant, the collimation energy as a function of the jet angle scales linearly with the engine time. Analyzing the simulations, we report our fit in Fig. 2 as a black dashed line over the colored collimation curve lines.', '3.2. Breakout criterion': 'Gottlieb & Nakar (2022) estimate the breakout criterion for a jet surrounded by anisotropic, fast-moving ejecta. Their condition for a jet of luminosity 𝐿 j reads:  \n𝐸 j ≡ 𝐿 j × 𝑡 e > 150 " GLYPH<18> 𝑡 d 𝑡 e GLYPH<19> 2 + 2 # 𝐸 ej ( < 𝜃 j ) × 𝜃 2 j , (4)  \nwhere 𝑡 d is the delay time between the merger and the jet launch, 𝑡 e is the engine time of the jet, 𝜃 j is the jet opening angle, and 𝐸 ej ( < 𝜃 j ) is the energy of the ejecta within 𝜃 j. We see immediately that if 𝑡 d ≫ 𝑡 e, it becomes energetically difficult for a jet to break out from the ejecta. If we now express 𝐸 ej ( < 𝜃 ) as a function of 𝜃 , then we fit the data from the BNS simulation as 𝐸 ej ( < 𝜃 j ) ≃ 𝐸 0 × 𝜃 𝑘 j (where 𝑘 ≈ 2 ) , and if we express 𝐸 j as 𝐿 j × 𝑡 e, we get a third-order polynomial in 𝑡 e, which gives us an approximate estimate of the minimal engine time for a jet to break out as a function of 𝐿 j, 𝜃 j and 𝑡 d. This estimate, for typical jet parameters, reads:  \n𝑡 e ≃ 0 . 1 s GLYPH<16> 𝐿 -1 j , 50 𝜃 4 j , 5 + 𝐿 -1 / 3 j , 50 𝜃 4 / 3 j , 5 𝑡 2 / 3 d , 0 GLYPH<17> , (5)  \nwhere 𝐿 j , 50 = 𝐿 j /( 10 50 erg / s ) , 𝜃 j , 5 = 𝜃 j / 5 · , and 𝑡 d , 0 = 𝑡 d /( 1 s ) . Since 𝑡 e ≲ 1 s for sGRB jets and the jet is likely to be narrowly collimated, the first term of Eq. (5) can be neglected. This second term agrees with the estimates given by Eq. (31) of Nakar (2019), even though it was derived under some restrictive conditions that do not hold here.  \nWereport in Fig. 3 the estimate of the engine time required for the jet to escape from the ejecta using Eq. (5). It is also worthwhile to notice how the prefactor in Eq. (4) quickly increases when 𝑡 e ≪ 𝑡 d, so we will mainly focus on engine times that have a duration between 0 . 05 and 0 . 3 s. We take  \n𝑡 d [s] 0.1 1 𝐿 j [erg/s] 10 50 10 51 𝜃 j [deg] 5 10  \nTable 1. Parameters used in the suite of simulations presented in Sec. 3.3  \nEq. (5) as the minimal estimate for the engine time required by the jets to escape.  \nFig. 2 depicts the escape and collimation conditions for the jet, based on the simulation analysis of Pais et al. (2023). The points in the plots show simulated jets: the triangles represent jets that emerged from the ejecta, while the stars represent jets that choked inside the ejecta. We stress that these curves (from Pais et al. 2023) represent only estimates of the energy required for the jet to escape and that the actual energy might differ by a factor of a few. Moreover, Gottlieb & Nakar (2022) develop their analytic theory for a spherical configuration, and here we used an adaptation (Pais et al. 2023) of this estimate to anisotropic ejecta, given by Eq. (4). Still, as we show later, the numerical results are consistent with the analytic estimates.', '3.3. Simulation': "Since the computations are at high resolution and computationally expensive (around 150000 CPU hours), we run the following simulations, changing one parameter each time, as reported in Table 1, giving a total of 8 simulations, 4 for each 𝑡 d. To test the collimation effect of the ambient pressure on the jet, we performed two further simulations with a wide opening angle, the first with 𝜃 j = 35 · , the second with 𝜃 j = 20 · , and 𝐿 j = 10 51 erg / s, and 𝑡 d = 1 s, for both cases. Except for the case study of 𝐿 j = 10 50 erg / s, 𝜃 j = 5 · and 𝑡 d = 1 s (where we let the injection run for 𝑡 e = 550 ms), we run our jets for a maximum time of 𝑡 e = 300 ms (100 ms) in the case of 𝑡 d = 1 s (0.1 s). Although the BH formation happens quite early after the merger in our imported ICs, we explore the case of a late launch at one second to maximize the effect of the ejecta collimation over the jet propagation. We also run shorter-duration jets for 𝑡 d = 0 . 1 s since the ejecta did not expand too far from the center, and the jet is likely to break out earlier. With a continuous and constant energy injection in our system, we can stop the injection at any time and follow the evolution of the jet-cocoon system. This way, we can determine whether the jet engine time has been sufficient for the jet to break out of the ejecta.  \nFurthermore, we are also interested in the jets' evolution until their breakout from the BNS outermost ejecta. For this reason, we run several simulations where we stop the injection of the jet at different moments, i.e., 𝑡 e = [ 0 . 01 , 0 . 03 , 0 . 05 , 0 . 1 , 0 . 12 , 0 . 3 ] s. We performed a lowerresolution re-grid of the simulations (see Appendix A.2) to speed up the jet's evolution after the engine is shut down.  \n0  \nFigure 5. Proper velocity ( Γ 𝛽 ) and jet scalar tracer maps of a jet initialized with a delay of 𝑡 d = 1 s, a luminosity of 𝐿 j = 10 51 erg / s, and an opening angle 𝜃 j = 5 · at an early stage of its injection (taken after injecting for 0 . 025 s). At this stage, the jet is still in the ejecta's innermost (and densest) part of the ejecta and drilled through the polar funnel. We can see the shape of the collimation shock, which divides the unshocked jet material from the shocked one.  \n<!-- image -->", '3.3.1. Jet evolution for 𝑡 d = 1 s': "In Fig. 4, we show jet simulations with 𝑡 d = 1 s and 𝑡 e = 300 ms. Note that 𝑡 d = 1 s is approximately the time delay that took place in 170817A (Kasliwal et al. 2017; Gottlieb et al. 2018). The image shows three cases: one with a luminosity of 𝐿 j = 10 50 erg / s and 𝜃 j = 5 · ∼ 0 . 1 rad (first row), one with the same opening angle but a luminosity of 𝐿 j = 10 51 erg / s (second row), one with 𝐿 j = 10 51 erg / s and 𝜃 j = 10 · ∼ 0 . 2 rad (third row). The flow is injected at 𝑧 0 = 100 km and 𝑧 0 = 50 km for 𝜃 j = 5 · and 10 · , respectively. From those, the combination of 𝐿 j = 10 50 erg and 𝜃 j = 5 · ∼ 0 . 1 rad is very similar to the jet properties inferred in GRB 170817A from the afterglow observations (see e.g. Nakar & Piran 2020, for a summary of the afterglow analysis). Other simulations, not shown here, are qualitatively similar in their behavior. Here and in the subsequent figures, we show the proper velocity  \nFigure 6. Top : Differential energy distribution with respect to the proper velocity Γ 𝛽 of a relativistic jet with 𝐿 j = 10 50 erg / s and 𝜃 j = 5 · injected into the ejecta. The different colors show the jet time evolution during the injection (cold colors for early times, hot colors for late times). The dotted colored lines represent the material associated with the jet. The black-dashed line is the energy of the initial conditions ejecta. Note that we prolonged the simulation just in this case to see the evolution of the distribution till the breakout time 𝑡 bo ≃ 500 ms after the start of the injection. Bottom : Different regions of the energy-velocity distribution after 550 ms of injection with respect to different thresholds of the jet scalar tracer 𝑄 j .  \n<!-- image -->  \nΓ 𝛽 , the density 𝜌 , the pressure 𝑃 , and the jet tracer 𝑄 j, which is a scalar marker that flags the jet material fraction and is passively advected with the fluid, explicitly: 𝜕 𝑡 ( 𝜌𝑄 j ) + ∇ · ( 𝜌𝑄 j v ) = 0, where 𝜕 𝑡 denotes the partial derivative with respect to time.  \nThe jet starts to carve out of the ejecta, filling in ≲ 30 ms the low density-low pressure polar cavities above the compact object. As stated in Pais et al. (2023), this is a crucial step to let the jet escape: the anisotropic distribution of the ejecta, dense at the equator and dilute at the poles, allows the jetted material to quickly escape the central densest ejecta near the BH ( ≃ 10 4 km ), entering in a less dense external region, where it is easier for the jet to drill its way out up to the ejecta outer radius 𝑅 ej. As it propagates, the jet forms a cocoon, resulting from its interaction with the merger ejecta. The cocoon grows in size as the jet proceeds its vertical evolution. Alarge proper velocity characterizes the unshocked jet region. It resides in the low-density, low-pressure axial cavity and spans an angle roughly the size of the original opening angle. The white oblique lines represent the initial opening angle of the expanding jets. These lines roughly coincide with the opening angle covered by the highest proper velocity material close to the jet axis. Inspection of the early phases of the jet propagation (see Fig. 5) reveals a collimation shock. The ambient pressure quickly collimates the jet (in a few tens  \nFigure 7. Relativistic jets after we switched off the jet engine running for 𝑡 e = 100 ms (top) and 𝑡 e = 120 ms (bottom) and let the jet evolve until it reached 𝑅 ej (post-injection phase) for 𝐿 j = 10 50 erg / s and 𝜃 j = 5 · (top row), and 𝐿 = 10 51 erg/s and 𝜃 j = 10 · (bottom row). The jet's high-velocity material reaches the ejecta's outer boundary after the jet engine shuts down, leaving jetted material that mixes with the BNS ejecta matter (rightmost panels). In both cases, the breakout happens at 𝑡 j ≃ 0 . 5 s.  \n<!-- image -->  \nof milliseconds). The collimated jet maintains roughly its original opening angle.  \nIn the top panel of Fig. 6, we show the time evolution of the differential energy-proper velocity distribution of the jet for the case 𝐿 j = 10 50 erg / s and 𝜃 j = 5 · . The hot injected jet accelerates quickly, and a high-velocity energetic tail forms within a few tens of ms. The energy distribution changes very little after a few hundreds of ms. In the case shown, the terminal Lorentz factor equals Γ ∞ = Γ j , 0 ℎ j ≃ 800. We compare the high-velocity tail to the original ejecta energy distribution (black dashed line in the figure), which has a cutoff around Γ 𝛽 = 0 . 2. Analyzing the jet energy-velocity distribution as a function of the jet scalar 𝑄 j we can associate each region of the curve to a different structure in the jet-cocoon system (bottom panel of Fig. 6). We identify three main peaks in this distribution: a first peak at low velocities, associated with the original ejecta (purple line), a second peak at around Γ 𝛽 ≃ 40, associated with the jet head and shocked jet ma-  \nre with the Γ 𝛽 maps in Fig. 4), a third peak at Γ 𝛽 ≃ 400, associated with unshocked jet material (blue line). Between the first (low velocity) and second peaks, the energy distribution shows a long plateau associated mostly with the turbulent region of the high-pressure cocoon surrounding the jet (yellow and green lines).  \nAfter an injection time of 𝑡 e = 0 . 3 s, the jet and the cocoon enshrining it are still inside the ejecta in all four cases with 𝜃 = 5 · -10 · and 𝐿 j = 10 50 -10 51 erg / s; however, it has already reached the outermost and shallowest region of the ejecta, which is expanding at ∼ 0 . 1 c. If we extend the simulation further, continuing our injection, at around 𝑡 e = 0 . 5 s, the jet breaks out of the ejecta. However, the jet engine can be switched off earlier, and the jet can still break out. As we switch off the engine, the last jet material (the jet tail ) is launched and moves close to the speed of light in the lowdensity, low-pressure interior of the jet, trying to catch up with the head, which is causally disconnected and continues to propagate as if the engine were still active. When the tail reaches the head, the jet is choked. A good marker of jet choking is the disappearance of the central region with the high proper velocity (i.e., Γ 𝛽 ≃ Γ ∞ ) and the sideways spreading of the material associated with the jet, which starts to dilate in vortexes. If the jet head continues to be causally disconnected from the tail while reaching the end of the ejecta (located, for 𝑡 d = 1 s, at around 𝑧 = 120 , 000 km), the jet is unchoked, and it escapes at a high proper velocity along the axis.  \nFig. 7 shows the evolution of two of the simulations shown in the first and second rows after 100 ms and 120 ms of injection, respectively. In both cases, we see a jet breakout from the outer ejecta. This was largely expected from the results of Eq. (5). To carry out a more detailed investigation, we stopped the jets of also the other simulations initialized with 𝑡 d = 1 s earlier than an injection time of 300 ms to see whether the results of Eq. (5) are consistent.  \nWe stopped the jets earlier according to the selected set of parameters ( 𝐿 j, 𝜃 j) shown in Table 1 and reported our results in the third and fourth columns of Fig. 2. Each panel shows a case with the same configuration parameters ( 𝐿 j, 𝑡 e) and two different values of 𝜃 j (5 · , 10 · ) . In the third column, we select the minimal engine time such that the jet with the narrower angle escapes; in the adjacent fourth column, we select the minimal engine time such that the wider jet escapes. The triangles indicate jets that could break out from the ejecta, while the stars indicate jets choked inside the ejecta. We can see that all points match the expected results for the jet breakout pretty well.", '3.3.2. Jet evolution for 𝑡 d = 0 . 1 s': 'We now explore the case where the jet is injected at 𝑡 d = 0 . 1 s after the merger. This allows us to study the  \nFigure 8. Jet injection for the model 𝑡 d = 0 . 1 s at 𝑡 e = 100 ms after the injection (end of injection) for two cases. Top : 𝐿 = 10 50 erg/s and 𝜃 = 5 · and bottom : 𝐿 = 10 51 erg/s and 𝜃 j = 10 · . White lines and plotted physical quantities are the same as the previous multipanel figures depicting simulation snapshots. For space reasons, we omitted plotting the density map. The top panels correspond to the top row, second column, of Fig. 2.  \n<!-- image -->  \njet breakout much earlier (which requires less injection time to burst out) but, conversely, from a more dense environment, which already shows the presence of the polar funnels from the compact object. We run again four cases with the same luminosities and opening angles as those with 𝑡 d = 1 s (as described in Table 1). Fig. 8 shows the status of the jets injected at 𝑡 d = 0 . 1 s after 100 ms of injection. In these cases, we also see that the jets drilled their way out of the ejecta, maintaining their opening angle. As predicted by Eq. (5), a much shorter engine time is required for the jet to escape since the outer shock of the ejecta 𝑅 ej is located much closer to the center at ≃ 21 , 000 km. Echoing the 𝑡 d = 1 simulations, we decided to test if the engine time required 𝑡 e to break out from the ejecta is consistent with Eq. (5). A very short burst ( 𝑡 e < 0 . 1 s) is sufficient in these cases. We stopped the injection time in all four cases at 𝑡 e = [ 0 . 01 , 0 . 03 , 0 . 05 , 0 . 1 ] s. The results are reported in the first and second columns of Fig. 2. Again, we express the choked jets as starred dots and the escaping ones  \n]  \nm  \nk  \n[  \nz  \nFigure 9. Jet injection with initial conditions 𝑡 d = 1 s for a jet with a wide initial opening angle of 𝜃 j = 35 · and 𝐿 j = 10 51 erg/s. The jet is shown in two different moments: (top row) at 𝑡 j = 25 ms and (bottom row) at 𝑡 j = 150 ms. The solid white line is the initial opening angle at which the jet is injected, and the dashed line in the bottom row shows the actual - approximately - opening angle of the jet ≃ 10 · . The external pressure could collimate the flow to a narrower angle in a few tens of milliseconds. Again, the four panels, from left to right, represent the Γ 𝛽 factor, the density 𝜌 , the pressure 𝑃 , and the jet tracer (flagged cells with jet material advected with the fluid)  \n<!-- image -->  \nas triangles in the plot. Also, in these cases, we verify that the predicted engine times are sufficient to break out from the ejecta.', '3.4. Initially wide jets': "As a last step in our analysis, we address the collimation effect of the ambient medium on a jet initialized with a very wide angle, i.e., 𝜃 j = 35 · , and a luminosity of 𝐿 j = 10 51 erg/s. The results are reported in Fig. 9. We plotted two different moments of the simulation: 25 ms after injection and 150 ms after injection. The collimating effect of the surrounding matter is clearly seen now. The ambient pressure of the ejecta can redirect the material's flux going sideways and collimate the jet to a much smaller opening angle in a few tens of milliseconds. A characteristic cusp-shaped collimation shock is seen at 𝑡 = 25ms. As expected, this leads to an initially  \nFigure 10. Jet injection with initial conditions 𝑡 d = 1 s for a jet with a wide initial opening angle of 𝜃 j = 20 · and 𝐿 j = 10 51 erg/s (top row), and 𝐿 j = 5 . 4 × 10 51 erg / s (bottom row), at 𝑡 j = 150 ms.  \n<!-- image -->  \ncylindrical jet in the region past the collimation shock (see also at 𝑡 = 25 ms). These results are consistent with Eq. (2) and subsequent Fig. 2 from which we see that the luminosity of the jet is not sufficient to maintain a wide opening angle; thus, the resulting jet will be collimated to a smaller angle.  \nThe jet later expands sideways to ∼ 10 · , filling the postmerger environment's low density-low pressure polar cavity. This cavity has an opening angle of ≃ 10 · , resulting in the jet quickly filling it and propagating with a reduced opening angle, presenting a quasi-cylindrical structure in its evolution (second row).  \nIn Fig. 10, we show two simulations with an initial jet opening angle of 𝜃 j = 20 · with two different luminosities: 10 51 erg / s and 5 . 4 × 10 51 erg / s. The latter value is chosen according to Eq. (3) for the collimation energy for the selected opening angle. The first case presents again a wide jet collimated to a smaller angle by the ejecta pressure; the second is a more energetic jet that expands almost conically at later times.  \nThis has a nontrivial implication: the only way to have very wide jets in this post-BNS merger environment is to inject a large amount of energy in a very short amount of time, i.e., the jet would require a high luminosity ( > 10 52 erg / s) to achieve a wide opening angle and counteract the presence of the strong  \nambient pressure. Thus, it is reasonable to suppose that the jet coming out from the post-BNS merger is relatively narrow ( 𝜃 j ≤ 10 · ) and in line with the observations of the afterglow of GW 170817 (Nakar & Piran 2020).", '3.5. The angular distribution': 'The angular distribution of the total energy of the outflow (the jet, the cocoon, and the ejecta) has significant implications for the EM signals: the prompt gamma rays that arose from the cocoon in the case of GRB 170817A, the observed afterglow that arose first from the cocoon and then from the jet, and a possible late radio afterglow arising from the rest of the ejecta.  \nFig. 11 depicts the differential isotropic equivalent energy distribution for four different runs with 𝜃 j = [ 3 , 5 , 6 . 8 , 10 , 20 ] · , and two luminosities 𝐿 j = 10 50 erg / s and 10 51 erg / s, taken at 𝑡 e = 300 ms. For each luminosity, the runs have all the same total energy 𝐸 j ( 𝑡 e ) = 3 × 10 49 erg and 3 × 10 50 erg, respectively. As expected, the narrower the opening angle, the more collimated the energy distribution is. For the 10 50 erg / s runs, the effect of the collimation by the ejecta is evident since the fit shows a smaller opening angle than the initial one, and the curve looks closer to a standard Gaussian. Overall, we find a core with rather sharp edges surrounded by a lower energy cocoon material (from the jet up to ∼ 50 · -60 · and ejecta material from ∼ 50 · -60 · to the equatorial plane).  \nWe fit these curves with a generalized Gaussian function:  \n𝐸 Ω ( 𝑡, 𝜃 ) = d 𝐸 d Ω ≡ 1 2 𝜋 sin 𝜃 d 𝐸 d 𝜃 = 𝐸 0 ( 𝑡 ) exp " -1 2 GLYPH<18> 𝜃 𝜃 c ( 𝑡 ) GLYPH<19> 𝑠 ( 𝑡 ) # , (6)  \nwhere 𝜃 c is the angular spread of the core of the jet, 𝐸 0 ( 𝑡 ) is proportional to the total energy of the jet at the time 𝑡 , and 𝑠 is a shape parameter that depends on time. This fit is suitable for the core region that is dominated by the jet and has a rather steep profile. We notice that the more the jet evolves, the more 𝜃 c saturates to 𝜃 j, the shape parameter 𝑠 stabilizes to a constant value, and 𝐸 0 increases following the same time evolution of 𝐸 j = 𝐿 j × 𝑡 . The jet\'s core angle 𝜃 c is in excellent agreement with the jet\'s initial opening angle within the error, and the wider the opening angle, the higher the jet\'s shape parameter is (indicating a flatter distribution within 𝜃 ≤ 𝜃 c ≃ 𝜃 j). A more in-depth overview of the parameters regulating the jet shape and their correlations can be found in Appendix D.3.', '3.6. Afterglow light curves': 'Weusethesimulationresults shown in the previous sections to calculate the light curves of the afterglow emission. We implement a model similar to the one described in Takahashi &Ioka (2020) (see also van Eerten et al. 2010, and references therein). We then compare the result to the actual data points  \nFigure 11. Jet energy per solid angle for jets with constant luminosities of 𝐿 j = 10 50 erg / s (top panel) and 10 51 erg / s (bottom panel). In each figure, we report the distributions for 𝜃 j = 3 · (orange line), 𝜃 j = 5 · (blue line), 6 . 8 · (green line), 10 · (red line), and 20 · (purple) with an injection time of 𝑡 e = 300 ms, and their respective generalized Gaussian fits (solid, dashed, dotted, dashed-dotted, solid black lines). The dashed colored lines for each run indicate the energy per solid angle of the sole jet material. The fit parameters refer to Eq. (6). The grey line represents the initial ejecta\'s energy per solid angle. The minima in the calculated distribution are most likely due to the artificial axial symmetry that produces a \'plug" ahead of the jet. The total energy for each jet model is 3 × 10 50 erg ) .  \n<!-- image -->  \nof the afterglow curves of GW170817. For a more detailed description of the methodology implemented to calculate the afterglow emission, including formulas, we redirect the reader to Appendix D.  \nFor the afterglow curve, the viewing angle is set at 𝜃 obs = 18 · as this is approximately the angle found by Hotokezaka et al. (2019); Ghirlanda et al. (2019) who take also into account the VLBI observations of the centroid motion (see also Nakar & Piran 2020, for a comparison of different estimates of the opening and viewing angles). We fix the spectral index at 𝑝 = 2 . 17, consistent with models of diffusive shockaccelerated supra-thermal electrons. The ISM density 𝑛 0, 𝜀 e (fraction of accelerated electrons energy density over the ki-  \nFig. 12 depicts the afterglow curves at 5 . 5 GHz for two values of the luminosity 𝐿 j = [ 10 50 , 10 51 ] erg / s and four lower resolution jet models with 𝜃 j = [ 3 , 5 , 6 . 8 , 10 ] · and 𝐿 j = 10 51 erg / s for a viewing angle of 𝜃 obs = 18 · . The afterglow curves are obtained using the 𝐸 Ω distribution at 𝑡 e = 300 ms for all models when the jet engine is still active. Thus, the energy deposited by the jet with 𝐿 j = 10 50 erg / s  \n<!-- image -->  \nT [day]  \nFigure 12. Afterglow light curves at 5 . 5 GHz for 𝜃 obs = 18 · of several jet models with different opening angle, 𝐿 j = 10 50 erg / s (top) and 𝐿 j = 10 51 erg / s (bottom), that use the 𝐸 Ω modeled after the fitting curves shown Fig. 11 (same color scheme). On the upper panel, we report only three synthetic light curves out of four since the model for 𝜃 j = 10 · with 𝐿 j = 10 50 erg / s has unrealistic values for 𝑛 0 and 𝜀 B . The radio points and relative fit (black dashed line) are taken from Troja et al. (2019). The dashed lines are the afterglow from the generalized Gaussian fits only, while the solid lines are the light curves of the jet, including the energy distribution outside the Gaussian fits. For comparison, alongside the value of the jet opening angle, we report 𝜃 c from the fit.  \nnetic energy density at the shock), and 𝜀 B (the fraction of the magnetic energy density over the kinetic energy density at the shock), are linked by a degeneracy relation, which gives the dependency of one of the three parameters of our choice concerning the other two, for a fixed jet model and a fixed 𝜃 obs.  \nTable 2. Afterglow parameters used in Fig. 12 to calculate the light curves at 5 . 5 GHz. For brevity, we omitted the units for 𝐿 j (in erg / s), 𝐸 j (in erg), and 𝑛 0 (in cm -3 ).  \n(10 51 erg / s) is 𝐸 j = 3 × 10 49 erg (3 × 10 50 erg). Switching off the injection for the jets at this time will not change their 𝐸 Ω profile too much as long as the jet has already crossed most of the ejecta and breaks out. The curves are compared to the fits obtained by Mooley et al. (2018a) (black dashed curve) and Takahashi & Ioka (2020), which uses the average afterglow parameters of Troja et al. (2019).  \nFig. 12 depicts both the afterglow light curves obtained considering the full 𝐸 Ω distribution (solid curves) and the generalized Gaussian fit (dashed curves) that considers only the jet material. In the latter case, we have excluded the slow material associated with the ejecta, selecting only the fast matter with a jet tracer value 𝑄 j > 10 -5 and Γ 𝛽 ≥ 1.  \nThe first thing we notice is that, while there was no attempt to guide the calculations to fit a given observed afterglow (apart from setting 𝑛 0 and 𝜀 𝐵 to fit the time of the peak and its magnitude), the model with 𝜃 j = 6 . 8 · for the 10 50 erg/s jet gives a good fit to the observed light curve. The effective opening angle of the jet after collimation is 𝜃 c = 3 . 8 · . Thus, the ratio of the jet opening angle to the viewing angle, 𝜃 j / 𝜃 obs, is consistent with the one obtained for synthetic jets estimates (see Nakar & Piran 2020). The 10 51 erg/s jet is less collimated, and as it is too wide, 𝜃 c = 6 . 8 · , the resulting afterglow is somewhat inconsistent on the declining part of the afterglow.  \nTherising portion of the light curves for the jets with 𝜃 j = 3 · and 5 · are below the observed data. Additionally, these jets have a narrow peak. This is due to the sharp edges of the energy distribution that decline faster than a Gaussian. For 𝜃 j = 10 · and 𝐿 = 10 51 erg/s, the peak becomes much too wide and overestimates both the rising and the declining phase of the curve at early and later times. The parameters of these curves are found in Table 2, in which, for simplicity, we keep 𝜀 e = 0 . 1 constant for each model and we vary 𝑛 0 and 𝜀 B to adjust the position of the peak time 𝑇 peak and the height of the peak flux F 𝜈, peak.', '4. SUMMARY AND CONCLUSIONS': 'We presented a set of 2.5D special relativistic simulations of unmagnetized relativistic jets launched into a realistic postBNSmergerenvironment based on the neutrino radiation GRMHDsimulation by Kiuchi et al. (2023) that was evolved for more than 1 second after the merger. We explored a variety of configurations for relativistic jets launched in proximity of a compact object and expanded in a post-merger environment. We verified the semi-analytical prediction presented in Pais et al. (2023) about the jet escape and collimation and validated our results for a jet launched at 0 . 1 s and 1 s after the merger. We summarise our findings as follows:  \n- · The low-pressure, low-density polar funnel that arises in the BNS simulations makes it easier for the jet to emerge from the ejecta.\n- · For wide jets, the ambient pressure of the ejecta counterbalances the pressure of the jet while it propagates and collimates it.\n- · We confirm the predictions derived by Gottlieb & Nakar (2022) for a minimum engine operative time for the jet to escape. The numerical simulations verified that the predicted engine time 𝑡 e allows the jet to escape from the system after the shutdown of the central engine driving it.\n- · We explored the conditions of a late jet launch at 𝑡 d = 1 s and an early jet launch at 𝑡 d = 0 . 1 s, finding a significantly shorter engine time was required to escape in case of an early launch at 𝑡 d = 0 . 1 s, i.e., 0 . 01 -0 . 05 s vs. 0 . 1 -0 . 3 s for a launch with a delay of 𝑡 d = 1 s.\n- · We verified that the choice of the enthalpy (and thus of the terminal Lorentz factor) does not influence our results as long as the jet is hot.\n- · We tested the effect of the ejecta thermal pressure in collimating the jet injecting two very wide jets ( 𝜃 j = 20 · ≃ 0 . 35 rad and 35 · ≃ 0 . 6 rad) and following their evolution in the ejecta. We showed how the dense inner ejecta can collimate the jet to an opening angle of about 10 · in a few tens of milliseconds while the jet is still located close to the compact object where it is surrounded by dense ejecta matter.\n- · Wecompared our synthetic light curves generated from the energy profile of our simulations with the observed afterglow data of GW170817 at 5 . 5 GHz. Our isotropic equivalent energy distributions are characterized by sharp edges around the jet\'s opening angle and decline faster than a Gaussian. Nevertheless, the resulting synthetic light curves can match radio data for collimated jets with 𝐸 j = 3 × 10 49 erg / s, 𝜃 j = 6 . 8 · (corresponding to 𝜃 c = ( 3 . 8 ± 0 . 3 ) · ) and a viewing angle of  \n𝜃 obs = 18 · . In contrast, for smaller (larger) jets\' opening angles, the light curve\'s peak becomes too narrow (wide). These values are roughly in agreement with those obtained by fitting jet models to the light curve (Hotokezaka et al. 2019; Ghirlanda et al. 2019; Lamb et al. 2019b; Takahashi & Ioka 2020). The more powerful 𝐸 j = 3 × 10 50 erg / s jets are uncollimated (i.e., they roughly maintain their original opening angles). These jets decline or rise too steeply and do not match the radio data well.  \n- · These calculations that follow the jet propagation in the ejecta enable us to break the degeneracy between the total energy of the jet, 𝐸 j, and the external density, 𝑛 0 that exists in ad hoc fittings of jet profiles to the light curve (e.g. Nakar & Piran 2020).  \nThe novelty of this work is that we can study the jet launch much later than previous works due to the 1-second-long evolved initial conditions, which tracked the ejecta\'s expansion into the circumbinary environment. As we have seen, it is more difficult for the jet to break out from the ejecta in case of a late launch, thus putting stronger constraints on the jet engine. This is important, at least in the context of modeling GW170817A, as it has been suggested that there was about a 1-sec delay between the merger and the jet launch (Kasliwal et al. 2017; Gottlieb et al. 2018; Nakar et al. 2018) in that event.  \nBecause of this long duration, there is an external decelerating shock wave propagating outwards at a velocity of 0 . 1 c in the artificially imposed \'numerical vacuum" that has a floor value for the density of ≈ 0 . 17 g / cm 3 . This does not influence our main concerns here: jet collimation and breakout. However, it influences the fine details of the shock breakout from the ejecta that is needed in case one wishes to calculate the signature of the cocoon breakout - which was, most likely, the cause of the observed 𝛾 -rays in GRB 170817A (Kasliwal et al. 2017; Gottlieb et al. 2018; Nakar et al. 2018).  \nModeling the observed sGRB detected in GRB 170817A wouldrequire following up the shocked cocoon until it catches up with the low mass very-high velocity component of the ejecta so that the optical depth of the shock is of order 𝜏 ≃ 1. This would require a longer run of our simulations, which  \nwould take at least 2 seconds. However, for a realistic simulation, the external \'numerical vacuum" in such a case should be muchlowersothatitdoesnotslowdownthefastcomponentof the ejecta. Namely, the circumbinary atmosphere surrounding the post-merger ejecta has to have a mass smaller than the fast-moving ejecta (i.e., 𝑀 atmosphere < 10 -4 M ⊙ ) within the simulation box of about 10 11 -10 12 cm. We leave this analysis to future work. Additional future work would involve exploring the propagation of a Poynting flux-dominated jet in this ejecta.  \nIn this work, we have demonstrated that for jet luminosities of 10 50 -10 51 erg/s an engine time 𝑡 e of 0 . 1-0 . 3 s are required for a late launch ( 𝑡 d = 1 s) and 𝑡 e ∼ 0 . 01-0 . 05 s in the case of an early launch with a time delay of 𝑡 d = 0 . 1 s. Recalling that the jet in GW 170817 was less powerful compared to typical sGRBs, we note that a longer duration would be required for less powerful jets to break out. The jet coming out from the system is narrowly collimated. The collimation is supported by the narrow funnel in the initial outflow configuration and does not exceed an opening of 𝜃 j ≲ 10 · . A wide jet injected into the ejecta is quickly collimated to a smaller opening angle.', 'ACKNOWLEDGEMENTS': "We thank the anonymous referee for their constructive report, which helped improve this work's quality. MP thanks Riccardo Ciolfi for his valuable comments. This work was supported by the ERC grants Multĳets and by the SimonsCollaboration on Extreme Electrodynamics of Compact Sources (SCEECS) (TP, MP), and by Grant-in-Aid for Scientific Research (Grant Nos. JP20H00158 and JP23H04900) of Japanese MEXT/JSPS (MS), and by Grant-in-Aid for Scientific Research (Grant Nos. JP23K25869, and JP23H01172) of Japanese MEXT/JSPS (KK). This work used the computational resources of the supercomputer Fugaku, which RIKEN provided through the HPCI System Research Project (Project ID: hp220174). The simulation was also performed on Sakura, Cobra, and Raven clusters at the Max Planck Computing and Data Facility and the Cray XC50 at CfCA of the National Astronomical Observatory of Japan.", 'DATA AVAILABILITY': 'The data that support the findings of this study are available upon request from the authors.', 'REFERENCES': 'Abbott, B., Abbott, R., Abbott, T., et al. 2017a, Physical Review Letters, 119, 161101, doi: 10.1103/PhysRevLett.119.161101 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, The Astrophysical Journal, 848, L12, doi: 10.3847/2041-8213/aa91c9  \nAguilera-Miret, R., Palenzuela, C., Carrasco, F., Rosswog, S., &  \nVigan ò , D. 2024, Delayed jet launching in binary neutron star mergers with realistic initial magnetic fields, arXiv.  \nhttp://arxiv.org/abs/2407.20335  \n- Taylor, G. 1950, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 201, 175. https://www.jstor.org/stable/98396\n- Tominaga, N. 2009, The Astrophysical Journal, 690, 526, doi: 10.1088/0004-637X/690/1/526\n- Troja, E., Piro, L., van Eerten, H., et al. 2017, Nature, 551, 71, doi: 10.1038/nature24290\n- Troja, E., van Eerten, H., Ryan, G., et al. 2019, Monthly Notices of the Royal Astronomical Society, 489, 1919,  \ndoi: 10.1093/mnras/stz2248  \n- van Eerten, H., Zhang, W., & MacFadyen, A. 2010, The Astrophysical Journal, 722, 235,\n- doi: 10.1088/0004-637X/722/1/235\n- Wang, P., Abel, T., & Zhang, W. 2008, The Astrophysical Journal Supplement Series, 176, 467, doi: 10.1086/529434\n- Yang, B., Jin, Z.-P., Li, X., et al. 2015, Nature Communications, 6, 7323, doi: 10.1038/ncomms8323\n- Zenati, Y., Krolik, J., Werneck, L., et al. 2024, The Dynamics of Debris Disk Creation in Neutron Star Mergers,  \ndoi: 10.48550/arXiv.2404.03156  \n- Zenati, Y., Krolik, J. H., Werneck, L. R., et al. 2023, The Astrophysical Journal, 958, 161, doi: 10.3847/1538-4357/acf714  \nZhang, W., Woosley, S. E., & Heger, A. 2004, The Astrophysical Journal, 608, 365, doi: 10.1086/386300', 'A.1. Environment data': 'The simulation domain of the initial conditions based on Kiuchi et al. (2023) comprises 13 levels of a fixed mesh refinement (FMR) embedded in a Cartesian geometry. The finest FMR domain has the size of 𝐿 ∈ [-37 . 875 , 37 . 875 ] km with the grid spacing of Δ 𝑥 = 0 . 15 km, while the coarsest has a grid size of 𝐿 ∈ [-1 . 55136 , 1 . 55136 ] × 10 5 km with a grid spacing of Δ 𝑥 = 614 . 4 km.  \nTo import the initial conditions data into /p.pc/l.pc/u.pc/t.pc/o.pc and run the simulation, we proceeded as follows:  \n- 1. for each refinement level, we create a cylindrical average around the 𝑧 -axis for each cell, obtaining a 2D cylindrically symmetric Cartesian grid. To run the code in 2.5D we retain the 𝑣 𝜙 components (which are also averaged);\n- 2. we merge the 13 refinement levels into one, interpolating linearly between the regions where the grid resolution doubles in the original simulation grid;\n- 3. we import the data into /p.pc/l.pc/u.pc/t.pc/o.pc creating a mock grid file that reflects the actual center of each cell;  \nAfter the merger, we imported the initial conditions at 𝑡 d = 0 . 1 s and 𝑡 d = 1 s. The first snapshot at 0 . 1 s represents a case where the ejecta did not expand too further into the dense circumbinary environment ( 𝜌 floor = 0 . 167 g / cm 3 ) of the simulation box. The ejecta outer radius is about 2 × 10 4 km from the center of the simulation box; the second snapshot represents the system at the end of the original simulation, with the outer radius at 1 . 2 × 10 5 km from the center.', 'A.2. Computational grid': 'For the simulation in /p.pc/l.pc/u.pc/t.pc/o.pc we use a grid of size 3549 × 5541 cells, with the radial cylindrical coordinate /two.sup extending within the range 𝑟 = [ 0 , 1 . 5 × 10 5 ] km and the vertical coordinate extending within the range 𝑧 = [ 𝑧 0 , 1 . 5 × 10 5 ] km. /three.sup . Our choice of a very low initial height to inject the jet is dictated by the fact that we want to explore the effect of the collimation in the innermost region of the ejecta, which retains most of the cavity. We start our jet to a distance comparable to the light cylinder radius of the compact object ( ≃ 40 km).  \nWeuseacombination of a uniform and a non-uniform mesh grid in 𝑟 -𝑧 coordinates with a decreasing resolution from the inner region of the simulation box to the outer boundaries. The grid mesh is uniform in the inner part to maintain a high resolution of the jet injection and the formation of the resulting high-pressure cocoon. The uniform mesh has 1000 × 4000 grid points extending in the ranges 𝑟 = [ 0 , 100 ] km and 𝑧 = [ 𝑧 0 , 5000 ] km with a resolution along both coordinates of Δ ( 𝑟, 𝑧 ) unif . = [ 0 . 1 , 1 . 2375 ] km, allowing enough mesh coverage for the nozzle region. Next to the uniform mesh, we placed a stretched mesh with 2549 × 1541 grid points extending along both coordinates within the range ( 𝑟, 𝑧 ) = [ 100 , 1 . 5 × 10 5 ] × [ 5000 , 1 . 5 × 10 5 ] km with a stretching ratio of ≃ 1 . 0033. The number of grid points for this mesh is chosen such that its initial grid spacing is the same as the adjacent uniform mesh Δ ( 𝑟, 𝑧 ) s , init = Δ ( 𝑟, 𝑧 ) unif . = 0 . 1kmand its final grid spacing is Δ ( 𝑟, 𝑧 ) s , final = 500 km. This ensures a smooth resolution increase without jumps for the entire simulation grid.  \nAfter switching off the jet injection, we lower the resolution of our computational domain using a 𝑟 × 𝑧 grid with a combination of uniform and geometrically stretched patches. The uniform grid 1000 × 900 grid extends in the range ( 𝑟, 𝑧 ) = [ 0 , 5000 ] km × [ 𝑧 0 , 5000 ] km, with a resolution Δ 𝑥 = 5 km, while the geometrically stretched grid with 1352 2 cells extends in the range ( 𝑟, 𝑧 ) = [ 5000 , 1 . 5 × 10 5 ] km in both directions.', 'B.1. Jet injection': 'The jet is characterized by four parameters: luminosity 𝐿 j, opening angle 𝜃 j, enthalpy ℎ j, and the jet engine duration 𝑡 e. In our setup, we inject the jet cylindrically through a nozzle such that the initial Lorentz factor Γ 0 , j is related to the initial choice of the jet opening angle via Γ 0 , j ≃ 1 /( 1 . 4 𝜃 j ) (Mizuta & Ioka 2013). The following stress-energy tensor gives the energy density of a relativistic, unmagnetized fluid:  \n𝑇 𝜇𝜈 = 𝜌ℎ𝑐 2 𝑢 𝜇 𝑢 𝜈 + 𝑃𝑔 𝜇𝜈 (B1)  \nwhere 𝑢 𝜇 is the 4-velocity and 𝑔 𝜇𝜈 is the metric tensor (which coincides with the Minkowski space-time in our case). In the laboratory frame, where 𝑢 𝜇 = Γ ( 1 , β ) , the energy density of the fluid assumes the form /four.sup :  \n𝑇 00 = 𝜌ℎ𝑐 2 Γ 2 -𝑃 (B2)  \nThe jet is injected through a nozzle with radius 𝑟 j and area Σ j = 𝜋𝑟 2 j . At the nozzle, we inject a flux of energy per unit time such that:  \n𝐿 j = ∫ 𝑟 j 0 𝑇 00 j 𝛽 jd Σ j = ( 𝜌 j ℎ j 𝑐 2 Γ 2 0 , j -𝑃 j ) 𝛽 0 , j Σ j , (B3)  \nwhere we imposed a top-hat profile of the jet (i.e., the properties of the jet do not vary with the position within the nozzle radius). From Eq. (B3), we can derive the value for 𝜌 j. The expression of the enthalpy links pressure and density: ℎ = 1 + 𝜀 + 𝑃 / 𝜌 , where 𝜀 ( 𝑃, 𝜌 ) is the internal specific energy density of the fluid. The explicit expression for 𝜀 ( 𝑃, 𝜌 ) depends on the choice of the EOS. In the case of a perfect fluid, we have 𝜌𝜀 = 𝑃 /( 𝛾 -1 ) , where the adiabatic factor 𝛾 = 4 / 3 for a relativistic fluid, such that 𝜌ℎ𝑐 2 = 𝜌𝑐 2 + 4 𝑃 . In the case of a Taub-Matthews EOS for a Synge gas, we use the approximation showed in Mignone & McKinney (2007) where the product 𝜌𝑐 2 ℎ is expressed by  \n𝜌ℎ𝑐 2 = 5 2 𝑃 + √︂ 𝜌 2 𝑐 4 + 9 4 𝑃 2 (B4)  \nInverting Eq. (B4) and since we fix the enthalpy when we inject the jet, we find the following expression for the pressure:  \n𝑃 = 1 8 GLYPH<16> 5 ℎ -√︁ 9 ℎ 2 + 16 GLYPH<17> 𝜌𝑐 2 = ℎ eff 𝜌𝑐 2 (B5)  \nwhich reduces to 𝜌ℎ𝑐 2 / 4 for a relativistic gas. Since we inject our jet vertically through the nozzle, the only nonzero component of 𝑢 𝜇 is 𝛽 𝑧 . Since we fix the value of the luminosity 𝐿 j and the enthalpy ℎ j the quantities injected at the excision height 𝑧 0 and for 𝑟 ≤ 𝑟 j are:  \n𝛽 𝑟, j = 𝛽 𝜙, j = 0 (B6)  \n𝛽 𝑧, j ≡ 𝛽 0 , j = √︃ 1 -1 / Γ 2 0 , j (B7)  \n𝜌 j = 𝐿 𝑗 Σ j ( ℎ j Γ 2 0 , j -ℎ eff ) 𝛽 𝑧, j 𝑐 3 (B8)  \n𝑃 j = ℎ eff 𝜌 j 𝑐 2 (B9)  \nEqs. (B6)-(B9) describe the jet injection at the nozzle and are constant in time by the entire duration of the engine activity as long as 𝐿 j = const . , otherwise 𝜌 j and 𝑃 j will vary according to the time dependency of 𝐿 j.  \nAdditionally, we enforce that, after 𝑞 time steps, the increase in the total energy over time follows the temporal law we specified for the jet luminosity. To do so, we measure the total energy added to the system up to that point on top of the ejecta energy of the initial conditions (i.e., 𝐸 added = 𝐸 measued -𝐸 ej ) , and we compare it to the total energy we expect up to that moment integrating the luminosity function of our choice, i.e. 𝐸 ( 𝑡 ) = ∫ 𝑡 0 𝐿 j ( 𝑡 ) d 𝑡 . For a constant luminosity, the energy added to the system is trivial  \n𝐸 ( 𝑡 ) = 𝐿 j × 𝑡 . If the energy in the system does not match the expected energy, we re-adjust iteratively the value of the jet pressure (and thus the jet density) at the injection point according to the following scheme:  \n𝑃 j ,𝑛 + 𝑞 =            ( 1 + 𝑥 ) 𝑃 j ,𝑛 if 𝐸 added < 𝐸 ( 𝑡 ) 𝑃 j ,𝑛 if 𝐸 added = 𝐸 ( 𝑡 ) 𝑃 j ,𝑛 /( 1 + 𝑥 ) if 𝐸 added > 𝐸 ( 𝑡 ) (B10)  \nwhere 𝑃 j ,𝑛 is the jet pressure at the nozzle at the current time step 𝑛 , while 𝑛 + 𝑞 indicates the 𝑞 -th time step after 𝑛 . The quantity 𝑥 = 10 -3 is a small correction we give to the pressure at each time step to match the desired total energy of the system up to that point. In our simulations, we found that 𝑞 = 500 was sufficient to control the pressure value over time.', 'B.2. Energy calculation': "Eq. (B1) expresses the fluid's total energy density. If we plug the definition of specific dimensionless enthalpy ℎ = 1 + 𝜀 / 𝑐 2 + 𝑃 / 𝜌𝑐 2 , we can rewrite it as:  \n𝑇 00 = 𝜌𝑐 2 Γ 2 +( 𝜌𝜀 + 𝑃 ) Γ 2 -𝑃 = 𝐷𝑐 2 ( Γ -1 )+ 𝐷𝑐 2 +( 𝜌𝜀 + 𝑃 ) Γ 2 -𝑃 (B11)  \nwhere we introduced the laboratory density 𝐷 = 𝜌 Γ . Writing the total energy density in this way allows us to distinguish different contributions:  \n𝑒 kin = 𝐷𝑐 2 ( Γ -1 ) (B12)  \n𝑒 int = ( 𝜌𝜀 + 𝑃 ) Γ 2 -𝑃 = ( ℎ -1 ) 𝜌𝑐 2 Γ 2 -𝑃 = GLYPH<18> 𝛾 𝛾 -1 Γ 2 -1 GLYPH<19> 𝑃 (B13)  \nwhich are the kinetic and internal energy of the fluid in the lab frame, respectively. In case of a Taub-Matthews EOS we can still find an effective adiabatic factor 𝛾 that allows writing Eq. (B13) in the form shown and, generalizing, we can write 𝛾 /( 𝛾 -1 ) = ( ℎ -1 ) 𝜌𝑐 2 / 𝑃 (Mignone & McKinney 2007). Integrating Eqs. (B12) and (B13) over the entire volume gives the total energy of the relativistic fluid in the simulation measured in the laboratory frame.", 'C. JETS WITH SAME TERMINAL VELOCITY': 'Theprevious runs we presented for different opening angles were initialized so that the initial opening angle 𝜃 j determines the initial Lorentz factor via Γ 0 , j ≃ 1 /( 1 . 4 𝜃 j ) ≃ [ 8 , 4 ] (for 𝜃 0 , j = [ 5 · , 10 · ] , respectively), and keeping the enthalpy of the jet fixed at ℎ j = 100. Since Γ ∞ ≃ Γ 0 , j ℎ j, the terminal Lorentz factor of the system depends on the initial Lorentz factor and the specific enthalpy used to initialize the jet. We show here runs where we keep the terminal Lorentz factor constant, varying the specific enthalpy according to the choice of the initial Lorentz factor, i.e., using ℎ j = Γ ∞ / Γ 0 , j, where Γ ∞ =  \nFigure 13. Differential energy distribution with respect to the proper velocity Γ 𝛽 of relativistic jets with 𝐿 j = 10 51 erg / s and 𝜃 j = [ 3 , 5 , 6 . 8 , 10 ] · after 300 ms of injection. The solid lines represent runs with ℎ j = 100 and Γ j , 0 = [ 12 , 8 , 6 , 4 ] , respectively, such that the final Lorentz factor Γ ∞ = Γ j , 0 ℎ j is not fixed. The dashed lines represent runs with the same opening angles but fixed Γ ∞ = 500.  \n<!-- image -->  \nj = 3  \n,  \n=1200  \nj = 3  \n,  \n=500  \nFigure 14. Comparison, after 𝑡 j = 300 ms, between the velocity maps of two jets with the same opening angle 𝜃 j = 3 · , 𝐿 j = 10 51 erg / s, but different Γ ∞ . We modified the color map and scale with respect to the other velocity maps to enhance the difference between the two runs. Notice the similarity in the overall structure.  \n<!-- image -->  \n500, which is value high enough to avoid cold jets. We run three cases with the same luminosity, fixed at 𝐿 j = 10 51 erg / s, and different initial opening angles 𝜃 j = [ 3 , 5 , 6 . 8 , 10 ] · . The choice of the angles is such that they correspond to initial Lorentz factors of Γ j , 0 = [ 12 , 8 , 6 , 4 ] , respectively. In Fig. 13, weshowthedifferential distribution of the energy with respect to the proper velocity where we compared the cases with ℎ j = 100 (solid lines) with the cases with ℎ j = Γ ∞ / Γ 0 , j where Γ ∞ = 500 (dashed lines). We notice that the energy distribution is essentially the same between the two types of runs except for the ending portion of the tail, where energy accumulates.  \nIn Fig 14, we show two jets with the same luminosity, opening angle 𝜃 j = 3 · but different terminal velocity (i.e., different enthalpy). We chose a small opening angle to emphasize the difference in 𝑢 ∞ = 1200 vs. 500, which translates into ℎ j = 100 vs. ≃ 42. Both jets are hot enough, so most of the initial energy is injected into the thermal channel. The jet expansion, head velocity, and drilling through the ejecta are very similar in both cases. We conclude that choosing a fixed Γ ∞ , as long as the jet is hot, does not affect our results and runs.', 'D.1. Framework': "The synchrotron emission of the GRB afterglow is modeled after Sari et al. (1998) and Granot et al. (1999).  \nThe received flux is given by the volume integral:  \nF 𝜈 ( 𝑇, 𝜃 obs ) = 1 4 𝜋𝐷 2 ∫ 2 𝜋 0 d 𝜙 ∫ 𝜃 max 0 d 𝜃 sin 𝜃 𝑅 2 Δ 𝑅 𝜀 ' 𝜈 ' Γ 2 ( 1 -𝛽𝜇 ) 2 GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝑡 = 𝑇 + 𝜇𝑅 / 𝑐 ≃ 1 4 𝜋𝐷 2 ∫ 2 𝜋 0 ∫ 𝜃 max 0 d 𝜙 d 𝜃 sin 𝜃𝑅 3 𝜀 ' 𝜈 ' 12 Γ 4 ( 1 -𝛽 sh 𝜇 )( 1 -𝛽𝜇 ) 2 GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> 𝑡 = 𝑇 + 𝜇𝑅 / 𝑐 (D14)  \nwhere we assumed that the observed radiation is produced by a thin shell of width Δ 𝑅 ≃ 𝑅 /[ 12 Γ 2 ( 1 -𝛽 sh 𝜇 )] ≪ 𝑅 sh, smaller than the size of the shock radius. In the formula, 𝐷 is the distance from the source, 𝜀 ' 𝜈 ' is the comoving frame emissivity, and 𝜇 ( 𝜃, 𝜙 ) is the angle between observer direction and local velocity, defined as:  \n𝜇 = cos 𝜃 cos 𝜃 obs + sin 𝜃 sin 𝜃 obs cos 𝜙 (D15)  \nFurthermore, the integral is calculated over different emission times, which depend on the position ( 𝜃, 𝜙 ) . The implicit formula  \n𝑡 = 𝑇 + 𝜇 𝑐 𝑅 ( 𝑡 ) (D16)  \nrelates the laboratory time 𝑡 at each position ( 𝜃, 𝜙 ) to a given observer time 𝑇 . Solving it we obtain 𝑡 = 𝑡 ( 𝑇, 𝜃, 𝜙 ) for fixed 𝜃 obs. The functional form of 𝑅 ( 𝑡 ) is found taking some approximations: we assume that the jet expands in a cold and  \nFigure 15. Comparison between the formula used in van Eerten et al. (2010) for the shock proper velocity and our broken power law expressed by Eq. (D17). The curves are plotted in dimensionless time 𝑡 ∗ and are scale-free. The dashed vertical lines indicate the two transition times 𝑡 1 and 𝑡 2 between the relativistic, trans-relativistic, and the trans-relativistic and Newtonian regimes, respectively. The horizontal dotted lines represent the values of 𝑢 sh at the transition times. Our formula interpolates the two asymptotic regimes for 𝑡 ∗ → 0 and 𝑡 ∗ →∞ with a softer power law between the two with 𝑢 ∝ 𝑡 -0 . 9 .  \n<!-- image -->  \nhomogeneous interstellar medium; the jet is assumed to be axisymmetric, the energy of the blast wave is angle-dependent, and each segment expands spherically with a portion of the total energy given by the local isotropic equivalent energy 𝐸 Ω ≡ d 𝐸 ( 𝜃 )/ d Ω = 𝐸 0 𝑓 ( 𝜃 ) .  \nEach shock segment is described by an ultra-relativistic blast wave, which gradually decelerates. We distinguish three regimes: i) R: relativistic ( Γ 𝛽 > 1), in which the evolution of the (proper) velocity of the shock is described by the selfsimilar solution of Blandford & McKee (1976); ii) TR: transrelativistic (0 . 1 ≲ Γ 𝛽 ≤ 1), which marks the gradual passage from a relativistic blast wave to a Newtonian one and. iii) NR: non-relativistic ( Γ 𝛽 ≲ 0 . 1), in which the shock dynamics is described by the Sedov-Taylor solution (Taylor 1950; Sedov 1959).  \nUnlike van Eerten et al. (2010), we do not combine the two self-similar solutions by summing their squares under the square root. Instead, we use a broken power law for the proper velocity to describe the three phases of the shock radius evolution. The self-similar solution for Γ 𝛽 reads  \n( Γ 𝛽 ) sh ≡ 𝑢 sh =            𝐶 BM 𝑡 -3 / 2 ( R ) , if Γ 𝛽 > 1 𝐶 3 / 5 BM 𝑡 -9 / 10 ( TR ) , if 0 . 12 ≲ Γ 𝛽 ≤ 1 𝐶 ST 𝑡 -3 / 5 ( NR ) , if Γ 𝛽 ≲ 0 . 12 (D17)  \nwhere  \n𝐶 BM = √︄ 17 𝐸 Ω 8 𝜋𝑛 0 𝑚 p 𝑐 5 , 𝐶 ST = 2 5 × 1 . 15 GLYPH<18> 𝐸 Ω 𝑛 0 𝑚 p 𝑐 5 GLYPH<19> 1 / 5 (D18)  \nare the coefficients of the Blandford-McKee and the SedovTaylor self-similar solutions, respectively.  \nThe trans-relativistic (TR) regime is obtained by modifying the Blandford-McKee solution raised to the power of 3 / 5, which imposes a different slope to make the transition to a non-relativistic regime smoother. The transition times between the different regimes are found imposing ( Γ 𝛽 ) R = ( Γ 𝛽 ) TR, and ( Γ 𝛽 ) TR = ( Γ 𝛽 ) NR. We find  \n𝑡 R → TR = 𝐶 2 / 3 BM , 𝑡 TR → NR = 𝐶 2 BM 𝐶 -10 / 3 ST (D19)  \nWe can rewrite the formulas more concisely, introducing a characteristic time scale for the blast wave segment such that  \n𝑡 = 𝛼𝑡 ∗ , 𝑅 = 𝛼𝑐𝑅 ∗ , (D20)  \nwhere  \n𝛼 ( 𝜃 ) = GLYPH<18> 𝐸 Ω 𝑛 0 𝑚 p 𝑐 5 GLYPH<19> 1 / 3 = 337 d GLYPH<18> 𝐸 0 𝑓 ( 𝜃 ) 10 51 erg GLYPH<19> 1 / 3 GLYPH<16> 𝑛 0 1 cm -3 GLYPH<17> -1 / 3 , (D21)  \nwhich depends on 𝜃 via the local isotropic energy. We can rewrite the coefficients in Eq. (D18) as 𝐶 BM , ∗ = √︁ 17 / 8 𝜋 and 𝐶 ST , ∗ = 0 . 46, and the transition times as:  \n𝑡 1 , ∗ = 𝑡 R → TR , ∗ = GLYPH<18> 17 8 𝜋 GLYPH<19> 1 / 3 ≃ 0 . 878 , 𝑡 2 , ∗ = 𝑡 TR → NR , ∗ = 17 8 𝜋 GLYPH<18> 2 5 × 1 . 15 GLYPH<19> -10 / 3 ≃ 9 . (D22)  \nSubstituting, these values in Eq. (D17) correspond to the values of the proper velocity:  \n𝑢 R → TR = 1 , 𝑢 TR → NR = 𝐶 3 ST 𝐶 -6 / 5 BM ≃ 0 . 123 (D23)  \nA comparison between the proper velocity formula used by van Eerten et al. (2010) and our approximation is shown in Fig. 15.  \nTheLorentzfactor of the fluid is related to the Lorentz factor of the shock via the jump conditions, giving 𝑢 = 𝑢 sh / √ 2 in the relativistic and trans-relativistic regimes, and 𝑢 = 3 𝑢 sh / 4 in the non-relativistic regime. The shock radius is given by  \n𝑅 sh ( 𝑡 ) = 𝑐 ∫ 𝑡 0 𝛽 sh ( ˜ 𝑡 ) d˜ 𝑡 -→ 𝑅 sh , ∗ ( 𝑡 ∗ ) = ∫ 𝑡 ∗ 0 𝛽 sh ( ˜ 𝑡 ) d˜ 𝑡 (D24)  \n<!-- image -->  \n<!-- image -->  \nFigure 16. Peak position, peak height and peak spread of the afterglow curves as a function of the viewing angle 𝜃 obs , for four different opening angles: 3 · (blue points), 5 · (orange), 6 . 8 · (green), 10 · (red). We superposed to the points the power-law fitting lines with the same color as the points. The viewing angle is normalized with respect to the jet opening angle of the corresponding run. In contrast, 𝑇 peak and F 𝜈, peak , for simplicity, are normalized with respect to their values for 𝜃 obs = 60 · (the last point of the series).  \n<!-- image -->  \nwith 𝛽 sh = 𝑢 sh / √︃ 1 + 𝑢 2 sh . For each regime, the integration is analytical and yields to  \n𝑅 sh , R , ∗ = 2 𝐹 1 1 3 , 1 2 , 4 3 ; -𝑡 3 ∗ 𝐶 2 BM ! 𝑡 ∗ , 𝑅 sh , TR , ∗ = 10 𝐶 3 / 5 BM 2 𝐹 1 -1 18 , 1 2 , 17 18 , -𝐶 6 / 5 BM 𝑡 9 / 5 ∗ ! 𝑡 1 / 10 ∗ 𝑅 sh , NR , ∗ = 𝐶 ST 𝑡 2 / 5 ∗ (D25)  \nwhere 2 𝐹 1 ( . . . ) is Gauss' hypergeometric function. The integral in Eq. (D24) can be broken into the individual contributions of the terms expressed by Eq. (D25) for the different regimes:  \n𝑅 sh , ∗ =              𝑅 sh , R , ∗ ( 𝑡 ∗ ) 𝑡 ∗ ≤ 𝑡 1 , ∗ 𝑅 sh , R , ∗ ( 𝑡 1 , ∗ ) + 𝑅 sh , TR , ∗ GLYPH<12> GLYPH<12> 𝑡 ∗ 𝑡 1 , ∗ 𝑡 1 , ∗ < 𝑡 ∗ ≤ 𝑡 2 , ∗ 𝑅 sh , R , ∗ ( 𝑡 1 , ∗ ) + 𝑅 sh , TR , ∗ GLYPH<12> GLYPH<12> 𝑡 2 , ∗ 𝑡 1 , ∗ + 𝑅 sh , NR , ∗ GLYPH<12> GLYPH<12> 𝑡 ∗ 𝑡 2 , ∗ 𝑡 2 , ∗ < 𝑡 ∗  \nwhere 𝑥 ( 𝑡 )| 𝑏 𝑎 = 𝑥 ( 𝑏 ) -𝑥 ( 𝑎 ) . Substituting the constants, we find:  \n𝑅 sh , ∗ =           𝑅 sh , R , ∗ ( 𝑡 ∗ ) 𝑡 ∗ ≤ 𝑡 1 , ∗ 𝑅 sh , TR , ∗ ( 𝑡 ∗ ) -8 . 176 𝑡 1 , ∗ < 𝑡 ∗ ≤ 𝑡 2 , ∗ 𝑅 sh , NR , ∗ ( 𝑡 ∗ ) + 1 . 8 𝑡 2 , ∗ < 𝑡 ∗ (D26)  \nThis method helps to quickly find the numerical root for 𝑡 in Eq. (D16), which can be easily rewritten as 𝑡 ∗ = 𝑇 ∗ + 𝜇𝑅 ∗ ( 𝑡 ∗ ) , with 𝑇 ∗ = 𝑇 / 𝛼 . When the scale 𝛼 ( 𝜃 ) becomes very small (i.e. for 𝜃 ≫ 𝜃 j and 𝑓 ( 𝜃 ) ≪ 1), then, at later times, 𝑇 ∗ becomes larger than 𝑅 ∗ ( 𝑡 ∗ ) and we can impose 𝑡 ≃ 𝑇 to speed up the calculations.  \nTheshockedmediumischaracterized by its comoving number density 𝑛 ' , internal energy density 𝑒 ' 𝑖 , and magnetic field strength 𝐵 ' . Those are given as follows Sari et al. (1998):  \n𝑛 ' = 4 Γ 𝑛 0 (D27)  \n𝑒 ' 𝑖 = ( Γ -1 ) 𝑛 ' 𝑚 p 𝑐 2 = 4 Γ ( Γ -1 ) 𝑛 0 𝑚 p 𝑐 2 (D28)  \n𝐵 ' = √︃ 8 𝜋𝜀 B 𝑒 ' 𝑖 (D29)  \nwhere 𝜀 B is a parameter that expresses the energy fraction of the shocked matter going into the magnetic field.  \nThe non-thermal electrons generating the synchrotron emission are described by an energy spectrum with a single power-law 𝑝 , and a constant fraction 𝜀 e of the shock energy is assumed to go into these accelerated electrons. The synchrotron spectrum is approximated by a broken power law, with a characteristic frequency 𝜈 ' m and a cooling frequency 𝜈 ' c . These are given by:  \n𝜈 ' m = 3 16 GLYPH<20> 𝜀 e GLYPH<18> 𝑝 -2 𝑝 -1 GLYPH<19> 𝑚 p 𝑚 e ( Γ -1 ) GLYPH<21> 2 𝑞 e 𝐵 ' 𝑚 e 𝑐 (D30)  \n𝜈 ' c = 3 16 GLYPH<20> 3 𝑚 e 𝑐 Γ 4 𝜎 T 𝜀 B 𝑒 ' 𝑖 𝑡 GLYPH<21> 2 𝑞 e 𝐵 ' 𝑚 e 𝑐 (D31)  \nwhere 𝑞 e is the electron charge and 𝜎 T = 8 𝜋𝑞 4 e /( 3 𝑚 2 e 𝑐 4 ) is the Thomson cross section. The emissivity 𝜀 ' 𝜈 ' , in case of slow cooling ( 𝜈 ' m < 𝜈 ' c ) is given by:  \n𝜀 ' 𝜈 ' = 𝜀 ' 𝜈 ' , peak            ( 𝜈 ' / 𝜈 ' m ) 1 / 3 𝜈 ' < 𝜈 ' m ( 𝜈 ' / 𝜈 ' m ) ( 1 -𝑝 )/ 2 𝜈 ' m ≤ 𝜈 ' < 𝜈 ' c ( 𝜈 ' c / 𝜈 ' m ) ( 1 -𝑝 )/ 2 ( 𝜈 ' / 𝜈 ' c ) -𝑝 / 2 𝜈 ' c ≤ 𝜈 ' (D32)  \nwhere 𝜀 ' 𝜈 ' , peak is the peak emissivity (Granot et al. 1999):  \n𝜀 ' 𝜈 ' , peak ≃ 0 . 88 × 256 27 GLYPH<18> 𝑝 -1 3 𝑝 -1 GLYPH<19> 𝑞 2 e 𝑚 𝑒 𝑐 2 𝑛 ' 𝐵 ' . (D33)  \nFinally, the rest-frame frequency 𝜈 ' is connected to the observer frequency 𝜈 by:  \n𝜈 ' = Γ ( 1 -𝛽𝜇 ) 𝜈 (D34)  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 17. Time evolution of the parameters of the generalized Gaussian fit for four jet simulations with 𝜃 j = [ 3 , 5 , 6 . 8 , 10 ] · and 𝐿 j = 10 51 erg / s. Aspecific color and marker identify each run in each plot. The top left panel shows 𝐸 0 multiplied by an integration factor that depends upon 𝜃 c and the shape parameter normalized to the jet energy 𝐸 j , namely 𝑓 ( 𝜃 c , 𝑠 ) = 2 𝜋𝜃 2 c [ 𝑐 1 ( 𝑠 ) -𝑐 2 ( 𝑠 ) 𝜃 2 c ] . The top right panel shows the ratio 𝜃 c / 𝜃 j vs time. Analogously, the bottom left panel shows the shape parameter 𝑠 vs. time, and the bottom right parameter shows the linear correlation between 𝜃 c and the shape parameter.  \n<!-- image -->", 'D.2. Degeneracy relations': "The afterglow light curves are subject to degeneracy relations for the peak position and the peak height, derived from the following scaling laws (Nakar et al. 2002; Nakar & Piran 2020):  \n𝑇 peak ∝ 𝐸 1 / 3 0 𝑛 -1 / 3 0 ( 𝜃 obs -𝜃 j ) 2 (D35)  \nF 𝜈, peak ∝ 𝐸 0 𝑛 𝑝 + 1 4 𝜀 𝑝 -1 e 𝜀 𝑝 + 1 4 B 𝜃 -2 𝑝 obs 𝜈 1 -𝑝 2 𝐷 -2 (D36)  \nfor a fixed observer time 𝑇 . Here, 𝐸 0 = 𝐸 Ω ( 𝜃 = 0 ) . Additionally, the scaling for 𝜃 obs becomes more precise as long as 𝜃 obs ≫ 𝜃 j. We verified the scaling relations by numerically integrating the light curves for several viewing angles, from 𝜃 obs = 14 · to 60 · at steps of 1 · , and for four different opening angles of the jet, i.e. 𝜃 j = [ 3 , 5 , 6 . 8 , 10 ] · . We show the results of the peak time and the peak flux as a function of the viewing angle in Fig. 16. From the fitting lines we see that, for a ratio 𝜒 = 𝜃 obs / 𝜃 j sufficiently larger than unity, the scaling relation is a power law with slope ≃ 2 . 2 (instead of 2 in Eq. (D35))  \nfor 𝑇 peak and a slope of ≃ -2 . 4 𝑝 for F 𝜈, peak (instead of -2 𝑝 from Eq. (D36)). Furthermore, the afterglow peak tends to be wide for 𝜒 ≃ 1 -3 and narrower for 𝜒 ≫ 1 until the shape converges to a certain curve. Combining Eqs. (D35) and (D36), and taking into account the corrections above for the scaling with respect to 𝜃 obs, we get the generalized scaling relations for the peak time and the peak flux:  \n𝑝 + 5 4 log GLYPH<18> ˜ 𝑛 0 𝑛 0 GLYPH<19> + 𝑝 + 1 4 log GLYPH<18> ˜ 𝜀 B 𝜀 B GLYPH<19> + ( 𝑝 -1 ) log GLYPH<18> ˜ 𝜀 e 𝜀 e GLYPH<19> = = ( 3 𝑙 1 + 𝑙 2 𝑝 ) log GLYPH<18> ˜ 𝜃 obs 𝜃 obs GLYPH<19> + 3 𝑙 1 log GLYPH<18> 1 -1 / ˜ 𝜒 1 -1 / 𝜒 GLYPH<19> (D37)  \nwith 𝑙 1 ≃ 2 . 2 and 𝑙 2 ≃ 2 . 4. If the observer's angle does not change, the second term goes to 0, and the previous equation reduces to the scaling relation of Takahashi & Ioka (2021). We remark that, even if the peak time and peak flux coincide for different viewing angles, the previous equation does not define a family of perfectly superposing curves since the af-  \nght curve tends to spread the more the observer gets closer to the jet's opening angle. This is visible from the third panel of Fig. 16 where we measured the time Δ 𝑇 between the peak time 𝑇 peak and the asymptotic decay of the light curve (i.e., the observer time when the flux decreases by roughly 10% from its peak value), which roughly follows the power law scaling of a spreading jet (Granot 2007).", 'D.3. Jets described by a generalized Gaussian and their time evolution': 'The integral of the energy per solid angle of a jet described by Eq. (6) is rather complicated and is analytical only for specific integer or fractional values of 𝑠 ; thus, it is better to approximate the integrated function using a Taylor series expansion of the result around 𝜃 c = 0, as it follows:  \n𝐸 j = 2 𝜋𝐸 0 ∫ 𝜃 max 0 d 𝜃 sin 𝜃 exp GLYPH<20> -1 2 GLYPH<18> 𝜃 𝜃 c GLYPH<19> 𝑠 GLYPH<21> ≃ 2 𝜋𝜃 2 c ( 𝑐 1 -𝑐 2 𝜃 2 c ) 𝐸 0 (D38)  \nVarying 𝑠 and calculating the integral numerically, we can extrapolate the behavior of the coefficients as a function of the shape parameter 𝑠 . The coefficients 𝑐 1 ( 𝑠 ) and 𝑐 2 ( 𝑠 ) have the following approximate forms:  \n𝑐 1 ( 𝑠 ) ≃ 1 2! ( 2! × 4 ) 𝑠 -𝜋 / 2 𝑐 2 ( 𝑠 ) ≃ 1 4! ( 4! × 4 2 ) 𝑠 -𝜋 / 2 . (D39)  \nFor 𝑠 →∞ , the edges of the jet become infinitely sharp. The jet shape degenerates into a Heaviside step function 𝐸 0 Θ ( 𝜃 𝑐 -𝜃 ) , thus the integral of Eq. (D38), dropping a factor 2 𝜋 , will give exactly 1 -cos ( 𝜃 c ) = " ∞ 𝑗 = 1 (-1 ) 𝑗 + 1 1 2 𝑗 ! 𝜃 2 𝑗 c = 1 2! 𝜃 2 c -1 4! 𝜃 4 c + . . . , which are the asymptotic limits of the coefficients 𝑐 1 and 𝑐 2.  \nIn Fig. 17, we show the time evolution of the fitting parameters for four values of the opening angle 𝜃 j = [ 3 , 5 , 6 . 8 , 10 ] · and fixed 𝐿 j = 10 51 erg / s. We took snapshots from 25 ms to 450 ms of continuous jet injection, at Δ 𝑡 = 25 ms from each other. We cut off the rest of the slow ejecta in each snapshot and fitted the on-axis peak with a generalized Gaussian. We used Eqs. (D38) and (D39) to calculate the total energy under the fitting function and compared it to the jet energy for each given time. We see that for all the runs but one ( 𝜃 j = 10 · ), the result matches the expected jet energy very well. The discrepancy at larger angles can be explained considering that there is a non-negligible amount of energy outside the peak, thus giving less than expected energy. The value of 𝜃 c grows in the first 250 ms of injection to roughly stabilize around the value of the jet\'s opening angle 𝜃 j. Also, the shape parameter of the Gaussian seems to converge to a certain value at later times, around 5. However, due to its wide angular spread, the run with 𝜃 j = 10 · seems to jump to higher values (flatter central distribution). Not surprisingly, 𝜃 c and the shape parameters are linearly correlated at later times, as shown in the bottom-right panel of the figure. From these data, we  \ncan conclude that taking the energy distribution 𝐸 Ω too early might strongly affect the final result when the synthetic light curves are calculated since its parameters tend to stabilize at later times.'}