# The propagation of relativistic jets in expanding media

Ore Gottlieb1,2, Ehud Nakar1
1School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel
2Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA),
Physics & Astronomy, Northwestern University, Evanston, IL 60201, USA
<EMAIL_ADDRESS>
###### Abstract

We present a comprehensive analytic model of relativistic jet propagation in
expanding homologous media (ejecta). This model covers the entire jet
evolution as well as a range of configurations that are relevant to binary
neutron star mergers. These include low and high luminosity jets, unmagnetized
and mildly magnetized jets, time-dependent luminosity jets, and Newtonian and
relativistic head velocities. We also extend the existing solution of jets in
a static medium to power-law density media with index $\alpha<5$. Our model
provides simple analytic formulae (calibrated by 3D simulations) for the jet
head propagation and breakout times. We find that the system evolution has two
main regimes: strong and weak jets. Strong jets start their propagation
immediately within the ejecta. Weak jets are unable to penetrate the ejecta at
first, and breach it only after the ejecta expands significantly, thus their
evolution is independent of the delay between the onset of the ejecta and the
jet launching. After enough time, both strong and weak jets approach a common
asymptotic phase. We find that a necessary, but insufficient, criterion for
the breakout of unmagnetized [weakly magnetized] jets is $E_{j,{\rm
iso,tot}}\gtrsim 3[0.4]\,{E_{ej,{\rm
tot}}}\left({\,{\theta_{j,0}}}/{0.1{\rm~{}rad}}\right)^{2}$, where $E_{j,{\rm
iso,tot}}$ is the jet total isotropic equivalent energy, $\,{\theta_{j,0}}$ is
its opening angle, and $\,{E_{ej,{\rm tot}}}$ is the ejecta energy. Applying
our model to short GRBs, we find that there is most likely a large diversity
of ejecta mass, where mass $\lesssim 10^{-3}~{}{\rm M}_{\odot}$ (at least
along the poles) is common.

###### keywords:

gamma-ray bursts — neutron star mergers — methods: analytical — methods:
numerical

††pubyear: 2021††pagerange: The propagation of relativistic jets in expanding
media–B

##  1 Introduction

Relativistic jets appear to be common in a wide range of astrophysical
systems. The jet forms following accretion of infalling mass onto a compact
object through an accretion disk or via fast rotation of a highly magnetized
neutron star. In a variety of astrophysical environments the jet encounters
dense media during its propagation. These can be the infalling mass, outflows
from the compact object or the disk, or some other mass that surrounds the
compact object. The interaction between the jet and the media is an essential
part of the jet evolution, affecting the jet velocity, collimation and
stability. This encounter ultimately dictates whether the jet breaks out from
the media and it strongly affects its structure in case that it does.

The first works that studied the jet-medium interaction considered the
simplest setup of hydrodynamic (unmagnetized) jets propagating in a static
medium, addressing it both analytically (e.g., Blandford & Rees, 1974;
Begelman & Cioffi, 1989; Mészáros & Waxman, 2001; Matzner, 2003; Lazzati &
Begelman, 2005; Bromberg et al., 2011), and numerically (e.g., Marti et al.,
1995; Martí et al., 1997; Aloy et al., 2000; MacFadyen et al., 2001; Zhang et
al., 2004; Mizuta et al., 2006; Morsony et al., 2007; Wang et al., 2008; Rossi
et al., 2008; Lazzati et al., 2009; Mizuta & Aloy, 2009; Morsony et al., 2010;
Nagakura et al., 2011; López-Cámara et al., 2013, 2016; Matsumoto & Masada,
2013, 2019; Mizuta & Ioka, 2013; Ito et al., 2015; Toma et al., 2017; Harrison
et al., 2018; Gourgouliatos & Komissarov, 2018; Gottlieb et al., 2019, 2020a,
2021a). These studies found that the interaction of the jet with the medium
leads to the formation of a hot cocoon that envelops the jet, and has a
fundamental role in the collimation and stability of the jet.

In recent years additional numerical works also included magnetic fields in
order to study how they alter the hydrodynamic picture (e.g., Meliani &
Keppens, 2009; Mizuno et al., 2009, 2012; Tchekhovskoy & Bromberg, 2016;
Bromberg & Tchekhovskoy, 2016; Bromberg et al., 2019; Kathirgamaraju et al.,
2018, 2019; Gottlieb et al., 2020b, 2021b; Gottlieb et al., 2022a; Gottlieb et
al., 2022c; Matsumoto et al., 2021; Nathanail et al., 2021). In general, the
effect of the magnetization depends on the magnetic field configuration (e.g.,
poloidal, toroidal etc.), and its strength compared to the internal energy,
$\sigma\equiv\frac{B^{\prime 2}}{4\pi h\rho c^{2}}$, where $B^{\prime}$ is the
proper magnetic field, $h$ is the specific enthalpy and $\rho$ is the proper
mass density. It was shown that while low $\sigma$ jets
($10^{-2}\lesssim\sigma\lesssim 10^{-1}$) are still hydrodynamically
dominated, their magnetization might be sufficient to stabilize the jet
against hydrodynamic instabilities, thereby supporting faster jets and less
energetic cocoons (Gottlieb et al., 2020b; Gottlieb et al., 2021a; Matsumoto
et al., 2021). In high $\sigma$ jets ($\sigma\gtrsim 1$), current-driven
instabilities such as kink emerge, rendering the jet structure globally
unstable and considerably altering the hydrodynamical picture (e.g., Bromberg
& Tchekhovskoy, 2016; Gottlieb et al., 2022a; Gottlieb et al., 2022c).

The first detection of a neutron star merger, GW170817, by gravitational waves
(GW) was followed by detections across the entire electromagnetic spectrum
(see reviews in Nakar, 2020; Margutti & Chornock, 2020, and references
therein). The kilonova optical signal, powered by radioactive decay of
$r$-process elements, revealed the ejection of subrelativistic massive ($\sim
5\times 10^{-2}~{}\,{\rm M_{\odot}}$) component from the merger, slightly more
massive than the predictions of numerical relativity simulations which find
mass ejection of $10^{-3}-10^{-2}~{}\,{\rm M_{\odot}}$. Months later, the
multi-band afterglow light curve and VLBI observations revealed a relativistic
jet that is similar to the one expected in short gamma-ray burst (sGRB).
According to any reasonable chain of events, at least a non-negligible
fraction of the subrelativistic outflow was ejected before the launching of
the jet. Thus, the observations of GW170817 have shown that the jet has
propagated through the expanding subrelativistic ejecta that generated the
kilonova emission, and it stands to reason that a similar propagation also
takes place in sGRBs that are detected via their $\gamma$-rays. Additional
astrophysical system where there may be a jet that propagates in expanding
medium is hydrogen-poor superluminous supernovae (e.g., Margalit et al.,
2018).

The propagation of relativistic jets in expanding medium started to draw
special attention in anticipation of the first GW signal from binary neutron
star (BNS) mergers (Murguia-Berthier et al., 2014; Murguia-Berthier et al.,
2017; Lazzati et al., 2017; Gottlieb et al., 2018a), and it became the focus
of many studies following GW170817 (Kasliwal et al., 2017; Gottlieb et al.,
2018b; Gottlieb et al., 2021a; Duffell et al., 2018; Kathirgamaraju et al.,
2018; Lazzati & Perna, 2019; Geng et al., 2019; Gottlieb & Loeb, 2020;
Gottlieb & Globus, 2021; Gottlieb et al., 2022b; Ito et al., 2021; Klion et
al., 2021; Murguia-Berthier et al., 2021; Pavan et al., 2021; Urrutia et al.,
2021; Lamb et al., 2022; Nativi et al., 2022). Most of these studies used
numerical simulations showing that the jet-ejecta interplay shapes the
emerging jet-cocoon structure, and thus has important observational
implications on the prompt emission and the afterglow light curve, and
possibly also on the kilonova emission.

First attempts of deriving an analytic model for the propagation of
relativistic jets in expanding media have started recently with the works of
Duffell et al. (2018); Matsumoto & Kimura (2018); Lyutikov (2020); Margalit et
al. (2018); Nakar (2020); Hamidani et al. (2020); Hamidani & Ioka (2021).
These studies have provided a major step-forward in understanding the effect
of expanding medium on the jet propagation. However, the analytic description
in these studies is incomplete and in some of the studies it is also not fully
accurate. Moreover, these studies focused only on part of the evolutionary
regimes and considered only a limited range of the jet-medium properties phase
space, leaving many aspects of the jet propagation unresolved, some of which
are directly relevant for BNS mergers and sGRBs. In this paper we derive an
analytic description of jet propagation within expanding homologous media that
covers all the phases of the evolution and a wide range of initial conditions.
These include: (i) strong jets that breach the ejecta immediately after launch
as well as weak jets that are initially stalled by the ejecta, and can only
penetrate it following a considerable expansion; (ii) the entire evolution of
the jet, starting on timescales that are much shorter than the delay time
between the subrelativistic mass ejection and the jet launching; (iii) the
initial collimation phase of the jet; (iv) relativistic jet head; (v) various
jet magnetizations; (vi) time dependent luminosity. We calibrate our analytic
formula using 3D relativistic magnetohydrodynamic simulations. Previous
analytic models either did not use numerical simulations or made use of
axisymmetric simulations. It was shown in the past that in significant part of
the parameter phase space 3D models produce different results than those of 2D
ones (Gottlieb et al., 2021a), as we also find for our models here. Thus, even
in regimes that were explored by previous studies we obtain different
calibrations.

The paper starts with an overview of the jet dynamics (§2) followed by an
analytic model of the jet propagation in the Newtonian regime through all the
various phases and the entire range of initial conditions (§3). We derive a
very simple breakout criterion and a formula for the breakout time in the
Newtonian regime in §4. We generalize our solution to a time-dependent
luminosity in §5 and to the relativistic head regime in §6. We test our
analytic model using numerical simulations and find calibration coefficients
in §7. In §8 we examine the effect of the density in the cavity, which may
form in the ejecta before the jet is launched, on the jet propagation. In §9
we test numerically the effect of non-negligible sub-dominant magnetization
(i.e., $\sigma\sim 0.01-0.1$) on the jet propagation velocity showing that the
presence of magnetic fields might be equivalent to an increase of an order of
magnitude in the jet luminosity. In §10 we review the results of previous
studies and compare them with ours. In §11 we discuss the implications of our
results on the understanding of the sGRB population in general and on GW170817
in particular. We summarize and conclude in §12.

##  2 An overview of the jet dynamics

###  2.1 Initial conditions

Consider a bi-polar relativistic jet. The jet is launched with a half opening
angle $\theta_{j,0}$ and the jet power is assumed to be constant over the
entire engine work-time (in §5 we relax this assumption) with a total (two-
sided111Note that in many studies of jet propagation (including Bromberg et
al. 2011) the jet luminosity is defined to include only the energy launched in
one direction.) luminosity $L_{j}$. The jet launching starts with a delay
$t_{d}$ with respect to the merger time, which is also assumed to be the time
at which the subrelativistic ejecta is launched. The subrelativistic ejecta is
assumed to be homologous (i.e., $r=vt)$ with a power-law density profile,

$\rho_{ej}=Kr^{-\alpha}t^{\alpha-3}=Kt^{-3}v^{-\alpha}~{}~{}~{};~{}~{}~{}\,{v_{\rm
min}}<v<\,{v_{\rm max}}$ (1)

where $K$, $\,{v_{\rm min}}$ and $\,{v_{\rm max}}$ are constants and $t$ is
the time since the merger. The power-law index at $\,{v_{\rm min}}<v<\,{v_{\rm
max}}$ is constant $\alpha<5$, so that most of the energy is carried by the
fast part of the ejecta. We assume that if there is ejecta with $v>\,{v_{\rm
max}}$ and/or $v<\,{v_{\rm min}}$, then its density is significantly lower
than the power-law extrapolation of the density at $\,{v_{\rm
min}}<v<\,{v_{\rm max}}$. We denote the mass in the velocity range given by
Eq. 1 as the “bulk” of the ejecta and consider the jet as propagating inside
the ejecta only when the ejecta velocity at the jet head location is within
this velocity range. The values of $\,{v_{\rm min}}$, $\,{v_{\rm max}}$ and
$\alpha$ that represent the ejecta in astrophysical explosions may vary
significantly between different systems. Observations of the BNS merger
GW170817, as well as many numerical simulations, indicate that in BNS mergers
$\,{v_{\rm min}}\approx 0.05$c and $\,{v_{\rm max}}\approx 0.3$c (Nakar, 2020,
and references therein). The value of $\alpha$ in BNS mergers is less certain,
although there are some observational and numerical results that suggest that
at least in some mergers $\alpha\approx 3$ gives a decent description of the
ejecta density profile (see §11).

Note that we approximate the ejecta as being spherical. However, also If the
ejecta is not spherically symmetric, then this approximation is good as long
as the angular profile of the ejecta at the poles does not vary strongly over
an angle that is significantly larger than $\,{\theta_{j,0}}$. The adujstment
needed for the application of our model to a non-spherical ejecta are
discussed in §4.4.

###  2.2 General properties of the jet evolution

Before deriving an analytic approximated solution for the propagation of a jet
in such a setup, we describe the properties of the solution in the various
regimes based on general considerations. Here and the following several
sections we restrict our discussion to Newtonian head velocities222Note that
while the jet is launched at relativistic velocities the head itself under a
wide range of conditions propagates at Newtonian velocities.. The propagation
of a jet with a relativistic head is discussed in §6. A key parameter of the
problem is the ratio between the velocity of the jet head as seen in the local
ejecta frame and the ejecta velocity at the location of the head:

$\eta\equiv\frac{v_{h}-\,{v_{ej,h}}}{\,{v_{ej,h}}}~{},$ (2)

where $v_{h}$ is the head velocity and $\,{v_{ej,h}}$ is the ejecta velocity
at the location of the head, both measured in the lab frame. Below we show
that the evolution of $\eta$ with time can be approximated using very general
considerations, and once $\eta(t)$ is known, one can derive the head location
$r_{h}(t)$ by integrating the following equation:

$\frac{d\log(r_{h})}{d\log(\,{\tilde{t}})}=\frac{dr_{h}}{d\,{\tilde{t}}}\frac{\,{\tilde{t}}}{r_{h}}=\frac{v_{h}}{\,{v_{ej,h}}}=\eta+1~{},$
(3)

where here and elsewhere we use normalized time

$\,{\tilde{t}}\equiv\frac{t}{t_{d}}.$ (4)

Then, for a given initial condition of $r_{h}$ at some time
$\,{\tilde{t}}_{0}$, the head location is:

$r_{h}(\,{\tilde{t}})=r_{h}(\,{\tilde{t}}_{0})\exp\left[\int_{\,{\tilde{t}}_{0}}^{\,{\tilde{t}}}\frac{\eta(\,{\tilde{t}}^{\prime})+1}{\,{\tilde{t}}^{\prime}}d\,{\tilde{t}}^{\prime}\right]~{}.$
(5)

To understand the evolution of $\eta$ with time, it is instructive to look at
the two extreme regimes of the parameter $\eta$. When $\eta\gg 1$ the
expansion of the ejecta can be neglected and the propagation is similar to the
one obtained for a jet in a static medium with a power-law density gradient
$\rho\propto r^{-\alpha}$ (Bromberg et al., 2011). In this solution the
evolution of the shock in the Newtonian regime satisfies
$v_{h}\propto(t-t_{d})^{\frac{\alpha-2}{5-\alpha}}$ and
$r_{h}\propto(t-t_{d})^{\frac{3}{5-\alpha}}$, so $v_{h}\propto
r_{h}^{\frac{\alpha-2}{3}}$. At the same time, homologous expansion dictates
$\,{v_{ej,h}}\propto r_{h}$. Thus, as long as $\eta\gg 1$, for $\alpha<5$ the
value of $\eta$ decreases as the head propagates,
$\eta\approx\frac{v_{h}}{\,{v_{ej,h}}}\propto r_{h}^{-\frac{5-\alpha}{3}}$. On
the other hand, when $\eta\ll 1$, the head can be approximated as being at
rest in the local frame of the ejecta. Thus, the density seen by the head of
the jet evolves as $\rho_{ej}(r_{h})\propto r_{h}^{-3}$. The velocity of the
head is determined by finding the frame in which the magnitudes of the
momentum fluxes from the jet and from the ejecta are equal. This condition was
derived by Matzner (2003) for a static medium and modified for an expanding
one by moving to the forward shock local upstream frame, namely, the frame of
the ejecta at the head location (e.g., Ioka & Nakamura, 2018; Hamidani et al.,
2020):

$\frac{v_{h}-\,{v_{ej,h}}}{c}\approx\frac{1}{1+\,{\tilde{L}}^{-1/2}}\approx\,{\tilde{L}}^{1/2}~{}~{}~{}(v_{h},\,{v_{ej,h}}\ll
c),$ (6)

where

$\,{\tilde{L}}\simeq\frac{L_{j}}{\Sigma_{j}\rho_{ej}(r_{h})c^{3}}~{},$ (7)

and $\Sigma_{j}\approx 2\pi r_{h}^{2}\theta_{h}^{2}$ is the jet cross-section
at the location of the head (the factor of 2 accounts for both sides of the
jet to match the two-sided luminosity $L_{j}$) and $\theta_{h}$ is the jet
head opening angle. Thus, for a constant head opening angle,
$\Sigma_{j}\propto r_{h}^{2}$ and if $\eta\ll 1$ then $\,{\tilde{L}}\propto
r_{h}$ so $v_{h}-\,{v_{ej,h}}$, and thus $\eta$, must grow with the radius.
Even if the head opening angle is not constant, it is straight forward to
show333 $\theta_{h}$ depends on the cocoon pressure. If the pressure is high
enough the jet is collimated and $\theta_{h}<\,{\theta_{j,0}}$. Otherwise, the
jet is roughly conical up to the head and $\theta_{h}\approx\,{\theta_{j,0}}$.
Thus, in principle, $\,{\tilde{L}}$ may grow faster or slower than $r_{h}$ if
the jet collimation by the cocoon varies with time. If the collimation
increases, $\,{\tilde{L}}$ grows faster and so does $\eta$. If, instead, the
collimation decreases, then there must be a time where the jet become
uncollimated and $\theta_{h}\approx\theta_{j}$ cannot grow anymore. After this
time $\,{\tilde{L}}$ and $\eta$ are guaranteed to increase with time until
$\eta\sim 1$. that as long as $\eta\ll 1$ it must start growing at some point
until $\eta\sim 1$.

From the discussion above we learn that if $\eta\ll 1$ it must grow with time,
whereas for $\eta\gg 1$ it is decreasing with time. Thus, we expect that
regardless of the initial conditions, after enough time $\eta$ should converge
to a constant asymptotic value, $\eta_{a}$. As we show next, the value of
$\eta_{a}$ can be found by simple dimensional considerations.

###  2.3 The temporal evolution during the asymptotic phase

During the asymptotic phase $t_{d}\ll t$ and therefore the initial timescale,
$t_{d}$, is forgotten and the system has only a single timescale, $t$.
Additionally, during this phase the jet has propagated a significant way in
the ejecta, so $\,{v_{\rm min}}$ is forgotten and there is no velocity scale
in the system. Thus, the only dimensional parameters of the system that are
relevant during the asymptotic phase are $L_{j}$ and the ejecta density
normalization constant $K$. Using dimensional analysis we see that the only
combination that gives a length scale dictates

$r_{h}\propto\left(\frac{L_{j}}{K}\right)^{\frac{1}{5-\alpha}}t^{\frac{6-\alpha}{5-\alpha}}~{}.$
(8)

When $\eta=\eta_{a}$, Eq. 5 reads $r_{h}\propto t^{1+\eta_{a}}$. We therefore
conclude

$\eta_{a}=\frac{1}{5-\alpha}$ (9)

###  2.4 The evolution before the asymptotic regime

Had we been interested only in the asymptotic regime, all we should have done
to provide a solution is finding the normalization factor for Eq. 8 (as we do
in §3.3), which depends only on the dimensionless parameters
$\,{\theta_{j,0}}$ and $\alpha$. However, there may be systems where the
interesting part of the evolution takes place before the asymptotic phase. For
example, there are observational indications suggesting that in GW170817 the
delay time, $t_{d}$, was of order of a second (e.g., Metzger et al., 2018;
Nakar, 2020). In this event, the delay between the GW signal and the gamma-ray
emission was 1.7s. Now, the propagation time of the jet in the merger ejecta
until it breaks out, $t_{bo}$, must satisfy $t_{bo}<1.7s-t_{d}$. Therefore, if
indeed $t_{d}$ was of order of a second then in GW 170817 $t_{bo}\lesssim
t_{d}$. These indications are far from being conclusive but, if correct, then
in GW170817 the jet breakout from the ejecta took place before its propagation
reached the asymptotic phase. Also, it is possible that in some BNS mergers
the jets are choked before ever reaching the asymptotic phase. Therefore, we
devote a significant part of this paper to study the evolution at early times,
before it reaches the asymptotic phase.

The evolution before the asymptotic phase is also divided into two regimes -
the collimation phase, during which the jet is being collimated by the cocoon
pressure, and the evolution after the jet is fully collimated by the cocoon.
As we show here, the collimation phase is, by itself, a transient and its
properties (which may vary from one system to another) do not have a strong
effect on the jet evolution after a full collimation is achieved. On the other
hand, we find that the collimation phase can have an effect on the propagation
of a jet in BNS merger ejecta over a significant part of the jet propagation,
and therefore we need to solve for the evolution during this phase as well.
Unfortunately, this phase may depend on the density in regions close to the
jet launching site where the ejecta may not be homologous and the density may
not follow Eq. 1. Yet, we find a description that constrains the evolution
during this phase for a wide range of configurations (including a non-
homologous ejecta at small radii). Luckily, as we show in §11, for typical
parameters the breakout from the ejecta in BNS mergers takes place near the
end or after the completion of the collimation phase. Thus, the criterion for
a successful breakout that we derive, which is applicable only after the
collimation phase is completed, is largely independent of the unknown
conditions near the jet launching site. Next we describe our strategy to
constrain the evolution during the collimation phase and the conditions under
which it is valid.

####  2.4.1 The effect of the density at $r<\,{r_{\rm min}}$ and non-
homologous ejecta

If there is a significant delay between the time at which the bulk of the
ejecta is launched and the time that the jet is launched, then the jet may
encounter at first a low-density region. In BNS mergers such a delay can be a
result of a delayed collapse of the central object to a black hole, as
observations suggest for GW170817 (Nakar, 2020, and references therein) . We
denote this low-density region as the “cavity” and within our framework it
extends at the time that the jet launch starts, $t_{d}$, between the jet
launching site and $\,{r_{\rm min}}\equiv\,{v_{\rm min}}t_{d}$.

The definition of $\,{v_{\rm min}}$ is the minimal velocity at which a
significant amount of mass is ejected, and therefore the density at
$r<\,{r_{\rm min}}$ is often negligible. However, there are realistic
scenarios where it is not. For example, it is possible that there is ejecta
with $v<\,{v_{\rm min}}$ although its density is low compared to extrapolation
of $\rho\propto v^{-\alpha}$ to low velocities. Another possibility is that
there is a continuous mass ejection so a non-homologous outflow that was
ejected at $t\sim t_{d}$ fills the cavity. Here we assume that if the latter
is the case then the mass ejection rate at $t\sim t_{d}$ is such that the
density in the cavity does not exceed an extrapolation of the power-law of the
bulk of the ejecta to small radii, namely $\rho_{ej}(r<\,{r_{\rm
min}},t_{d})\leq\rho_{ej}(r=\,{r_{\rm min}},t_{d})(r/\,{r_{\rm
min}})^{-\,{\alpha}}$.

The exact solution at early times, while the jet propagates at $r<\,{v_{\rm
min}}t$ and also for a limited time after it starts its propagation in the
bulk of the ejecta, depends on the exact density profile at $r<\,{v_{\rm
min}}t$. This may change from one system to another and therefore we cannot
provide a general solution for that phase. However, we can provide solutions
in two extreme cases, such that all other solutions must lie between them. The
first is an empty cavity where we assume that $\rho(r<\,{r_{\rm
min}},t=t_{d})=0$. This solution is a very good approximation for all
scenarios where the mass ejection rate drops significantly at $t\ll t_{d}$ and
there is a negligible amount of ejecta at $v<\,{v_{\rm min}}$. It is also the
solution where the jet propagates at the highest possible velocity to
$\,{r_{\rm min}}$ and starts its interaction with the bulk of the ejecta
without being collimated at all. The analytic solution that we provide in the
next section assumes an empty cavity and it follows the jet evolution from the
first interaction of the jet with the base of the ejecta at $\,{r_{\rm min}}$
through the collimation phase and up to the asymptotic phase. The second
extreme possibility is the propagation of a jet in a “full cavity”, where the
same power-law density of the ejecta is extended at $t_{d}$ to $r<\,{r_{\rm
min}}$, namely $\rho(r<\,{r_{\rm min}},t=t_{d})=\rho(r=\,{r_{\rm
min}},t=t_{d})(r/\,{r_{\rm min}})^{-\,{\alpha}}$. This solution can be
obtained simply by the empty cavity solution in the limit of $\,{v_{\rm
min}}\rightarrow 0$. Any solution of a partially filled cavity should lie
between these two extreme solutions.

Another assumption of our model that may not be satisfied at $r<\,{r_{\rm
min}}$ is that of a homologous ejecta. This may happen if ejecta continues to
be launched at $t\sim t_{d}$ or later. We expect the effect of a non-
homologous outflow close to the jet and ejecta launching sites to be minor, as
long as the density in the cavity does not exceed our assumed upper limit (an
extension of the power-law to $r<\,{r_{\rm min}}$) and the velocity of freshly
launched ejecta does not increase over time (so there are no internal shocks
that can generate high pressure within the ejecta). The reason is that if the
jet is strong ($\eta\gg 1$), then the motion of the ejecta does not play a
significant role and it is unimportant whether it is homologous or not. If
however the jet is weak compared to the ejecta, then it is stalled by the
ejecta until it expands significantly, at which point the ejecta expansion
becomes homologous and our solution is valid.

Finally, we stress that the ejecta density at $r<\,{r_{\rm min}}$ affects only
the initial phase of the jet evolution. After the jet propagates a significant
way in the bulk of the ejecta, the initial density in the cavity is forgotten.
This happens roughly when $r_{h}\gtrsim\theta_{j,0}^{-2/3}\,{v_{\rm min}}t$,
as we show in §8, where we discuss the effect of a non-empty cavity on the jet
evolution based on analytic considerations as well as numerical simulations.

##  3 Analytic propagation model of a jet with a Newtonian head

We derive an analytic model of a Newtonian jet head propagation in expanding
media and then (§4) deduce its breakout time and a breakout criterion. For
convenience we concentrate in Table 1 the main propagation equations for
$\,{\alpha}=2$ and $\,{\alpha}=4$, including the appropriate numerical
coefficients (which are presented later in §7). In §11 we provide the
propagation equations for $\,{\alpha}=3$, normalized for canonical parameters
of BNS mergers. For convenience we also assemble the main notations used
throughout the paper in Table 2.

Quantity | Regime | $\alpha=2$ | $\alpha=4$
---|---|---|---
$r_{h}$ | $\,{\eta_{\rm crit}}>1$ |  $\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm col}}$: |  $\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm col}}$: |  $\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm col}}$: |  $\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm col}}$:
$\,{v_{\rm min}}t_{d}\left[1+\frac{\theta_{j,0}^{-2/3}-1}{\,{\,{\tilde{t}}_{\rm col}}-1}\left(\,{\tilde{t}}-1\right)\right]$ | $\theta_{j,0}^{-2/3}\left(\frac{\,{\tilde{t}}-1}{\,{\,{\tilde{t}}_{\rm col}}-1}\right)\left(\frac{\,{\tilde{t}}}{\,{\,{\tilde{t}}_{\rm col}}}\right)^{1/4}$ | $\,{v_{\rm min}}t_{d}\left[1+\frac{\theta_{j,0}^{-2/3}-1}{\,{\,{\tilde{t}}_{\rm col}}-1}\left(\,{\tilde{t}}-1\right)\right]$ | $\theta_{j,0}^{-2/3}\left(\frac{\,{\tilde{t}}-1}{\,{\,{\tilde{t}}_{\rm col}}-1}\right)^{3}\left(\frac{\,{\tilde{t}}}{\,{\,{\tilde{t}}_{\rm col}}}\right)^{-3/2}$
$\,{\,{\tilde{t}}_{\rm col}}=2.5[1.3]\,{\eta_{\rm crit}}^{-2/3}$ | $\,{\,{\tilde{t}}_{\rm col}}=9.4[4.7]\,{\eta_{\rm crit}}^{-2/3}\,{\theta_{j,0}}^{-4/9}$
$\,{\eta_{\rm crit}}<1$ |  $\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm a}}$: |  $\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm a}}$: |  $\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm a}}$: |  $\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm a}}$:
$\,{v_{\rm min}}t$ | $\,{v_{\rm min}}t^{4/3}t_{a}^{-1/3}$ | $\,{v_{\rm min}}t$ | $\,{v_{\rm min}}t^{2}/t_{a}$
$\,{\,{\tilde{t}}_{\rm a}}=16[2](\,{\theta_{j,0}}/\,{\eta_{\rm crit}})^{2}$ | $\,{\,{\tilde{t}}_{\rm a}}=844[105](\,{\theta_{j,0}}/\,{\eta_{\rm crit}})^{2}$
Asymptotic $r_{h}$ |  | $0.25[0.5]{\cal L}^{\frac{1}{3}}\,{v_{\rm max}}t^{4/3}$ | $10^{-3}[7\times 10^{-3}]{\cal L}\,{v_{\rm max}}t^{2}$
$\tilde{E}_{d}$ |  | $0.014[0.11]{\cal L}t_{d}$ | $0.006[0.05]{\cal L}t_{d}$
$t_{bo}$ | $\tilde{E}_{d}\gtrsim 1$ | $5[2.5]t_{d}^{2/3}{\cal L}^{-1/3}$ | $5[3]t_{d}^{2/3}{\cal L}^{-1/3}$
$\tilde{E}_{d}\lesssim 1$ | $287[36]{\cal L}^{-1}$ | $662[87]{\cal L}^{-1}$
Breakout criterion |  | $E_{j,{\rm iso,tot}}>150[20]\left((t_{d}/t_{e})^{2}+2\right)\,{E_{ej,{\rm tot}}}\theta_{j,0}^{2}$
Table 1: The main equations that describe the evolution of a subrelativistic head in the cases of $\alpha=2,4$ (including the appropriate numerical coefficients for each $\,{\alpha}$ value), using $\,{\eta_{\rm crit}}=\frac{1}{t_{d}}\bigg{(}\frac{\,{L_{j,\rm iso}}}{8\rho_{ej}(t_{d},\,{r_{\rm min}})\,{v_{\rm min}}^{5}}\bigg{)}^{1/2}$ and ${\cal L}=\frac{\,{L_{j,\rm iso}}}{\,{E_{ej,{\rm tot}}}\theta_{j,0}^{2}}$. If coefficients differ between unmagnetized and mildly magnetized jets, the numbers for magnetized jets are shown in square brackets. The same equations for $\,{\alpha}=3$, normalized for canonical parameters of BNS mergers, are given in §11. For application of these formula to ejecta that is not spherically symmetric see §4.4.  Notation | Definition | Equations
---|---|---
| Jet notations |
$\theta_{j,0}$ | jet launching half-opening angle |
$\,{L_{j,\rm iso}}$ | jet isotropic equivalent luminosity |
$L_{j}$ | jet total (two-sided) luminosity $(\,{L_{j,\rm iso}}\theta_{j,0}^{2}/2)$ |
$\,{E_{j,\rm iso}}(t)$ | $\int_{t_{d}}^{t}\,{L_{j,\rm iso}}~{}dt$ |
$E_{j}(t)$ | $\int_{t_{d}}^{t}L_{j}~{}dt$ |
$\,{E_{j,{\rm tot}}}$ | $\int_{t_{d}}^{t_{d}+t_{e}}L_{j}~{}dt$ |
$t_{e}$ | working time of the jet engine |
$r_{h}$ | jet head location |  5,25,3.3,29,8
$v_{h}$ | jet head velocity |  6,32
| Ejecta notations |
$\,{v_{ej,h}}$ | ejecta velocity at the location of the jet head |
$\,{v_{\rm min}}[\,{\beta_{\rm min}}]$ | minimal ejecta velocity [in units of c] |
$\,{v_{\rm max}}[\,{\beta_{\rm max}}]$ | maximal ejecta velocity [in units of c] |
$\,{r_{\rm min}}$ | $\,{v_{\rm min}}t_{d}$ |
$\,{r_{\rm max}}$ | $\,{v_{\rm max}}t_{d}$ |
$w$ | $\,{v_{\rm max}}/\,{v_{\rm min}}$ |
$\rho_{ej}$ | ejecta density | 1
$K$ | normalization of the ejecta density profile | 1
$A_{\rho}$ | normalization of the ejecta density profile at $t_{d}$ | 10
$\,{m_{\rm min}}$ | $\rho_{ej}(t_{d},\,{r_{\rm min}})\,{r_{\rm min}}^{3}$ |
$\alpha$ | power-law of the ejecta density profile |
$\,{E_{ej}}$ | energy carried by ejecta with $v<\,{v_{ej,h}}$ | 33
$\,{E_{ej,{\rm tot}}}$ | total ejecta energy |
| Jet-ejecta relation notations |
$\eta$ | $(v_{h}-\,{v_{ej,h}})/\,{v_{ej,h}}$ |  2,22,28,33,36
$\eta_{0}$ |  $\eta$ at $\,{r_{\rm min}}$ if the jet is fully collimated | 12
$\,{\eta_{\rm uncol}}$ |  $\eta$ at $\,{r_{\rm min}}$ if the jet is uncollimated | 13
$\,{\eta_{\rm cav}}$ |  $\eta$ when the jet collimation is cavity dominated | 16
$\,{\eta_{\rm crit}}$ | critical $\eta$ to determine if the jet is strong or weak | 17
$\eta_{a}$ |  $\eta$ during the asymptotic regime | 9
$\,{r_{h,{\rm col}}}$ | jet head location when the jet is collimated | 18
$\,{v_{h,{\rm col}}}$ | jet head velocity when the jet is fully collimated | 19
$\,{\eta_{\rm col}}$ |  $\eta$ when the jet becomes fully collimated | 20
$\tilde{E}$ | dimensionless parameter proportional to $\frac{E_{j}}{\,{E_{ej}}}\theta_{j,0}^{-4}$ |  35,38
$\,{\tilde{E}}_{a}$ |  $\tilde{E}$ during the asymptotic regime | 37
$\,{\tilde{E}_{\rm tot}}$ |  $\,{\tilde{E}}$ with the total jet and ejecta energies |  39,40,42
$\tilde{E}_{d}$ |  $\,{\tilde{E}_{\rm tot}}$ where $\,{E_{j,{\rm tot}}}$ is replaced by $L_{j}t_{d}$ | 44
${\cal L}$ | $\,{L_{j,\rm iso}}\,{\theta_{j,0}}^{-2}/\,{E_{ej,{\rm tot}}}$ |
| Time notations |
$t$ | time from the merger |
$t_{d}$ | delay time between the merger and the jet launch |
$\,{\tilde{t}}$ | $t/t_{d}$ | 4
$\delta$ | $t/(t-t_{d})=\,{\tilde{t}}/(\,{\tilde{t}}-1)$ | 23
$\,{\,{\tilde{t}}_{\rm col}}$ | time for jet to reach full collimation | 26
$\,{\,{\tilde{t}}_{\rm a}}$ | time of the transition to the asymptotic regime $\eta_{a}$ |  24,30
$\,{t_{d,\rm eff}}$ | effective $t_{d}$ of weak jets | 31
$t_{bo}$ | breakout time since jet launching |  43,45
$t_{e,bo}$ | minimal engine work-time needed for breakout |
Table 2: Notations throughout the paper.

###  3.1 Critical parameters of the jet evolution

Altogether, under the assumptions described above, the evolution of the system
depends on six parameters: $L_{j}$, $\,{\theta_{j,0}}$, $\,{v_{\rm min}}$,
$t_{d}$, $\,{\alpha}$, and $K$. However, we find that most of the jet
evolution depends only on three parameters, $\,{\theta_{j,0}}$, $\,{\alpha}$
and another dimensionless parameter that we denote $\eta_{0}$, which is
related to the initial value of $\eta$. At very early times there may also be
a minor dependence on $\,{v_{\rm min}}$. Below, we first define $\eta_{0}$
using the fact that for $\eta\gg 1$ the motion of the ejecta can be neglected
and the solution is approaching that of a static medium. Later, we show that
the various phases of the jet evolution depend on various combinations of
$\eta_{0}$, $\,{\theta_{j,0}}$ and $\,{\alpha}$.

Consider a jet that propagates in a hypothetical static medium with a density
profile that is a snapshot of the ejecta density at $t_{d}$, with a single
modification: the density profile of the hypothetical static medium is a
single power-law, $\rho_{stat}\propto r^{-\alpha}$ not only at $r>\,{r_{\rm
min}}$ but also at $r<\,{r_{\rm min}}$. We can therefore write the density
profile of the static medium as444For some values of $\alpha$ the density
profile must have a cut-off at small radius. We assume that this radius is
much smaller than $\,{r_{\rm min}}$ and then its value has a negligible effect
on the hypothetical jet dynamics that we consider in order to define
$\eta_{0}$ when its head reaches $\,{r_{\rm min}}$.
$\rho_{stat}=A_{\rho}r^{-\alpha}$, where this equation is applicable starting
at $r\ll\,{r_{\rm min}}$ and $A_{\rho}$ is related to the expanding medium
density (Eq. 1) via

$A_{\rho}=Kt_{d}^{\alpha-3}=\rho_{ej}(t_{d},\,{r_{\rm min}})\,{r_{\rm
min}}^{\alpha}.$ (10)

Under the assumption that $\rho_{stat}\propto r^{-\alpha}$ is starting at
$r\ll r_{min}$, the jet that propagates in the hypothetical static medium is
fully collimated by the cocoon when the head reaches $\,{r_{\rm min}}$, and
its evolution is described well by Bromberg et al. (2011) and Harrison et al.
(2018). We define the velocity $v_{h,0}$ as the velocity of the jet head that
propagates in the static medium when it reaches $\,{r_{\rm min}}$. $v_{h,0}$
can be found using equations A2 and A3 of Harrison et al. (2018):

$v_{h,0}\approx\frac{0.4}{(5-\alpha)^{1/3}}\left(\frac{L_{j}}{A_{\rho}\theta_{j,0}^{4}}\right)^{1/3}\,{r_{\rm
min}}^{\frac{\alpha-2}{3}}~{}.$ (11)

Now, consider a moving ejecta, with a profile that is similar to that of the
original system, except for the modification that the density profile
$\rho_{ej}=Kt_{d}^{-3}v^{-\alpha}$ is extended also to $v\ll\,{v_{\rm min}}$,
such that $\rho_{ej}(r,t_{d})=\rho_{stat}(r)$. If $v_{h,0}\gg\,{v_{\rm min}}$,
then the jet arrives to $\,{r_{\rm min}}$ at a time which is much shorter than
the ejecta dynamical time, $t_{d}$, namely at $t-t_{d}\ll t$. Thus, as far as
the propagation of the jet head is concerned, the expanding ejecta can be
approximated as static up to the time that the jet head arrives to $\,{r_{\rm
min}}$ with a velocity that is approximately $v_{h,0}$. The ejecta velocity at
this location is approximately $\,{v_{\rm min}}$ so the corresponding value of
$\eta$ at this location, in case that $\eta\gg 1$, can be approximated by Eq.
11. We use this result to define

$\eta_{0}\equiv\left(\frac{L_{j}t_{d}^{\alpha-2}}{A_{\rho}\theta_{j,0}^{4}\,{v_{\rm
min}}^{5-\alpha}}\right)^{1/3}=\left(\frac{L_{j}t_{d}}{\,{m_{\rm min}}v_{\rm
min}^{2}}\right)^{1/3}\theta_{j,0}^{-4/3}.$ (12)

which satisfies, up to a factor of order unity, $\eta_{0}\approx
v_{h,0}/\,{v_{\rm min}}$. In the second expression we rearranged the
parameters to highlight that the dimensionless factor $\eta_{0}$ can be
expressed using the ratio between the injected jet energy over $t_{d}$,
$L_{j}t_{d}$, and the kinetic energy carried by ejecta with velocity
$\sim\,{v_{\rm min}}$, $\,{m_{\rm min}}v_{\rm min}^{2}$, where we define
$\,{m_{\rm min}}\equiv\rho_{ej}(t_{d},\,{r_{\rm min}})\,{r_{\rm min}}^{3}$.

In the real system the value of $\eta$ at $\,{r_{\rm min}}$ is not necessarily
$\eta_{0}$. First, in the derivation of $\eta_{0}$ we neglected the
collimation phase, which depends on the density in the region $r<\,{r_{\rm
min}}$ at $t_{d}$. Second, if $\eta\lesssim 1$, then the static approximation
used to estimate $\eta_{0}$ is not valid and $\eta(\,{r_{\rm
min}})\neq\eta_{0}$. Thus, next we explore the evolution of the jet during the
collimation phase for all values of $\eta_{0}$

####  3.1.1 The collimation phase

Figure 1:  Illustration in arbitrary scales of the different phases in strong
(a,b) and weak (c,d) jets, using the logarithmic specific entropy meridional
maps from simulations (models $\,{\it\alpha_{2}-190}$ and
$\,{\it\alpha_{2}-{0.15}}$, see §7 for details). Left: the jet is partially
collimated and its head is inside the cavity moving at $\,{v_{\rm min}}$.
Right: the jet is collimated and its head is breaching the ejecta from the
cavity. Black dashed lines mark the radius of the cavity, $\,{v_{\rm min}}t$.

The jet evolution depends strongly on the collimation by the cocoon pressure.
However, it takes time for this pressure to build up after the jet launching
starts. Thus, the early jet evolution goes through a phase of collimation,
during which its propagation is affected by the collimation process. The
collimation process depends on the density of the medium that lies between the
launching site and $\,{r_{\rm min}}$. Here we consider a setup where there is
a negligible ejecta mass with $v<\,{v_{\rm min}}$ and that the ejection of
subrelativitic material drops significantly at $t\ll t_{d}$, so that at
$t_{d}$ there is effectively an empty cavity at $r<\,{r_{\rm min}}$ and we can
approximate $\rho_{ej}(t_{d},r<\,{r_{\rm min}})=0$. In §8 we discuss the jet
propagation when the cavity is not empty and the density at $r<\,{r_{\rm
min}}$ affects its evolution.

In an empty cavity, the jet propagates to $\,{r_{\rm min}}$ at a velocity
close to the speed of light. Since there is no pressure in the cavity, the jet
is conical (i.e., uncollimated) up to $\,{r_{\rm min}}$, and its cross-section
when it starts interacting with the medium is $\Sigma_{j}=2\pi r_{\rm
min}^{2}\theta_{j,0}^{2}$ (the factor of 2 accounts for the two sides of the
jet). Plugging this cross-section into Eq. 7, while taking $\rho(\,{r_{\rm
min}},t_{d})$ from Eq. 1, we obtain the value of $\,{\tilde{L}}$. Then we use
Eq. 6 to find the value of $\,{\eta_{\rm uncol}}(t_{d})$ while the head starts
propagating through the dense medium at $\,{r_{\rm min}}$ (we neglect the head
propagation time in the empty cavity). A comparison to $\eta_{0}$ shows that:

$\,{\eta_{\rm uncol}}(\,{\tilde{t}}=1)\approx\beta_{\rm
min}^{1/2}\eta_{0}^{3/2}\,{\theta_{j,0}}~{},$ (13)

where $\,{\beta_{\rm min}}=\,{v_{\rm min}}/c$. If $\,{\eta_{\rm
uncol}}(\,{\tilde{t}}=1)\gg 1$, the jet promptly starts propagating in the
medium. An example of a jet in that phase is shown in Fig. 1a. However, if
$\,{\eta_{\rm uncol}}(\,{\tilde{t}}=1)\ll 1$, the jet is pounding at first on
the slow end of the ejecta without penetrating it almost at all. We focus
first on the latter case, in which the location of the jet head at any given
time is roughly $\,{v_{\rm min}}t$. Before the jet collimation, the jet is
conical and the cross-section of the uncollimated jet head increases as
$\Sigma_{j}\propto r_{h}^{2}\propto\,{\tilde{t}}^{2}$ while the density of the
ejecta at the head location drops as
$\rho_{ej}(r_{h}))\propto\,{\tilde{t}}^{-3}$. Therefore, we obtain that as
long as the jet is uncollimated $\,{\eta_{\rm
uncol}}(\,{\tilde{t}})\propto(\Sigma_{j}\rho_{ej}(r_{h}))^{-1/2}\propto\,{\tilde{t}}^{1/2}$.

As the jet tries unsuccessfully to propagate into the ejecta, the gas that
goes into the jet head (mostly through the reverse shock) spills sideways and
fills the spherical cavity at $r<\,{v_{\rm min}}t$. The pressure in the cavity
increases as the energy deposited by the jet accumulates:

$P\approx\frac{L_{j}(t-t_{d})}{4\pi\,{v_{\rm min}}^{3}t^{3}}.$ (14)

When the pressure is high enough, it starts collimating the jet so the head
cross-section is (Bromberg et al., 2011):

$\Sigma=\frac{L_{j}\,{\theta_{j,0}}^{2}}{4cP}.$ (15)

Plugging Eqs. 15 and 14 into Eqs. 7 and 6, and using the fact that
$\,{v_{ej,h}}\approx\,{v_{\rm min}}$ at this stage, one obtains that once the
jet collimation starts, $\eta$ satisfies,

$\,{\eta_{\rm
cav}}(\,{\tilde{t}})\approx\eta_{0}^{3/2}\,{\theta_{j,0}}(\,{\tilde{t}}-1)^{1/2}.$
(16)

The jet propagates at first at $\eta=\,{\eta_{\rm uncol}}$ and once
$\,{\eta_{\rm cav}}=\,{\eta_{\rm uncol}}$, the collimation starts being
significant and the jet propagates at $\eta=\,{\eta_{\rm cav}}$. Comparing
Eqs. 13 and 16, while taking into consideration the relations $\,{\eta_{\rm
uncol}}(\,{\tilde{t}})\propto\,{\tilde{t}}^{1/2}$ and $\,{\beta_{\rm min}}\ll
1$, one finds that by $\,{\tilde{t}}=2$ all jets with $\,{\eta_{\rm
uncol}}(\,{\tilde{t}}=1)\ll 1$ maintain $\eta=\,{\eta_{\rm cav}}$. This
implies that within one dynamical time of the ejecta since the jet launching
starts, the pressure in the cavity is always high enough to start collimating
the jet. The discussion above implies that we can define a critical value for
the system parameters that dictates its evolution:

$\,{\eta_{\rm crit}}\equiv\,{\eta_{\rm cav}}(\,{\tilde{t}}=2)=\,{N_{\rm
crit}}\eta_{0}^{3/2}\,{\theta_{j,0}}=\,{N_{\rm
crit}}\left(\frac{L_{j}\theta_{j,0}^{-2}t_{d}}{\,{m_{\rm min}}v_{\rm
min}^{2}}\right)^{1/2}~{},$ (17)

where $\,{N_{\rm crit}}$ is a numerical coefficient of order unity that we
extract from the simulations described in §7. If $\,{\eta_{\rm crit}}\ll 1$,
the jet is unable to propagate into the ejecta for many dynamical times. In
that case $\eta$ grows slowly until it reaches $\eta_{a}$. We refer to this
regime as the “weak jet” regime. If, however, $\,{\eta_{\rm crit}}\gg 1$, then
$\eta$ becomes much larger than unity on a time that is much shorter than the
ejecta dynamical time ($t_{d}$) and the jet starts its propagation through the
ejecta almost immediately after it is launched. In such a case, at early time
$\eta\gg 1$ and it starts decreasing a short time after the jet launch starts,
until it reaches $\eta_{a}$. We refer to this case as the “strong jet” regime.
We note that although we use the terms ’strong’ and ’weak’ to describe the
initial propagation regime of a jet, the determination of the regime of a
specific system depends not only on the jet luminosity but also on its opening
angle, ejecta properties and the delay time. Thus, even in case that we
observe the jet and the ejecta, e.g., via the afterglow and kilonova emission,
we may be still unable to determine whether the jet propagation was in the
strong or weak regime unless we know what $t_{d}$ is. In the following we
provide an analytic solution to the jet evolution in the strong and the weak
regimes.

###  3.2 Strong jet - $\,{\eta_{\rm crit}}\gg 1$

When $\,{\eta_{\rm crit}}\gg 1$ the jet starts propagating in the ejecta at
$\,{\tilde{t}}-1\ll 1$ with $\eta\gg 1$ (see Fig. 1a). Thus, $\eta$ approaches
its asymptotic value, $\eta_{a}$, from above, and as long as $\eta>1$, we can
use the approximation of a jet that propagates in a static medium. However, at
first the size of the cavity at $r<\,{r_{\rm min}}$ affects the collimation of
the jet. Therefore, the propagation approaches the solution of Bromberg et al.
(2011) only when the jet head reaches

$\,{r_{h,{\rm col}}}\approx\,{r_{\rm min}}\theta_{j,0}^{-2/3}~{},$ (18)

and the volume of the bubble inflated by the cocoon, $\sim 2\pi
r_{h}^{3}\theta_{j,0}^{2}$ (Bromberg et al., 2011), is comparable to the
volume of the cavity, $(4/3)\pi r_{\rm min}^{3}$. Only at $r\gtrsim r_{h}$ the
volume of the cavity has a minor effect, so the collimation process is
complete and we can refer to the jet as being “fully collimated”. An example
of a jet in that phase is shown in Fig. 1b. Using equations A2 and A3 from
Harrison et al. (2018) we estimate that at $\,{r_{h,{\rm col}}}$,

$\,{v_{h,{\rm col}}}\approx v_{h,0}(\,{r_{h,{\rm col}}}/\,{r_{\rm
min}})^{\frac{\,{\alpha}-2}{3}}~{}.$ (19)

Neglecting the ejecta motion between $t_{d}$ and the time that the jet head
reaches $\,{r_{h,{\rm col}}}$, we find that $\,{v_{ej,h}}(\,{r_{h,{\rm
col}}})\approx\,{v_{\rm min}}\,{r_{h,{\rm col}}}/\,{r_{\rm min}}$ and
therefore we can approximate,

$\,{\eta_{\rm col}}\equiv\eta(\,{r_{h,{\rm col}}})\approx\,{N_{\rm
col}}(\alpha)\frac{0.4}{(5-\alpha)^{1/3}}\eta_{0}\theta_{j,0}^{\frac{2(5-\,{\alpha})}{9}}~{},$
(20)

where $\,{N_{\rm col}}(\alpha)$ is a numerical coefficient that we find from
the simulations. At $r_{h}>\,{r_{h,{\rm col}}}$, as long as $\eta\gg 1$, the
solution of a jet in a static media provides a good approximation and
$v_{h}\propto(t-t_{d})^{\frac{\alpha-2}{5-\alpha}}$. We can therefore relate
the distance that the head traveled up to time $t$ and the distance traveled
by the ejecta at the same location

$\,{v_{ej,h}}t\approx\int_{t_{d}}^{t}v_{h}(t^{\prime})dt^{\prime}\approx\frac{5-\alpha}{3}v_{h}(t)(t-t_{d})~{},$
(21)

where the last approximation is accurate for $r\gg\,{r_{h,{\rm col}}}$.
Rearranging Eq. 21 and adding a constant so
$\eta(\,{\tilde{t}}\to\infty)\to\eta_{a}$ we obtain

$\eta\approx\frac{3}{5-\alpha}(\delta-1)+\eta_{a}=\frac{3\delta-2}{5-\alpha}~{}~{}~{};~{}~{}~{}\eta<\,{\eta_{\rm
col}},$ (22)

where

$\delta\equiv\frac{t}{t-t_{d}}=\frac{\,{\tilde{t}}}{\,{\tilde{t}}-1}.$ (23)

Eq. 22 has two interesting implications. First, at $r>\,{r_{h,{\rm col}}}$
there is a universal solution for $\eta(\,{\tilde{t}})$ which is applicable
for all strong jets. The difference between different jets is the time and the
value of $\eta$ at which they join this solution (given roughly by Eq. 20).
Second, in all strong jets $\eta$ approaches its asymptotic value on a
timescale that is not much larger than $t_{d}$ after the jet launching starts.
We denote the time that the jet approaches $\eta_{a}$ as $t_{a}$. From Eq. 22
we see that $\eta_{a}<\eta(\,{\tilde{t}}>4)<2\eta_{a}$, implying that in
strong jets

$\,{\,{\tilde{t}}_{\rm a}}\approx 4~{}~{}~{};~{}~{}~{}\,{\eta_{\rm crit}}\gg
1.$ (24)

The evolution of $\eta$ before the jet becomes fully collimated is complex and
among other things, it depends on the density in the cavity (whether it is
negligible, as in an empty cavity that we consider in this section or not). We
therefore take here the simplest possible approximation during that phase, a
constant velocity. The entire evolution can be then approximated based on Eqs.
5, 9, 18 and 22:

$r_{h}(\,{\tilde{t}})=\,{r_{\rm
min}}\begin{cases}1+\frac{\theta_{j,0}^{-2/3}-1}{\,{\,{\tilde{t}}_{\rm
col}}-1}(\,{\tilde{t}}-1)&\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm col}}\\\ &\\\
\theta_{j,0}^{-2/3}\bigg{(}\frac{\,{\tilde{t}}-1}{\,{\,{\tilde{t}}_{\rm
col}}-1}\bigg{)}^{\frac{3}{5-\alpha}}\Big{(}\frac{\,{\tilde{t}}}{\,{\,{\tilde{t}}_{\rm
col}}}\Big{)}^{\frac{3-\alpha}{5-\alpha}}&\,{\,{\tilde{t}}_{\rm
col}}<\,{\tilde{t}}\end{cases}~{},$ (25)

where

$\,{\,{\tilde{t}}_{\rm col}}\approx\frac{3}{(5-\,{\alpha})\,{\eta_{\rm
col}}}+1~{},$ (26)

so Eq. 22 is approximately satisfied at $t>\,{\,{\tilde{t}}_{\rm col}}$. As we
show in §7, using numerical simulations that scan a wide range of parameters,
Eq. 25 provides a good approximation for $r_{h}(t)$ of jets with $\,{\eta_{\rm
crit}}>\eta_{a}$. For very strong jets with $\,{\eta_{\rm crit}}\gg 1$ the
approximation is very good, to within a factor of order unity in all our
simulated range, while for jets with $\,{\eta_{\rm crit}}\sim 1$ it provides a
fairly good approximation, to within a factor of 2.

###  3.3 The asymptotic phase

Regardless of whether the jet is strong or weak, after enough time the jet
must approach the asymptotic phase. Moreover, since the evolution in this
phase is independent of $t_{d}$ and $\,{v_{\rm min}}$, all jets, strong and
weak, with the same value of $L_{j}/K$, $\,{\theta_{j,0}}$ and $\,{\alpha}$
must eventually converge to the same late time evolution. Thus, we can use the
late time evolution of strong jets (at $t\gg t_{d}$) to find the location of
the head during the asymptotic phase. Note, however, that in reality some jets
break out before reaching the asymptotic phase. This is expected for example
in BNS mergers where $\,{v_{\rm max}}/\,{v_{\rm min}}\approx 6$ and strong
jets are expected to break out on a timescale $\lesssim t_{d}$.

Taking the limit of $t\gg t_{d}$ in Eq. 25 we find that during the asymptotic
phase

$\displaystyle r_{h}$
$\displaystyle=\left(\frac{2(5-\,{\alpha})^{\frac{2}{3}}\,{N_{\rm
col}}}{15}\right)^{\frac{3}{5-\,{\alpha}}}\left(\frac{L_{j}}{K}\theta_{j,0}^{-4}\right)^{\frac{1}{5-\,{\alpha}}}t^{\frac{6-\,{\alpha}}{5-\,{\alpha}}}$
$\displaystyle=\left(\frac{2(2\pi(5-\,{\alpha}))^{\frac{1}{3}}\,{N_{\rm
col}}}{15}\right)^{\frac{3}{5-\,{\alpha}}}\left(\frac{L_{j}}{\,{E_{ej,{\rm
tot}}}}\theta_{j,0}^{-4}\right)^{\frac{1}{5-\,{\alpha}}}\,{v_{\rm
max}}t^{\frac{6-\,{\alpha}}{5-\,{\alpha}}}~{},$ (27)

where $\,{E_{ej,{\rm tot}}}$ is the total ejecta kinetic energy. As expected
this expression follows the scaling of Eq. 8 and it is independent of $t_{d}$
and $\,{v_{\rm min}}$. It is a good approximation for strong jets at $t\gtrsim
4t_{d}$. An example of a strong jet that approaches the asymptotic phase can
be seen in the left panels of Fig. 2. The time at which it becomes a good
approximation for weak jets is found next.

###  3.4 Weak jet - $\,{\eta_{\rm crit}}<\eta_{a}$

As long as $\eta\ll 1$ the jet is too weak to breach the ejecta. Instead, the
head propagates roughly at $\,{v_{\rm min}}$ and the volume of the cocoon is
dominated by the cavity at $r<\,{v_{\rm min}}t$ so $\eta\approx\,{\eta_{\rm
cav}}$. Thus, $\,{\eta_{\rm crit}}=\,{\eta_{\rm cav}}(\,{\tilde{t}}=2)\ll 1$
and the jet starts to propagate into the ejecta only at $\,{\tilde{t}}\gg 1$
(i.e., $t\gg t_{d}$). An example of a weak jet that is stalled at the base of
the ejecta on a timescale that is longer than $t_{d}$ is shown in Fig. 1c. We
can use Eq. 16 to approximate

$\eta(\,{\tilde{t}})\approx\,{\eta_{\rm
crit}}(\,{\tilde{t}}-1)^{1/2}~{}~{}~{};~{}~{}~{}\,{\tilde{t}}\ll\,{\,{\tilde{t}}_{\rm
a}}~{}.$ (28)

As $\eta$ approaches $\eta_{a}$, the jet starts its propagation inside the
ejecta and the collimation is not determined by the cavity alone so Eq. 28 is
no longer accurate. An example of a weak jet that starts breaching the ejecta
is shown in Fig. 1d. We know that before the jet starts propagating in the
ejecta the head velocity satisfies $v_{h}\approx\,{v_{\rm min}}$, while after
it starts to propagate in the ejecta it approaches the asymptotic solution
given in Eq. 3.3. Therefore we can approximate the head location of a weak jet
as

$r_{h}(\,{\tilde{t}})\approx\,{r_{\rm
min}}\,{\tilde{t}}\begin{cases}1&\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm a}}\\\
&\\\ \Big{(}\frac{\,{\tilde{t}}}{\,{\,{\tilde{t}}_{\rm
a}}}\Big{)}^{\frac{1}{5-\alpha}}&\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm a}}\\\
\end{cases}~{},$ (29)

where

$\,{\tilde{t}}_{a}=\left(\frac{15}{2(5-\,{\alpha})^{\frac{2}{3}}\,{N_{\rm
col}}}\right)^{3}\eta_{0}^{-3}~{}~{}~{};~{}~{}~{}\,{\eta_{\rm crit}}\ll 1~{}$
(30)

is chosen so the propagation during the asymptotic phase satisfies Eq 3.3. An
example of a weak jet that approaches the asymptotic phase is shown on the
right panels of Fig. 2. Our numerical simulations show that Eq. 29 provides a
very good approximation of $r_{h}$ for weak jets, to within a factor of order
unity during the entire simulated range.

####  3.4.1 Effective delay time of weak jets

Since a weak jet is unable to propagate a significant way within the ejecta up
to $t_{a}$, which is much longer than $t_{d}$, there is a little difference
between a jet that starts being launched at $t_{d}\ll t_{a}$ and a jet that
starts being launched at $t_{d}\approx t_{a}$. To see that consider two
similar jets (in terms of luminosity and opening angle), one is a weak jet
that starts being launched at $t_{d}\ll t_{a}$ and the other starts being
launched at $t_{d}=t_{a}$. In the first jet the pressure in the cavity builds
up and collimate the jet on a timescale comparable to $t_{d}$, but the ejecta
is too dense so the jet is stalled at first. Only at $t\sim t_{a}$ the ejecta
density drops enough so the jet can start its propagation within the ejecta.
In the second jet, the cavity pressure starts building up only at $t_{a}$, but
at $t\approx 2t_{a}$ it becomes comparable to that of the first jet at the
same time, so the propagation of the two jets becomes similar. This implies
that the propagation of weak jets is only weakly dependent on $t_{d}$ and it
depends mostly on $t_{a}$ and all weak jets behave as if $t_{d}\sim t_{a}$. We
therefore define $t_{a}$ as an effective delay time for weak jets. Following
Eq. 30

$\,{t_{d,\rm
eff}}=t_{a}=\left(\frac{15}{2(5-\,{\alpha})^{\frac{2}{3}}\,{N_{\rm
col}}}\right)^{3}\frac{\,{m_{\rm min}}v_{min}^{2}}{L_{j}}\theta_{j,0}^{4}~{}.$
(31)

This result has an interesting implication for BNS mergers. As $t_{d}$ depends
on the post merger evolution, it is an important quantity that may teach us
about the physics of mergers. Unfortunately, the discussion above shows that
in case of weak jets the observations are almost independent of $t_{d}$.

###  3.5 The structure of the jet-cocoon

Figure 2:  Meridional maps of our canonical strong $\,{\it\alpha_{2}-190}$
(left) and marginal $\,{\it\alpha_{2}-1}$ (right) jet models (taken from the
low resolution simulations, see §7 for details) close to the asymptotic phase,
at $\,{\tilde{t}}=5.7$ and $\,{\tilde{t}}=41$, respectively. Top panels: the
logarithmic mass density (in arbitrary units) map, where a dashed black line
is stretched from the origin at the jet opening. It shows that in both
simulations the cocoon opening angle is
$\theta_{c}=r_{h}/l_{c}\approx\,{\theta_{j,0}}$. Bottom panels: logarithmic
proper-velocity.

The solution of Bromberg et al. (2011) to a jet in a static medium includes
the full structure of the jet-cocoon system during the propagation. This
includes the cocoon opening angle, $\theta_{c}$, and the Lorentz factor of jet
material after it is shocked by the collimation shock and before it reaches
the head, $\Gamma_{1}$ (we use here the notation as in Bromberg et al. 2011).
$\theta_{c}$ is of interest since it defines the part of the ejecta that
affects the jet propagation. $\Gamma_{1}$ is of interest since it can affect
the emission from GRB jets (e.g., the efficiency of the photospheric emission,
Gottlieb et al. 2019). Here we extend the solution of these quantities to a
jet in an expanding medium.

As in most other properties, the discussion of the cocoon opening angle
depends on the value of $\eta$. Here we focus on the asymptotic phase since in
the regime of $\eta\ll 1$, the cocoon is confined to the cavity. If the cavity
is empty, then the cocoon structure is trivial (approximately a sphere with
radius $\,{v_{\rm min}}t$), otherwise the structure depends on the cavity
details. On the other hand, for $\eta\gg 1$ the solution can be approximated
by the static case and therefore the structure derived by Bromberg et al.
(2011) is applicable.

Similarly to the value of $\eta_{a}$, the evolution of the jet-cocoon
structure during the asymptotic phase can be deduced from dimensional
analysis. Since there is only a single combination of the system parameters
that can give a length scale during the asymptotic phase (see §2.3), all
length scales must have the same dependence on time. This includes the cocoon
cylindrical radius, $l_{c}$, which evolves as $l_{c}\propto r_{h}\propto
t^{\frac{6-\,{\alpha}}{5-\,{\alpha}}}$. Thus, all length scale ratios,
including $\theta_{c}=(l_{c}/r_{h})$, remain constant during the asymptotic
phase, at least until the head becomes relativistic. The way to evaluate
$\theta_{c}$ is by noting that one of the pathways to the asymptotic phase is
via a strong jet that starts with $\eta\gg 1$. According to the static
solution, during this phase the cocoon opening angle is always constant in
time and it satisfies $\theta_{c}\approx\,{\theta_{j,0}}$. Thus, assuming that
$\theta_{c}$ does not change much during the transition from $\eta\gg 1$ to
$\eta=\eta_{a}$ we can approximate $\theta_{c}\approx\,{\theta_{j,0}}$ also
during the asymptotic phase. We verify that this relation holds in the three
simulations that we have that reach this phase (strong, marginal and weak), as
can be seen in Fig. 2 in two examples. This result is of special interest in
the context of BNS mergers where the ejecta is most likely not isotropic. It
implies that during the entire evolution the jet-cocoon structure is sensitive
only to the density within an opening angle $\sim\,{\theta_{j,0}}$ around the
jet symmetry axis. Given that typical jet angles are considered to be of the
order of $0.1$ rad (e.g., in GW170817 it is most likely $\approx 0.05$ rad),
the solution presented here is applicable even if the ejecta is highly
aspherical as long as its density and velocity distributions do not vary much
over an angle of $\,{\theta_{j,0}}$ around the jet axis.

Another interesting property of the jet-cocoon structure is the Lorentz factor
above the first collimation shock. Bromberg et al. (2011) showed that in
static media $\Gamma_{1}\approx 1/\,{\theta_{j,0}}$. This relation arises from
the requirement that after the collimation the pressure in the shocked jet is
similar to the cocoon pressure. The same equality holds also in expanding
media in all of the various regimes, and therefore also here
$\Gamma_{1}\approx 1/\,{\theta_{j,0}}$. In Fig. 2 we show maps of two jets
that are approaching the asymptotic phase, one very strong and the other
marginally strong. Those maps show that in both cases
$\theta_{c}\approx\,{\theta_{j,0}}$ and $\Gamma_{1}\approx
1/\,{\theta_{j,0}}$.

##  4 A Breakout criterion and the dimensionless parameter $\,{\tilde{E}}$

The ability of a jet to break out from the merger ejecta has profound
observational implications. Therefore, it is most interesting to characterise
the necessary conditions for breakout. One way to find analytically if a jet
breaks out of the ejecta is to use the approximated formula for $r_{h}$ given
above (Eqs. 25 and 29) and check whether at the time that the engine stops the
jet has broken out or not. However, as we show below if we are only interested
in whether the jet breaks out or not, then there is a much simpler criterion
that provides an answer. To find this criterion we define a new dimensionless
parameter that we denote $\,{\tilde{E}}$, which depends on the energy ratio
between the jet and the ejecta and on the jet opening angle. Below we find
first the relation between $\,{\tilde{E}}$ and $\eta$.

###  4.1 Definition of $\,{\tilde{E}}$ and the relation with $\eta$

####  4.1.1 Fast jet head - $\eta\gg 1$

When $\,{\eta_{\rm crit}}\gg 1$, the jet head is initially fast, namely
$v_{h}\gg\,{v_{ej,h}}$ and $\eta\gg 1$. As long as this condition is
satisfied, and after the jet is fully collimated
($\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm col}}$) we can use the equations of
propagation in a static medium as a good approximation. We start from equation
A3 of Harrison et al. (2018) written in the instantaneous ejecta rest frame:

$v_{h}-\,{v_{ej,h}}=\bigg{(}\frac{\tilde{N}_{s}^{5}}{(5-\alpha)^{3-\alpha}}\frac{L_{j}}{A_{\rho}\theta_{j,0}^{4}}\bigg{)}^{\frac{1}{5-\alpha}}(t-t_{d})^{\frac{\alpha-2}{5-\alpha}}~{},$
(32)

where $\tilde{N}_{s}$ is a numerical coefficient that is related to, but
different than, the coefficient $N_{s}$ that was found for the case of a
static medium (Harrison et al., 2018) and later verified for an expanding
ejecta (Hamidani et al., 2020). We discuss its role and calculate its value
for various $\,{\alpha}$ in Appendix A. Dividing both sides by $\,{v_{ej,h}}$
yields:

$\eta=\bigg{(}\frac{\tilde{N}_{s}^{5}}{(5-\alpha)^{3-\alpha}}\frac{L_{j}(t-t_{d})^{\alpha-2}t^{5-\alpha}}{A_{\rho}\theta_{j,0}^{4}r^{5-\alpha}}\bigg{)}^{\frac{1}{5-\alpha}}=\bigg{(}\frac{2\pi\tilde{N}_{s}^{5}}{(5-\alpha)^{4-\alpha}}\frac{\delta^{3-\alpha}E_{j}}{\,{E_{ej}}\theta_{j,0}^{4}}\bigg{)}^{\frac{1}{5-\alpha}}~{},$
(33)

where $E_{j}=L_{j}(t-t_{d})$ is the energy injected into the jet up to time
$t$ and

$\,{E_{ej}}(v<\,{v_{ej,h}})=\frac{2\pi}{5-\alpha}A_{\rho}r_{h}^{5-\alpha}t^{-2}$
(34)

is the kinetic energy carried by the ejecta555Note that we assume here
spherically symmetric ejecta. For application to an ejecta that is not
spherically symmetric see §4.4 within $v<\,{v_{ej,h}}=r_{h}/t$.

We then define a dimensionless parameter of the system that evolves as the jet
propagates and depends mostly on the jet to ejecta instantaneous energy ratio
and the jet opening angle:

$\,{\tilde{E}}(\,{\tilde{t}},v)\equiv\frac{2N_{E}^{5}(\alpha)}{(5-\alpha)^{4-\alpha}}\frac{E_{j}(\,{\tilde{t}})}{\,{E_{ej}}(\,{v_{ej,h}})}\theta_{j,0}^{-4}~{},$
(35)

where $N_{E}(\alpha)$ is a numerical coefficient (obtained from the
simulations), which includes $\tilde{N}_{s}$. Then Eq. 33 can be written in a
compact form:

$\eta(\,{\tilde{t}})\approx\big{[}\delta(\,{\tilde{t}})^{3-\alpha}\,{\tilde{E}}(\,{\tilde{t}},v)\big{]}^{\frac{1}{5-\alpha}}~{}~{}~{};~{}~{}~{}\eta\gg
1~{},~{}\,{\,{\tilde{t}}_{\rm col}}<\,{\tilde{t}}<\,{\,{\tilde{t}}_{\rm
a}}~{}.$ (36)

This approximation is valid up to the time where $\eta\sim 1$ and
$\,{\tilde{E}}\sim 1$. Note that Eqs. 33-36 are all depending on the value of
$\frac{E_{j}}{\,{E_{ej}}\theta_{j,0}^{4}}$, which shows a non-trivial
dependence on $\,{\theta_{j,0}}$. The reason for this dependence is that the
relevant parameter on which $\eta$ depends is the ratio of the jet and ejecta
energies per unit of area along the path of the jet. Thus, if the collimation
of the jet by the cocoon is (wrongfully) ignored, one obtains that Eqs. 33-36
depend on
$\frac{E_{j}}{\,{E_{ej}}\theta_{j,0}^{2}}\left(\approx\frac{E_{j,iso}}{\,{E_{ej}}}\right)$,
in which case the dependence on $\,{\theta_{j,0}}$ corresponds to the energy
per unit of area of an uncollimated jet. When the collimation, which by itself
depend on $\,{\theta_{j,0}}$, is taken into account, one obtains the
dependence on $\frac{E_{j}}{\,{E_{ej}}\theta_{j,0}^{4}}$.

####  4.1.2 The asymptotic regime - $\eta=\eta_{a}$

During the asymptotic phase the jet evolution is independent of whether it
started strong or weak, and thus also the value of $\,{\tilde{E}}$ during the
asymptotic phase is similar for both types of jets. The evolution of
$\,{\tilde{E}}$ during this phase can be found analytically. When
$\eta=\eta_{a}$ the head radius increases as $t^{1+\eta_{a}}$ implying that
$\,{v_{ej,h}}\propto t^{\eta_{a}}$ and therefore
$\,{E_{ej}}(\,{v_{ej,h}})\propto\,{v_{ej,h}}^{5-\alpha}\propto
t^{(5-\alpha)\eta_{a}}\propto t$. Since $E_{j}\propto t$ as well, we obtain
that $\,{\tilde{E}}$ is constant during the asymptotic phase:

$\,{\tilde{E}}=\,{\tilde{E}}_{a}~{}~{}~{};~{}~{}~{}t\gg t_{a}.$ (37)

Below, we use numerical simulations to find the value of $\,{\tilde{E}}_{a}$,
which is indeed of order unity with some dependence on $\alpha$ (see Table 4).

####  4.1.3 Slow jet head - $\eta<\eta_{a}$

When $\,{\eta_{\rm crit}}\ll\eta_{a}$, the jet head propagates very slowly up
to $t\sim t_{a}$ without propagating a significant distance inside the ejecta.
For our purpose, it is not very interesting to follow the value of
$\,{\tilde{E}}$ during this phase. At $t\gtrsim t_{a}$ the jet propagation
merges with the asymptotic regime discussed above.

###  4.2 Breakout criterion

The conclusion of the previous section is that after the jet propagates a
significant way in the ejecta (i.e., $\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm
col}}$ in strong jets and $\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm a}}$ in weak
jets), we can use the approximation

$\,{\tilde{E}}(\,{\tilde{t}})\approx\delta^{\alpha-3}(\delta-1)^{5-\alpha}+\,{\tilde{E}}_{a}\approx(\,{\tilde{t}}-1)^{-2}+\,{\tilde{E}}_{a}~{}.$
(38)

In strong jets this approximation is obtained from Eqs. 22, 36 and 37 when we
neglect the coefficient that depends on $\alpha$ in Eq. 22 and use the fact
that $\delta(\,{\tilde{t}}\gg 1)-1\ll 1$. In weak jets it is derived using the
fact that the jet propagates in the ejecta only at
$\,{\tilde{t}}\gtrsim\,{\,{\tilde{t}}_{\rm a}}\gg 1$ when $\delta-1\ll 1$ and
$\,{\tilde{E}}\approx\,{\tilde{E}}_{a}$.

Based on that approximation to $\,{\tilde{E}}(\,{\tilde{t}})$ we can find a
breakout criterion. We define the following parameters:

$\,{\tilde{E}_{\rm
tot}}\equiv\frac{2N_{E}^{5}}{(5-\alpha)^{4-\alpha}}\frac{\,{E_{j,{\rm
tot}}}}{\,{E_{ej,{\rm
tot}}}}\theta_{j,0}^{-{4}}=\frac{N_{E}^{5}}{(5-\alpha)^{4-\alpha}}\frac{\,{\tilde{E}}_{j,{\rm
iso,tot}}}{\,{E_{ej,{\rm tot}}}}\theta_{j,0}^{-{2}}~{},$ (39)

where $\,{E_{j,{\rm tot}}}=L_{j}t_{e}$ is the total jet energy,
$\,{\tilde{E}}_{j,{\rm iso,tot}}=2\,{E_{j,{\rm tot}}}\theta_{j,0}^{-2}$ is the
total jet isotropic equivalent energy, and $t_{e}$ is the jet engine work-
time. The jet must propagate a significant way within the ejecta in order to
break out and therefore Eq. 38 holds at the time of the breakout, when it can
be written as
$\,{\tilde{E}}(\,{\tilde{t}}=\,{\tilde{t}}_{bo}+1)\approx\,{\tilde{t}}_{bo}^{-2}+\,{\tilde{E}}_{a}~{}$
where $\,{\tilde{t}}_{bo}\equiv t_{bo}/t_{d}$, and $t_{bo}$ is the duration
that the engine works until the jet breakout (the jet breaks out at time
$t_{d}+t_{bo}$ after the merger). A successful jet requires $t_{e}>t_{bo}$ and
$\,{\tilde{E}}_{\rm tot}>\,{\tilde{E}}(\,{\tilde{t}}=\,{\tilde{t}}_{bo}+1)$.
Therefore a criterion for a successful breakout is:

$\displaystyle\,{\tilde{E}}_{\rm tot}>\,N_{\tilde{E}_{\rm
tot}}(\alpha)\,{\tilde{t}}_{e}^{-2}+\,{\tilde{E}}_{a}(\alpha)~{},$ (40)

where $\,N_{\tilde{E}_{\rm tot}}(\alpha)$ is another numerical coefficient
obtained from the simulations. When all the numerical calibration factors are
included we find a criterion that is accurate to within an order of magnitude
to all the density profiles that we explored.

$E_{j,{\rm
iso,tot}}\gtrsim\left[\left(\frac{t_{d}}{t_{e}}\right)^{2}+2\right]\begin{cases}1.5\,{E_{ej,{\rm
tot}}}\left(\frac{\,{\theta_{j,0}}}{0.1{\rm~{}rad}}\right)^{2}~{}&unmagnetized\\\
&\\\ 0.2\,{E_{ej,{\rm
tot}}}\left(\frac{\,{\theta_{j,0}}}{0.1{\rm~{}rad}}\right)^{2}~{}&magnetized\end{cases}~{}.$
(41)

The different factors of magnetized and unmagnetized jets are discussed in §9.
This criterion is applicable only in strong jets with $t_{e}>t_{\rm
col}-t_{d}$ and in weak jets with $t_{e}>t_{a}$. Any weak jet that breaks out
must have $t_{e}>t_{a}$, since at $t<t_{a}$ the jet is unable to make a
significant way within the ejecta. In strong jets the collimation takes place
at $r_{h,col}$ (given in Eq. 18), and therefore the breakout takes place at
$t>t_{\rm col}$ (i.e., after the jet is fully collimated) if $\,{v_{\rm
max}}/\,{v_{\rm min}}>\theta_{j,0}^{-2/3}$, which is expected in BNS mergers
(see discussion at §11). If $\,{v_{\rm max}}/\,{v_{\rm
min}}<\theta_{j,0}^{-2/3}$, then the breakout may take place before a full
collimation is achieved, and then Eq. 41 provides a lower limit for the jet
energy needed for breakout, owing to the faster propagation of collimated jets
compared to uncollimated ones. Note that Eq. 41 is independent of the exact
ejecta density distribution and as we show in the next section it is also
applicable to jets with luminosity that evolves with time, as long as the
total jet energy can be approximated by $E_{j,tot}\sim L_{j}(t_{e})t_{e}$.
Therefore this criterion can be applied to a very wide range of setups.

In case that $t_{e}\gtrsim t_{a}$, as in strong jets with $t_{e}>t_{d}$, and
weak jets with $t_{e}\gtrsim\,{t_{d,\rm eff}}$, the term $(t_{d}/t_{e})^{2}$
in Eq. 41 can be neglected and the criterion can be simplified:

$E_{j,{\rm iso,tot}}\gtrsim\begin{cases}3\,{E_{ej,{\rm
tot}}}\left(\frac{\,{\theta_{j,0}}}{0.1{\rm~{}rad}}\right)^{2}~{}&t_{e}>t_{a},~{}~{}unmagnetized\\\
&\\\ 0.4\,{E_{ej,{\rm
tot}}}\left(\frac{\,{\theta_{j,0}}}{0.1{\rm~{}rad}}\right)^{2}~{}&t_{e}>t_{a},~{}~{}magnetized\end{cases}~{}.$
(42)

If $t_{e}$ is unknown then this simplified criterion is necessary, but it may
be insufficient.

We stress that Eq. 41 and 42 apply to only subrelativistic heads, and given
all the approximations we used, they are accurate to within an order of
magnitude. Thus, if the two sides of Eq. 41 are comparable to within an order
of magnitude, a numerical simulation is probably required to find if the jet
breaks out or not.

###  4.3 Breakout time

In case of a successful breakout, the previous section can be used to estimate
the breakout time. It is important to note that the breakout that we consider
here is the emergence of the jet from the region where the ejecta velocity is
$\,{v_{\rm max}}$, which is the maximal velocity up to which the ejecta
distribution satisfies $\rho\propto v^{-\alpha}$ (we assume $\,{v_{\rm
max}}\ll c$). If there is ejecta at higher velocities than $\,{v_{\rm max}}$,
we assume that the density at these velocities drops faster than $v^{-5}$, so
if the jet breaks out of $\,{v_{\rm max}}$ it is likely to break out of the
faster velocity material, even if the engine stops. We stress this point since
it is likely that in BNS mergers in addition to the bulk of the ejecta there
is also a low-mass fast tail (see Nakar 2020 and references therein). Thus,
while the observations of GW170817, as well as numerical simulations of BNS
mergers, suggest $\,{v_{\rm max}}\approx 0.3$c, the numerical simulations also
find a very steep density profile that continues up to mildly or even ultra-
relativistic velocities. Thus, the breakout time that we find here (from
$\,{v_{\rm max}}$) should not be confused with the breakout time from the fast
tail which can be significantly longer. Yet, the breakout time found here is
approximately the minimal time that the engine should work for the jet to
break out of $\,{v_{\rm max}}$, and thus also most likely from the fast tail.

We start from Eq. 38 which at the time of the breakout can be written as

$\,{\tilde{E}}_{d}\,{\tilde{t}}_{bo}\approx\,N_{\tilde{E}_{\rm
tot}}\,{\tilde{t}}_{bo}^{-2}+\,{\tilde{E}}_{a}~{},$ (43)

where $\,{\tilde{t}}_{bo}\equiv t_{bo}/t_{d}$, and $t_{bo}$ is the duration
that the engine works until the jet breakout (the jet breaks out at time
$t_{d}+t_{bo}$ after the merger). $\,{\tilde{E}}_{d}$ is defined using Eq. 39
by replacing $\,{E_{j,{\rm tot}}}$ with $L_{j}t_{d}$ (or
$\,{\tilde{E}}_{j,{\rm iso,tot}}$ with $\,{L_{j,\rm iso}}t_{d}$):

$\,{\tilde{E}}_{d}\equiv\frac{N_{E}^{5}}{(5-\alpha)^{4-\alpha}}\frac{\,{L_{j,\rm
iso}}~{}t_{d}}{\,{E_{ej,{\rm tot}}}}\theta_{j,0}^{-{2}}~{},$ (44)

Eq. 43 can be solved directly to find $t_{bo}$, or without loss of much
accuracy, the solution to $t_{bo}$ of a subrelativistic jet head can be
approximated as:

$t_{bo}\approx\begin{cases}\,N_{\tilde{E}_{\rm
tot}}^{\frac{1}{3}}\,{\tilde{E}}_{d}^{-\frac{1}{3}}t_{d}=\left(\frac{\,N_{\tilde{E}_{\rm
tot}}(5-\alpha)^{4-\alpha}}{N_{E}^{5}}\frac{\,{E_{ej,{\rm
tot}}}\theta_{j,0}^{2}}{\,{L_{j,\rm
iso}}}t_{d}^{2}\right)^{\frac{1}{3}},&{\scriptstyle\,{\tilde{E}}_{d}>\,{\tilde{E}}_{a}~{}(t_{bo}<t_{d})}\\\
&\\\
\frac{\,{\tilde{E}}_{a}}{\,{\tilde{E}}_{d}}t_{d}=\frac{\,{\tilde{E}}_{a}(5-\alpha)^{4-\alpha}}{N_{E}^{5}}\frac{\,{E_{ej,{\rm
tot}}}}{\,{L_{j,\rm
iso}}}\theta_{j,0}^{2}~{}~{},&{\scriptstyle\,{\tilde{E}}_{d}<\,{\tilde{E}}_{a}~{}(t_{bo}>t_{d})}\end{cases}~{}.$
(45)

Eq. 45 shows two interesting properties of $t_{bo}$. First, in strong jets
with $\,{\tilde{E}}_{d}>\,{\tilde{E}}_{a}$ the dependence of $t_{bo}$ on the
jet and ejecta parameters is relatively weak. Second, in weak jets where
$\,{\tilde{E}}_{d}<\,{\tilde{E}}_{a}$, there is a relatively strong dependency
on the jet parameters (linear or stronger), but no dependency on $t_{d}$. This
implies that for a given ejecta and jet parameters where only $t_{d}$ can
vary, the breakout time has a minimal value which is obtained for
$\,{\tilde{E}}_{d}\lesssim 1$ where $t_{d}<t_{bo}$. The physical reason for
that is that when $t_{d}$, and thus $\,{r_{\rm min}}=\,{v_{\rm min}}t_{d}$, is
small enough, the ejecta density is too high for the jet to propagate through,
and the ejecta has to expand to radii that are much larger than $\,{r_{\rm
min}}$ before the jet head can start propagating. This expansion time provides
an effective delay, which sets the breakout time (see discussion in §3.4.1).

###  4.4 Propagation in aspherical ejecta

In many systems, including binary mergers, the ejecta is expected to be
nonspherical. Our model, which assumes spherically symmetric ejecta, is
applicable to aspherical ejecta (with the minor adjustments listed below), as
long as the properties of the ejecta around the jet axis do not considerably
vary over an angle that is larger than $\,{\theta_{j,0}}$. The reason is that
the entire jet-cocoon structure is confined to an angle $\,{\theta_{j,0}}$
during the jet propagation in the ejecta (§3.5). As a result the jet
propagation is insensitive to the ejecta properties at
$\theta\gg\,{\theta_{j,0}}$. If a nonspherical ejecta satisfy this criterion,
then our analytic model can be used when the ejecta properties are taken as
the isotropic equivalent properties of the ejecta near the poles.
Specifically, one should use the following mapping:

$\begin{array}[]{ccc}\rho_{ej}(v,t)&\to&\rho_{ej}(v,t,\theta=0)\\\
v_{ej}&\to&v_{ej}(\theta=0)\\\
\,{E_{ej}}(v)&\to&\,{E_{ej}}(v,\theta<\,{\theta_{j,0}})\times
4\theta_{j,0}^{-2}\\\ \,{E_{ej,{\rm tot}}}&\to&\,{E_{ej,{\rm
tot}}}(\theta<\,{\theta_{j,0}})\times 4\theta_{j,0}^{-2}\end{array}$ (46)

where $\,{E_{ej}}(v,\theta<\,{\theta_{j,0}})$ and $\,{E_{ej,{\rm
tot}}}(\theta<\,{\theta_{j,0}})$ are, respectively, the kinetic energy at
velocities lower than $v$, and the total kinetic energy which is carried by
material that is confined to $\theta<\,{\theta_{j,0}}$ (one sided).

##  5 Time-dependent luminosity

In the previous sections, and in most previous studies, the jet luminosity is
assumed to be constant. In GRBs the jet luminosity may vary on short time
scales, but time averaged energy output seems to be relatively constant over
the duration of the prompt emission phase. Yet, there are scenarios where the
jet luminosity may vary continuously with time. It is straightforward to
generalize the theory derived above for a luminosity that evolves as a power-
law, as long as most of the jet energy at any given time is injected by recent
(rather than early) engine activity, namely

$L\propto(t-t_{d})^{k}~{}~{}~{};~{}~{}~{}k>-1~{}.$ (47)

We do not give here a complete solution for this case, nor do we test it
numerically. Instead we provide the power-law dependence of $v_{h}$ and an
order of magnitude breakout criterion. Note that if the luminosity drops
faster than Eq. 47 (e.g., $k<-1$) then as long as $\,{\alpha}<5$ (and possibly
for larger values of $\,{\alpha}$) the jet cannot support the propagation of
the head and it is eventually stalled, even if it had been propagating in the
ejecta at first.

For a strong jet, as long as $\eta\gg 1$ we use the static approximation.
Following the derivation of Bromberg et al. (2011), we find that the head
propagation in the Newtonian regime (their equations B2-B11) can be
generalized to accommodate a time evolving luminosity simply by replacing the
constant $L_{j}$ with $L_{j}(t)$. This is accurate up to a constant correction
of the normalization factors that depend on $k$. This corresponds to a head
velocity that evolves as:

$v_{h}\propto(t-t_{d})^{\frac{2+k-\,{\alpha}}{5-\,{\alpha}}}~{}~{}~{};~{}~{}~{}\eta\gg
1~{}.$ (48)

The evolution of the velocity during the asymptotic phase can be derived using
dimensional analysis. Following similar lines as in §3.3 we find (see also
Margalit et al. 2018):

$v_{h}\propto
t^{\frac{6+k-\,{\alpha}}{5-\,{\alpha}}}~{}~{}~{};~{}~{}~{}t>t_{a}$ (49)

and

$\eta_{a}=\frac{1+k}{5-\,{\alpha}}$ (50)

We can now replace the time dependence in Eq. 32 with the one derived in Eq.
48 and following the same line of arguments as in §4, we find that Eq. 36 is
applicable (up to a constant factor that depends on $k$) to a time-dependent
luminosity where we approximate $E_{j}(t)\approx L_{j}(t)(t-t_{d})$.
Similarly, using Eq. 50 we see that Eq. 37 is also applicable to a time-
dependent luminosity. Thus, both the breakout criterion (Eq. 41) and breakout
time given in Eq. 45 are applicable to a time-dependent luminosity (as long as
the $\int_{t_{d}}^{t}Ldt\sim L(t)(t-t_{d})$), where in the definition of
$\,{\tilde{E}}_{d}$ in Eq. 44 the luminosity $\,{L_{j,\rm iso}}$ is taken as
$\,{L_{j,\rm iso}}(t=2t_{d})$.

##  6 Relativistic head

In this paper we consider only Newtonian ejecta. Thus, if the head is
relativistic then $\eta\gg 1$, and the static solution is an excellent
approximation, so we can readily apply the results of Bromberg et al. (2011)
and Harrison et al. (2018). Note that since there is no study of the weakly
magnetized jets in the relativistic regime, our results in this section are
applicable to unmagnetized jets. First we find the condition for the jet head
to be relativistic at the time of the breakout from $\,{v_{\rm max}}$. The
transition from a Newtonian to a relativistic head takes place roughly when
the solution of the Newtonian head (§3) gives $v_{h}=c$. Since the jet head
reaches the location where $\,{v_{ej,h}}=\,{v_{\rm max}}$ after it is fully
collimated, we can estimate the head velocity upon breakout from Eq. 11 where
$\,{r_{\rm min}}$ is replaced with $\,{r_{\rm max}}=\,{v_{\rm max}}t_{d}$.
Doing that and using the definition of $\eta_{0}$ (Eq. 12) we obtain that the
breakout is relativistic if

$\eta_{0}\gtrsim 3\frac{w^{\frac{2-\alpha}{3}}}{\,{\beta_{\rm
min}}}~{}~{}~{};~{}~{}~{}{\rm relativistic~{}breakout}~{},$ (51)

where we ignore the weak dependence of the coefficient in Eq. 11 on $\alpha$
and use the definition

$w\equiv\frac{\,{v_{\rm max}}}{\,{v_{\rm min}}}~{}.$ (52)

The jet head is relativistic when it starts to propagate in the ejecta, even
before it is collimated, if $\,{\eta_{\rm uncol}}>c/\,{v_{\rm min}}$ or:

$\eta_{0}\gtrsim\beta_{\rm
min}^{-1/2}\theta_{j,0}^{-2/3}~{}~{}~{};~{}~{}~{}{\rm
relativistic~{}start}~{}.$ (53)

When both criteria are satisfied the head is relativistic during its entire
propagation through the ejecta. If only one is satisfied, then the jet starts
in one regime and crosses to the other during the propagation in the ejecta.
Below we find the breakout criterion for a jet head that is relativistic at
all times. Usually this criterion is applicable also when only the breakout is
relativistic since most of the propagation time is spent at larger radii
(after collimation). Yet, in this case it is safer to solve the entire
evolution, starting at the collimation phase, using the static approximation.

The location of a jet head when it is relativistic at all times is simply

$r_{h}\approx c(t-t_{d})~{}.$ (54)

Unlike the case of a Newtonian head, a successful breakout is possible also in
cases where the engine stops long before the breakout takes place. The
criterion for a successful breakout of a relativistic head is that the jet
engine time will be long enough to allow for the jet tail (the fluid elements
ejected last) to break out from the ejecta before crossing the reverse shock
at the jet head. Defining $t_{e,bo}$ as the engine work-time at the point that
the fluid element that crosses the reverse shock upon breakout is launched,
the breakout criterion is $t_{e}>t_{e,bo}$. Since the tail Lorentz factor is
much larger than that of the head, the relative velocity between them at any
time is $c/2\Gamma_{h}^{2}(t)$. Assuming that the head Lorentz factor evolves
as $\Gamma_{h}\propto t^{m}$, we obtain

$t_{e,bo}=\frac{R_{bo}}{2c(1-m)\Gamma_{h,bo}^{2}}=\frac{\beta_{\max}}{2(1-m)\Gamma_{h,bo}^{2}}t_{d}~{},$
(55)

where $\Gamma_{h,bo}$ is the head Lorentz factor at the breakout radius,
$R_{bo}\approx\,{v_{\rm max}}t_{d}$.

The dependence of $\,{\Gamma_{h}}$ on the system parameters depends on whether
the jet is collimated or not666Note that we consider breakout that takes place
at $r>\,{r_{h,{\rm col}}}$, so a jet with a relativistic head is uncollimated
only if it is too strong for the cocoon pressure to collimate it (see Bromberg
et al. 2011). This is not to be confused with a jet with a Newtonian head that
may be uncollimated only at very early times, before the pressure in the
cocoon builds up.. Following equation 2 in Harrison et al. (2018) we find that
the criterion for the jet to be collimated upon breakout is

$\eta_{0}\lesssim\frac{w^{\frac{2-\alpha}{3}}}{\,{\beta_{\rm
min}}}\theta_{j,0}^{-10/9}~{}~{}~{};~{}~{}~{}{\rm collimated~{}jet}~{}.$ (56)

In the collimated regime we use equation A11 of Harrison et al. (2018) to
obtain:

$\Gamma_{h,bo}^{2}\approx\left(\frac{1}{10\pi}\frac{(3+\,{\alpha})^{2}}{7-\,{\alpha}}\right)^{\frac{1}{5}}\left(\frac{L_{j}\theta_{j,0}^{-4}}{A_{\rho}c^{5-\,{\alpha}}}\right)^{\frac{1}{5}}(\beta_{\max}t_{d})^{\frac{\,{\alpha}-2}{5}}$
(57)

Using this relation together with Eqs. 55 and 34 we obtain an expression for
relation between $t_{e,bo}$ and $t_{d}$. The requirement that $t_{e}>t_{e,bo}$
gives the following breakout criterion,

$t_{e}\gtrsim 0.7\left(\frac{\,{E_{ej,{\rm tot}}}\theta_{j,0}^{2}}{\,{L_{j,\rm
iso}}}\right)^{\frac{1}{5}}\beta_{\max}^{\frac{2}{5}}~{}~{}t_{d}^{\frac{4}{5}}~{}~{};~{}~{}{\rm
collimated~{}breakout}~{},$ (58)

where we neglect a very weak dependence on $\,{\alpha}$. The same criterion
can be written, using Eq. 55, in terms of the total jet energy,

$E_{j,{\rm
iso,tot}}\gtrsim\frac{3(5-\,{\alpha})}{(3+\,{\alpha})^{2}}\,{E_{ej,{\rm
tot}}}\left(\frac{\,{\theta_{j,0}}}{0.1{\rm~{}rad}}\right)^{2}\left(\frac{\beta_{\max}}{0.3}\right)^{-2}\Gamma_{h,bo}^{8}\\\
~{}~{};~{}~{}{\rm collimated~{}breakout}~{}.$ (59)

In the uncollimated regime equation A19 of Harrison et al. (2018) reads (after
correcting a typo777Note that there is a typo in Eq. A19 of Harrison et al.
(2018) (as well as in Eq. B26 of Bromberg et al. 2011), the dependence on the
opening angle should be $\Gamma_{h}^{2}\propto\theta_{j,0}^{-1}$.)

$\Gamma_{h,bo}^{2}\approx\left(\frac{L_{j}\theta_{j,0}^{-2}}{4\pi
A_{\rho}c^{5-\,{\alpha}}}\right)^{\frac{1}{2}}(\beta_{\max}t_{d})^{\frac{\,{\alpha}-2}{2}}~{}.$
(60)

Following the same steps as in the collimated regime we obtain the breakout
criterion on $t_{e}$

$t_{e}\gtrsim\frac{5(5-\,{\alpha})^{\frac{1}{2}}}{7-\,{\alpha}}\left(\frac{\,{E_{ej,{\rm
tot}}}}{\,{L_{j,\rm
iso}}}\right)^{\frac{1}{2}}\beta_{\max}^{-\frac{1}{2}}~{}~{}t_{d}^{\frac{1}{2}}~{}~{};~{}~{}{\rm
uncollimated~{}breakout}~{},$ (61)

and the equivalent criterion on $E_{j,{\rm iso,tot}}$:

$E_{j,{\rm iso,tot}}\gtrsim
100\frac{(5-\,{\alpha})}{(7-\,{\alpha})}\,{E_{ej,{\rm
tot}}}\left(\frac{\beta_{\max}}{0.3}\right)^{-2}\Gamma_{h,bo}^{2}\\\
~{}~{};~{}~{}{\rm uncollimated~{}breakout}~{}.$ (62)

Note that in the uncollimated regime, the solution is applicable only for
$\alpha<4$ (for $\alpha>4$ the head can propagate without support of the jet).

A comparison of the breakout criteria in the relativistic regimes to the
breakout criterion in the Newtonian head regime (Eq.41) shows that, for a
given ejecta and initial jet opening angle, a breakout of a relativistic head
always requires more energy than the breakout of the Newtonian head. This
implies that the necessary (but possibly insufficient) breakout criterion, Eq.
42, is valid also when a relativistic head is considered.

##  7 Numerical simulations

Model | $\,{\eta_{\rm crit}}$ | $\eta_{0}$ | $\,{\eta_{\rm uncol}}$ | $\,{\eta_{\rm col}}$ | $\theta_{j,0}$ | $\alpha$ | $\,{r_{\rm min}}~{}[z_{\rm beg}]$ | $w_{f}$ | $\tilde{t}_{f}$
---|---|---|---|---|---|---|---|---|---
$\,{\it\alpha_{2}-190}$ | 190 | 240 | 26 | 10 | 0.1 | 2 | 3 | 42 | 2.2
$\,{\it\alpha_{2}n-240}$ | 240 | 450 | 26 | 12 | 0.05 | 2 | 4 | 61 | 2.2
$\,{\it\alpha_{2}w-390}$ | 390 | 290 | 40 | 16 | 0.18 | 2 | 4 | 52 | 2.2
$\,{\it\alpha_{2}-1.7}$ | 1.7 | 11 | 0.4 | 0.4 | 0.1 | 2 | 3 | 2.3 | 16
$\,{\it\alpha_{2}-1}$ | 1.0 | 7.4 | 0.2 | 0.3 | 0.1 | 2 | 3 | 3.6 | 53
$\,{\it\alpha_{2}-0.5}$ | 0.5 | 4.6 | 0.1 | 0.2 | 0.1 | 2 | 3 | 1.8 | 11
$\,{\it\alpha_{2}-{0.15}}$ | 0.15 | 2.1 | 0.04 | 0.08 | 0.1 | 2 | 3 | 1.7 | 72
$\,{\it\alpha_{2}-{0.05}}$ | 0.05 | 1.0 | 0.01 | 0.04 | 0.1 | 2 | 3 | 1.1 | 29
$\,{\it\alpha_{1}-370}$ | 370 | 380 | 60 | 11 | 0.1 | 1 | 8 | 18 | 1.7
$\,{\it\alpha_{1}-1}$ | 1.0 | 7.4 | 0.2 | 0.2 | 0.1 | 1 | 3 | 2.2 | 17
$\,{\it\alpha_{3}-39}$ | 39 | 84 | 4.5 | 4.8 | 0.1 | 3 | 3 | 16 | 2.3
$\,{\it\alpha_{3}-1}$ | 1.0 | 7.4 | 0.2 | 0.4 | 0.1 | 3 | 3 | 2.3 | 9
$\,{\it\alpha_{3}-0.05}$ | 0.05 | 1.0 | 0.01 | 0.05 | 0.1 | 3 | 3 | 1.1 | 27
$\,{\it\alpha_{4}-50}$ | 50 | 100 | 5.8 | 11 | 0.1 | 4 | 3 | 27 | 1.6
$\,{\it\alpha_{4}-1}$ | 1.0 | 7.3 | 0.2 | 0.8 | 0.1 | 4 | 3 | 3.3 | 12
Table 3:  Table of numerical models. Shown are $\,{\eta_{\rm crit}}$ (Eq. 17),
$\eta_{0}$ (Eq. 12), $\,{\eta_{\rm uncol}}$ (Eq. 13) and $\,{\eta_{\rm col}}$
(Eq. 20), as well as the jet opening angle $\theta_{j,0}$, the ejecta density
power-law index $\alpha$, and the minimal radius of the ejecta, $\,{r_{\rm
min}}$ in units of the jet injection height, $z_{\rm beg}$, $w_{f}$ is the
ratio between the ejecta velocity at the jet head location at the end of the
simulation and the minimal ejecta velocity, and $\tilde{t}_{f}\equiv
t_{f}/t_{d}$.

We calibrate the analytic model by performing a set of relativistic
hydrodynamic simulations with pluto (Mignone et al., 2007). The simulations
also test the model within the parameter phase space that they cover. In the
simulations we use a relativistic ideal gas equation of state, both for the
jet and for the ejecta, as appropriate for radiation dominated gas. We use
piecewise parabolic reconstruction method and an HLL Riemann solver. As
pointed out in the previous sections, the general jet evolution is dictated by
three parameters: $\theta_{j,0}$, $\alpha$ and $\,{\eta_{\rm crit}}$ (one can
use $\eta_{0}$ instead of $\,{\eta_{\rm crit}}$), and thus they compose our
parameter space. The full list of simulations and their parameters is given in
Table 3. Once those three parameters are set, the specific choice for the rest
of the parameters is degenerate (e.g., the degeneracy between
$L_{j},A_{\rho},\,{v_{\rm min}},t_{d}$ in $\,{\eta_{\rm crit}}$), and only
affect the scaling of the problem.

In reality the engine has a finite working time, $t_{e}$, and the ejecta have
a maximal velocity, $\,{v_{\rm max}}$, above which the density drops sharply
(sharper than $\rho\propto v^{-5}$). In our simulations the jet engine
operates throughout the entire time of the simulations and the jet never
reaches the edge of the ejecta. Thus, the simulation ends at time $t_{f}$
where the head reached ejecta with velocity $\,{v_{ej,h}}(t_{f})$. This setup
enables us to find out for each simulation whether the jet breaks out
successfully for any value of $t_{e}\leq t_{f}-t_{d}$ and/or $\,{v_{\rm
max}}\leq\,{v_{ej,h}}(t_{f})$. To see that, let us first consider $\,{v_{\rm
max}}$. The evolution of the system at any given time is independent of ejecta
at $v>\,{v_{ej,h}}$, i.e., at any time of the simulation $t$, one can define
$\,{v_{\rm max}}=\,{v_{ej,h}}(t)$ and find the breakout time for this value of
$\,{v_{\rm max}}$. Next, consider an engine that stops working at time
$t<t_{f}$ (so $t_{e}=t-t_{d}$). The jet is choked once the last fluid element
launched by the engine crosses the reverse shock at the jet head. The relative
velocity between the last fluid element and the jet head is $\approx
c-\,{v_{ej,h}}(t)$. Since we simulate only Newtonian heads, $r_{h}(t)\ll
c(t-t_{d})$ and the last fluid element that is launched at time $t$ reaches
the head before it doubles the radius at which it was at time $t$, so the jet
is choked roughly at $r_{h}(t)$. Therefore, for a given engine time $t_{e}$
the jet is choked approximately at $r_{h}(t_{e}+t_{d})$. Thus, our simulations
find the breakout time for all the systems with $\theta_{j,0}$, $\alpha$ and
$\,{\eta_{\rm crit}}$ that appear in Table 3 and with $\,{v_{\rm
max}}/\,{v_{\rm min}}\leq\,{v_{ej,h}}(\,{\tilde{t}}_{f})/\,{v_{\rm min}}$ and
$t_{e}/t_{d}<\,{\tilde{t}}_{f}-1$.

###  7.1 Numerical setup

We recently demonstrated in Gottlieb et al. (2021a) that the propagation
velocity of jets in dense media, such as the ejecta from the merger, is more
accurate in 3D simulations than in 2D simulations. The two main differences
between axisymmetric 2D and 3D models lie in the structure of the jet-cocoon
interface and the jet head. Axisymmetric hydrodynamic 2D jets do not exhibit
the development of local hydrodynamic instabilities along the jet-cocoon
interface, which disrupt the jet’s spine and subsequently slow its head down.
At the same time, those 2D simulations are subject to a numerical artifact of
heavy material that is accumulated on top of the jet head due to the
axisymmetry imposition (known as the “plug”). Before carrying all our
simulations we first tested whether 2D simulations are accurate enough for our
purposes. Unfortunately, we found that axisymmetric simulations yield
significantly different jet propagation velocities than 3D ones, primarily in
jets with low values of $\,{\eta_{\rm crit}}$. The reason is that in weaker
jets the plug plays a much more dominant role. Therefore, all our simulations
are carried out in 3D. We note that all previous works that studied the jet
propagation in BNS ejecta, verified and calibrated their results by comparing
them with axisymmetric 2D simulations (§10).

In our 3D Cartesian grid, the jet is injected axisymmetrically as a top-hat
jet along the $\hat{z}$ axis from the center of the lower boundary. The jet is
injected through a nozzle with a cylindrical radius $r_{j}=10^{7}\,{\rm cm}$,
at the lower boundary at $z=z_{beg}$, where $z_{beg}$ is set by the jet
opening angle such that $z_{\rm beg}=r_{j}\theta_{j,0}^{-1}$. The top-hat jet
is injected hot (initial specific enthalpy, $h=100\gg 1$), with initial
Lorentz factor $\Gamma_{0}$ such that it expands sideways to an angle of
$\theta_{j,0}=1.4\Gamma_{0}^{-1}$ (there is no difference between this mode of
injection and the injection of a conical cold jet with the same
$\,{\theta_{j,0}}$; Mizuta & Ioka 2013; Harrison et al. 2018).

All simulation grids are constituted by three patches on the
$\hat{x}-\hat{y}-$ axes, a central patch of a uniform cell distribution with
two outer logarithmic patches. The $\hat{z}$-axis includes one uniform patch
from $z_{\rm beg}$ to $4\times 10^{8}$ cm, after which the cell distribution
becomes logarithmic. The number of cells and the size of each patch vary
between the cases of $\,{\eta_{\rm crit}}<1$, in which the jet evolution takes
place on small scales, and $\,{\eta_{\rm crit}}>1$, where the jet reaches
farther distances. For $\,{\eta_{\rm crit}}<1$ ($\,{\eta_{\rm crit}}>1$), the
$\hat{x}-\hat{y}-$ axes have one central patch inside the inner $|6\times
10^{7}|$ cm ($|3\times 10^{7}|~{}\theta_{j,0}/0.1$ cm) with a uniform
distribution of 100 ($50~{}\theta_{j,0}/0.1$) cells, and 160 (180) logarithmic
cells up to $|6\times 10^{9}|\,{\rm cm}$. The uniform and logarithmic patches
on the $\hat{z}$-axis have 200 and 400 (500) cells, respectively. The grid
boundary on $\hat{z}$-axis is at $8\times 10^{9}$ cm ($4\times 10^{10}$ cm).
In Appendix B we show that these grids are more than sufficient to reach
convergence.

###  7.2 Numerical results

Coefficient | Equation | $\alpha=1$ | $\alpha=2$ | $\alpha=3$ | $\alpha=4$ | Magnetized
---|---|---|---|---|---|---
$\,{N_{\rm crit}}$ | 17 | 0.5 | $\times 1$
$\,{N_{\rm col}}$ | 20 | 1.0 | 0.9 | 0.5 | $\times 2$
$N_{E}$ | 35 | 0.7 | 0.57 | 0.31 | $\times 1.5$
$\,{\tilde{E}}_{a}$ | 37 | 1.0 | 2.0 | 1.0 | $\times 1$
$\,N_{\tilde{E}_{\rm tot}}$ | 40 | 0.8 | 0.5 | $\times 1$
Table 4: The calibrated numerical coefficients from the simulations. The last
column shows the factor by which one needs to multiply the coefficient to
account for the effect of a non-negligible magnetic field $\sigma\gtrsim
10^{-2}$ (see §9).

Figs. 3, 4, and 5 depict a comparison of the analytic expressions with the
numerical results. The upper panel of each model shows how the numerical value
of $\eta$ (solid blue) compares with analytic approximations (solid red) as
obtained by Eqs. 22 ($\,{\eta_{\rm crit}}>\eta_{a}$, where the analytic fit
starts at $\,{\,{\tilde{t}}_{\rm col}}$) and 28 ($\,{\eta_{\rm
crit}}<\eta_{a}$). The lower panel compares the distance covered by the jet
head, $r_{h}$, as given by Eqs. 25 ($\,{\eta_{\rm crit}}>\eta_{a}$) and 29
($\,{\eta_{\rm crit}}<\eta_{a}$) with its location in the numerical
simulations. The asymptotic curve that all the jets eventually converge to
(Eq. 3.3) is shown as well (dashed red). For very strong jet models (Fig. 3)
we also show $\eta(\,{\tilde{E}})$ (yellow) from Eq. 36, which is the basis
for the breakout criterion (Eq. 40) in those jets. The simulations presented
in those figures were used to find the calibrating factors given in Table 4,
and the analytic curves in the figures include these factors. The comparison
of the simulations to the analytic approximation shows a good agreement. Not
only does the general behavior of the jets agree with the analytic
expectation, but also the numerical coefficients are in good agreement as all
the calibrating factors that we find numerically are of order unity (Table 4).

Figure 3:  A comparison of the analytic expressions with the numerical results
for $\,{\eta_{\rm crit}}\gg 1$ models, with our canonical power-law index,
$\alpha=2$ (top three models), and with different $\alpha$ values (bottom
three models). For each model we show in the upper panel how the numerical
value of $\eta$ (blue) compares with Eq. 22 (red), and with
$\eta(\,{\tilde{E}})$ (yellow, Eq. 36). The lower panel of each model shows
the jet head radius, $r_{h}$, from simulations (blue) and as calculated
analytically (solid red, Eq. 25). We also show $r_{h}$ in the asymptotic phase
(Eq. 3.3, dashed red line) although none of the simulations shown in this
figure reach this phase.

Figure 4:  Same as Fig. 3, but for $\,{\eta_{\rm crit}}\sim 1$ models. Top:
models with different $\,{\eta_{\rm crit}}$ and $\alpha=2$. Bottom: models
with $\,{\eta_{\rm crit}}=1$ and different $\alpha$.  Figure 5:  Same as Figs.
3 and 4, but for $\,{\eta_{\rm crit}}\ll 1$ models. The analytic curves of
$r_{h}$ are obtained from Eq. 29.  Figure 6:  Top: A comparison between
numerical curves of $\,{\tilde{E}}$ (Eq. 35, grey) and the r.h.s. of the
analytic expression in Eq. 40 (blue). The dotted grey lines are weak jet
models $\,{\it\alpha_{2}-{0.15}}$ and $\,{\it\alpha_{2}-{0.05}}$, where
$\,{\eta_{\rm crit}}\ll 1$. The light grey lines represent models with cavity
of uniform density, with solid and dashed lines reflect magnetized and
hydrodynamic jets, respectively (see discussion in §9). Bottom: The same
comparison as in the top panel, but for $\alpha=1$ (red), $\alpha=3$ (yellow)
and $\alpha=4$ (purple). Dashed lines are the numerical curves (dotted is the
weak jet model $\,{\it\alpha_{3}-0.05}$), and solid lines are the analytic
expressions.

Fig. 3 depicts very strong jets with $\,{\eta_{\rm crit}}\gg 1$. The top three
models are with our canonical power-law index of $\alpha=2$ and varying value
of $\,{\theta_{j,0}}$. The three bottom models are with our canonical opening
angle, $\,{\theta_{j,0}}=0.1$ rad and a variety of $\alpha$ values. All models
show a remarkable agreement between the numerical results and both estimates
for $\eta$, with the main deviation taking place with $\eta(\tilde{E})$ as
$\eta$ is approaching unity, where Eq. 36 for $\eta(\tilde{E})\gg 1$ can no
longer be applied. Similarly, all the numerical and the analytic head
locations are consistent with each other to within a factor of order unity.
The simulations shown in Fig. 3 do not reach the asymptotic phase. In all
models $\eta$ shows a universal temporal evolution during which it drops to
$\eta\sim 1$ by $\,{\tilde{t}}=2$ and we expect that after several ejecta
dynamical times $\eta$ converges to its asymptotic value. However, our
computational resources prevent us from continuing most of the strong jet
simulations to later times. We carried out a single strong jet simulation for
longer time, finding that $\eta\approx\eta_{a}$ by $\,{\tilde{t}}=5$ (Fig. 7).
At that time also the location of the jet head is roughly consistent with the
analytic prediction for the asymptotic phase (Eq. 3.3).

Fig. 4 shows marginally strong jets with $\,{\eta_{\rm crit}}\approx 1$. The
top three models are with $\alpha=2$ and $\,{\eta_{\rm crit}}=0.5,1$ and
$1.7$. The three bottom models are all with $\,{\eta_{\rm crit}}=1$ and
different $\alpha$ values. The analytic approximation in these jets is good,
but not as good as in the very strong jet case. This is expected since the
strong jet model is based on the static approximation which is more accurate
for larger values of $\,{\eta_{\rm crit}}$. Nevertheless, even for these
marginal values of $\,{\eta_{\rm crit}}$ the approximation of $r_{h}$ is
better than a factor of $2$ during the entire range of the simulations. Fig. 5
depicts weak jets with $\,{\eta_{\rm crit}}<1$. The agreement between the
analytic and numerical curves of $r_{h}$ is reasonably good (better than about
a factor of 2 at all times).

Fig. 6 tests the accuracy of the calibrated analytic criterion for breakout
(Eq. 40). It compares the analytic estimate of $\,{\tilde{E}}(\,{\tilde{t}})$
(Eq. 38 with calibrated coefficients), which is the basis for the breakout
criterion, with the instantaneous numerical value of $\,{\tilde{E}}$ in each
of our simulations. The numerical value of $\,{\tilde{E}}(\,{\tilde{t}})$ is
extracted by finding $E_{j}$ and $\,{E_{ej}}$ at every time step of the
simulations and plugging them into Eq. 35. In the top panel we present the
cases of $\alpha=2$, with $\tilde{E}$ of all $\alpha=2$ models (grey),
compared with the analytic expression shown in blue. One can see that all
$\,{\eta_{\rm crit}}\gtrsim 1$ models lie roughly on the curve of the analytic
model, including the cases of full cavity (light grey lines) and magnetized
jets (solid light grey, see §9), showing that the breakout criterion is
accurate to all strong jet models. The dotted lines represent numerical models
$\,{\it\alpha_{2}-{0.15}}$, $\,{\it\alpha_{2}-{0.05}}$ in which the jet is
weak. In these models the curve of $\tilde{E}$ begins as the jet first punches
through the ejecta, which corresponds to the sharp drop as $\,{E_{ej}}$
increases from zero. This brief stage is followed by a continuous gradual rise
as the jet is still in the collimation phase, and thus is yet to reach its
asymptotic $\eta_{a}$. Similar to strong jets, the two simulated weak jets
also show a very good agreement with the breakout criterion. The bottom panel
shows the same comparison between analytic and numerical (dashed lines)
results for $\alpha\neq 2$ models. Note that the scatter around the analytic
estimate is larger for larger values of $\alpha$, especially in the simulation
of a weak jet that we carry for $\alpha=3$ (model $\,{\it\alpha_{3}-0.05}$).
Yet, for all our models, the analytic model of $\,{\tilde{E}}$ and the
simulations are in agreement to within an order of magnitude or better. This
agreement implies that our breakout criterion is applicable, at least in the
whole phase-space that we studied numerically.

Finally, we run three additional simulations of strong, marginal and weak jets
(models $\,{\it\alpha_{2}-190}$, $\,{\it\alpha_{2}-1}$ and
$\,{\it\alpha_{2}-{0.15}}$, respectively) for a longer time. Due to
computational limitations these simulations were carried out in half of the
resolution compared to the rest of the simulations. We verify that this
resolution is high enough by a comparison to the full resolution simulations
of the same models. The goal of these simulations is to see the transition to
the asymptotic phase and the convergence of $\eta$ to its asymptotic value,
which for $\,{\alpha}=2$ is predicted to be $\eta_{a}=1/3$ (Eq. 9). In Fig. 7
we present the behavior of $\eta$ (top) and $r_{h}$ (bottom) of the three long
duration simulations. For the very strong jet (blue) the simulation is carried
out up to $\,{\tilde{t}}=5$ at which time $\eta\approx\eta_{a}$. The marginal
(red) and weak (yellow) jets are simulated up to $\,{\tilde{t}}\approx 50$ and
$\,{\tilde{t}}\approx 70$, respectively. The high resolution simulations are
depicted in cyan, magenta and green, showing a good agreement at all times
with the low resolution simulations. In all marginal jets we see that within
several dynamical times the value of $\eta$ drops at first slightly below
$\eta_{a}$ (see Fig. 4). In Fig. 7 we see that following this drop, $\eta$ of
model $\,{\it\alpha_{2}-1}$ rises slowly until it converges to $\eta_{a}$. The
effect of these small variations of $\eta$ (where $\eta$ never drops far below
$\eta_{a}$) on the analytic model of $r_{h}$ is not significant, as the jet
head location is in a good agreement with the analytic estimate at all times.
The weak jet ($\,{\it\alpha_{2}-{0.15}}$) features a slow monotonic increase,
up to fluctuations in the jet head velocity, reaching $\eta_{a}$ over a longer
timescale of $\,{\tilde{t}}\gtrsim 70$.

Figure 7: Long duration simulations of very strong jet
($\,{\it\alpha_{2}-190}$; blue), marginally strong jet ($\,{\it\alpha_{2}-1}$;
red) and weak jet ($\,{\it\alpha_{2}-{0.15}}$; yellow) that extend to the
asymptotic phase. Solid lines mark the numerical results and dashed lines mark
the analytic model (Eq. 25). High resolution $\,{\it\alpha_{2}-190}$ (cyan),
$\,{\it\alpha_{2}-1}$ (magenta) and $\,{\it\alpha_{2}-{0.15}}$ (green) models,
which are available only for a limited duration, are also shown in the bottom
panel for comparison.

##  8 Non-empty cavity

Figure 8:  Comparison of model $\,{\it\alpha_{2}-190}$ with two different
types of cavity. Shown are the original model with an empty cavity (blue) and
the same model with a cavity of uniform density that is equal to
$\rho(\,{r_{\rm min}},t_{d})$ (red). The figure also shows a model with a
similar jet and a larger cavity ($\,{r_{\rm min}}=3\times 10^{9}$ cm compared
to $\,{r_{\rm min}}=3\times 10^{8}$ cm in $\,{\it\alpha_{2}-190}$, yellow). At
$t_{d}$ the density at $r>3\times 10^{9}$ cm is similar to the two other
models, while at smaller radii $\rho(r<3\times 10^{9}~{}\,{\rm
cm},t_{d})=\rho(3\times 10^{9}~{}\,{\rm cm},t_{d})$.

In our analytic solution we assumed that at $t=t_{d}$ there is an empty cavity
at $r<\,{r_{\rm min}}$. However, in a realistic scenario it is possible that
the density in the cavity is not negligible. As discussed in §2.4.1, it is
possible, for example, that there is ejecta with $v<\,{v_{\rm min}}$ with
density that is low compared to extrapolation of $\rho\propto v^{-\alpha}$.
Another possibility is that there is a continuous mass ejection so mass that
was ejected at $t\sim t_{d}$ fills the cavity. In this section we consider the
effect of a non-empty cavity on the propagation of the jet.

For a strong jet a non-negligible density at $r<\,{r_{\rm min}}$ is expected
to have two opposite effects on the head velocity. The obvious one is that the
density can be high enough to slow down the head velocity to be
subrelativistic in the cavity, thereby delaying the time that the head reaches
$\,{r_{\rm min}}$ compared to the empty cavity case. The second is that the
higher density reduces the volume that the shocked jet can occupy, thereby
increasing the pressure that the cocoon applies on the jet. The result is that
the jet is at least partially collimated by the time that it reaches
$\,{r_{\rm min}}$, and therefore its velocity at this point is faster than
that of a jet that propagates in an empty cavity. Thus, when comparing the
head propagation in empty and non-empty cavities we expect the former to
propagate faster at first, but the latter should catch up at $r>\,{r_{\rm
min}}$. The propagation in a non-empty cavity is bounded from above by the
solution of empty cavity (Eq. 25), and it can be bounded from below by a jet
in a full cavity where at $t_{d}$ the density in the cavity is
$\rho(r<\,{r_{\rm min}})=\rho(\,{r_{\rm min}})(r/\,{r_{\rm
min}})^{-\,{\alpha}}$. Our analytic solution already includes this scenario,
since it is equivalent to taking Eq. 25 in the limit $\,{v_{\rm min}}\to 0$:

$\displaystyle r_{h}=$
$\displaystyle\left(\frac{2(2\pi(5-\,{\alpha}))^{\frac{1}{3}}\,{N_{\rm
col}}}{15}\right)^{\frac{3}{5-\,{\alpha}}}\left(\frac{L_{j}}{\,{E_{ej,{\rm
tot}}}}\theta_{j,0}^{-4}\right)^{\frac{1}{5-\,{\alpha}}}\,{v_{\rm max}}$
$\displaystyle\times~{}(t-t_{d})^{\frac{3}{5-\,{\alpha}}}t^{\frac{3-\,{\alpha}}{5-\,{\alpha}}}~{}~{}~{};~{}~{}~{}{\rm
strong~{}jet~{}in~{}full~{}cavity},$ (63)

The two solutions (in empty and in full cavities) merge at
$\,{\,{\tilde{t}}_{\rm col}}$ of the empty cavity jet (since at
$\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm col}}$ both solutions follow the
approximation of a fully collimated jet is a static medium), implying that at
$\,{\tilde{t}}>\,{\,{\tilde{t}}_{\rm col}}$ the jet location is independent of
the density in the cavity.

The above expectation is confirmed by numerical simulations888Here we are
interested in comparing the effect of different cavities rather than the
absolute behavior and calibration of the jet head motion. Thus, the
simulations in this section are conducted in 2D, where $r_{h}$ is smaller by
$\sim 1.5$ compared to their 3D counterparts.. Fig. 8 shows $\eta$ and $r_{h}$
of two simulations with similar setups, where the only difference is that in
one the cavity is empty whereas in the other it is filled with a uniform
density $\rho(r<\,{r_{\rm min}},t_{d})=\rho(\,{r_{\rm min}},t_{d})$. As
expected, at first $r_{h}$ of the head in the empty cavity is larger, but
quickly the velocity of the head in the filled cavity becomes larger, until at
$\,{\tilde{t}}\approx 1.15$ it catches up, and from that time on the two heads
propagate with similar $r_{h}(t)$. The jet in the empty cavity becomes fully
collimated roughly at $\,{\,{\tilde{t}}_{\rm col}}=1.1$.

Fig. 8 also shows a third model with a uniform density inside the cavity,
$\rho(r<\,{r_{\rm min}},t_{d})=\rho(\,{r_{\rm min}},t_{d})$. It has the same
jet parameters, $t_{d}$ and density normalization $A_{\rho}$ as the first two
simulations, but its $\,{r_{\rm min}}=3\times 10^{9}$ cm rather than
$\,{r_{\rm min}}=3\times 10^{8}$ cm of the first two (equivalent to $\,{v_{\rm
min}}$ that is larger by a factor of 10). The head in the third model
propagates at first at a faster velocity than the other two, up to a radius
larger than $10^{9}$ cm. But, as soon as it starts propagating inside the bulk
of the ejecta, it significantly decelerates such that the three models
converge at $\,{\tilde{t}}\sim 1.3$. At later times all the models have
similar $r_{h}(t)$. Our conclusion is, that as expected, the density in the
cavity affects the head location of strong jets up to a few times $\,{r_{\rm
min}}$ and it has only a small effect on $r_{h}(t)$ once the head propagates a
significant distance in the bulk of the ejecta (at $r>\,{r_{\rm min}}$).
Naturally, the density in the cavity should have no effect on the propagation
during the asymptotic phase.

For a weak jet a non-empty cavity may prevent the head from propagating in the
cavity. Nevertheless, once the ejecta expands, the jet starts propagating in
the cavity at $t\ll t_{a}$. Since at $t\gg t_{a}$ the jet propagates in the
asymptotic phase where the initial conditions (e.g., $t_{d}$ and $\,{r_{\rm
min}}$) are forgotten, it is expected that by that time the head of the jet
with the non-empty cavity approaches the same $r_{h}(t)$ as that of the jet
with the empty cavity. Fig. 9 shows (in addition to the magnetic jet) two
unmagnetized jets. One with an empty cavity and $\,{\eta_{\rm crit}}=1$ and
the other with a constant density cavity. The latter propagates at first very
slowly inside the dense cavity, but within several dynamical times it
approaches the head location of the former.

We conclude that in both strong and weak jets, the jet propagation is affected
by the cavity density as long as the jet head is still at radii that are
comparable to the cavity size. As the jet head enters the ejecta, the strong
(weak) jet head is approaching Eq. 25 (29), regardless of the cavity density,
and when $r_{h}\gg\,{v_{\rm min}}t$ it forgets the initial conditions in the
cavity entirely.

Figure 9: $\eta$ (top) and $r_{h}$ (bottom) plots of numerical (solid lines)
and analytic (dashed lines) models of magnetized (blue) and hydrodynamic (red)
jets in model $\,{\it\alpha_{2}-1}$ with cavity of uniform density, compared
with the canonical unmagnetized $\,{\it\alpha_{2}-1}$ model (yellow) with an
empty cavity. Magnetized jets with $\sigma=0.1$ and $\sigma=0.01$ yield the
same results, and thus the blue curve represents the entire range in between.

##  9 Effect of magnetization

An important jet property, which has never been addressed in the context of
analytic description of jet propagation in expanding medium, is the jet’s
magnetization, $\sigma\equiv\frac{B^{\prime 2}}{4\pi h\rho c^{2}}$, where
$B^{\prime}$ is the proper magnetic field. Even a subdominant999Studies of
jets in stars found that even if the jets are launched as Poynting-flux
dominated, magnetic dissipation at the collimation nozzle may reduce their
magnetization to subdominant above the nozzle (Bromberg & Tchekhovskoy, 2016;
Gottlieb et al., 2022c). magnetization can greatly alter the jet propagation
in the medium, e.g., by a stabilizing the jet against the development of
hydrodynamic instabilities on the jet-cocoon boundary (Gottlieb et al., 2020b,
2021b; Matsumoto et al., 2021). In a recent study, Gottlieb et al. (2020b)
performed numerical simulations to compare hydrodynamic and magnetized jets
with toroidal magnetic fields. They found that for a hydrodynamic jet which
propagates at $v\approx 0.2$c in a star, introducing an initial magnetization
of $\sigma\sim 10^{-2}-10^{-1}$ to the jet increases its velocity by a factor
of three. This result, albeit obtained in a static medium and for a specific
set of jet and medium parameters, can be utilized for our results in the
regime of strong jets with $\,{\eta_{\rm crit}}\gg 1$.

Here we carry out two simulations of a magnetized jet in the regime of
$\,{\eta_{\rm crit}}\sim 1$, to complement the work of Gottlieb et al. (2020b)
on propagation in static media, which is applicable to strong jets. We use a
similar magnetic configuration of a toroidal field with initial
$\sigma=10^{-1}$ and $\sigma=10^{-2}$ which we apply to model
$\,{\it\alpha_{2}-1}$101010For the magnetic jet simulations we make use of the
constrained transport method in pluto.. To avoid numerical instability, we
fill the cavity with uniform density such that $\rho(r<\,{r_{\rm
min}},t_{d})=\rho(\,{r_{\rm min}},t_{d})$. We also conduct a similar filled
cavity simulation with a hydrodynamic jet for comparison (see §8 for a
discussion about the effect of dense cavity on the jet propagation).

In Fig. 9 we compare the numerical (solid lines) and analytic (dashed lines)
results of magnetized (blue) and hydrodynamic (red) jets with filled cavity,
and a hydrodynamic jet with an empty cavity (yellow). The analytic models are
derived assuming an empty cavity, and therefore they are expected to converge
with the filled cavity simulation only when $r_{h}$ reaches a few times
$\,{v_{\rm min}}t$. We find that the magnetized jets with $\sigma=0.1$ and
$\sigma=0.01$ feature the same jet evolution, similar to what was found in
Gottlieb et al. (2020b), and thus the blue curve represents both values of
$\sigma$. The effect of the magnetic stabilization of the instabilities along
the jet boundary is that the jet head velocity is $\sim 2-3$ faster than that
of hydrodynamic jets. This implies, in turn, that the jet head location is
also about 2-3 times farther than its hydrodynamic counterpart. Therefore, the
jet crosses the filled cavity and becomes collimated faster. We find that the
required correction factor is doubling the value of $\,{N_{\rm col}}$. As can
be seen in the bottom panel, taking $\,{N_{\rm col}}$ to be twice larger sets
the analytic jet head location to fit the simulations well. The difference in
the head velocity that we find here (a factor of 2-3) is similar to the one
found in the case of a static medium (Gottlieb et al., 2020b), which is
applicable for strong jets as long as $\eta\gg 1$. We conclude that both in
the regime of $\,{\eta_{\rm crit}}\sim 1$ and $\,{\eta_{\rm crit}}\gg 1$ the
jet head velocity propagates 2-3 times faster, and taking $\,{N_{\rm col}}$ to
be twice larger than in hydrodynamic jets is expected to provide a reasonable
approximation to magnetized jets with $\sigma\sim 10^{-2}-10^{-1}$.

By the requirement that the jet head velocity is $\sim$ 2 faster than its
hydrodynamic counterpart, Eqs. 32 and 35 dictate that $N_{E}$ is 1.5 times
larger than the value obtained for hydrodynamic jets. By applying the
correction factor to $N_{E}$, the magnetized jet model is found to coincide
with the universal asymptotic behavior of $\tilde{E}$ in Fig. 6 (solid light
grey line in the top panel). This follows a short delayed phase of
$\tilde{E}\gg 1$ which comes as a result of the filled cavity which delays the
onset of the jet propagation in the ejecta. For comparison, we also present in
Fig. 6 a hydrodynamic jet with a filled cavity which exhibits the same delayed
behavior (dashed light grey line in the top panel). Finally, Eq. 44 shows that
the increase in $N_{E}$ for mildly-magnetized jets is equivalent to an
increase of almost an order of magnitude in the jet luminosity. Similarly, Eq.
45 implies that the presence of a non-negligible magnetic field shortens the
breakout time by a factor of two if the jet breaks out at $t_{bo}<t_{d}$, and
by almost an order of magnitude if $t_{bo}>t_{d}$.

##  10 Comparison with previous studies

The propagation of relativistic jets in expanding media was studied by many
authors over the past few years (e.g., Murguia-Berthier et al., 2014; Murguia-
Berthier et al., 2017, 2021; Duffell et al., 2018; Margalit et al., 2018;
Matsumoto & Kimura, 2018; Gill et al., 2019; Lazzati & Perna, 2019; Nakar,
2020; Beniamini et al., 2020; Lyutikov, 2020; Salafia et al., 2020; Hamidani
et al., 2020; Hamidani & Ioka, 2021). Here we summarize the findings of
previous studies and compare them to our results of the jet head propagation,
the breakout time and the breakout criterion. We focus on studies that derived
analytic expressions that can be compared to ours.

The most prominent difference between our work and previous studies is that
our model was derived (and tested and calibrated numerically) for all types of
jets during the entire evolution (weak and strong, unmagnetized and
magnetized, Newtonian and Relativistic). That is, we formulate the jet
propagation for both strong and weak jets and we do that i) during the pre-
collimation and the collimation phases; ii) during the collimated phase of a
strong jet when $v_{h}\gg\,{v_{ej,h}}$; and iii) during the asymptotic phases
as $v_{h}-\,{v_{ej,h}}<\,{v_{ej,h}}$. All previous analytic works considered
only a single regime, either a strong jet with a Newtonian head and
$v_{h}\gg\,{v_{ej,h}}$ or a Newtonian head in the asymptotic phase. As we
discuss below, even in these phases there are some differences between our
results and previous studies.

First attempts for deriving a breakout criterion for a jet in the Newtonian
regime in expanding media were provided by Murguia-Berthier et al. (2014);
Murguia-Berthier et al. (2017); Duffell et al. (2018). These studies (and also
Lyutikov 2020 later on) ignored the jet collimation. Namely, they assume
(explicitly or implicitly) that the jet propagates in the ejecta with an
opening angle that is similar to the one it was injected with, thus missing
the factor of $\,{\theta_{j,0}}^{2}$ in the breakout criterion (e.g., in Eq.
42). This is inconsistent with any of the jet regimes, since when the head is
Newtonian the pressure, which builds up in the cocoon during the propagation,
is always strong enough to collimate the jet. Simulations show that the only
case where a jet with a Newtonian head remains uncollimated is when the
initial jet opening angle is large, $\,{\theta_{j,0}}\gtrsim 0.5$ rad (e.g.,
Gottlieb et al., 2018b). In such jets the shocked ejecta is accumulated on top
of the jet head since it does not have enough time to be pushed aside. The
evolution of such jets is more similar to that of a segment from a spherical
explosion,where Eq. 6, which is used by all these studies, is not applicable.
Therefore, the various results derived in studies that neglect the collimation
by the cocoon of a jet with a Newtonian head, would not be applicable to a
broad range of regimes of physical jets.

All other studies carried out some type of generalization of the model derived
by Bromberg et al. (2011). Some used the simplification of a single typical
jet head velocity (e.g., Gill et al., 2019; Lazzati & Perna, 2019; Nakar,
2020), which is similar to taking $\alpha=2$. Others, considered a range of
$\alpha$ values (e.g., Margalit et al., 2018; Matsumoto & Kimura, 2018;
Hamidani et al., 2020; Hamidani & Ioka, 2021). Nakar (2020) considered only
the strong jet regime where $\eta>1$ obtaining analytic results that are
consistent with ours in jets with $\,{\eta_{\rm crit}}\gg 1$ at $t\lesssim
2t_{d}$. Margalit et al. (2018) and Matsumoto & Kimura (2018) provided an
analytic expression for the evolution during the asymptotic phase, when $t\gg
t_{d}$. The power-law of the temporal evolution that they found (equations 22
of Margalit et al. 2018 and A10 of Matsumoto & Kimura 2018) is the same as the
one that we derive (Eq. 8), but their normalizations of the asymptotic phase,
which were not tested numerically, are different (e.g., in Matsumoto & Kimura
2018 it is about an order of magnitude larger than in Eq. 3.3 for
$\,{\alpha}=2$). They also derived a breakout criterion for the asymptotic
phase that has the same functional form as our Eq. 42. However, also here
their normalizations are significantly different. Their breakout criterion are
lower than ours by about 2 orders of magnitude, which is a result of their
overestimate of the jet head location. Hamidani & Ioka (2021) carried out a
careful analytic study with a comparison to 2D numerical simulations. Their
analytic expression to the jet head in expanding media (their Eq. 42) includes
some dependence on $t_{d}$ although it is derived assuming that $t\gg t_{d}$.
At late times the temporal dependence of their expression is $r_{h}\propto
t^{\frac{6-\,{\alpha}}{5-\,{\alpha}}}\ln(\,{\tilde{t}})^{-\frac{2}{5-\,{\alpha}}}$,
with a very weak dependence on $t_{d}$. This is compared to our expression
$r_{h}\propto t^{\frac{6-\,{\alpha}}{5-\,{\alpha}}}$. Comparing the predicted
value of $r_{h}$ for typical jet and ejecta parameters we find that the
difference between the predictions depends mostly on $\,{\alpha}$ and on
whether the jet is strong or weak, where in some of the regimes their equation
42 predicts $r_{h}$ that agree with ours to within an order of magnitude while
in others their $r_{h}$ is larger by more than an order of magnitude. Finally,
for $\,{\alpha}\geq 3$ their expression is undefined and when $\,{\alpha}$
approaches 3 from below, their $r_{h}$ diverges and so does the ratio between
their solution and ours. In Fig. 10 we demonstrate the ratio between our
solution to the jet head evolution in unmagnetized strong and weak jet models
with different values of $\,{\alpha}\leq 2.5$, and those of Hamidani & Ioka
(2021) (HI21) and Matsumoto & Kimura (2018) (MK18). These are the only works
that derived an analytic solution to $r_{h}(t)$. It shows that at early times
($\,{\tilde{t}}\lesssim 1$) the solutions of HI21 and MK18 often deviate from
our solution by $\sim 2$ orders of magnitude. This is expected since their
models were designed to fit only the asymptotic phase. Indeed at later times,
($\,{\tilde{t}}\gg 1$), the difference between the models is lower. The
solution of HI21 deviates from ours by $\sim 0.5-1$ order of magnitude during
the asymptotic phase of most of the presented scenarios, while MK18 model
yields a similar ratio in weak jets, but features a larger ratio in strong
jets, with a jet head location that is between 1 to 2 orders of magnitude
farther than ours.

Figure 10: The ratio between the analytic solution to $r_{h}$ of Hamidani &
Ioka (2021) (HI21; blue) and Matsumoto & Kimura (2018) (MK18; red), for strong
(top) and weak (bottom) jets, and our analytic solution for unmagnetized jets.
Each set of models is shown for $\alpha=1$ (dashed), $\alpha=2$ (solid) and
$\alpha=2.5$ (dash-dotted). The model parameters that are used for the
analytic solutions in the strong (weak) jets with $\alpha=1,2,2.5$ are the
same as those of $\,{\it\alpha_{1}-370}$ ($\,{\it\alpha_{1}-1}$),
$\,{\it\alpha_{2}-190}$ ($\,{\it\alpha_{2}-{0.05}}$), and
$\,{\it\alpha_{3}-39}$ ($\,{\it\alpha_{3}-0.05}$), respectively. Note that
since the models of HI21 and MK18 are undefined at $\alpha=3$, we use here the
parameters of models $\,{\it\alpha_{3}-39}$ and $\,{\it\alpha_{3}-0.05}$ with
the only modification of taking $\alpha=2.5$. For the solution of HI21 we use
their table 1 to set the following values for their parameters
$\langle\eta^{\prime}\rangle=1/4$ and $\langle\chi\rangle=5/4$, as those are
obtained in their model with an opening angle similar to ours,
$\,{\theta_{j,0}}=0.1$ rad.

The results of Duffell et al. (2018) deserve some additional discussion. They
report that for successful jets they find two regimes of jet breakout from the
ejecta: early breakout on timescales shorter than the engine duration, and
late breakout well after the engine shuts off. This seems to be in
disagreement with the result that any jet fluid element that crosses the
reverse shock at the head and spills into the cocoon is mixed with material
from medium and cannot be part of any ultra-relativistic outflow. This result,
which was found in a large number of simulations by different groups and
different codes (e.g., Mizuta & Aloy, 2009; Mizuta & Ioka, 2013; Harrison et
al., 2018; Gottlieb et al., 2018a), implies that when the jet head is
Newtonian it is choked quickly after the engine stops (more precisely a time
$r_{h}/c$ after the engine stops). A careful look at the setup used by Duffell
et al. (2018) finds that they performed simulations of jets that are launched
with a constant luminosity for a characteristic time, $T$, where at $t\gg T$
the jet has a long lasting tail with a luminosity that drops as $t^{-2}$.
Namely, their engine does not really shut off at any time. Thus, if their
simulations are consistent with those of other groups, then in their
simulations in which $t_{bo}\gg T$, it is the tail of the ejecta that is
launched at $t\gtrsim t_{bo}$ that breaks out of the ejecta. Here it is
important to note that the density distribution considered by Duffell et al.
(2018) has two characteristic regimes. Their ejecta is homologous and it has a
leading edge which at the beginning of the simulation is located at $R_{0}$.
The initial distribution at $r\ll R_{0}$ is approximately a power-law in $r$
with $\alpha=2.5$. Near the leading edge, the distribution drops sharply as a
power-law of the distance from the edge (instead of a power-law in $r$),
$\rho_{ej}\propto(1-r/R_{0})^{3/4}$. Our results show that a jet can propagate
in a medium with $\alpha<5$ only if the head is supported by a continuous flow
of jet material through the reverse shock with a luminosity that drops slower
than $t^{-1}$ (see also111111As was recently shown by Govreen-Segal et al.
(2021) a spherical Newtonian shock in expanding media always dies out when
$\,{\alpha}<8$, unless there is a continuous energy injection into the shock.
Since jets expand sideways, as the energy injection stops they cannot be
sustained without a continuous energy injection even at $\,{\alpha}$ values
that are larger than 8. Govreen-Segal et al. 2021). The combination of their
density profile and jet luminosity evolution leads to a head that can
propagate at $t\lesssim T$ through the inner parts of the ejecta, and then be
supported by a jet tail where $L_{j}(j)\propto t^{-2}$ as it propagates near
the edge of the ejecta, where the density gradient is very sharp. This implies
that the ’late breakout’ is actually a breakout of the tail of the jet and it
is not a generic property of jets in expanding ejecta in the sense that it is
a result that depends of the specific setup chosen by Duffell et al. (2018).
Therefore, the various relations that they find in their simulations that
refer to the ‘late breakout’ are relevant only to their setup. For example, if
the setup is changed so the jet launching is truly terminated at $T$ (i.e.,
$L_{j}(t>T)=0$), then there is no late jet breakout.

This result raises a more general interesting question. Is it possible that
many of the observed GRBs are generated by jets that break out while their
luminosity decays strongly with time (faster than $t^{-1}$) after propagating
in sharp density media ($\frac{d\ln\rho}{d\ln v}\ll-5$), for which the
breakout time that we derive here is incorrect? Current observations suggest
that typical sGRBs are inconsistent with this picture. There are hundreds of
sGRBs for which we have detailed observations of the prompt gamma-ray
emission. The temporal properties of this emission were studied by many
authors and it is found to be typically variable, with a complex temporal
structure which is composed of many pulses, some of them overlapping and some
separated (e.g., McBreen et al., 2001; Nakar & Piran, 2002; Bhat et al., 2012;
Li et al., 2020). Interestingly, these studies did not find an evolution of