Two steps forward and one step sideways: the propagation of relativistic jets in realistic binary neutron star merger ejecta

The association of GRB170817A with GW170817 has confirmed the long-standing hypothesis that binary neutron star (BNS) mergers are the progenitors of at least some short gamma-ray bursts (SGRBs). This connection has ushered in an era in which broadband observations of SGRBs, together with measurements of the time delay between the gravitational waves and the electromagnetic radiation, allow to probe the properties of the emitting outflow and its engine to an unprecedented detail. Since the structure of the radiating outflow is molded by the interaction of a relativistic jet with the binary ejecta, it is of paramount importance to study the system in a realistic setting. Here we present a three-dimensional hydrodynamic simulation of a relativistic jet propagating in the ejecta of a BNS merger, which were computed with a general relativistic magnetohydrodynamic simulation. We find that the jet's centroid oscillates around the axis of the system, due to inhomogeneities encountered in the propagation. These oscillations allow the jet to find the path of least resistance and travel faster than an identical jet in smooth ejecta. In our setup the breakout time is ~0.6 sec, comparable to the expected central engine duration in SGRBs and possibly a non-negligible fraction of the total delay between the gravitational and gamma-ray signals. Our simulation also shows that energy is carried in roughly equal amounts by the jet and by the cocoon, and that about 20 per cent of the injected energy is transferred to the ejecta via mechanical work.


INTRODUCTION
The simultaneous detection of GW170817 and GRB170817A (Abbott et al. 2017a) has been a milestone event for many fields, from constraining the difference between the speed of gravitational waves and that of light (Abbott et al. 2017b), to the confirmation of the production of very heavy elements in a binary neutron (NS) star merger (Abbott et al. 2017a), to constraints on the equation of state of very dense matter (Abbott et al. 2018). In addition, and especially relevant to this work, it has confirmed the prediction (Eichler et al. 1989;Narayan et al. 1992) that at least some short gamma-ray bursts (SGRBs) are produced in binary neutron star (BNS) mergers (e.g. Abbott et al. 2017a,c;Goldstein et al. 2017;Hallinan et al. 2017;Margutti et al. 2017;Troja et al. 2017;Lazzati et al. 2017aLazzati et al. , 2018Mooley et al. 2018;Ghirlanda et al. 2019). More specifically, the presence of a relativistic jet, which characterizes SGRBs, was uncovered also in GRB170817A, despite its lower luminosity compared to that of the standard SGRBs.
The side-emission of the jet, which is the one more likely to be seen in GW-detected BNS mergers, plays an important role for the interpretation of the observations. The jet shape, and especially its wings, are molded by the BNS ejecta in which it propagates. A top-hat jet will become a structured outflow by the time it emerges, as a result of the interaction between the relativistic jet and the ejecta material. 1 The importance of the dynamical ejecta in molding the properties of the outflow has already been recognized in a number of works, with special emphasis in the modeling of the multi-wavelength spectrum of GRB170817A (e.g. Murguia-Berthier et al. 2014;Nagakura et al. 2014;Lazzati et al. 2017bLazzati et al. ,a, 2018Xie et al. 2018;Bromberg et al. 2018;van Eerten et al. 2018;Lamb & Kobayashi 2018;Wu & Mac-Fadyen 2018;Gottlieb et al. 2018;Beniamini & Nakar 2019;Gottlieb et al. 2019;Kathirgamaraju et al. 2019;Geng et al. 2019). Independently of the connection to the SGRBs, but also very important, the ejecta mass and its profile contain information on the equation of state of NSs, as well as being the primary site for the production of the heaviest elements in the Universe via the r-process (Lattimer & Schramm 1974). A recent study by Lazzati & Perna (2019) has shown how, for a given engine, the distribution of the Lorentz factor and the energy per unit solid angle with viewing angle has a noticeable dependence on the ejecta mass, its speed, and on the time delay between the merger and the jet onset.
Besides one very recent exception (Pavan et al. 2021), SGRB jet propagation studies have so far used analytical, axially symmetric (and often even isotropic) density profiles for the ejecta, constructed to broadly reproduce what found in general relativistic simulations of BNS mergers (e.g. Lazzati et al. 2018;Xie et al. 2018;Geng et al. 2019;Kathirgamaraju et al. 2019;Nathanail et al. 2020;Salafia et al. 2020;Lazzati et al. 2020;Hamidani & Ioka 2021;Gottlieb et al. 2021;Murguia-Berthier et al. 2021;Urrutia et al. 2021). However, those ejecta profiles miss relevant small scale structure and, more generally, lack a direct connection with the realistic three-dimensional environment of the BNS merger. In this work, we present a three-dimensional relativistic hydrodynamic simulation in which a relativistic jet propagates through the ejecta produced by a BNS merger, as pre-computed via a fully general relativistic magnetohydrodynamic (GRMHD) simulation. In particular, we focus on a physical model where the mass of the environment around the central engine at the jet launching time is as high as ∼ 0.1 M , without a pre-evacuated funnel along the rotation axis. This corresponds to a very different physical regime compared to Pavan et al. (2021), where the initial environment is much lighter and contains a low density funnel along the rotation axis. The breakout time found in Pavan et al. (2021) is indeed much shorter (∼ 30 ms in their fiducial case), and their focus is not on the details of the pre-breakout jet's head propagation, but rather on the post-breakout evolution.. We explore the dynamics and energetics of the propagation up to breakout, occurring after ∼ 0.6 s. In addition, we compare the first 0.2 s of jet evolution with an equivalent case where the environment is isotropic and its density distribution is set via an analytical prescription (with the same total mass and kinetic energy). Based on the emerging differences, we find that for precision better than ∼ 20 per cent in the jet's head velocity inside the ejecta, realistic initial data from actual BNS merger simulations need to be used (see also Pavan et al. 2021). The paper is organized as follows: Sec. 2 describes our numerical methods and initial conditions. The simulations results are presented in Sec. 3, and we summarize and conclude in Sec. 4.

METHODS
We carry out two 3D special relativistic hydrodynamic (SRHD) simulations to study the effect of realistic BNS ejecta environments on the propagation of a SGRB relativistic jet. Our fiducial simulation is based on the direct import of density, pressure, and 3-velocity distributions from a BNS merger simulation from Ciolfi et al. (2019). Our control simulation is identical in all aspects to the fiducial simulation with the exception of the BNS ejecta, which are prescribed with an analytic profile that matches the total mass and average velocity of the fiducial simulations, but removing the polar dependence as well as all the small scale density, pressure, and/or velocity inhomogeneities. Both simulations are carried out with the public code FLASH (Fryxell et al. 2000) in 3D Cartesian coordinates. The inner boundary is set at a distance y 0 = 5 × 10 7 cm = 500 km from the central engine and the domain extends to ±6.25 × 10 8 cm in both the x and z directions (equatorial plane) and to 1.25 × 10 9 cm in the y direction, along which the jet flows. A nested grid of refinement levels is adopted (analogously to Morsony et al. 2007), with a maximum resolution at the inner boundary where the jet is injected of ∆x = ∆y = ∆z = 1.22 × 10 6 cm= 12.2 km.
A relativistic jet is injected as an inflow boundary condition. The jet has an injection half opening angle θ j = 10 • , a constant luminosity L j = 10 50 erg s −1 for the one-sided jet, an injection Lorentz factor Γ 0 = 5, and an internal energy that allows for an asymptotic Lorentz factor η = 300. The jet base is resolved with approximately 164 elements, which was shown by López-Cámara et al. (2013) to provide sufficient resolution for convergence (see their Figures 10 and  14). These injection parameters are consistent with the ones Figure 1. Pseudocolor volume rendering of the 3D rest-mass comoving density distribution of the merger ejecta at six different times in the simulation (green-yellow shading). Semi-transparent isocontours of the jet Lorentz factor are overlaid in blue (cocoon material, β = 0.1), red (jet material, Γ = 2), and black (jet core, Γ = 20). A bar of size 1000 km is shown in each panel for scale. The simulation -which covers only one hemisphere -was reflected across the equator to produce this figure. Apparent north-south asymmetries are only due to the viewing perspective and not intrinsic to the simulation. derived by Lazzati et al. (2020) for GW170817 2 (see their  Table 3). Since the BNS ejecta have a net outward velocity (for both the fiducial and control cases), the entire equatorial plane is set as an inflow boundary. Outside of the jet base a non-relativistic (NR) inflow is implemented according to an analytical fit based on the GRMHD simulation from which the initial condition was taken. In practice, the density, pressure, and velocity field in the BNS merger simulation domain that did not overlap with the FLASH domain were used to create an analytical model to predict the density, pressure, and velocity of the inflow in the FLASH domain at later times. This NR inflow has density and velocity that decline with time and the velocity vectors are oriented diverging from a single point where the central engine is located. Even though this additional inflow is rather simplified, the head of the jet never comes in contact with it and we con-sider its details irrelevant for the jet propagation dynamics. All other boundaries are set as outflowing.
The specific BNS merger simulation under consideration was carried out with the WhiskyMHD code (Giacomazzo & Rezzolla 2007;Giacomazzo et al. 2011) and is the one presented in Ciolfi et al. (2019), reproducing an equal mass 1.35-1.35 M system 3 and employing as equation of state (EOS) a piece-wise polytropic approximation of APR4 (Akmal et al. 1998;Endrizzi et al. 2016). We refer to Ciolfi et al. (2019) for details on the numerical and physical setup. This is one of the GRMHD simulations with the longest post-merger massive NS remnant evolution to date ( 100 ms; see also Ciolfi 2020). The initial condition for the FLASH simulation is taken 71 ms after merger, corresponding to the jet launching time plus the propagation time to the inner boundary of the simulation, which would be ∼ 1.6 ms at the speed of light and ∼ 11 ms at the average propagation speed within 3 Corresponding to a total mass only ∼ 1 per cent lower than that of GW170817 (Abbott et al. 2019). . Pseudocolor image of the region of the jet around the head (very low density, red colors) at a time close to the breakout from the ejecta (high density, blue color). False colors show different levels of the logarithm of comoving density. The left panel is displayed with color cuts chosen to emphasize the structure of the ejecta, which are crossed by multiple shock systems caused by the inhomogeneities within the ejecta themselves. The right panel, instead, shows the same area but with color cuts that emphasize the internal structure of the jet and the fact that the low-density, fast material is propagating along a slightly off-axis channel. the ejecta measured in the simulation. In Ciolfi et al. (2019), it was shown that inducing the collapse of the NS remnant into a BH after a survival time of about 70 ms leads to a BH mass and dimensionless spin of 2.5 M and 0.5, respectively, surrounded by an accretion disk of mass of the order of ∼ 0.1 M , i.e. suitable to power a SGRB jet. We note that the two NSs composing the merging BNS were initially endowed with a strong magnetic field (maximum field strength of 10 16 G and magnetic energy of ∼ 10 48 erg), which eventually contributed to enhancing the post-merger matter outflow and thus the total mass of the baryon-loaded cloud surrounding the central object. At 71 ms after merger, this mass is ∼ 0.1 M . The original GRMHD dynamics was embedded in a highdensity artificial atmosphere (6.18 × 10 6 g/cm 3 ), unsuitable for the propagation of a relativistic jet on the larger domain of our FLASH simulation. To allow for a lower density at larger distances without altering the GRMHD dynamics we substituted the constant density floor with a static atmosphere with the same central density but with an exponential cutoff: ρ(r) = 6.18 × 10 6 exp − r 3 × 10 8 cm g cm 3 . (1) The total mass contained within our FLASH domain is 0.6 M , of which 0.56 M are due to the artificial atmosphere and 0.042 M are due to the ejecta. In the ejecta region (initially located within a distance r < 2100 km), however, the ejecta dominates, since the artificial atmosphere contributes only 0.027 M . The jet propagation out to the breakout is therefore driven predominantly by the ejecta and not by the artificial atmosphere in which they are embedded.
Regarding the angular distribution of the material, a distinctive feature is that in this case no lower density funnel is present along the spin axis of the remnant, due to the nearly isotropic magnetically driven matter outflow (Ciolfi et al. 2019). This implies a more difficult initial propagation of the jet, with a slower and irregular jet's head advancement (see Section 3). We also note that we are using a matter distribution above 500 km imported from a non-collapsing simulation, which corresponds to the implicit assumption that, if the jet is launched by a central BH, the time between the collapse and the jet emergence at 71 ms is short enough that the environment is not yet affected by the newly-formed BH at distances equal or larger than 500 km.
3. RESULTS Figure 1 shows pseudocolor volume renderings of the logarithmic density at various times in the simulation. Semitransparent isocontours of the jet's Lorentz factor are also shown. Our simulation shows that the jet does reach the leading edge of the ejecta, and breaks out at about 0.6 seconds after injection. Since the typical duty cycle of a SGRB engine is of the order of one second or more, we conclude that SGRB jets can reach breakout while the central engine is still active, even in presence of rather massive ejecta like the ones that we have considered here.
By looking at the results in more detail, it is interesting to follow the jet propagation as it encounters -and navi- gates around -higher density regions in the ejecta. Figure 2 shows the X-Y plane projection of the jet head region in log density pseudocolors. In the figure, the blue colors show the BNS ejecta, while the red colors show the jet. The left-side panel is optimized for showing differences of density in the ejecta, while the right-side panel is optimized for showing the jet internal structure. It is clear by visual inspection of the figure that multiple shock systems are present in the ejecta, driven by the jet propagation, which proceeds discontinuously. The distance that the various components of the system (ejecta, jet, reverse shock in the jet) travel since the injection time is shown in Figure 3, together with the centroid of the jet computed with respect to the jet axis at the initial time.
Here and in the following, we call reverse shock the location at which the jet material is first decelerated after injection. It is the shock along the jet spine at the smallest distance from the injection boundary. The centroid is seen to oscillate by ∼ 1 • − 2 • as it propagates, as a result of inhomogeneities in the ejecta. A striking feature of Figure 3 is the almost one-to-one correspondence between setbacks in the jet head propagation and deviation of the jet centroid form its propagation direction. This correspondence (highlighted in the figure with black vertical dashed lines) suggests a scenario in which the jet propagates along a direction until it impacts onto an overdensity. When that happens, the jet head shifts towards a region of lower density. On the short time-scale the shift halts the jet propagation. However, it quickly recovers and eventually propagates faster than an identical jet that propagates in ejecta with the same mass and velocity but no inhomogeneities.
This effect is emphasized in Figure 4 that shows a comparison between our fiducial and control simulations. Unfortunately, due to limited resources, it was impossible to compute the evolution of the control simulation all the way to breakout. A comparison for the first 200 ms is, however, very informative nonetheless. The solid lines in Figure 4 show the same distances shown in Figure 3 but for the control simulation, compared to the fiducial simulation in thin dashed lines. The control simulation appears to propagate in a continuous and monotonic way, confirming that the head oscillations and propagation setbacks seen in the fiducial simulation are caused by the small scale overdensities. Overall, despite the setbacks encountered, the jet in the fiducial simulation propagates ∼20 per cent faster than the control simulation jet, at least in the initial 200 ms of the simulation for which a direct comparison can be performed. The jet centroid is also seen to oscillate in the control case (red line in the figure, see also Zhang et al. 2004;López-Cámara et al. 2013). However, the oscillation has smaller amplitude and has a much longer timescale, compared to the sudden changes seen in Figure 3.
We finally investigated the energy balance in the various components of the simulated system. For this aim, we divide the simulation domain into jet material, cocoon material, and ejecta material. This selection is based on the local asymptotic Lorentz factor: We define as jet all the material with Γ ∞ ≥ 20, as cocoon all the material with 1.007 ≤ Γ ∞ < 10 and as ejecta all the material with Γ ∞ < 1.007 (note that for Γ ∞ = 1.007 the maximum attained speed is v = 3.5 × 10 9 cm s −1 ∼0.1c). The time evolution of the energy content of the various components in our fiducial simulation is displayed in Figure 5, together with the energy injected by the central engine as the jet inflow. The jet energy is initially very low, showing that for the initial 200 ms after the injection the jet is struggling to form and all the energy is used to evacuate the funnel and inflate the cocoon. Notice also how a significant fraction of the engine energy is used to evacuate the channel, and so the cocoon energy is only about half of the total injected one in the first ∼ 200 ms. After that initial time, a large scale jet forms and its energy grows quickly, eventually settling to about 30 per cent of the  engine energy. The cocoon energy also grows slightly, settling on ∼ 50 per cent of the injected energy. The missing ∼ 20 per cent was used as mechanical energy for evacuating the path for the jet propagation. Finally, we note that the growth of the energy in the NR ejecta is due to the inflow of NR material at the equatorial boundary throughout the simulation (see Sect 2).

DISCUSSION AND CONCLUSIONS
The evolution of a relativistic jet driven by a BNS merger is molded into a structured outflow by the interaction with its immediate environment, i.e. the ejecta produced by the merger itself. The emission generated by the outflow, when compared with theoretical models such as in the case of GRB170817A, allows us to infer physical properties not only of the large-scale outflow (e.g., Lazzati et al. 2018) but also of the jet and its engine (Lazzati et al. 2020). Therefore, it is of paramount importance to refine as much as possible our understanding of the jet-ejecta interaction and of the ensuing outflow. One way of gaining better insight into these processes is by using BNS ejecta from merger simulations as an initial condition in jet simulations, rather than simple analytical parameterizations of the ejecta environment (see also Pavan et al. 2021).
In this paper we have presented a 3D special relativistic hydrodynamic simulation of a jet propagating in the ejecta as obtained from a GRMHD simulation of an equal-mass BNS merger. As a control case, we also present a simulation of an identical jet propagating in an analytical rendering of the same ejecta, removing small scale inhomogeneities and any polar dependency of the ejecta density and velocity field. Our simulations are similar in scope to those recently presented in Pavan et al. (2021). However, they differ in the code used (Pluto vs. FLASH), the BNS simulation physics (low vs. high magnetization), the mass of the ejecta (low vs. high mass) and the size of the domain, allowing for an interesting comparison that furthers the scope of our investigation.
Our fiducial simulation (with realistic BNS merger ejecta), shows that the jet propagation in a inhomogeneous medium is non-uniform. The head of the jet oscillates around the system symmetry axis (red line in Figure 3) and its outward propagation is characterized by numerous sudden slow-down instances, which occur almost any time the centroid shifts its position angle. While some of these features have been seen also in 3D simulations of jets propagating in smooth material (e.g., Zhang et al. 2004;López-Cámara et al. 2013), the behaviour is enhanced by the presence of inhomogeneities, as seen via the comparison between our fiducial and control simulations, both run in 3D (Figure 4). Despite the somewhat erratic behavior of the jet in the fiducial simulation, its advancement is faster than in the control case. We attribute this to the fact that the centroid oscillations allow the jet to find a path of least resistance through the ejecta. This optimizes the overall jet propagation, even if the immediate, short time-scale effect is to slow down the jet head.
The overall comparison between the fiducial and control simulations shows that differences exist, even though for the case at hand they are not large. The jet propagation in the analytic environment captures quite well the overall properties of the phenomenon, and the jet's head velocity is correct to within 20 per cent from the more time-consuming case of realistic ejecta. Our fiducial simulation also shows that the jet and its surrounding cocoon share the injected energy almost equally. More interesting, a significant fraction of the injected energy is lost to mechanical work on the ejecta by cleaning a path for the jet to propagate and the cocoon to expand. This loss of energy is an important effect that should be considered in analytic models such as those of Lazzati & Perna (2019); Salafia et al. (2020).
Among the important results of our simulation is the time that the jet takes to break out, ∼ 0.6 sec in our 0.1 M ejecta. This, together with the time delay between merger and jet launch, and the (angular-dependent) time delay for the radiation to emerge from the jet, constitutes an important component of the measurable time delay between the gravitational wave signal and the electromagnetic one, and hence is a crucial quantity to compute. Comparing our results to those of Pavan et al. (2021) we find general agreement. Not surprisingly, we also find that their jets break off the ejecta much faster, given the lower mass of the ejecta considered in their case and the presence of pre-existing even lower density funnel along the jet propagation path. These differences are mainly brought about by the different magnetization of the progenitor NSs in the merger simulation, and jet dynamics could therefore be used to constrain the stars' original magnetic fields.
Looking into the future, while employing ejecta distributions from BNS merger simulations already constitutes a considerable improvement towards modeling a realistic setting for the production of SGRBs, further improvements will require moving from a pure hydrodynamical to a magnetohydrodynamical jet, as well as allowing the incipient luminosity of the jet to have an angular dependence (rather than being a simple top-hat as assumed here) and a temporal dependence derived from the merger simulations.
More generally, simulating a variety of jet models propagating within a variety of ejecta structures will allow for a more direct mapping of input physics to observations, hence deepening our understanding of the astrophysics of BNS mergers and the powerful jets that they produce.
We thank the anonymous referee for useful suggestions that led to a much improved version of the original manuscript. DL acknowledges support from NASA grants 80NSSC18K1729 (Fermi) and NNX17AK42G (ATP), Chandra grant TM9-20002X, and NSF grant AST-1907955. RP acknowledges support by NSF award AST-2006839 and from NASA (Fermi) award 80NSSC20K1570. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.