The propagation of choked jet outflows in power-law external media

TL;DR

This paper models the evolution of choked relativistic jets and their cocoons in power-law media, providing analytical formulas and comparing them with simulations to understand their shape and observational signatures.

Contribution

It introduces an analytical framework for the shape and evolution of cocoons from choked jets in power-law environments, validated against numerical simulations.

Findings

Analytical formulas for cocoon shape and evolution in different regimes.

Agreement with simulations within 15% up to 45 degrees from the axis.

Cocoon shape diagnostics can inform jet properties long after activity ceases.

Abstract

Observations of both gamma-ray bursts (GRBs) and active galactic nuclei (AGNs) point to the idea that some relativistic jets are suffocated by their environment before we observe them. In these "choked" jets, all the jet's kinetic energy is transferred into a hot and narrow cocoon of near-uniform pressure. We consider the evolution of an elongated, axisymmetric cocoon formed by a choked jet as it expands into a cold power-law ambient medium $\rho \propto R^{-\alpha}$, in the case where the shock is decelerating ($\alpha<3$). The evolution proceeds in three stages, with two breaks in behaviour: the first occurs once the outflow has doubled its initial width, and the second once it has doubled its initial height. Using the Kompaneets approximation, we derive analytical formulae for the shape of the cocoon shock, and obtain approximate expressions for the height and width of the outflow versus time in each of the three dynamical regimes. The asymptotic behaviour is different for flat ($\alpha \le 2$) and steep ($2 < \alpha < 3$) density profiles. Comparing the analytical model to numerical simulations, we find agreement to within $\sim15$ per cent out to 45 degrees from the axis, but discrepancies of a factor of 2-3 near the equator. The shape of the cocoon shock can be measured directly in AGNs, and is also expected to affect the early light from failed GRB jets. Observational constraints on the shock geometry provide a useful diagnostic of the jet properties, even long after jet activity ceases.

Figures

Tables

Table 1. Table 1: Asymptotic cocoon properties for the regimes discussed in Section  3 .

	$\alpha \le 2b^2/a^2$		$\alpha =2$
	$\left(k_{\alpha }\ge 1-b^2/a^2\right)$		$\left(k_{\alpha }=0\right)$
$f_{\alpha }$	1			${\left[\sin \left(\pi /\alpha \right)\right]}^{\alpha /\left(\alpha -2\right)}$
$A_{\alpha }$	$\left(2/k_{\alpha }^2\right){\sin }^2\left(k_{\alpha }\pi /4\right)$		${\pi }^2/8$	$\left(1/2k_{\alpha }^2\right){\cot }^2\left(\pi /\alpha \right)$
$\zeta$	$\left(1/k_{\alpha }\right)\left[{\left(z_c/a\right)}^{k_{\alpha }}-1\right]$		$\ln \left(z_c/a\right)$	$\left(-1/k_{\alpha }\right)\left[{\left(z_c/a\right)}^{-k_{\alpha }}-1\right]$
bulge forms?	no	yes
asymptotic shape	sphere			‘peanut’

Table 2. Table 2: Glossary of symbols.

Symbol	Meaning	Introduced
Coordinates
$R$, $\theta$	polar coordinates
$r$, $z$	cylindrical coordinates
Initial conditions
$E_0$	energy of the cocoon	Section 2
$a$	initial height of the cocoon	Section 2
$b$	initial width of the cocoon	Section 2
$P_0$	initial pressure in the cocoon	Section 2
$V_0$	initial volume in the cocoon	Section 2
${\rho }_0$	ambient density at $R=a$	eq. 11
$\alpha$	power-law index of the density profile	eq. 11
$k_{\alpha }$	frequently appearing quantity defined by $k_{\alpha }\equiv 1-\alpha /2$	eq. 29
$R_i\left({\theta }_i\right)$	curve describing the shock’s initial shape	eq. 13
${\chi }_i\left({\theta }_i\right)$	angle between the axis and the initial direction of motion	eqs. 12 and 14
${\omega }_i\left({\theta }_i\right)$	angle between the initial direction of motion and the radial direction	eq. 22
Time-scales
$t_0$	time when the jet is choked	Section 2
$t_b$	time it takes for the cocoon’s width to double	eq. 1
$t_a$	time it takes for the cocoon’s height to double	eq. 2
$t_{\mathrm{sph}}$	time for the cocoon’s width to become 90 per cent of its height	eqs. 78 and 79
$t_{\mathrm{eq}}$	time to reach pressure equilibrium with the ambient medium	eqs. 81 and 82
$t_w$	time when the solution breaks down in an expanding medium	eqs. 83 and 84
$t_{\mathrm{bo}}\left(R_{\ast }\right)$	time when the cocoon breaks out of an ambient medium extending to $R_{\ast }$	eqs. 85 and 86
Cocoon properties
$V_c\left(t\right)$	volume of the cocoon	Section 2
$P_c\left(t\right)$	average pressure in the cocoon	Section 2
$\lambda$	dimensionless ratio of postshock pressure to volume-averaged pressure	eq. 16
$x\left(t\right)$	dimensionless time parameter used to solve for the shock’s shape	eq. 18
$z_c\left(t\right)$	height of the cocoon	eq. 25
$r_c\left(t\right)$	width of the cocoon	eq. 26
$r_{\mathrm{eq}}\left(t\right)$	width of the cocoon at the equator	see Fig. 5
${ϵ}_{\mathrm{eff}}\left(t\right)$	effective eccentricity, i.e. ${ϵ}_{\mathrm{eff}}\equiv \sqrt{1-r_c^2/z_c^2}$	Sections 3.1 and 3.2
$R({\theta }_i,t)$, $\theta ({\theta }_i,t)$	Parametric equations for the shape of the shock at time $t$	Sections 3.1 and 3.2
Quantities related to the bulge
$R_b\left(t\right)$, ${\theta }_b\left(t\right)$	polar coordinates of the bulge	eq. 36 and 37
$r_b\left(t\right)$, $z_b\left(t\right)$	cylindrical coordinates of the bulge	see Fig. 5
${\theta }_{\mathrm{b0}}\left(t\right)$, $R_{\mathrm{b0}}\left(t\right)$, ${\chi }_{\mathrm{b0}}\left(t\right)$	initial values of ${\theta }_i$, $R_i$, and ${\chi }_i$ for a trajectory passing through ($R_b,{\theta }_b$)	see Fig. 6
Asymptotic quantities
$\zeta \left(t\right)$	quantity important in the expansion of $r_c\left(z_c\right)$ at infinity	eq. 43
$A_{\alpha }$	coefficient appearing in the expansion of $r_c\left(z_c\right)$ at infinity	eq. 44 and  63
$f_{\alpha }$, $g_{\alpha }$	ratios of $r_c/z_c$ and $r_{\mathrm{eq}}/z_c$ respectively, as $t\to \infty$	eq. 39
$C_{\alpha }$	ratio of $V_c$ to $4\pi z_c^3/3$ as $t\to \infty$	eq. 70 and 71
${\theta }_{b\infty }$	value of ${\theta }_b\left(t\right)$ as $t\to \infty$	eq. 40
${\theta }_{\mathrm{b0},\min }$	value of ${\theta }_{\mathrm{b0}}\left(t\right)$ as $t\to \infty$	Section 3.2
Jet properties
$L_j$	jet luminosity	Section 4
$t_j$	jet duration	Section 4
${\theta }_{\mathrm{op}}$	jet opening angle at injection	Section 4
${\beta }_h$, ${\Gamma }_h$	velocity and Lorentz factor of the jet head upon choking	Section 4
${\beta }_c$	sideways velocity of the cocoon upon choking	Section 4
$\tilde{L}$	dimensionless parameter governing jet propagation	Section 4

Equations

$$t_b\sim t_0\sim {\left(\frac{24\pi \rho b^{4}a}{7E_0}\right)}^{1/2}.$$

$$t_a\sim {\left(\frac{48\pi \rho a^5}{7E_0}\right)}^{1/2}.$$

$$z_{\mathrm{c}}\sim \{\begin{array}{ll}a,& t-t_0\ll t_b\\ a\left[1+{\left(\frac{t}{t_a}\right)}^{1/2}\right]& t_b\ll t-t_0\ll t_a\\ 2a{\left(\frac{t}{t_a}\right)}^{2/(5-\alpha )},& t-t_0\gg t_a\end{array}$$

$$r_{\mathrm{c}}\sim \{\begin{array}{ll}b,& t-t_0\ll t_b\\ a{\left(\frac{t}{t_a}\right)}^{1/2},& t_b\ll t-t_0\ll t_a\\ a{\left(\frac{t}{t_a}\right)}^{2/(5-\alpha )},& t-t_0\gg t_a\end{array}.$$

$${\theta }_{\mathrm{b}}\sim \{\begin{array}{ll}{\left(\frac{t}{t_a}\right)}^{1/2},& t_b\ll t-t_0\ll t_a\\ 1,& t-t_0\gg t_a\end{array}$$

$$z_{\mathrm{c}}-r_{\mathrm{c}}\propto z_{\mathrm{c}}^{(\alpha -2)/2}\propto t^{(\alpha -2)/(5-\alpha )}$$

$$a\sim \{\begin{array}{ll}{z}_{\mathrm{c}}-r_{\mathrm{c}},& z_{\mathrm{c}}-r_{\mathrm{c}}\text{ga}{z}_{\mathrm{c}}/2\\ [z_{\mathrm{c}}^{-\alpha }(z_{\mathrm{c}}-r_{\mathrm{c}}{)}^2{]}^{1/(2-\alpha )},& z_{\mathrm{c}}-r_{\mathrm{c}}\ll z_{\mathrm{c}}/2.\end{array}$$

$$L_{\mathrm{j}}{t}_{\mathrm{j}}=E_0$$

$$t_{\mathrm{j}}{\theta }_{\mathrm{op}}^{-2}\sim {\left(\frac{\rho (z_{\mathrm{c}})z_{\mathrm{c}}^5}{E_0}\right)}^{1/2}{\left(\frac{z_{\mathrm{c}}-r_{\mathrm{c}}}{z_{\mathrm{c}}}\right)}^{q_{\alpha }},$$

$$q_{\alpha }\equiv \{\begin{array}{ll}(5-\alpha )/2,& z_{\mathrm{c}}-r_{\mathrm{c}}\text{ga}{z}_{\mathrm{c}}/2\\ (5-\alpha )/(2-\alpha ),& z_{\mathrm{c}}-r_{\mathrm{c}}\ll z_{\mathrm{c}}/2\end{array}.$$

$$\rho ={\rho }_0(R/a{)}^{-\alpha }.$$

$${\chi }_i({\theta }_i)={\tan }^{-1}\left(\frac{\tan {\theta }_i-\frac{\mathrm{d}\ln R_i}{\mathrm{d}{\theta }_i}}{1+\frac{\mathrm{d}\ln R_i}{\mathrm{d}{\theta }_i}\tan {\theta }_i}\right).$$

$$R_i({\theta }_i)=\frac{ab}{\sqrt{a^2{\sin }^2{\theta }_i+b^2{\cos }^2{\theta }_i}}$$

$${\chi }_i({\theta }_i)={\tan }^{-1}\left(\frac{a^2}{b^2}\tan {\theta }_i\right).$$

$$\sqrt{\frac{\gamma +1}{2}\frac{P^{\prime }(t)}{\rho (R)}}=v(R,t)=\frac{\left|\partial \theta /\partial t\right|}{\sqrt{1/R^2+(\partial \theta /\partial R{)}^2}},$$

$$P^{\prime }(t)=\frac{\lambda (\gamma -1)E_0}{V_{\mathrm{c}}(t)},$$

$$\lambda \equiv \left(\frac{E_{\mathrm{th}}}{E_0}\right)\left(\frac{P^{\prime }}{P_{\mathrm{c}}}\right)$$

$$x(t)\equiv {\int }_{t_0}^t\sqrt{\frac{\gamma +1}{2}\frac{P_{\mathrm{c}}(t^{\prime })}{{\rho }_0{a}^2}}\mathrm{d}{t}^{\prime }.$$

$$x(t)={\int }_{t_0}^t\sqrt{\frac{P_{\mathrm{c}}(t^{\prime })}{P_0}}\frac{v_0\mathrm{d}{t}^{\prime }}{a}.$$

$$\frac{1}{a^2}{\left(\frac{\partial \theta }{\partial x}\right)}^2={\left(\frac{R}{a}\right)}^{\alpha }\left[{\left(\frac{\partial \theta }{\partial R}\right)}^2+\frac{1}{R^2}\right],$$

$$R(x,{\theta }_i)={\left(R_i^2+a^2{x}^2+2ax{R}_i\cos {\omega }_i\right)}^{1/2},$$

$${\omega }_i\equiv {\chi }_i-{\theta }_i.$$

$$\theta (x,{\theta }_i)={\theta }_i+{\tan }^{-1}\left(\frac{x\sin {\omega }_i}{R_i/a+x\cos {\omega }_i}\right).$$

$$R(\theta ,{\theta }_i)=\frac{R_i}{\cos (\theta -{\theta }_i)-\sin (\theta -{\theta }_i)\cot {\omega }_i}.$$

$$z_{\mathrm{c}}(x)\equiv R(x,{\theta }_i=0),$$

$$r_{\mathrm{c}}(x)\equiv \max \{R(x,{\theta }_i)\sin [\theta (x,{\theta }_i)]\}$$

$$\frac{r_{\mathrm{c}}}{z_{\mathrm{c}}}=\frac{x+b/a}{x+1}=1-\frac{a}{z_{\mathrm{c}}}\left(1-b/a\right).$$

$${ϵ}_{\mathrm{eff}}\approx {\left[2\left(1-b/a\right)a/z_{\mathrm{c}}\right]}^{1/2}$$

$$k_{\alpha }\equiv \frac{2-\alpha }{2},$$

$$\frac{1}{k_{\alpha }^2}{\left(\frac{\partial \varphi }{\partial x}\right)}^2={\left(\frac{\partial \varphi }{\partial \xi }\right)}^2+\frac{1}{{\xi }^2},$$

Full text

††pubyear: 2019††pagerange: The propagation of choked jet outflows in power-law external media–LABEL:lastpage

1 Introduction

Relativistic jets are ubiquitous in high-energy astrophysics and essential to our understanding of phenomena such as short and long gamma-ray bursts (GRBs) and active galactic nuclei (AGNs). Although these systems have markedly different origins—AGNs are powered by accretion on to a supermassive black hole (e.g., Rees, 1984; Blandford et al., 2018, and references therein); long GRBs (LGRBs) are connected to the death of high-mass stars (see, e.g., Woosley & Bloom, 2006, for a review); and, thanks to GW 170817, short GRBs (SGRBs) have at long last been unambiguously linked to the merger of two neutron stars (e.g., Abbott et al., 2017)—the underlying physics is none the less similar, as ultimately they all involve a compact object driving a jet into an external medium. The relevant medium is the intracluster medium in the case of AGNs, the host star and surrounding circumstellar gas in the case of LGRBs, or the ejecta thrown out by the merger in the case of SGRBs.

While the jet-launching process remains murky, and the emission mechanism of relativistic jets is a long-standing problem, the dynamics of jet propagation are comparatively well-understood, both analytically (e.g., Begelman & Cioffi, 1989; Mészáros & Waxman, 2001; Matzner, 2003; Lazzati & Begelman, 2005; Bromberg et al., 2011) and numerically (e.g., 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; Mizuta & Aloy, 2009; Lazzati et al., 2009; Nagakura et al., 2011; López-Cámara et al., 2013; Mizuta & Ioka, 2013; Ito et al., 2015; Harrison et al., 2018). Whether describing an AGN or a GRB, the dynamics are captured by the same physical quantities: the kinetic power $L_{\rm j}$ and Lorentz factor $\Gamma_{j}\gg 1$ of the jet, the opening angle $\theta_{\rm op}$ into which the jet is injected, and the external density $\rho$. In AGN environments, the ambient pressure $P_{a}$ and magnetic field may also be important. As the jet drills through its surroundings, a forward and reverse shock structure known as the jet head develops at the interface between the jet ejecta and the ambient gas. Due to a strong pressure gradient in the jet head, matter flowing in through the forward and reverse shocks is squeezed out to the sides, forming a hot ‘cocoon’ full of shocked ejecta. This cocoon, which is an unavoidable consequence of jet propagation, surrounds and exerts pressure on the jet. Because this pressure can be significant enough to collimate the jet, the dynamics of the jet-cocoon system must be solved self-consistently, as described by Bromberg et al. (2011).

The cocoon contains both an inner part comprised of relativistic jet ejecta that entered through the reverse shock, and an outer part containing non-relativistic ambient matter swept up by the forward shock. The two components have different densities and temperatures—the relativistic gas is light and hot, whereas the non-relativistic gas is heavy and cold—and are likely mixed together to some degree (Nakar & Piran, 2017; Gottlieb et al., 2018). Importantly, however, the different parts of the cocoon are roughly in pressure equilibrium, provided that the cocoon expands non-relativistically, because the near-relativistic sound speed in the lighter regions acts to smooth out pressure differences on time-scales shorter than a dynamical time (Bromberg et al., 2011). The assumption of near-uniform pressure, which is crucial to our model, is supported by numerical jet simulations (e.g., Mizuta & Aloy, 2009; Mizuta & Ioka, 2013; Harrison et al., 2018).

The above description is valid while the jet remains active and matter continues to flow into the reverse shock. However, there is reason to believe that some GRB jets fail to penetrate their immediate surroundings. In these ‘choked’ jets, the reverse shock crosses the flow before the jet breaks out, and all of the jet energy is dumped into the cocoon. Several lines of evidence point towards the existence of choked jets. First, the duration distribution of LGRBs has a plateau, with few objects having an observed duration much shorter than the typical breakout time, indicating a significant population of failed jets (Bromberg et al., 2012). A similar plateau has been observed in the duration distribution of SGRBs, as well (Moharana & Piran, 2017). Second, in low-luminosity GRBs ($ll$GRBs), a peculiar faint and long-lived class of long GRBs, early optical emission suggests the presence of an extended, optically thick envelope (Nakar, 2015; Irwin & Chevalier, 2016). The inferred mass ($\sim 10^{-2}M_{\odot}$) and radius ($\ga 100R_{\odot}$) of the envelope are more than sufficient to choke a GRB of typical luminosity and duration (Nakar, 2015; Irwin & Chevalier, 2016). Mildly relativistic shock breakout from such an extended envelope could explain the unusual prompt emission of $ll$GRBs (Kulkarni et al., 1998; Campana et al., 2006; Nakar & Sari, 2012; Nakar, 2015). Interestingly, although they are more difficult to observe, $ll$GRBs may be more common per cosmic volume than standard GRBs (Soderberg et al., 2006), again hinting that choked jets may be common. Third, early spectroscopy of Type Ib and Ic supernovae (SNe) has unveiled a distinct high-velocity component in several cases (e.g., Piran et al., 2019; Izzo et al., 2019). The inferred energy ($\sim 10^{51}$ erg s*-1*) and velocity ($\sim 0.1c$) of this component are consistent with the expectations for a GRB jet’s cocoon. Finally, in AGNs, the evidence for choked jets is even more direct: some galaxies contain ‘relic bubbles’ that were likely formed by past jet activity (e.g. Churazov et al., 2000; McNamara et al., 2000; Fabian et al., 2002; McNamara et al., 2005; Tang & Churazov, 2017). In recently quenched systems where the bubbles have not yet been deformed by buoyancy, the shape of the bubbles can be used to constrain the jet properties (Irwin et al., in preparation).

Compared to the case where the jet is active, the evolution of the outflow after the jet is choked is less clear. The expectation is that the outflow will eventually become self-similar, but our understanding of the transition to this scale-free regime is lacking. While several authors have explored the process of jet choking through numerical simulations (e.g. Lazzati et al., 2012; Gottlieb et al., 2018), there is not yet a firm analytic theory underpinning the results. We aim to rectify this situation by developing an analytical model for the propagation of a choked jet outflow in a power-law external density profile, thereby extending the solutions of Bromberg et al. (2011) to beyond the moment of choking. To do so, we employ the well-known Kompaneets approximation (hereafter, KA; Kompaneets, 1960), which takes advantage of the fact that the cocoon pressure is nearly uniform (as discussed above).

Before further discussing aspherical outflows, we briefly review important results in the spherical case. In the scale-free limit, a spherical explosion expands according to the well-known Sedov–Taylor (ST) blast wave solution, with the radius growing in time as $R\propto t^{2/(5-\alpha)}$ if the density obeys $\rho\propto R^{-\alpha}$ (Taylor, 1950; Sedov, 1959). The ST solution applies when the outflow is decelerating, which is true for $\alpha<3$. In media of finite mass ($\alpha>3$), the shock accelerates with time and eventually becomes causally disconnected from the swept up mass (Koo & McKee, 1990). The case $3<\alpha<5$ can be described by a second type of self-similar solution (Koo & McKee, 1990; Waxman & Shvarts, 1993), although the ST formula remains a valid approximation until causal connection is lost. For $\alpha>5$, however, the ST solution breaks down completely, since the radius goes to infinity in a finite time, and a different class of self-similar solution applies (Waxman & Shvarts, 1993).

While aspherical explosions have been considered in the past, the literature has mainly focused on point explosions. Early work on this topic was carried out by Kompaneets (1960), who developed the eponymous approximation and obtained a solution for the case of an exponentially stratified density. Later authors applied the KA to a wide variety of other problems, as reviewed by Bisnovatyi-Kogan & Silich (1995) (see also Lyutikov, 2011; Bannikova et al., 2012; Rimoldi et al., 2015, and references therein). The work of Korycansky (1992) (hereafter K92), who investigated off-centre point explosions in a power-law ambient medium, is of particular relevance to our study. Also closely related is the case of an initially spherical explosion with a velocity depending on polar angle, which was treated by Bisnovatyj-Kogan & Blinnikov (1982) for a homogeneous medium.

We consider a different generalization of the Kompaneets problem to power-law ambient media, in which the explosion is centred at the origin but elongated along the symmetry axis. This configuration arises naturally when the energy is injected by a relativistic jet. We restrict our discussion to the $\alpha<3$ case where the shock is decelerating and the flow evolves like an ST blast wave at late times; the case of an accelerating shock will be treated in a separate paper. Applying the KA, we derive general analytical expressions for the evolution of the shock’s shape, and examine how the outflow transitions from an initially elongated shape to a quasi-spherical one. As we will show, although the evolution eventually comes to resemble a point explosion, deviations from spherical symmetry persist long after the jet has been quenched, especially when the initial shape is narrow. Meaningful information about the jet can therefore be extracted by measuring the shape of the shock, even when jet activity has already ended (as in AGN with relic bubbles) or was hidden from view (as in failed GRB jets).

The KA is powerful, but our idealized model is not without its limitations. First, as mentioned above, the model depends on the cocoon pressure being uniform, and will not give reliable results if this is not the case. Second, we only consider the case where the cocoon propagates non-relativistically, and the speed of the ambient gas is negligible. This assumption is accurate for LGRBs and AGNs, and is somewhat reasonable for $ll$GRBs, which are at most mildly relativistic, but is not suitable for SGRBs, where the ejecta may have an expansion speed comparable to the cocoon’s speed. Third, we ignore gravity. This is fine for GRBs where the gravitational time-scale is much longer than a dynamical time; however, in AGNs, our model breaks down once buoyancy becomes important. In spite of these drawbacks, the idealized model presented here is an important first step towards understanding the ultimate fate of choked jet explosions.

The paper has a somewhat unconventional structure in that we begin with a summary of key results in Section 2, where we overview the important time-scales and the role of the density profile in shaping the evolution of the choked jet outflow. Then, in Section 3, we use the KA to derive analytical solutions for the shape of the shock. We start by considering the case of uniform density in Section 3.1, then generalize to a power-law external density in Section 3.2. We find different behaviour in each of the cases $\alpha<2$, $\alpha=2$, and $2<\alpha<3$, which are discussed respectively in subsections 3.2.1–3.2.3. We calculate the volume of the expanding cocoon in Section 3.3, and provide approximate formulae for the cocoon’s height and width as functions of time in Section 3.4. We then wrap up this section with a discussion of the long-term breakdown of the solution in Section 3.5. Next, in Section 4, we consider how the initial conditions of the cocoon problem depend on the properties of the injected jet, and conversely how observationally-inferred cocoon properties can constrain the underlying jet. After that, we compare our analytical results to numerical simulations of choked jets in Section 5. Finally, we present our conclusions and discuss some possible applications in Section 6.

2 Overview

Consider a relativistic jet propagating in a power-law external medium, $\rho\propto R^{-\alpha}$, with a head velocity $c\beta_{\rm h}$. As the jet drills into the ambient gas, shocked material is pushed to the side to form a hot cocoon around the jet. Assuming that the ambient pressure is negligible, this cocoon expands sideways supersonically, driving a shock into the ambient medium at a speed $c\beta_{\rm c}$. Since the sideways expansion is driven by the cocoon pressure alone, while the jet head is pushed forward by both internal pressure and ram pressure from the jet, we necessarily have $\beta_{\rm h}>\beta_{\rm c}$.

After the central engine shuts off, jet material continues to flow into the cocoon until the last-emitted material catches up with the head. Once all of the jet ejecta have entered the cocoon, at the time $t_{0}$, we consider the jet to be ‘choked.’ The time $t_{0}$ is the endpoint of the Bromberg et al. (2011) jet propagation model, and is the starting point for our choked jet model. The evolution of the system after the jet is choked can be described by three parameters: the height $a$ and width $b$ of the outflow upon choking, and the total energy $E_{0}$ injected into the cocoon. (These quantities are directly related to the jet properties in the Bromberg et al. (2011) model, as discussed in Section 4). From the condition $\beta_{\rm h}>\beta_{\rm c}$, it follows that $a>b$.

Elongated, axisymmetric explosions are best understood through analogy with the familiar spherical case. In the case of a spherical explosion of energy $E_{0}$, with initial size $R_{0}$, detonated at time $t_{0}$, there is one length scale ($R_{0}$), and one characteristic time-scale, ${t_{\rm ST}\sim(\rho R_{0}^{5}/E_{0})^{1/2}}$, which is roughly the time for the outflow to double its initial radius. This divides the evolution into two dynamical regimes: the planar phase ($t-t_{0}\ll t_{\rm ST}$), and the self-similar phase ($t-t_{0}\gg t_{\rm ST}$). During the planar phase, outflow properties such as the volume, pressure, and expansion speed are nearly constant in time, with values that depend on $R_{0}$. Conversely, in the self-similar phase, these quantities vary in time, but the dependence on the initial conditions is lost.

The axisymmetric case, on the other hand, has two inherent length scales: the initial width ($b$), and the initial height ($a$). Consequently, there are two characteristic time-scales, and three dynamical regimes. The relevant time-scales are the time for the initial width to double ($t_{b}$), and the time for the initial height to double ($t_{a}$). The early phase of evolution ($t-t_{0}\ll t_{b}$) is analogous to the planar phase of a spherical explosion. Also, like in the spherical case, the expansion becomes self-similar when $t-t_{0}\gg t_{a}$. However, in the axisymmetric solution there is a transitional regime ($t_{b}\ll t-t_{0}\ll t_{a}$) which does not appear in the spherical case. During this period, the outflow’s width changes significantly, while its height remains about the same. The more elongated the explosion, the more pronounced is the transitional regime.

To put it another way, spherical evolution corresponds to a special case of the axisymmetric solution, with $a=b$ and $t_{b}=t_{a}$, where pressure changes become important at around the same time that initial conditions become unimportant. In the more general case of an elongated axisymmetric explosion, however, changes in pressure start to affect the outflow while the initial size is still relevant.

In Fig. 1, we illustrate the shape of the shock in each dynamical regime for several different density profiles, taking as an initial condition an ellipsoidal cocoon shock with $b/a=0.1$. The system is initially in the planar phase, with $t-t_{0}\ll t_{b}$ (upper left panel). In this regime, the shape remains roughly ellipsoidal. The evolution is not very sensitive to the density profile, although minor differences appear towards the equator. Once the width has roughly doubled, the system transitions to the sideways expansion phase (upper right panel), for which $t_{b}\ll t-t_{0}\ll t_{a}$. The shock shape in this case is strongly influenced by the initial height $a$, but only weakly dependent on the initial width $b$. The effect of the density profile also starts to become apparent, with the cocoon being wider towards the base in shallower density gradients, and wider towards the tip in steeper gradients. However, the aspect ratio of the shock is similar regardless of density profile. Over a time-scale $\sim t_{a}$, the height and width of the shock become comparable, and the outflow enters the quasi-spherical regime (lower left panel). In this phase, $t-t_{0}\gg t_{a}$, and the outflow’s size is much larger than its initial size, so the evolution is approximately scale-free. The differences in shape are mainly due to the density profile. Different density profiles no longer have a similar aspect ratio; outflows in shallow density profiles are closer to becoming spherical. Finally, as $t\rightarrow\infty$, the system approaches an asymptotic self-similar solution, as shown in the lower right panel.

The asymptotic behaviour can be divided into two different regimes depending on $\alpha$. For $\alpha\leq 2$, the shock becomes fully spherical and the system is described by the classical spherical ST blast wave solution at late times. Interestingly, however, non-sphericity persists indefinitely in the KA when $\alpha>2$. Although the shock still becomes self-similar in this case, it does not become spherical: the shock radius remains smaller at the equator than at the poles, resulting in a peanut-like shape. The reasons for the different shapes will be explored further in Section 3.2.3.

Our model assumes that the pressure in the ambient medium is negligible compared to the cocoon pressure, and that the ambient medium is effectively static, with an expansion velocity much smaller than the shock velocity. If either of these conditions are violated, the solution breaks down. The range of validity of the model is discussed further in Section 3.5.

Because the sideways expansion is driven by pressure in the cocoon, both before and after choking, the width of the outflow doubles over a time-scale comparable to the choking time, i.e. $t_{b}\sim t_{0}$. Alternatively, we can write $t_{b}\sim b/c\beta_{\rm c}$, with the sideways shock expansion speed $\beta_{\rm c}$ determined by balancing the upstream and downstream momentum flux at the shock. Applying the shock jump conditions gives $\beta_{\rm c}=[(\gamma+1)P_{0}/2\rho c^{2}]^{1/2}$, where $P_{0}=(\gamma-1)E_{0}/V_{0}$ is the initial pressure in the cocoon, $V_{0}\approx\frac{4}{3}\pi ab^{2}$ is the cocoon’s initial volume, and $\gamma$ is the adiabatic index. Adopting $\gamma=4/3$, we arrive at an expression for $t_{b}$ in terms of the cocoon parameters $a$, $b$, and $E_{0}$:

[TABLE]

The time-scale $t_{a}$ is also straightforward to derive. After the jet is choked, both the forward and sideways expansion are driven by the cocoon’s pressure. Most of the early expansion occurs near the tip of the jet where the density is lowest. Forward-moving material near the tip encounters nearly the same density as sideways-moving material, so the expansion speed is comparable in both directions. Therefore, in the time it takes for the jet’s height to increase from $a$ to $2a$, its width also increases from $b$ to $\sim b+a$. The volume at that time is roughly $V_{\rm c}\sim\frac{8}{3}\pi a^{3}$, so the cocoon pressure is $P_{\rm c}\sim E_{0}/8\pi a^{3}$. The time-scale $t_{a}$ is then given by $\sim a/v_{z}$, where $v_{z}\sim[(\gamma+1)P_{\rm c}/2\rho]^{1/2}$ is the on-axis expansion speed, leading to

[TABLE]

The fact that $t_{a}\sim(a/b)^{2}t_{0}$ implies that an outflow takes much longer to become quasi-spherical if it is initially narrow (i.e., if $b\ll a$). This is due to a combination of two effects which cause narrow jets to be more strongly decelerated. First, narrow jets undergo a larger drop in velocity when the jet ram pressure vanishes. After the jet is choked, the on-axis velocity decreases from its initial value, $\beta_{\rm h}$, to a new value, $\sim\beta_{\rm c}$, which is comparable to the sideways expansion speed. Since $a\sim c\beta_{\rm h}t_{0}$ and $b\sim c\beta_{\rm c}t_{0}$, the velocity is reduced by a factor of $\beta_{\rm c}/\beta_{\rm h}\sim b/a$. Second, narrow jets are more strongly affected by sideways expansion. As discussed above, the width of the outflow grows from $b$ to $\sim a$ in the time it takes the height to double. This increases the volume by a factor of $\sim(a/b)^{2}$, and decreases the pressure by a factor of $\sim(b/a)^{2}$. Since $v_{z}\propto P_{\rm c}^{1/2}$, the sideways expansion reduces the velocity by an additional factor of $(P_{\rm c}/P_{0})^{1/2}\sim b/a$, so by the time the shock reaches a height of $2a$, the velocity is only $\sim(b/a)\beta_{\rm c}\sim(b/a)^{2}\beta_{\rm h}$.

For $b\ll a$, the cocoon’s height $z_{\rm c}$ and width $r_{\rm c}$ in the three dynamical regimes evolve approximately as

[TABLE]

and

[TABLE]

Once the cocoon’s width has increased significantly (i.e., for $t\gg t_{b}$), the angle $\theta_{\rm b}$ between the axis and the point where the cocoon is widest is of order

[TABLE]

(For the derivation of equations 3–5, including the accurate determination of order-unity prefactors, see Section 3.) If $\alpha\leq 2$, the ratio $r_{\rm c}/z_{\rm c}$ approaches 1 as $t\rightarrow\infty$, and the angle $\theta_{\rm b}$ approaches $\pi/2$. Otherwise, $r_{\rm c}/z_{\rm c}\rightarrow[\sin(\pi/\alpha)]^{\alpha/(\alpha-2)}$ and $\theta_{\rm b}\rightarrow\pi/\alpha$.

In cases where the outflow becomes spherical, we find that the convergence to a sphere is rather gradual: as $t\rightarrow\infty$, the difference between the height and width goes to zero as (see Section 3)

[TABLE]

for $\alpha<2$, or as $(\ln z_{\rm c})^{-1}\propto(\ln t)^{-1}$ for $\alpha=2$. Because this quantity decreases so slowly, the difference $z_{\rm c}-r_{\rm c}$ can be used to constrain the radius at which the jet was choked, even at times considerably larger than $t_{a}$. If the both the height and width of the outflow are measured through observations, the choking radius can be estimated via

[TABLE]

Knowing the choking radius places strong constraints on the properties of the jet, as we discuss in Section 4. In particular, if the explosion energy and the density at $z=z_{\rm c}$ are also known, then the jet’s luminosity $L_{\rm j}$, opening angle $\theta_{\rm op}$ and duration $t_{\rm j}$ satisfy

[TABLE]

and

[TABLE]

where we have defined

[TABLE]

3 Analytical Solutions for Cocoon Evolution

As discussed above, the choking of a relativistic jet produces a narrow cocoon filled with hot gas. We suppose that the cocoon propagates in an ambient medium with a negligible pressure and a radially stratified density,

[TABLE]

The forward shock bounding the cocoon can be described by a curve $\theta=\theta(R,t)$ at time $t$, where $\theta$ and $R$ are the usual polar coordinates.

To help keep track of the many definitions employed throughout this paper, we provide a list of the symbols used in Appendix A, along with their meaning and the place in the text where they are defined. In naming variables, we make use of the following subscripts to simplify our notation:

• •

The subscript ‘0’ refers to to values measured at the choking time, $t=t_{0}$.

• •

Global properties of the cocoon (such as its height, volume, and pressure) which depend only on time are marked with the subscript ‘c’.

• •

Quantities pertaining to the ‘bulge’ (a special point on the cocoon surface where the instantaneous velocity is parallel to the equator; see Section 3.2) are denoted by the subscript ‘b’.

• •

Parameters characterizing the jet prior to choking are subscripted with a ‘j’.

• •

Finally, various constants that depend on the density profile are given the subscript ‘$\alpha$’.

Let $(R_{i},\theta_{i})$ be the set of coordinates describing the cocoon surface at the initial time $t=t_{0}$ when the jet was choked. We can then parametrize the initial shape of the forward shock by a curve $R_{i}(\theta_{i})$. Additionally, we define $\chi_{i}(\theta_{i})$ as the angle between the surface normal at $(R_{i},\theta_{i})$ and the axis, which is given by

[TABLE]

A schematic of the initial conditions is shown in Fig. 2. Particles starting with different values of $\theta_{i}$ each follow a unique trajectory as the outflow expands, such that any point along the shock surface can be associated with a particular value of $\theta_{i}$.

For concreteness, we take the initial shape of the cocoon to be a prolate ellipsoid of semi-major axis $a$ and semi-minor axis $b$, centred at the origin and oriented along the $z$-axis. In that case, geometry dictates

[TABLE]

and

[TABLE]

However, the model is general and applicable to arbitrary initial shapes, as long as 1) the shape has axisymmetry and reflective symmetry about the equator; 2) $R_{i}$111Hereafter, to reduce clutter, we refer to $R_{i}(\theta_{i})$ and $\chi_{i}(\theta_{i})$ as simply $R_{i}$ and $\theta_{i}$. These and all other quantities subscripted with ‘$i$’ should be understood as having an implicit dependence on $\theta_{i}$. decreases smoothly with $\theta_{i}$; and 3) there is no cusp on the axis (cusps at the equator are permissible since, as we will see, they form even when not present initially). The choice of shape mainly affects the early evolution of the outflow, and does not significantly impact our conclusions.

Now, suppose that at $t=t_{0}$ the gas within the cocoon has a pressure $P_{0}$, distributed uniformly over a volume $V_{0}$. At each point along the shock surface, the internal pressure applies a force in the direction normal to the surface. The speed $v$ of a surface element is set by balancing the mass and momentum flux in a frame comoving with the shock, leading to $v=[(\gamma+1)P^{\prime}/2\rho]^{1/2}$, where $P^{\prime}$ is the postshock gas pressure. The ambient pressure is assumed to be negligible compared to the ambient ram pressure. The speed in the direction normal to the surface can alternatively be found by direct differentiation of $\theta(R,t)$, yielding $v=|\partial\theta/\partial t|/|\nabla\theta|$. Equating the two expressions for $v$ results in a partial differential equation for $\theta(R,t)$,222Eq. 15 is undefined for a spherical flow since in that case $\partial\theta/\partial t\rightarrow\infty$ and $\partial\theta/\partial R\rightarrow\infty$. The equation is more well-behaved in the spherical limit when $\theta$ is chosen as the independent variable, i.e. when the shock is described by a curve $R=R(\theta,t)$. However, in the general case of an aspherical shock in a spherically symmetric density, it is more convenient to use $R$ as the independent variable, because the density is a function of $R$, not $\theta$.

[TABLE]

where we have assumed a strong, non-relativistic shock and applied the jump condition $\rho^{\prime}/\rho=v/v^{\prime}={(\gamma+1)/(\gamma-1)}$ for adiabatic index $\gamma$.

The KA allows us to replace the postshock pressure $P^{\prime}(t)$ in equation 15 by the volume-averaged pressure of the cocoon $P_{\rm c}(t)$. If losses are unimportant, the total cocoon energy, $E_{0}$, can be treated as constant, so that $P_{\rm c}(t)\propto E_{0}/V_{\rm c}(t)$. We can then write

[TABLE]

where $V_{\rm c}(t)$ is the volume of the cocoon at time $t$, and

[TABLE]

is an order-unity dimensionless parameter which depends on the ratio of the thermal energy $E_{\rm th}$ to the total energy $E_{0}$, as well as the ratio of the postshock pressure to the average pressure. The essential assumption of the KA is that $\lambda$ is constant in time. In this case the dynamics are governed solely by the change in pressure due to the expanding volume, and the evolution of the cocoon can be determined to a reasonable approximation without the need to calculate the structure of the postshock region in detail. For simplicity, we adopt $\lambda=1$ and $P^{\prime}=P_{\rm c}=E_{0}/(3V_{\rm c})$ in subsequent calculations. The actual value of $\lambda$ can be estimated from numerical simulations (see Section 5).

Equation 15 can be simplified by noting that $P^{\prime}$ depends only on $t$, while $\rho$ depends only on $R$. Following K92, we define a dimensionless parameter $x=x(t)$:

[TABLE]

The parameter $x$ serves as a dimensionless replacement for the time. It is initialized to $x=0$ at $t=t_{0}$ and thereafter increases monotonically with $t$. Using the fact that the initial shock velocity along the axis is $v_{0}=[(\gamma+1)P_{0}/2\rho_{0}]^{1/2}$, equation 18 can be rewritten as

[TABLE]

We then see that $x$ is essentially the product of two dimensionless ratios: the ratio $[P_{\rm c}(t)/P_{0}]^{1/2}$ comparing the pressure at time $t$ to the initial pressure, and the ratio $t/(a/v_{0})$ comparing the age to the initial lengthwise sound-crossing time, $a/v_{0}\sim\sqrt{t_{b}t_{a}}$.

It is possible for $x$ to be bounded or unbounded, depending on whether integral in eq. 18 converges or diverges as ${t\rightarrow\infty}$ (see also Section 3.4). Since we expect to recover the ST solution at late times, with $R\propto t^{2/(5-\alpha)}$ and $P_{\rm c}\propto R^{-3}\propto t^{-6/(5-\alpha)}$, we find that the integral diverges for $\alpha\leq 2$. The transition to the sideways expansion phase discussed in Section 2 occurs at $x\sim b/a$, while the transition to the quasi-spherical phase occurs when $x\sim 1$.

Substituting equations 11 and 18 into equation 15, we obtain

[TABLE]

which can be solved to find the cocoon shape as a function of $R$ and $x$ (see Sections 3.1 and 3.2). Once the shape is known, the volume can be calculated (see Section 3.3), and the pressure can be found via ${P_{\rm c}=P_{0}(V_{0}/V_{\rm c})}$. The time can then be determined by inverting equation 18 and integrating (see Section 3.4).

As we will see, the evolution of the cocoon depends strongly on $\alpha$, the power-law index of the external density profile. We begin by discussing the simple case of a constant ambient density as an illustrative example. Then, we examine how the evolution changes as the external density profile steepens.

3.1 A constant external density

When the external density is uniform with $\rho(R)=\rho_{0}$, each point along the shock surface travels with the same velocity, $v(t)=[(\gamma+1)P_{\rm c}(t)/2\rho_{0}]^{1/2}$, since all points are driven by the same uniform pressure into the same density. Furthermore, since there is no differential velocity between adjacent points on the surface, each small patch of surface maintains its original orientation. Therefore, a particle on the surface travels to infinity along a straight line in the direction normal to the surface. Examining equation 19, we see that it reduces to $x(t)=\int v(t){\rm d}t/a$ in the constant density case. In other words, when $\alpha=0$, the distance traveled along each trajectory at time $t$ is simply $ax(t)$. Thus, a particle beginning at $(R_{i},\theta_{i})$ at time $t_{0}$ travels a distance $ax$ along a straight line to reach a new position $(R,\theta)$ at time $t$, as illustrated in Fig. 3.

Because the trajectories are straight lines in this case, $R$ and $\theta$ can be determined simply from geometry. A triangle with side lengths $R$, $R_{i}$, and $ax$ is formed by connecting the points $(0,0)$, $(R_{i},\theta_{i})$, and $(R,\theta)$, as seen in Fig. 3. From the Law of Cosines,

[TABLE]

where we have defined $\omega_{i}$ as the angle between the initial direction of motion and the radial direction, i.e.

[TABLE]

The cylindrical coordinates $r$ and $z$ are obtained via ${r=R_{i}\sin\theta_{i}+ax\sin\chi_{i}}$ and ${z=R_{i}\cos\theta_{i}+ax\cos\chi_{i}}$. After some manipulation using angle sum and difference identities, ${\theta=\tan^{-1}(r/z)}$ becomes

[TABLE]

If $\theta_{i}$ is held fixed, equations 21 and 23 give the trajectory of a particle that started at $\theta_{i}$, while if $x$ is fixed, the equations describe the shape of the cocoon surface at the time $t=t(x)$. It is also possible to combine equations 21 and 23 to eliminate $x$, which gives $R$ as a function of $\theta$ along a trajectory stemming from $\theta_{i}$:

[TABLE]

As expected, this describes a straight line that passes through $(R_{i},\theta_{i})$ and forms an angle $\chi_{i}$ with the axis.

In order to understand the long term evolution of the outflow, it is informative to compare the expansion along the axis to the expansion perpendicular to the axis. We denote the height and width of the cocoon as $z_{\rm c}$ and $r_{\rm c}$, respectively. The height is defined as the radius at the pole,

[TABLE]

while the width is defined as the maximum distance between the shock and the axis, i.e.

[TABLE]

for a given $x$. In the constant-density case, since the trajectories are linear and the expansion velocity is equal in magnitude at each point on the surface, the width is ultimately set by the ejecta which have an initial velocity perpendicular to the axis, with $\chi_{i}=\pi/2$.

We now turn to the specific case of an initially ellipsoidal shape. For an ellipsoid with semi-major axis $a$ and semi-minor axis $b$, we have $(R_{i},\chi_{i})=(a,0)$ at $\theta_{i}=0$, and $(R_{i},\chi_{i})=(b,\pi/2)$ at $\theta_{i}=\pi/2$. Then, from equation 21, $z_{\rm c}=a(x+1)$ and $r_{\rm c}=a(x+b/a)$, so that the ratio of the cocoon’s width to its height is

[TABLE]

Defining an effective eccentricity, $\epsilon_{\rm eff}=\left(1-r_{\rm c}^{2}/z_{\rm c}^{2}\right)^{1/2}$, we have

[TABLE]

to leading order in $a/z_{\rm c}$.

As $x\rightarrow\infty$, the outflow becomes gradually wider and ${r_{\rm c}/z_{\rm c}\rightarrow 1}$ while $\epsilon_{\rm eff}\rightarrow 0$. In fact, we see from equation 21 that for large $x$, $R\approx ax$ regardless of $\theta_{i}$, indicating that the cocoon becomes spherical. Taking $x\propto R$ and $P_{\rm c}\propto V_{\rm c}\propto R^{-3}$ in equation 18, we see that the usual $R\propto t^{2/5}$ behaviour for a blast wave in a constant density is recovered.

The solution for $\alpha=0$ is plotted in Fig. 4 at the times when the cocoon’s height has reached 2, 4, 6, 8, and 10 times its initial value. The initial shape was taken to be an ellipsoid with aspect ratio $b/a=0.2$. Several particle paths are also indicated. A few key properties of the constant-density solution are worth emphasizing. First, all of the trajectories are straight lines that do not cross. There are no self-intersections of the cocoon surface. Second, the width of the outflow is always a maximum in the equatorial plane. Finally, the solution gradually approaches a self-similar spherical blast wave as $t$ becomes large. The process of becoming spherical is relatively rapid: by the time the cocoon’s height has doubled to $z_{\rm c}=2a$, eq. 27 shows that its width is already $r_{\rm c}>z_{\rm c}/2$, regardless of how narrow it was to start. As we will see, these basic features do not necessarily hold true at higher $\alpha$.

3.2 A power-law external density

In a non-uniform, spherically symmetric density profile, there is a density gradient in the radial direction that is not aligned with the surface normal. Although each small patch of the surface is driven from behind by a uniform pressure, adjacent points encounter different densities, feel different ram pressures, and therefore move with different speeds. As a result, elements of the surface change their orientation as they travel outwards, following trajectories that curve towards higher density. The strength of this effect depends on the misalignment of the density gradient and the surface normal, which is expressed by $|\chi_{i}-\theta_{i}|$. The cocoon becomes more deformed in regions where $|\chi_{i}-\theta_{i}|$ is larger, eventually developing a bulge near the location where $|\chi_{i}-\theta_{i}|$ is maximal.

Because the bulge forms at a larger radius where the density is lower, it expands sideways faster than material in the equatorial plane, and eventually overtakes it. Thus, unlike the constant-density case, in decreasing density profiles the cocoon’s width is ultimately set by the behaviour of the bulge, not the expansion near the equator. Once the bulge has formed, we can describe the cocoon’s shape with four quantities: the extent along the $z$-axis, $z_{\rm c}$; the width in the equatorial plane, $r_{\rm eq}$; and the coordinates $R_{\rm b}$ and $\theta_{\rm b}$ indicating the location of the bulge, which we define to be the point where the cocoon surface is locally parallel to the axis. This point is located at a height $z_{\rm b}=R_{\rm b}\cos\theta_{\rm b}$ above the equatorial plane and at a distance $r_{\rm b}=R_{\rm b}\sin\theta_{\rm b}$ from the axis. The cocoon width is then $r_{\rm c}=\text{max}(r_{\rm b},r_{\rm eq})$. This configuration is illustrated in Fig. 5.

The determination of $r_{\rm eq}$ comes with an important caveat. Because the points on the surface follow curved paths, self-intersections can occur in a non-uniform density profile: when any part of the cocoon reaches the equator, it encounters material coming from the opposite side due to the reflective symmetry of the system. Collisions at the equator are expected to generate a high-pressure region and drive a shock outwards in the equatorial plane. Our simple KA model does not take equatorial interactions into account, and therefore underestimates the true equatorial width. In addition, there is a cusp at the equator in the idealized model which is unlikely to exist in reality. These discrepancies are clearly seen in a comparison with numerical simulations (see Section 5). While the analytically-obtained value of $r_{\rm eq}$ is only accurate to within a factor of a few, it may none the less be useful as a lower bound.

For $\alpha\neq 2$, K92 showed that equation 20 can be solved by means of a conformal coordinate transformation. We define

[TABLE]

and then introduce new dimensionless coordinates ${\xi=(R/a)^{k_{\alpha}}}$ and $\phi=k_{\alpha}\theta$. Applying this transformation to equation 20 results in

[TABLE]

which has the same form as the Kompaneets equation 20 for a constant density medium. Essentially, this converts the problem of solving equation 20 for an arbitrary density power-law to the much easier problem of solving the constant-density case for different initial conditions. For convenience, we also define $\xi_{i}=(R_{i}/a)^{k_{\alpha}}$ and $\phi_{i}=k_{\alpha}\theta_{i}$.

The angle $\psi_{i}$ between the axis and the initial direction of motion in $\xi$–$\phi$ coordinates can be derived as follows. We first observe that the quantity ${\rm d}\ln R/{\rm d}\theta={\rm d}\ln\xi/{\rm d}\phi$ is preserved by the coordinate transformation. Then, we note that in order for the slope of a trajectory to be $\cot\chi_{i}$ at the point $(R_{i},\theta_{i})$, the initial condition ${\rm d}\ln R/{\rm d}\theta=\cot(\chi_{i}-\theta_{i})$ must be satisfied at $\theta=\theta_{i}$. By the same logic, ${\rm d}\ln\xi/{\rm d}\phi=\cot(\psi_{i}-\phi_{i})$ when evaluated at $\phi=\phi_{i}$. Thus, we have $\cot(\psi_{i}-\phi_{i})=\cot(\chi_{i}-\theta_{i})=\cot\omega_{i}$, and therefore $\psi_{i}=\phi_{i}+\omega_{i}$ (up to a phase that we choose to be zero).

Having determined the initial conditions in $\xi$–$\phi$ space, the solution of equation 30 now proceeds in the same fashion as in the $\alpha=0$ case. Following the procedure outlined in the previous section yields

[TABLE]

and

[TABLE]

Transforming back to $R$–$\theta$ coordinates, and using the fact that ${\psi_{i}-\phi_{i}=\omega_{i}}$ (see above) we obtain:

[TABLE]

[TABLE]

Just as before, we can combine the parametric equations 33 and 34 to eliminate $x$ and obtain an expression describing the trajectory of a particle that began at $\theta_{i}$:

[TABLE]

Note that when $R_{i}=a$, $\theta_{i}=0$, and $\omega_{i}=\chi_{i}$, equations 33–35 reduce to the expressions given in K92 for an on-axis point explosion.

Equation 35 can be used to derive parametric equations for the coordinates of the bulge. To do so, we first note that since $R_{\rm b}$ and $\theta_{\rm b}$ describe a particular point on the cocoon’s surface, they must be functions of time alone, i.e. $R_{\rm b}=R_{\rm b}(x)$ and $\theta_{\rm b}=\theta_{\rm b}(x)$. The trajectory passing through the point ($R_{\rm b}(x)$, $\theta_{\rm b}(x)$) at a given time $x$ can be traced back to a specific initial angular position $\theta_{\rm b0}(x)$, as shown in Fig. 6. That is to say, for each time $x$, there exists a corresponding initial angle $\theta_{\rm b0}(x)$ such that $R(x,\theta_{\rm b0}(x))=R_{\rm b}(x)$ and $\theta(x,\theta_{\rm b0}(x))=\theta_{\rm b}(x)$. The relationship between $x$ and $\theta_{\rm b0}$ is derived in Appendix B, but the exact functional form is not important for the present discussion. For now, it suffices to say that there is a one-to-one mapping between $\theta_{\rm b0}$ and $t$. Therefore, $\theta_{\rm b0}$ can be thought of as an alternative dimensionless time parameter, which has an initial value of $\theta_{\rm b0}=\pi/2$ and decreases monotonically to a final value $\theta_{\rm b0,min}$ as $t\rightarrow\infty$.

Now, since every trajectory meets the surface at a right angle, the trajectory which intersects the surface at ($R_{\rm b}(x),\theta_{\rm b}(x)$) must be parallel to the equator at that point (see Fig. 6). This means that $d\ln R/d\theta=\tan\theta_{\rm b}$ must be satisfied when evaluated at $R=R_{\rm b}$ and $\theta=\theta_{\rm b}$. Differentiating equation 35 and setting $\theta_{i}=\theta_{\rm b0}(x)$ leads to the equations

[TABLE]

[TABLE]

where we defined $R_{\rm b0}(x)\equiv R_{i}(\theta_{i}=\theta_{\rm b0}(x))$ and $\chi_{\rm b0}(x)\equiv\chi_{i}(\theta_{i}=\theta_{\rm b0}(x))$. As the cocoon expands, the position of the bulge traces out a path described by the parametric curve ($R_{\rm b}(\theta_{\rm b0}),\theta_{\rm b}(\theta_{\rm b0}))$.333From here on, we write $\theta_{\rm b0}(x)$ as simply $\theta_{\rm b0}$.

To understand the transition to spherical flow, our ultimate goal is to derive an expression akin to equation 27 which relates the width of the cocoon $r_{\rm c}$ to its height $z_{\rm c}$. To make a direct comparison between $r_{\rm c}$ and $z_{\rm c}$, we must first express them in terms of the same parameter. For example, in the constant density case, we derived expressions for both $r_{\rm c}$ and $z_{\rm c}$ in terms of $x$, which we then solved to obtain $r_{\rm c}$ as an explicit function of $z_{\rm c}$. In the general case, however, things are not so straightforward. While it is still possible to determine the height in terms of $x$ by setting $\theta_{i}=\chi_{i}=0$ and $R_{i}=a$ in equation 33,

[TABLE]

in general there is not a way to express $r_{\rm c}(x)$ explicitly in a similar fashion. The underlying reason that $z_{\rm c}$ can be described by a simple function of $x$, but $r_{\rm c}$ can not, is that the location where the height is greatest always occurs along the axis, while the place where the width is greatest changes its relative position on the surface as the cocoon evolves. Put another way, whereas the height is always greatest along a particular trajectory, the point where the width is greatest is associated with different trajectories at different times. (This also explains why an explicit form of $r_{\rm c}(x)$ exists in the constant-density case, since in that case the width is always determined by the ejecta traveling along a certain trajectory parallel to the equator.)

We thus conclude that $x$ is not, in general, a suitable choice of parameter for comparing $r_{\rm c}$ and $z_{\rm c}$ in a non-uniform density. What, then, would be a better choice? The discussion surrounding equations 36–37 suggests another option: parametrizing the system in terms of $\theta_{\rm b0}$. This turns out to be a good choice, in the sense that the important quantities characterizing the cocoon’s shape ($z_{\rm c}$, $r_{\rm b}$, and $z_{\rm b}$; see Fig. 5) can be written explicitly in terms of $\theta_{\rm b0}$ and the quantities $R_{\rm b0}$ and $\chi_{\rm b0}$ derived from the initial shape. For clarity of presentation, we reserve the derivation of the full system of parametric equations describing the shape for Appendix B.

The new parametrization allows a direct comparison between the height $z_{\rm c}(\theta_{\rm b0})$ and the width $r_{\rm c}(\theta_{\rm b0})$ of the cocoon at a given time $x$. However, a drawback compared to the constant-density case is that the equations cannot be solved to obtain an explicit function $r_{\rm c}(z_{\rm c})$ relating the height and width. None the less, it is possible to simplify the exact equations to obtain an approximate relation between $r_{\rm c}$ and $z_{\rm c}$ in each of the dynamical regimes discussed in Section 2 (see Appendix C). Expressions relating $r_{\rm eq}$ and $z_{\rm c}$ can also be derived in a similar manner (see Appendix D), although there are issues with the KA near the equatorial plane as discussed above.

Here, we give only a qualitative overview of the results, assuming an initially ellipsoidal shape as given by equation 13. In this case, the cocoon width is initially set by the ejecta in the equatorial plane, with $r_{\rm c}=r_{\rm eq}$. We find that the evolution of $r_{\rm c}$ depends on the relationship between $\alpha$ and the outflow’s initial aspect ratio, $a/b$. If $\alpha\leq 2b^{2}/a^{2}$, no bulge develops, and the outflow remains widest at the equator indefinitely, similar to the constant-density case. In this scenario, there are no trajectories which pass through the equator, so there is no interaction between ejecta from opposite sides of the equator. On the other hand, if $\alpha>2b^{2}/a^{2}$, the bulge eventually forms and overtakes the equatorial ejecta, and collisions at the equator occur as well.

We are mainly interested in the scenario where ${a/b\gg 1}$, in which case $\alpha\gg 2b^{2}/a^{2}$ also holds unless the density profile is very flat. In this limit, the bulge overtakes the equatorial material early in the planar phase, once the width has increased by $r_{\rm c}-b\approx(2b^{2}/\alpha a^{2})b\ll b$. The cocoon’s width is therefore governed by the behaviour of the bulge, with $r_{\rm c}=r_{\rm b}$ for the majority of the evolution. The motion of the bulge during each dynamical phase can be summarized as follows:

During the planar phase, the bulge initially remains close to the equator. At first, particles have not had sufficient time to significantly change their direction of motion, so the cocoon width is determined by the ejecta from $\theta_{\rm b0}\gg b/a$, which had an initial velocity almost perpendicular to the axis, with $\chi_{\rm b0}\approx\pi/2$. As the cocoon expands, material at larger $z$ (which started closer to the axis, but with a higher speed due to the lower density) begins to overtake material at smaller $z$, causing the bulge to march along the surface towards smaller $\theta$. By the time the cocoon’s width has doubled, the bulge has crossed most of the surface, with $\theta_{\rm b}$ decreasing from $\sim\pi/2$ to $\sim b/a$. 2. 2.

In the sideways expansion phase, the particles at the location of the bulge originated from ${b^{2}/a^{2}\ll\theta_{\rm b0}\ll b/a}$. On this part of the surface, the $r$- and $z$-components of the initial velocity are comparable, with $\chi_{\rm b0}\ga\pi/4$. The particles in the bulge have had to change their direction of motion by $\la 45$ degrees, moving farther from the axis in the process. By the time a particle starting from $\theta_{\rm b0}$ passes through the bulge, its angular position has increased to $\theta_{\rm b}\gg\theta_{\rm b0}$. Because the width of the outflow increases in this phase, while the height stays about the same (see Section 2), $\theta_{\rm b}$ now grows over time. This phase ends when the height and width of the cocoon become comparable, at which point $\theta_{\rm b}\sim 1$. 3. 3.

Finally, during the quasi-spherical phase, the width is governed by material that came from $\theta_{\rm b0}\ll b^{2}/a^{2}$. The material in the bulge was initially moving almost parallel to the axis, with $\chi_{\rm b0}\ll\pi/4$, and had its direction of motion changed by $\sim 90$ degrees. As in the previous phase, $\theta_{\rm b}\gg\theta_{\rm b0}$. As time goes on, the value of $\theta_{\rm b}$ continues to gradually increase. Eventually, the evolution becomes scale-free, with all length scales $\propto z_{\rm c}$. However, the outflow does not necessarily become spherical in the KA, as we discuss in Section 3.2.3. To capture the shape at infinity, we define two constants of proportionality depending on $\alpha$, $f_{\alpha}$ and $g_{\alpha}$, such that

[TABLE]

In addition, we define the maximum angle reached by the bulge as

[TABLE]

If the outflow becomes spherical, then ${f_{\alpha}=g_{\alpha}=1}$, and the bulge moves to the equator so $\theta_{\rm b\infty}=\pi/2$. Otherwise, we have $f_{\alpha}<1$ and $g_{\alpha}<1$, with $f_{\alpha}\neq g_{\alpha}$, and $\theta_{\rm b\infty}<\pi/2$.

In both the sideways expansion and quasi-spherical phases, we find that $\theta_{\rm b0}$ is negligible compared to $\chi_{\rm b0}$ and $\theta_{\rm b}$. This property is not unique to the bulge point. In fact, as long as $b\ll a$ and $\alpha\gg 2b^{2}/a^{2}$, the conditions $\theta_{i}\ll\theta$ and $\theta_{i}\ll\chi_{i}$ hold over all parts of the surface where $\theta\la x$. Particles with $\theta_{i}\ll 1$ also satisfy $R_{i}\approx a$, since they started out near the tip of the cocoon. Therefore, we can neglect $\theta_{i}$ and set $R_{i}/a\approx 1$ in equations 33 and 34 to obtain an approximate formula for the shape,

[TABLE]

which is valid for initially narrow cocoons out to angles of $\theta\la x$. During the sideways expansion phase, it can be shown that $x\sim\theta_{\rm b}$ (see Appendix B); therefore, in this regime, eq. 41 describes the shape between the axis and the bulge, for ${0\leq\theta\la\theta_{\rm b}}$. On the other hand, during the quasi-spherical phase $x\ga 1$ and eq. 41 is approximately valid over most of the cocoon surface.

The approximation used to generate equation 41 is equivalent to assuming that the shock was generated by a point explosion placed at $z=a$. Although K97 did not explicitly give an equation for $R(\theta,x)$ in the point explosion case, equation 41 can be reproduced by combining the two expressions given in their equation 23. As with other results derived from the KA, equation 41 is expected to be less accurate towards the equator.

While the evolution during the planar and sideways expansion phases is not affected much by the density gradient, the asymptotic behaviour of the outflow depends strongly on the density profile. Two different scenarios are possible, depending on whether the density profile is flat ($\alpha<2$) or steep ($\alpha>2$). We discuss each regime in turn below, along with the special case $\alpha=2$. As before, we treat only the specific case of an initially ellipsoidal shock.

3.2.1 A flat external density profile $(0<\alpha<2)$

For $\alpha<2$, $k$ is positive, and the integral in equation 18 diverges so $x\rightarrow\infty$ as $t\rightarrow\infty$. Inspecting equation 38, we see that $z_{\rm c}\propto x^{1/k_{\alpha}}\propto x^{2/(2-\alpha)}$ for large $x$. Assuming that the pressure drops as $P_{\rm c}\propto V_{\rm c}\propto z_{\rm c}^{-3}$ asymptotically, equation 18 implies $x\propto t^{(2-\alpha)/(5-\alpha)}$. Thus, as expected, we recover a blast wave evolution with $z_{\rm c}\propto t^{2/(5-\alpha)}$ at late times.

As discussed above, interactions at the equator occur in a non-uniform external density when $\alpha>2b^{2}/a^{2}$. Which regions on the cocoon surface eventually experience a collision? For $2b^{2}/a^{2}<\alpha<2$, it turns out that $\theta_{\rm b0}$ has a minimum value $\theta_{\rm b0,min}>0$. Trajectories stemming from $\theta_{\rm b0,min}$ are asymptotically parallel to the equator, i.e. they have $\theta\rightarrow\pi/2$ as $x\rightarrow\infty$. Particles originating from $\theta_{i}<\theta_{\rm b0,min}$ follow paths that never become parallel to the equator, so they can never be located at the bulge’s location. These particles stream freely to infinity, without interacting with any other ejecta. On the other hand, particles with $\theta_{i}>\theta_{\rm b0,min}$ eventually undergo a collision at the equator. The value of $\theta_{\rm b0,min}$ is derived in Appendix C.

As time goes on, $\theta_{\rm b0}$ decreases (see Fig. 6) and becomes progressively closer to $\theta_{\rm b0,min}$. To understand the sideways spreading at late times, we can therefore study the behaviour of equations 36 and 37 as $\theta_{\rm b0}\rightarrow\theta_{\rm b0,min}$. For $z_{\rm c}\gg a$, we find (see Appendix C) that the ratio of the cocoon’s width to its height is approximated by

[TABLE]

where

[TABLE]

and

[TABLE]

is an order-unity factor that depends on the density power-law. The effective eccentricity in this case is

[TABLE]

Since $r_{\rm c}/z_{\rm c}\rightarrow 1$ and $\theta_{\rm b}\rightarrow\pi/2$ as $z_{\rm c}\rightarrow\infty$, we infer that the outflow becomes spherical in the $0<\alpha<2$ case, and therefore $f_{\alpha}=g_{\alpha}=1$. However, the process of spherization is slower than in the constant-density case, with $r_{\rm c}/z_{\rm c}$ approaching unity as $\zeta^{-1}\propto z_{\rm c}^{-k_{\alpha}}$, rather than as $z_{\rm c}^{-1}$ as in equation 27.

The shape of the cocoon at $z/a=2$, 4, 6, 8, and 10 is shown in Fig. 7, for $\alpha=1$ and $b/a=0.2$. The trajectory stemming from $\theta_{\rm b0,min}$ and the evolution of the bulge’s location are also shown, respectively as dotted and dashed lines. To summarize, as in the constant density solution, the flow becomes more spherical with time, eventually becoming a self-similar ST blast wave. However, in contrast with the constant-density case, if $\alpha>2b^{2}/a^{2}$ a portion of cocoon material with $\theta_{i}>\theta_{\rm b0,min}$ eventually reaches the equator and experiences a collision. Furthermore, whereas for a constant density the cocoon is always widest at the equator, for a decreasing density profile, the cocoon bulges out and its width is maximum at some height $z_{b}>0$. But, because $z_{b}$ grows more slowly than $z_{\rm c}$, $z_{b}/z_{\rm c}$ steadily decreases and the relative position of the bulge gradually approaches the equator at late times.

3.2.2 A wind-like external density profile $(\alpha=2)$

When $\alpha=2$, equation 15 cannot be solved by the coordinate transformation used in Section 3.2.1, because $k_{\alpha}=0$. Instead, we go back to the original differential equation 20, which for $\alpha=2$ becomes

[TABLE]

This equation admits a solution of the form $\theta^{2}(x,R)=\theta_{x}(x)+\theta_{R}(R)$ (Kompaneets, 1960). By separation of variables, we find

[TABLE]

where $\kappa_{1}=\kappa_{1}(\theta_{i})$ and $\kappa_{2}=\kappa_{2}(\theta_{i})$ depend on the original shape. The initial conditions $\theta(x=0,\theta_{i})=\theta_{i}$, $R(x=0,\theta_{i})=R_{i}$, and ${\rm d}\ln R_{i}/{\rm d}\theta_{i}=-\tan\omega_{i}$ give $\kappa_{2}=\theta_{i}\cot\omega_{i}-\ln R_{i}$ and $\kappa_{1}^{2}=[\kappa_{2}^{2}+\theta_{i}^{2}\csc^{2}\omega_{i}]/2$. Rewriting, we have

[TABLE]

To derive an equation for the trajectories, we note that since the shock is a surface of constant $x$, $\nabla x$ is normal to the cocoon surface and parallel to the particle paths. Therefore, the trajectories satisfy

[TABLE]

Integrating subject to the initial condition $R=R_{i}$ at $\theta=\theta_{i}$ yields

[TABLE]

The $\alpha=2$ case is unique in the sense that $d\ln R/d\theta=\cot\omega_{i}$ is constant along particle paths.

As in Section 3.2, $R$ and $\theta$ can be written as parametric functions of $x$ and $\theta_{i}$. After substitution of equation 48 into 50, we eliminate either $\theta$ or $R$ to obtain, respectively,

[TABLE]

and

[TABLE]

As before, we obtain an approximate expression for the shape of the shock in the sideways expansion and quasi-spherical regimes by neglecting $\theta_{i}$ and adopting $R_{i}\approx a$. This gives

[TABLE]

Setting $\theta_{i}=\chi_{i}=0$ and $R_{i}=a$ in equation 51 gives the cocoon height,

[TABLE]

Repeating the procedure of Section 3.2, we find that the bulge is located at

[TABLE]

[TABLE]

Equations 50–56 agree with equations 33–37 and 41 in the limit of $\alpha\rightarrow 2$ and $k_{\alpha}\rightarrow 0$.

Overall, the late-time evolution for $\alpha=2$ is similar to the $\alpha<2$ case. One notable difference is that ${\theta_{\rm b0,min}=0}$ for $k_{\alpha}=0$, so $\theta_{i}>\theta_{\rm b0,min}$ holds everywhere except on the axis. Therefore, all of the cocoon material eventually reaches the equator and interacts. The asymptotic ratio of width to height in this case is

[TABLE]

and the eccentricity is

[TABLE]

Equations 57 and 58 can also be derived by taking the limit of equations 42–45 as $k_{\alpha}\rightarrow 0$. Although the cocoon once again becomes spherical, the logarithmic dependence on $z_{\rm c}$ makes the transition to sphericity relatively slow in a wind-like medium. As such, significant deviations from spherical symmetry are expected even after the cocoon has expanded to much larger than its original size.

Cocoon shapes and particle trajectories for $\alpha=2$ are shown in Fig. 8, using the same values as before for $z_{\rm c}$ and $b/a$. The dashed line, which tracks the position of the bulge in space, would become horizontal for a spherical flow. The extremely gradual change in the slope of this line captures the slow transition to sphericity.

3.2.3 A steep external density profile $(2<\alpha<3)$

When $2<\alpha<3$ ($k_{\alpha}<0$), the integral 18 converges, implying that $x$ has a finite value as $t\rightarrow\infty$. From equation 38, we see that $z_{\rm c}\propto(1+k_{\alpha}x)^{1/k_{\alpha}}$. Thus, in order to have $z_{\rm c}\rightarrow\infty$ as $t\rightarrow\infty$, we require $x\rightarrow 1/|k_{\alpha}|$. Evaluating both sides of equation 18 as $t\rightarrow\infty$ gives

[TABLE]

Then, subtracting equation 18 from equation 59 and multiplying through by $|k_{\alpha}|$, we obtain

[TABLE]

Now, if we now choose $t$ to be very large, we must have $P_{\rm c}\propto z_{\rm c}^{-3}\propto(1+k_{\alpha}x)^{-3/k_{\alpha}}$ since the flow will be self-similar. Plugging this into equation 60, we find ${(1+k_{\alpha}x)}\propto t^{(2-\alpha)/(5-\alpha)}$. Therefore the usual blast wave scaling $z_{\rm c}\propto t^{2/(5-\alpha)}$ is also recovered for $\alpha>2$.

A peculiar feature of the solution for $\alpha>2$ is that, while the flow eventually becomes self-similar, it does not become spherical. This is most easily seen from equation 36. Recall that as $t\rightarrow\infty$, we have ${\theta_{\rm b0}\rightarrow\theta_{\rm b0,min}}$ and $\theta_{\rm b}\rightarrow\theta_{\rm b\infty}$. Since $\theta_{\rm b0,min}=0$ in this case, we take $\theta_{\rm b0}\rightarrow 0$ in equation 36 to obtain $\theta_{\rm b\infty}=\pi/\alpha<\pi/2$, implying that the bulge stops short of reaching the equator. As we show in Appendix C, the width and the height in this case are related by

[TABLE]

where $\zeta$ is defined as before, and

[TABLE]

In this regime the parameter $A_{\alpha}$ takes on a different form,

[TABLE]

In the limit $\alpha\rightarrow 2$ and $k\rightarrow 0$, $A_{\alpha}\rightarrow\pi^{2}/8$ and $f_{\alpha}\rightarrow 1$, so we recover equation 57. The effective eccentricity is now given by

[TABLE]

Unlike the previous cases, $\epsilon_{\rm eff}$ no longer goes to zero, but instead approaches $(1-f_{\alpha}^{2})^{1/2}$ as $z_{\rm c}$ becomes large.

Considering equation 41 in the limit $x\rightarrow 1/|k_{\alpha}|$ leads to an asymptotic expression for the cocoon’s shape: as $z_{\rm c}$ tends to infinity, the curve bounding the cocoon is increasingly well-approximated by

[TABLE]

The overall shape of the shock front in the scale-free limit is the same as the shape that would be obtained from setting off two point explosions at $z=\pm a$ (see, e.g., K92).444Interestingly, the asymptotic shape often takes a familiar mathematical form when $k_{\alpha}$ is a rational fraction. For example, when $k_{\alpha}=-1/2$, the asymptotic shape is a cardioid. Setting $\theta=\pi/2$ in equation 65 and comparing with formula 39, we find

[TABLE]

Whereas $g_{\alpha}=f_{\alpha}=1$ for density profiles flatter than $r^{-2}$ because the cocoon ultimately becomes a sphere, steeper density profiles have $g_{\alpha}<f_{\alpha}<1$ instead. Therefore, $r_{\rm eq}<r_{\rm c}<z_{\rm c}$ at late times; rather than a sphere, the outflow becomes a self-similar ‘peanut.’

In Fig. 9, we plot the cocoon’s shape and several trajectories for $\alpha=3$. Compared to the previous cases, the trajectories are more curved. Also, because of the steep density gradient, the density near the tip is much lower than the density near the base. As a result, material originating near the axis wraps around and actually reaches the equator before material on the sides. The path traced out by the bulge (dashed line) straightens out and never becomes parallel to the equator, as expected for an outflow which does not become spherical.

We emphasize again that the above results are based on an idealized model that ignores heating due to collisions. The effects of equatorial interactions will tend to make the outflow more spherical than what the analytical model predicts. Interactions may be particularly important for steeper density profiles, because whereas in flat density profiles the equatorial collisions are glancing, in steep density profiles they can occur head-on (see Appendix D). This can be seen in Fig. 9: the first ejecta to reach the equator follow the dot-dashed path and collide fully head-on with material from the other side. The other paths in Fig. 9 also intersect the equator at a steep angle $\ga 45$ degrees. In comparison, the trajectories in the $\alpha=1$ case (Fig. 7) approach the equator at a much shallower angle and interact more weakly.

As a final note, we point out that for sufficiently steep $\alpha$, there can be a region inside the cocoon which is never shocked in the analytical model. This is simply a mathematical feature of the idealized KA model, which is not present in numerical simulations (see Section 4), and which most likely does not exist in reality. In a more realistic model taking interaction into account, we expect that a pressure gradient would develop due to collisions on either side of the unshocked region, which would then drive material into this region and shock it. In any case, because the unshocked region is engulfed by the rest of the cocoon, its presence or absence does not affect our conclusions about the late-time behaviour.

Gathering results obtained so far, it is possible to succinctly express the evolution of the cocoon’s width versus its height in each dynamical phase. In the planar phase, the width and height remain roughly constant. Once the width has doubled ($r_{\rm c}\gg 2b$), the velocity becomes comparable in the forwards and sideways directions, and the change in width ($r_{\rm c}-b\approx r_{\rm c}$) is about the same as the change in height ($z_{\rm c}-a$). This approximation remains valid throughout the sideways expansion phase, i.e. while $z_{\rm c}-a\ll a$. Finally, in the quasi-spherical phase, the evolution depends on $\alpha$, with the ratio of width to height is given by equation 42, 57, or 61, respectively for $\alpha<2$, $\alpha=2$, or $\alpha>2$. Combining these results, and assuming $b\ll a$, we then have

[TABLE]

Note that the intermediate regime $b\ll z_{\rm c}-a\ll a$ may not be realized if $b/a$ is not particularly small. In Table 1, we give the appropriate values of $f_{\alpha}$, $A_{\alpha}$, and $\zeta$ in the various regimes, and also summarize some important asymptotic characteristics of the cocoon.

Fig. 10 compares the limiting expressions in 67 with the full solution for each of the values of $\alpha$ discussed above. We see that when $\alpha$ is not too close to 2, the asymptotic approximation agrees well with the exact solution for $z_{\rm c}\ga 10a$. For $\alpha\approx 2$, however, the leading-order approximation given in equation 67 becomes less accurate. The reason for the slow convergence is that, since $\zeta\propto\ln z_{\rm c}$ for $\alpha=2$, the contribution from higher-order terms in $\zeta$ is non-negligible. If greater accuracy is required, the parametric equations presented in Appendix B can be used to calculate $r_{\rm c}/z_{\rm c}$ instead.

3.3 The volume of the cocoon

In order to determine the expansion speed of the cocoon, or the time that has elapsed since energy deposition, it is necessary to know the cocoon’s volume, $V_{\rm c}$. However, due to the complicated shape of the boundary, it is not always possible to obtain an analytical expression for the volume. To make the problem tractable, we consider the limiting behaviour of $V_{\rm c}$ for $t\ll t_{a}$ and $t\gg t_{a}$. Assuming the initial shape is an ellipsoid, the initial volume is $V_{0}=(4/3)\pi ab^{2}$. Subsequently, the volume grows as $V_{\rm c}\propto r_{\rm c}^{2}z_{\rm c}$. At early times ($t\ll t_{a}$), when the cocoon is still roughly elliptical, we approximate the volume as

[TABLE]

At late times ($t\gg t_{a}$), we have $r_{\rm c}\propto z_{\rm c}$ and therefore $V_{\rm c}\propto z_{\rm c}^{3}$. In this case we write

[TABLE]

where the constant of proportionality $C_{\alpha}$ is defined by

[TABLE]

For $\alpha\leq 2$, $C_{\alpha}=1$ since the outflow becomes spherical, while for steeper density profiles, $C_{\alpha}$ can be determined by integrating equation 65 over the enclosed volume. We then have

[TABLE]

The integral in equation 70 can be evaluated analytically in limited cases, but it is generally more convenient to evaluate it numerically. In Fig. 11, we plot $C_{\alpha}$ as a function of $\alpha$, and compare it to the value of $3V_{\rm c}/(4\pi z_{\rm c}^{3})$ found by direct numerical integration of equations 33 and 34 at $z_{\rm c}/a=10$, 100, and 1000. As expected from the results of Section 3.2, the $\alpha=2$ case is slowest to converge to the limiting volume, while for $\alpha\la 1$, the asymptotic expression is already a reasonably good approximation for $z_{\rm c}\ga 10a$. Additionally, we plot the estimate $C_{\alpha}\simeq f_{\alpha}^{2}$ obtained by treating the cocoon as an ellipsoid with volume $(4/3)\pi r_{\rm c}^{2}z_{\rm c}$. Replacing $C_{\alpha}$ by $f_{\alpha}^{2}$ is reasonably accurate, and has the advantage of being easy to compute via eq. 62.

3.4 The age of the cocoon

We now turn to the question of how $z_{\rm c}$ and $r_{\rm c}$ evolve in time. To start, we rewrite equation 18 as

[TABLE]

where we used $P(x)V(x)=P_{0}V_{0}$, $x(t_{0})=0$, $\gamma=4/3$, and $t_{0}\simeq t_{b}$ as given by equation 1. Now, we differentiate equation 38 to find $dx=(z_{\rm c}/a)^{k-1}d(z_{\rm c}/a)$ and change variables in the integral 72:

[TABLE]

Then, substituting for $V_{\rm c}(z_{\rm c})$ using equation 68 (for $z_{\rm c}-a\ll a$) or equation 69 (for $z_{\rm c}-a\gg a$), equation 73 becomes

[TABLE]

Finally, we replace $r_{\rm c}$ using equation 67 and perform the integration to obtain

[TABLE]

The leading constant term in the second expression was chosen to ensure a smooth transition at $z_{\rm c}-a=b$.

Inverting equation 75, we find the height as a function of time:

[TABLE]

As expected, the evolution becomes like an ST blast wave at late times with $z_{\rm c}\propto t^{2/(5-\alpha)}$. Finally, combining equations 67 and  76 gives an approximation for the cocoon width,

[TABLE]

In Fig. 12, we plot the analytical estimates in equations 76 and 77 respectively, compared to the accurate result obtained from numerical integration of the dynamical equations. We choose a small value of ${b/a=0.03}$ so that the intermediate regime shows up clearly. In general, the formulae in 76 and 77 do a good job at reproducing the full solution.

We conclude this section with an estimate for the time-scale for the cocoon expansion to become spherical. We define this time-scale, $t_{\rm sph}$, as the time when the width equals 90 per cent of the height. Note that $t_{\rm sph}$ is only defined for $\alpha\leq 2$, since for $\alpha>2$ the outflow does not become spherical (as discussed in Section 3.2.3). The spherization time can be written as the product of the time-scale for the height to double, $t_{a}$, and a scaling factor $Q_{\alpha}$ which depends on the density profile. Since $t_{a}\simeq\sqrt{2}(a/b)^{2}t_{0}$ (equations 1 and 2), we define $Q_{\alpha}$ via

[TABLE]

The numerically-computed value of $Q_{\alpha}$ as a function $\alpha$ is plotted in Fig. 13. For reference, we also show the times when the width becomes 60 and 80 per cent of the height, as well as an analytical estimate for $Q_{\alpha}$ (see below). The value of $Q_{\alpha}$ grows extremely fast with $\alpha$, with $Q_{\alpha}=59$, 441, and $4.4\times 10^{6}$ respectively for ${\alpha=0}$, 1, and 2. Outflows in steep density profiles therefore take a considerably longer time to become spherical.

To obtain an analytical estimate for $Q_{\alpha}$, we start from equation 67 and observe that $r_{\rm c}/z_{\rm c}=0.9$ holds when ${\zeta\simeq 10A_{\alpha}}$. Using the definition of $\zeta$ (equation 43), we then determine that this occurs at ${z_{\rm c}\simeq a(1+10k_{\alpha}A_{\alpha})^{1/k_{\alpha}}}$ when $\alpha<2$, or $z_{\rm c}\simeq a{\rm e}^{5\pi^{2}/4}$ when $\alpha=2$. Finally, we use equation 76 to translate the limit on $z_{\rm c}$ into a limit on the age, arriving at

[TABLE]

This approximation overestimates the true value of $Q_{\alpha}$, because equation 67 underestimates the true value of $r_{\rm c}/z_{\rm c}$ at late times. For $\alpha\la 1.5$, the analytical estimate of $Q_{\alpha}$ is roughly twice the numerically-determined value. However, the approximation is less accurate when $\alpha$ is close to 2, because the series expansion of $r_{\rm c}(z_{\rm c})$ converges slowly in that case, and we have used only the lowest-order term. For $\alpha=2$, equation 79 overestimates $Q_{\alpha}$ by a factor of 9. None the less, Fig. 13 and equation 79 both clearly demonstrate that $t_{\rm sph}$ is much larger than $t_{a}$, especially for $\alpha\sim 2$. The long-lingering asphericity provides a useful means to constrain the jet properties, as we discuss in Section 4.

3.5 Eventual breakdown of the solution

In this section, we estimate the time at which the KA solution breaks down. The KA model makes two key assumptions which may eventually become invalid. First, the shock is assumed to be strong. In a static medium, this condition holds until the outflow reaches pressure equilibrium with the surrounding medium, on a time-scale $t_{\rm eq}$. In an expanding medium, on the other hand, the shock becomes weak once the shock velocity becomes comparable to the expansion velocity of the ambient gas, which takes a time $t_{\rm w}$. Second, the model assumes that the ambient medium extends to infinity. If the extent is finite, then the solution is valid only up until the time $t_{\rm bo}$ when the shock breaks out. The KA model applies up until $t=\min(t_{\rm eq},t_{\rm w},t_{\rm bo})$.

To make for an easy comparison, we will write each of the time-scales as a product of $t_{a}$ and a dimensionless factor $\eta$. Also, we neglect order-unity constants which depend on the density profile. We first consider $t_{\rm eq}$. Pressure equilibrium is attained once the cocoon pressure, $P_{\rm c}$, becomes comparable to the ambient pressure, $P_{a}$. We suppose that the ambient pressure obeys

[TABLE]

where $\gamma_{a}$ is the adiabatic index of the ambient medium.555We note that pressure equilibrium can only be achieved if the cocoon pressure falls off more steeply than the ambient pressure; this is true only for $\alpha<3/\gamma_{a}$. Furthermore, we assume that pressure upon choking satisfies $P_{0}\gg P_{a0}$ (otherwise, the solution was not valid in the first place). In this case the shock remains strong throughout the planar phase.

At the time $t_{a}$, when the outflow transitions from sideways expansion to quasi-spherical expansion, the cocoon pressure is $P_{\rm c}(t_{a})\sim E_{0}/(4\pi a^{3})$. Therefore, there are two possible scenarios: if $P_{a0}\ga P_{\rm c}(t_{a})$, pressure equilibrium is reached during the sideways expansion phase, while if $P_{a0}\la P_{\rm c}(t_{a})$, it happens in the quasi-spherical regime. In the former case, the cocoon’s height is given by $z_{\rm c}\sim a$ and its volume is $V_{\rm c}\sim 4\pi ar_{\rm c}^{2}/3$. The expansion occurs near the tip of the cocoon, where $P_{a}\sim P_{a0}$. Pressure equilibrium is attained when the width is $r_{\rm c}\sim a(E_{0}/4\pi P_{a0}a^{3})^{1/2}$. In the latter case, the outflow is roughly spherical, and pressure balance is achieved when the shock’s size satisfies $z_{\rm c}\sim a(E_{0}/4\pi P_{a}(z_{\rm c})a^{3})^{1/3}$. Substituting equation 80 for $P_{\rm c}(z_{\rm c})$ then leads to $z_{\rm c}\sim a(E_{0}/4\pi P_{a0}a^{3})^{1/(3-\gamma_{a}\alpha)}$. As a final step, we note that age of the system in the sideways expansion phase ($t_{b}\ll t\ll t_{a}$) is given by $t\sim t_{0}(r_{\rm c}/b)^{2}\sim t_{a}(r_{\rm c}/a)^{2}$ via equation 77, and the age in the quasi-spherical phase ($t\gg t_{a}$) is given by $t\sim t_{a}(z_{\rm c}/a)^{(5-\alpha)/2}$ via equation 76. We then arrive at the pressure equilibrium time-scale

[TABLE]

where

[TABLE]

The condition $\eta_{\rm eq}\gg(b/a)^{2}$ comes from our assumption that $P_{0}\gg P_{a0}$.

The above analysis assumes a static medium. An alternative possibility is that the cocoon propagates in a supersonically expanding medium with an expansion velocity $v_{\rm w}$ exceeding the ambient sound speed. Then, once $v\sim v_{\rm w}$, the shock becomes weak and our solution no longer applies. Let us assume a steady wind with $\rho\propto r^{-2}$ and $v_{\rm w}=const$. We also require that $v\gg v_{\rm w}$ initially. Setting $v\sim(P_{\rm c}/\rho)^{1/2}\sim v_{\rm w}$, the calculation then proceeds as before. We ultimately arrive at

[TABLE]

with $\eta_{\rm w}$ given by

[TABLE]

where we used the relation $\dot{M}=4\pi R^{2}v_{\rm w}\rho(R)$ for a steady wind with a mass loss rate of $\dot{M}$.

Finally, we consider the break out of the shock from a density profile extending out to a radius of $R_{*}$. The breakout time is sensitive to location where the jet is choked relative to $R_{*}$. It is important to note that, due to the rapid sideways expansion and deceleration of a choked outflow, failed jets take considerably longer to break out than successful jets, unless the choking occurs close to the edge. We suppose that the jet is choked sufficiently far from the edge (i.e., $R_{*}-a\gg b$), so that the breakout takes place in either the sideways expansion regime (if $b\ll R_{*}-a\ll a$) or in the quasi-spherical regime (if $R_{*}-a\gg a$). Then, taking $z_{\rm c}=R_{*}$ in equation 75, the breakout time is

[TABLE]

where

[TABLE]

The eventual fate of the outflow depends on the ordering of the time-scales discussed above. We investigate three astrophysically relevant scenarios: a GRB jet choked inside a WR star (i.e., a failed long GRB); a GRB jet choked within an extended, low-mass stellar envelope (as may be the case in $ll$GRBs); and an AGN jet in a cluster environment choked by the intracluster medium. In each case, the ambient medium can be treated as static. Therefore, the relevant time-scales are $t_{a}$, $t_{\rm eq}$, and $t_{\rm bo}$.

In a typical long GRB, it takes the jet $\sim 15$ seconds to drill through the progenitor WR star (Bromberg et al., 2012). Therefore, we expect jets with a duration of $\sim$ a few seconds to be choked at an appreciable fraction of $R_{*}$, i.e. we have $a/R_{*}\sim$ a few tenths. In this case $\eta_{\rm bo}$ is order-unity and $t_{\rm bo}\sim t_{a}$. A WR star can be approximated by an $n=3$ polytrope, with the ambient density and pressure scaling as $\rho\propto R^{-5/2}$ and ${P_{a}\propto\rho^{4/3}\propto R^{-10/3}}$ in the outer parts of the star. As $P_{a}R^{3}\propto R^{-1/3}$ changes slowly with radius, the quantity $P_{a0}a^{3}$ appearing in equation 82 is roughly constant, regardless of the choking location. Numerical modelling of pre-supernova stars (e.g., Heger et al., 2005) suggests a pressure on the order of $2\times 10^{18}\,\text{erg cm}^{-3}$ at $R=10^{10}$ cm. We then have $\eta_{\rm eq}\simeq 80(E_{0}/10^{51}\,\text{erg})(P_{a0}a^{3}/10^{48}\,\text{erg})^{-1}\gg 1$. Thus, $t_{a}\sim t_{\rm bo}\ll t_{\rm eq}$, and the breakout occurs while the KA model is still valid.

The progenitors of $ll$GRBs, on the other hand, may have an extended optically thick envelope containing a mass $\sim 10^{-2}M_{\odot}$ within a radius $\sim 100-1000R_{\odot}$ (Campana et al., 2006; Nakar, 2015; Irwin & Chevalier, 2016). The pressure in this envelope is unknown, but it must be smaller than in the WR star case, since the envelope is both cooler and less dense. We then expect $\eta_{\rm eq}\ga 100$, as before. However, because of the larger radius compared to the WR star case, a typical GRB jet is choked deep within the extended envelope, with $a/R_{*}\ll 1$. We therefore have $\eta_{\rm bo}\gg 1$ and a quasi-spherical breakout occurs at $t_{\rm bo}\gg t_{a}$. Assuming that the jet successfully escapes the star, $a$ is at least $\sim 10R_{\odot}$, so $\eta_{\rm bo}$ is at most $\sim 100$. We then find $t_{a}\ll t_{\rm bo}<t_{\rm eq}$, so in this case the KA model also applies up until breakout.

In AGNs, bubbles with a size of $\sim 10$ kpc are sometimes observed in the centre of galaxy clusters (e.g., Bîrzan et al., 2004; Diehl et al., 2008). If these bubbles were inflated by jets, the jets must have been choked at $a<10$ kpc. However, the dense cores in clusters extend to much larger radii of ${R_{*}\sim 100}$ kpc, so we have $a/R_{*}\la 0.1$ and $\eta_{\rm bo}\ga 10$. Thus, like in the $ll$GRB case, the breakout time is $t_{\rm bo}\gg t_{a}$. The difference is that, compared to the GRB case, the ambient pressure in an AGN environment is much more significant. Taking a typical jet energy $E_{0}\sim 10^{59}$  erg and a typical ambient pressure of $P_{a}(10\,\text{kpc})\sim 0.1$ keV cm*-3* at $R=10$ kpc, we obtain $\eta_{\rm eq}\ga 1.6(E_{0}/10^{59}\,\text{erg})[P_{a}(10\,kpc)/0.1\,\text{keV\,cm}^{-3}]^{-1}$ for $a\la 10$ kpc. Assuming a flat ($\alpha\simeq 1$) and isothermal ($\gamma_{a}=1$) ambient medium, the lower limits on $\eta_{\rm bo}$ and $\eta_{\rm eq}$ imply $t_{\rm bo}\ga 100t_{a}$ and $t_{\rm eq}\ga 2t_{a}$, respectively. Furthermore, in the isothermal case, the ratio $t_{\rm bo}/t_{\rm eq}$ is independent of $a$. Dividing equation 85 by equation 81 and applying $P_{a0}=P(10\,\text{kpc})(a/10\,\text{kpc})^{-1}$, we obtain

[TABLE]

Therefore, AGN jets choked in the core of galaxy clusters satisfy $t_{a}\la t_{\rm eq}\ll t_{\rm bo}$. Pressure equilibrium is achieved early in the quasi-spherical phase, and the KA approximation breaks down long before the flow crosses the cluster core.

In AGNs, buoyancy will also become important at late times. The effects of buoyancy are discussed in more detail in a separate paper (Irwin et al., in prep.), where we apply the model to young bubbles inflated by jets in galaxy clusters. Here, we simply remark that if the ambient medium is in hydrostatic balance, and the size of the bubbles is comparable to their distance from the central source, the buoyancy time-scale $t_{buoy}$ is comparable to $t_{\rm eq}$ (Churazov et al., 2001, see also Irwin et al., in prep).

4 Inferring the jet parameters from the cocoon properties

For a cocoon formed by a choked jet expanding with a non-relativistic velocity, the initial conditions of the cocoon ($a$, $b$, and $E_{0}$) are related to the parameters of the jet (the luminosity $L_{\rm j}$, the duration $t_{\rm j}$, and the initial opening angle $\theta_{\rm op}$) in a straightforward way. In this section, we consider how to work backwards from the observed characteristics of the cocoon to obtain useful constraints on the jet properties. We will assume that it is possible to estimate the cocoon’s height ($z_{\rm c}$), width ($r_{\rm c}$), and energy ($E_{0}$) from observations, and additionally that the density profile is known. The calculation then proceeds in two steps. First, we use the KA model to estimate the age of the system and determine the height and width upon choking ($a$ and $b$). Then, we apply the jet propagation model of Bromberg et al. (2011) to find the jet properties needed to match this initial shape. Throughout this discussion, we ignore order-unity prefactors for simplicity.

We caution that our analysis applies specifically to outflows which are initially narrow with $b\ll a$, as expected when the cocoon is inflated by a relativistic jet. When the observed aspect ratio is $z_{\rm c}/r_{\rm c}\gg 1$, the association with a jet is not controversial. However, when $z_{\rm c}/r_{\rm c}$ is of order unity, the alternative possibility that the outflow was powered by a wide-angle wind cannot be ruled out.

The first step is to determine which dynamical phase is appropriate. Roughly speaking, the outflow is in the quasi-spherical phase when $r_{\rm c}/z_{\rm c}\ga 0.5$, and in the planar or sideways expansion phase when $r_{\rm c}/z_{\rm c}\la 0.5$. From equations 77 and 76, along with the definitions of $t_{b}$ and $t_{a}$, we deduce that the age in each case is approximately

[TABLE]

where $\rho(z_{\rm c})$ is the observed density at $z_{\rm c}$.

Next, we wish to estimate $a$ and $b$. Because the outflow is slow to become spherical (see Section 3.4), the dependence of the aspect ratio on $a$ persists until well after $t_{a}$. The observed value of $r_{\rm c}/z_{\rm c}$ is then directly related to the ratio $z_{\rm c}/a$, as shown in Fig. 10. In the quasi-spherical case there is also a significant dependence on the density profile. An analytical approximation for $a$ is obtained by inverting equation 67:

[TABLE]

For elongated outflows in the planar or sideways expansion phase, we find that the choking occurred closer to the edge of the outflow than the centre, while in quasi-spherical outflows the opposite holds true.

The initial width, however, is more difficult to determine. In the quasi-spherical regime, the shape has little to no dependence on the initial width, so $b$ cannot be recovered. In the sideways expansion phase, there is a weak dependence on $b$, but there is also another issue, which is that the sideways expansion and planar cases are difficult to distinguish observationally. One indication that the cocoon is in the sideways expansion phase is a bulging out towards the tip, as seen in Fig. 1. However, without knowing the initial shape, this method is unreliable, and moreover it only applies when the density profile is fairly steep. With good observations, a better way to test which case is relevant is to compare the strength of the shock along the axis to the strength near the widest part of the cocoon. If the shock is significantly stronger towards the axis, it would be good evidence for a recently quenched jet, suggesting a planar evolution with $t\sim t_{0}$ and $b\sim r_{\rm c}$.

Statistically speaking, though, it is much more likely to catch the system during the sideways expansion phase, since it lasts roughly $(a/b)^{2}$ times longer than the planar phase. In this case it is hard to tell if the jet was choked somewhat recently, with a width comparable to the observed width, or if it was choked farther in the past, but was initially narrower. All we can say for certain is that $b<r_{\rm c}$ and $t_{0}<t$. For a known choking radius $a$ given by equation 89, the degeneracy between the initial width and the choking time is captured by

[TABLE]

where in the second equality we used $\rho_{0}a^{\alpha}=\rho(z_{\rm c})z_{\rm c}^{\alpha}$ to replace $\rho_{0}$ by observable quantities.

Now that we have established how $a$ and $b$ depend on the observables, let us consider how they relate to the jet properties. Suppose the jet is injected relativistically with luminosity $L_{\rm j}$ and opening angle $\theta_{\rm op}$, and operates for a duration $t_{\rm j}$. Bromberg et al. (2011) showed that the dynamics of the jet are governed by a dimensionless quantity $\tilde{L}$. At the time $t=t_{0}$ when the jet is choked, $\tilde{L}=\text{max}\{(\frac{L_{\rm j}}{\rho_{0}t_{0}^{2}\theta_{\rm op}^{4}c^{5}})^{2/5},\frac{L_{\rm j}}{\rho_{0}t_{0}^{2}\theta_{\rm op}^{2}c^{5}}\}$. The jet is collimated with a non-relativistic head when $\tilde{L}\ll 1$; collimated with a relativistic head when $1\ll\tilde{L}\ll\theta_{\rm op}^{-4/3}$; and uncollimated with a relativistic head when $\tilde{L}\gg\theta_{\rm op}^{-4/3}$. $\tilde{L}$ is closely related to the head velocity $\beta_{\rm h}$, with $\beta_{\rm h}\approx\tilde{L}^{1/2}$ for a non-relativistic jet, and $2\Gamma_{h}^{2}=\tilde{L}^{1/2}$ for a relativistic head, where $\Gamma_{h}\approx(1-\beta_{\rm h}^{2})^{-1/2}$.

The case of an uncollimated jet can be further subdivided depending on the relative values of $\tilde{L}$ and $\theta_{\rm op}^{-4}$. If $\tilde{L}\ll 4\theta_{\rm op}^{-4}$, the jet and cocoon are in causal contact and the pressure in the cocoon remains more or less uniform. However, if $\tilde{L}\gg 4\theta_{\rm op}^{-4}$, causal connection is lost and the approximation of uniform pressure breaks down. We therefore restrict our discussion to the $\tilde{L}\ll 4\theta_{\rm op}^{-4}$ case in which the KA remains valid.

After the jet is switched off at $t_{\rm j}$, how long does it take to be choked? Material continues to flow into the cocoon until $t_{0}$, when the last bit of material launched by the jet catches up to the head. For a non-relativistic head, this time is negligible compared to $t_{\rm j}$, and therefore $t_{0}\simeq t_{\rm j}$. If the head is relativistic instead, choking takes longer, and $t_{0}=t_{\rm j}/(1-\beta_{\rm h})\approx 2\Gamma_{h}^{2}t_{\rm j}$ (Nakar, 2015). Thus we have

[TABLE]

The distance the head travels in time $t_{0}$ is $\beta_{\rm h}ct_{\rm j}$ in the former case, and $2\Gamma_{h}^{2}ct_{\rm j}$ in the latter; in either regime, we can write

[TABLE]

Then, estimating the cocoon width at the time of choking as ${b\simeq(\beta_{\rm c}/\beta_{\rm h})a}$, where the sideways expansion speed is $\beta_{\rm c}\simeq\tilde{L}^{1/2}\theta_{\rm op}$ ($\tilde{L}\ll\theta_{\rm op}^{-4/3}$), or $\beta_{\rm c}\simeq\tilde{L}^{1/8}\theta_{\rm op}^{1/2}$ ($\theta_{\rm op}^{-4/3}\ll\tilde{L}\ll 4\theta_{\rm op}^{-4}$)(Bromberg et al., 2011), we arrive at

[TABLE]

We now have three equations (91, 92, and 93) relating the initial conditions of the cocoon ($a$, $b$, and $t_{0}$) to three quantities describing the jet ($\tilde{L}$, $t_{\rm j}$, and $\theta_{\rm op}$). When the choking radius satisfies $a\ll ct_{0}$, we find that the head was non-relativistic, with $\tilde{L}$ given by

[TABLE]

In this case, the three jet parameters are uniquely determined by the cocoon properties:

[TABLE]

When the jet head is relativistic, however, a degeneracy arises because $a$ and $t_{0}$ are related trivially by $a=ct_{0}$. In this case only two of the cocoon parameters are independent and $\tilde{L}$ is not uniquely determined by the cocoon properties. Instead, $\tilde{L}$ and $\theta_{\rm op}$ follow a closure relation given by

[TABLE]

None the less, the relationship between the jet parameters can still be constrained. In either case, we have the same constraint on the total injected energy as before, i.e. $L_{\rm j}t_{\rm j}\simeq E_{0}$. For a collimated jet with a relativistic head ($1\ll\tilde{L}\ll\theta^{-4/3}$), we obtain the additional constraint

[TABLE]

while for an uncollimated jet in causal contact with its cocoon ($\theta^{-4/3}\ll\tilde{L}\ll 4\theta^{-4}$), we find

[TABLE]

In both GRBs and AGNs, the usual case is that the jet head is Newtonian (see, for example, the discussion in Bromberg et al., 2011). The jet properties are then given by equation 95. However, a precise calculation requires knowledge of the initial width, which may be difficult to extract from observations for the reasons discussed above. In cases where $E_{0}$ and $a$ are known, but $b$ is unconstrained, we can only learn about the quantities $L_{\rm j}t_{\rm j}$ and $t_{\rm j}\theta_{\rm op}^{-2}$, which are independent of $b$. From equation 95 we have $L_{\rm j}t_{\rm j}=E_{0}$ (as expected), and $t_{\rm j}\theta_{\rm op}^{-2}\simeq(a/b)^{2}t_{0}\simeq t_{a}$. We then find that $t_{\rm j}\theta_{\rm op}^{-2}$ satisfies a relation resembling equation 90:

[TABLE]

The reason this degeneracy comes about is that the choking radius scales as $a\propto\tilde{L}^{1/2}t_{\rm j}$ via equation 92. In the non-relativistic regime, this becomes $a\propto E_{0}^{1/5}(t_{\rm j}/\theta_{\rm op}^{2})^{1/5}$. Therefore, two jets with the same energy and the same $t_{\rm j}\theta_{\rm op}^{-2}$ are choked at the same location. Even if the two jets had different widths upon choking, after several $t_{0}$ the initial width becomes unimportant, and the resulting outflows look about the same.

To conclude this section, we consider how the jet properties influence the late-time evolution. The characteristic time-scale for the cocoon’s height to double is related to the jet parameters by

[TABLE]

Unsurprisingly, $t_{a}$ is larger for narrow or long-lived jets. Interestingly, however, for jets of a given opening angle and duration, $t_{a}$ has a minimum with respect to $\tilde{L}$, with $\min(t_{a})\simeq\theta_{\rm op}^{-4/3}t_{\rm j}$ occurring at $\tilde{L}\simeq\theta_{\rm op}^{-4/3}$. This suggests that barely collimated jets become spherical on the shortest time-scale compared to the jet working time, which makes sense given that the time to become spherical scales as $t_{\rm sph}\propto(a/b)^{2}t_{0}$. As $\tilde{L}$ increases, $a/b$ decreases while $t_{0}$ grows. A tightly collimated jet is easy to suffocate, but takes a longer time to spherize because it leaves behind a narrow cocoon. On the other hand, a powerful uncollimated jet produces a wider cocoon, but takes a longer time to choke in the first place because the jet head is highly relativistic. Mildly collimated jets balance these two effects to minimize the spherization time. The radius where the flow becomes spherical, however, is $\propto a\propto\tilde{L}^{1/2}t_{\rm j}$ and is always larger for more powerful jets.

5 Comparison with Numerical Simulations

To check the validity of the analytical results, we compared them to numerical simulations carried out by a collaborator, O. Gottlieb, using the publicly available hydrodynamics code PLUTO (Mignone et al., 2007). The axisymmetric simulation setup is akin to previous studies of GRB jet propagation (e.g., Gottlieb et al., 2018; Harrison et al., 2018). The jet is injected with a two-sided luminosity of $10^{51}$ erg s*-1* and a Lorentz factor of 5 into a nozzle of width $8\times 10^{7}$ cm. An adiabatic index of 4/3 is used in all cases. The lower boundary of the simulation is placed at $z_{inj}=10^{9}$ cm, and the external density is defined by $\rho(R)=6\times 10^{4}(R/z_{inj})^{-\alpha}$ g cm*-3*. The ambient pressure is set to a negligibly low value.

We compared our results with three simulations, for three values of $\alpha$ relevant to astrophysical jets: $\alpha=1$, as is typical for the dense core of a galaxy cluster; $\alpha=2$, as expected for a dense circumstellar wind; and $\alpha=2.5$, as appropriate in the outer parts of a WR star. Tracer particles are included in the unshocked jet to determine when all of the jet material has flowed into the cocoon. Once all the jet has been shocked, the maximum speed of the tracers drops. In each simulation, we locate the time step where this drop occurs and call this the choking time. We then find the height, width, and volume-averaged pressure of the outflow at the choking time, and use these as initial conditions for the analytical model. From there we evolve the analytical cocoon model according to Section 3, assuming a nominal value of $\lambda=1$ in equation 16.

The simulation results are shown in figures 14, 15, and 16, respectively for $\alpha=1$, 2, and 2.5. The analytical solution obtained from the KA at each time is also shown in a heavy black line. In the $\alpha=1$ and $\alpha=2$ cases, there is good agreement (within 10–15 per cent) between the analytical and numerical results out to about 45 degrees from the axis. The slight disagreement is due to the fact that the pressure behind the shock in the simulations is slightly larger than assumed by the KA. The shapes can be made to agree by adjusting the parameter $\lambda$, with $\lambda\approx 1.2$ giving an improved fit.

Beyond 45 degrees, the shapes start to deviate, and as expected the agreement becomes worse near the equator. Typically, the analytical model underestimates the equatorial radius by a factor of 2–3. As discussed above, one possible source of the discrepancy at the equator is the increased pressure due to collisions. However, another reason that the equatorial radius is larger in the simulations is that the analytical model assumes, following the KA, that the initial pressure is uniform. In the simulations, this is not the case; rather, there is a clear gradient in pressure and velocity at the time of choking, with both being larger towards the equator. One possible reason why the shock might be stronger at large angles is that the ejecta near the equator were shocked at earlier times when the sideways expansion was faster and the density was higher. An alternative possibility is that some interaction at the equator already took place while the jet was active.

In the $\alpha=2.5$ case, there are also discrepancies near the equator, for the same reasons discussed above. In addition, the extent along the axis is somewhat larger in the numerical simulation compared to before. This comes about because the analytical model assumes that the outflow has already decelerated to a steady state where the velocity is given by $\sim(P/\rho)^{1/2}$, or in other words, that the ram pressure of the jet suddenly becomes unimportant at the choking time. It seems that this assumption is reasonable when $\alpha=1$ or 2, but not for $\alpha=2.5$. This makes sense considering that, prior to choking, the jet head is accelerating in the $\alpha=2.5$ case, whereas in the other cases it is is decelerating ($\alpha=1$) or has a constant velocity ($\alpha=2$). Therefore, for $\alpha=2.5$, it takes longer for the on-axis ejecta to decelerate and approach the KA prediction. In fact, a shock is clearly visible in Fig. 16, indicating that deceleration is ongoing. Despite the differences on the axis and at the equator, the analytical solution still does a good job at predicting the cocoon’s width to within roughly 10 per cent.

In addition to the jet simulations described above, we also tested our conclusion that the outflow does not become spherical for $\alpha>2$. To do so, we compare with the results of a different simulation, in which an ellipsoid filled with a uniform high pressure is placed into a $\rho\propto R^{-2.5}$ density profile and allowed to evolve. Since there is no need to resolve a jet in this case, the flow can be tracked to much larger radii, until it becomes approximately self-similar. The initial aspect ratio of the outflow is $a/b=70/3\gg 1$, and the ambient pressure is set to a small value ($10^{-8}$ times the initial cocoon pressure) so that it does not affect the evolution.

In Fig. 17, we compare the shape of the cocoon seen in the simulation to the prediction of the analytical model. The figure shows a snapshot of the simulation at $t/t_{0}=840$, when the height is $z_{\rm c}/a=10.4$. As expected, we find that the outflow does not become spherical when $t$ grows large. The shape in the analytic model, using the same values of $z_{\rm c}/a$ and $b/a$, is drawn as a solid black line. The shape observed in the simulations is slightly narrower towards the axis and wider towards the base, but otherwise agrees well with the analytical prediction. In particular, the ratio of the cocoon’s width to its height is extremely close to the analytically-predicted value of 0.62.

However, although the agreement at a given height is very good, the agreement at a given time is noticeably worse. Although the cocoon’s size grows as $\propto t^{2/(5-\alpha)}$ in the simulations (as expected), the simulation reaches the quasi-spherical phase more quickly than anticipated analytically. As a result, the size of the cocoon in the simulation is larger than the analytically-estimated size at a given time. At late times, the ratio of the sizes approaches a constant value of $\approx 2.5$. For example, whereas the simulations reach a height of $z_{\rm c}/a=10.4$ at a time $t/t_{0}=840$, the analytical model predicts a height of only $z_{\rm c}/a=4.2$ at this time, as shown by the dashed line in Fig. 17. One possible reason for the discrepancy is the fact that, as discussed above, there is a clear pressure gradient along the cocoon surface in the simulations, with the pressure increasing by a factor of a few towards the equator. The higher pressure towards the sides enhances the lateral expansion at early times, causing the outflow to transition to the quasi-spherical regime sooner than it would if the pressure were constant across the surface.

From these considerations, we conclude that our results describing the shock’s shape as a function of height (i.e., eq. 67 and Fig. 10) are rather accurate, while our expressions relating the height and width to time (i.e., eqs. 76–77, and Fig. 12) are only reliable to within a factor of a few. An accurate calibration of the analytical model to more extensive numerical simulations is reserved as an interesting topic for future study.

6 Discussion and Conclusions

Jets which fail to penetrate the surrounding medium are expected to leave behind a narrow cocoon filled with hot gas of almost uniform pressure. To describe the evolution of these systems, we consider an axisymmetric generalization of the ST solution in a power-law ambient density profile $\rho\propto R^{-\alpha}$, in which the explosion energy is deposited throughout an elongated region along the polar axis. Applying the Kompaneets approximation (KA), we derive a differential equation governing the shock motion, and solve it to yield analytical expressions for the shape of the shock. For a given curve $R_{i}(\theta_{i})$ describing the initial shock shape, we obtain exact parametric solutions for the shock coordinates $R(\theta_{i},z_{\rm c})$ and $\theta(\theta_{i},z_{\rm c})$ as functions of $z_{\rm c}$, the shock height along the axis. The solutions vary in form depending on $\alpha$. In addition, we provide simplified expressions for the height, width, and volume of the shock versus time in the case of an initially ellipsoidal cocoon with an aspect ratio $a/b\gg 1$, as expected for a jet. Our solution extends the previous work of B11, which describes the dynamics of the system while the jet is active, to times after energy injection has ended.

We find that the evolution of the outflow is characterized by three dynamical regimes, which we call the planar regime, the sideways expansion regime, and the quasi-spherical regime. Two important dynamical time-scales are the time for the outflow’s width to double ($t_{b}$) and the time for its height to double ($t_{a}$). For a jet choked at time $t_{0}$, the planar regime persists until the age is $\sim 2t_{0}$. During this time, dynamical quantities remain roughly constant, and the outflow more or less keeps its original shape. Once $t\sim 2t_{b}\sim 2t_{0}$, sideways expansion starts to become important; the pressure in the cocoon drops and the shock starts to decelerate. The shape also changes, with the cocoon bulging out towards the tip. This is especially pronounced in steep density gradients where there is a large density contrast between the base and the tip. Most of the sideways expansion occurs over the time for the cocoon’s height to double, which is $t\sim t_{a}\sim(a/b)^{2}t_{0}$, where $a/b$ is the aspect ratio of the outflow upon choking. After this, the evolution becomes quasi-spherical, the dependence on the initial shape becomes weak, and the shock’s shape is governed primarily by the density profile. As $a/b$ approaches unity, the two characteristic time-scales ($t_{b}$ and $t_{a}$) become the same and the usual spherical blast wave solution is recovered.

Although the height ($z_{\rm c}$) and width ($r_{\rm c}$) of the cocoon become comparable around $t_{a}$, the time for them to become equal (to within, say, 10 per cent) can be considerably longer. In general we find that the asymptotic behaviour depends strongly on the density profile. For $\alpha\leq 2$, the flow becomes spherical. The convergence to a sphere is slower for steeper density profiles, with the quantity $1-r_{\rm c}/z_{\rm c}$ going to zero as $\propto z_{\rm c}^{(\alpha-2)/2}$ for $\alpha<2$, or as $\propto(\ln z_{\rm c})^{-1}$ for $\alpha=2$. An approximation for the time it takes for the outflow to become spherical in different density profiles is given in eq. 79. We find that in a constant-density medium, the flow spherizes after about $\sim 10t_{a}$, while in a wind-like medium it takes $\sim 1000t_{a}$.

In steep density profiles with $\alpha>2$, the flow does not become spherical (at least in the idealized case where the KA applies). Instead, the ratio $r_{\rm c}/z_{\rm c}$ converges to a constant value of $[\sin(\pi/\alpha)]^{\alpha/(\alpha-2)}$ as ${\propto z_{\rm c}^{(2-\alpha)/4}}$. As before, the convergence becomes slower as $\alpha\rightarrow 2$. For $2<\alpha<3$, the shape of the shock at infinity is narrower towards the equator, like a peanut.

Comparing the analytical solution with numerical hydrodynamics simulations, we find a generally good agreement for $z_{\rm c}/a\la$ a few, even without calibration. We compared our results with realistic simulations of a choked jet for three different density profiles of astrophysical interest, with $\alpha=1$, 2, and 2.5. In each case, the height and the width of the outflow estimated analytically match the numerically determined value to within about 15 per cent. The agreement is not so good near the equator, but the KA is known to have issues in that region. In addition, we explored the late-time behaviour in the $\alpha=2.5$ case through comparison with a simplified simulation in which an elliptical region is filled with an initially high pressure and left to evolve for $840t_{0}$. For $z_{\rm c}/a\gg 1$, we find that the size of the simulated outflow is larger than predicted analytically by a factor of $\approx 2.5$. However, when rescaled to be the same size, the analytical and numerical shapes are in close agreement.

The fact that the time-scale to become quasi-spherical is long compared to the working time of the jet is good news when it comes to observational applications, since it means that it is not necessary to catch the system shortly after jet activity ceases in order to meaningfully constrain the jet. The relationship between the parameters of the choked jet model (initial height $a$, initial width $b$, and the energy $E_{0}$) and the parameters characterizing the jet (luminosity $L_{\rm j}$, duration $t_{\rm j}$, and opening angle $\theta_{\rm op}$) is discussed in Section 4. By measuring the size and shape of a choked jet outflow, along with its energy, it is in principle possible to work backwards to obtain the initial conditions at the moment of choking, and therefore derive the jet properties. However, because the system quickly becomes insensitive to the initial width, constraining $b$ is difficult in practice, unless we happen to catch the system early on. The situation where only $E_{0}$ and $a$ can be estimated from observations is expected to be more typical. In this case, the jet parameters are subject to degeneracy and cannot be determined uniquely. None the less, knowing the location where the jet was choked is a powerful constraint, and one which is not attainable in spherically symmetric models.

We conclude by discussing several possible future applications of our results. The first involves quenched AGN jets in galaxy clusters. These jets inflate bubbles filled with shocked jet debris that initially expand supersonically, driving a shock into the ambient medium, and later drift away due to buoyancy. Our model can be used to describe the shape of the forward shock surrounding the bubbles during the early evolution when the shocked jet material has not yet reached pressure equilibrium with the intracluster medium or been deformed by buoyancy. The size, shape, and speed of the forward shock can be directly measured from X-ray observations, and then compared with the predictions of the model to derive constraints on the jet properties. This idea is developed further and applied to 5 AGNs with suitably young bubbles in a separate paper (Irwin et al., in preparation).

Another possible application deals with observations of SNe that may have harboured a jet. We envision a scenario where the jet is choked before breaking out of its progenitor star, at a radius $a<R_{*}$ where $R_{*}$ is the star’s radius. The cocoon continues to propagate, eventually breaking out of the star when its height is $z_{\rm c}=R_{*}$. If the choking does not occur too deep inside the star, the width of the outflow upon breakout, $r_{\rm c}$, will be small compared to $R_{*}$, and the mass swept up by the cocoon will also be small compared to the star’s mass $M_{*}$. Once outside the star, the cocoon is expected to quickly expand sideways and engulf the star due to the much lower density. However, since a typical GRB jet imparts an energy of $10^{51}$ erg which is comparable to the SN energy, the cocoon will expand faster than the SN due to its lower mass. Therefore, a prediction for the choked jet scenario is the presence of a small amount of mass ($0.01$–$0.1M_{*}$) surrounding the supernova that is moving with a speed several times greater than the bulk SN velocity.

This fast-moving ejecta, which is only visible at early times, may have already been observed in some objects (e.g., Piran et al., 2019; Izzo et al., 2019). Early-time spectroscopy makes it possible to estimate the total mass ($M_{\rm ej}$) and energy ($E_{\rm ej}$) contained in this high-velocity component, and can even constrain the distribution of mass with velocity. $E_{\rm ej}$ corresponds to the cocoon energy $E_{0}$ in our model, while $M_{\rm ej}$ serves as a proxy for the width at breakout, with $M_{\rm ej}\sim M_{*}(r_{\rm c}/R_{*})^{2}$. By applying our results, this information can be translated into an estimate of the choking radius, which provides significant insight into the jet properties. Whether the choked jet model can reproduce the observed mass-velocity distribution in Piran et al. (2019) is an interesting question for future study.

The asphericity of the cocoon also has important implications for shock breakout, where even a small departure from spherical flow can lead to oblique effects that dramatically alter the observed signature (e.g., Matzner et al., 2013). Therefore, we typically expect the breakout of a choked jet to look markedly different than a fully spherical shock breakout. At the least, an aspherical breakout will last longer and be fainter than a spherical breakout of the same energy, because it takes some time for the breakout to wrap around from the pole to the equator. The expected changes in the breakout emission when it is powered by a choked jet, rather than a spherical SN explosion, will be addressed in future work.

As these examples show, a 2D analytical model for the shape of the shock in choked jet outflows has a wide variety of astrophysical applications, and offers significant advantages over the standard simplifying assumption of spherical symmetry. The Kompaneets approximation model developed here is a useful tool which is straightforward and inexpensive to implement, without sacrificing too much accuracy. Bridging the gap between the jet phase and the spherical phase of evolution is an important step towards a complete picture of jet dynamics, which will remain important as observations continue to reveal evidence for suffocated jets in diverse astrophysical systems.

Acknowledgements

We thank O. Gottlieb, who conducted the numerical simulations discussed in Section 5 and graciously shared his data with us, for his generous assistance; his contributions greatly strengthened the final paper. This research is supported by the CHE-ISF I-Core Center for Excellence in Astrophysics. TP is supported by an advanced ERC grant TReX. This work was supported in part by the Zuckerman STEM Leadership Program. CI and EN were partially supported by consolidator ERC grant JetNS and by an ISF grant 1114/17.

Appendix A Glossary of symbols

Many of the quantities appearing in our problem evolve in time. Additionally, at a given time, certain quantities such as the shock velocity vary along the cocoon surface. Since each point on this surface had a unique initial angular position, $\theta_{i}$, at the moment the jet was choked (see Section 3), we use $\theta_{i}$ to parametrize the relative position on the shock front. In Table 2, we list the symbols defined throughout the paper, explicitly including any dependence on $t$ or $\theta_{i}$ in parentheses following each symbol.

Appendix B Parametric equations for the cocoon shape

In this appendix, we derive parametric equations for the height and width of the cocoon as a function of the parameter $\theta_{\rm b0}$. As discussed in Section 3.2, a trajectory originating from $\theta_{\rm b0}$ becomes parallel to the equator at the point ($R_{\rm b},\theta_{\rm b}$), with $\theta_{\rm b}$ and $R_{\rm b}$ given by equations 36 and 37 respectively. How long does it take a particle at $(R_{\rm b0},\theta_{\rm b0})$ to travel to $(R_{\rm b},\theta_{\rm b})$? By substituting $\theta(x(\theta_{\rm b0}),\theta_{\rm b0})=\theta_{\rm b}$ into equation 34 or 52, as appropriate, we obtain

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

to make the formulae more compact. The values of $R_{\rm b0}$ and $\chi_{\rm b0}$ are given by setting $\theta_{i}=\theta_{\rm b0}$ in equations 13 and 14, respectively.

If the initial shape is an ellipsoid with $R_{\rm b0}$ and $\chi_{\rm b0}$ given by equations 13 and 14 at $\theta_{i}=\theta_{\rm b0}$, then equation 101 only has a solution when $x$ is greater than a critical value,

[TABLE]

Physically, this corresponds to the time when the bulge width becomes equal to the equatorial width, i.e. we have $r_{\rm b}=r_{\rm eq}$ at the time $t=t(x_{b})$. The relation between $r_{\rm c}$ and $z_{\rm c}$ is different before and after this transition, with $r_{\rm c}=r_{\rm eq}$ up to $x=x_{b}$, and $r_{\rm c}=r_{\rm b}$ afterwards. Note that if $\alpha\leq 2b^{2}/a^{2}$, $x_{b}$ is negative or undefined. Since $x$ can only be positive and finite, this means that when $\alpha\leq 2b^{2}/a^{2}$, the transition at $x=x_{b}$ never occurs. The bulge width never exceeds the equatorial width in this case, and $r_{\rm c}=r_{\rm eq}$ at all times. However, as long as $b/a$ is small, this is only relevant to very flat density profiles with $\alpha\approx 0$. In the case where $b\ll a$ and $\alpha\gg 2b^{2}/a^{2}$, the equatorial phase lasts only briefly, until the height has increased by $ax_{b}\approx(2/\alpha)(b/a)^{2+k_{\alpha}}a\ll a$, and the bulge-dominated case is the relevant one for most of the evolution.

When the width is determined by equatorial expansion (i.e. when $\alpha\leq 2b^{2}/a^{2}$, or when $\alpha>2b^{2}/a^{2}$ but $x\leq x_{b}$), the cocoon shape parameters are simple functions of the parameter $x$. The height is given by equation 38 or 54, and the equatorial radius is found by taking $\theta_{\rm b0}=\chi_{\rm b0}=\pi/2$ and $R_{\rm b0}=b$. We then have $r_{\rm c}=r_{\rm eq}$, and

[TABLE]

In this regime, it is possible to express $r_{\rm c}$ as an explicit function of $z_{\rm c}$:

[TABLE]

For $k_{\alpha}=1$ (constant density), we reproduce equation 27.

When $\alpha>2b^{2}/a^{2}$ and $x>x_{b}$, however, the width is governed by the bulge, and the relation between $r_{\rm c}$ and $z_{\rm c}$ becomes more nuanced. It is no longer possible to express $r_{\rm c}$ as an explicit function of $z_{\rm c}$, but the width and height can be parametrized in terms of $\theta_{\rm b0}$. The height of the cocoon when $x=x(\theta_{\rm b0})$ is, combining equations 38 or 54 with 101,

[TABLE]

with $\theta_{\rm b}$ given by equation 36. The width, on the other hand, is $r_{\rm c}=r_{\rm b}=R_{\rm b}\sin\theta_{\rm b}$; from equations 36 and 37, we have

[TABLE]

Similarly, the height of the bulge above the equatorial plane is $z_{\rm b}=R_{\rm b}\cos\theta_{\rm b}$, or

[TABLE]

The equatorial width in this case is discussed separately, in Appendix D. In the first part of equations 107–109, we simplified the denominator by using the identity $\omega_{\rm b0}-k_{\alpha}\Delta_{\rm b}=\pi/2-\theta_{\rm b}$, which can be derived from equations 102–103.

Since we expect $r_{\rm c}\propto z_{\rm c}$ to hold approximately at late times, it is sometimes more convenient to work with the dimensionless ratio $r_{\rm c}/z_{\rm c}$. Dividing equation 108 by equation 107, we obtain

[TABLE]

Appendix C Evolution of the cocoon width versus height

Approximate formulae for $r_{\rm c}$ versus $z_{\rm c}$ can be obtained in certain limiting situations where equations 107 and 108 simplify. To understand the interesting limiting cases, let us investigate how the three angular scales appearing in equation 110 ($\theta_{\rm b}$, $\omega_{\rm b0}$, and $\Delta_{\rm b}$) vary with $\theta_{\rm b0}$. (We restrict the discussion here to an initially ellipsoidal cocoon, but the conclusions can be generalized to other initial shapes.) To guide the discussion, we plot $\chi_{\rm b0}$, $\omega_{\rm b0}\equiv\chi_{\rm b0}-\theta_{\rm b0}$, $\theta_{\rm b}$, and $\Delta_{\rm b}\equiv\theta_{\rm b}-\theta_{\rm b0}$ versus $\theta_{\rm b0}$ in Fig. 18 for $\alpha=1$ (left panel) or 3 (right panel) and $a/b=10$. Focusing now on the $\theta_{\rm b}$ and $\omega_{\rm b0}$ curves, the key features are as follows:

• •

Differentiation of equation 36 implies that $\theta_{\rm b}$ has a minimum if ${\rm d}\chi_{\rm b0}/{\rm d}\theta_{\rm b0}=\alpha/2$ is satisfied. Applying equation 14, we see that this condition is only fulfilled when $\frac{2b^{2}}{a^{2}}<\alpha<\frac{2a^{2}}{b^{2}}$. If we assume that $b\ll a$, and that $\frac{2b^{2}}{a^{2}}\ll\alpha\ll\frac{2a^{2}}{b^{2}}$ holds as well,666The case of $\alpha\leq\frac{2b^{2}}{a^{2}}$ describes the scenario where the cocoon always remains widest in the equatorial plane, as discussed above. The case of $\alpha\geq\frac{2a^{2}}{b^{2}}$ we do not discuss further, because if $a\gg b$, it requires an extreme value of $\alpha$. then $\theta_{\rm b}$ has a minimum value of $\approx\frac{2b}{a}\sqrt{\frac{2}{\alpha}}$ at $\theta_{\rm b0}\approx\frac{b}{a}\sqrt{\frac{2}{\alpha}}$. Away from the minimum, $\theta_{\rm b}$ rises to $\pi/2$ at $\theta_{\rm b0}=\pi/2$, and to $\pi/\alpha$ at $\theta_{\rm b0}=0$.777If $\theta_{\rm b}>\pi/2$ for some $\theta_{\rm b0}$, the physical interpretation is that trajectories stemming from $\theta_{\rm b0}$ never become parallel to the equator. If $\alpha<2$, $\theta_{\rm b}$ becomes equal to $\pi/2$ at $\theta_{\rm b0}=\theta_{\rm b0,min}$ (see below).

• •

$\omega_{\rm b0}$ is determined from equation 14, and its behaviour is straightforward since there is no dependence on $\alpha$: $\omega_{\rm b0}$ is 0 when $\theta_{\rm b0}=0$, rises to a maximum value of $\omega_{\rm b0}\approx\pi/2-b/a$ at $\theta_{\rm b0}\approx b/a$, and then goes to zero again as $\theta_{\rm b0}\rightarrow\pi/2$.

• •

The $\theta_{\rm b}$ and $\omega_{\rm b0}$ curves have two crossings, one at $\theta_{\rm b0}\sim b^{2}/a^{2}$, and another at $\theta_{\rm b0}\sim\pi/4$. Over the range $b^{2}/a^{2}\la\theta_{\rm b0}\la\pi/4$, $\theta_{\rm b}<\omega_{\rm b0}$; otherwise $\theta_{\rm b}>\omega_{\rm b0}$.

As discussed in Section 3.2.1, $\theta_{\rm b0,min}$ corresponds to a special trajectory separating the ejecta that eventually interact at the equator from those that do not. To find its value, we look for a trajectory that has $\theta\rightarrow\pi/2$ as $x\rightarrow\infty$. Then, from equation 34, we get an implicit equation,

[TABLE]

where $R_{\rm b0,min}$ is related to $\theta_{\rm b0,min}$ by equation 14. Alternatively, formula 111 may be obtained by taking $\theta_{\rm b}\rightarrow\pi/2$ in  36. Note that when $\alpha=2b^{2}/a^{2}$, the solution to equation 111 is $\theta_{\rm b0,min}=R_{\rm b0,min}=\pi/2$. This is another way to see that there is no interaction in the $\alpha\leq 2b^{2}/a^{2}$ case. As $\alpha$ increases from this value, $\theta_{\rm b0,min}$ decreases, and a larger portion of the cocoon becomes subject to collisions. When $\alpha\rightarrow 2$, $\theta_{\rm b0,min}\rightarrow 0$, and all parts of the surface will interact given sufficient time. As long as $\alpha\gg 2b^{2}/a^{2}$, $\theta_{\rm b0,min}$ is small and can be approximated by

[TABLE]

The outflow’s shape evolves in a different way depending on the relative size of the angular parameters. There are four distinct evolutionary phases (labeled by Roman numerals in the figure) as $\theta_{\rm b0}$ is decreased from $\pi/2$ to $\theta_{\rm b0,min}$ (when $\alpha<2$) or to [math] (when $\alpha\geq 2$). These four phases are closely related to the dynamical phases discussed in Section 2. The relationship between the angular scales in each phase is as follows:

1. I.

$\theta_{\rm b0}\approx\pi/2$ and $\chi_{\rm b0}>\theta_{\rm b0}\approx\theta_{\rm b}\gg\omega_{\rm b0}\gg\Delta_{\rm b}$. This corresponds to the early planar phase, when $x\sim x_{b}$ (see also Appendix B), and the bulge is still located close to the equator. 2. II.

$b/a\ll\theta_{\rm b0}\ll\pi/4$ and $\chi_{\rm b0}>\omega_{\rm b0}\gg\theta_{\rm b0}\approx\theta_{\rm b}\gg\Delta_{\rm b}$. This applies during the later part of the planar phase, once $x\gg x_{b}$ and the bulge is closer to the axis than to the equator. 3. III.

$b^{2}/a^{2}\ll\theta_{\rm b0}\ll b/a$ and $\chi_{\rm b0}\simeq\omega_{\rm b0}\gg\theta_{\rm b}\simeq\Delta_{\rm b}\gg\theta_{\rm b0}$. This holds during the sideways expansion phase. 4. IV.

$\theta_{\rm b0}\rightarrow\theta_{\rm b0,min}$ or $\theta_{\rm b0}\rightarrow 0$ (respectively for $\alpha<2$ and $\alpha\geq 2$), and $\theta_{\rm b}\approx\Delta_{\rm b}>\chi_{\rm b0}\approx\omega_{\rm b0}\gg\theta_{\rm b0}$. This is appropriate in the quasi-spherical phase.

We consider each phase in turn below.

Phase I ($\theta_{\rm b0}\approx\pi/2$:) Let $\delta\equiv\pi/2-\theta_{\rm b0}$ be a small quantity. Then, expanding equations 14 and 36 in terms of $\delta$, we find $\chi_{\rm b0}\approx\pi/2-(b^{2}/a^{2})\delta$, $\omega_{\rm b0}\approx(1-b^{2}/a^{2})\delta$, $\theta_{\rm b}\approx\pi/2-(1-2b^{2}/\alpha a^{2})\delta$, and $\Delta_{\rm b}\approx(2b^{2}/\alpha a^{2})\delta$. In this case, the term in brackets in equation 110 is unity, up to a small correction of order $b^{2}/a^{2}$. Therefore we can write $r_{\rm c}/z_{\rm c}\approx(R_{\rm b0}/a)\sin\theta_{\rm b}\approx b/a$.

Phase II ($b/a\la\theta_{\rm b0}\la\pi/4$): In this limit $\delta=\pi/2-\theta_{\rm b0}\geq\pi/4$ is no longer a small quantity. However, $(b^{2}/a^{2})\delta$ is still small. So, to a reasonable approximation $\theta_{\rm b}\approx\theta_{\rm b0}$, $\chi_{\rm b0}\approx\pi/2$, and $\omega_{\rm b0}\approx\pi/2-\theta_{\rm b0}$. Thus the bracketed part of equation 106 is still order-unity, and $r_{\rm c}/z_{\rm c}\approx(R_{\rm b0}/a)\sin\theta_{\rm b}$ remains valid. Although $R_{\rm b}$ and $\theta_{\rm b}$ both vary with $\theta_{\rm b0}$ in this case, one decreases while the other increases. The effects roughly cancel so that $r_{\rm b}$ turns out to be nearly constant with respect to $\theta_{\rm b0}$. When phase II ends, $\theta_{\rm b}$ is near its minimum value of $\approx(2b/a)\sqrt{2/\alpha}$. At that time, $R_{\rm b}\approx R_{\rm b0}\approx a/\sqrt{1+2/\alpha}$ via equations 13 and 37, and $r_{\rm c}/z_{\rm c}\approx R_{\rm b}\theta_{\rm b}/a\approx(2b/a)/\sqrt{1+\alpha/2}$. We therefore find that $r_{\rm c}/z_{\rm c}\simeq b/a$ throughout this phase, to within a factor of two. This suggests that the end of phase II roughly coincides with the end of the planar phase.

Phase III ($b^{2}/a^{2}\la\theta_{\rm b0}\la b/a$): In phases III and IV, $\theta_{\rm b0}$ is the smallest angular scale in the problem. In the limit $\theta_{\rm b0}\rightarrow 0$, the shock’s shape is similar to that of an on-axis point explosion at $z=a$. Since $\theta_{\rm b0}$ can be neglected, we have $\theta_{\rm b}\approx\Delta_{\rm b}$ and $\chi_{\rm b0}\approx\omega_{\rm b0}\approx\pi/2-(\alpha/2)\theta_{\rm b}$. Equations 106 and 107 then simplify to functions of $\theta_{\rm b}$ alone:

[TABLE]

and

[TABLE]

these expressions are valid for both phase III and phase IV, since so far we have only assumed $\chi_{\rm b0}\gg\theta_{\rm b0}$ and $\theta_{\rm b}\gg\theta_{\rm b0}$. In what follows it will also be useful to know how the parameter $\zeta$ varies with $\theta_{\rm b}$. Combining equation 113 with the definition of $\zeta$ from equation 43, we obtain

[TABLE]

For phase III specifically, we want $\theta_{\rm b}\ll\chi_{\rm b0}$, and since $\chi_{\rm b0}\leq\pi/2$ this implies $\theta_{\rm b}$ must be small. With the approximation $R_{\rm b0}\approx a$, expanding 113 and 114 to leading order for $\theta_{\rm b}\ll 1$ yields simply $r_{\rm c}/z_{\rm c}\simeq\theta_{\rm b}$ and $z_{\rm c}\simeq a(1+\theta_{\rm b})$. Therefore,

[TABLE]

We caution that this approximation is coarse and is only relevant for sufficiently large $a/b$. The $r_{\rm c}\propto z_{\rm c}-a$ behaviour can be understood through analogy with a point explosion at $z=a$. Within a small neighbourhood of size $\ll a$ near the tip of the cocoon, the density is nearly constant. Since parts of the cocoon moving in different directions encounter the same density, the shocked region near $z_{\rm c}=a$ is a roughly sphere of radius $z_{\rm c}-a$. As phase III ends when $\theta_{\rm b}\sim 1$ and $z_{\rm c}\sim a$, it can be identified with the sideways expansion phase discussed in Section 2.

Phase IV ($\theta_{\rm b0}\la b^{2}/a^{2}$): Finally, we address the case where $\theta_{\rm b}>\chi_{\rm b0}\gg\theta_{\rm b0}$. This limit describes the behaviour as the outflow’s size tends to infinity. In contrast to the previous phases, the evolution in phase IV depends strongly on $\alpha$. Below, we separately discuss the cases $\alpha<2$, $\alpha=2$, and $\alpha>2$. The basic solution method is the same in each case. First, we create some small quantity $\delta\equiv\theta_{\rm b\infty}-\theta_{\rm b}$, with $\theta_{\rm b\infty}$ defined as in equation 40. Next, we replace $\theta_{\rm b}$ by $\delta$ in equation 114, and expand to leading order in $\delta$. We then do the same for equation 115 to obtain $\zeta$ in terms of $\delta$. Finally, we combine the two expressions to determine $r_{\rm c}/z_{\rm c}$ to first order in $\zeta$.

For $\alpha<2$, we have $\theta_{\rm b}\rightarrow\pi/2$ as $\theta_{\rm b0}\rightarrow\theta_{\rm b0,min}$, so we let $\delta=\pi/2-\theta_{\rm b}$. Then, ignoring terms of order $\delta^{2}$ or higher, we find $\cos\theta_{\rm b}\approx\delta$; $\sin\theta_{\rm b}\approx 1$; $\cos(\alpha\theta_{\rm b}/2)\approx\sin(k_{\alpha}\pi/2)[1+(\alpha/2)\delta\cot(k_{\alpha}\pi/2)]$; and $\sin(k_{\alpha}\theta_{\rm b})\approx\sin(k_{\alpha}\pi/2)[1-k_{\alpha}\delta\cot(k_{\alpha}\pi/2)]$. Applying equation 115, we obtain

[TABLE]

where $R_{\rm b0,min}\equiv R_{i}(\theta_{i}=\theta_{\rm b0,min})$. Now, expanding equation 114 and using 117 to replace $\delta$ by $\zeta$, we find

[TABLE]

If we further assume that $\theta_{\rm b0,min}$ is small (i.e., that $\alpha\gg 2b^{2}/a^{2}$), then $R_{\rm b0,min}\approx a$ and equation 118 reduces to the form in equations 42–44.

For $\alpha=2$, we can repeat the same procedure using the second expressions in 114 and 115, or simply take the limit $k_{\alpha}\rightarrow 0$ in formula 118. Either way, we recover equation 57.

Lastly, we treat the case of $\alpha>2$, for which $\theta_{\rm b0}\rightarrow 0$ and $\theta_{\rm b}\rightarrow\theta_{\rm b\infty}=\pi/\alpha$. This time, we suppose $\delta=\pi/\alpha-\theta_{\rm b}$ is small. In this limit, it can be shown from equations 21 and 36 that $R_{\rm b0}/a=1-\mathcal{O}(\delta^{2})$, so we ignore corrections arising from the $R_{\rm b0}/a$ factors. Using $k_{\alpha}=1-\alpha/2$, the denominator appearing in equations 114 and 115 can be rewritten as:

[TABLE]

The term in brackets can be neglected, since it vanishes as $\delta^{2}$ when $\delta\rightarrow 0$. Therefore, equations 114 and 115 simplify to

[TABLE]

and

[TABLE]

When $\theta_{\rm b}\rightarrow\pi/\alpha$, we see that $r_{\rm c}/z_{\rm c}$ approaches the value $f_{\alpha}$ given by equation 62.

Expanding once more in terms of $\delta$ yields $\sin\theta_{\rm b}\approx\sin(\pi/\alpha)[1-\delta\cot(\pi/\alpha)]$; $\sin(k_{\alpha}\theta_{\rm b})\approx-\cos(\pi/\alpha){[1-k_{\alpha}\delta\cot(k_{\alpha}\pi/\alpha)]}$; and $\cos(\alpha\theta_{\rm b}/2)\approx(\alpha/2)\delta$. So, to leading order,

[TABLE]

Expanding formula 120 up to first order in $\delta$, and then substituting $\delta$ for $\zeta$ according to equation 122, we arrive at equations 61–63.

Appendix D Collisions and the evolution of the equatorial width versus height

Just like the other shape parameters, the equatorial radius $r_{\rm eq}$ can also be parametrized in terms of $\theta_{\rm b0}$. Before continuing, we stress again that the presence of shocks near the equator (formed either by the jet’s passage or the collision of cocoon material) may significantly alter the evolution of $r_{\rm eq}$. As compared to numerical simulations (see Section 5), the equations presented here are only accurate to within a factor of a few. However, the true equatorial radius cannot be less than the analytical estimate which ignores interactions.

Consider the ejecta that are just arriving to the equator at the time $x=x(t)$. These ejecta must satisfy $\theta(x(\theta_{\rm eq0}),\theta_{\rm eq0})=\pi/2$, where $\theta_{\rm eq0}$ is their initial angle. Similar to the way that the width of the bulge can be parametrized in terms of $\theta_{\rm b0}$ (see Appendix B), the width at the equator can be parametrized in terms of $\theta_{\rm eq0}$. Equation 34 then leads to the condition

[TABLE]

where (in analogy with Appendix B) we defined

[TABLE]

and

[TABLE]

and $R_{\rm eq0}$ and $\chi_{\rm eq0}$ are now given by taking $\theta_{i}=\theta_{\rm eq0}$ in equations 13 and 14.

Depending on the density profile and how much time has elapsed since choking, there can be 0, 1, 2, or 3 values of $\theta_{\rm eq0}$ which satisfy $x(t)=x(\theta_{\rm eq0})$ on the interval $0<\theta_{\rm eq0}<\pi/2$. The number of solutions for $\theta_{\rm eq0}$ corresponds to the number of places where the cocoon surface intersects the equator (not counting the point $\theta_{\rm eq0}=\pi/2$). If there are no solutions for $\theta_{\rm eq0}$, the equatorial radius $r_{\rm eq}$ is determined by the ejecta which started at $\theta_{i}=\pi/2$, and $r_{\rm eq}$ is described by equation 106. Otherwise, the value of $r_{\rm eq}$ is determined by the intersection point farthest from the axis, which is set by the smallest solution for $\theta_{\rm eq0}$. We then have $r_{\rm eq}=R(x(\theta_{\rm eq0}),\theta_{\rm eq0})$. Inserting $x=x(\theta_{\rm eq0})$ into equation 33 and simplifying, we obtain

[TABLE]

To compare $r_{\rm eq}$ to the values of $z_{\rm c}$, $r_{\rm c}$, and $z_{\rm b}$ given in Appendix B, it is necessary to determine the relationship between $\theta_{\rm eq0}$ and $\theta_{\rm b0}$. This is done by equating the two expressions for $x$,  101 and 123, and solving numerically to obtain $\theta_{\rm eq0}$ in terms of $\theta_{\rm b0}$, or vice versa.

Equation 123 is useful for determining if and when different parts of the surface reach the equator and undergo a collision. In general, there are two ways that collisions can occur. The first way involves material just adjacent to the equator that gets pushed down into the equatorial plane as the cocoon expands. This material tends to impact the equator at a relatively shallow angle; therefore, we call collisions that happen in this way ‘glancing collisions.’ Another type of collision can also occur in steep density profiles, because the trajectories are highly curved, and due to the large density contrast, material at large radii travels considerably faster than material at small radii. If $\alpha$ is sufficiently large, material originating near the tip of the cocoon wraps around along a curved trajectory and impacts the equator before material that started near the equator. The first material to reach the equator in this way has a velocity perpendicular to the equator at the moment of impact. We therefore call this a ‘head-on collision.’ As we will see, the relationship between $\alpha$ and the ratio $b/a$ governs whether or not collisions occur, and also which type of collisions is relevant. In addition, the behaviour is different for $\alpha\leq 2$ and $2<\alpha<3$. We discuss the various possibilities in detail below.

Let us first consider the case $\alpha\leq 2$. In this scenario, we find that there is no solution to equation 123 when $\alpha\leq 2b^{2}/a^{2}$, indicating that no collisions occur in this case. When $\alpha>2b^{2}/a^{2}$, we find that glancing collisions start to occur at $x=x_{eq}$, where

[TABLE]

for an initially ellipsoidal cocoon. Just as $r_{\rm c}$ has a transition at $x=x_{b}$ (as discussed in Appendix B), $r_{\rm eq}$ breaks to a different behaviour at $x=x_{eq}$. While $x<x_{eq}$, there is no solution for $\theta_{\rm eq0}$ and $r_{\rm eq}$ is given by equation 106. Once $x>x_{eq}$, ejecta from $\theta_{i}<\pi/2$ start to arrive at the equator, and there is one solution for $\theta_{\rm eq0}$. Thereafter, $r_{\rm eq}$ is determined by equation 126. Head-on collisions never occur for $\alpha\leq 2$ (see below).

When $\alpha>2$, the situation is more complicated, because both glancing and head-on collisions are possible. To determine if and when a head-on collision occurs, we look for trajectories that intersect the equatorial plane at right angles. Such trajectories satisfy $dR/d\theta=0$ at the point $(R,\theta)=(r_{\rm eq},\pi/2)$. From equation 35, we find that ejecta originating from $\theta_{i}$ undergo a head-on collision if the condition

[TABLE]

is satisfied. Solving for $\alpha$ yields

[TABLE]

For an initially ellipsoidal shape, the right hand side of equation 129 has a minimum value. Therefore, equation 129 only has a solution if $\alpha\geq\alpha_{col,ho}$, where

[TABLE]

to leading order in $b/a$. If $\alpha<\alpha_{col,ho}$, equation 128 does not have a solution, and no head-on collision occurs. (Because $\alpha_{col,ho}$ is always greater than 2, this also implies that head-on collisions are impossible for $\alpha<2$.)

In the $\alpha\geq\alpha_{col,ho}$ case, equation 128 admits two solutions. The relevant one for determining $r_{\rm eq}$ is approximated by

[TABLE]

for $b/a\ll 1$. Eq. 131 pertains to the path taken by the first ejecta to experience a head-on collision (e.g., the dot-dashed line in Fig. 9). If we set $\theta_{i}=\theta_{\rm eq0}$, equation 128 implies $\omega_{\rm eq0}\approx\pi(1-\alpha/4)$, $\Delta_{\rm eq}\approx\pi/2$, and $R_{\rm eq0}\approx a$. Inserting these values into equations 123 and 126, we find that the first head-on collision occurs at a time

[TABLE]

and a radius

[TABLE]

The evolution of $r_{\rm eq}$ in the $2<\alpha<3$ case depends on how $\alpha$ compares with $\alpha_{col,ho}$. If $2<\alpha<\alpha_{col,ho}$, material just adjacent to the equator starts to undergo collisions at $x=x_{eq}$, and from then on there are continuous glancing collisions, as in the case $2b^{2}/a^{2}<\alpha\leq 2$ discussed above. In this case, the equatorial radius $r_{\rm eq}$ grows continuously with $x$. However, if $\alpha_{col,ho}\leq\alpha<3$, there is a head-on collision that occurs at $x=x_{col}$. The head-on collision causes the equatorial radius to suddenly increase to $\approx r_{col}$, so that $r_{\rm eq}$ has a discontinuity at $x_{col}$. The value of $r_{\rm eq}$ increases smoothly up for $x<x_{col}$, then jumps to a higher value at $x=x_{col}$, and afterwards continues to increase gradually.

The results of this section are summarized by Fig. 19, where we plot the critical values of $\alpha$ as functions of $b/a$. The phase space of $\alpha$ versus $b/a$ can be divided into the four regions labeled A through D in the figure, each with a distinct behaviour:

• •

A ($\alpha\leq 2b^{2}/a^{2}$): There are no collisions.

• •

B ($2b^{2}/a^{2}<\alpha\leq 2$): The part of the surface with $\theta_{i}>\theta_{\rm b0,min}$ eventually experiences glancing collisions. The rest of the surface never reaches the equator.

• •

C ($2<\alpha<\alpha_{col,ho}$): All the ejecta ultimately undergo a glancing collision.

• •

D ($\alpha_{col,ho}\leq\alpha<3$): All the ejecta eventually experience a collision. For material starting out near the equator, the collision is glancing; for material originating near the axis, the collision is head-on.