Impact of anisotropic ejecta on jet dynamics and afterglow emission in binary neutron-star mergers

Binary neutron stars mergers widely accepted as potential progenitors of short gamma-ray bursts. After the remnant of the merger has collapsed to a black hole, a jet is powered and may breakout from the the matter expelled during the collision and the subsequent wind emission. The interaction of the jet with the ejecta may affect its dynamics and the resulting electromagnetic counterparts. We here examine how an inhomogeneous and anisotropic distribution of ejecta affects such dynamics, dictating the properties of the jet-ejecta cocoon and of the afterglow radiated by the jet upon deceleration. More specifically, we carry out general-relativistic hydrodynamical simulations of relativistic jets launched within a variety of geometrically inhomogeneous and anisotropic distributions of ejected matter. We find that different anisotropies impact the variance of the afterglow light-curves as a function of the jet luminosity and ejected mass. A considerable amount of the jet energy is deposited in the cocoon through the jet-ejecta interaction with a small but important dependence on the properties of the ejecta. Furthermore, all configurations show a two-component behaviour for the polar structure of the jet, with a narrow core at large energies and Lorentz factors and a shallow segment at high latitudes from the jet axis. Hence, afterglows measured on off-axis lines of sight could be used to deduce the properties of the ejected matter, but also that the latter need to be properly accounted for when modelling the afterglow signal


INTRODUCTION
The association between binary neutron-star (BNS) mergers and short gamma-ray bursts (SGRBs) was recently consolidated in August 2017 via the coincident detection of the gravitational-wave signal from the inspiral of a BNS (GW170817, The LIGO Scientific Collaboration & The Virgo Collaboration 2017) and of a short, hard gamma-ray burst (GRB170817, Goldstein et al. 2017;LIGO Scientific Collaboration et al. 2017).Although within the population of SGRBs observed so far this event is characterised by a very low luminosity (Matsumoto et al. 2019), its multi-messenger nature provides strong evidence in support of the interpretation of BNS mergers as progenitors of SGRBS (Paczynski 1986;Eichler et al. 1989).In the case of GW178717, in particular, the afterglow observations suggest that a regular, bright SGRB would have been observed if it has been seen from an on-axis line of sight (Salafia et al. 2019).
After these signals were first detected, a follow up campaign allowed to localise the source via the detection of the kilonova optical transient (Abbott et al. 2017 and references therein).Spectroscopic observations of this transient suggests the presence of rprocess nucleosynthesis in the merger ejecta (Smartt et al. 2017;Watson et al. 2019).The discovery and long-term sampling of the afterglow signal across the electromagnetic spectrum opened a rich window on the outflow from the merger (Troja et al. 2017;Hallinan et al. 2017;Lyman et al. 2018).The combined photometry and VLBI image of the afterglow gave clear indication of a relativistic jet emerging from the source, being viewed with a significantly off-axis line-of-sight (Mooley et al. 2018a;Ghirlanda et al. 2019).
The multi-wavelength light-curves of GW170817 were characterized by a late-time power-law decay.In addition, observations using very long baseline interferometry revealed an apparent superluminal motion (Mooley et al. 2018a;Lamb et al. 2019;Ghirlanda et al. 2019).This motion helped constrain the range of allowed observational angles to θ obs = 15 ± 1.5 deg at 1-σ and hinted to a very energetic jet with a core energy of Ejet ≈ 10 49−50 erg; other works (Nathanail et al. 2021;Wang et al. 2023) extended the inclination range to even higher values, and close to 30 deg.The quantity of ejecta, inferred from the kilonova transient suggests that the remnant object collapsed at a time of ≈ 1sec after the merger (Gill et al. 2019;Murguia-Berthier et al. 2021;Beniamini et al. 2020b).
Considerable theoretical work has been invested over the years to study the evolution of SGRB jets, starting from the launching of the incipient jet in the immediate environment of the mergers (see, e.g., Aloy et al. 2005;Rezzolla et al. 2011;Murguia-Berthier et al. 2016;Most et al. 2019;Hayashi et al. 2022;Most et al. 2021), continuing with its interaction with the ejecta and breakout (Gottlieb & Globus 2021;Urrutia et al. 2021;Gottlieb & Nakar 2022), and up to the exploration of the various processes that shape the large-scale structure of the jet and that eventually lead to the afterglow signal.Among these processes, the interaction with the merger ejecta is expected to be particularly interesting as the jet can potentially acquire (or modify) its structure, hence determining how the Lorentz factor and energy of the material in the jet at a large distance from the source depend on the latitude with respect to the jet axis.These dependencies on the latitude in the jet are usually captured in two functions Γ∞(θ) and dEjet/dΩ(θ) denoting the Lorentz factor and distribution of total energy of the jet material per unit solid angle at latitude θ in the jet, which is assumed to be axis-symmetric.Together, these functions as referred to as the "jet structure" (see, e.g., Salafia & Ghirlanda 2022, for a review).
Hence, understanding the formation of the jet structure in this chain of events is not only interesting, but also important when making use of the afterglow emission to infer the structure of the relativistic jet and, in turn, to explore the lower-scale physics of jet launching.
This connection between the large-scale structure of the jet and the afterglow signal has been explored a number of times in the literature and employing a variety of methods: • by prescribing a functional form for the structure and by inferring its parameters through light-curve fitting (see, e.g., Gill & Granot 2018;Troja et al. 2019, for the case of GW170817); • by deriving a standard inversion procedure to determine the structure starting from the light-curve (see, e.g., Takahashi & Ioka 2021); • by discovering a set of jet structures all consistent with the data from GW170817 (see, e.g., Takahashi & Ioka 2020); • by studying analytically the diverse light-curve morphologies that can arise from a given form of the jet structure (see, e.g., Beniamini et al. 2020a); • by focusing on the mimicking of jet-structure effects by other jet dynamics, such as its lateral expansion (see, e.g., Lamb et al. 2021).What all of these methods have in common is the realisation that degeneracies exist in the theoretical modelling and hence all the available information needs to be incorporated when constraining the jet structure from afterglow observations.The connection between the BNS merger and the large-scale structure of the jet is tightly linked with the dynamical evolution that takes place in the post-merger.More specifically, the violent merger produces an immediate dynamical ejection of mass (see, e.g., Rezzolla et al. 2010;Sekiguchi et al. 2015;Lehner et al. 2016;Radice et al. 2018;Ruiz et al. 2020), accompanied by secular mechanisms that further eject mass from the system through magnetically-driven or and neutrino-driven winds from the accretion disk and from the remnant before its eventual collapse (Dessart et al. 2009;Siegel et al. 2014;Fernández et al. 2015;Foucart et al. 2016;Fujibayashi et al. 2018).Subsequently, and on much larger scales away from the central object, the jet-ejecta interaction leads to the formation of a jet-cocoon system (Bromberg et al. 2011) and imprints the jet structure upon breakout, before it propagates to larger distances and produces gamma-ray burst prompt and afterglow radiation.General-relativistic simulations have been instrumental in the study of this early phase of the jet dynamics and has allowed, among other things, to study the conditions for the jet to breakout (Duffell et al. 2018;Gottlieb et al. 2022), the effect of baryon mixing between the jet and the cocoon (Gottlieb et al. 2020a;Gottlieb et al. 2020b), the possible formation of a hollowcore structure and its dependence on magnetic fields (Nathanail et al. 2019(Nathanail et al. , 2020)), the role of the structure of the incipient jet in its interaction with the ejecta (Urrutia et al. 2021), or the possible asymmetries and oscillations of the jet imprinted by the interaction with the ejecta due to three-dimensional (Lazzati et al. 2021).Furthermore, by performing simulations of the jet-ejecta interaction and comparing spherically-symmetric ejecta with more realistic ejecta imported from BNS merger simulations, Pavan et al. (2021) have pointed out the considerable impact that even a small anisotropy in the ejecta rest-mass density profile can have on the emerging jet structure.These numerical studies use various prescriptions for the ejecta and incipient jet properties, and include various sets of physical processes in the simulations.
In this work, we extend the exploration of the impact of the jetejecta interactions on the large-scale properties of the jet to the case of an inhomogeneous and anisotropic distribution of matter in the merger ejecta.These radial and angular gradients can be produced by a number of processes, such as, unequal composition and heating rates by nucleosynthetic processes stemming from anisotropic neutrino irradiation (Combi & Siegel 2023), or by sporadic baryon loading of the base of the jet, to name a few (Dionysopoulou et al. 2015;Gottlieb et al. 2019).Furthermore, the different components of the merger ejecta (dynamical ejecta, winds, etc.) are expected to mix and form shocks (Bovard et al. 2017;Radice et al. 2018), resulting in an overall inhomogeneous envelope which the jet will have to penetrate.
In particular, we carry out two-dimensional (2D) generalrelativistic hydrodynamic simulations of relativistic jets propagating through an envelope of ejected material following a variety of prescribed rest-mass density profiles featuring strong gradients in both the radial and polar direction.We study the interaction of the jets for up to 3 s after jet launching, focusing on the effects of these gradients on the jet-cocoon system, on the propagation of the jet and on the large-scale structure adopted by the jet after breakout.Overall, we find that the cocoon is generally robust to different anisotropies in the ejecta and that the more massive the ejecta are relative to the jet energy, the more the jet structure is affected.As a result, a strong observational imprint is in principle present on off-axis afterglow observations of the most massive ejecta.
The structure of the paper is organised as follows.In Sec. 2 we describe the numerical methods employed for both the numerical simulations and the calculation of the afterglow, accompanied by the different setups for the ejected matter and the injection of the jet.We instead report the results of the simulations and afterglow observations in Sec. 3, together with comparison with the observations and how they can be used to infer the properties of the ejected material.The conclusions and prospects for future work are presented in Sec. 4.

Numerical setup
Our simulations are performed employing the numerical code BHAC (Porth et al. 2017;Olivares et al. 2019), which solves the equations of general-relativistic hydrodynamics or magnetohydrodynamics using high-resolution, shock capturing methods in one, two and three spatial dimensions and in a variety of coordinates (we have employed spherical polar coordinates here).Adopting the Schwarzschild solution as a fixed curved background metric, our grid starts from an inner radius of Rin = 10 9 cm and extends to Rout = 6 × 10 10 cm.While for the outer boundary we impose standard outgoing boundary conditions, a more involved prescription is used for the inner ones as long as the jet is injected (see discussion below for when the jet is progressively shut-down).More specifically, on the first radial shell we examine whether the polar angle θ is smaller than the opening angle of our jet θjet.If so, we set an initial Lorentz factor and fix density and pressure in order to match a top-hat jet, for a given luminosity Ljet, and asymptotic Lorentz factor Γ∞.For larger angles, we set a smooth solution that follows the prescriptions detailed below for the ejected matter, that extends from the inner boundary radius rin = Rin = 10 9 cm out to rout = 5 × 10 10 cm.
The plasma is modelled as a simple ideal-gas fluid (Rezzolla & Zanotti 2013) with adiabatic index γ = 4/3 and, as anticipated above, is characterised by a large degree of inhomogeneity and anisotropy, as expected for the matter ejected on dynamical and secular timescales.Because the actual degree depends sensitively on the mass ratio of the binary, the equation of state and the magnetisation, it is difficult to derive from first principles analytic expressions that capture this large diversity.Of course, one could resort to importing the distribution of the ejecta directly from numerical simulations, but this would limit significantly the scenarios considered to those for which numerical simulations are performed and essentially prevent the systematic exploration that wish to carry out in our work.In view of these considerations, we have opted for a simplified prescription of the inhomogeneities and anisotropies in the ejected matter in terms of a reference "homogeneous" and isotropic configuration in which the rest-mass density follows a power-law fall in radius (see, e.g., Nagakura et al. 2014) where ρ0 is a constant.On the other hand, to study the impact of the inhomogeneous and anisotropic ejecta on the dynamics and emission of the jet, the rest-mass distribution is taken to have also a polar dependence of the type where we vary the dimensionless function ξ to explore three different anisotropic configurations.More precisely, we further define (3) and ensure that low-density patches exist in the ejecta by setting the function ϕ(r, θ) to: sin (6θ) sin (6π(r−rin)/rout) ; cos (6θ) sin (6π(r−rin)/rout) ; cos (6θ + 2π(r−rin)/(rout − rin)) sin (6π(r−rin)/rout) .
(4) which we refer to as: "on-axis", "off-axis", and "mixed", respectively; clearly, the "homogeneous" profile (2) corresponds to the case in which ϕ(r, θ) = 0. We note that while the profiles given by the expressions (4) are simplified, their spectral distribution in wavelengths resembles that observed in numerical simulations; in addition, the ability to place in precise locations regions of over-or under-density allows us to ascertain in much more controlled way the impact that inhomogeneities have on the jet propagation and hence on the afterglow signal (see discussion in Sec.3.4).The initial conditions in the ejecta in the three cases on-axis, off-axis or mixed can be found in the three panels of Fig. 1.In all of the cases considered, the constant ρ0 is determine such that the total rest-mass in the ejecta, denoted by Mej, varies in the range Mej/M⊙ ∈ [0.001, 0.01], and corresponds to the range expected from BNS mergers (Rezzolla et al. 2010;Rezzolla et al. 2011;Bauswein et al. 2013;Hotokezaka et al. 2013;Radice et al. 2016;Dietrich et al. 2017;Radice et al. 2018;Fujibayashi et al. 2023).
For a given choice of ejecta distribution and total rest-mass Mej, the jet will sweep on its path is a rest-mass column defined as which does not change by more than 5% across the various profiles and is at least 0.4% of Mej.
By construction, the matter in the ejecta is all gravitationally unbound and the corresponding initial velocity profile is assumed to be isotropic and purely radial, increasing linearly with radius, i.e., where er is the unit-vector in the radial direction.Expression (6) leads to values that are close to those measured in numerical simulations (see, e.g., Baiotti & Rezzolla 2017;Paschalidis 2017;Shibata & Hotokezaka 2019, for some reviews), and within the ranges inferred from the observations (see, e.g., Arcavi et al. 2017;Drout et al. 2017;Shappee et al. 2017;Pian et al. 2017;Kilpatrick et al. 2017;Kasliwal et al. 2017;Kasen et al. 2017;Evans et al. 2017;Smartt et al. 2017).
The jet is assumed to be launched by the "central engine" at the innermost radial shell and with an initial opening angle θjet,i = 0.1.Its structure is assumed to be that of a "top-hat" jet, that is, with sharp cut-offs for all quantities at angles larger than the opening angle.For each simulation, we inject a jet with constant luminosity, Ljet, using three different reference values 10 49 erg/s, 5 × 10 49 erg/s, and 10 50 erg/s.The jet is injected in the grid for a duration of tjet = 1 s and has an initial Lorentz factor Γ0 = 10.We set the asymptotic value of the Lorentz factor to Γ∞ = 100, which is translated to an initial specific enthalpy for the injected plasma of h0 = Γ∞/Γ0 since h∞ = 1.In practice, at every timestep and for the grid cells in the first radial slice within the jet opening angle we set After a time t = tjet, which can be smaller or larger than the breakout time depending on its energy and mass of the ejecta, the central engine is not turned-off abruptly but following an exponential fall-off for a duration of τ rel = 0.1 tjet, which, in turn, produces a sufficiently large decay of the asymptotic Lorentz factor.This smooth turning-off of the jet is implemented to avoid the generation of artificial shocks in the numerical domain, especially in the tail of the jet.After the transition time, we treat the inner radial boundary inside the jet-angle in the same way as we do for the ejected material.In other words, we set the rest-mass density to follow a power-law in time with exponent −6.A summary of the physical properties of the 27 different simulations carried out in our investigations can be found in Table 1.
Finally, our 2D grid has a number of base cells in the radial and polar direction given by Nr × N θ = 1120 × 576 and we employ three refinement levels, thus reaching a peak effective resolution of Nr × N θ = 8960 × 4608.This is a rather high resolution that provides us with sufficiently detailed information about the dynamics of the jet during its propagation and breakout.Furthermore, a discussion of the consistency and robustness of our results when considering different numerical resolutions can be found in Appendix A.

Large-scale jet structure and afterglow emission
In order to determine the astronomical observables that would allow us to probe the jet-ejecta interaction, we calculate the afterglow radiation produced as a result of the jet morphology computed in the various simulations.To this scope, we first need to compute the angular structure in total jet energy dEjet/dΩ(θ) and in Lorentz factor Γ∞(θ), which are extracted from the simulation, at a time of t = 2.25 s after launching of the jet.To compute the energy density e in the inertial frame for a perfect fluid moving with Lorentz factor Γ we can make use of the fact that it possesses a conserved quantity, namely, the Bernoulli constant given by (see, e.g., Rezzolla & Zanotti 2013, Secs. 3.6.2and 11.9.2) Hence, for a perfect fluid obeying an ideal-fluid equation of state with adiabatic index γ = 4/3, expression (8) can be inverted to obtain the energy density of the fluid in the inertial frame To obtain the distribution of jet energy per solid angle at latitude θ, we integrate Eq. ( 9) over the radial direction: which we will refer to as the "E-structure" 1 , where r ej,front is the ejecta front, a choice made to avoid including plasma from unshocked ejecta in the jet structure and is calculated by : where we set vej,max = 0.6 [Eq.( 6)] and tex = 2.25 s is the extraction time of the jet structure in the simulation domain (see also Sec. 3.3).The choice of the definition of the energy structure as 1 Indeed, according to Eq. ( 10), the total energy contained in a 3D volume r ∈ [r ej,front , rout] and (θ, ϕ) ∈ Ω 1 , where Ω 1 is a solid angle, will be: Etot =  1.Properties of the various simulations performed in our analysis, each of which is marked by the corresponding jet luminosity L jet and by the ejected mass M ej .
the energy per unit volume is guided by the fact that this function is constant in a top-hat jet, or generally in an isotropic expansion, and because it is usually considered in prior numerical work and analytical work.In this integral, the choice of bounds in radius r1 and r2 is done to consider only relativistic material in the domain, as explained below.
We define the jet structure in Lorentz factor via the asymptotic Lorentz factor Γ∞, which can be easily calculated using again the Bernoulli constant to obtain: that is, the Lorentz factor that the material would reach if all its internal energy were converted to kinetic energy.Then, we determine the average asymptotic Lorentz factor in the solid angle around latitude θ that we refer to as the "Γ-structure".We note that when performing the volume integrals to determine both the Eand Γ-structures we need to avoid taking into account plasma that is not relativistic a simply part of the moving ejecta.To this scope, we define different criteria to examine whether a cell should be excluded or not from the domain.In essence, they consist of cut-offs in asymptotic velocity, which are considered to only include in the relativistic jet material with Lorentz-factor high enough to accelerate a population of non-thermal particles in the forward shock susceptible to produce synchrotron radiation (Sironi et al. 2013).(see Sec. 3.1 for more details on these Lorentz-factor cutoffs and the robustness of our results with respect to different choices).
The procedure outlined above allows us to determine the structure of the jet, i.e., the dependence of the Lorentz factor and energy of the outflow on the polar angle with respect to the jet axis.Once emerged from the jet, the outflow will then decelerate by colliding with the interstellar medium leading to the afterglow emission.Given the likely achromaticity of afterglows of jets viewed off axis, we restrict our study in the emission in the radio band (Duque et al. 2019;Beniamini et al. 2020a) and model the afterglow light-curves in the radio by following the procedure described by Granot et al. (2002) and Gill & Granot (2018), which we briefly review below.
First, we assume that each angular slice of the relativistic jet evolves independently of the others; this amounts to ignoring the lateral spreading (i.e., considering v θ = 0 = v ϕ ) and tends to an over-estimation of dEjet/dΩ and Γ∞.The interaction of the jet with the ambient medium obviously leads to what is normally referred to as the "forward" shock and which is responsible for dominant part of the emission (see Granot & Sari 2002 for the validity of this assumption).The dynamics of the forward shock proceeds through two different phases in time.First, during the socalled "coasting phase", the shock front propagates with a constant Lorentz factor Γ(r) = Γ0 until a deceleration radius for uniform circum-burst medium, where n is the number density of the ambient medium the jet penetrates after the breakout.We determine the deceleration radius and subsequent deceleration dynamics for all the angular slices according to the Eand Γ-structure as described previously.
Following the case of GW170817 and in line with the low densities expected in media around BNS mergers, we adopted n = 3 × 10 −3 cm −3 with a uniform circum-burst medium for all our afterglow calculations.Once sufficient material has been swept up and energised in the forward shock, the jet decelerates significantly, such that, for any polar angle θ, the Lorentz factor evolution with radius follows (Panaitescu & Kumar 2000): where ξ := r/R d .After the coasting phase, the following relation is recovered for the radial dependence of the Lorentz factor of the forward shock: where s is the power-law index of the circum-burst medium, equal to s = 0 in our case of a uniform medium.This corresponds to the "deceleration phase".To capture the entire evolution of the afterglow light-curve, we consider ξ in the range from r = 0 up to the Newtonian radius where the jet becomes non-relativistic.
For a given observation angle θ obs and a given angular slice of the jet θ, we calculate the afterglow flux assuming synchrotron emission from the forward shock formed as the relativistic jet decelerates (Nakar et al. 2002).According to the semi-analytical solution for the deceleration of the jet where the forward shock forms, the forward shock has a comoving internal energy density of u ′ = (Γ − 1)n ′ mp, where n ′ (r) ≈ 4Γ(r)n is the comoving particle number density in the shock and Γ(r) is the Lorentz factor of the forward shock once it has reached radius r (Blandford & McKee 1976).
The fraction of this energy that is deposited to relativistic electrons is denoted as ϵe and essentially determines the mean Lorentz factor of these electrons in the frame comoving with the shock.We assume that the electrons obey a power-law distribution n ′ (γe) ∝ γ −q e , as it is expected from acceleration processes, and shown by particle-in-cell simulations (see, e.g., Sironi & Spitkovsky 2011;

Meringolo et al. 2023
).The power law index q (> 2) of the nonthermal electrons in the afterglow emission is chosen to be q = 2.2, as suggested by the ratio of afterglow fluxes in various bands in GW170817 (e.g.The LIGO Scientific Collaboration et al. 2017) and other SGRB afterglows (Cenko et al. 2011).Furthermore, in the shocked material, a small-scale turbulent magnetic field is produced and we assume that a fraction ϵ B of the comoving internal energy is transformed into magnetic energy in the shocked matter to power the synchrotron emission, such that the magnetic-field strength is Granot et al. 1999).Hereafter, and unless stated otherwise, we set ϵe = 0.112 and ϵ B = 10 −3 ; these values are in agreement with the posterior bounds given by Ghirlanda et al. (2019).
Once the distribution of non-thermal electrons and the magnetic-field strength is set, we can determine the comoving spectral luminosity L ′ ν ′ due to the synchrotron emission for each angular slice of the jet over the region of deceleration.We use the expression for L ′ ν ′ from Gill & Granot (2018), equations 13 and following.This expression covers all the segments of the standard synchrotron spectrum as a function of the basic physical quantities of the shock that we previously presented, namely, n ′ , u ′ , ϵe, ϵ B , and q.
Finally, we calculate the total flux from the jet at radio frequency ν = 3 GHz for a distant observer at viewing angle θ obs and observing time t by integrating the comoving spectral luminosity over the jet structure on equal-arrival-time surfaces, as in Gill & Granot (2018).In the end, the jet-structure integral is written as an integration over angles where δ D := [Γ(1 − vjet cos χ)] −1 is the so-called Doppler factor (vjet is the jet velocity), so that the two frequencies are related through ν = δ D ν ′ , and χ is the angle between the line of sight and the angular slice at (θ, ϕ) at hand in the integration over the structure.For definiteness, we chose a luminosity distance of d L = 40 Mpc, small enough to allow us to ignore redshift corrections.The choice of this distance is not very important, as we are mostly interested in the morphology of the afterglow light-curves.Given the small external density we consider (n = 3 × 10 −3 cm −3 ), the effects of synchrotron self-absorption can be neglected (Duque et al. 2020).

Jet propagation and breakout
We start by studying the effect of the anisotropies in the rest-mass density distribution on the propagation of the jet in the ejecta.
Figure 2 shows a snapshot of the jet-ejecta evolution at the time of the jet breakout for different ejecta density profiles from the Lj.5.49.Me.010 runs (see Table 1).In particular, Fig. 2  the breakout region).For each panel, in addition, the left part refers to a simulation with an anisotropic and inhomogeneous distribution ("off-axis" on the top row and "on-axis" on the bottom row), while the right part provides the evolution in the reference case of an isotropic but inhomogeneous distribution.
When comparing the different parts of the figure, it becomes apparent that the jet propagation is only mildly affected by the presence of the anisotropic perturbations in the ejecta distribution.In turn, the breakout time and energies of the jet at breakout are rather similar among all cases considered, including the homogeneous one (see discussion below for details).At the same time, differences do appear in the morphology of the matter that interacts with the jet and that is set in turbulent motion by the propagation of the jet.This is more evident in the case of the "off-axis" perturbations, which lead to a significant lateral expansion of the jet as it encounters regions in the ejecta with significantly smaller densities.However, differences are present also for the "on-axis" perturbations but, interestingly, these mostly cancel out.In particular, the acceleration following the interaction with an underdense region along the jet direction are compensated by the deceleration when the jet encounters an overdense region.As a result, the jet breaks-out almost at the same time as in the case of a propagation across an isotropic ejecta distribution (compare the left and right parts of the bottom-left panel of Figs. 2 and 3).
The qualitative description offered by Figs. 2 and 3 can be complemented by a more quantitative description summarised in Fig. 4, which reports, in the left panel, the jet breakout time t jet,bo , and in the middle panel, the opening angle of the structured jet θjet.More specifically, we define t jet,bo as the time when the material shocked by the propagating shock first overtakes the radially expanding ejected matter.
As we will show in more detail in Fig. 7, the Lorentz factor experiences a steep decrease outside of inner region (core) that allows us to robustly define the angular size of the ultra-relativistic core of the jet, denoted by θjet and defined as the smallest angle for which Γ∞ < 10.We extract this data from the profiles presented in Fig. 7.We plot this angle on the middle panel of Fig. 4. We find that this angle varies by at most 1.5 deg when changing the pattern of the perturbations in the ejecta density and, in cases where the ejecta mass is low, this variation is much less than 10%.This suggests that, under the condition of a successful jet breakout, the jet core is largely unaffected by the ejecta profile.
The jet break-out and core opening angles reported in Fig. 4 are shown as function of the three luminosities assumed for the jet (horizontal axis) and of the rest-mass density contained in the ejected material (black, red and blue symbols, respectively).Furthermore, the variance shown with vertical lines across the solid symbols reflects the different geometry of the anisotropies (i.e., "off-axis", "on-axis", or "mixed") and indicates one of the most important results of our work, namely that the impact of the geometric anisotropies on jet break-out and jet core properties is not very large and that these jet properties are mainly set by the energy in the jet and by the mass in the ejecta.
Starting from the left in Fig. 4, it is apparent that the breakout time decreases as the energy of the jet is increased, but also that it increases as the ejected matter is more massive.The interpretation of this behaviour is straightforward: the jet travel-time across the ejecta will depend on its energy and on the matter resistance it will encounter, so that the most massive ejecta may actually prevent the breakout of the weakest jet, while the most powerful jet will cross the least massive ejecta in only 0.4 s.At the same time, for a given mass of the ejecta, the breakout time can even dou-  ble between the least and the most powerful jets.Interestingly, the jet with Ljet,50 = 0.10 erg/s (simulation Lj.1.49.Me.010) is indeed chocked when going across the massive ejecta with Mej = 0.010 M⊙.This behavior is in agreement with the predictions for the "early-breakout" criterion of (Duffell et al. 2018) [Eq.( 20) there] but in tension with the "late-breakout" criterion, which instead predicts a successful jet; as pointed out by Gottlieb & Nakar (2022), this tension may be due to the way that the jet is being "turned off".An abrupt shutdown of a jet, or an extended transition phase will have an impact on the cases of "late-breakout".We also note that the variability of the break-out time induced by the anisotropies is larger for the intermediate mass ejecta (see, e.g., Fig. 4).This behaviour, could be due to a number of potential effects but we conjecture it is the result of the fact that in the case of intermediate mass ejecta the system is intrinsically closer to a near-cancellation of effects, which is not present if the mass of the ejecta is very large or very small.Under these conditions, small fluctuations, even numerical in nature, can lead to large variations in quantities such as the break-out time or the jet opening angle.Also, the 2D nature of our simulations, inevitably implies the presence of the so-called "plug-instability" (Zhang et al. 2004;Mizuta & Ioka 2013;López-Cámara et al. 2013), which is not really an instability but rather the artificial piling up of matter ahead of the propagating jet (in 3D this matter has more opportunities to move in the direction orthogonal to that of the jet) and hence affects in part the dynamics of the jet propagation.Fully three-dimensional simulations have shown that the plug-instability tends to over-estimate the breakout time by introducing a systematic bias in the breakout time (see Appendix of Urrutia et al. 2021), which is likely to affect also our results.However, we expect that the basic phenomenological trend of the breakout time as a function of Ljet and Mej portrayed in the left panel of Fig. 4 to be robust and hence to be present also in fully 3D simulations.
The physically intuitive behaviour of the breakout time is shared also by the jet opening angles θjet, shown in the middle panel of Fig. 4, respectively.In essence, what these panels highlight is that this angle increases as the energy of the jet is increased and decreases as the mass in the ejecta is increased.Once again, this is not surprising as it is natural to expect that more powerful jets will deposit systematically larger amounts of energy in the ejecta, leading to larger opening angles at the breakout.Conversely, a larger ejecta mass will exert more pressure on the propagating jet and close the core angle before breakout.Less obvious is that this behaviour should be essentially linear in Ljet, while it is clearly nonlinear in terms of Mej (see Fig. 4).

Cocoon evolution
As mentioned in the previous section, the propagation of the jet within the ejected envelope leads to a turbulent shear layer that is normally referred to as the "cocoon" for its peculiar shape.At the breakout time, Figs. 2 and 3 illustrate the detailed morphology of the cocoon, which can be readily identified by the steep gradient in density.More precisely, we distinguish the matter in the cocoon from the jet material by following the conditions proposed by Hamidani & Ioka (2023), that is, we mark as filled with cocoonmatter cells for which the following conditions are met: Γ∞ ≤ 10, Γ < 5, and |v θ | > 0, where the last criterion is imposed to avoid the mis-identification of ejected material or interstellar medium (for which |v θ | = 0) with the matter in the cocoon.Clearly, Fig. 2 shows that the cocoon varies at different vertical heights and hence it changes during the evolution.More importantly, it is easy to realise that the different geometric properties of the ejected material do have an impact on the morphology of the cocoon.
During its propagation across the ejected material, the jet will deposit a part of its kinetic energy increasing the internal energy of the ejecta E th .Since the Ljet is the source of energy, we can define the efficiency ϵ with which the jet energy is converted into the thermal energy of the cocoon as where ⟨E th ⟩ is the time-averaged thermal energy deposited in the ejecta from the start and up to the breakout time t jet,bo .The unknown efficiency ϵ (and its time-averaged counterpart ⟨ϵ⟩) takes into account the work done both on the jet head surface via the forward shock and that done in generating a lateral expansion of the cocoon [the (kinetic) energy in the turbulent motion in the cocoon is not included in expression ( 17)].Clearly, the value of such an efficiency will depend on the properties of the ejected material and, in our setup, on the properties of the anisotropies introduced.Using our simulations and the measurements of ⟨E th ⟩, we have estimated Figure 5.The same as in Fig. 4 but for the diagnostics of the thermal energy in the cocoon during the jet evolution.In particular, shown are: the average efficiency of jet-cocoon energy transfer [Eq.( 17), left panel], the time-averaged thermal energy of the cocoon up until the jet breakout (middle panel), and the peak thermal energy of the cocoon during the jet evolution (right panel).In the last two panels, the dashed lines indicate the square-root fit with the jet luminosity.
the average efficiency ⟨ϵ⟩ and summarised the results in the left panel of Fig. 5 for the configurations considered.It is then interesting to note that the time-averaged efficiency ⟨ϵ⟩ decreases linearly with increasing energy in the jet and increases with the mass in the ejecta.Stated differently, more energetic jets are less efficient in converting the kinetic energy into thermal energy, especially if the ejecta are not very massive.This behaviour is rather natural and indeed one would expect a vanishing efficiency for extremely tenuous ejecta and very energetic jets, which would propagate in them essentially ballistically.Note also how the efficiency is also upper-limited by a value of about 1.5%, which corresponds to the configuration Lj.1.49.Me.004 that we will later-on see, is that of a marginally successful breakout2 .Considering the amount of energy that is needed to explain the early blue component of the kilonova of GW170817, Duffell et al. (2018) concluded that the exchange of thermal energy in the jet-cocoon system is subdominant.To validate this conclusion, we report in Fig. 5 the time-averaged thermal energy deposited in the ejecta ⟨E th ⟩ (middle panel) and the peak thermal energy E th,peak (right panel), that is, the largest value of E th recorded for the cocoon during the jet propagation through the ejecta, as a function of the energy in the jet and mass in the ejecta.The behaviour of the two quantities is similar even though they really refer to distinct quantities as ⟨E th ⟩ contains information on the overall "heating" of the ejecta, while E th,peak is mostly indicative of observational properties, especially when the observations are performed with sensitivity-limited instruments.Overall, both panels reveal that the energy deposited in the ejecta to produce the cocoon increases with the mass of the ejecta and the luminosity of the jet, and scales like the square root of the latter, as also anticipated by (Duffell et al. 2018).This is shown by the dashed lines that show a fitting function with power-law index equal to 1/2.
In addition, note how the variance introduced by the different anisotropies plays a very modest role in determining the amount of energy deposited and that the impact is larger for the most energetic jets.Still within the context of the GW170817 event, early emission from the cocoon has been studied (see, e.g., Nakar & Piran 2017;Gottlieb et al. 2018) and associated with early observations of an electromagnetic counterpart.Because of the delay of about 10 hours in the localisation of the source, the cocoon contribution to the emission was likely unobserved.Nonetheless, Piro & Kollmeier (2018) claimed the detection of prompt emission in association with the GW170817 event.More recently, Hamidani & Ioka (2023) have built an analytical model to examine the differences between jets produced in collapsars and in SGRBs environment, the main difference being the static or homologously expanding ejecta.In this way, they found that at most 10% of the cocoon mass can escape from the ejecta, and therefore only a small portion of the cocoon energy can escape and contribute to the early emission.Indeed, as long as the cocoon is buried in the ejecta and these are optically thick, it cannot be observed.However, since the ejecta are expanding and cooling, the optical depth will reach, on a timescale of about one day, values below unity revealing the cocoon.Hence, it is important to compare the timescale for the emergence of the emission from the cocoon with the other timescales, namely t jet,bo or the delay between breakout and the γ-ray emission.
In view of these considerations, and given the wealth of information provided by our simulations, it is interesting to assess what impact the anisotropies in the ejecta have on the amount of energy in the cocoon that is trapped in the ejecta.To this scope, and since all of the ejecta is already unbound, we need to mark the part of the "untrapped" cocoon as the material of the cocoon that is moving faster than the ejecta; since the velocity in the ejecta increases linearly with radius [Eq.( 6)] and reaches a maximum of v ej,edge at the forward edge, matter in the untrapped cocoon should satisfy the additional condition where the terminal velocity is calculated as vcoc,∞ = hcocvcoc.For this matter, we then measure the corresponding energy polar distribution dEcoc/dΩ(θ) of the untrapped cocoon in analogy with what was done for the energy in the jet [Eq.( 10)].Also, we compute the total energy in the untrapped cocoon Ecoc as the volume integral of the energy density for matter satisfying the condition (18).The results are summarised in Fig. 6, whose left panel reports the angular distribution of the energy of the "untrapped" cocoon for three different cases of ejecta, namely, "on-axis", "off-axis" and "homogeneous".Note how the polar distribution in the co- Right panel: total energy in the untrapped cocoon shown as a function of the jet luminosity and masses in the ejecta (different coloured symbols).Also in this case, the vertical bars mark the variation that comes from the different density profiles.
coon energy is rather similar in the case of the homogeneous and on-axis ejecta, while it differs in the case of off-axis anisotropies, where more energy is concentrated at higher latitudes.While this may appear surprising, it reflects the analogy already remarked between the propagation across homogeneous and on-axis perturbations, where the sequence of over-and under-densities along the jet path compensate and cancel out, hence giving a very similar cocoon and breakout time as in the homogeneous case (cf., Figs. 2 and 3).The right panel of Fig. 6, on the other hand, shows Ecoc as a function of the two main parameters of our setup, i.e., Ljet and Mej, and the corresponding variance with the initial anisotropies of the ejecta shown again as vertical bars.Rather unsurprisingly, the energy in the cocoon that is able to escape grows with both the energy in the jet and with the mass in the ejecta.In turn, this implies that even if parts of the ejecta are energised more significantly because of their lower density, follow-up interactions with a denser medium do not produce significant changes on the overall dynamics of the cocoon.

Large-scale jet polar structure
Using the procedure described in Sec.2.2, we next discuss the large-scale polar structure of the jet.For the sake of brevity, we introduce the short-hand notation EΩ(θ) for the distribution of jet energy per unit solid angle dEjet/dΩ(θ).In particular, for all of the runs performed Fig. 7 reports the polar dependence of the jet Lorentz factor Γ∞(θ) (red curves) and energy EΩ(θ) (blue curves).More specifically, for each pair of Ljet, Mej, the solid lines in the various panels report respectively the average EΩ(θ) and Γ∞(θ) across the three perturbation distributions ("on-axis", "off-axis", and "mixed"); the shaded areas, on the other hand, span the ranges ∆EΩ(θ) and ∆Γ∞(θ) among the three distributions, as due to the different choices in profiles of the density in the inhomogeneous ejecta.Note that in each panel of Fig. 7 we present the angular profiles assuming two different cuts in velocity when computing the energy and Lorentz-factor integrals [Eqs.( 10) and ( 13)].The first cut (solid lines), aims at representing both the cocoon and the jet and corresponds to matter with a velocity v∞ > 0.3.The second cut (dotted line), shows instead only the mildly-relativistic to relativistic outflow and corresponding to a cut Γ∞v∞ ≥ 3.This cut is meant to refer to material that is expected to contribute the most to the afterglow radiation through the forward shock.Indeed, below this velocity, the efficiency of the forward shock in accelerating the population of non-thermal electrons necessary for the synchrotron emission from the shocked region, is considerably smaller (see, e.g., Morlino & Caprioli 2012;Crumley et al. 2019).
Figure 7 reveals that, except for the particular case of Lj.1.49.Me.004, which we discuss separately later, all the simulations show a similar behaviour for the polar structure of the jet.In particular, the polar profiles of the Lorentz factor consistently show a "core" with large (i.e., Γ∞ ∼ 60 − 90) and near-constant values of Γ∞ within an angle of θ ≈ 5 deg, which is closely related to the initial jet opening angle θjet,i.The Lorentz factor then exhibits a sharp decrease down to Γ∞ ≲ 10 for θ ≳ 7 deg; this sharp decrease, which we have indicated with θjet in Fig. 4, can be well described by a power-law dependence Γ∞(θ) ∝ θ −β with β > 7 for all configurations considered.
Similarly, the polar profiles of the energy also consistently exhibit a core region, with an energy level that is largely independent of the choice of the perturbation pattern in the ejecta.Outside of this core, the energy polar structures then present a decreasing segment where, just like for the Lorentz factor, the profiles are largely independent of the density anisotropies in the ejecta.However, in contrast to the Γ-structure, the E-structure does not have a sharp fall-off for θ ≳ θjet,i, but a much slower decay ranging from a flat distribution (e.g., Lj.5.49.Me.001) to slopes with a power-law index β ∼ 3 (e.g., Lj.5.49.Me.004) up to large latitudes up to θ ≈ 25 deg.This region of slower decay contains a significant reservoir of energy and, remarkably does depend on the actual pattern of anisotropies in the ejecta and, as we will show below, strongly impacts the afterglow emission, most notably at early times.Figure 7. Large-scale polar structures at t = 2.25 s of the nine jets simulated and shown in terms of the total energy (blue lines) and of the Lorentz factor (red lines).Each panel is marked by the name of the corresponding case reported in Table 1 and the top-right panel refers to the case when the jet is unable to break out from the ejecta (i.e., "chocked jet").Solid lines refer to polar structures computed considering the entire plasma outside the ejecta (i.e., "jet and cocoon"), while the dashed lines refer to structures after imposing a a velocity cut (i.e., u∞ > 3).In all cases, the blue/red shaded areas report the variance due to the different types of anisotropy in the ejecta.Shown with a black solid vertical line is the transition angle θtrans (cf., Fig. 4).
To discuss this original feature of the E-structures, we introduce the angle θtrans, defined as the transition angle from the steep segment to the shallow segment in the energy structure of the jet.To capture this transition from the ejecta-pattern-independent to the pattern-dependent segment of the structure we define θtrans as the angle marking a relative difference larger than 50%, between the average structure and the structure range i.e., which depends on the injected jet energy and on the mass of the ejecta.This angle is reported as a vertical black solid for all cases shown in Fig. 7 3 .In addition, in the right panel of Fig. 4, we report θtrans as a function of Ljet for different values of Mej.We find that, analogously to θjet, θtrans increases with Ljet and decreases with Mej, underlining the larger or lesser capacity for the jet to expand without being affected by the anisotropies in the ejecta.Before concluding this section, we concentrate on the case of Lj.1.49.Me.004, which, as already anticipated, deviates from this general trend.While this configuration does produce a largescale jet with a significant Lorentz factor in the core (on average in the core, Γ∞ ≳ 30), there are some density anisotropies in the ejecta that do not lead to an ultra-relativistic jet emerging on a large scale (indeed Γ∞ ≲ 10 in the core for a "off-axis" anisotropic distribution).In addition, the average Lorentz factor in the core is significantly lower than that of the injected jet, and is very low with respect to constraints on SGRB jets.Considering that a similar configuration with a more massive set of ejecta, i.e., Lj.1.49.Me.010, does not lead to a successful breakout but to a chocked jet, we interpret the phenomenology associated with Lj.1.49.Me.004 as to that of a marginally successful breakout.In this configuration, therefore, the energy injected is only marginally sufficient to produce a successful motion across the ejecta that, when slightly more massive, are instead capable of preventing the progression of the jet.We can use these two cases, namely Lj.1.49.Me.004 and Lj.1.49.Me.010, to set a minimum threshold for the "critical efficiency" necessary for a successful breakout given an ejecta mass and a jet propagating at a given luminosity.In particular, using the data for Lj.1.49.Me.010 we determine such an efficiency to be The robustness of this result is provided by the fact that all configurations have efficiencies larger than ηcrit and the same behaviour discussed so far is measured when using other extraction times, e.g., tex = 2.0, 2.2 s, or when varying the threshold for our velocity cut, e.g., Γ∞v∞ = 2.0, 2.5.Finally, the estimate in Eq. ( 20) is in reasonable agreement with similar but distinct estimates (Duffell et al. 2018).

Afterglow predictions and observations
Starting from the jet structures shown in Fig. 7, we can produce the radio afterglow light-curves emission following the procedure discussed in Sec.2.2.These light-curves are reported in Fig. 8 for three different viewing angles θ obs = 4 (magenta lines), 20 (orange lines) and 27 deg (green lines), representing, respectively lines of sight that are either "on-axis" (magenta line) or "off-axis" (orange and green lines); the shading around these lines follows the convention already employed in Fig. 7, namely, marking the variance when considering different anisotropies in the ejecta.Postponing again the discussion of the marginally successful breakout Lj.1.49.Me.004, we note that the afterglows for an observer whose line of sight is on-axis, exhibit a short brightening phase followed by a monotonous power-law decay with a by a monotonous power-law decay with a temporal slope of Fν ∝ t −1.8 , which coincides with the theoretical prediction Fν ∝ t −3(q−1)/2 for the afterglow from a non-expanding jet observed below the synchrotron cooling break, where q = 2.2 is the spectral index of the population of non-thermal electrons in the forward shock (see Sec. 2.2) (Sari et al. 1998;Lamb et al. 2018).On these lines of sight, the afterglow is fully dominated by the emission from the core jet which, as discussed in the previous section, does not depend on the type of anisotropies in the ejecta.On the other hand, for observers with off-axis lines of sight, the afterglows exhibit light-curves as first revealed by GRB170817A, that is, with a long-lasting brightening phase lasting for up to hundreds of days, before reaching a broad peak and connecting to the power-law decay seen by on-axis observers (Mooley et al. 2018c;Gill & Granot 2018).This general behaviour is due to the delayed opening of the beaming cone of the jet core as a result of the large difference between the line of sight and the direction of propagation of the jet.Indeed, given the estimates of the core opening angle θjet ∼ 5 deg (Fig. 4), the off-axis afterglows of Fig. 8 have θ obs /θjet > 4, well inside the regime of significantly off-axis lines of sight.On these off-axis lines of sight, and as time proceeds, the detected radiation is dominated progressively by decelerating material which is moving with angles from θ ≲ θ obs down to θ ∼ 0, i.e., the jet core.
As discussed above, the jet structures at high latitudes are most affected by the ejecta density perturbations, explaining the variation (shown by the shading in Fig. 8) in the afterglow light-curves at early times.As the jet propagates, the variance in the light-curves reduces as the emission becomes dominated by the core structure, which is insensitive to the ejecta pattern.
Figure 8 is also useful to appreciate the impact that the different anisotropies have on the variance of the light-curves as a function of the jet luminosity and ejected mass.Recalling that the smaller the impact, the smaller the variance in the light-curvesand hence the width of the shaded regions -it is apparent that the variance is reduced as the injected energy in the jet is increased Ljet (i.e., from top row to bottom row), but it increases with the mass in the ejecta Mej (i.e., from left column to right column).As a result, the impact of the anisotropies in the ejecta effectively increases with the ratio Mej/Ljet.It is also particularly interesting to note the special behaviour of the afterglow light-curves for the case Lj.5.49.Me.010 (middle row, right column), which exhibit two local maxima.The origin of these two maxima in the afterglows is probably to be found in the two local maxima that are present in the energy distribution for this configuration with a highmass ejecta (cf., panel in middle row, right column of Fig. 7) and, in particular to the significant amount of energy present at highlatitudes, i.e., at θ ≃ 15 deg ≫ θtrans ∼ 3 deg.The presence of high-Lorentz factor material at these high latitudes also plays a role in this morphology of the afterglow light curves.In particular, it is reasonable that first (earlier) of these peaks is produced exactly by the shallow portion of the jet energy structure which dominates at early times and progressively decays.A few tens or hundreds of days later depending on the observer line of sight, the emission is instead dominated by the steep jet core and this leads to the second and more traditional emission observed also in the other configurations.
Interestingly, Beniamini et al. (2020a) have analysed in which conditions a jet structure with a single power-law in energy and Lorentz factor can give rise to a double-peaked afterglow, just as observed for our Lj.5.49.Me.010 system.They found that such afterglows are a generic feature of light-curves from off-axis lines of sight resulting from the interplay between the progressive debeaming of the core and the deceleration of the material on the line of sight.While this is indeed a possibility, our results here show that double-peaks in the light-curves can also arise from a simple angular exploration of the energy structure of the jet, provided it has two portions, with a steep inner core and a shallow exterior with a local maximum at large latitudes.Hence, the system Lj.5.49.Me.010 represents a further illustration of the degeneracies that can be present when inferring jet structures from afterglow light-curves; we will come back to this point in Sec.3.5.
Finally, the presence of an extended shallow segment in the jet structure for systems with moderate Ljet/Mej ratios results in long-lasting phases of slow brightening of the afterglow, as observed in the case of GRB170817A, as the emission is dominated by material progressively ranging up the shallow segment (see.e.g., configurations Lj.5.49.Me.004,Lj.1.49.Me.004).Indeed, the very luminous jets emerge with a steep structure (e.g., Figure 8. Afterglow emission of the nine jets simulated; each panel is marked by the name of the corresponding case reported in Table 1 and the top-right panel refers to the chocked jet.For each panel, the different lines refer to the different viewing angles of the observers the corresponding coloured shaded areas report the variance due to the different types of anisotropy in the ejecta.Note that in the case of model Lj.1.49.Me.004 the light-curve is double peaked as a result of complex distribution of the jet energy in the shallow segment (see main text for a discussion).
Lj.1.50.Me.001), and display a rapidly brightening afterglow, close to Fν ∝ t 3 which is the case for top-hat (Nakar et al. 2002).On the contrary, systems such as Lj.1.50.Me.010 display a slow phase before the peak, signalling the shallow segment.Furthermore, in the specific case of Lj.1.49.Me.004, we have shown in Sec.3.3 that the relativistic nature of the jet after breakout was very sensitive to the perturbation pattern.The structure of the outflow after a marginally successful or a failed breakout could then be better described as a nearly spherically symmetric outflow with a radial stratification in velocity (Kasliwal et al. 2017;Gottlieb et al. 2018), rather than as a jet with angular dependence of energy and Lorentz factor.In this case, the afterglow radiation is better captured by a model of continuous energy injection in the mildly relativistic forward shock, as the early afterglow of GRB170817 suggested (Lazzati et al. 2018;D'Avanzo et al. 2018), rather than the ultra-relativistic structured jet afterglow model used here.
A couple of remarks are useful before concluding this section.First, Fig. 7 shows that the central regions of the jet structure for case Lj.1.49.Me.004 can contain material at the edge of our velocity cut considered for the structure, at Γ∞v∞ ∼ 3.This material, in turn, corresponds to the lower bound of the shaded area for the afterglow light-curve in Fig. 8 relative to an on-axis observer.This material, however, is likely not to be able to accelerate electrons to the level of efficiency of ϵe = 0.1 considered in this work; as a result, the lower bound on the afterglow light-curve for this system is likely overestimated.Furthermore, given the high sensitivity of the outcome of this system on the ejecta anisotropies, it is reasonable to consider the shading in the afterglow of in Fig. 8 as extending to much lower fluxes.Second, for the configuration Lj.5.49.Me.004, we examine how different cases match the radio observations of the afterglow of GW170817 (data collected from Margutti et al. 2017;Alexander et al. 2017;Hallinan et al. 2017;Mooley et al. 2018b;Dobie et al. 2018;Nynka et al. 2018;Troja et al. 2018;Lamb et al. 2019;Fong et al. 2019;Hajela et al. 2019;Makhathini et al. 2021).More specifically, in Fig. 9 we overlay the data from GW170817 (black circles) with the afterglow light-curves determined from the jet structures found in the case of either homogeneous ejecta (blue solid line) or for anisotropic ejecta that are either on-axis (black solid line) or off-axis (red solid line).To this scope, we used the jet polar-structure profiles obtained by applying the velocity cut Γ∞v∞ > 3, and afterglow parameters n = 5 × 10 −4 cm −3 , ϵe = 10 −0.32 , and ϵ B = 10 −3 , while the viewing angle is assumed to be that of an off-axis observer with θ obs = 21 deg.The values of the afterglow parameters n, ϵe and ϵB are different than those adopted in Fig. 8, as they were determined to best-fit the GW170817 data assuming the structure from the Lj.5.49.Me.004.
The figure clearly illustrates that the afterglows differ mostly at earlier times (a result we have already anticipated) and that these differences are rooted both in the anisotropies of the ejecta and in the emission from high latitudes.In view of these results, and while there may be some uncertainty in the structure at these high angles, it is clear that earlier photometry data would be the sharpest way to study jet-ejecta interaction through its signature in the afterglow (see also Sec. 3.5 for a discussion of afterglows to infer jet-ejecta interaction physics and the particular case of GW170817).

On the use of afterglows to study jet structures
The afterglows of relativistic outflows can obviously be used to infer the structure of these jets.This structure can, in turn, provide insight on the physics at smaller scales, on the launching of the jet by the central engine and its interaction with the merger ejecta.Overall, our simulations show that the interaction of the jet with the ejecta, quite independently of their geometric anisotropic properties, leads to an energy polar structure of the jet that is characterised by a steep segment followed by a shallow segment can arise from jet-ejecta interaction (cf., Fig. 7), and that the extent of the shallow segment depends on the strength of this interaction as measured by the ratio of the ejecta mass and of the jet luminosity, Mej/Ljet ratio.Interestingly, to our knowledge, this double steep-shallow-segment jet structure in energy has not been found in previous study of jet-ejecta interaction, and has not been considered as a functional form for the fitting of afterglow data.It does, however, appear as a possible jet structure reconstructed directly from the afterglow data of GRB170817A (Takahashi & Ioka 2021).In particular, as shown in Fig. 8, this structure can give rise to a variety of afterglow morphologies for off-axis observers, from a fast-rise-fast-decline for very small interactions resulting in near top-hat structures (e.g., Lj.1.50.Me.001) to slow-rise-fast-decline for moderate interaction resulting in a reduced shallow segment (e.g., Lj.1.50.Me.010) and a doublepeaked afterglow for significant interaction resulting in a large energy reservoir at high-latitudes (e.g., Lj.5.49.Me.010).Analytical studies have shown that this variety of light-curve morphologies can also arise from simple power-law structures (Beniamini et al. 2020b), underlining once more the large degeneracy (and risks) in reconstructing structure from afterglow light-curves.As an additional caveat, our simulations also highlight that differences in the distribution of the rest-mass density of the ejecta can have a strong impact on the large-scale structure of the jet and of the cocoon.This makes the inferring of the ejecta properties more challenging even when assuming a perfect knowledge of the jet structure.At the same time, the simulations also show that the impact of the ejecta properties depends systematically on both the jet luminosity Ljet to the ejecta mass Mej, and has a regular behaviour with their ratio.
We should also remark that the effect of the jet energy structure has also been recently considered on slightly misaligned lines of sight, i.e., with θ obs ≲ 2 × θjet, in relation to the origin of plateaus (Eichler & Granot 2006;Beniamini et al. 2020a) or of flares in SGRB afterglows (Duque et al. 2022).The models explored in these studies generally require a steep jet structure in this near-core range of angles, in apparent tension with the shallow structure on large angles revealed in the jet of GRB170817A when fit with single power-law structures (Ghirlanda et al. 2019;Gill & Granot 2018).However, the double steep-shallow structure revealed by our simulations naturally reconciles these two observational properties, with the steep segment allowing for these plateaus and flare models on slightly misaligned lines of sight and for a GRB170817A-like afterglows on significantly misaligned ones.
Finally, we note that in the case of GRB170817, the ejecta mass inferred from the kilonova observations is ≳ 0.01M⊙ (Villar et al. 2017;Metzger et al. 2018;Coughlin et al. 2019) and the iso-equivalent energy of the core jet inferred from the afterglow photometry and image is of ∼ 10 52 erg, within a core angle of θjet ∼ 6 deg (Ghirlanda et al. 2019).Assuming a jet duration of ∆t = 1 s and negligible losses to SGRB prompt emission and jetejecta interactions in the core jet, this would translate into a jet luminosity of Ljet ∼ 5 × 10 49 erg/s.This places GRB170817 in the regime of our Lj.5.49.Me.010 simulation, where we expect significant interaction of the jet with the ejecta and a strong effect on the high-latitude energy structure of the jet.On the one hand, this would be compatible with the GRB170817 afterglow and its extended slowly-increasing phase (Mooley et al. 2018c) as due, in our simulations, to a strong jet-ejecta interaction after launch.On the other hand, this conclusion would cast doubts about the possibility of using the afterglow in infer ejecta properties, since in this regime we have shown the sensitivity to ejecta details.
In summary, our simulations show that a great diversity of afterglow light-curve morphologies can arise due to the jet-ejecta interactions.This suggests that the afterglow of GW170817 could not be a standard for BNS merger afterglows, as future GW follow-up observations could reveal.As a result, our systematic investigation reveals that the role of the ejecta in shaping this afterglow provides an additional caveat on inference of the jet structure from the afterglow observations, which fortunately is less important for lower Mej/Ljet ratios.Indeed, if this ratio could be determined to be small using a independent method (e.g., based on the kilonova signal), then reconstructing the jet from the afterglow data will likely be more robust and carry less uncertainty.

CONCLUSIONS
Using general-relativistic hydrodynamical simulations, we have examined the propagation of an ultra-relativistic jet through an envelope of matter ejected by a BNS merger having a variety of anisotropic distributions of rest-mass.Overall, our results can be summarised as follows.
• The geometry and degree of anisotropy of the ejecta impacts the propagation of the jet and influences its breakout time as well as the energy in the cocoon produced.
• Contrary to expectations, the propagation of a jet encountering several anisotropies along its path (i.e., on-axis perturbations) only marginally affects its breakout time, which is comparable with that in the absence of anisotropies.
• The breakout time decreases as the energy of the jet is increased and increases as the ejected matter is more massive.This is simply because the jet travel-time across the ejecta will depend on its energy and on the matter resistance it will encounter.As a result, the most massive ejecta actually prevent the breakout of the weakest jet.
• The time-averaged efficiency ⟨ϵ⟩ in the conversion of the jet kinetic energy into the thermal energy of the cocoon decreases linearly with increasing energy in the jet and with the mass in the ejecta.This behaviour reflects the fact that more energetic jets are less efficient in converting the kinetic energy into thermal energy, especially if the ejecta are not very massive.
• In agreement with Duffell et al. (2018), ⟨ϵ⟩, but also the average thermal energy ⟨E th ⟩ and the peak internal energy E th,peak , show a square-root scaling with the jet energy and a linear scaling with the mass in the ejecta.The variance produced by the different types of anisotropies is marginal and most important at the largest jet luminosities.
• The injection of internal energy into the matter in the cocoon by the propagating jet is sufficiently large to allow for the escape of this energy from the ejecta.The amount of energy in the cocoon that is able to escape Ecoc grows with both the energy in the jet and with the mass in the ejecta.
• Given an ejecta mass and a jet propagating at a given luminosity, a critical efficiency of ηcrit ∼ Ljet,50 (t jet,bo /1 s)(Mej/10 −3 M⊙) −1 ∼ 0.03 exists for a successful jet breakout after a time t jet,bo .
• Essentially all configurations considered show a similar twocomponent behaviour for the polar structure of the jet.In particular, the polar profiles of the Lorentz factor consistently show a "core" with Γ∞ ∼ 60 − 90 within an angle of θ ≈ 5 deg.The Lorentz factor then exhibits a sharp decrease down to Γ∞ ≲ 10 for θ ≳ 7 deg with a power-law dependence Γ∞(θ) ∝ θ −β with β > 7. The polar structure of the energy distribution, on the other hand, does not have a sharp fall-off but a much slower decay ranging from a flat distribution to power-law decays with index β ∼ 3 up to large latitudes up to θ ≈ 25 deg.
• Within the two-components structure of the jet, the core is essentially unaffected by the anisotropies in the ejecta, while the shallow segment displays the greatest variance, which is stronger for lower ratios of Ljet/Mej.In the case of the marginally successful jet, on the other hand, the two-components structure is essentially not present and the polar distributions of energy and Lorentz factor are far more irregular.
• Different anisotropies impact the variance of the afterglow light-curves as a function of the jet luminosity and ejected mass.In particular, the variance is reduced as the injected energy in the jet is increased Ljet, but it increases with the mass in the ejecta Mej.As a result, the impact of the anisotropies in the ejecta effectively increases with the ratio Mej/Ljet.
• The anisotropies in the afterglow observations are most dominant in the first tens of days, for off-axis observations of systems with generally high ejecta mass.In the case of BNS mergers, strong candidates for detection of these effects, are the systems with low mass ratio q := M2/M1, where M1 is the mass of the primary and M2 that of the secondary.The reason for this is that these systems generally produce larger amounts of dynamically ejected matter.
• Double-peaks in the light-curves can also arise from the polar structure of the energy distribution of the jet, provided it has two portions, with a steep inner core and a shallow exterior with a local maximum at large latitudes.This represents a further illustration of the degeneracies that can be present when in inferring jet structures from afterglow light-curves.
• The presence of an extended shallow segment in the jet structure for systems with moderate Ljet/Mej ratios results in longlasting phases of slow brightening of the afterglow as the emission is dominated by material progressively ranging up the shallow segment (cf., GRB170817A).
• In the case of GRB170817A, assuming a jet duration of ∆t = 1 s and negligible losses to GRB prompt emission and jet-ejecta interactions in the core jet, leads to a jet luminosity of Ljet ∼ 5 × 10 49 erg/s.This is in the regime where we expect significant interaction of the jet with the ejecta and a strong effect on the high-latitude energy structure of the jet.
While the present work arguably represents the most extended and systematic investigation of the interaction of a relativistic jet with the ejecta from a BNS merger with non-smooth features in the density distribution, it can be improved in at least three different directions.First, by including the presence of a magnetic field of different topologies and strengths (Nathanail et al. 2020(Nathanail et al. , 2021)), which could influence the propagation of the jet across the ejecta and that, in the presence of turbulence-driven reconnection, could lead to an increase of the internal energy of the cocoon.Second, by extending to three the number of spatial dimensions, which would improve the description of the dynamics of the jet head and remove some of the artefacts discussed here.Finally, by moving away from the simple equation of state considered in our work and by considering instead the same equations of state describing the neutronstar matter in full numerical-relativity calculations.This final step, which is less trivial than it may appear given the complexity of the equations of state employed in BNS mergers, will finally allow for a consistent and faithful import of realistic ejecta distributions from merger simulations.
As a concluding remark, we note how future gravitationalwave observing runs will see more afterglows of jets with significantly off-axis lines of sight, both with current facilities (e.g., Duque et al. 2019;Colombo et al. 2022) and with next-generation instruments (Ronchini et al. 2022), thus allowing for deeper in-sights into the jet structure.In view of our results, the combination of this data with kilonovae signals (see, e.g., Mochkovitch et al. 2021) could provide a very powerful tool to understand the properties of the sources.Indeed, if dim kilonova signals were to be more often associated with rapidly-increasing afterglows, they could provide the signature of a small ejecta mass weakly imprinting the incipient jet, which therefore retained its steep initial structure as described here.In this case, the incipient jet structures could be considered a faithful probe of the jet-launching physics at the scales of the central engine.If, on the contrary, dim kilonovae were to be most regularly associated with slowly-increasing afterglows, they could hint to the fact that the associated shallow jet structures are acquired directly at jet launch.Finally, the inclusion of gravitational-wave constraints on the progenitor binaries can help bridge the gap between the kilonova and afterglow jet signals with the theoretical predictions obtained by numerical simulations on the ability of magnetised merging BNSs to launch relativistic jets (see, e.g., Rezzolla et al. 2011;Hotokezaka et al. 2013;Palenzuela et al. 2015;Ruiz et al. 2016).Me.010 at t = 2.25 s after jet launching.The solid and dashed lines refer to simulations performed with a numerical resolution that is either "high" or "standard", respectively.Clearly, the qualitative features of the polar structures are robust and no new features are introduced when increasing the resolution.

Figure 1 .
Figure 1.Initial rest-mass density distributions of the ejecta considered in our set of simulations.Note that the distributions are both inhomogeneous and anisotropic, with voids that are either "on-axis" (left panel), " off-axis" (middle panel), or "mixed" (right panel) [see Eq. (2) for details].

Figure 2 .
Figure 2. Comparison of jet-breakout dynamics for the case Lj.5.49.Me.010 (i.e., L jet = 5 × 10 49 erg/s and M ej = 0.10 M ⊙ ).The panels in the left column show the rest-mass density, while those in the right column report the distributions of the Lorentz factor.To highlight the contrast among different distributions, the left part of each panel reports the anisotropic distribution of rest-mass density ("off-axis" in the top row and "on-axis" in the bottom row), while the right part the dynamics across a homogeneous distribution.All panels refer to a time t = 0.7 s since the launch of the jet.

Figure 3 .
Figure 3.The same as in Fig. 2 but on a scale highlighting the breakout region.

Figure 4 .
Figure 4. Left panel: Jet-breakout time from the ejecta as a function of the jet luminosity (horizontal axis) and of the mass in the ejecta (different coloured lines and symbols).Middle panel: The same as in the left panel but for the jet opening angle.Right panel: The same as in the left panel but for the transition angle from the steep to the shallow segment in the energy polar structure of the jet [see Eq. (19) for details].In all panels, the vertical bars mark the variation that comes from the different rest-mass density distributions. 50.Me.01Lj.1.50.Me.01

Figure 9 .
Figure 9.Comparison of the radio afterglow from GW170817 (black circles with corresponding error bars) with the afterglow light-curves relative to model Lj.5.49.Me.004 for an observer with viewing angle θ obs = 21 deg.The different lines refer to the different anisotropies in the ejecta (red line for off-axis, blue line for on-axis and black line for the homogeneous distribution).The data is taken at a frequency of ν = 6 GHz and up until 940 d after the burst.
Figure A1.The same as in Fig.7for model Lj.1.50.Me.010 at t = 2.25 s after jet launching.The solid and dashed lines refer to simulations performed with a numerical resolution that is either "high" or "standard", respectively.Clearly, the qualitative features of the polar structures are robust and no new features are introduced when increasing the resolution.