# Jet Propagation in Expanding Medium for Gamma-Ray Bursts

Hamid Hamidani1 and Kunihito Ioka1
1 Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502,
Japan
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The binary neutron star (BNS) merger event GW170817 clearly shows that a BNS
merger launches a short Gamma-Ray Burst (sGRB) jet. Unlike collapsars, where
the ambient medium is static, in BNS mergers the jet propagates through the
merger ejecta that is expanding outward at substantial velocities ($\sim
0.2c$). Here, we present semi-analytic and analytic models to solve the
propagation of GRB jets through their surrounding media. These models improve
our previous model by including the jet collimation by the cocoon self-
consistently. We also perform a series of 2D numerical simulations of jet
propagation in BNS mergers and in collapsars to test our models. Our models
are consistent with numerical simulations in every aspect (the jet head
radius, the cocoon’s lateral width, the jet opening angle including
collimation, the cocoon pressure, and the jet-cocoon morphology). The energy
composition of the cocoon is found to be different depending on whether the
ambient medium is expanding or not; in the case of BNS merger jets, the cocoon
energy is dominated by kinetic energy, while it is dominated by internal
energy in collapsars. Our model will be useful for estimating electromagnetic
counterparts to gravitational waves.

###### keywords:

gamma-ray: burst – hydrodynamics – relativistic processes – shock waves – ISM:
jets and outflows – stars: neutron – gravitational waves

††pubyear: 2020††pagerange: Jet Propagation in Expanding Medium for Gamma-Ray
Bursts–C.3

##  1 Introduction

Observation of the gravitational wave signal from the binary neutron star
(BNS) merger event GW170817 by the Laser Interferometer Gravitational-Wave
Observatory (LIGO) and the Virgo Consortium (LVC) [Abbott et al. 2017a], and
the follow-up observation campaign across the electromagnetic (EM) spectrum
marked the dawn of the era of multi-messenger astronomy (Abbott et al. 2017b).
One of the most important findings was the association of GW170817 with the
prompt emission of the short Gamma-Ray Burst (sGRB) sGRB 170817A $\sim 1.7$ s
after the GW signal (Abbott et al. 2017b). Clear evidences of a relativistic
jet have also been obtained from radio observations at later times (Mooley et
al. 2018). As GRB observations show, relativistic jets are very important in
time-domain astronomy, especially because of the EM emission over a wide
spectrum (i.e., the prompt and the afterglow emission).

This scenario, that the BNS merger powers a relativistic jet, was
theoretically suggested for sGRBs in the past (Paczynski 1986; Goodman 1986;
Eichler et al. 1989). Recent studies of numerical relativity show that after
the merger, a central engine is formed and is surrounded by $\sim
10^{-3}-10^{-2}M_{\odot}$ of matter that has been ejected during the merger,
referred to as “the dynamical ejecta”, that expands at a speed of $\sim 0.2c$,
where $c$ is the light speed (Shibata 1999; Shibata & Uryū 2000; Hotokezaka et
al. 2013; Bauswein et al. 2013; Radice et al. 2018a; etc.). Also, later on
($<10$ s after the merger), matter gets ejected from the torus surrounding the
merger remnant in the form of wind (Siegel & Metzger 2017; Radice et al.
2018b; Fujibayashi et al. 2018, 2020; Nedora et al. 2020; Ciolfi & Vijay
Kalinani 2020; etc.).

Multi-messenger observations of the BNS merger event GW170817 indicates that
the relativistic jet was launched from the central engine within $\sim 1.7$ s
after the merger (more precisely, within $\sim 1.3$ s after the merger
according to Hamidani et al. 2020, see their Figure 9; within $0.4$ s
according to Lazzati et al. 2020). Therefore, the jet must have propagated
through the dense surrounding medium (i.e., the dynamical ejecta), and
successfully broke out of it for the sGRB 170817A to be emitted (as it has
been observed). This is because the jet outflow is not observable unless it
propagates up to the outer edge of the medium and eventually breaks out of it;
as it is the case for long GRBs, where the relativistic jet propagates through
the stellar envelope of a massive star.

During its propagation, the jet continuously injects energy into the expanding
ejecta material. This produces the hot cocoon. The cocoon immediately
surrounds the jet, interacts with it, and collimates it. Although this phase,
where the jet is confined inside the ejecta, is short, it is critical, as it
shapes the jet (and the cocoon) structure (see Gottlieb et al. 2020). After
the breakout, the jet (and the cocoon) is the source of different EM
counterparts over a wide band, and it is the key to interpreting them (Nakar &
Piran 2017; Lazzati et al. 2017b; Gottlieb et al. 2018a; Lazzati et al. 2018;
Piro & Kollmeier 2018; Nakar et al. 2018; Ioka & Nakamura 2018; etc.). Also,
this connection between the jet and the EM counterparts allows us to make use
of observational data to extract crucial information and better understand the
phenomenon of sGRBs (e.g., the jet angular structure, Troja et al. 2019;
Takahashi & Ioka 2020; the property of the jet outflow, Ioka & Nakamura 2019;
the central engine, Gill et al. 2019; Hamidani et al. 2020; Lazzati et al.
2020; Salafia & Giacomazzo 2020; the physics of neutron density matter,
Lazzati & Perna 2019; the viewing angle, Nakar & Piran 2020; etc.). Therefore,
the jet propagation through the ejecta surrounding the BNS merger remnant,
until the breakout, is a key process in sGRBs (as it is in collapsars and long
GRBs).

The propagation of astrophysical jets through dense ambient media has been the
subject of intensive theoretical works; mostly, in the context of Active
Galactic Nuclei (AGNs) and collapsars (Begelman & Cioffi 1989; Martí et al.
1997; MacFadyen & Woosley 1999; Matzner 2003; Bromberg et al. 2011; Mizuta &
Ioka 2013; Harrison et al. 2018; etc.). One critical difference in the context
of BNS mergers is that the ambient medium expands at substantial velocities
(Hotokezaka et al. 2013), while it is static in AGNs and collapsars (MacFadyen
& Woosley 1999). This further complicates the problem of modeling the jet
propagation in BNS mergers.

There have been an increasing number of studies dedicated to solving jet
propagation in BNS mergers through numerical simulations, especially after the
discovery of GW170817 (Nagakura et al. 2014; Murguia-Berthier et al. 2014;
Duffell et al. 2015; Bromberg et al. 2018; Duffell et al. 2018; Lazzati et al.
2017b; Gottlieb et al. 2018a; Gottlieb et al. 2018b; Xie et al. 2018;
Nathanail et al. 2020; Gottlieb et al. 2020; etc.). However, the subject is
still far from being well understood.

Using ideas from the modeling of the jet-cocoon in collapsars (e.g., Bromberg
et al. 2011; Mizuta & Ioka 2013; Harrison et al. 2018), several studies
presented analytic modeling of the jet-cocoon in an expanding medium (Lazzati
et al. 2017a; Margalit et al. 2018; Duffell et al. 2018; Matsumoto & Kimura
2018; Gill et al. 2019; Salafia et al. 2020; Hamidani et al. 2020; Lyutikov
2020; Beniamini et al. 2020; etc.). Although some of these works offer
promising results, many of them overlooked important aspects, such as the jet
collimation, the expansion of the ejecta and its effect on the cocoon (e.g.,
on the cocoon pressure and on the cocoon radius), etc. And still, there is no
analytic model for the jet propagation in an expanding medium that is simple
to use, and robust at the same time (which is necessary to investigate other
related topics, such as the emission from the cocoon).

This work presents analytic modeling of jet propagation in an expanding
medium. This model is an upgrade of the model presented in Hamidani et al.
(2020). The main improvement is that the jet collimation (i.e., the jet
opening angle) is calculated in a self-consistent manner by using the cocoon
pressure; rather than relying on the assumption of a constant opening angle,
and on the free parameter $f_{j}$. Another addition is the semi-analytic model
presented here, where results are found through numerical integration of
ordinary differential equations, relying less on approximations.

The aim of our work is to present a physical model that accurately describes
the jet-cocoon system in an expanding medium, in consistency with numerical
simulation. A crucial point in our study is that our modeling is based on
rigorous analysis of the jet-cocoon system in numerical simulations (in both
expanding and static media), which shows that the expansion of the medium does
intrinsically affect the jet-cocoon (e.g., the energy composition of the
cocoon, the expansion velocity of the cocoon, etc.). We show that the jet-
cocoon system, in an expanding medium, can be described by a set of equations
that can be solved numerically (referred to as the “semi-analytic" solution).
We also show that, with some reasonable approximations, the system of
equations can be simplified and solved analytically (the “analytic" solution).
Both solutions are rigorously compared to the results from the numerical
simulations and found consistent.

This paper is organized as follows. In Section 2, physical modeling of the
jet-cocoon system in both expanding and static media is presented, and two
(semi-analytic and analytic) solutions are derived. In Section 3, numerical
simulations are presented and compared to both solutions. A conclusion is
presented in Section 4.

##  2 The jet-cocoon physical model

The jet-cocoon model presented here is an upgrade of previous models, in
particular in Bromberg et al. (2011); Mizuta & Ioka (2013); Harrison et al.
(2018); Hamidani et al. (2020); etc. Unlike previous models, this model allows
us to treat the jet collimation by the cocoon in expanding media. Therefore,
this model can be applied not only to the case of collapsar jets (where the
jet propagates through the static stellar envelope of a massive star) but also
to the case of BNS merger jets (where the jet propagates through the expanding
dynamical ejecta). This model is also an upgrade of the model presented
Hamidani et al. (2020); it takes into account the cocoon and its pressure on
the jet (i.e., collimation), hence allowing to derive the jet opening angle in
a self-consistent manner.

###  2.1 Jump conditions

The jet head dynamics in a dense ambient medium can be determined by the shock
jump conditions at the jet’s head (e.g., Begelman & Cioffi 1989; Martí et al.
1997; Matzner 2003):

$\displaystyle
h_{j}\rho_{j}c^{2}(\Gamma\beta)_{jh}^{2}+P_{j}=h_{a}\rho_{a}c^{2}(\Gamma\beta)_{ha}^{2}+P_{a},$
(1)

where $h$, $\rho$, $\Gamma=(1-\beta^{2})^{-1/2}$, and $P$ are enthalpy,
density, Lorentz factor (with $\beta$ being the velocity normalized by the
light speed), and pressure of each fluid element, all measured in the fluid’s
rest frame. The subscripts $j$, $h$, and $a$ refer to three domains: the
relativistic jet, the jet head, and the cold ambient medium, respectively.
Typically, both $P_{a}$ and $P_{j}$ in equation (1) are negligible terms;
hence, we can write the jet head velocity as (for more details see Ioka &
Nakamura 2018; Hamidani et al. 2020):

$\displaystyle\beta_{h}=\frac{\beta_{j}-\beta_{a}}{1+\tilde{L}^{-1/2}}+\beta_{a},$
(2)

where $\tilde{L}$ is the ratio of energy density between the jet and the
ejecta,
$\tilde{L}=\frac{h_{j}\rho_{j}\Gamma_{j}^{2}}{h_{a}\rho_{a}\Gamma_{a}^{2}}$,
which can be approximated as:

$\displaystyle\tilde{L}\simeq\frac{L_{j}}{\Sigma_{j}(t)\rho_{a}c^{3}\Gamma_{a}^{2}},$
(3)

where $L_{j}$ is the (true) jet luminosity (per one jet),
$\Sigma_{j}(t)=\pi\theta_{j}^{2}(t)r_{h}^{2}(t)$ is the jet head cross
section, with $r_{h}(t)$ being the jet head radius
[$r_{h}(t)=\int^{t}_{t_{0}}\beta_{h}dt+r_{0}$], and $\theta_{j}(t)$ being the
jet opening angle (Matzner 2003; Bromberg et al. 2011).

###  2.2 Main approximations

Here, in our analytic modeling, we consider a similar set of approximations to
that in Hamidani et al. (2020). In summary:

  1. 1.

The analytic treatment presented here is limited to the case of a non-
relativistic jet head where:

$\tilde{L}\ll(1-\beta_{a})^{2}.$ (4)

  2. 2.

In BNS mergers, during the merger, matter is dynamically ejected from the
system (i.e., the dynamical ejecta) and surrounds the later formed central
engine. After the merger, another plausible source of mass ejection is the
wind. The mass of matter driven by the wind could be substantial in early time
(e.g., $\sim 0.01M_{\odot}$ in the case of magnetically driven wind; see
Ciolfi & Vijay Kalinani 2020), however, the launch time (of the wind) may be
later than the jet depending on the type of wind (see Fujibayashi et al. 2018
for the case of viscous wind). For simplicity, we consider the case where the
central engine is only surrounded by the dynamical ejecta. Note that, here,
the ambient medium (often referred to as “the ejecta”) is defined as the
medium surrounding the central engine through which the jet propagation takes
place; regardless of whether this inner region is gravitationally bound to the
central engine or not. Hence, for simplicity, the total mass of the ambient
medium, $M_{a}$, is defined by the ambient medium’s density through which the
jet head propagates $\rho_{a}(r,t)$ [i.e., in the polar direction]111Note
that, in the case of BNS mergers, the density throughout the dynamical ejecta
is angle dependent, with the density in the polar region being much lower than
that near the equatorial region (see Figure 8 in Hamidani et al. 2020). This
effect results in equation (5) giving a mass, $M_{a}$, of $\sim
0.002M_{\odot}$ if accounting for a total dynamical ejecta mass of $\sim
0.01M_{\odot}$ (for more details see Hamidani et al. 2020). Note that, this
value for the mass $M_{a}$ can be scaled-up to account for contribution form
the wind.:

$M_{a}=\int_{r_{0}}^{r_{m}(t)}4\pi r^{2}\rho_{a}(r,t){\rm d}r,$ (5)

where $r_{0}$ is the inner boundary of the ambient medium, which is of the
order of $10^{6}-10^{7}$ cm, and $r_{m}(t)$ is the outer radius of the ambient
medium. This definition of $M_{a}$ is also used for the collapsar case.

  3. 3.

In the case of an expanding medium (BNS merger), the approximation of a
homologous expansion is used (see Figure 8 in Hamidani et al. 2020). Hence,
the radial velocity of the ambient medium as a function of radius $r$ and time
$t$ is approximated as:

$v_{a}(r,t)=\left[\frac{r}{r_{m}(t)}\right]v_{m},$ (6)

with $r_{m}(t)$ being the outer radius of the ambient medium, and $v_{m}$
being the maximum velocity of the ambient medium at the radius $r_{m}(t)$. In
the case of collapsars, the velocity is negligible, i.e., $v_{m}=0$.

  4. 4.

The ambient medium’s density profile is approximated to a power-law function,
with “$n$" being its index. Hence, considering the homologous expansion of the
ambient medium, the density can be written as:

$\rho_{a}(r,t)=\rho_{0}\left[\frac{r_{0}}{r}\right]^{n}\left[\frac{r_{m,0}}{r_{m}(t)}\right]^{3-n},$
(7)

where $\rho_{0}=\rho_{a}(r_{0},t_{0})=\left[\frac{M_{a}}{4\pi
r_{0}^{n}}\right]\left[\frac{3-n}{r_{m,0}^{3-n}-r_{0}^{3-n}}\right]$, and
$r_{m,0}=r_{m}(t_{0})$, with $t_{0}$ being the jet launch time. This
expression is simpler in the collapsar case, where $r_{m}(t)=r_{m,0}$ ($\equiv
r_{m}$). Also, the value $n=2$ is assumed for the density profile of the
ambient medium, for both the BNS merger case and the collapsar case. It should
be noted that this is a simplification as, ideally, $n\sim 2-3.5$ in the BNS
merger case (see Figure 8 in Hamidani et al. 2020), and $n\sim 1.5-3$ in the
collapsar case (see Figure 2 in Mizuta & Ioka 2013). Also, it should be noted
that the analytic modeling presented here is limited for the case $n<3$
(Bromberg et al. 2011; Hamidani et al. 2020).

  5. 5.

The pressure in the cocoon, $P_{c}$, is dominated by radiation pressure.
Hence, it can be written as:

$P_{c}=\frac{E_{i}}{3V_{c}},$ (8)

where $E_{i}$ is the cocoon’s internal energy and $V_{c}$ is the cocoon’s
volume.

  6. 6.

Based on rigorous analysis of the cocoon in numerical simulations, we suggest
that the cocoon’s shape is better approximated to an ellipsoidal (see Hamidani
et al. 2020; also see Figure 5 below); where the ellipsoid’s semi-major axis
and semi-minor axis at a time $t$ are $\frac{1}{2}r_{h}(t)$ and $r_{c}(t)$,
respectively, with $r_{c}(t)$ being the cocoon’s lateral width (from the jet
axis) at the radius $\frac{1}{2}r_{h}(t)$ [see also equation (16)]. Hence, the
volume of the cocoon (in one hemisphere) can be written as:222Ideally, the jet
volume should be subtracted from the above expression of $V_{c}$ to give a
more accurate expression of the cocoon volume. However, as long as the jet
opening angle is not very large (as it is the case here), the jet volume can
be neglected.

$V_{c}=\frac{2\pi}{3}r_{c}^{2}(t)r_{h}(t).$ (9)

Note that this presents one of the differences compared to previous works –
typically assuming a cylindrical cocoon shape (e.g., Bromberg et al. 2011;
Mizuta & Ioka 2013; Salafia et al. 2020).

  7. 7.

As previously explained in Mizuta & Ioka (2013), and more rigorously in
Harrison et al. (2018), the analytic description of $\tilde{L}$ in equation
(3) needs to be calibrated by numerical simulations. The parameter $N_{s}$ is
introduced to calibrate the analytic value of $\tilde{L}$ so that:

$\tilde{L}_{c}=N_{s}^{2}\tilde{L},$ (10)

where $\tilde{L}_{c}$ is the calibrated counterpart of $\tilde{L}$. Here, the
value of $N_{s}$ for the analytic (or semi-analytic) solution is chosen so
that the analytic (or semi-analytic) breakout time is calibrated to the
breakout time measured in numerical simulations (see Table 1).

As previously noted in Hamidani et al. (2020) (also see Mizuta & Ioka 2013;
Harrison et al. 2018), $N_{s}$ accounts for the part of meandering energy,
without contributing to the forward jet head motion. However, it should be
noted that value of $N_{s}$ should be dependent on the parameter space (in
particular on the value of $\tilde{L}$, see Figure 12 in Harrison et al. 2018;
and $v_{m}$, see Appendix C). Therefore, the values of $N_{s}$ used here
should be limited to the parameter space used here, and should not be taken at
face value (for more details on $N_{s}$ refer to Appendix C).

Following the introduction of the calibration coefficient $N_{s}$, $\tilde{L}$
is substituted by $\tilde{L}_{c}$ in equation (2), and with $\beta_{j}\simeq
1$, the jet head velocity can be written as:

$\displaystyle\beta_{h}=\left[(1-\beta_{a})(1+\tilde{L}_{c}^{1/2})^{-1}\right]\tilde{L}_{c}^{1/2}+\beta_{a},$
(11)

where $\tilde{L}_{c}$ can be found from equations (3) and (10). Given the jet
luminosity $L_{j}$, the ambient medium’s velocity $\beta_{a}$
[$=v_{a}(r,t)/c$; see equation (6)], and the density $\rho_{a}(r,t)$ [see
equation (7)], the only unknown quantity for $\tilde{L}_{c}$ (i.e.,
$\beta_{h}$) to be determined is the jet head cross-section $\Sigma_{j}(t)$.
The jet head opening angle $\theta_{j}(t)$; and hence $\Sigma_{j}(t)$; will be
determined in Section 2.3 by considering the collimation of the jet by the
cocoon.

The jet head velocity [i.e., $\beta_{h}$ in equation (11)] will be solved in
two different ways: Semi-analytically and analytically (details are given in
Sections 2.4 and 2.5, respectively). In the semi-analytic solution, the
expression of $\beta_{h}$ in equation (11) is used as it is, and is solved
through numerical integration. In the analytic solution, the above expression
of $\beta_{h}$ is further approximated so that it is solved analytically (see
Section 2.5.1).

###  2.3 The cocoon and jet collimation

####  2.3.1 The system of equations

We follow the same treatment of Bromberg et al. (2011). The unshocked jet’s
height $\hat{z}$ can be written as a function of the jet luminosity $L_{j}$
and the cocoon’s pressure $P_{c}$:

$\displaystyle\hat{z}=\sqrt{\frac{L_{j}}{\pi cP_{c}}}+z_{*}.$ (12)

With $r_{in}$ being the radius at which the jet is injected into the medium,
$z_{*}=\max[r_{in},z(P_{c}=P_{j0})]$ is the radius at which the pressure of
the injected jet and the pressure of the cocoon are balanced; beyond $z_{*}$
the pressure of the cocoon is higher than the pressure of the injected jet
(Bromberg et al. 2011). In our simulations, $z_{*}$ is typically of the same
order of $r_{in}$, hence, for simplicity, we take $z_{*}\approx r_{in}$.

At a certain time $t$, the jet is uncollimated if the jet head’s radius,
$r_{h}(t)$, is below $\hat{z}/2$, and collimated if it is beyond $\hat{z}/2$
(see Figure 2 in Bromberg et al. 2011). Hence, the jet head’s cross-section
can be found for the two modes as follows:

$\Sigma_{j}(t)=\begin{cases}\pi r_{h}^{2}(t)\theta_{0}^{2}&\text{if
$r_{h}(t)<\hat{z}/2$ (uncollimated jet)},\\\ \pi
r_{h}^{2}(t)\theta_{j}^{2}(t)&\text{if $r_{h}(t)>\hat{z}/2$ (collimated
jet)},\end{cases}$ (13)

where $\theta_{0}$ is the initial opening angle of the jet333The initial
opening angle is given by $\theta_{0}\approx\theta_{inj}+1/\Gamma_{0}$ where
$\theta_{inj}$ is the opening angle of the injected jet at $t=t_{0}$ and
$r=r_{in}$, and $\Gamma_{0}$ is its initial Lorentz factor., and
$\theta_{j}(t)$ is the opening angle of the jet head at a given time $t$.

Since the cocoon shape is approximated to an ellipsoidal [see (vi) in Section
2.2], $r_{c}(t)$ is the cocoon’s lateral width at the radius
$\frac{1}{2}r_{h}(t)$. $r_{c}(t)$ is determined by integrating the lateral
velocity, $\beta_{\perp}$, with which the cocoon expands into the ambient
medium at the radius
$\sqrt{[r_{h}(t)/2]^{2}+r_{c}^{2}(t)}\approx\frac{1}{2}r_{h}(t)$ [since
$r_{h}(t)\gg r_{c}(t)$; see Figures 2, 3, and 5]. At this radius, since the
ambient medium’s velocity $v_{a}(r_{h}/2)$ is $\leq v_{m}/2$, and considering
the value of $v_{m}$ (see Table 1), $\Gamma_{a}(r_{h}/2)$ is $\approx 1$ and a
non-relativistic treatment is reasonable. $\beta_{\perp}$ is therefore
determined by the ram pressure balance between the cocoon and the ambient
medium at the radius $\frac{1}{2}r_{h}(t)$, giving:

$\displaystyle
P_{c}\approx\rho_{a}(r_{h}/2,t)c^{2}[\beta_{\perp}-\beta_{a,\perp}]^{2},$ (14)

where $\beta_{a,\perp}$ is the ambient medium’s expansion velocity [see
equation (6)] in the lateral direction:

$\displaystyle\beta_{a,\perp}=\left[\frac{r_{c}(t)}{r_{m}(t)}\right]\frac{v_{m}}{c}.$
(15)

In summary, the equations describing the jet-cocoon system can be found as
follows:

$\displaystyle\frac{dr_{c}(t)}{dt}=$ $\displaystyle c\beta_{\perp},$ (16)
$\displaystyle\beta_{\perp}=$
$\displaystyle\sqrt{\frac{P_{c}}{{\rho}_{a}(r_{h}/2,t)c^{2}}}+\left[\frac{r_{c}(t)}{r_{m}(t)}\right]\frac{v_{m}}{c},$
(17) $\displaystyle P_{{c}}=$
$\displaystyle\>\>\>\>\>\>\frac{E_{i}}{3\>V_{{c}}}\>\>\>\>\>\>\>=\eta\frac{L_{j}\left(1-\langle{\beta_{h}}\rangle\right)\>(t-t_{0})}{2\pi
r_{c}^{2}(t)r_{{h}}(t)},$ (18) $\displaystyle\Sigma_{j}(t)=$ $\displaystyle\pi
r_{h}^{2}(t)\theta_{j}^{2}(t)=\frac{L_{j}\theta_{0}^{2}}{4cP_{c}},$ (19)

where $P_{c}$ and $V_{c}$ are defined as in equations (8) and (9),
respectively; and
$\left<\beta_{h}\right>=\frac{1}{c}\frac{r_{h}(t)-r_{0}}{t-t_{0}}$ is the
time-averaged jet head velocity, which is a term that takes into account the
fact that a part of the injected energy
[$=L_{j}\langle{\beta_{h}}\rangle(t-t_{0})$] is contained in the jet and does
not make its way into the cocoon instantly. The last equation (19) is
determined by the pressure balance between the post-collimated jet and the
cocoon (Bromberg et al. 2011).

The expression of $\beta_{\perp}$ [and eventually
$r_{c}(t)=\int_{t_{0}}^{t}c\beta_{\perp}dt+r_{0}\theta_{0}$] here is different
from the original collapsar case where the medium is static (Bromberg et al.
2011; Harrison et al. 2018); it is instead applicable to both the case of
static medium and the case of expanding medium. The term
$\left[\frac{r_{c}(t)}{r_{m}(t)}\right]\frac{v_{m}}{c}$ in equation (17) is
new and is the result of the homologous expansion of the medium. It is worth
mentioning that the term
$\left[\frac{r_{c}(t)}{r_{m}(t)}\right]\frac{v_{m}}{c}$ is far more dominant
(in $\beta_{\perp}$) over the term
$\sqrt{\frac{P_{c}}{\rho_{a}(r_{h}/2,t)c^{2}}}$ in the case of an expanding
medium as in BNS mergers, and hence it is important.

####  2.3.2 The parameters $\eta$ and $\eta^{\prime}$

$\eta$ in equation (18) is a parameter that expresses the fraction of internal
energy in the total energy delivered into the cocoon (by the engine and
through the jet) at a given time $t$ (Bromberg et al. 2011; Mizuta & Ioka
2013). It takes values between 0 and 1, and it can be expressed as:

$\eta=\frac{3P_{c}V_{c}}{L_{j}\left(1-\langle{\beta_{h}}\rangle\right)\>(t-t_{0})},$
(20)

with $E_{i}=3P_{c}V_{c}$ [see equation (8)]. For convenience, we define the
parameter $\eta^{\prime}=\eta[1-\langle{\beta_{h}}\rangle]$; it relates to the
fraction of internal energy in the cocoon out of the total energy delivered by
the central engine, at a given time $t$. Hence:

$\eta^{\prime}=\frac{3P_{c}V_{c}}{L_{j}\>(t-t_{0})}.$ (21)

$\eta$ and $\eta^{\prime}$ can be easily deduced from numerical simulations by
measuring both $P_{c}$ and $V_{c}$, or by measuring the internal energy in the
cocoon $E_{i}$. In Section 3.2, using numerical simulations’ results, we will
show that, on average, $\langle{\eta^{\prime}}\rangle\sim 1/2$ for the
collapsar case, and $\langle{\eta^{\prime}}\rangle\sim 1/4$ for the BNS merger
case [see Figure 1 and equation (53)], where:

$\langle{\eta^{\prime}}\rangle=\frac{1}{t_{b}-t_{0}}\int_{t_{0}}^{t_{b}}\eta^{\prime}dt.$
(22)

These fiducial values will be adopted to solve the jet head motion (see Table
1).

###  2.4 The semi-analytic solution

Here, the system of equations [equations (3), (10), (11), (12), (13), (16),
(17), (18), and (19)] is solved though numerical integration. At every time
step, the time is updated (from $t$ to $t+dt$, where $dt$ is sufficiently
small). The density ${\rho}_{a}(r_{h}/2,t)$ in equation (17) is calculated
using equation (7). Then, using equation (18) the pressure is calculated; the
parameter $\eta^{\prime}$ [as defined in equation (21)] is represented by its
time-averaged value $\langle{\eta^{\prime}}\rangle$ [see equation (22)] as
measured in numerical simulations (see Table 1 for the values of
$\langle{\eta^{\prime}}\rangle$ used). Next, $\beta_{\perp}$ is derived using
equation (17). The jet head’s cross-section and opening angle are found by
calculating $\hat{z}$ first, using (12), and then determining the collimation
mode and the opening angle of the jet, using equation (13) together with
equation (19). $\tilde{L}_{c}$ is then calculated using equations (3) and
(10). Finally, at the end of each time step, the jet head radius $r_{h}(t)$,
the cocoon’s lateral width $r_{c}(t)$, and the cocoon’s volume $V_{c}$, for
the next time step are calculated using equations (11), (16) and (9),
respectively. These processes are repeated until the jet breaks out of the
ambient medium [i.e., the following condition is met: $r_{h}(t)\geqslant
r_{m}(t)$].

###  2.5 The analytic solution

Here, the system of equations (3), (10), (11), (13), (16), (17), (18), and
(19) [in Sections 2.1, 2.2, and 2.3] is simplified using several additional
approximations, and then solved analytically.

In summary, the jet head’s velocity [equation (11)] is simplified to equation
(23), which can be written as a function of $t$, $r_{h}(t)$, and
$\theta_{j}(t)$ using equations (3), (10), and (13) [see Section 2.5.1]. The
expression of the cocoon’s lateral width, $r_{c}(t)$, is simplified from
equation (16) to equation (30) [$\langle{\chi}\rangle$ can be found with
equations (26) and (28); see Section 2.5.2], and with equation (32) the
expression of the cocoon pressure, $P_{c}$, is derived analytically in Section
2.5.3 [in equation (34) as a function of $t$ and $r_{h}(t)$]. Next, equation
(19) is used to find the analytic expression of the jet opening angle
$\theta_{j}(t)$ [equation (35) as a function of $r_{h}(t)$ and $t$], which
allows us to derive an analytically solvable equation of motion of the jet
head [equation (36)], and to determine the solution, $r_{h}(t)$, as a function
of the initial parameters and $t$ (see Section 2.5.4).

The same logic can be used in the collapsar, and the equation of motion of the
jet head can be found accordingly [equation (46); see Section 2.5.5].

For reference, Table 1 presents a summary of the relevant parameters and the
values they take.

####  2.5.1 Approximated jet head velocity $\beta_{h}$

In the analytic solution, two additional approximations are used for the jet
head velocity. Firstly, in the case of BNS mergers where the medium is
expanding (i.e., $\Gamma_{a}>1$), the term $\frac{1}{\Gamma_{a}^{2}}$ in the
expression of $\tilde{L}_{c}$ in equation (3) is considered as constant and is
absorbed into $N_{s}$. Secondly, in the analytic solution, the term
$\left[(1-\beta_{a})(1+\tilde{L}_{c}^{1/2})^{-1}\right]$ in equation (11) is
also approximated as, roughly, constant over time and is also effectively
absorbed into the calibration coefficient $N_{s}$. The result is the following
expression:

$\displaystyle\beta_{h}\approx\tilde{L}_{c}^{\frac{1}{2}}+\beta_{a}.$ (23)

In the case of BNS mergers, and for typical parameters ($\beta_{a}\sim 0.2$
and ${\tilde{L}}_{c}\sim 0.1$–$0.4$), these approximations would result in a
factor of $\sim 0.5$ being absorbed in $N_{s}$ [values of $N_{s}$ are given in
the caption of Table 1; for details refer to Appendix C.2 and equation (59)].
In the case of collapsars ($\beta_{a}=0$) the above expression is even simpler
[see equation (45) in Section 2.5.5], and this approximation results in a
factor of $\sim 0.7$ being absorbed in $N_{s}$ [for details refer to Appendix
C.2 and equation (58)].

Harrison et al. (2018) showed that $N_{s}$ depends on the actual value of
$\tilde{L}$ (i.e., $\tilde{L}_{c}$), but overall $N_{s}\sim 0.3-0.4$ for the
case of a non-relativistic collapsar jet. As a remark, since $N_{s}$ here is
used to absorb the above two approximations, its value differs depending on
the type of the jet (BNS merger case or collapsar case) and on the type of the
solution (semi-analytic or analytic; see Sections 2.4 and 2.5). Even for the
case of a collapsar jet, the values of $N_{s}$ here do differ slightly from
those in Harrison et al. (2018) [see the caption of Table 1 for the values of
$N_{s}$]. This is because additional difference in $N_{s}$ emerges as a result
of the difference in the modeling [e.g., difference in the modeling of the
cocoon’s lateral width, volume, and in the value of $E_{i}/E_{c}$ (or $\eta$)
compared to Harrison et al. 2018].

####  2.5.2 The approximated cocoon’s lateral width $r_{c}$

Here, the expression of $r_{c}(t)$ is simplified based on approximation that
the term $\sqrt{\frac{P_{c}}{\rho_{a}(r_{h}/2,t)c^{2}}}$ in the expression of
$\beta_{\perp}$ [in equation (17)] is considered as roughly constant over
time. This approximation is justified later by comparison with numerical
simulations. This allows us to write equation (16) as:

$\frac{dr_{c}(t)}{dt}+\left[-\frac{v_{m}}{r_{m}(t)}\right]r_{c}(t)=\sqrt{\frac{P_{c}}{\rho_{a}(r_{h}/2,t)}}.$
(24)

This is integrated, with $r_{c}(t=t_{0})$ being very negligible, as

$r_{c}(t)\approx\chi(t)\sqrt{\frac{P_{c}}{\rho_{a}(r_{h}/2,t)}}(t-t_{0}),$
(25)

where $\chi(t)$ is given by:

$\chi(t)=\frac{r_{m}(t)}{r_{m}(t)-r_{m,0}}\ln{\left[\frac{r_{m}(t)}{r_{m,0}}\right]}\approx\frac{t-t_{m}}{t-t_{0}}\ln{\left[\frac{t-t_{m}}{t_{0}-t_{m}}\right]},$
(26)

with $t_{m}$ being the time of the merger in the case of BNS mergers. The
value of $\chi(t)$ in BNS mergers depends on the time since the merger, the
time since the jet launch, and the time delay between the merger and the jet
launch. Typically $\chi(t)$ is found to take values as follows (also, see
Table 1 for the average values):

$\chi(t)\begin{cases}=1&\text{if $\beta_{a}=0$ (Collapsar case)},\\\ \sim
1-2&\text{if $\beta_{a}\sim 0.2-0.3$ (BNS merger case)}.\end{cases}$ (27)

Since $\chi(t)\propto\ln{t}$, its evolution over time is very limited.
Therefore, in order to further simplify the expression of $r_{c}(t)$, we
consider the time-averaged value of $\chi(t)$:

$\langle{\chi}\rangle=\frac{1}{t_{b}-t_{0}}\int_{t_{0}}^{t_{b}}\chi(t)dt,$
(28)

so that $r_{c}(t)$ is simplified to:

$r_{c}(t)\approx\langle{\chi}\rangle\sqrt{\frac{P_{c}}{\rho_{a}(r_{h}/2,t)}}(t-t_{0}).$
(29)

See Table 1 for the typical values of $\langle{\chi}\rangle$.

When deriving the breakout time $t_{b}$, $\langle{\chi}\rangle$ and $t_{b}$
depend on each other [see equation (44) for the expression of $t_{b}$, with
$A_{1}\propto 1/\sqrt{\langle{\chi}\rangle}$ as in equation (39)]. However,
this dependency is very weak, and small variation in the value of $t_{b}$
hardly affects the value $\langle{\chi}\rangle$. Therefore, both $t_{b}$ and
$\langle{\chi}\rangle$ can be determined iteratively 444Initially, a typical
value is assumed for $t_{b}$ and $\langle{\chi}\rangle$ based on a guess on
$t_{b}$ which can be guided by numerical simulations. Then $t_{b}$ is found
using equation (44) and a new value of $\langle{\chi}\rangle$ is found by
inserting $t_{b}$ in equation (28). This new value of $\langle{\chi}\rangle$
results in a slightly different $t_{b}$, which is used again (to find a more
accurate $\langle{\chi}\rangle$). This process is repeated $\sim 2-3$ times
until the values of $t_{b}$ and $\langle{\chi}\rangle$ converge..

####  2.5.3 The system of equations and the analytic solution

The system of equations (16), (17), (18), and (19) can be simplified to the
following:

$\displaystyle r_{c}(t)\approx$
$\displaystyle\langle{\chi}\rangle\sqrt{\frac{P_{c}}{\rho_{a}(r_{h}/2,t)}}(t-t_{0}),$
(30) $\displaystyle\beta_{\perp}=$
$\displaystyle\sqrt{\frac{P_{c}}{{\rho}_{a}(r_{h}/2,t)c^{2}}}+\left[\frac{r_{c}(t)}{r_{m}(t)}\right]\frac{v_{m}}{c},$
(31) $\displaystyle P_{{c}}=$
$\displaystyle\>\>\>\>\>\>\frac{E_{i}}{3\>V_{{c}}}\>\>\>\>\>\>\>=\langle{\eta^{\prime}}\rangle\frac{L_{j}\>(t-t_{0})}{2\pi
r_{c}^{2}(t)r_{{h}}(t)},$ (32) $\displaystyle\Sigma_{j}(t)=$ $\displaystyle\pi
r_{h}^{2}(t)\theta_{j}^{2}(t)=\frac{L_{j}\theta_{0}^{2}}{4cP_{c}}.$ (33)

From equation (7), ${\rho}_{a}(r_{h}/2,t)$ can be found. Then, replacing
$r_{c}(t)$ in equation (32), and using $r_{h}(t)\gg r_{0}$ and $r_{m,0}\gg
r_{0}$, $P_{c}$ can be written as:

$P_{c}=\sqrt{\frac{\langle{\eta^{\prime}}\rangle}{\langle{\chi}\rangle^{2}}\frac{(3-n)L_{j}M_{a}}{2^{3-n}\>\pi^{2}(t-t_{0})}\frac{r_{m}^{n-3}(t)}{r_{h}^{n+1}(t)}}.$
(34)

Finally, substituting equation (34) in equation (33) gives the expression of
the opening angle of the collimated jet as:

$\frac{\theta_{j}(t)}{\theta_{0}}=\left[\frac{\langle{\chi}\rangle^{2}}{\langle{\eta^{\prime}}\rangle}\frac{L_{j}}{M_{a}c^{2}}\frac{1}{(3-n)2^{n+1}}\right]^{\frac{1}{4}}\>\left[\frac{r_{h}(t)}{r_{m}(t)}\right]^{\frac{n-3}{4}}\>[t-t_{0}]^{\frac{1}{4}}.$
(35)

Notice the weak dependence of the jet opening angle on time, which has already
been pointed out in Hamidani et al. (2020).

The opening angle of the jet depends on the two parameters
$\langle{\chi}\rangle$ and $\langle{\eta^{\prime}}\rangle$. In the BNS merger
case, the ratio $\langle{\chi}\rangle^{2}/\langle{\eta^{\prime}}\rangle$ can
take values up to $\sim 10$; hence, these two parameters are important and
should not be overlooked.

The expression of $P_{c}$ and $\theta_{j}(t)$ [equations (34) and (35)] is
valid for both a BNS merger jet case and a collapsar jet case [where,
$r_{m}(t)\equiv r_{m}$ and $\langle{\chi}\rangle=1$].

####  2.5.4 Analytic solution for the BNS merger case

The jet head velocity [equation (23)] with equations (3), (6), (7), (10), and
(35), with further simplifications, gives the following differential equation:

$\displaystyle\frac{dr_{h}(t)}{dt}+\left[-\frac{v_{m}}{r_{m}(t)}\right]r_{h}(t)={A(t)}\>{r_{m}(t)}^{\frac{3-n}{2}}{r_{h}(t)}^{\frac{n-2}{2}},$
(36)

where:

$A(t)=N_{s}\sqrt{\left(\frac{r_{m,0}^{3-n}-r_{0}^{3-n}}{(3-n)\>r_{m,0}^{3-n}}\right)\left(\frac{4\>L_{j}}{\theta_{0}^{2}M_{a}\>c}\right)}\times\left[\frac{\theta_{0}}{\theta_{j}(t)}\right].$
(37)

The jet is initially uncollimated until the jet head’s radius, $r_{h}(t)$,
crosses the radius $\hat{z}/2$ [see equation (13)]. Since this initial phase
is very short (relative to the jet propagation timescale until the breakout;
see Figure 2), we consider as if the jet is in the collimated mode from the
start ($t=t_{0}$). Therefore, the jet opening angle can be found using
equation (35), and after inserting it in the expression of $A(t)$ [equation
(37)], the equation of motion [equation (36)] can be found as:

$\left[\frac{r_{h}(t)}{r_{m}(t)}\right]^{\frac{5-n}{4}}=A_{1}\frac{5-n}{4}\int{r_{m}^{-3/4}(t)[1-r_{m,0}/r_{m}(t)]^{-1/4}dt},$
(38)

where $A_{1}$, a constant, can be found as follows:

$\displaystyle A_{1}=$ $\displaystyle
N_{s}\left[\frac{\langle{\eta^{\prime}}\rangle}{\langle{\chi}\rangle^{2}}\right]^{\frac{1}{4}}\left[\left(\frac{r_{m,0}^{3-n}-r_{0}^{3-n}}{r_{m,0}^{3-n}}\right)^{2}\frac{2^{n+5}}{3-n}\frac{\>L_{j}\>v_{m}}{\theta_{0}^{4}M_{a}}\right]^{\frac{1}{4}},$
(39) $\displaystyle\approx$ $\displaystyle
N_{s}\left[\frac{\langle{\eta^{\prime}}\rangle}{\langle{\chi}\rangle^{2}}\right]^{\frac{1}{4}}\left[\frac{2^{n+5}}{3-n}\frac{\>L_{j}\>v_{m}}{\theta_{0}^{4}M_{a}}\right]^{\frac{1}{4}}.$
(40)

In the case where the delay between the merger time and the jet launch time,
$t_{0}-t_{m}$, is significantly smaller in comparison to the breakout time
$t_{b}-t_{m}$: $t_{0}-t_{m}\ll t_{b}-t_{m}$, we have $r_{m}(t)\gg r_{m,0}$,
hence the following approximation can be made:555The approximation
$[1-r_{m,0}/r_{m}(t)]\approx 1$ is not good in the early phase of jet
propagation where $r_{m,0}\simeq r_{m}(t)$. Still, since this early phase’s
timescale is very short (relative to the whole jet propagation timescale; see
Figure 2) this approximation is reasonable as long as $t_{b}-t_{m}\gg
t_{0}-t_{m}$.

$\int{r_{m}^{-3/4}(t)[1-r_{m,0}/r_{m}(t)]^{-1/4}dt}\simeq\int{r_{m}^{-3/4}(t)dt}.$
(41)

With the boundary conditions $r_{m}(t_{0})=r_{m,0}$ and $r_{h}(t_{0})=r_{0}$
at $t=t_{0}$, the integration gives:

$r_{h}(t)=\left\\{\frac{(5-n)A_{1}}{v_{m}}(r_{m}^{\frac{1}{4}}(t)-r_{m,0}^{\frac{1}{4}})+\left[\frac{r_{0}}{r_{m,0}}\right]^{\frac{5-n}{4}}\right\\}^{\frac{4}{5-n}}r_{m}(t).$
(42)

The jet head velocity can be deduced from equation (42) as:

$v_{h}(t)=v_{m}\left[\frac{r_{h}(t)}{r_{m}(t)}\right]+A(t)\left[\frac{r_{h}(t)}{r_{m}(t)}\right]^{\frac{n-2}{2}}[r_{m}(t)^{\frac{1}{4}}(r_{m}(t)-r_{m,0})^{\frac{1}{4}}].$
(43)

Finally, the breakout time can be derived by taking
$r_{h}(t_{b})/r_{m}(t_{b})=1$ in equation (42):

$t_{b}-t_{0}=\left\\{\frac{v_{m}^{\frac{3}{4}}}{(5-n)A_{1}}\left[1-\left[\frac{r_{0}}{r_{m,0}}\right]^{\frac{5-n}{4}}\right]+\left(\frac{r_{m,0}}{v_{m}}\right)^{\frac{1}{4}}\right\\}^{4}-\frac{r_{m,0}}{v_{m}}.$
(44)

####  2.5.5 Analytic solution for the collapsar case

This case is a special case from the previous one (in Section 2.5.4) where
$v_{m}=0$ [i.e., the ambient medium is static: $\beta_{a}=0$, $\chi(t)=1$ and
$r_{m}(t)\equiv r_{m}$]. Therefore, the equation of motion [equation (11)],
after being approximated to equation (23) (see Section 2.5.1), can be further
simplified to the following:

$\displaystyle\beta_{h}\approx\tilde{L}_{c}^{\frac{1}{2}}.$ (45)

Hence, the equation of motion for the jet head can be found as:

$\displaystyle\frac{dr_{h}(t)}{dt}={A(t)}\>{r_{m}}^{\frac{3-n}{2}}{r_{h}(t)}^{\frac{n-2}{2}},$
(46)

where $A(t)$ here is:

$A(t)=N_{s}\sqrt{\left(\frac{r_{m}^{3-n}-r_{0}^{3-n}}{(3-n)\>r_{m}^{3-n}}\right)\left(\frac{4\>L_{j}}{\theta_{0}^{2}M_{a}\>c}\right)}\times\left[\frac{\theta_{0}}{\theta_{j}(t)}\right],$
(47)

which is the same expression as in equation (37) [where here $r_{m,0}\equiv
r_{m}$]. As in Section 2.5.4, the initial uncollimated phase is overlooked for
simplicity. Then the expression of $\theta_{j}(t)/\theta_{0}$ can be found
using equation (35) [with $r_{m}(t)=r_{m}$ in the collapsar case]. Inserting
$\theta_{j}(t)/\theta_{0}$ in the above expression of $A(t)$ [equation (47)],
integrating equation (46), and using the boundary condition
$r_{h}(t_{0})=r_{0}$, gives the following expression for the jet head radius:

$r_{h}(t)=\left\\{\frac{(5-n)A_{1}^{\prime}}{3}(t-t_{0})^{\frac{3}{4}}+\left[\frac{r_{0}}{r_{m}}\right]^{\frac{5-n}{4}}\right\\}^{\frac{4}{5-n}}r_{m},$
(48)

and the jet head velocity:

$v_{h}(t)=A_{1}^{\prime}r_{m}\left\\{\frac{(5-n)A_{1}^{\prime}}{3}(t-t_{0})^{\frac{3}{4}}+\left[\frac{r_{0}}{r_{m}}\right]^{\frac{5-n}{4}}\right\\}^{\frac{n-1}{5-n}}(t-t_{0})^{-\frac{1}{4}},$
(49)

where $A_{1}^{\prime}$ is a constant that can be written as:

$\displaystyle A_{1}^{\prime}=$ $\displaystyle
N_{s}\left[\frac{\langle{\eta^{\prime}}\rangle}{\langle{\chi}\rangle^{2}}\right]^{\frac{1}{4}}\left[\left(\frac{r_{m}^{3-n}-r_{0}^{3-n}}{r_{m}^{4-n}}\right)^{2}\frac{2^{n+5}}{3-n}\frac{L_{j}}{\theta_{0}^{4}M_{a}}\right]^{\frac{1}{4}},$
(50) $\displaystyle\approx$ $\displaystyle
N_{s}\left[\frac{\langle{\eta^{\prime}}\rangle}{\langle{\chi}\rangle^{2}}\right]^{\frac{1}{4}}\left[\frac{1}{r_{m}^{2}}\frac{2^{n+5}}{3-n}\frac{L_{j}}{\theta_{0}^{4}M_{a}}\right]^{\frac{1}{4}}.$
(51)

The breakout time can be found for $r_{h}(t_{b})=r_{m}$ as:

$t_{b}-t_{0}=\left\\{\frac{3}{(5-n)A_{1}^{\prime}}\left[1-\left(\frac{r_{0}}{r_{m}}\right)^{\frac{5-n}{4}}\right]\right\\}^{\frac{4}{3}}.$
(52)

##  3 Comparison with numerical simulations

###  3.1 Numerical simulations

In addition to the analytic (and semi-analytic) modeling presented above, we
carried out a series of 2D relativistic hydrodynamical simulations. In total,
the series includes a total of over a hundred models covering a wide parameter
space (see Table 1 in Hamidani et al. 2020). The essential aim of carrying out
numerical simulations here is to test the semi-analytic (Section 2.4) and the
analytic (Section 2.5) solutions, and calibrate them if necessary. These tests
and calibrations are presented for both, the case of BNS merger jet, and the
case of collapsar jet.

We pick up four models, as a subsample, out of our sample of numerical
simulations. As presented in Table 1, two out of four are BNS merger models,
with different initial opening angles (T03-H and T13-H); and the other two are
collapsar models, with different initial opening angles as well (A and B). The
parameters of the stellar envelope in collapsar simulations (models A and B)
follows the widely used model 16TI in Woosley & Heger (2006). However, for
simplicity, the density profile is approximated to a power-law function with
an index $n=2$ [see (iv) in Section 2.2]. This allows the analytic results to
be tested fairly with simulations.

The injection radius, $r_{in}$, is set at $1.2\times 10^{8}$ cm for the BNS
merger case, and $10^{9}$ cm for the collapsar case (see Table 1). This might
seem quite large; ideally the injection radius should be of the order of
$10^{7}$ cm. However, since the density profile in the inner region can be
approximated to a power-law function with an index $n<3$ ($n=2$ for the
dynamical ejecta of BNS mergers, see Figure 8 in Hamidani et al. 2020;
$n\approx 1.5$ for the collapsar case, see Figure 2 in Margutti et al. 2013),
the mass contained in the inner region ($<10^{8}$ cm in BNS mergers; $<10^{9}$
cm in collapsars) is negligible, in comparison to the total ambient medium
mass, as long as $r_{in}\ll r_{m,0}$ (or $r_{in}\ll r_{m}$), which is the case
in our simulations [see equation (5)]. Hence, this inner region is expected to
have a very limited effect on the overall jet dynamics. For an estimation of
the effect of these values on jet dynamics, please refer to the analytic model
[in particular refer to equations (44), and (52), for the effect of the value
of the inner boundary, $r_{0}$, on the jet breakout time]. Note that the
motivation for taking such large values for $r_{in}$ is because smaller
injection radii make numerical simulation extremely expensive in terms of
computation time.

Further details about the numerical code are presented in Hamidani et al.
(2017). For more information about the setup of the numerical simulations,
refer to Section 3.1 in Hamidani et al. (2020).

###  3.2 Measurement of the internal energy in the cocoon and $\eta^{\prime}$

Figure 1 shows the time evolution of the fraction of internal energy in the
cocoon $E_{i}/E_{c}$, and the two parameters $\eta$ and $\eta^{\prime}$
[previously defined in equations (20) and (21); see Section 2.3.2], from the
jet launch time ($t=t_{0}$) to the jet breakout time ($t=t_{b}$), as measured
from numerical simulation 666The data was deduced by measuring $E_{i}$ (or the
cocoon pressure $P_{c}$ as presented in Figure 4, combined with the total
volume of the cocoon $V_{c}$), the total energy contained in the cocoon
$E_{c}$, and the jet head radius [$r_{h}(t)$ or $\langle{\beta_{h}}\rangle$],
all from numerical simulations.. For comparison, both the case of collapsar
jet and the case of BNS merger jet are shown. In the collapsar jet case (model
A and B), the fraction of internal energy in the cocoon, $E_{i}/E_{c}$, is
high ($\sim 0.7-0.8$). The values of $\eta$ and $\eta^{\prime}$ (in the range
$\sim 0.5-1$) are also high. On the other hand, in the BNS merger jet case
(model T03-H and T13-H), where the medium is expanding, the values of
$E_{i}/E_{c}$, $\eta$, and $\eta^{\prime}$ are significantly lower ($<0.5$).

This contrast is related to the adiabatic expansion of the cocoon, which is
very effective in the case where the medium is expanding; as the medium’s
expansion velocity (i.e., the dynamical ejecta’s radial velocity) is
comparable to the jet head velocity, up to the breakout time, this expansion
enhances the volume of the (over-pressurized) cocoon significantly, while
depleting its internal energy. Through this process, in the BNS merger case,
the inner region of the cocoon, where initially the density is high and the
expansion velocity is small, is propelled further outward within the cocoon,
up to velocities of the order of the homologous expansion of the medium,
gaining kinetic energy at the expense of internal energy. Note that in the
collapsar case, the same process happens but to a lesser extent; the inner
cocoon is propelled outward, but as a result of the initially static medium,
the gained velocity (and fraction of kinetic energy) is much less significant.

Another reason is the high jet head velocity ($\langle{\beta_{h}}\rangle$) in
the BNS merger jet case (roughly $\sim 2v_{m}/c$; see Ioka & Nakamura 2018),
implying that the fraction of the injected energy that reaches the cocoon
[$\propto(1-\langle{\beta_{h}}\rangle)$] is less significant [relative to the
case of a collapsar jet; see equation (55)]. For more details refer to
Appendix A.

To our best knowledge, this is the first time that the fraction of internal
energy of the cocoon, and the parameter $\eta$ (and $\eta^{\prime}$), has been
measured (from simulations), and such a significant difference, between the
case of collapsar jet and the case of BNS merger jet, has been found. The
parameter $\eta$ has been discussed in the literature in the case of collapsar
jets, and it has been suggested to take a value of $\sim 1$ (e.g., in Bromberg
et al. 2011). Our results show that this assumption is quite reasonable. On
the other hand, several recent works naively assumed the same value, $\eta\sim
1$, for the case of BNS merger jet (e.g., Matsumoto & Kimura 2018; Gill et al.
2019; Salafia et al. 2020). Here, we show that such assumption is not
reasonable by a factor of $\sim 2$; $\eta$ should rather be smaller in the
case of BNS merger jet, as it can be seen in Figure 1 (unless $t_{b}-t_{m}\sim
t_{0}-t_{m}$).

In summary, $\eta^{\prime}$ is found to take values as follows:

$\eta^{\prime}\begin{cases}\sim 0.5-1&\text{if $\beta_{a}=0$ (Collapsar jet
case)},\\\ \sim 0.1-0.5&\text{if $\beta_{a}\gg 0$ (BNS merger jet
case)}.\end{cases}$ (53)

As a remark, for typical cases, naively assuming $\eta\sim 1$ for the case of
BNS merger tends to incorrectly give a factor of $\sim\sqrt{2}$ times more
collimation [i.e., $\sim\sqrt{2}$ times more jet head velocity, and hence much
shorter breakout times; see equations (35); (39); and (44)].

Figure 1: The fraction of the internal energy to the total energy of the
cocoon (top); and the parameters $\eta$ (middle; the fraction of the cocoon’s
internal energy to the energy injected into the cocoon) and $\eta^{\prime}$
(bottom; the fraction of the cocoon’s internal energy to the injected jet
energy), as measured in our 2D simulations [see equations (20) and (21)]. The
red and dark red lines are for collapsar jet models (models A and B in Table
1). The blue and dark blue lines are for BNS merger models (models T03-H and
T13-H in Table 1).

Figure 2: Results for the case of BNS mergers showing the jet’s and the
cocoon’s evolution over time, as measured in numerical simulations (black
dotted line with filled squares), and as inferred with the analytic (solid
blue line) and the semi-analytic (solid red line) solutions [see Sections 2.4
and 2.5.4]. The black dotted line in the top two panels shows the outer radius
of the expanding ejecta. From top to bottom, the jet head radius, the cocoon’s
lateral width, and the jet opening angle relative to the initial opening
angle, are shown respectively. From left to right, results for the narrow jet
case (model T03-H) and the wide jet case (model T13-H) are shown respectively.
We take $N_{s}=0.46$ in the analytic solution, and $N_{s}=0.75$ in the semi-
analytic solution (see Section 2.5.1 and Table 1).

Figure 3: Same as Figure 2 for the collapsar case, showing the results from numerical simulations (black dotted line with filled squares), and results of the analytic (solid blue line) and the semi-analytic (solid red line) solutions [see Sections 2.4 and 2.5.5]. The horizontal black dotted line in the top two panels shows the radius of the stellar envelope. From left to right, the results for the narrow jet case (model A) and the results for the wide jet case (model B) are shown, respectively. We take $N_{s}=0.38$ in the analytic solution, and $N_{s}=0.53$ in the semi-analytic solution (see Section 2.5.1 and Table 1).  Table 1: A subsample showing the simulated models and the corresponding parameters. From the left: The model name; the type of jet (BNS merger or collapsar); the ambient medium’s mass; the jet initial opening angle; the engine’s isotropic equivalent luminosity [$L_{iso,0}=\frac{2L_{j}}{1-\cos\theta_{0}}\simeq\frac{4L_{j}}{\theta_{0}^{2}}$] where $L_{j}$ is the jet true luminosity (one sided); the inner radius at which the jet is injected in simulations; the ambient medium’s outer radius at the start of the simulation; the maximum expansion velocity of the ambient medium; the time-averaged value of $\eta^{\prime}$ [see equation (22)] estimated from simulations (see Figure 1); the time-averaged value of $\chi(t)$ used in the analytic solution [using equation (28); see Section 2.5.2]; the breakout time measured in numerical simulations; the inferred breakout time by the analytic solution [using equation (44) for the BNS merger jet case, and equation (52) for the collapsar jet case], and by the semi-analytic solution (see Section 2.4); The values of the calibration coefficient $N_{s}$ are: $0.46$ and $0.75$ in the BNS merger case (for the analytic and the semi-analytic solution, respectively), $0.38$ and $0.53$ in the collapsar case (for the analytic and the semi-analytic solution, respectively). The density profile of the ambient medium in all models is approximated to power-law with the index $n=2$ [see (iv) in Section 2.2]. |  | $M_{a}$ | $\theta_{0}$ | $L_{iso,0}$ | $r_{in}$ | $r_{m}(t_{0})$ | $v_{m}$ | $\langle{\eta^{\prime}}\rangle$ | $\langle{\chi}\rangle$ |  $t_{b}-t_{0}$ [s] |  $t_{b}-t_{0}$ [s] |  $t_{b}-t_{0}$ [s]
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  |  |  |  |  |  | (Semi-
Model | Type | [$M_{\odot}$] | [deg] | [erg s-1] | [cm] | [cm] | [c] |  |  | (Simulation) | (Analytic) | analytic)
T03-H | BNS | $0.002$ | 6.8 | $5\times 10^{50}$ | $1.2\times 10^{8}$ | $1.67\times 10^{9}$ | $\frac{\sqrt{3}}{5}$ | $1/4$ | $1.25$ | 0.221 | 0.203 | 0.231
T13-H | BNS | $0.002$ | 18.0 | $5\times 10^{50}$ | $1.2\times 10^{8}$ | $1.67\times 10^{9}$ | $\frac{\sqrt{3}}{5}$ | $1/4$ | $1.48$ | 0.429 | 0.456 | 0.408
A | Collapsar | $13.950$ | 9.2 | $7.83\times 10^{52}$ | $10^{9}$ | $4\times 10^{10}$ | 0 | $1/2$ | $1.00$ | 3.804 | 3.282 | 3.947
B | Collapsar | $13.950$ | 22.9 | $1.27\times 10^{52}$ | $10^{9}$ | $4\times 10^{10}$ | 0 | $1/2$ | $1.00$ | 9.930 | 11.137 | 9.681

###  3.3 Time evolution of the jet propagation

####  3.3.1 BNS merger’s case

In Figure 2 we show results for the two models of jet propagation in the BNS
merger ejecta [T03-H and T13-H with $\theta_{0}=6.8^{\circ}$ (left) and
$18.0^{\circ}$ (right), respectively; see Table 1]. Three different quantities
are shown (from top to bottom): the jet head radius $r_{h}(t)$, the cocoon’s
lateral width $r_{c}(t)$, and the opening angle of the jet head
$\theta_{j}(t)$.

The calibration coefficient $N_{s}$ has been used to calibrate the analytic
solution (with $N_{s}=0.46$) and the semi-analytic solution (with
$N_{s}=0.75$); the value of $N_{s}$ is set so that the breakout time in the
analytic (or semi-analytic) solution matches the breakout time in numerical
simulations [refer to equation (10) and the explanation that follows]. Also,
as noted in Section 2.5.1, the different values of $N_{s}$ are due to the
additional approximations in the analytic solution. It should be noted that
the value of $N_{s}$ for the analytic model, here (0.46), is slightly
different from the one in the analytic model presented in Hamidani et al.
(2020) (0.40). This difference is due to the main difference between the two
models; in Hamidani et al. (2020) the jet opening angle is fixed using the
parameter $f_{j}$ (measured from numerical simulations; see Figure 3 in
Hamidani et al. 2020), while here, the jet opening angle is determined self-
consistently (by calculating the jet collimation by the cocoon), and should be
more robust (see the bottom two panels in Figure 2). For more details on
$N_{s}$ refer to Appendix C.

The time evolution of the analytic and the semi-analytic jet head radius,
$r_{h}(t)$, in Figure 2 shows a very good agreement with simulations (within
$\sim 10\%$). Analytic and semi-analytic results hold fairly well for both
models (T03-H and T13-H) showing that the jet-cocoon model here works well
regardless of the initial jet opening angle $\theta_{0}$.

The time evolution of the analytic and semi-analytic cocoon’s lateral width
$r_{c}(t)$ is also consistent with simulations, especially for the case of
small $\theta_{0}$ (within $\sim 10\%$; see T03-H in Figure 2). For the case
of large $\theta_{0}$ (T13-H), the agreement with simulations is less
significant, but roughly within $\sim 30\%$, where $r_{c}(t)$ is slightly
underestimated in the analytic and semi-analytical model (relative to
numerical simulations).

In simulations, the jet head’s opening angle has been estimated by taking the
average opening angle from $r=\frac{1}{2}r_{h}(t)$ to $r=r_{h}(t)$ [see
equation (34) in Hamidani et al. 2020]. This average opening angle is compared
with the analytic and the semi-analytic jet opening angles in Figure 2. Except
the early time evolution of the jet opening angle, during which the jet-cocoon
is highly inhomogeneous in simulations [in particular, in terms of entropy and
Lorentz factor which are used to discriminate the jet from the cocoon: for
more details see Section 3.1.2 in Hamidani et al. (2020)], analytic and semi-
analytic jet opening angles are consistent with the average opening angle in
simulations within $\sim 30\%$.

####  3.3.2 Collapsar’s case

Figure 3 shows a comparison of the analytic and the semi-analytic results with
simulations for three quantities $r_{h}(t)$, $r_{c}(t)$, and $\theta_{j}(t)$,
in the same manner as in Figure 2 (and Section 3.3.1), but for the collapsar
jet case. We present two models with different initial opening angles [A
(left) and B (right), with $\theta_{0}=9.2^{\circ}$ and $22.9^{\circ}$,
respectively; see Table 1]. The calibration coefficient is found as
$N_{s}=0.38$ for the analytic solution, and $N_{s}=0.53$ for the semi-analytic
solution (see Section 2.5.1 for more details about the origin of this
difference). Although slightly larger, these values of $N_{s}$ are fairly
consistent with those found by Harrison et al. (2018), despite several
differences in the jet-cocoon modeling (such as for the cocoon’s lateral
width, volume, $\langle{\eta^{\prime}}\rangle$, etc.).

The analytic and the semi-analytic solutions for the time evolution of the jet
head radius, $r_{h}(t)$, show a clear agreement with simulations (within $\sim
10$–$20\%$) for both models (A and B). The same can be said about the time
evolution of the cocoon’s lateral width $r_{c}(t)$ (within $\sim 10-20\%$).

The time evolution of the average opening angle of the jet head
$\theta_{j}(t)$ in simulation can be divided into two phases; with the first
phase showing relatively large opening angles and unstable behavior, and the
second phase showing collimated opening angles, and relatively a stable
behavior. In the first phase, the effects of the initial conditions are still
present. However, since this phase is relatively short, its contribution to
the jet structure and propagation, up to the breakout, is limited. In the
second phase, which represents most of the jet propagation time, the jet
head’s opening angle in the analytic and semi-analytic solutions is, overall,
consistent with the average opening angle in simulations (well within $\sim
50\%$).

###  3.4 The cocoon pressure $P_{c}$

Figure 4 shows the average cocoon pressure in simulations (measured throughout
the cocoon’s grid in numerical simulations and volume-averaged) compared with
the cocoon pressure as inferred from our analytic and semi-analytic solutions
[see Section 2.4; and equation (34) in Section 2.5]. For both cases (BNS
merger jet and collapsar jet), and for both the analytic and the semi-analytic
solutions, the time evolution of the cocoon pressure up to the breakout is
well consistent with numerical simulations. This agreement indicates that the
modeling presented here is a good representation of the cocoon and its
interaction with the jet and the ambient medium.

Note that in the analytic solution, as the early uncollimated jet phase is not
taken into account, the cocoon pressure diverges at $t\sim t_{0}$ due to the
approximation in equation (41). The semi-analytic solution shows no such
anomaly.

Figure 4: The pressure of the cocoon $P_{c}$ from the jet launch time $t_{0}$
to the breakout time. The top two panels show two BNS merger models (T03-H and
T13-H, from left to right), and the bottom two panels show two collapsar
models (A and B, from left to right; see Table 1). The black dotted line with
filled squares shows the average pressure in the cocoon as measured in
numerical simulation. The blue line shows the cocoon pressure according to the
analytic model [equation (34)]. The red solid line shows the cocoon pressure
according to the semi-analytic model (calculated numerically; see Section
2.4).

Figure 5: Density maps showing the jet-cocoon system inside the ambient medium
just before the jet breakout. Four models are shown (see Table 1), where the
top two models are for jet propagation in BNS merger ejecta (T03-H and T13-H),
and the bottom two are for jet propagation in a stellar envelope (collapsar
jets; model A and B). The black filled square shows the inferred jet head
radius $r_{h}$ by our semi-analytic model, and the two black filled circles
show the inferred cocoon’s lateral width $r_{c}$. The ellipsoidal shape (solid
black line) shows the jet-cocoon’s shape as predicted by our modeling (the
semi-analytic solution).

###  3.5 Morphology of the jet-cocoon system

Figure 5 presents snapshots from our numerical simulations showing the density
map of the jet-cocoon system just before its breakout out of the ambient
medium. The four models in Table 1 are shown for BNS mergers (top) and
collapsars (bottom).

The jet head radius and the cocoon’s lateral width, as inferred from our semi-
analytic solution, are also shown for comparison with simulations (with a
black filled square, and two black filled circles, respectively). Also, the
inferred jet-cocoon morphology using the approximation of an ellipsoidal shape
is shown (with a solid black line) where the semi-major axis is
$\frac{1}{2}r_{h}$ and the semi-minor axis is $r_{c}$.

In Figure 5, we see a clear similarity between the morphology of the whole
jet-cocoon system as inferred from our modeling and numerical simulations.
With $r_{h}$ being well consistent with simulations, and the error on the
cocoon’s lateral width of the order of $\sim 20\%$ [at
$r\sim\frac{1}{2}r_{h}(t)$], it can be claimed that our modeling can robustly
give the cocoon volume (within $\sim 50\%$). Together, with the modeling of
the cocoon pressure (see Section 3.4), we can conclude that all aspects of the
jet-cocoon system are fairly well reproduced with our modeling.

##  4 Conclusion

In this paper we present a new jet-cocoon model. The model is based on
previous works of collapsar jet-cocoon modeling, in particular, models in
Matzner (2003); Bromberg et al. (2011); Mizuta & Ioka (2013); Harrison et al.
(2018). From the analysis of jet propagation in numerical simulations over a
wide parameter range, the model has been generalized to enable proper
treatment of jet propagation in the case of BNS merger where the ambient
medium is expanding. For each jet case, equations have been solved through
numerical integration (semi-analytic solution), or analytically (analytic
solution) after adding some approximations (see Section 2.5.1 and Section
2.5.2). Table 2 presents an overview of previous works on the modeling of GRB-
jet propagation, and a comparison with this work. In summary, our results can
be summarized as follows:

  1. 1.

Comparisons with numerical simulations show that our model’s results are in a
clear agreement with numerical simulations (overall, within $\sim 20\%$). The
time evolution of the following quantities has been shown to be well
consistent with measurements from numerical simulations: the jet head radius
$r_{h}(t)$, the cocoon radius $r_{c}(t)$, the jet head opening angle
$\theta_{j}(t)$, and the cocoon pressure $P_{c}$ (see Figures 2, 3 and 4). The
cocoon’s morphology, and volume, as inferred with our model are also
consistent with numerical simulations (see Figure 5). This is the first time
that results from the modeling of jet propagation in an expanding medium (as
in BNS mergers) has been compared with numerical simulations over such a large
set of parameters, and found consistent to such extent (see Table 1).

  2. 2.

The results of our jet-cocoon model are proven to be consistent with numerical
simulations regardless of the jet case (collapsar jet or BNS merger jet), and
regardless of the jet initial opening angle (see Figures 2, 3, 4, and 5).

  3. 3.

In addition to the semi-analytic solution, where equations are solved through
numerical integration, our model offers an analytic solution, where equations
are solved analytically after being approximated and simplified (see Section
2.5). Still, we showed that, within certain conditions (e.g., $t_{0}-t_{m}\ll
t_{b}-t_{m}$), the analytic solution’s results are, almost, as consistent with
numerical simulations as the semi-analytic model’s results, despite being much
simpler. This analytic modeling, with its simplified but fairly robust
equations [e.g., equation (44) for the breakout time], is very useful for
further investigations; for instance, on the cocoon’s cooling emission
(Hamidani et al. in prep).

  4. 4.

The composition of the cocoon energy has been measured for both jet cases
(collapsar and BNS merger), thanks to numerical simulations. Results show a
clear contrast between the two cases; the cocoon energy in the case of a
collapsar jet is overwhelmingly dominated by internal energy, while in the
case of a BNS merger jet it is, rather, overwhelmingly dominated by kinetic
energy. This is the first time that such difference has been revealed. As a
result of this difference in internal energy, we showed that the parameter
$\eta$ [see equation (18)] is much smaller ($\sim 2$ times) in the case of BNS
mergers than in the case of collapsars (see Figure 1). Such difference has not
been taken into account in previous works although it substantially affects
every aspect of the jet propagation. Also, such difference in internal energy
(in the BNS merger case) is very important when estimating the cooling
emission of the cocoon; hence the importance of this result (previously
mentioned in Kimura et al. 2019).

It should be noted that the analytic modeling presented here includes several
limitations. The most important limitation is that jets here are assumed as
unmagnetized. Other notable limitations are that effects from neutrinos,
r-process, viscous wind, general relativity, stellar rotation, stellar
magnetic field, etc. have been overlooked for the sake of simplicity. Future
works are likely to update our results.

Finally, it should also be noted that due to the limited computational
resources, the numerical simulations presented here use the approximation of
axial symmetric jets (2D) and jets are injected at relatively larger radii.
This may result in some numerical artifacts. Therefore, results such as; the
value of $N_{s}$; and the overall agreement of $\sim 20\%$ between analytic
results and numerical simulations; should not be taken at face value. We
expect these values to be updated in the future once more refined numerical
simulations are available.

Table 2: Comparison of jet-cocoon models in the literature.  | Context (medium): |  | Consistency |
---|---|---|---|---
| BNS mergers | Collapsars | Analytic | with |
Work | (expanding) | (static) | solution | simulations | Comment
Bromberg et al. (2011) | No | ✓ | ✓ | ✓ | Limited to the collapsar case
Mizuta & Ioka (2013) | No | ✓ | ✓ | ✓ | Limited to the collapsar case
Margalit et al. (2018) | No | ✓∗ | ✓ | ? |  ∗Jet propagation in SLSN ejecta.
Harrison et al. (2018) | No | ✓ | ✓ | ✓ | Limited to the collapsar case
Duffell et al. (2018) | ✓ | No | ✓ | ✓ | No treatment for jet collimation.
Matsumoto & Kimura (2018) | ✓ | No | ✓ | No | Overlooks $\eta$ and $\chi$.
Lazzati & Perna (2019) | ✓ | No | ✓ | ? | Describes the jet-wind interaction.
Gill et al. (2019) | ✓ | No | ✓ | ? | Overlooks $\eta$ and $\chi$.
Salafia et al. (2020) | ✓ | ✓ | No | ✓ | The effect of the expansion was not included
|  |  |  |  | [in $\beta_{\perp}$, second term in equation (17); in $\eta$; etc.].
Hamidani et al. (2020) | ✓ | ✓ | ✓ | ✓ | Simplified, after showing that $\theta_{j}(t)\sim\rm{Constant}$.
Lyutikov (2020) | ✓ | No | ✓ | No | No treatment for jet collimation.
Beniamini et al. (2020) | ✓ | No | ✓ | No | No treatment for jet collimation.
This work | ✓ | ✓ | ✓ | ✓ |

## Acknowledgements

We thank Amir Levinson, Atsushi Taruya, Bing Zhang, Christopher M. Irwin,
Hendrik van Eerten, Hirotaka Ito, Kazumi Kashiyama, Kazuya Takahashi, Kenta
Kiuchi, Kohta Murase, Koutarou Kyutoku, Masaomi Tanaka, Masaru Shibata, Ore
Gottlieb, Tomoki Wada, Toshikazu Shigeyama, Tsvi Piran, and Yudai Suwa, for
fruitful discussions and comments. We thank the participants and the
organizers of the workshops with the identification number YITP-T-19-04,
YITP-W-18-11 and YITP-T-18-06, for their generous support and helpful
comments.

Numerical computations were achieved thanks to the following: Cray XC50 of the
Center for Computational Astrophysics at the National Astronomical Observatory
of Japan, and Cray XC40 at the Yukawa Institute Computer Facility.

This work is partly supported by JSPS KAKENHI nos. 20H01901, 20H01904,
20H00158, 18H01213, 18H01215, 17H06357, 17H06362, 17H06131 (KI).

##  5 Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2017a) Abbott B. P., et al., 2017a, Physical Review Letters, 119, 161101
  * Abbott et al. (2017b) Abbott B. P., et al., 2017b, ApJ, 848, L13
  * Bauswein et al. (2013) Bauswein A., Goriely S., Janka H. T., 2013, ApJ, 773, 78
  * Begelman & Cioffi (1989) Begelman M. C., Cioffi D. F., 1989, ApJ, 345, L21
  * Beniamini et al. (2020) Beniamini P., Duran R. B., Petropoulou M., Giannios D., 2020, ApJ, 895, L33
  * Bromberg et al. (2011) Bromberg O., Nakar E., Piran T., Sari R., 2011, ApJ, 740, 100
  * Bromberg et al. (2018) Bromberg O., Tchekhovskoy A., Gottlieb O., Nakar E., Piran T., 2018, MNRAS, 475, 2971
  * Ciolfi & Vijay Kalinani (2020) Ciolfi R., Vijay Kalinani J., 2020, arXiv e-prints, p. arXiv:2004.11298
  * Duffell et al. (2015) Duffell P. C., Quataert E., MacFadyen A. I., 2015, ApJ, 813, 64
  * Duffell et al. (2018) Duffell P. C., Quataert E., Kasen D., Klion H., 2018, ApJ, 866, 3
  * Eichler et al. (1989) Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126
  * Fujibayashi et al. (2018) Fujibayashi S., Kiuchi K., Nishimura N., Sekiguchi Y., Shibata M., 2018, ApJ, 860, 64
  * Fujibayashi et al. (2020) Fujibayashi S., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M., 2020, arXiv e-prints, p. arXiv:2007.00474
  * Gill et al. (2019) Gill R., Nathanail A., Rezzolla L., 2019, ApJ, 876, 139
  * Goodman (1986) Goodman J., 1986, ApJ, 308, L47
  * Gottlieb et al. (2018a) Gottlieb O., Nakar E., Piran T., 2018a, MNRAS, 473, 576
  * Gottlieb et al. (2018b) Gottlieb O., Nakar E., Piran T., Hotokezaka K., 2018b, MNRAS, 479, 588
  * Gottlieb et al. (2020) Gottlieb O., Nakar E., Bromberg O., 2020, arXiv e-prints, p. arXiv:2006.02466
  * Hamidani et al. (2017) Hamidani H., Takahashi K., Umeda H., Okita S., 2017, MNRAS, 469, 2361
  * Hamidani et al. (2020) Hamidani H., Kiuchi K., Ioka K., 2020, MNRAS, 491, 3192
  * Harrison et al. (2018) Harrison R., Gottlieb O., Nakar E., 2018, MNRAS, 477, 2128
  * Hotokezaka et al. (2013) Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., Shibata M., Taniguchi K., 2013, Phys. Rev. D, 87, 024001
  * Ioka & Nakamura (2018) Ioka K., Nakamura T., 2018, Progress of Theoretical and Experimental Physics, 2018, 043E02
  * Ioka & Nakamura (2019) Ioka K., Nakamura T., 2019, MNRAS, 487, 4884
  * Kimura et al. (2019) Kimura S. S., Murase K., Ioka K., Kisaka S., Fang K., Mészáros P., 2019, ApJ, 887, L16
  * Lazzati & Perna (2019) Lazzati D., Perna R., 2019, ApJ, 881, 89
  * Lazzati et al. (2017a) Lazzati D., Deich A., Morsony B. J., Workman J. C., 2017a, MNRAS, 471, 1652
  * Lazzati et al. (2017b) Lazzati D., López-Cámara D., Cantiello M., Morsony B. J., Perna R., Workman J. C., 2017b, ApJ, 848, L6
  * Lazzati et al. (2018) Lazzati D., Perna R., Morsony B. J., Lopez-Camara D., Cantiello M., Ciolfi R., Giacomazzo B., Workman J. C., 2018, Phys. Rev. Lett., 120, 241103
  * Lazzati et al. (2020) Lazzati D., Ciolfi R., Perna R., 2020, arXiv e-prints, p. arXiv:2004.10210
  * Lyutikov (2020) Lyutikov M., 2020, MNRAS, 491, 483
  * MacFadyen & Woosley (1999) MacFadyen A. I., Woosley S. E., 1999, ApJ, 524, 262
  * Margalit et al. (2018) Margalit B., Metzger B. D., Thompson T. A., Nicholl M., Sukhbold T., 2018, MNRAS, 475, 2659
  * Margutti et al. (2013) Margutti R., et al., 2013, MNRAS, 428, 729
  * Martí et al. (1997) Martí J. M., Müller E., Font J. A., Ibáñez J. M. Z., Marquina A., 1997, ApJ, 479, 151
  * Matsumoto & Kimura (2018) Matsumoto T., Kimura S. S., 2018, ApJ, 866, L16
  * Matzner (2003) Matzner C. D., 2003, MNRAS, 345, 575
  * Mizuta & Ioka (2013) Mizuta A., Ioka K., 2013, ApJ, 777, 162
  * Mooley et al. (2018) Mooley K. P., et al., 2018, Nature, 561, 355
  * Murguia-Berthier et al. (2014) Murguia-Berthier A., Montes G., Ramirez-Ruiz E., De Colle F., Lee W. H., 2014, ApJ, 788, L8
  * Nagakura et al. (2014) Nagakura H., Hotokezaka K., Sekiguchi Y., Shibata M., Ioka K., 2014, ApJ, 784, L28
  * Nakar & Piran (2017) Nakar E., Piran T., 2017, ApJ, 834, 28
  * Nakar & Piran (2020) Nakar E., Piran T., 2020, arXiv e-prints, p. arXiv:2005.01754
  * Nakar et al. (2018) Nakar E., Gottlieb O., Piran T., Kasliwal M. M., Hallinan G., 2018, ApJ, 867, 18
  * Nathanail et al. (2020) Nathanail A., Gill R., Porth O., Fromm C. M., Rezzolla L., 2020, MNRAS, 495, 3780
  * Nedora et al. (2020) Nedora V., et al., 2020, arXiv e-prints, p. arXiv:2008.04333
  * Paczynski (1986) Paczynski B., 1986, ApJ, 308, L43
  * Piro & Kollmeier (2018) Piro A. L., Kollmeier J. A., 2018, ApJ, 855, 103
  * Radice et al. (2018a) Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., Roberts L. F., 2018a, ApJ, 869, 130
  * Radice et al. (2018b) Radice D., Perego A., Hotokezaka K., Bernuzzi S., Fromm S. A., Roberts L. F., 2018b, ApJ, 869, L35
  * Salafia & Giacomazzo (2020) Salafia O. S., Giacomazzo B., 2020, arXiv e-prints, p. arXiv:2006.07376
  * Salafia et al. (2020) Salafia O. S., Barbieri C., Ascenzi S., Toffano M., 2020, A&A, 636, A105
  * Shibata (1999) Shibata M., 1999, Phys. Rev. D, 60, 104052
  * Shibata & Uryū (2000) Shibata M., Uryū K. ō., 2000, Phys. Rev. D, 61, 064001
  * Siegel & Metzger (2017) Siegel D. M., Metzger B. D., 2017, Phys. Rev. Lett., 119, 231102
  * Takahashi & Ioka (2020) Takahashi K., Ioka K., 2020, MNRAS,
  * Troja et al. (2019) Troja E., et al., 2019, MNRAS, 489, 1919
  * Woosley & Heger (2006) Woosley S. E., Heger A., 2006, ApJ, 637, 914
  * Xie et al. (2018) Xie X., Zrake J., MacFadyen A., 2018, ApJ, 863, 58

##  Appendix A Energy composition of the cocoon

The total energy of the cocoon can be written as:

$E_{c}=E_{i}+E_{k,in}+E_{k,0}.$ (54)

$E_{i}$ is the internal energy remaining in the cocoon, out of the energy
injected by the central engine into the cocoon (through the jet head).
$E_{k,in}$ is the part of the kinetic energy of the cocoon that originates
from the energy injected by the engine into the cocoon (even if the injected
energy by the engine is in the form of internal energy, a part of it gets
converted into kinetic energy; inside the shocked jet, before reaching the
cocoon; or later, inside the cocoon, due to the adiabatic expansion of the
cocoon). These two energies ($E_{i}$ and $E_{k,in}$) summed up equal the total
energy delivered by the engine into the cocoon:

$E_{i}+E_{k,in}=(t-t_{0})L_{j}(1-\langle{\beta_{h}}\rangle).$ (55)

Finally, $E_{k,0}$ is the kinetic energy initially carried by the part of the
ambient medium that became the cocoon. It can be calculated with following
integration:

$E_{k,0}=\int_{cocoon}(\Gamma_{a}-1)\rho_{a}dV.$ (56)

Note that, $E_{k,0}=0$ in the collapsar case as the medium is static.

The fraction of the internal energy in the cocoon has been presented in Figure
1. Figure 6 offers more details by showing the fraction of each of the three
energies [in equation (54)], as measured from numerical simulations, and for
each of the four models (see Table 1 for the parameters of each model).

As it can be seen from Figure 6, numerical simulations show that, in the case
of BNS merger jets, the energy contribution of the initial kinetic energy of
the ambient medium ($E_{k,0}$) in the cocoon energy increases with time and is
substantial at later times ($t\sim t_{b}$). This is because of the increasing
volume of the cocoon with time (since, $E_{k,0}\propto V_{c}$). In the
collapsar jet case, the situation is simpler as $E_{k,0}=0$.

In the case of BNS merger jets, the energy contribution of both the internal
energy ($E_{i}$) and the kinetic energy that originates from the central
engine ($E_{k,in}$) declines with time (while it is constant over time in the
collapsar jet case). This is due to two main reasons. First, the increasing
contribution of $E_{k,0}$ in $E_{c}$ with time, which ends up lowering the
share of $E_{i}$ and $E_{k,in}$ in $E_{c}$. Second, the higher jet head
velocity (in the central engine frame) in the BNS merger jet case (relative to
the collapsar jet case) [see equations (43) and (49)]; this implies that in
the BNS merger jet case, and relative to the collapsar jet case, the fraction
of the energy injected by the central engine that ends up in the cocoon is
smaller, and the energy that is contained inside the jet and never reaches the
cocoon is higher [$\langle{\beta_{h}}\rangle$ is higher in the BNS merger case
Ioka & Nakamura 2018 $\sim{2v_{m}}/{c}$; see equation (55)].

Another notable difference between the BNS merger jet case and the collapsar
jet case is the ratio $E_{i}/E_{k,in}$. In the BNS merger jet case, $E_{i}$ is
comparable to $E_{k,in}$ (and even slightly smaller at $t\sim t_{b}$) while in
the collapsar jet case $E_{i}$ is much higher than $E_{k,in}$ ($\sim 80\%$ and
$\sim 20\%$, respectively). This can be explained by the continues expansion
of the medium in the BNS merger jet case, resulting in the conversion of a
substantial part of the internal energy $E_{i}$ into kinetic energy (into
$E_{k,0}$). Roughly, it is estimated that about half of the cocoon’s internal
energy is lost with this process.

To conclude, we showed, for the first time (to our best knowledge), that the
fraction of internal energy in the cocoon, $E_{i}/E_{c}$, is much smaller in
the case of BNS merger jet (relative to the case of collapsar jet) due to
following three reasons: i) the increasing contribution of the initial kinetic
energy of the ambient medium, $E_{k,0}/E_{c}$, in the cocoon of a BNS merger
jet case; ii) the high jet head velocity (in the frame of the central engine)
in the BNS merger jet case; and iii) the adiabatic expansion process in the
cocoon of a BNS merger jet case (due to continues expansion of the ambient
medium).

The fraction $E_{i}/E_{c}$ is an important quantity as it strongly affects the
jet propagation by affecting the value of $\eta$ (and $\eta^{\prime}$)
[Bromberg et al. 2011], hence the importance of these results.

Figure 6: Energy composition of the cocoon, from the jet launch time
$t=t_{0}$, to the jet breakout time $t=t_{b}$, as measured from numerical
simulations. $E_{i}/E_{c}$ is the fraction of the internal energy remaining in
the cocoon (same as the previously presented in the top panel of Figure 1),
$E_{k,in}/E_{c}$ is the fraction of the kinetic energy originating from the
central engine in the cocoon (mostly injected as internal energy, but later
converted into kinetic energy), and $E_{k,0}/E_{c}$ is the fraction of the
kinetic energy that was initially carried by the part of the ambient medium
that became the cocoon (unlike $E_{k,in}$, it did not originate from the
central engine) in the cocoon [see equation (54) and the explanation that
follows]. The top two panels are for BNS merger jet models (from left to
right, T03-H and T13-H), and the bottom two panels are for collapsar jet
models (from left to right, A and B). See Table 1 for more details about each
of the four models.

##  Appendix B Energy composition of the cocoon: testing BNS models with a
hypothetical static ambient medium

Here, our explanation of the origin of the difference in the energy
composition of the cocoon, between the collapsar case and the BNS merger case
(being different due to the adiabatic expansion of the cocoon, and due to the
contribution from the initial kinetic energy of the ambient medium, in the BNS
merger case, hence being due to the expanding aspect of the ambient medium;
see Section 3.2; and Appendix A) is tested. We present two test models
T03-H-static, and T13-H-static. These two models share the same parameters as
the models T03-H, and T13-H (respectively), except that the ambient medium
velocity is set as $v_{m}=0$ (hence the naming “-static”). This test is
carried out to make sure that the overwhelming difference between the BNS
merger case and the collapsar case, in terms of the cocoon’s energy
composition, is not caused by other parameters (such as the parameters of the
jet or the ambient medium; e.g., $L_{j}$, or $M_{a}$).

In Figure 7, the two models T03-H-static and T13-H-static (with $v_{m}=0$),
compared with the two models T03-H and T13-H ($v_{m}=0.345c$), show that the
expanding aspect of the ambient medium is, indeed, the origin of the large
difference in the cocoon’s energy composition (mostly internal for a static
medium and mostly kinetic for an expanding medium). This is in perfect
agreement with our explanation in Section 3.2 and in Appendix A.

Figure 7: Same as the top panel in Figure 1, with the addition of two test
models (T03-H-static and T13-H-static) having the same parameters as the BNS
merger case models (T03-H and T13-H, respectively) expect that their ambient
medium is set as static ($v_{m}=0$, while in T03-H and T13-H $v_{m}=0.345c$;
see Table 1).

##  Appendix C The calibration coefficient $N_{s}$

The calibration coefficient, $N_{s}$, was first introduced by Harrison et al.
(2018) after the realization (first made by Mizuta & Ioka 2013) that the
analytic model by Bromberg et al. (2011) needs to be corrected to match the
simulation results. This parameter is used to correct the analytic (or semi-
analytic) value of $\tilde{L}$ (to $\tilde{L}_{c}$), in accordance with
numerical simulations, as follows:

$\tilde{L}_{c}=N_{s}^{2}\tilde{L}.$ (57)

$N_{s}$ was studied by Harrison et al. (2018) for the case of collapsar jets.
By measuring $\tilde{L}$ in numerical simulations, it was found that $N_{s}$
should take values in the range $\sim\frac{1}{2.5}-\frac{1}{3}$ in the non-
relativistic (Newtonian) regime, and $\sim 1$ in the relativistic regime (see
Figure 12 in Harrison et al. 2018).

Here, the value of $N_{s}$ is deduced based on the jet breakout time; the
value of $N_{s}$ calibrates $\tilde{L}$ (to $\tilde{L}_{c}$) so that the semi-
analytic (or analytic) breakout time matches the numerical simulation’s
breakout time. Hence, with this procedure, the values of $N_{s}$ found here
can be understood as average values – throughout the jet propagation from the
jet launch time ($t_{0}$) to the jet breakout time ($t_{b}$).

In this study, the jet propagation is solved for the case of BNS mergers where
the ambient medium is expanding, in addition to the case of collapsars. We
adopt the analytic model of Bromberg et al. (2011), with a few notable
modifications. In particular, the modeling of the cocoon’s lateral width
($r_{c}$ and $\beta_{\perp}$), the cocoon morphology and volume ($V_{c}$), and
the fraction of internal energy in the cocoon $E_{i}/E_{c}$ (or
$\eta^{\prime}$) have been modified [see equations (9); (17); and Table 1;
respectively]. These modifications are expected to affect the value of
$N_{s}$.

The system of equations of the jet propagation is solved in two methods:

  * •

Semi-analytic solution: the equations are solved through numerical integration
(see Section 2.4).

  * •

Analytical solution: approximations are added to the semi-analytic solution to
simplify it to a fully analytic form (see Section 2.5).

Due to the approximations introduced in the analytic solution, $N_{s}$ is not
exactly the same in these two methods (see Appendix C.2).

###  C.1 $N_{s}$ in the semi-analytic solution

From the semi-analytic modeling (see Section 2.4), we get $N_{s}=0.53$ for the
collapsar case ($\tilde{L}\sim 0.08-1$), and $N_{s}=0.75$ for the BNS merger
case ($\tilde{L}\sim 0.2-0.7$).

First, the value of $N_{s}$ in the collapsar case is slightly different from
the value suggested by Harrison et al. (2018) ($N_{s}\approx 0.33-0.40$ in the
non-relativistic case; see their Figure 12 in the range $\tilde{L}_{a}\lesssim
1$). This can be explained by differences in the modeling [of the cocoon’s
lateral width ($\beta_{\perp}$ and $r_{c}$), the cocoon’s volume ($V_{c}$),
the fraction of internal energy in the in the cocoon ($E_{i}/E_{c}$), and
$\eta$ (or $\eta^{\prime}$)].

Second, $N_{s}$ in the BNS merger case is substantially different than in the
collapsar case. This difference is mainly related to the expanding aspect of
the ambient medium, which is found to strongly affect parameters such as
$E_{i}/E_{c}$ (or $\eta^{\prime}$; see Figure 1), and hence affecting the
cocoon pressure on the jet.

###  C.2 $N_{s}$ in the analytic solution

Using the analytic solution (see Section 2.5) we get $N_{s}=0.38$ for the
collapsar case, and $N_{s}=0.46$ for the BNS merger case.

In the collapsar case, the analytic solution is overall very similar to the
semi-analytic solution. One major difference is the absorption of the term
$(1+\tilde{L}_{c}^{1/2})^{-1}\lesssim 1$ into $N_{s}$ [see Section 2.5.1 and
equation (45)]. This is why $N_{s}$ here differs from $N_{s}$ in the semi-
analytic solution. $N_{s}$ can be estimated for our analytic solution based on
our results in the semi-analytic solution ($\tilde{L}\sim 0.08-1$, and
$N_{s}=0.53$, for models A and B; see Appendix C.1) as follows:

$N_{s}\approx 0.53\>(1+0.53\>\tilde{L}^{1/2})^{-1}.$ (58)

Using these values in the above equation, it can be confirmed that, overall,
$N_{s}\sim 0.3-0.4$.

In the BNS merger case the difference in $N_{s}$ is more substantial;
$N_{s}=0.75$ for the semi-analytic solution, and $N_{s}=0.46$ for the analytic
solution. This is because, in addition to the term
$(1+\tilde{L}_{c}^{1/2})^{-1}$, the term $(1-\beta_{a})$ [usually $\sim
0.6-0.9$] and the term $\frac{1}{\Gamma_{a}}$ have also been absorbed into
$N_{s}$ [see Section 2.5.1 and equation (23); also see equations (3) and
(10)]. Hence, $N_{s}$ for the analytic solution (BNS merger case) can be
estimated using the results of the semi-analytic solution $N_{s}$,
$\beta_{a}$, and $\tilde{L}$, as follows:

$N_{s}\approx
0.75\left[\frac{1-\beta_{a}}{(1+0.75\>\tilde{L}^{1/2})(1-\beta_{a}^{2})^{1/2}}\right],$
(59)

where $\tilde{L}\sim 0.2-0.7$ and $\beta_{a}\sim 0.2$ gives roughly $N_{s}\sim
0.4$ (for T03-H and T13-H; see Appendix C.1) 777It should be noted that more
difference is caused by the approximation introduced in equation (41);
although, this is largely (but not entirely) mitigated by focusing on cases
where $t_{b}-t_{m}\gg t_{0}-t_{m}$..

###  C.3 $N_{s}$ dependence on the parameter space

It should be noted that $N_{s}$ can take different values beyond the parameter
space studied here (i.e., $\tilde{L}\sim 10^{-1}-1$; and $v_{m}=0.35c$).
Therefore, the values of $N_{s}$ given here should not be taken at face value;
$N_{s}$ should always be associated with its parameter space, and calibrated
with simulations (if possible) if used outside of this parameter space.

Nevertheless, we believe that at the close vicinity of our parameter space,
and when using the semi-analytic (or analytic) model presented here, the
values of $N_{s}$ given here should be reliable.