Jet-cocoon outflows from neutron star mergers: structure, light curves, and fundamental physics

The discovery of GW170817, the merger of a binary neutron star (NS) triggered by a gravitational wave detection by LIGO and Virgo, has opened a new window of exploration in the physics of NSs and their cosmological role. Among the important quantities to measure are the mass and velocity of the ejecta produced by the tidally disrupted NSs and the delay - if any - between the merger and the launching of a relativistic jet. These encode information on the equation of state of the NS, the nature of the merger remnant, and the jet launching mechanism, as well as yielding an estimate of the mass available for r-process nucleosynthesis. Here we derive analytic estimates for the structure of jets expanding in environments with different density, velocity, and radial extent. We compute the jet-cocoon structure and the properties of the broadband afterglow emission as a function of the ejecta mass, velocity, and time delay between merger and launch of the jet. We show that modeling of the afterglow light curve can constrain the ejecta properties and, in turn, the physics of neutron density matter. Our results increase the interpretative power of electromagnetic observations by allowing for a direct connection with the merger physics.


INTRODUCTION
The discovery of GRB170817A in association with the gravitational wave (GW) source GW170817 produced by a double neutron star (NS) merger (Abbott et al. 2017), has opened a new era in the study of short gamma-ray bursts (SGRBs) as well as of NSs, both separately, but especially when the information from GWs and the electromagnetic counterparts are combined. Prior to the association with this GW event, only the brightest SGRBs were followed at longer wavelengths. However, even under these conditions, the afterglows were generally dim, and rarely full broadband coverage was possible (Fong et al. 2015).
GRB170817, despite having an isotropic luminosity ∼ 10 4 times lower than that of 'standard' SGRBs, was very bright due to its relatively small distance of 40 Mpc, having the initial trigger been in GWs rather than in γ-rays as for the standard cosmological SGRBs (Goldstein et al. 2017;Savchenko et al. 2017). A massive broadband follow-up campaign to GW170817 allowed an unprecedented data coverage, revealing a kilonova/macronova (e.g., Arcavi et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Smartt et al. 2017;Soares-Santos et al. 2017), as well as the signatures of a relativistic jet (Alexander et al. 2017Margutti et al. 2017Margutti et al. , 2018Lazzati et al. 2018;Mooley et al. 2018;Ghirlanda et al. 2019). The mass of the ejecta was inferred to be ≈ 0.065 M by Kasen et al. (2017) and between 0.042 and 0.077 M by Perego et al. (2017) from modeling of the the kilonova light curve. The former analysis assumed two components, while the latter assumed three components. As pointed out by Barnes et al. (2016), key uncertainties in kilonova modeling include the emission profiles of the radioactive decay products -non-thermal β particles, α particles, fission fragments, and γ rays -and the efficiency with which their kinetic energy is absorbed by the ejecta.
The ejecta mass, on the hand, is an important quantity to measure in BNS mergers: it contains information on the equation of state (EoS) of the NS (e.g. Hotokezaka et al. (2013)), and is a favorable site for the production of the heaviest elements in the Universe via the r-process (Lattimer & Schramm 1974). As LIGO and Virgo have begun a new observing run at higher sensitivity, more NS-NS events will be expected, and hence any additional diagnostics of the ejecta mass will be extremely valuable.
The role of the ejecta in determining the properties of the jet has been studied under different angles by a number of investigators, especially within the context of explaining the multi-wavelength observations of GRB170817A (e.g. Lazzati et al. 2017bLazzati et al. , 2018Xie et al. 2018;Bromberg et al. 2018;van Eerten et al. 2018;Lamb & Kobayashi 2018;Wu & MacFadyen 2018;Gottlieb et al. 2018Gottlieb et al. , 2019Kathirgamaraju et al. 2019;Geng et al. 2019). These data were found to be incompatible with a top-hat jet, but rather required a 'structured jet', as produced by the interaction of a relativistic outflow with the high-velocity ejecta generated during the merger process.
The role of the ejecta was also recently discussed (Barbieri et al. 2019) in the context of NS-BH (black hole) mergers, together with the role of the mass and spin of the BH, in determining the observable electromagnetic signatures accompanying the coalescence. Additionally, the interaction of the ejecta themselves with the surrounding interstellar medium is expected to produce a long-lasting radio flare (Hotokezaka & Piran 2015).
Given the importance of the dynamical ejecta in molding the properties of the relativistic outflow that propagates in it, in this paper we perform a detailed study aimed at connecting the properties of this matter (and in particular its mass and velocity), with the jet-cocoon structure that develops in it, and the resulting afterglow emission in representative bands. The sensitivity of this radiation to the ejecta properties (as will be shown here) will provide an additional element linking observables to fundamental physics.
Our paper is organized as follows: We begin (Sec. 2) by deriving analytic expressions for the structured-jet properties, connecting the jet morphology to the properties of the merger ejecta in which it propagates. Aided by the results of hydrodynamic numerical simulations, we derive energy and Lorentz factor profiles as a function of the off-axis angle. We feed these profiles to an external shock synchrotron code and compute the corresponding afterglow light curves. Finally, we discuss the connection between this observables and the NS EoS (Sec. 4). We summarize and conclude in Sec. 5.

THE STRUCTURED JET MODEL
In this section we discuss an analytic framework to compute the properties of the structured jet from a binary NS merger, whithin the framework of the Kompaneets approximation (e.g., Begelman & Cioffi 1989;Matzner 2003;Lazzati & Begelman 2005;Bromberg et al. 2011Bromberg et al. , 2014. We first discuss the input of hydrodynamic simulations in setting up an analytic profile for the energy and Lorentz factor as a function of the off-axis angle. Subsequently we derive the scale angles and energy of the jet and cocoon components. (a) The energy and Lorentz factor profiles To describe the angular distribution of the relativistic outflow at large radii (well after the jet interaction with the NR ejecta) we adopt an exponential model, as supported by the results of numerical simulations ( Fig. 1): The constants A and B are determined by the conditions that A Ω e − θ θ j dΩ is the total energy in the relativistic jet and B Ω e − θ θc dΩ is the total energy in the cocoon. The variables θ j and θ c indicate the scale angles for the jet and cocoon, respectively. Note that θ j is in general smaller than the injected jet opening angle θ j,inj due to hydrodynamic recollimation.
A similar profile is adopted for the Lorentz factor, but the value is not allowed to be less than unity: where Γ core is a free parameter of the model. The cocoon maximum Lorentz factor Γ c is instead the result of the complex interaction between the jet and the ejecta, and is mainly controlled by the mixing between the two components. This cannot be modeled analytically and we assume here that Γ c = Γ core /10, a value that fits well our numerical results (see the bottom panel of Figure 1). Given the fact that we only have one simulation to compare to, this value should be considered tentative. It should be remembered, however, that the initial Lorentz factor does not affect significantly the afterglow light curve at times longer than the deceleration time.
We now focus on calculating the energy profile (top panel of Figure 1). To evaluate the values of the four structure constants A, B, θ j , and θ c we model the propagation of a jet through the non-relativistic ejecta. We adopt the Kompaneets approximation (Begelman & Cioffi 1989;Matzner 2003;Lazzati & Begelman 2005;Bromberg et al. 2011Bromberg et al. , 2014 and identify the scale angles θ j and θ c as the opening angles of the jet and cocoon at break-out.

(b) The NS ejecta
The ejecta structure described above is created by the interaction of a relativistic jet with non-relativistic ejecta from the merger. We model the ejecta as a non-relativistic wind with constant velocity v w and mass loss rateṁ w that is released during the merger and by the remnant a time δt before the jet is launched. The ejecta are assumed to be spherically symmetric within a solid angle Ω w . The density of the ejecta is therefore given by: (c) The jet At time δt after the NR ejecta are launched, a jet is launched from the central engine. We consider a top-hat jet of luminosity L j with uniform injection Lorentz factor Γ j,inj within an injection opening angle θ j,inj . The jet is injected very simple, its "structure" is acquired during the interaction with the ejecta.
One can qualitatively understand the jet propagation as follows. The jet head propagates through the ejecta driving a bow shock that inflates a cocoon of hot material that enshrouds the jet. The cocoon pressure acts on the jet, collimating it into a narrower angle. The cocoon is trapped as well inside the ejecta and drives a shock into them. The model has three fundamental unknowns: the velocity of the head of the jet v jh , the opening angle of the jet at breakout θ j , and the opening angle of the cocoon at breakout θ c . The system can be solved because it is constrained by three pressure balances as will be discussed in the part (d) below.
The jet propagates inside the ejecta and eventually breaks out when the outer radius of the ejecta coincides with the position of the jet head, i.e. when v w (t bo + δt) = v jh t bo . (4) This yields a break out time (measured from the launching of the jet) and a breakout radius (Murguia-Berthier et al. 2014) (d) The pressure balances The structures of the jet and cocoon are defined by three pressure balances: at the head of the jet, the ram pressure of the jet on the contact discontinuity is balanced by the ram pressure of the merger ejecta on the same contact discontinuity. At the jet-cocoon transition, the cocoon thermal pressure is balanced by the jet pressure (be either the jet internal pressure or the ram pressure due to the deflection of the jet material). Finally, at the cocoon-ejecta transition, the thermal pressure of the cocoon is balanced by the ram pressure of the merger ejecta. These conditions yield the equations: This set of equation allows us to solve for the three unknowns, v jh , θ j , and θ c . The first pressure balance (Eq. 7) has been studied in detail previously (e.g. Marti et al. 1994;Matzner 2003;Lazzati & Begelman 2005;Bromberg et al. 2011) and yields the speed with which the jet head advances in the merger ejecta. In static ejecta, the velocity of the jet head reads where we have made the assumption that the jet material moves with v j ≈ c and that the enthalpy density is dominated by the pressure component. If the jet propagates in a moving medium, the velocity in Eq. (10) Eq. (10) can be written in terms of unknown quantities and model parameters. The wind density is given by Eq. (3) and, using 4p j Γ 2 j r 2 = L j /(Ω j c), we obtain: For the second pressure balance, we assume that the cocoon energy is dominated by radiation so that P c = E c /(3V c ), where V c = Ω c r 3 is the cocoon volume. Combining the above, we obtain an expression for the cocoon pressure, The ram pressure of the wind material on the cocoon is given by p ram,wc = ρ w v 2 ⊥c , where v 2 ⊥c = Ω c r 2 /(πδt 2 ) is the velocity squared of the shock front driven by the cocoon into the ejecta. We obtain: Finally, let us consider the jet pressure, which has two components. First, the internal pressure of the jet, which is a relativistic invariant, and given by:  Figure 3. The jet parameters are: Lj = 10 50 erg/s, θj,inj = 10 • , r0 = 10 7 cm, Γ0 = 1. The wind is assumed to be isotropic (Ωw = 4π).
Additionally, one must consider the ram pressure of the jet material that is deflected from its initially radial velocity into a cylindrical flow in which the velocity is predominantly axial. This pressure component is obtained by correcting the jet ram pressure with a geometrical factor sin 2 θ j,inj , yielding p ram,jc = L j Ω j cr 2 sin 2 θ j,inj . Off-axis angle (degrees) In most cases, unless the jet is not significantly recollimated, the ram pressure dominates and hence we will use Eq. (16) for the jet pressure on the cocoon hereafter.
Using Eq. 13 and 14 in Eq. 8 we obtain an equation for the cocoon solid angle as Alternatively, an expression for the cocoon solid angle can be derived using Eq. 13 and 16 in Eq. 9 to yield Equating Eq. (17) to the square of (18) and solving for the jet solid angle we obtain: The right hand side of this equation can now be inserted for the jet opening angle in Eq. (11) and solved for the jet head velocity, by using Eq. (6) for the radius r.
Finally, the cocoon energy can be derived from and, thanks to energy conservation, the remaining jet energy E j = E j,inj − E c , where E j,inj is the energy injected at the base of the jet. To validate this model we show in Figure 1 a comparison between the numerical results of Lazzati et al. (2017) and the energy and gamma profiles derived from the above equations. While there are some differences due to the complexity of the interaction between a relativistic jet and the relatively slow and dense ejecta, the main features are well represented in the analytic model. We stress that the simulation was used as an input for the functional form in Equations 1 and 2, but the scale angles and cocoon energy used to draw the green curve in the figure are obtained form the analytic model. The general agreement between the numerical result (in blue) and the model (in green) is therefore not a circular argument. Figure 2 shows the jet and cocoon opening angles (top and middle panels), and cocoon energies (bottom panels) as a function of the total mass in the ejecta at the time of the jet injection and of the delay between the launching of the Time since merger (days)   Figure 3 shows the energy (top panels) and Lorentz factor (bottom panels) profiles for a few selected cases. The cases shown in Figure 3 are indicated in Figure 2 with colored star symbols.

AFTERGLOW LIGHT CURVES
We compute the afterglow light curves by feeding the energy and Lorentz factor profiles into a semi-analytic external shock code. We adopt the same version of the radiative code that was previously used by Lazzati et al. (2017bLazzati et al. ( ,a, 2018 and Perna et al. (2019). Afterglow light curves are calculated assuming a Blandford-McKee self-similar profile (Blandford & McKee 1976) for the blast wave structure behind the shock. The radiation reaches the observer at a time which depends on both the emission radius R em , as well as on the position of the emitting element with respect to the line of sight. For an emitting element at radius r em located at an angle θ los from the line of sight, the arrival time is given by t = t lab − r em /c cos θ los where t lab = c −1 β −1 dr (Panaitescu & Mészáros 1998;Granot et al. 1999).
The external material is assumed to be a low-density, uniform interstellar medium of density n ISM , and the shock microphysics is parameterized, as customary, with the equipartition parameters e and B and with the non-thermal electron's energy distribution index p el . All electrons are assumed to be accelerated into the non-thermal distribution. We perform our light curves calculations adopting typical average values for the shock microphysics parameters, as inferred in the afterglow modeling of the Swift-detected short GRBs (Fong et al. 2015) and of GRB170817 (Lazzati et al. 2017a). These are: e = 0.1, B = 0.01, and p el = 2.3. The interstellar medium is assumed to be uniform and tenuous, as expected in the surrounding of a BNS merger, with density n ISM = 10 −2 cm −3 . The observer angle with respect to the jet axis is set to θ jet = 30 • , the typical viewing angle expected for Ligo/Virgo detected events (Schutz 2011) and the distance to the burst is set to 40 Mpc as a representative case to relate to GRB170817.
Examples of light curves in X-ray (1 keV) and radio (6 GHz) are shown in the top and bottom panels of Figure 4, respectively. Interestingly, the light curves are noticeably different from one , giving hope that the ejecta mass and velocity, as well as the time delay δt, can be constrained with a high quality dataset. In particular, the slope of the slowrising (or, in some cases, slowly decaying) phase before the maximum and the prominence of the feature at the peak are characteristic of the jet-cocoon configuration. A larger ejecta mass results in a more energetic cocoon confined in a narrower scale angle θ c (Figure 2) that moves with lower Lorentz factor (Figure 3). The lower Lorentz factors cause a late onset of the afterglow at both X-ray and radio frequencies, which is therefore a signature of significant mass ejection. Large ejecta masses also cause the jet scale angle θ j to narrow significantly, resulting in a clearer jet feature at the light curve peak, especially at radio frequencies (see, e.g., the green light curve in the top-center panel of Figure 4). Very low mass ejecta produce instead a very weak but fast moving cocoon and wide angle jet. These light curves peak early, especially in X-rays, and more closely resemble those of canonical, onaxis jets. Note that in the complete absence of ejecta, the light curves peak late, since no cocoon is present, as shown with the case of a pure jet by Perna et al. (2019). This might seem contradictory, but there is indeed a big difference between the complete absence of ejecta, which leaves the tophat jet unaltered, and the presence of light ejecta, which produce a fast, low-energy cocoon. The former configuration gives a peak time that corresponds with the time at which the jet emission comes into the line of sight, at about 100 s in Figure 4. The latter, instead, peaks early, since the deceleration radius scales as r dec ∝ E c 1/3 Γ −2/3 c , yielding a deceleration time for material along the line of sight t dec ∝ E c 1/3 Γ −8/3 c . This yields a scaling for the cocoon-afterglow peak luminosity L pk ∝ E c 2/3 Γ 8/3 c , showing that the peak luminosity of the cocoon grows quickly for increasing Γ c , unless its energy vanishes.
In some cases, however, there are some degeneracies. For example, a careful comparison of the leftmost and rightmost panels of Figure 4 reveals similar light curves at both radio and X-ray frequencies. This is not surprising since both the delay δt and the ejecta velocity v w control the radius of the ejecta in which the jet is propagating. Breaking these degeneracies will require the use of additional data. For example, one could use the shock-breakout model of Kasliwal et al. (2017) to constrain the outer radius of the ejecta from the brightness of the prompt γ-ray emission. In addition, constraints on δt can be set since it has to be smaller than the time delay between the detection of gravitational waves and of γ-ray photons. Additional constraints could come form the kilonova properties and from the gravitational wave detection (e.g., constraints on the inclination of the system, the chirp mass, the individual masses of the two NS). Thus one can envisage that afterglow modeling with the appropriate priors from the complete dataset will be able to set meaningful constraints on the ejecta properties.

EJECTA PROPERTIES AND THE EQUATION OF STATE OF NEUTRON STARS
As already mentioned in §1 to motivate this work, the connection between the amount and velocity of tidally disrupted mass and observable properties of the mergers (in particular in the electromagnetic spectrum) is especially interesting in light of the fact that it contains information on the EoS of the NS, as demonstrated by several groups via GRMHD simulations (e.g. Shibata et al. 2005;Kiuchi et al. 2010;Rezzolla et al. 2010;Hotokezaka et al. 2011;Bauswein et al. 2013;Hotokezaka et al. 2013;Sekiguchi et al. 2016;Lehner et al. 2016;Ruiz et al. 2016;Kawamura et al. 2016;Radice et al. 2016;Ciolfi et al. 2017;Radice et al. 2018;Foucart et al. 2019) Most of the tidally disrupted mass remains bound and circularizes into a hot torus, whose accretion onto the remnant compact object formed from the merger powers a relativistic jet. The total energy in the jet, proportional to the total accreted mass, carries information on the EoS; this connection was exploited by  with the sample of the short GRBs available at the time, to draw inferences on the likelihood of various EoS.
Here we are focusing on signatures of the ejecta mass in the light curves. This high-velocity material constitutes a (small) fraction of the total mass which is tidally disrupted during the merger. More specifically, this dynamical ejecta is partly due to matter ejected (and unbound) by the tidal disruption of the less massive NS, and partly by matter ejected (and unbound) during merger due to the formation of shocks. Analogously to the total matter which is tidally disrupted, this small fraction of matter, and its velocity, carries the imprint of the EoS of the NS, in addition to being dependent on the masses and mass ratios of the NSs.
Several studies have specifically focused on the dynamical ejecta (e.g. Bauswein et al. 2013;Hotokezaka et al. 2013;Sekiguchi et al. 2016;Lehner et al. 2016;Radice et al. 2018), and computed, via GRMHD simulations, its properties for a set of representative EoS, NS masses, and mass ratios. The simulations show that larger mass ratios tend to result in the partial disruption of the smaller star during the merger, and hence a smaller amount of shocked material, since the stars merge less violently. On the hand, the mass asymmetery also produces more massive tidal out- flows. A larger amount of ejecta further tends to correlate with softer EoSs.
Dietrich & Ujevic (2017) collected a large set of numerical relativity simulations performed by different groups, and combined them to determine the best fits for the ejecta mass, velocity and kinetic energy. In particular, they find for the ejecta mass where M i and M * i indicate the gravitational and baryonic mass of star 'i', respectively, and C i = 2 × 10 33 × GM i /c 2 R i the compactness, with R i being the radius. All masses are in units of solar masses. The symbol (1 ↔ 2) means that the expression in brackets should be repeated exchanging the subscripts. The best fit parameters were found to be: a = −1. 35695, b = 6.11252, c = −49.43355, d = 16.1144, n = −2.5484.
Similarly, the ejecta velocity was found to be best fit by the expression where the best fits values for the constants are: a = −0.219479, b = −2.67385, c = 0.444836 . Note that this same parametrization was used Radice et al. (2018) al. to fit the results of their simulations; however they found systematically lower ejecta mass compared to the best fits by  with the best fit values for the parameters differing by a factor of ∼ 50% − 100%. These variations among different groups should be kept in mind in interpreting results at this stage. However, it is conceavable to expect that simulation results among different groups will eventually converge, and hence the dynamical ejecta will become an additional strong probe of the EoS.
A visual representation of the dependence of the ejecta properties on the EOS of NS is offered in Figure 5, which displays the ejecta mass (top panels) and velocity (bottom panels) for three representative equations of state of neutron star matter 2 . These quantities are shown as a function of the primary mass M 1 and the chirp mass M, which is extremely well constrained by GW observations alone (e.g., The LIGO Scientific Collaboration et al. (2018)). In the figures, white areas correspond to regions that are either unphysical (M > M 1 /2 1/5 ), for which Eq. 22 or 23 yield negative results, or for which the mass-radius relationship is not allowed by that equation of state. For the sake of example, a white horizontal line shows the parameter region allowed by the constraints on GW170817: M = 1.186 ± 0.001 and 1.36 ≤ M 1 ≤ 1.89 M . In this case, being able to set a constraint on the ejecta velocity would allow us to determine the correct EoS, since the three cases shown predict different velocities for the ejecta. We note that the relatively large ejecta masses inferred from the kilonova modeling (Kasen et al. 2017) are consistent with the late rise of the afterglow radiation, both in X-rays and in radio observed in the followup to GRB170817A (Abbott et al. 2017). Therefore, early multi-band monitoring will be especially important in order to exploit the sensitivity of the afterglow light curves to the ejecta properties.

SUMMARY AND CONCLUSIONS
The detection of electromagnetic radiation from a binary NS merger triggered by GWs has allowed an unprecedented view into the inner workings of a relativistic jet, likely powered by accretion onto the remnant object, and propagating onto the ejecta from the merger.
Since the LIGO horizon is much smaller than the Swift one, events triggered in GWs will be generally much closer than the ones triggered in γ-rays. Hence their afterglows, especially in the radio band which is not contaminated by the kilonova emission, will generally be much brighter than for the more distant short GRBs triggered in γ-rays. As shown by the case of GW170817, this allows for a detailed monitoring of the afterglow light curves in several bands by a variety of instruments, and hence a detailed shape reconstruction of the light curve.
In this paper we have studied the structure of the jet-cocoon which forms as a result of the jet expanding into the fastmoving ejecta produced by the NS-NS merger. More specifically, we have explored how this structure varies with the density, velocity, and radial extent of the ejecta. We have then used this structure at the radiative stage to compute the shape of the afterglow in the X-ray and radio bands (the optical is likely dominated by the kilonova, hence may have less diagnostic power), and explored its dependence on the ejecta properties, as well as on the time delay between merger and onset of the jet.
We identified features in the light curves which are especially distinctive of a jet-cocoon configuration, where the cocoon is generated by the interaction of the jet with the dynamical ejecta. A more pronounced cocoon yields an additional, earlier bump in the light curve, compared to the singly-peaked light curve of a pure relativistic jet. The rise time of the light curves is also very sensitive to the amount of ejecta, with larger values yielding delayed rises.
Therefore, a detailed monitoring of the afterglow light curve, especially in the cleaner X-ray and radio bands, can help constrain the properties of the dynamical ejecta from the merger event. Since this material is sensitive to the equation of state of dense matter (in addition to being important for production of heavy elements), detailed monitoring of afterglow light curves of GW-detected binary NS mergers (and NS-BH alike) has the power to help probing the physics of the merger events and the equation of state of dense nuclear matter.
We thank Riccardo Ciolfi, Brian Metzger, David Radice, and Om Sharam Salafia for friutful discussions and constructive comments on an early draft of this manuscript. DL acknowledges support from NASA ATP grant NNX17AK42G, NASA Fermi GI grant 80NSSC19K0330, and Chandra grant TM9-20002X. RP acknowledges support from the NSF under grant AST-1616157. The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.