the pulses properties with time. Instead, the prompt emission seems to last
for a duration, $T$, during which the pulses’ properties such as luminosity
and duration do not show any specific evolution. This continues until at some
point the gamma-ray emission stops over a time-scale that is comparable to the
duration of a single pulse, which is often much shorter than the duration of
the burst. Since the prompt emission reflects the activity of the central
engine after the breakout, this result suggests that the central engine does
not change its properties considerably for the duration of the burst. The
gamma-ray emission stops either because the engine shuts off abruptly or the
efficient generation of gamma-rays stops for some unknown reason. Otherwise,
we would have seen the marks of the evolution of the jet properties on the
prompt emission. For example, if the jet luminosity would have been decaying
monotonically and significantly over the duration of the burst, we would
expect that on average the pulses seen at the beginning of the burst would be
brighter than those seen at the end of the burst. When combined with the fact
that the engine is causally disconnected from the head and is unaware of the
jet breakout, the observational properties of the prompt emission support a
picture where the jet luminosity is roughly constant (up to possibly a short
duration variability) during the entire propagation in the ejecta, through the
breakout and the production of the gamma-rays. Thus, these observations are in
tension with the late breakout scenario suggested by Duffell et al. (2018).

Duffell et al. (2018) also found a breakout criterion that they calibrate
numerically. The functional form of their criterion for breakout of jets with
$t_{bo}<T$ (their equations 20 and 21) differs from ours (and from that of
Margalit et al. 2018 and Matsumoto & Kimura 2018). Their criterion is
equivalent to our Eq. 42 without the dependence on the r.h.s. on
$\theta_{j,0}^{2}$. The reason that they did not find this difference in their
numerical tests is likely a result of the narrow range of $\,{\theta_{j,0}}$
(factor of 2) with which they tested the criterion and the effect of the late
breakout on the test results. When we set the jet opening angle in our
breakout criterion to the range in which they used to calibrate the numerical
coefficient $\kappa$ of their breakout criterion, $\,{\theta_{j,0}}=0.07-0.14$
rad, their results are similar to ours to within an order of magnitude. In
fact, we find that when our Eq. 42 is applied to their simulations it predicts
that the transition from successful to failed jets should take place right at
the phase-space region that Duffell et al. (2018) found a transition from
’early’ to ’late’ breakouts. This is consistent with our understanding that
late breakouts corresponds to jets that almost broke free at $t\approx T$ and
that are sustained longer by the support of the jet tail and/or the
propagation in the very steep density gradient near the leading edge of the
ejecta.

Finally, Beniamini et al. (2020) studied the impact of the propagation time on
observations of sGRBs. They use it to obtain constraints on $t_{d}$. They
study two regimes, static and dynamic, corresponding to $t_{bo}<t_{d}$ and
$t_{bo}>t_{d}$, respectively. For the static regime they applied a similar
analysis to that of Murguia-Berthier et al. (2017); Gill et al. (2019) and for
the dynamic regime they adopted the solution by Duffell et al. (2018). In §11
we discuss their results in view of our findings.

##  11 Application to BNS mergers and short GRBs

The results of the jet propagation have some important implications for BNS
mergers and short GRBs. For this discussion we summarize first the main
propagation parameters normalized to those inferred by observations of
GW170817 (see Nakar 2020 and references therein), $\,{E_{j,{\rm tot}}}\approx
10^{50}$ erg, $\,{\theta_{j,0}}\approx 0.05$ rad, $\,{E_{ej,{\rm tot}}}\approx
10^{51}$ erg, $\,{v_{\rm min}}\approx 0.05$c and121212Note that $\,{v_{\rm
max}}$ is the maximal velocity of the bulk of the ejecta. Namely, it is the
velocity where most of the kinetic energy of the ejecta is. It is believed
that in GW170817 there was ejecta with velocities higher than $\sim 0.3$ c,
but it carried a negligible fraction of the energy and did not have a
significant effect on the jet propagation. See discussion in §4.3. $\,{v_{\rm
max}}\approx 0.3$c, so $w\approx 6$. While the duration of the engine,
$t_{e}$, and thus the jet isotropic equivalent luminosity, $\,{L_{j,\rm
iso}}$, are not constrained, we take a canonical value of $\,{L_{j,\rm
iso}}=10^{52}~{}{\rm erg~{}s^{-1}}$. We use $\alpha=3$ since the amount of
mass in low ($v\approx 0.1c$) and high ($v\approx 0.3c$) velocities in
GW170817 was found to be comparable, and under the assumption of isotropic
mass distribution, one obtains $\alpha=3$. There are also numerical studies
that find that $\alpha\approx 3-4$ provide a reasonable description of the
ejecta density profile (e.g., Nagakura et al., 2014). When the parameters of
unmagnetized and weakly magnetized jets are different we give both where the
weakly magnetized values are in ‘[]’.

The dimensionless parameters for BNS mergers are:

$\eta_{0}=80\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{-\frac{1}{3}}\left(\frac{w}{6}\right)^{\frac{2}{3}}\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{\frac{1}{3}}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{-\frac{2}{3}}\left(\frac{t_{d}}{1{\rm~{}s}}\right)^{\frac{1}{3}}~{},$
(64) $\,{\eta_{\rm crit}}=17\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{-\frac{1}{2}}\left(\frac{w}{6}\right)\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{\frac{1}{2}}\left(\frac{t_{d}}{1{\rm~{}s}}\right)^{\frac{1}{2}}~{},$
(65) $\,{\eta_{\rm col}}=4[8]\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{-\frac{1}{3}}\left(\frac{w}{6}\right)^{\frac{2}{3}}\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{\frac{1}{3}}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{-\frac{2}{9}}\left(\frac{t_{d}}{1{\rm~{}s}}\right)^{\frac{1}{3}}~{},$
(66) $\,{\tilde{E}}=8[60]\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{-1}\left(\frac{\,{E_{j,\rm
iso}}}{10^{52}{\rm~{}erg}}\right)\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{-2}~{},$
(67)

and

$\,{\tilde{E}}_{d}=8[60]\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{-1}\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{-2}\left(\frac{t_{d}}{1{\rm~{}s}}\right)~{}.$
(68)

The corresponding breakout time in the Newtonian regime is:

$t_{bo}\approx\begin{cases}0.4[0.2]~{}{\rm s}~{}\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm erg}}\right)^{\frac{1}{3}}\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm
erg~{}s^{-1}}}\right)^{-\frac{1}{3}}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm
rad}}\right)^{\frac{2}{3}}\left(\frac{t_{d}}{1{\rm
s}}\right)^{\frac{2}{3}}&{\scriptstyle\,{\tilde{E}}_{d}>1~{}(t_{bo}<t_{d})}\\\
&\\\ 0.1[0.02]~{}{\rm s}~{}\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{-1}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{2}~{}&{\scriptstyle\,{\tilde{E}}_{d}<1~{}(t_{bo}>t_{d})}\end{cases}~{}.$
(69)

The minimal energy work-time needed for breakout in the Newtonian regime is
$t_{e,bo}\approx t_{bo}$.

The transition to the relativistic regime can be estimated from the condition
for a relativistic breakout (Eq. 51):

$t_{d,{\rm rel}}=0.1~{}{\rm s}~{}\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{-1}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{2}\left(\frac{\beta_{\max}}{0.3}\right)^{3}~{}.$
(70)

For $t_{d}\gtrsim t_{d,{\rm rel}}$ the head becomes relativistic before it
breaks out of the ejecta. In that case, as long as the jet is collimated, the
minimal engine work-time required for a successful breakout is (Eq. 58):

$t_{e,bo}\approx 0.1~{}{\rm s}~{}\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{\frac{1}{5}}\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{-\frac{1}{5}}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{\frac{2}{5}}\left(\frac{\beta_{\max}}{0.3}\right)^{\frac{2}{5}}\left(\frac{t_{d}}{1{\rm~{}s}}\right)^{\frac{4}{5}}~{}.$
(71)

The transition to an uncollimated jet takes place at (see Eq. 56)

$t_{d,{\rm uncol}}=50~{}{\rm s}~{}\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{-1}\left(\frac{\,{\theta_{j,0}}}{0.05{\rm~{}rad}}\right)^{-\frac{4}{3}}\left(\frac{\beta_{\max}}{0.3}\right)^{-3}~{},$
(72)

where for $t_{d}>t_{d,{\rm uncol}}$ the jet is uncollimated (and surely the
head is ultra-relativistic). The minimal engine work-time required for a
successful breakout in the uncollimated regime is (Eq. 61):

$t_{e,bo}\approx 2~{}{\rm s}~{}\left(\frac{\,{E_{ej,{\rm
tot}}}}{10^{51}{\rm~{}erg}}\right)^{\frac{1}{2}}\left(\frac{\,{L_{j,\rm
iso}}}{10^{52}{\rm~{}erg~{}s^{-1}}}\right)^{-\frac{1}{2}}\left(\frac{\beta_{\max}}{0.3}\right)^{-\frac{1}{2}}\left(\frac{t_{d}}{50{\rm~{}s}}\right)^{\frac{1}{2}}~{}.$
(73)

Note that there is a slight disagreement by a factor of a few in the predicted
value of $t_{e,bo}$ by Eqs. 69 and 71 in the mildly relativistic region where
$t_{d}\sim t_{d,{\rm rel}}$. This is not surprising since both equations are
inaccurate in this regime. The reason for the difference is that the numerical
correction factor $N_{E}$ is different between the Newtonian and relativistic
regimes. For the relativistic regime we use the equations from Harrison et al.
(2018) with no numerical correction for both unmagnetized and magnetized jets
(see Harrison et al. 2018 for a discussion of numerical corrections in the
relativistic regime).

Figure 11: Required jet engine times for breakout as a function of $t_{d}$,
for $\alpha=2$ (top panel) and $\alpha=3$ (bottom panel). Shown are different
values of ${\cal L}\equiv\,{L_{j,\rm iso}}\,{\theta_{j,0}}^{-2}/\,{E_{ej,{\rm
tot}}}$ and magnetizations, while keeping $\,{\beta_{\rm min}}=0.05$ and
$\,{\beta_{\rm max}}=0.3$ fixed. The dashed black lines reflect the minimal
engine time for breakout in the asymptotic phase, assuming it takes place at
$\,{\,{\tilde{t}}_{\rm a}}=4$. Five regimes are present (from small to large
$t_{d}$): i) at small $t_{d}$, $\tilde{E}_{d}<\tilde{E_{a}}$; $t_{e}$ is
independent of $t_{d}$. ii) $\tilde{E}_{d}>\tilde{E_{a}}$; $t_{e}\propto
t_{d}^{2/3}$. iii) transition to relativistic velocities (dashed lines). iv)
the ejecta expanded sufficiently to allow the jet reaching a relativistic
breakout; $t_{e}\propto t_{d}^{4/5}$. v) ejecta is too dilute to collimate the
jet, and a relativistic uncollimated jet breaks out; $t_{e}\propto
t_{d}^{1/2}$.

###  11.1 GW170817

In GW170817 there was a delay of 1.7 s between the merger and the $\gamma$-ray
flash. The afterglow observations show that a jet with the properties given
above broke out of the ejecta successfully and the $\gamma$-rays where
presumably generated by the breakout of the cocoon driven by this jet from the
ejecta or by the jet itself. In either case the observed delay of the gamma-
rays sets an upper limit on the delay in the jet launching and the time it
takes the jet to break out of the ejecta. If the jet is Newtonian, then the
observed delay implies simply $t_{d}+t_{bo}<1.7$ s. If the head is
relativistic, then the limit is $t_{d}+t_{bo}/2\Gamma_{bo}^{2}<1.7$ s. Since
the minimal required engine work-time for breakout is $t_{e,bo}\approx t_{bo}$
in the Newtonian regime and $t_{e,bo}\approx t_{bo}/2\Gamma_{bo}^{2}$ in the
relativistic regime, we conclude that the observations of GW170817 imply that
$t_{d}+t_{e,bo}<1.7$ s. Eqs. 64-73 show that the jet properties, which are
inferred from the afterglow emission ($\,{E_{j,{\rm tot}}}$,
$\,{\theta_{j,0}}$), and the ejecta properties which are measured from the
kilonova ($\,{E_{ej,{\rm tot}}}$, $\,{v_{\rm min}}$, $\,{v_{\rm max}}$), are
consistent with this limit.

However, an interesting question is whether we can use the observations to
obtain a better constraint on $t_{d}$ than the obvious one ($t_{d}<1.7$ s).
Fig. 11 depicts the minimal engine work-time needed for a successful breakout
of the jet from the ejecta, $t_{e,bo}$, as a function of the parameter ${\cal
L}\equiv\,{L_{j,\rm iso}}\,{\theta_{j,0}}^{-2}/\,{E_{ej,{\rm tot}}}$ and the
delay time $t_{d}$. For GW170817 ${\cal L}\approx 10^{3}-10^{4}~{}{\rm
s^{-1}}$, and Fig. 11 shows that for every value of $t_{d}\lesssim 1$ s the
minimal engine work-time is in the range of $t_{e,bo}\sim 0.1-1$ s. This
implies that the observed delay is consistent with being dominated by any of
the various processes involved. Namely, it can be dominated either by $t_{d}$
(e.g., $t_{d}\approx 1.5$ s and $t_{e,bo}\ll 1$ s), or by $t_{e,bo}$ (e.g.,
$t_{d}\ll 1$ s and $t_{e,bo}\approx 1.5$ s), or by other factors such as
geometrical effects, so $t_{d},t_{e,bo}\ll 1$ s, or by any combination of
these processes. Therefore, every value $t_{d}<1.7$ is consistent with the
observed delay between the GWs and $\gamma$-rays and we cannot constrain
$t_{d}$ from the observation of GW170817, apart for the obvious upper limit.

In the shock breakout model for the $\gamma$-ray emission from GW170817, the
breakout emission is released when the cocoon driven by the jet breaks out of
a low mass fast tail that runs ahead of the core of the ejecta (Kasliwal et
al., 2017; Gottlieb et al., 2018b; Beloborodov et al., 2020). The core of the
ejecta is relatively massive ($\sim 0.05{\rm~{}M_{\odot}}$) and it has
$\alpha<5$. In contrast, the fast tail density is expected to be much lower
and most likely it has $\alpha$ much larger than 5, so if the jet breaks out
of the core of the ejecta then it almost certain to break out from the fast
tail as well. However, the breakout from the tail takes place at a much larger
radius and at a much later time than the breakout from the core. In the
breakout model the observed $\gamma$-rays determine the radius and Lorentz
factors of the ejecta and of the shock at the time of the breakout from the
fast tail. This, in turn, determines the time at which the jet should break
out of the core of the ejecta so the shock driven by the cocoon can overtake
the fast tail at the correct radius (Nakar, 2020). In GW170817 the breakout
from the core of the ejecta should be $0.1-1$ s after the merger for the
$\gamma$-rays to be consistent with the shock breakout emission. This range is
consistent with the predicted breakout time derived above based on the jet and
ejecta properties.

###  11.2 Limits on the jet and ejecta parameters from sGRB observations

Fig. 11 has also interesting implications to sGRBs. The engine work-time of
sGRBs cannot be much longer than the observed duration of the prompt emission
(Bromberg et al., 2013), implying that $t_{e,bo}\lesssim 1$ s in sGRBs.
Additionally, there is also an observational evidence that typically
$t_{e,bo}\approx 0.1-1$ s (Moharana & Piran, 2017). Therefore, Fig. 11 shows
that a successful sGRB requires ${\cal L}\gtrsim 10^{3}~{}[10^{2}]~{}{\rm
s^{-1}}$ for unmagnetized [magnetized] jets. The jet of GW170817 is strong and
narrow enough to easily satisfy this criterion with ejecta mass of $\sim
0.05~{}{\rm M_{\odot}}$. However, typical sGRB jets that point towards Earth
are much weaker and most likely also wider. The energy and opening angle of
typical sGRB jets can be estimated based on the properties of their prompt
$\gamma$-rays. The prompt emission is efficient and the total energy carried
by the jet is typically thought to be comparable to the one seen in
$\gamma$-rays, or at most larger by a factor of a few (Nakar, 2007, and
references therein). Studies of the sGRB luminosity function show that the
most common sGRBs, by volumetric rate, are much fainter than the jet seen in
GW170817. In fact, the volumetric rate131313Note that the luminosity function
constrains the rate of sGRBs whose jets are pointing towards Earth. Thus, if
the typical opening angle of more luminous jets is narrower than the opening
angle of less luminous jets then the total rate of luminous and faint jets,
including the beaming correction, may be comparable. The fact that the first
BNS merger that was detected by GWs had a jet with unusually high luminosity
compared to typical sGRBs suggests that this is the case. See (Nakar, 2020)
for discussion. of sGRBs with a peak isotropic equivalent $\gamma$-ray
luminosity $L_{\gamma,{\rm iso}}^{\rm peak}\sim 10^{50}~{}{\rm erg~{}s^{-1}}$
is roughly $\sim 5{\rm~{}Gpc^{-3}yr^{-1}}$ (Wanderman & Piran, 2015, and
references therein). This rate is at least 100 times more frequent than that
of sGRBs with $L_{\gamma,{\rm iso}}^{\rm peak}\sim 10^{52}~{}{\rm
erg~{}s^{-1}}$. Given that the jet luminosity is expected to be larger by a
factor of a few than the $\gamma$-ray luminosity, and that the average
luminosity is a factor of a few lower than the peak luminosity we can
approximate $\,{L_{j,\rm iso}}\sim L_{\gamma,{\rm iso}}^{\rm peak}$. Since the
sGRB luminosity function counts only jets that point towards us, rate
considerations show that the opening angle of these jets cannot be very small.
Given that the total rate of BNS mergers, as found by their gravitational
waves, is $\sim 1000{\rm~{}Gpc^{-3}yr^{-1}}$ (Abbott et al., 2019), jets with
$\,{L_{j,\rm iso}}\sim 10^{50}~{}{\rm erg~{}s^{-1}}$ should have typically
$\,{\theta_{j,0}}\gtrsim 0.1$ rad. Combining these results with the lower
limit on ${\cal L}$ of sGRBs implies that the energy in the ejecta of sGRBs
with $L_{\gamma,{\rm iso}}\sim 10^{50}~{}{\rm erg~{}s^{-1}}$ cannot exceed
$10^{49}~{}[10^{50}]$ erg for unmagnetized [magnetized] jets. Since the jet is
sensitive only to the ejecta in its opening angle, it sets limits on the
isotropic equivalent energy of the ejecta along the poles. Since the
velocities of the ejecta in various mergers are expected to be rather similar,
this result indicates that the isotropically equivalent mass ejected along the
poles in those events is $\lesssim 10^{-3}M_{\odot}$, significantly lower than
the ejecta mass of GW170817, unless the ejecta along the pole in GW170817 was
significantly lower than the ejecta closer to the equator (see also discussion
in Gottlieb et al., 2021a).

Fig. 11 also enables us to learn during which phase the breakout in sGRBs
typically takes place. We have seen that strong jets approach the asymptotic
phase at $t\gtrsim 4t_{d}$. Weak jets break out at $t\gg t_{d}$ and always
after reaching the asymptotic phase. Therefore every jet with
$t_{e}>t_{e,bo}\gtrsim 3t_{d}$ breaks out successfully after approaching the
asymptotic phase, if its head is Newtonian at the breakout (note that a
Newtonian breakout takes place at $t_{d}+t_{e,bo}$). The dashed black line in
Fig. 11 marks $t_{e,bo}\gtrsim 3t_{d}$. Thus, only events with $t_{d}$ and
$t_{e,bo}$ that puts it to the left and above this line break out during the
asymptotic phase. In GW070817 this certainly requires $t_{d}\lesssim 0.1$ s
and possibly even $t_{d}<20$ ms if the jet is weakly magnetized or if its
${\cal L}\approx 10^{4}{\rm~{}s^{-1}}$. For typical sGRBs with lower
luminosities, the breakout upper limit of $\sim 1$ s dictates that the
breakout can take place in the asymptotic regime only if $t_{d}\lesssim 0.2$
s. We therefore conclude that it is possible that many BNS mergers launch jets
that break out or choke before reaching the asymptotic phase.

Another interesting question is whether the breakout takes place before or
after the jet is fully collimated. In §3 we showed that if the cavity is
empty, then a full collimation is obtained roughly when
$\frac{\,{v_{ej,h}}}{\,{v_{\rm min}}}\approx
5\left(\frac{\,{\theta_{j,0}}}{0.1{\rm~{}rad}}\right)^{-2/3}$. If the cavity
is not empty, the collimation is obtained earlier. The expectation, which is
also supported by the observation of GW170817, is that in BNS mergers
$\,{v_{\rm max}}/\,{v_{\rm min}}\approx 6$. This implies that the jet is
likely to spend a significant part of its propagation before being fully
collimated, but it approaches full collimation near the time of the breakout.
As we showed in §8, this result implies that the initial conditions inside the
cavity should have at most a minor effect on the breakout time (and thus on
the breakout criterion) in BNS mergers.

Finally, Fig. 11 can also be used to set an upper limit on $t_{d}$ in sGRBs.
It shows that a successful breakout with $t_{e,bo}<1$ s requires at least
${\cal L}\approx 10^{3}~{}[10^{2}]~{}{\rm s^{-1}}$ for unmagnetized
[magnetized] jet, and that for such values of ${\cal L}$ the delay time has to
satisfy $t_{d}\lesssim 5$ s (regardless of the jet magnetization). If ${\cal
L}$ is larger by an order of magnitude, then $t_{e,bo}<1$ s as long as
$t_{d}\lesssim 50$ s. Assuming that the ejected mass from typical BNS mergers
that produce sGRBs does not fall below $\sim 10^{-3}~{}\,{\rm M_{\odot}}$, we
conclude that the observations of sGRBs imply $t_{d}\lesssim 10$ s. Note that
recently Beniamini et al. (2020) used similar considerations to conclude that
the observations require $t_{d}\lesssim 0.1$s. The two main reasons for the
difference between our conclusions and those of Beniamini et al. (2020), are:
i) their study depends in part on Duffell et al. (2018) (which as we show in
§10, cannot be used for an accurate estimate of $t_{e,bo}$); ii) they ignored
the relativistic regime, which as we show here is the relevant one at large
delay times.

##  12 Summary

Jet-medium interaction plays a key role in the evolution of jets. There are
astrophysical systems, such as BNS mergers and possibly some types of
superluminous SNe, where a relativistic jet encounters an expanding medium.
The jet propagation is then different than in the simple case of a static
medium. Recently, primarily motivated by GW detection of BNS mergers, several
authors have addressed the jet evolution in expanding medium. However, these
studies explored a limited set of configurations and a large area of the
parameter phase space, some of which is directly relevant to BNS mergers,
remained unexplored. In this work we present an analytic model of the jet
propagation in homologous expanding media, which covers a wide range of jet
and ejecta parameters, including the first calculations of propagation in
expanding medium of low-luminosity jets, mildly magnetized jets, jets with a
time-dependent luminosity and jets with a relativistic head.

Our solution makes several assumptions. First, we assume a Newtonian expanding
homologous ejecta with a power-law density profile with index $\,{\alpha}<5$,
so most of the energy is carried by the fast material. The assumption of
homologous flow is certainly valid if all the ejecta is launched before the
jet injection starts. We show that it is valid also in cases where ejecta is
continued to be launched simultaneously with the jet, as long as most of the
ejecta energy is carried by material that was launched before the jet
injection starts. The ejecta does not have to be spherically symmetric, but
the solution assumes that it can be approximated as such over the region where
the jet is passing. Namely, that the ejecta density and velocity distributions
do not vary significantly inside an angle significantly larger than
$\,{\theta_{j,0}}$ around the jet symmetry axis. We further assume that the
jet is launched with an opening angle that does not vary with time and that
the energy launched into the jet grows with time (i.e., if $L_{j}$ drops with
time, it is slower than $t^{-1}$). In order to derive this solution we also
expand the solution of jet propagation in a static medium (see Appendix A),
which was until now limited to a density profile with $\,{\alpha}<3$, to
include steeper density profiles of $3\leq\,{\alpha}<5$. We calibrate our
model using a series of 3D RMHD simulations.

We find that the jet evolution is mainly dictated by a single dimensionless
parameter (that may evolve with time),
$\eta=(v_{h}-\,{v_{ej,h}})/\,{v_{ej,h}}$. Namely, the ratio between the jet
head velocity, as measured in the local ejecta frame, and the ejecta velocity
at the location of the jet head. Different values of $\eta$ define three main
regimes of the evolution:
1\. $\eta\ll 1$: The head of the jet is roughly at rest in the local reference
frame of the ejecta. The velocity of the head in this regime can be
approximated as a constant. The jet cannot remain in this regime indefinitely
and the value of $\eta$ must eventually grow with time.
2\. $\eta\gg 1$: The head velocity is much larger than the ejecta velocity at
the head location. The ejecta velocity can be neglected in the calculation of
the head velocity, which can be approximated very well by the theory of jet
propagation in static media (Bromberg et al., 2011). The value of $\eta$ in
this regime is dropping continuously.
3\. $\eta=\eta_{a}$ (the asymptotic regime): Regardless of the initial
conditions, the value of $\eta$ always converges after enough time to a
constant value141414If the luminosity varies with time as $L\propto t^{k}$
where $k>-1$ then $\eta_{a}=(1+k)/(5-\,{\alpha})$.
$\eta_{a}=1/(5-\,{\alpha})$.

This three regimes define two different types of jets, strong and weak. The
type of the jet is determined by the system initial conditions, and it depends
on another parameter that we define, $\,{\eta_{\rm crit}}$:
(i) Strong jet ($\,{\eta_{\rm crit}}>\eta_{a}$): The jet is strong enough to
start its propagation in the ejecta, either before or after being collimated
by the pressure that builds up in the cavity, within a timescale that is
shorter than $t_{d}$. Note that the delay between the merger and the jet
launching, $t_{d}$, is also the dynamical time of the ejecta when the jet
launching starts. Such jets propagate at first with $\eta\gtrsim 1$ and within
several dynamical times of the ejecta ($t\sim 4t_{d}$) $\eta$ approaches its
asymptotic value, regardless of how large $\eta$ is initially.
(ii) Weak jet ($\,{\eta_{\rm crit}}<\eta_{a}$): The jet is too weak to breach
the ejecta and it is stalled at first at the base of the ejecta. Only after
the ejecta expands significantly and its density becomes low enough, the jet
can break through and propagate inside it. In these jets, initially $\eta\ll
1$, and it grows over a timescale that is much longer than $t_{d}$ until
approaching $\eta_{a}$. As a result, the breakout time is independent of
$t_{d}$, as the effective delay time is the time at which the jet starts
propagating in the ejecta.

We provide simple analytic approximations to the time-dependent location of a
subrelativistic jet head in strong (Eq. 25) and weak (Eq. 29) jets. We derive
a simple breakout criterion that can be applied to any jet with known $t_{d}$
and $t_{e}$, as well as the jet and ejecta total energies and the jet opening
angle (Eq. 41). The same criterion can be applied to a wide range of density
profiles and to jets with time evolving luminosity (as long as
$E_{j,tot}\approx L_{j}(t=t_{d}+t_{e})\cdot t_{e}$), implying that this
criterion is largely independent of the exact density profile and the
luminosity temporal evolution. We also derive a criterion, Eq. 42, that is
necessary but may be insufficient for jet breakout, where $E_{j}$ is known,
but neither $t_{e}$, $L_{j}$ and/or $t_{d}$ are unknown (e.g., such as in
GW170817 where the afterglow reveals only $E_{j}$). If we know that
$t_{e}\gtrsim t_{d}$ in a strong jet or that $t_{e}\gtrsim t_{a}$ in a weak
jet , then satisfying Eq. 42 is sufficient for a successful breakout. All the
analytic formulae that we provide include calibration coefficients that were
obtained via numerical simulations. These coefficients, which depend on
$\,{\alpha}$, are given in Table 4.

Our approximation for the head location includes the early collimation phase
of the jet. This phase depends on the density between the jet launching site
and the slow end of the bulk of the ejecta (moving at $\,{v_{\rm min}}$). In
our analytic expressions we assume that there is an empty cavity at
$r<\,{r_{\rm min}}$. If this is not the case, the evolution at early times can
deviate from our solution for some time, but we show that after the jet is
collimated, the head location becomes independent of the initial density in
the cavity.

We explore the effect of subdominant magnetic field and find that if the jet
maintains a toroidal magnetic field with a magnetization of $\sigma\sim
0.01-0.1$, it propagates 2-3 times faster than an unmagnetized jet. That is
equivalent to an increase of about an order of magnitude in ${\cal
L}=\,{L_{j,\rm iso}}\theta_{j,0}^{-2}/\,{E_{ej,{\rm tot}}}$. Furthermore, such
a magnetic field considerably shortens the breakout time, by an order of
magnitude if the breakout time is $t_{bo}>t_{d}$, and by a factor of two if
$t_{bo}<t_{d}$. Lower magnetizations of $10^{-3}\lesssim\sigma\lesssim
10^{-2}$ are anticipated to be equivalent to a smaller correction in ${\cal
L}$. As the jet is likely to have some degree of magnetization, at least
before the collimation shock, this result implies that magnetic fields are
likely to play a central role in the breakout of jets from dense media.

Applying our model to GW170817 we find that the breakout time of the jet from
the ejecta was most likely $0.1-1$ s, with a very weak dependence on $t_{d}$.
Thus, any $t_{d}<1.7$ s is consistent with the observed delayed $\gamma$-ray
emission. When considering the sample of observed sGRBs, we find that in order
for an sGRB jet to break out within an engine work-time of $t_{e}\lesssim 1$ s
(as required by the observations), it must satisfy two conditions. First, the
system must maintain ${\cal L}\gtrsim 10^{3}~{}[10^{2}]~{}{\rm s^{-1}}$ for
unmagnetized [magnetized] jets. We show that this implies that for the most
common sGRBs with $\gamma$-ray luminosity of $10^{50}~{}{\rm erg~{}s^{-1}}$,
the ejecta mass must be lower than $\sim 10^{-3}{\rm~{}M_{\odot}}$. If the
ejecta is highly anisotropic then this limit is applicable for the isotropic
equivalent ejecta mass along the jet axis. Second, unless ${\cal L}\gg
10^{3}~{}{\rm s^{-1}}$, the jet must be launched after the merger within
$t_{d}\lesssim 10$ s.

## Acknowledgements

We thank Tsvi Piran, Hamid Hamidani, Kunihito Ioka, Paz Beniamini and Ariadna
Murguia-Berthier for useful comments. This research was partially supported by
a consolidator ERC grant (JetNS) and an ISF grant 1995/21.

## Data Availability

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

## References

  * Abbott et al. (2019) Abbott B. P., et al., 2019, Physical Review X, 9, 031040
  * Aloy et al. (2000) Aloy M. A., Müller E., Ibáñez J. M., Martí J. M., MacFadyen A., 2000, ApJ, 531, L119
  * Begelman & Cioffi (1989) Begelman M. C., Cioffi D. F., 1989, ApJ, 345, L21
  * Beloborodov et al. (2020) Beloborodov A. M., Lundman C., Levin Y., 2020, ApJ, 897, 141
  * Beniamini et al. (2020) Beniamini P., Duran R. B., Petropoulou M., Giannios D., 2020, ApJ, 895, L33
  * Bhat et al. (2012) Bhat P. N., et al., 2012, ApJ, 744, 141
  * Blandford & Rees (1974) Blandford R. D., Rees M. J., 1974, MNRAS, 169, 395
  * Bromberg & Tchekhovskoy (2016) Bromberg O., Tchekhovskoy A., 2016, MNRAS, 456, 1739
  * Bromberg et al. (2011) Bromberg O., Nakar E., Piran T., Sari R., 2011, ApJ, 740, 100
  * Bromberg et al. (2013) Bromberg O., Nakar E., Piran T., Sari R., 2013, ApJ, 764, 179
  * Bromberg et al. (2019) Bromberg O., Singh C. B., Davelaar J., Philippov A. A., 2019, ApJ, 884, 39
  * Duffell et al. (2018) Duffell P. C., Quataert E., Kasen D., Klion H., 2018, ApJ, 866, 3
  * Geng et al. (2019) Geng J.-J., Zhang B., Kölligan A., Kuiper R., Huang Y.-F., 2019, ApJ, 877, L40
  * Gill et al. (2019) Gill R., Nathanail A., Rezzolla L., 2019, ApJ, 876, 139
  * Gottlieb & Globus (2021) Gottlieb O., Globus N., 2021, ApJ, 915, L4
  * Gottlieb & Loeb (2020) Gottlieb O., Loeb A., 2020, MNRAS, 493, 1753
  * Gottlieb et al. (2018a) Gottlieb O., Nakar E., Piran T., 2018a, Monthly Notices of the Royal Astronomical Society, 473, 576
  * Gottlieb et al. (2018b) Gottlieb O., Nakar E., Piran T., Hotokezaka K., 2018b, MNRAS, 479, 588
  * Gottlieb et al. (2019) Gottlieb O., Levinson A., Nakar E., 2019, MNRAS, 488, 1416
  * Gottlieb et al. (2020a) Gottlieb O., Levinson A., Nakar E., 2020a, MNRAS, 495, 570
  * Gottlieb et al. (2020b) Gottlieb O., Bromberg O., Singh C. B., Nakar E., 2020b, MNRAS, 498, 3320
  * Gottlieb et al. (2021a) Gottlieb O., Nakar E., Bromberg O., 2021a, MNRAS, 500, 3511
  * Gottlieb et al. (2021b) Gottlieb O., Bromberg O., Levinson A., Nakar E., 2021b, MNRAS, 504, 3947
  * Gottlieb et al. (2022a) Gottlieb O., Liska M., Tchekhovskoy A., Bromberg O., Lalakos A., Giannios D., Mösta P., 2022a, arXiv e-prints, p. arXiv:2204.12501
  * Gottlieb et al. (2022b) Gottlieb O., Moseley S., Ramirez-Aguilar T., Murguia-Berthier A., Liska M., Tchekhovskoy A., 2022b, arXiv e-prints, p. arXiv:2205.01691
  * Gottlieb et al. (2022c) Gottlieb O., Lalakos A., Bromberg O., Liska M., Tchekhovskoy A., 2022c, MNRAS, 510, 4962
  * Gourgouliatos & Komissarov (2018) Gourgouliatos K. N., Komissarov S. S., 2018, Nature Astronomy, 2, 167
  * Govreen-Segal et al. (2021) Govreen-Segal T., Nakar E., Levinson A., 2021, ApJ, 907, 113
  * Hamidani & Ioka (2021) Hamidani H., Ioka K., 2021, MNRAS, 500, 627
  * Hamidani et al. (2020) Hamidani H., Kiuchi K., Ioka K., 2020, MNRAS, 491, 3192
  * Harrison et al. (2018) Harrison R., Gottlieb O., Nakar E., 2018, MNRAS, 477, 2128
  * Ioka & Nakamura (2018) Ioka K., Nakamura T., 2018, Progress of Theoretical and Experimental Physics, 2018, 043E02
  * Ito et al. (2015) Ito H., Matsumoto J., Nagataki S., Warren D. C., Barkov M. V., 2015, ApJ, 814, L29
  * Ito et al. (2021) Ito H., Just O., Takei Y., Nagataki S., 2021, ApJ, 918, 59
  * Kasliwal et al. (2017) Kasliwal M. M., et al., 2017, Science, 358, 1559
  * Kathirgamaraju et al. (2018) Kathirgamaraju A., Barniol Duran R., Giannios D., 2018, MNRAS, 473, L121
  * Kathirgamaraju et al. (2019) Kathirgamaraju A., Tchekhovskoy A., Giannios D., Barniol Duran R., 2019, MNRAS, 484, L98
  * Klion et al. (2021) Klion H., Duffell P. C., Kasen D., Quataert E., 2021, MNRAS, 502, 865
  * Lamb et al. (2022) Lamb G. P., Nativi L., Rosswog S., Kann D. A., Levan A., Lundman C., Tanvir N., 2022, arXiv e-prints, p. arXiv:2201.09796
  * Lazzati & Begelman (2005) Lazzati D., Begelman M. C., 2005, ApJ, 629, 903
  * Lazzati & Perna (2019) Lazzati D., Perna R., 2019, ApJ, 881, 89
  * Lazzati et al. (2009) Lazzati D., Morsony B. J., Begelman M. C., 2009, ApJ, 700, L47
  * Lazzati et al. (2017) Lazzati D., Deich A., Morsony B. J., Workman J. C., 2017, MNRAS, 471, 1652
  * Li et al. (2020) Li X.-J., Zhang Z.-B., Zhang C.-T., Zhang K., Zhang Y., Dong X.-F., 2020, ApJ, 892, 113
  * López-Cámara et al. (2013) López-Cámara D., Morsony B. J., Begelman M. C., Lazzati D., 2013, ApJ, 767, 19
  * López-Cámara et al. (2016) López-Cámara D., Lazzati D., Morsony B. J., 2016, ApJ, 826, 180
  * Lyutikov (2020) Lyutikov M., 2020, MNRAS, 491, 483
  * MacFadyen et al. (2001) MacFadyen A. I., Woosley S. E., Heger A., 2001, ApJ, 550, 410
  * Margalit et al. (2018) Margalit B., Metzger B. D., Thompson T. A., Nicholl M., Sukhbold T., 2018, MNRAS, 475, 2659
  * Margutti & Chornock (2020) Margutti R., Chornock R., 2020, arXiv e-prints, p. arXiv:2012.04810
  * Marti et al. (1995) Marti J. M. A., Muller E., Font J. A., Ibanez J. M., 1995, ApJ, 448, L105
  * Martí et al. (1997) Martí J. M., Müller E., Font J. A., Ibáñez J. M. Z., Marquina A., 1997, ApJ, 479, 151
  * Matsumoto & Kimura (2018) Matsumoto T., Kimura S. S., 2018, ApJ, 866, L16
  * Matsumoto & Masada (2013) Matsumoto J., Masada Y., 2013, ApJ, 772, L1
  * Matsumoto & Masada (2019) Matsumoto J., Masada Y., 2019, MNRAS, 490, 4271
  * Matsumoto et al. (2021) Matsumoto J., Komissarov S. S., Gourgouliatos K. N., 2021, MNRAS, 503, 4918
  * Matzner (2003) Matzner C. D., 2003, MNRAS, 345, 575
  * McBreen et al. (2001) McBreen S., Quilligan F., McBreen B., Hanlon L., Watson D., 2001, A&A, 380, L31
  * Meliani & Keppens (2009) Meliani Z., Keppens R., 2009, The Astrophysical Journal, Volume 705, Issue 2, pp. 1594-1606 (2009)., 705, 1594
  * Mészáros & Waxman (2001) Mészáros P., Waxman E., 2001, Phys. Rev. Lett., 87, 171102
  * Metzger et al. (2018) Metzger B. D., Thompson T. A., Quataert E., 2018, ApJ, 856, 101
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, The Astrophysical Journal Supplement Series, 170, 228
  * Mizuno et al. (2009) Mizuno Y., Lyubarsky Y., Nishikawa K.-I., Hardee P. E., 2009, The Astrophysical Journal, Volume 700, Issue 1, pp. 684-693 (2009)., 700, 684
  * Mizuno et al. (2012) Mizuno Y., Lyubarsky Y., Nishikawa K.-I., Hardee P. E., 2012, The Astrophysical Journal, Volume 757, Issue 1, article id. 16, 14 pp. (2012)., 757
  * Mizuta & Aloy (2009) Mizuta A., Aloy M. A., 2009, ApJ, 699, 1261
  * Mizuta & Ioka (2013) Mizuta A., Ioka K., 2013, The Astrophysical Journal, Volume 777, Issue 2, article id. 162, 17 pp. (2013)., 777
  * Mizuta et al. (2006) Mizuta A., Yamasaki T., Nagataki S., Mineshige S., 2006, ApJ, 651, 960
  * Moharana & Piran (2017) Moharana R., Piran T., 2017, MNRAS, 472, L55
  * Morsony et al. (2007) Morsony B. J., Lazzati D., Begelman M. C., 2007, ApJ, 665, 569
  * Morsony et al. (2010) Morsony B. J., Lazzati D., Begelman M. C., 2010, ApJ, 723, 267
  * Murguia-Berthier et al. (2014) Murguia-Berthier A., Montes G., Ramirez-Ruiz E., De Colle F., Lee W. H., 2014, ApJ, 788, L8
  * Murguia-Berthier et al. (2017) Murguia-Berthier A., et al., 2017, ApJ, 835, L34
  * Murguia-Berthier et al. (2021) Murguia-Berthier A., Ramirez-Ruiz E., De Colle F., Janiuk A., Rosswog S., Lee W. H., 2021, ApJ, 908, 152
  * Nagakura et al. (2011) Nagakura H., Ito H., Kiuchi K., Yamada S., 2011, ApJ, 731, 80
  * Nagakura et al. (2014) Nagakura H., Hotokezaka K., Sekiguchi Y., Shibata M., Ioka K., 2014, ApJ, 784, L28
  * Nakar (2007) Nakar E., 2007, Phys. Rep., 442, 166
  * Nakar (2020) Nakar E., 2020, Phys. Rep., 886, 1
  * Nakar & Piran (2002) Nakar E., Piran T., 2002, MNRAS, 330, 920
  * Nathanail et al. (2021) Nathanail A., Gill R., Porth O., Fromm C. M., Rezzolla L., 2021, MNRAS, 502, 1843
  * Nativi et al. (2022) Nativi L., Lamb G. P., Rosswog S., Lundman C., Kowal G., 2022, MNRAS, 509, 903
  * Pavan et al. (2021) Pavan A., Ciolfi R., Kalinani J. V., Mignone A., 2021, MNRAS, 506, 3483
  * Rossi et al. (2008) Rossi P., Mignone A., Bodo G., Massaglia S., Ferrari A., 2008, Astronomy and Astrophysics, Volume 488, Issue 3, 2008, pp.795-806, 488, 795
  * Salafia et al. (2020) Salafia O. S., Barbieri C., Ascenzi S., Toffano M., 2020, A&A, 636, A105
  * Tchekhovskoy & Bromberg (2016) Tchekhovskoy A., Bromberg O., 2016, Monthly Notices of the Royal Astronomical Society: Letters, Volume 461, Issue 1, p.L46-L50, 461, L46
  * Toma et al. (2017) Toma K., Komissarov S. S., Porth O., 2017, MNRAS, 472, 1253
  * Urrutia et al. (2021) Urrutia G., De Colle F., Murguia-Berthier A., Ramirez-Ruiz E., 2021, MNRAS, 503, 4363
  * Wanderman & Piran (2015) Wanderman D., Piran T., 2015, MNRAS, 448, 3026
  * Wang et al. (2008) Wang P., Abel T., Zhang W., 2008, ApJS, 176, 467
  * Zhang et al. (2004) Zhang W., Woosley S. E., Heger A., 2004, ApJ, 608, 365

##  Appendix A Jet propagation in static media with $3\leq\,{\alpha}<5$

Bromberg et al. (2011) limited their solution of jet propagation in static
medium to $\,{\alpha}<3$. As we show here the solution is actually valid (with
minor revisions) in the range $3\leq\,{\alpha}<5$ as well. We then complete
the numerical calibration of the jet propagation in a static medium (Harrison
et al., 2018) for $\alpha=3$ and $\alpha=4$.

The limit $\,{\alpha}<3$ that Bromberg et al. (2011) set on the solution
originates in their calculation of the average cocoon density under the
approximation of the cocoon shape as a cylinder. Then, if $\,{\alpha}\geq 3$
the mass in the cocoon diverges and there is no solution. This approximation,
however, is irrelevant for large values of $\,{\alpha}$. The reason is that
the shape of the cocoon depends on $\,{\alpha}$, where for $\,{\alpha}<2$ the
cocoon is wide at the base and narrow at the front (the shape of an arrowhead)
and for $\,{\alpha}>2$ it is wide at the front and narrow at the base (an
inverted arrowhead, see figure 3 of Harrison et al. 2018). Only for
$\,{\alpha}=2$ the cocoon looks like a barrel and can be approximated by a
cylinder. Physically, it is clear that the pressure in the cocoon, which
determines the propagation of the head, should not be affected by the density
distribution at radii that are much smaller than the smallest structure of the
system, which is the collimation shock. Therefore, there is no reason why the
solution should not be applicable for density profiles where the mass diverges
at small radii. In what follows we derive the propagation equations in such
density profiles.

Figure 12:  Calibration of the analytic expression of the jet head
dimensionless velocity in a static medium with $\alpha=1$ and $\alpha=2$
(top), and $\alpha=2.5$, $\alpha=3$ and $\alpha=4$ (bottom). Solid lines mark
numerical results and dotted lines are the analytic curves with best fit
$\tilde{N}_{s}$ from Table 5. Time is given in arbitrary units.

We start from Harrison et al. (2018), who derived the jet head dimensionless
velocity in the Newtonian regime (their equation A3):

$\beta_{h}=\frac{N_{s}^{\frac{5}{5-\alpha}}}{c}\left(\frac{\,{L_{j,\rm
iso}}}{A_{\rho}\theta_{j,0}^{2}}\right)^{\frac{1}{5-\alpha}}\left(\frac{3^{3-\alpha}}{3-\alpha}\right)^{\frac{1}{5-\alpha}}\left(\frac{4(5-\alpha)^{\alpha-3}}{\pi}\right)^{\frac{1}{5-\alpha}}t^{\frac{\alpha-2}{5-\alpha}}~{},$
(74)

where $N_{s}$ is a numerical coefficient that was calibrated by their
numerical simulations (Note that in Harrison et al. 2018 the notation $L_{j}$
refers to a one-sided luminosity rather than the two-sided luminosity used
here). The factor
$\left(\frac{3^{3-\alpha}}{3-\alpha}\right)^{\frac{1}{5-\alpha}}$ appears in
this equation due to the assumption of Bromberg et al. (2011) on the cocoon
shape described above and it reflects the divergence of the mass at small
radii for $\,{\alpha}>3$. Thus, during our derivation of Eq. 32, we omitted
this factor and introduced the notation $\tilde{N}_{s}$ instead of
$N_{s}=0.35$, to account for the missing factor:

$\beta_{h}=\frac{\tilde{N}_{s}^{\frac{5}{5-\alpha}}(5-\alpha)^{\frac{\alpha-3}{5-\alpha}}}{c}\left(\frac{\,{L_{j,\rm
iso}}}{A_{\rho}\theta_{j,0}^{2}}\right)^{\frac{1}{5-\alpha}}t^{\frac{\alpha-2}{5-\alpha}}~{}.$
(75)

This notation enables us to extend the solution to $\,{\alpha}<5$.

| $\alpha=1$ | $\alpha=2$ | $\alpha=2.5$ | $\alpha=3$ | $\alpha=4$
---|---|---|---|---|---
$\tilde{N}_{s}$ | 0.47 | 0.38 | 0.25 | 0.19 | 0.16
Table 5: The numerical coefficient $\tilde{N}_{s}$ for jet propagation in a
static medium with power-law index $\alpha$.

We carry out simulations of relativistic jets in static media with
$\,{\alpha}=1,2,2.5,3,4$ to calibrate $\tilde{N}_{s}$ (Table 5). We use
$\alpha=1,2,2.5$ to verify consistency with the results of Harrison et al.
(2018), and $\alpha=3,4$ to extend the calibration to $\tilde{N}_{s}$ of
$\,{\alpha}\geq 3$. The four simulations with integer values of $\alpha$ are
identical to $\,{\it\alpha_{1}-370}$, $\,{\it\alpha_{2}-190}$,
$\,{\it\alpha_{3}-39}$ and $\,{\it\alpha_{4}-50}$ (and the one with
$\alpha=2.5$ has $\,{\eta_{\rm crit}}=140,\eta_{0}=200$), but for the static
case, i.e., $v_{ej}=0$, $t_{d}=0$ and $\,{r_{\rm min}}=z_{\rm beg}$. Fig. 12
depicts the comparison of the analytic (dotted lines) and the numerical (solid
lines) velocities. When the power-law index is $\alpha=1~{}(2)$, the best fit
is obtained for $\tilde{N}_{s}=0.47~{}(0.38)$, indicating that
$N_{s}=0.33~{}(0.29)$ which is rather consistent with $N_{s}=0.35$ of Harrison
et al. (2018). Similarly, for $\alpha=2.5$ we find $\tilde{N}_{s}=0.25$ which
corresponds to $N_{s}=0.19$. This implies that $N_{s}$ (and $\tilde{N}_{s}$)
is not a constant when $\,{\alpha}$ is approaching 3, but decreasing
monotonically with $\alpha$. When $\alpha=3$, Eq. 75 shows that the jet head
velocity accelerates as $t^{0.5}$, in agreement with the numerical simulation.
The best fit coefficient is found to be $\tilde{N}_{s}=0.19$. When $\alpha=4$,
the head accelerates fast as $t^{2}$ to relativistic velocities, and thus we
can fit it over a limited range. Nevertheless, for $\tilde{N}_{s}=0.16$ the
analytic and the numerical curves coincide between $t=6$ and $t=10$ ($t$ is
shown in arbitrary units), including when the jet head becomes relativistic.

##  Appendix B Convergence tests

Figure 13:  Convergence tests for high $\,{\eta_{\rm crit}}$ model
$\,{\it\alpha_{2}-190}$ (left) and low $\,{\eta_{\rm crit}}$ model
$\,{\it\alpha_{2}-1}$ (right). Shown are comparisons of $\eta$ (top panels)
and $r_{h}$ (bottom panels) between original grid resolution (blue), 3/4 the
number of cells (red), and 1/2 the number of cells in the grid (yellow).

Highly demanding 3D grids are naturally run in lower grid resolutions compared
to 2D simulations. Thus, it is particularly important to verify that the
results are independent of the grid resolution chosen for the simulations. We
present two convergence tests for high and low $\,{\eta_{\rm crit}}$ values,
for which different grids were employed. For the numerical convergence tests
we choose our canonical $\alpha=2$ models: $\,{\it\alpha_{2}-190}$ and
$\,{\it\alpha_{2}-1}$ for high and low $\,{\eta_{\rm crit}}$ models,
respectively. For each of those models we run additional two simulations whose
grids are distributed in the same way of our original grids, but their number
of cells is 3/4 and 1/2 of that in the original grid. Fig. 13 depicts the
comparison of $\eta$ (top) and $r_{h}$ (bottom) for models
$\,{\it\alpha_{2}-190}$ (left) and $\,{\it\alpha_{2}-1}$ (right) between the
different simulations. Both quantities show a remarkable agreement between the
different grid resolutions. The $\eta$ plots show full convergence between all
resolutions. Model $\,{\it\alpha_{2}-190}$ shows no clear trend between one
resolution to another, and model $\,{\it\alpha_{2}-1}$ indicates that it takes
more time for lower resolutions to merge with the higher resolution curve. The
$r_{h}$ plots show very small deviations, with somewhat convergence trends
towards our original resolution which moves somewhat slower, as expected from
higher resolution simulations (Harrison et al., 2018).