Two Steps Forward and One Step Sideways: The Propagation of Relativistic Jets in Realistic Binary Neutron Star Merger Ejecta

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 (GWs) 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, 2017c; Goldstein et al. 2017; Hallinan et al. 2017; Margutti et al. 2017; Troja et al. 2017; Lazzati et al. 2017a, 2018; Mooley 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.

Detailed modeling of this source (Granot et al. 2017; Margutti et al. 2017; Murguia-Berthier et al. 2017; Alexander et al. 2018; Bromberg et al. 2018; Gottlieb et al. 2018; Hotokezaka et al. 2018; Ioka & Nakamura 2018; Lamb & Kobayashi 2018; Lazzati et al. 2018; Margutti et al. 2018; Mooley et al. 2018; Ruan et al. 2018; Troja et al. 2018; Xie et al. 2018; Beniamini et al. 2019; Beniamini & Nakar 2019; Geng et al. 2019; Kathirgamaraju et al. 2019) allowed for the discovery that the lower luminosity was due to a viewing angle effect, i.e., the observer being positioned at an angle ∼14°–28° with respect to the jet axis (Lazzati et al. 2018; Mooley et al. 2018; Ghirlanda et al. 2019). Therefore, this event alone has demonstrated that at least some of the diversity in SGRBs is due to different viewing angles, and that their gamma-ray emission drops off slowly at large angles.

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. ¹¹ 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. 2017b, 2017a; Bromberg et al. 2018; Gottlieb et al. 2018; Lamb & Kobayashi 2018; Lazzati et al. 2018; Troja et al. 2019; Wu & MacFadyen 2018; Xie et al. 2018; Beniamini et al. 2019; Beniamini & Nakar 2019; Geng et al. 2019; Gottlieb et al. 2019; Kathirgamaraju 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.

With 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; Lazzati et al. 2020; Nathanail et al. 2020; Salafia et al. 2020; Gottlieb et al. 2021; Hamidani & Ioka 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 (3D) environment of the BNS merger. In this work, we present a 3D 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% 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).

This Letter is organized as follows. Section 2 describes our numerical methods and initial conditions. The simulations results are presented in Section 3, and we summarize and conclude in Section 4.

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₀ = 5 × 10⁷ cm = 500 km from the central engine and the domain extends to ±6.25 × 10⁸ cm in both the x and z directions (equatorial plane) and to 1.25 × 10⁹ 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⁶ 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⁵⁰ erg s⁻¹ for the one-sided jet, an injection Lorentz factor Γ₀ = 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 derived by Lazzati et al. (2020) for GW170817 ¹² (see their Table 3). Because 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 consider 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 ¹³ and employing as equation of state 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 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 M_⊙, 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¹⁶ G and magnetic energy of ∼10⁴⁸ 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 high-density artificial atmosphere (6.18 × 10⁶ gcm⁻¹), 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:

$$\begin{eqnarray}&&\rho (r)=6.18\times {10}^{6}\,\exp \left(-\displaystyle \frac{r}{3\times {10}^{8}\,\mathrm{cm}}\right)\qquad \displaystyle \frac{{\rm{g}}}{{\mathrm{cm}}^{3}}.\end{eqnarray} \tag{ 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, as 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.

Figure 1 shows pseudocolor volume renderings of the logarithmic density at various times in the simulation. Semi-transparent 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 s after injection. As 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.

Zoom In Zoom Out Reset image size

By looking at the results in more detail, it is interesting to follow the jet propagation as it encounters and navigates 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 panel is optimized for showing differences of density in the ejecta, while the right 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.

Zoom In Zoom Out Reset image size

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 toward a region of lower density. On the short timescale 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.

Zoom In Zoom Out Reset image size

This effect is emphasized in Figure 4, which 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, however, is 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% 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 a much longer timescale compared to the sudden changes seen in Figure 3.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\infty }={\rm{\Gamma }}\left(1+4\displaystyle \frac{p}{\rho {c}^{2}}\right).\end{eqnarray} \tag{ 2 }$$

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⁹ cm s⁻¹ ∼ 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% of the engine energy. The cocoon energy also grows slightly, settling on ∼50% of the injected energy. The missing ∼20% 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 Section 2).

Zoom In Zoom Out Reset image size

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 Letter 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 versus FLASH), the BNS simulation physics (low versus high magnetization), the mass of the ejecta (low versus 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 behavior 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 timescale 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% 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) and Salafia et al. (2020).

Among the important results of our simulation is the time that the jet takes to break out, ∼0.6 s 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 GW 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 therefore jet dynamics could 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 toward modeling a realistic setting for the production of SGRBs, further improvements will require moving from a pure hydrodynamical to a magneto-hydrodynamical 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. D.L. acknowledges support from NASA grants 80NSSC18K1729 (Fermi) and NNX17AK42G (ATP), Chandra grant TM9-20002X, and NSF grant AST-1907955. R.P. 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.

Software:FLASH (Fryxell et al. 2000), VisIt (https://wci.llnl.gov/simulation/computer-codes/visit), Python (www.python.org).