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

Irwin, Christopher M; Nakar, Ehud; Piran, Tsvi

doi:10.1093/mnras/stz2268

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 ρ ∝ R^−α, in the case where the shock is decelerating (α < 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 shallow (⁠|$\alpha \leq 2$|⁠) and steep (2 < α < 3) density profiles. Comparing the analytical model to numerical simulations, we find agreement to within ∼15 per cent out to 45 deg 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.

hydrodynamics, shock waves, gamma-ray burst: general, galaxies: jets

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, Meier & Readhead 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. Other systems such as microquasars and isolated black holes are also expected to produce relativistic jets which interact with the surrounding interstellar environment.

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, Woosley & Heger 2001; Zhang, Woosley & Heger 2004; Mizuta et al. 2006; Morsony, Lazzati & Begelman 2007; Wang, Abel & Zhang 2008; Lazzati, Morsony & Begelman 2009; Mizuta & Aloy 2009; Nagakura et al. 2011; López-Cámara et al. 2013; Mizuta & Ioka 2013; Ito et al. 2015; Harrison, Gottlieb & Nakar 2018). Whether describing an AGN or a GRB, the dynamics are captured by the same physical quantities: the kinetic power L_j and Lorentz factor Γ_j ≫ 1 of the jet, the opening angle θ_op into which the jet is injected, and the external density ρ. 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, Nakar & Piran 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). Secondly, in low-luminosity GRBs (llGRBs), 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} \, \mathrm{M}_{\odot }$|⁠) and radius (≳100 R_|$\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 llGRBs (Kulkarni et al. 1998; Campana et al. 2006; Nakar & Sari 2012; Nakar 2015). Interestingly, although they are more difficult to observe, llGRBs may be more common per cosmic volume than standard GRBs (Soderberg et al. 2006), again hinting that choked jets may be common. Thirdly, early spectroscopy of Type Ib and Ic supernovae (SNe) has unveiled a distinct high-velocity component in several cases (e.g. Izzo et al. 2019; Piran et al. 2019). The inferred energy (∼10⁵¹ erg s⁻¹) and velocity (∼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. 2019).

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 processes of jet choking and spherization through numerical simulations (e.g. Ayal & Piran 2001; 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 ∝ t^{2/(5 − α)} if the density obeys ρ ∝ R^−α (Taylor 1950; Sedov 1959). The ST solution applies when the outflow is decelerating, which is true for α < 3. In media of finite mass (α > 3), the shock accelerates with time and eventually becomes causally disconnected from the swept-up mass (Koo & McKee 1990). The case 3 < α < 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 α > 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 α < 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. Secondly, 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 llGRBs, 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. Thirdly, 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 α < 2, α = 2, and 2 < α < 3, which are discussed, respectively, in Sections 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, ρ ∝ R^−α, with a head velocity cβ_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β_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 β_h > β_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₀, we consider the jet to be ‘choked’. The time t₀ 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₀ 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 β_h > β_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₀, with initial size R₀, detonated at time t₀, there is one length-scale (R₀), and one characteristic time-scale, |${t_{\rm ST} \sim (\rho (R_0) 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₀ ≲ t_ST), and the self-similar phase (t − t₀ ≳ t_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₀. 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). Assuming that a > b, t_a is always larger than t_b. The early phase of evolution (t − t₀ ≲ 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₀ ≳ t_a. However, in the axisymmetric solution there is a transitional regime (t_b ≲ t − t₀ ≲ 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₀ ≪ 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 ≪ t − t₀ ≪ 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 ∼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₀ ≫ 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 → ∞, the system approaches an asymptotic self-similar solution, as shown in the lower right panel.

Figure 1.

Comparison of the shape of the shock in each phase of the evolution, for a density profile ρ ∝ R^−α with α = 0, 1, 2, and 3. The axes are scaled to the height of the cocoon, z_c, to allow for an easy comparison. The black dashed line is a sphere of radius z_c. Upper left: The planar phase (t − t₀ < t_b). Upper right: The sideways expansion phase (t_b < t − t₀ < t_a). Lower left: The quasi-spherical phase (t_a < t − t₀). Lower right: The asymptotic shape as t → ∞; note that the lines for α = 0, 1, and 2 are overlapping because in each of these cases the shock becomes spherical.

Open in new tab Download slide

The asymptotic behaviour can be divided into two different regimes depending on α. 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 α > 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 ∼ t₀. Alternatively, we can write t_b ∼ b/cβ_c, with the sideways shock expansion speed β_c determined by balancing the upstream and downstream momentum flux at the lateral shock. Applying the shock jump conditions at R ≈ a gives β_c ≈ [(γ + 1)P₀/2ρ(a)c²]^1/2, where P₀ = (γ − 1)E₀/V₀ is the initial pressure in the cocoon, |$V_{0} \approx \frac{4}{3} \pi a b^2$| is the cocoon’s initial volume, and γ is the adiabatic index. Adopting γ = 4/3, we arrive at an expression for t_b in terms of the cocoon parameters a, b, and E₀:

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

(1)

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. Near the tip, forward-moving material 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 ∼b + a. The volume at that time is roughly |$V_{\rm c}\sim \frac{4}{3} \pi (2a)(a^2) \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 ∼a/v_z, where v_z ∼ [(γ + 1)P_c/2ρ(a)]^1/2 is the on-axis expansion speed, leading to

$$\begin{eqnarray*} t_a \sim \left(\dfrac{48 \pi \rho (a) a^5}{7E_0}\right)^{1/2}. \end{eqnarray*} $$

(2)

The fact that t_a ∼ (a/b)²t₀ implies that an outflow takes much longer to become quasi-spherical if it is initially narrow (i.e. if b ≪ a). This is due to a combination of two effects which cause narrow cocoons to be more strongly decelerated. First, narrow outflows 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, β_h, to a new value, ∼β_c, which is comparable to the sideways expansion speed. Since a ∼ cβ_ht₀ and b ∼ cβ_ct₀, the velocity is reduced by a factor of β_c/β_h ∼ b/a. Secondly, initially narrow cocoons are more strongly affected by sideways expansion. As discussed above, the width of the outflow grows from b to ∼a in the time it takes the height to double. This increases the volume by a factor of ∼(a/b)², and decreases the pressure by a factor of ∼(b/a)². Since |$v_z \propto P_{\rm c}^{1/2}$|⁠, the sideways expansion reduces the velocity by an additional factor of (P_c/P₀)^1/2 ∼ b/a, so by the time the shock reaches a height of 2a, the on-axis velocity is only v_z ∼ (b/a)cβ_c ∼ (b/a)²cβ_h.

For b ≪ a, the cocoon’s height z_c and width r_c in the three dynamical regimes evolve approximately as

$$\begin{eqnarray*} z_{\rm c}\sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}a, & t-t_0 \lesssim t_b \\ a\left[1+ \left(\dfrac{t}{t_a}\right)^{1/2} \right] & t_b \lesssim t-t_0 \lesssim t_a \\ 2a\left(\dfrac{t}{t_a}\right)^{2/(5-\alpha)}, & t-t_0 \gtrsim t_a \end{array}\right. \end{eqnarray*} $$

(3)

and

$$\begin{eqnarray*} r_{\rm c}\sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b, & t-t_0 \lesssim t_b \\ a\left(\dfrac{t}{t_a}\right)^{1/2}, & t_b \lesssim t-t_0 \lesssim t_a \\ a\left[2\left(\dfrac{t}{t_a}\right)^{2/(5-\alpha)}-1\right], & t-t_0 \gtrsim t_a \end{array}\right.. \end{eqnarray*} $$

(4)

The cocoon’s opening angle, θ_c = tan ⁻¹(r_c/z_c), is therefore of order

$$\begin{eqnarray*} \theta _{\rm c} \sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b/a, & t-t_0 \lesssim t_b \\ \left(\dfrac{t}{t_a}\right)^{1/2}, & t_b \lesssim t-t_0 \lesssim t_a \\ 1, & t-t_0 \gtrsim t_a \end{array}\right.. \end{eqnarray*} $$

(5)

(For the derivation of equations 3 and 4, including the accurate determination of order-unity prefactors, see Section 3.) If |$\alpha \leq 2$|⁠, the ratio r_c/z_c approaches 1 as t → ∞. Otherwise, we find that |$r_{\rm c}/z_{\rm c}\rightarrow [\sin (\pi /\alpha)]^{\alpha /(\alpha -2)}$|⁠.

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

$$\begin{eqnarray*} z_{\rm c}-r_{\rm c}\propto z_{\rm c}^{(\alpha -2)/2} \propto t^{(\alpha -2)/(5-\alpha)} \end{eqnarray*} $$

(6)

for α < 2, or as (ln z_c)⁻¹ ∝ (ln t)⁻¹ for α = 2. Because this quantity decreases so slowly, the difference z_c − r_c can be used to constrain the radius at which the jet was choked, even at times considerably larger than t_a. If both the height and width of the outflow are measured through observations, the choking radius can be estimated via

$$\begin{eqnarray*} a \sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}z_{\rm c}- r_{\rm c}, & z_{\rm c}- r_{\rm c}\gtrsim z_{\rm c}/2 \\ [z_{\rm c}^{-\alpha } (z_{\rm c}-r_{\rm c})^2]^{1/{(2-\alpha)}}, & z_{\rm c}- r_{\rm c}\lesssim z_{\rm c}/2. \end{array}\right. \end{eqnarray*} $$

(7)

Knowing the choking radius places strong constraints on the properties of the original jet, as we discuss in Section 4. In particular, if the explosion energy and the density at z = z_c are also known, then the initial jet’s luminosity L_j, opening angle θ_op and duration t_j satisfy

$$\begin{eqnarray*} L_{\rm j}t_{\rm j}= E_0 \end{eqnarray*} $$

(8)

and

$$\begin{eqnarray*} t_{\rm j}\theta _{\rm op}^{-2} \sim \left(\dfrac{\rho (z_{\rm c}) z_{\rm c}^5}{E_0}\right)^{1/2} \left(\dfrac{z_{\rm c}-r_{\rm c}}{z_{\rm c}}\right)^{q_\alpha }, \end{eqnarray*} $$

(9)

where we have defined

$$\begin{eqnarray*} q_\alpha \equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(5-\alpha)/2, & z_{\rm c}- r_{\rm c}\gtrsim z_{\rm c}/2 \\ (5-\alpha)/(2-\alpha), & z_{\rm c}- r_{\rm c}\lesssim z_{\rm c}/2 \end{array}\right.. \end{eqnarray*} $$

(10)

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,

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

(11)

The forward shock bounding the cocoon can be described by a curve θ = θ(R, t) at time t, where θ 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 values measured at the choking time, t = t₀.
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 ‘α’.

Let (R_i, θ_i) be the set of coordinates describing the cocoon surface at the initial time t = t₀ when the jet was choked. We can then parametrize the initial shape of the forward shock by a curve R_i(θ_i). Additionally, we define χ_i(θ_i) as the angle between the surface normal at (R_i, θ_i) and the axis, which is given by

$$\begin{eqnarray*} \chi _i(\theta _i) = \tan ^{-1}\left(\dfrac{\tan \theta _i - \frac{{\rm d}\ln R_i}{{\rm d}\theta _i}}{1+\frac{{\rm d}\ln R_i}{{\rm d}\theta _i} \tan \theta _i} \right). \end{eqnarray*} $$

(12)

A schematic of the initial conditions is shown in Fig. 2. Particles starting with different values of θ_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 θ_i.

Figure 2.

Problem setup. A high-pressure cocoon (purple) expands into an external medium (blue). In a power-law density profile, particles at (R_i, θ_i) travel along the curved pink path to reach the point (R, θ). The initial direction of motion is perpendicular to the shock surface; the angle between the axis and the surface normal is χ_i.

Open in new tab Download slide

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

$$\begin{eqnarray*} R_i(\theta _i) = \dfrac{ab}{\sqrt{a^2 \sin ^2 \theta _i+b^2 \cos ^2 \theta _i }} \end{eqnarray*} $$

(13)

and

$$\begin{eqnarray*} \chi _i(\theta _i) = \tan ^{-1}\left(\dfrac{a^2}{b^2} \tan \theta _i\right). \end{eqnarray*} $$

(14)

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¹ decreases smoothly with θ_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 first condition is necessary because it greatly simplifies the analytical calculation of the cocoon’s volume; the second condition ensures that there are no intersections between trajectories, except at the equator; and the third condition makes certain that trajectories on either side of the equator converge smoothly to a purely radial trajectory at the axis. 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₀ the gas within the cocoon has a pressure P₀, distributed uniformly over a volume V₀. 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 = [(γ + 1)P′/2ρ]^1/2, where P′ is the post-shock 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 derived geometrically. Consider a point on the shock surface with speed |$\boldsymbol {v}$| which travels from an initial angular position (R, θ(R, t)) to a new position (R + dR, θ(R + dR, t + dt)) in time dt. The displacement vector |${\rm d}\boldsymbol {s}$| connecting the initial and final points is normal to the surface and has a magnitude of ds = vdt. If the shock is inclined by an angle ω with respect to the radial direction, the radial component of the displacement is given by dR = dscos ω. We now decompose the displacement vector into two components: the first, |${\rm d}\boldsymbol {s}_1$|⁠, goes from (R, θ) to (R + dR, θ(R + dR, t)) with t held fixed, and the second, |${\rm d}\boldsymbol {s}_2$|⁠, goes from (R + dR, θ(R + dR, t)) to (R + dR, θ(R + dR, t + dt)) with R held fixed. On the first path, R changes by dR = dscos ω, and θ changes by |$(\partial \theta /\partial R) {\rm d}R$| (since t is fixed), so the total path-length is |${\rm d}s_1= {\rm d}s \cos \omega \sqrt{1+ R^2 (\partial \theta /\partial R)^2}$|⁠. On the second path, R is held constant and therefore the path-length is |${\rm d}s_2=R |\partial \theta /\partial t| {\rm d}t= (R/v) |\partial \theta /\partial t| {\rm d}s$|⁠. To complete the derivation, we note that the angle between |${\rm d}\boldsymbol {s}_1$| and |${\rm d}\boldsymbol {s}_2$| is ω, and that |${\rm d}\boldsymbol {s}_1$| is tangent to the surface and perpendicular to |${\rm d}\boldsymbol {s}$|⁠. Therefore, the vectors |${\rm d}\boldsymbol {s}$|⁠, |${\rm d}\boldsymbol {s}_1$|⁠, and |${\rm d}\boldsymbol {s}_2$| form a right triangle such that ds₁/ds₂ = cos ω. Substituting the values of ds₁ and ds₂ in this equation results in |$v = R |\partial \theta /\partial t| [1+ R^2 (\partial \theta /\partial R)^2]^{-1/2}$|⁠.

Equating the two expressions for v, we obtain a partial differential equation for θ(R, t),²

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

(15)

where we have assumed a strong, non-relativistic shock and applied the jump condition ρ′/ρ = v/v′ = (γ + 1)/(γ − 1) for adiabatic index γ.

The KA allows us to replace the post-shock pressure P′(t) in equation (15) by the volume-averaged pressure of the cocoon P_c(t). If losses are unimportant, the total cocoon energy, E₀, can be treated as constant, so that P_c(t) ∝ E₀/V_c(t). We can then write

$$\begin{eqnarray*} P^{\prime }(t) = \dfrac{\lambda (\gamma -1)E_0}{V_{\rm c}(t)}, \end{eqnarray*} $$

(16)

where V_c(t) is the volume of the cocoon at time t, and

$$\begin{eqnarray*} \lambda \equiv \left(\dfrac{E_{\rm th}}{E_0}\right)\left(\dfrac{P^{\prime }}{P_{\rm c}}\right) \end{eqnarray*} $$

(17)

is an order-unity dimensionless parameter which depends on the ratio of the thermal energy Eth to the total energy E₀, as well as the ratio of the post-shock pressure to the average pressure. The essential assumption of the KA is that λ 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 post-shock region in detail. For simplicity, we adopt λ = 1 and P′ = P_c = E₀/(3V_c) in subsequent calculations. The actual value of λ can be estimated from numerical simulations (see Section 5).

Equation (15) can be simplified by noting that P′ depends only on t, while ρ depends only on R. Following K92, we define a dimensionless parameter x = x(t):

$$\begin{eqnarray*} x(t) \equiv \int _{t_0}^t \sqrt{\dfrac{\gamma +1}{2} \dfrac{P_{\rm c}(t^{\prime })}{\rho _0 a^2}} {\rm d}t^{\prime }. \end{eqnarray*} $$

(18)

The parameter x serves as a dimensionless replacement for the time. It is initialized to x = 0 at t = t₀ and thereafter increases monotonically with t. Using the fact that the initial shock velocity along the axis is v₀ = [(γ + 1)P₀/2ρ₀]^1/2, equation (18) can be rewritten as

$$\begin{eqnarray*} x(t) = \int _{t_0}^t \sqrt{\dfrac{P_{\rm c}(t^{\prime })}{P_0}} \dfrac{v_0 {\rm d}t^{\prime }}{a}. \end{eqnarray*} $$

(19)

We then see that x is essentially the product of two dimensionless ratios: the ratio [P_c(t)/P₀]^1/2 comparing the pressure at time t to the initial pressure, and the ratio t/(a/v₀) 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 equation (18) converges or diverges as t → ∞ (see also Section 3.4). Since we expect to recover the ST solution at late times, with R ∝ t^{2/(5 − α)} and P_c ∝ R⁻³ ∝ t^{−6/(5 − α)}, we find that the integral diverges for |$\alpha \leq 2$|⁠. The transition to the sideways expansion phase discussed in Section 2 occurs at x ∼ b/a, while the transition to the quasi-spherical phase occurs when x ∼ 1.

Substituting equations (11) and (18) into equation (15), we obtain

$$\begin{eqnarray*} \dfrac{1}{a^2} \left(\dfrac{\partial \theta }{\partial x} \right)^2 = \left(\dfrac{R}{a}\right)^\alpha \left[ \left(\dfrac{\partial \theta }{\partial R} \right)^2 + \dfrac{1}{R^2} \right], \end{eqnarray*} $$

(20)

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_c = P₀(V₀/V_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 α, 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 ρ(R) = ρ₀, each point along the shock surface travels with the same velocity, v(t) = [(γ + 1)P_c(t)/2ρ₀]^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) = ∫v(t)dt/a in the constant density case. In other words, when α = 0, the distance travelled along each trajectory at time t is simply ax(t). Thus, a particle beginning at (R_i, θ_i) at time t₀ travels a distance ax along a straight line to reach a new position (R, θ) at time t, as illustrated in Fig. 3.

$As in Fig. 2, but for a uniform density ρ = ρ0. Particles follow a straight-line trajectory from (Ri, θi) to (R, θ), covering a distance ax. A triangle of side lengths Ri, R, and ax is formed by connecting the starting location, the new location, and the origin. The interior angle opposite to R is given by $\pi -(\chi _i-\theta _i)$.$

Figure 3.

As in Fig. 2, but for a uniform density ρ = ρ₀. Particles follow a straight-line trajectory from (R_i, θ_i) to (R, θ), covering a distance ax. A triangle of side lengths R_i, R, and ax is formed by connecting the starting location, the new location, and the origin. The interior angle opposite to R is given by |$\pi -(\chi _i-\theta _i)$|⁠.

Open in new tab Download slide

Because the trajectories are straight lines in this case, R and θ 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, θ_i), and (R, θ), as seen in Fig. 3. From the Law of Cosines,

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

(21)

where we have defined ω_i as the angle between the initial direction of motion and the radial direction, i.e.

$$\begin{eqnarray*} \omega _i \equiv \chi _i-\theta _i. \end{eqnarray*} $$

(22)

The cylindrical coordinates r and z are obtained via r = R_isin θ_i + axsin χ_i and z = R_icos θ_i + axcos χ_i. After some manipulation using angle sum and difference identities, θ = tan ⁻¹(r/z) becomes

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

(23)

If θ_i is held fixed, equations (21) and (23) give the trajectory of a particle that started at θ_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 θ along a trajectory stemming from θ_i:

$$\begin{eqnarray*} R(\theta ,\theta _i) = \dfrac{R_i}{\cos (\theta -\theta _i)-\sin (\theta -\theta _i)\cot \omega _i}. \end{eqnarray*} $$

(24)

As expected, this describes a straight line that passes through (R_i, θ_i) and forms an angle χ_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_c and r_c, respectively. The height is defined as the radius at the pole,

$$\begin{eqnarray*} z_{\rm c}(x) \equiv R(x,\theta _i=0), \end{eqnarray*} $$

(25)

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

$$\begin{eqnarray*} r_{\rm c}(x) \equiv \max \lbrace R (x, \theta _i) \sin [\theta (x,\theta _i)]\rbrace \end{eqnarray*} $$

(26)

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 semimajor axis a and semiminor axis b, we have (R_i, χ_i) = (a, 0) at θ_i = 0, and |$(R_i,\chi _i)=(b,\pi /2)$| at |$\theta _i=\pi /2$|⁠. Then, from equation (21), z_c = a(x + 1) and r_c = a(x + b/a), so that the ratio of the cocoon’s width to its height is

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} = \dfrac{x+b/a}{x+1} = 1-\dfrac{a}{z_{\rm c}}\left(1-b/a\right). \end{eqnarray*} $$

(27)

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

$$\begin{eqnarray*} \epsilon _{\rm eff}\approx \left[2 \left(1-b/a\right) a/z_{\rm c}\right]^{1/2} \end{eqnarray*} $$

(28)

to leading order in a/z_c.

As x → ∞, the outflow becomes gradually wider and r_c/z_c → 1 while ϵ_eff → 0. In fact, we see from equation (21) that for large x, R ≈ ax regardless of θ_i, indicating that the cocoon becomes spherical. Taking x ∝ R and P_c ∝ V_c ∝ R⁻³ in equation (18), we see that the usual R ∝ t^2/5 behaviour for a blast wave in a constant density is recovered.

The solution for α = 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. Secondly, 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_c = 2a, equation (27) shows that its width is already r_c > z_c/2, regardless of how narrow it was to start. As we will see, these basic features do not necessarily hold true at higher α.

Figure 4.

The shape of the cocoon (thick lines) when z_c/a = 1, 2, 4, 6, 8, and 10, for a constant density profile with α = 0 and a/b = 5. Several trajectories for different starting points on the cocoon surface are also shown as thin lines.

Open in new tab Download slide

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 |χ_i − θ_i|. The cocoon becomes more deformed in regions where |χ_i − θ_i| is larger, eventually developing a bulge near the location where |χ_i − θ_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_c; the width in the equatorial plane, r_eq; and the coordinates R_b and θ_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_b = R_bcos θ_b above the equatorial plane and at a distance r_b = R_bsin θ_b from the axis. The cocoon width is then r_c = max(r_b, r_eq). This configuration is illustrated in Fig. 5.

Figure 5.

The quantities z_c, r_b, z_b, and r_eq that characterize the cocoon shape. We define the cocoon width r_c to be the larger of r_b and r_eq (in this case, r_c = r_b).

Open in new tab Download slide

The determination of r_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_eq is only accurate to within a factor of a few, it may none the less be useful as a lower bound.

For α ≠ 2, K92 showed that equation (20) can be solved by means of a conformal coordinate transformation. We define

$$\begin{eqnarray*} k_\alpha \equiv \dfrac{2-\alpha }{2}, \end{eqnarray*} $$

(29)

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

$$\begin{eqnarray*} \dfrac{1}{k_\alpha ^2}\left(\dfrac{\partial \phi }{\partial x} \right)^2 = \left(\dfrac{\partial \phi }{\partial \xi } \right)^2 + \dfrac{1}{\xi ^2}, \end{eqnarray*} $$

(30)

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 ϕ_i = k_αθ_i.

The angle ψ_i between the axis and the initial direction of motion in ξ–ϕ coordinates can be derived as follows. We first observe that the quantity dln R/dθ = dln ξ/dϕ 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, θ_i), the initial condition |${\rm d}\ln R/{\rm d}\theta =\cot (\chi _i-\theta _i)$| must be satisfied at θ = θ_i. By the same logic, |${\rm d}\ln \xi /{\rm d}\phi =\cot (\psi _i-\phi _i)$| when evaluated at ϕ = ϕ_i. Thus, we have |$\cot (\psi _i-\phi _i)=\cot (\chi _i-\theta _i)=\cot \omega _i$|⁠, and therefore ψ_i = ϕ_i + ω_i (up to a phase that we choose to be zero).

Having determined the initial conditions in ξ–ϕ space, the solution of equation (30) now proceeds in the same fashion as in the α = 0 case. Following the procedure outlined in the previous section yields

$$\begin{eqnarray*} \xi (x,\phi _i) = \left[ \xi _i^2 + k_\alpha ^2 x^2 + 2 k_\alpha x \xi _i \cos (\psi _i-\phi _i) \right]^{1/2} \end{eqnarray*} $$

(31)

and

$$\begin{eqnarray*} \phi (x,\phi _i) = \phi _i + \tan ^{-1}\left(\dfrac{x \sin (\psi _i-\phi _i)}{\xi _i + x \cos (\psi _i-\phi _i)}\right). \end{eqnarray*} $$

(32)

Transforming back to R–θ coordinates, and using the fact that ψ_i − ϕ_i = ω_i (see above) we obtain:

$$\begin{eqnarray*} R(x,\theta _i) = \left[R_i^{2k_\alpha }+a^{2k_\alpha } k_\alpha ^2 x^2 + 2 k_\alpha x (R_i a)^{k_\alpha } \cos \omega _i \right]^{1/2k_\alpha } \end{eqnarray*} $$

(33)

$$\begin{eqnarray*} \theta (x,\theta _i) = \theta _i + \dfrac{1}{k_\alpha } \tan ^{-1} \left(\dfrac{ k_\alpha x \sin \omega _i}{ \left(R_i/a\right)^{k_\alpha } + k_\alpha x \cos \omega _i} \right). \end{eqnarray*} $$

(34)

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 θ_i:

$$\begin{eqnarray*} R(\theta ,\theta _i) = \dfrac{R_i}{\left\lbrace \cos \left[k_\alpha (\theta -\theta _i)\right] - \cot \omega _i \sin \left[k_\alpha (\theta -\theta _i)\right]\right\rbrace ^{1/k_\alpha }}. \end{eqnarray*} $$

(35)

Note that when R_i = a, θ_i = 0, and ω_i = χ_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_b and θ_b describe a particular point on the cocoon’s surface, they must be functions of time alone, i.e. R_b = R_b(x) and θ_b = θ_b(x). The trajectory passing through the point (R_b(x), θ_b(x)) at a given time x can be traced back to a specific initial angular position θ_b0(x), as shown in Fig. 6. That is to say, for each time x, there exists a corresponding initial angle θ_b0(x) such that R(x, θ_b0(x)) = R_b(x) and θ(x, θ_b0(x)) = θ_b(x). The relationship between x and θ_b0 is derived in Appendix B, but the exact functional form is not important for this discussion. For now, it suffices to say that there is a one-to-one mapping between θ_b0 and t. Therefore, θ_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 θ_{b0, min} as t → ∞.

The shape of the shock at three times – x1 (blue), x2 (red), and x3 (magenta) – taking α = 2 and b/a = 0.3 as an illustrative example. The coloured squares mark the location of the bulge where the velocity is parallel to the equator. The angular coordinates of the points A’, B’, and C’ are, respectively, θb(x1), θb(x2), and θb(x3). At each time, the dashed trajectory passing through the bulge can be traced to a unique point on the initial surface, as indicated by the filled coloured circles along the black curve. The points labelled A, B, and C have initial angular positions of θi = θb0(x1), θb0(x2), and θb0(x3) respectively.

Figure 6.

The shape of the shock at three times – x₁ (blue), x₂ (red), and x₃ (magenta) – taking α = 2 and b/a = 0.3 as an illustrative example. The coloured squares mark the location of the bulge where the velocity is parallel to the equator. The angular coordinates of the points A’, B’, and C’ are, respectively, θ_b(x₁), θ_b(x₂), and θ_b(x₃). At each time, the dashed trajectory passing through the bulge can be traced to a unique point on the initial surface, as indicated by the filled coloured circles along the black curve. The points labelled A, B, and C have initial angular positions of θ_i = θ_b0(x₁), θ_b0(x₂), and θ_b0(x₃) respectively.

Open in new tab Download slide

Now, since every trajectory meets the surface at a right angle, the trajectory which intersects the surface at (R_b(x), θ_b(x)) must be parallel to the equator at that point (see Fig. 6). This means that dln R/dθ = tan θ_b must be satisfied when evaluated at R = R_b and θ = θ_b. Differentiating equation (35) and setting θ_i = θ_b0(x) leads to the equations

$$\begin{eqnarray*} \theta _{\rm b}(x) = \theta _{\rm b0}(x) + \dfrac{2}{\alpha }\left[\pi /2-\chi _{\rm b0}(x)\right] \end{eqnarray*} $$

(36)

$$\begin{eqnarray*} R_{\rm b}(x) = R_{\rm b0}(x) \left(\dfrac{\sin [\chi _{\rm b0}(x)-\theta _{\rm b0}(x)]}{\cos [\theta _{\rm b}(x)]}\right)^{1/k_\alpha }, \end{eqnarray*} $$

(37)

where R_b0(x) ≡ R_i(θ_i = θ_b0(x)) and χ_b0(x) ≡ χ_i(θ_i = θ_b0(x)). As the cocoon expands, the position of the bulge traces out a path described by the parametric curve (R_b(θ_b0), θ_b(θ_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_c to its height z_c. To make a direct comparison between r_c and z_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_c and z_c in terms of x, which we then solved to obtain r_c as an explicit function of z_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 θ_i = χ_i = 0 and R_i = a in equation (33),

$$\begin{eqnarray*} z_{\rm c}(x) = a (1+k_\alpha x)^{1/k_\alpha }, \end{eqnarray*} $$

(38)

in general there is not a way to express r_c(x) explicitly in a similar fashion. The underlying reason that z_c can be described by a simple function of x, but r_c cannot, 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_c(x) exists in the constant-density case, since in that case the width is always determined by the ejecta travelling 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_c and z_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 θ_b0. This turns out to be a good choice, in the sense that the important quantities characterizing the cocoon’s shape (z_c, r_b, and z_b; see Fig. 5) can be written explicitly in terms of θ_b0 and the quantities R_b0 and χ_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_c(θ_b0) and the width r_c(θ_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_c(z_c) relating the height and width. None the less, it is possible to simplify the exact equations to obtain an approximate relation between r_c and z_c in each of the dynamical regimes discussed in Section 2 (see Appendix C). Expressions relating r_eq and z_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_c = r_eq. We find that the evolution of r_c depends on the relationship between α 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 α > 2b²/a², 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 ≫ 1, in which case α ≫ 2b²/a² 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_c − b ≈ (2b²/αa²)b ≪ b. The cocoon’s width is therefore governed by the behaviour of the bulge, with r_c = r_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 θ_b0 ≫ 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 θ. By the time the cocoon’s width has doubled, the bulge has crossed most of the surface, with θ_b decreasing from |$\sim \pi /2$| to ∼b/a.
In the sideways expansion phase, the particles at the location of the bulge originated from b²/a² ≪ θ_b0 ≪ b/a. On this part of the surface, the r- and z-components of the initial velocity are comparable, with |$\chi _{\rm b0}\gtrsim \pi /4$|⁠. The particles in the bulge have had to change their direction of motion by ≲45 deg, moving farther from the axis in the process. By the time a particle starting from θ_b0 passes through the bulge, its angular position has increased to θ_b ≫ θ_b0. Because the width of the outflow increases in this phase, while the height stays about the same (see Section 2), θ_b now grows over time. This phase ends when the height and width of the cocoon become comparable, at which point θ_b ∼ 1.
Finally, during the quasi-spherical phase, the width is governed by material that came from θ_b0 ≪ b²/a². 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 ∼90 deg. As in the previous phase, θ_b ≫ θ_b0. As time goes on, the value of θ_b continues to gradually increase. Eventually, the evolution becomes scale-free, with all length-scales ∝ z_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 α, f_α, and g_α, such that
$$\begin{eqnarray*} \begin{array}{ll}f_\alpha \equiv \lim _{t \rightarrow \infty } (r_{\rm c}/z_{\rm c}) \\ g_\alpha \equiv \lim _{t \rightarrow \infty } (r_{\rm eq}/z_{\rm c}) \end{array} . \end{eqnarray*} $$
(39)
In addition, we define the maximum angle reached by the bulge as
$$\begin{eqnarray*} \theta _{\rm b \infty }\equiv \lim _{t\rightarrow \infty } \theta _{\rm b}. \end{eqnarray*} $$
(40)
If the outflow becomes spherical, then f_α = g_α = 1, and the bulge moves to the equator so |$\theta _{\rm b \infty }= \pi /2$|⁠. Otherwise, we have f_α < 1 and g_α < 1, with f_α ≠ g_α, and |$\theta _{\rm b \infty }\lt \pi /2$|⁠.

In both the sideways expansion and quasi-spherical phases, we find that θ_b0 is negligible compared to χ_b0 and θ_b. This property is not unique to the bulge point. In fact, as long as b ≪ a and α ≫ 2b²/a², the conditions θ_i ≪ θ and θ_i ≪ χ_i hold over all parts of the surface where θ ≲ x. Particles with θ_i ≪ 1 also satisfy R_i ≈ a, since they started out near the tip of the cocoon. Therefore, we can neglect θ_i and set R_i/a ≈ 1 in equations (33) and (34) to obtain an approximate formula for the shape,

$$\begin{eqnarray*} R(\theta ,x) \approx a\left[\cos \theta + k_\alpha x \sqrt{1- \left(\dfrac{\sin k_\alpha \theta }{k_\alpha x}\right)^2}\right]^{1/k_\alpha }, \end{eqnarray*} $$

(41)

which is valid for initially narrow cocoons out to angles of θ ≲ x. During the sideways expansion phase, it can be shown that x ∼ θ_b (see Appendix B); therefore, in this regime, equation (41) describes the shape between the axis and the bulge, for |${0 \leq \theta \lesssim \theta _{\rm b}}$|⁠. On the other hand, during the quasi-spherical phase x ≳ 1 and equation (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 K92 did not explicitly give an equation for R(θ, 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 shallow (α < 2) or steep (α > 2). We discuss each regime in turn below, along with the special case α = 2. As before, we treat only the specific case of an initially ellipsoidal shock.

3.2.1 A shallow external density profile (0 < α < 2)

For α < 2, k is positive, and the integral in equation (18) diverges so x → ∞ as t → ∞. 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 ∝ t^{(2 − α)/(5 − α)}. Thus, as expected, we recover a blast wave evolution with z_c ∝ t^{2/(5 − α)} at late times.

As discussed above, interactions at the equator occur in a non-uniform external density when α > 2b²/a². Which regions on the cocoon surface eventually experience a collision? For 2b²/a² < α < 2, it turns out that θ_b0 has a minimum value θ_{b0, min} > 0. Trajectories stemming from θ_{b0, min} are asymptotically parallel to the equator, i.e. they have |$\theta \rightarrow \pi /2$| as x → ∞. Particles originating from θ_i < θ_{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 θ_i > θ_{b0, min} eventually undergo a collision at the equator. The value of θ_{b0, min} is derived in Appendix C.

As time goes on, θ_b0 decreases (see Fig. 6) and becomes progressively closer to θ_{b0, min}. To understand the sideways spreading at late times, we can therefore study the behaviour of equations (36) and (37) as θ_b0 → θ_{b0, min}. For z_c ≫ a, we find (see Appendix C) that the ratio of the cocoon’s width to its height is approximated by

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \approx 1-\dfrac{A_\alpha }{\zeta (z_{\rm c})}, \end{eqnarray*} $$

(42)

where

$$\begin{eqnarray*} \zeta (z_{\rm c}) \equiv \dfrac{(z_{\rm c}/a)^{|k_\alpha |}-1}{|k_\alpha |}, \end{eqnarray*} $$

(43)

and

$$\begin{eqnarray*} A_\alpha = \dfrac{2}{k_\alpha ^2} \sin ^2(k_\alpha \pi /4) \end{eqnarray*} $$

(44)

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

$$\begin{eqnarray*} \epsilon _{\rm eff}\approx \left(2 A_\alpha /\zeta \right)^{1/2}. \end{eqnarray*} $$

(45)

Since r_c/z_c → 1 and |$\theta _{\rm b}\rightarrow \pi /2$| as z_c → ∞, we infer that the outflow becomes spherical in the 0 < α < 2 case, and therefore f_α = g_α = 1. However, the process of spherization is slower than in the constant-density case, with r_c/z_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 α = 1 and b/a = 0.2. The trajectory stemming from θ_{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 α > 2b²/a² a portion of cocoon material with θ_i > θ_{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_c, z_b/z_c steadily decreases and the relative position of the bulge gradually approaches the equator at late times.

Figure 7.

As in Fig. 4, but for α = 1. The dotted line separates trajectories that eventually reach the equator and collide with cocoon material from the other side, from those that proceed freely to infinity without interacting. The dashed line shows the curve R_b(θ_b), which always passes through the cocoon surface at its widest point.

Open in new tab Download slide

3.2.2 A wind-like external density profile (α = 2)

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

$$\begin{eqnarray*} \left(\dfrac{\partial \theta }{\partial x}\right)^2 = R^2\left(\dfrac{\partial \theta }{\partial R}\right)^2+1. \end{eqnarray*} $$

(46)

This equation admits a solution of the form θ²(x, R) = θ_x(x) + θ_R(R) (Kompaneets 1960). By separation of variables, we find

$$\begin{eqnarray*} \theta ^2 = {(x+\kappa _1)^2 - (\ln R+\kappa _2)^2 + \kappa _1^2-\kappa _2^2}, \end{eqnarray*} $$

(47)

where κ₁ = κ₁(θ_i) and κ₂ = κ₂(θ_i) depend on the original shape. The initial conditions θ(x = 0, θ_i) = θ_i, R(x = 0, θ_i) = R_i, and dln R_i/dθ_i = −tan ω_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

$$\begin{eqnarray*} (x+\kappa _1)^2 = \theta ^2 + [\ln (R/R_i)+\theta _i \cot \omega _i]^2 -\kappa _1^2+ \kappa _2^2. \end{eqnarray*} $$

(48)

To derive an equation for the trajectories, we note that since the shock is a surface of constant x, |$\nabla x = (\partial x/\partial R) \hat{\boldsymbol {R}} + (1/R)(\partial x/\partial \theta) \hat{\boldsymbol {\theta }}$| is normal to the cocoon surface and parallel to the particle paths. (Here, |$\hat{\boldsymbol {R}}$| and |$\hat{\boldsymbol {\theta }}$| are the usual unit vectors in polar coordinates.) The infinitesimal displacement vector |${\rm d}\boldsymbol {s}=({\rm d}R) \hat{\boldsymbol {R}}+(R {\rm d}\theta) \hat{\boldsymbol {\theta }}$| is also parallel to the trajectories. Now, because ∇x and |${\rm d}\boldsymbol {s}$| are parallel vectors, the ratios of their R- and θ-components must be equal. We therefore obtain

$$\begin{eqnarray*} R \dfrac{{\rm d}\theta }{{\rm d}R} = \dfrac{(1/R)(\partial x/\partial \theta)}{\partial x/\partial R} = \dfrac{\theta }{\ln (R/R_i) + \theta _i \cot \omega _i}, \end{eqnarray*} $$

(49)

where the second equality comes from evaluating the partial derivatives of equation (48) and substituting. Integrating subject to the initial condition R = R_i at θ = θ_i, we find that the trajectories satisfy

$$\begin{eqnarray*} R(\theta ,\theta _i) = R_i {\rm e}^{(\theta -\theta _i) \cot \omega _i}. \end{eqnarray*} $$

(50)

The α = 2 case is unique in the sense that |${\rm d}\ln R/{\rm d}\theta =\cot \omega _i$| is constant along particle paths.

As in Section 3.2, R and θ can be written as parametric functions of x and θ_i. After substitution of equations (48) into (50), we eliminate either θ or R to obtain, respectively,

$$\begin{eqnarray*} R(x,\theta _i) = R_i {\rm e}^{x \cos \omega _i} \end{eqnarray*} $$

(51)

and

$$\begin{eqnarray*} \theta (x,\theta _i) = \theta _i + x \sin \omega _i. \end{eqnarray*} $$

(52)

As before, we obtain an approximate expression for the shape of the shock in the sideways expansion and quasi-spherical regimes by neglecting θ_i and adopting R_i ≈ a. This gives

$$\begin{eqnarray*} R(\theta ,x) \approx a {\rm e}^{\sqrt{x^2-\theta ^2}}. \end{eqnarray*} $$

(53)

Setting θ_i = χ_i = 0 and R_i = a in equation (51) gives the cocoon height,

$$\begin{eqnarray*} z_{\rm c}= a {\rm e}^{x}. \end{eqnarray*} $$

(54)

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

$$\begin{eqnarray*} \theta _{\mathrm{ b}}=\pi /2-(\chi _{\rm b0}-\theta _{\rm b0}) \end{eqnarray*} $$

(55)

$$\begin{eqnarray*} R_{\rm b}= R_{\rm b0}{\rm e}^{(\theta _{\rm b}-\theta _{\rm b0}) \cot (\chi _{\rm b0}-\theta _{\rm b0})} . \end{eqnarray*} $$

(56)

Equations (50)–(56) agree with equations (33)–(37) and (41) in the limit of α → 2 and k_α → 0.

Overall, the late-time evolution for α = 2 is similar to the α < 2 case. One notable difference is that θ_{b0, min} = 0 for k_α = 0, so θ_i > θ_{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

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \approx 1 - \dfrac{\pi ^2}{8\ln (z_{\rm c}/a)}, \end{eqnarray*} $$

(57)

and the eccentricity is

$$\begin{eqnarray*} \epsilon _{\rm eff}\approx \left[\pi ^2/4\ln (z_{\rm c}/a)\right]^{1/2}. \end{eqnarray*} $$

(58)

Equations (57) and (58) can also be derived by taking the limit of equations (42)–(45) as k_α → 0. Although the cocoon once again becomes spherical, the logarithmic dependence on z_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 α = 2 are shown in Fig. 8, using the same values as before for z_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.

Figure 8.

As in Figs 4 and 7, but for α = 2.

Open in new tab Download slide

3.2.3 A steep external density profile (2 < α < 3)

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

$$\begin{eqnarray*} \dfrac{1}{|k_\alpha |} = \int _{t_0}^{\infty }\sqrt{\dfrac{(\gamma +1)P_{\rm c}(t^{\prime })}{2\rho _0 a^2}} {\rm d}t^{\prime } . \end{eqnarray*} $$

(59)

Then, subtracting equation (18) from equation (59) and multiplying through by |k_α|, we obtain

$$\begin{eqnarray*} 1 + k_\alpha x = |k_\alpha | \int _t^{\infty }\sqrt{\dfrac{(\gamma +1)P_{\rm c}(t^{\prime })}{2\rho _0 a^2}} {\rm d}t^{\prime }. \end{eqnarray*} $$

(60)

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_αx) ∝ t^{(2 − α)/(5 − α)}. Therefore the usual blast wave scaling z_c ∝ t^{2/(5 − α)} is also recovered for α > 2.

A peculiar feature of the solution for α > 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 → ∞, we have θ_b0 → θ_{b0, min} and θ_b → θ_b∞. Since θ_{b0, min} = 0 in this case, we take θ_b0 → 0 in equation (36) to obtain |$\theta _{\rm b \infty }= \pi /\alpha \lt \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

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \approx f_\alpha \left(1-\dfrac{A_\alpha }{\zeta } \right), \end{eqnarray*} $$

(61)

where ζ is defined as before, and

$$\begin{eqnarray*} f_\alpha = \left[\sin \left(\dfrac{\pi }{\alpha }\right)\right]^{\alpha /(\alpha -2)}. \end{eqnarray*} $$

(62)

In this regime the parameter A_α takes on a different form,

$$\begin{eqnarray*} A_\alpha = \dfrac{1}{2k_\alpha ^2} \cot ^2(\pi /\alpha). \end{eqnarray*} $$

(63)

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

$$\begin{eqnarray*} \epsilon _{\rm eff}\approx \ \left(1-f_\alpha ^2 +2 f_\alpha ^2 A_\alpha /\zeta\right)^{1/2}. \end{eqnarray*} $$

(64)

Unlike the previous cases, ϵ_eff no longer goes to zero, but instead approaches |$(1-f_\alpha ^2)^{1/2}$| as z_c becomes large.

Considering equation (41) in the limit x → 1/|k_α| leads to an asymptotic expression for the cocoon’s shape: as z_c tends to infinity, the curve bounding the cocoon is increasingly well-approximated by

$$\begin{eqnarray*} R(\theta) \approx z_{\rm c}\left[\cos (k_\alpha \theta)\right]^{-1/k_\alpha }. \end{eqnarray*} $$

(65)

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 = ±a (see e.g. K92).⁴ Setting |$\theta =\pi /2$| in equation (65) and comparing with formula (39), we find

$$\begin{eqnarray*} g_\alpha =\left[\cos (k_\alpha \pi /2)\right]^{-1/k_\alpha }. \end{eqnarray*} $$

(66)

Whereas g_α = f_α = 1 for density profiles flatter than r⁻² because the cocoon ultimately becomes a sphere, steeper density profiles have g_α < f_α < 1 instead. Therefore, r_eq < r_c < z_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 α = 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.

Figure 9.

As in Figs 4, 7, and 8, but for α = 3. The dotted line marks the trajectory of the first particles to reach the equator.

Open in new tab Download slide

We emphasize again that the above results are based on an idealized model that ignores heating due to collisions between ejecta from opposite hemispheres at the equator. 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 relatively 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 ≳45 deg. In comparison, the trajectories in the α = 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 α, there can be a region inside the cocoon, located along the equator, 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 the interaction between ejecta colliding at the equator 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_c ≳ 2b), the velocity becomes comparable in the forwards and sideways directions, and the change in width (r_c − b ≈ r_c) is about the same as the change in height (z_c − a). This approximation remains valid throughout the sideways expansion phase, i.e. while z_c − a ≪ a. Finally, in the quasi-spherical phase, the evolution depends on α, with the ratio of width to height given by equations (42), (57), or (61), respectively, for α < 2, α = 2, or α > 2. Combining these results, and assuming b ≪ a, we then have

$$\begin{eqnarray*} r_{\rm c}\simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b, & z_{\rm c}-a \ll b \\ z_{\rm c}-a, & b \ll z_{\rm c}-a \ll a \\ f_\alpha z_{\rm c}\left(1-A_\alpha /\zeta \right), & z_{\rm c}\gg 2a \end{array}\right.. \end{eqnarray*} $$

(67)

Note that the intermediate regime b ≪ z_c − a ≪ a may not be realized if b/a is not particularly small. In Table 1, we give the appropriate values of f_α, A_α, and ζ in the various regimes, and also summarize some important asymptotic characteristics of the cocoon.

Table 1.

Open in new tab

Asymptotic cocoon properties for the regimes discussed in Section 3.

	\|$\alpha \leq 2b^2/a^2$\|	2b²/a² < α < 2	α = 2	2 < α < 3
	\|$(k_\alpha \geq 1-b^2/a^2)$\|	(0 < k_α < 1 − b²/a²)	(k_α = 0)	(− 1/2 < k_α < 0)
f_α	1			\|$\left[\sin (\pi /\alpha)\right]^{\alpha /(\alpha -2)}$\|
A_α	\|$\left(2/k_\alpha ^2\right) \sin ^2(k_\alpha \pi /4)$\|		\|$\pi ^2/8$\|	\|$\left(1/2k_\alpha ^2\right) \cot ^2(\pi /\alpha)$\|
ζ	\|$(1/k_\alpha)\left[(z_{\rm c}/a)^{k_\alpha }-1\right]$\|		ln (z_c/a)	\|$(-1/k_\alpha)\left[(z_{\rm c}/a)^{-k_\alpha }-1\right]$\|
Bulge forms?	No	Yes
Asymptotic shape	Sphere			‘Peanut’

	\|$\alpha \leq 2b^2/a^2$\|	2b²/a² < α < 2	α = 2	2 < α < 3
	\|$(k_\alpha \geq 1-b^2/a^2)$\|	(0 < k_α < 1 − b²/a²)	(k_α = 0)	(− 1/2 < k_α < 0)
f_α	1			\|$\left[\sin (\pi /\alpha)\right]^{\alpha /(\alpha -2)}$\|
A_α	\|$\left(2/k_\alpha ^2\right) \sin ^2(k_\alpha \pi /4)$\|		\|$\pi ^2/8$\|	\|$\left(1/2k_\alpha ^2\right) \cot ^2(\pi /\alpha)$\|
ζ	\|$(1/k_\alpha)\left[(z_{\rm c}/a)^{k_\alpha }-1\right]$\|		ln (z_c/a)	\|$(-1/k_\alpha)\left[(z_{\rm c}/a)^{-k_\alpha }-1\right]$\|
Bulge forms?	No	Yes
Asymptotic shape	Sphere			‘Peanut’

Table 1.

Open in new tab

Asymptotic cocoon properties for the regimes discussed in Section 3.

	\|$\alpha \leq 2b^2/a^2$\|	2b²/a² < α < 2	α = 2	2 < α < 3
	\|$(k_\alpha \geq 1-b^2/a^2)$\|	(0 < k_α < 1 − b²/a²)	(k_α = 0)	(− 1/2 < k_α < 0)
f_α	1			\|$\left[\sin (\pi /\alpha)\right]^{\alpha /(\alpha -2)}$\|
A_α	\|$\left(2/k_\alpha ^2\right) \sin ^2(k_\alpha \pi /4)$\|		\|$\pi ^2/8$\|	\|$\left(1/2k_\alpha ^2\right) \cot ^2(\pi /\alpha)$\|
ζ	\|$(1/k_\alpha)\left[(z_{\rm c}/a)^{k_\alpha }-1\right]$\|		ln (z_c/a)	\|$(-1/k_\alpha)\left[(z_{\rm c}/a)^{-k_\alpha }-1\right]$\|
Bulge forms?	No	Yes
Asymptotic shape	Sphere			‘Peanut’

	\|$\alpha \leq 2b^2/a^2$\|	2b²/a² < α < 2	α = 2	2 < α < 3
	\|$(k_\alpha \geq 1-b^2/a^2)$\|	(0 < k_α < 1 − b²/a²)	(k_α = 0)	(− 1/2 < k_α < 0)
f_α	1			\|$\left[\sin (\pi /\alpha)\right]^{\alpha /(\alpha -2)}$\|
A_α	\|$\left(2/k_\alpha ^2\right) \sin ^2(k_\alpha \pi /4)$\|		\|$\pi ^2/8$\|	\|$\left(1/2k_\alpha ^2\right) \cot ^2(\pi /\alpha)$\|
ζ	\|$(1/k_\alpha)\left[(z_{\rm c}/a)^{k_\alpha }-1\right]$\|		ln (z_c/a)	\|$(-1/k_\alpha)\left[(z_{\rm c}/a)^{-k_\alpha }-1\right]$\|
Bulge forms?	No	Yes
Asymptotic shape	Sphere			‘Peanut’

Fig. 10 compares the limiting expressions in (67) with the full solution for each of the values of α discussed above. We see that when α is not too close to 2, the asymptotic approximation agrees well with the exact solution for z_c ≳ 10a. For α ≈ 2, however, the leading-order approximation given in equation (67) becomes less accurate. The reason for the slow convergence is that, since ζ ∝ ln z_c for α = 2, the contribution from higher order terms in ζ is non-negligible. If greater accuracy is required, the parametric equations presented in Appendix B can be used to calculate r_c/z_c instead.

Figure 10.

The ratio r_c/z_c as a function of cocoon height, for b/a = 0.03 and α = 0, 1, 2, and 3. The dotted and dashed lines show the approximations given in equation (67), respectively, in the limits z_c − a ≪ a and z_c ≫ 2a.

Open in new tab Download slide

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_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_c for t ≪ t_a and t ≫ 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 ≪ t_a), when the cocoon is still roughly elliptical, we approximate the volume as

$$\begin{eqnarray*} V_{\rm c}\approx V_0 (z_{\rm c}/a) (r_{\rm c}/b)^2. \end{eqnarray*} $$

(68)

At late times (t ≫ t_a), we have r_c ∝ z_c and therefore |$V_{\rm c}\propto z_{\rm c}^3$|⁠. In this case we write

$$\begin{eqnarray*} V_{\rm c}\approx (4/3) \pi z_{\rm c}^3 C_\alpha , \end{eqnarray*} $$

(69)

where the constant of proportionality C_α is defined by

$$\begin{eqnarray*} C_\alpha \equiv \lim _{t\rightarrow \infty } \dfrac{3V_{\rm c}}{4 \pi z_{\rm c}^3}. \end{eqnarray*} $$

(70)

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

$$\begin{eqnarray*} C_\alpha = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1, & \alpha \leq 2 \\ \int _0^{\pi /2} \left[ \cos (k_\alpha \theta)\right]^{3/|k_\alpha |} \sin \theta {\rm d}\theta , & \alpha \gt 2 \end{array}\right.. \end{eqnarray*} $$

(71)

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_α as a function of α, and compare it to the value of |$3V_{\rm c}/\left(4\pi z_{\rm c}^3\right)$| found by direct numerical integration of equations (33) and (34) at z_c/a = 10, 100, and 1000. As expected from the results of Section 3.2, the α = 2 case is slowest to converge to the limiting volume, while for α ≲ 1, the asymptotic expression is already a reasonably good approximation for z_c ≳ 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_α by |$f_\alpha ^2$| is reasonably accurate, and has the advantage of being easy to compute via equation (62).

$The volume of the cocoon scaled to $4 \pi z_{\rm c}^3/3$, as a function of α for b/a = 0.1. The blue, red, and yellow curves depict the volume when the cocoon has expanded to 10a, 100a, and 1000a, respectively. The limiting value as zc → ∞ is shown in black: the solid line is the accurate value given by equation (70), while the dotted line is the simple approximation $C_\alpha \simeq f_\alpha ^2$.$

Figure 11.

The volume of the cocoon scaled to |$4 \pi z_{\rm c}^3/3$|⁠, as a function of α for b/a = 0.1. The blue, red, and yellow curves depict the volume when the cocoon has expanded to 10a, 100a, and 1000a, respectively. The limiting value as z_c → ∞ is shown in black: the solid line is the accurate value given by equation (70), while the dotted line is the simple approximation |$C_\alpha \simeq f_\alpha ^2$|⁠.

Open in new tab Download slide

3.4 The age of the cocoon

We now turn to the question of how z_c and r_c evolve in time. To start, we rewrite equation (18) as

$$\begin{eqnarray*} \dfrac{t-t_0}{t_0} \simeq \dfrac{a}{b} \int _0^x \left(\dfrac{V(x^{\prime })}{V_0} \right)^{1/2} {\rm d}x^{\prime }, \end{eqnarray*} $$

(72)

where we used P(x)V(x) = P₀V₀, x(t₀) = 0, γ = 4/3, and t₀ ≃ t_b as given by equation (1). Now, we differentiate equation (38) to find dx = (z_c/a)^{k − 1}d(z_c/a) and change variables in the integral (72)

$$\begin{eqnarray*} \dfrac{t-t_0}{t_0} \simeq \dfrac{a}{b}\int _a^{z_{\rm c}} \left(\dfrac{V(z_{\rm c}^{\prime })}{V_0}\right)^{1/2} \left(\dfrac{z_{\rm c}^{\prime }}{a}\right)^{k-1} \dfrac{{\rm d}z_{\rm c}^{\prime }}{a}. \end{eqnarray*} $$

(73)

Then, substituting for V_c(z_c) using equation (68) (for z_c − a ≪ a) or equation (69) (for z_c − a ≫ a), equation (73) becomes

$$\begin{eqnarray*} \dfrac{t-t_0}{t_0} \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\left(\dfrac{a}{b} \right)^2 \int _a^{z_{\rm c}} \left(\dfrac{r_{\rm c}(z_{\rm c}^{\prime })}{a}\right) \dfrac{{\rm d}z_{\rm c}^{\prime }}{a}, & z_{\rm c}-a \ll a \\ \left(\dfrac{a}{b} \right)^2 C_\alpha ^{1/2} \int _a^{z_{\rm c}} \left(\dfrac{z_{\rm c}^{\prime }}{a} \right)^{(3-\alpha)/2} \dfrac{{\rm d}z_{\rm c}^{\prime }}{a}, & z_{\rm c}-a \gg a \\ \end{array}\right.. \end{eqnarray*} $$

(74)

Finally, we replace r_c using equation (67) and perform the integration to obtain

$$\begin{eqnarray*} \dfrac{t-t_0}{t_0} \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(a/b)(z_{\rm c}/a-1), & z_{\rm c}-a \lesssim b \\ \frac{1}{2}+ \frac{1}{2}(a/b)^2(z_{\rm c}/a-1)^2, & b \lesssim z_{\rm c}-a \ll a \\ \frac{2}{5-\alpha } C_\alpha ^{1/2} (a/b)^2 (z_{\rm c}/a)^{(5-\alpha)/2}, & z_{\rm c}\gg 2a \end{array}\right.. \end{eqnarray*} $$

(75)

The leading constant term in the second expression was chosen to ensure a smooth transition at z_c − a = b.

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

$$\begin{eqnarray*} z_{\rm c}\simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}a+b(t/t_0-1), & t \lesssim 2 t_0\\ a+ b\left(2t/t_0-3\right)^{1/2}, & 2 t_0 \lesssim t \ll t_a \\ a \left(\frac{5-\alpha }{2} C_\alpha ^{-1/2} \frac{b^2}{a^2} \frac{t}{t_0}\right)^{2/(5-\alpha)}, & t \gg t_a \end{array}\right.. \end{eqnarray*} $$

(76)

As expected, the evolution becomes like an ST blast wave at late times with z_c ∝ t^{2/(5 − α)}. Finally, combining equations (67) and (76) gives an approximation for the cocoon width,

$$\begin{eqnarray*} r_{\rm c}\simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}b , & t \lesssim 2 t_0 \\ b\left(2t/t_0-3\right)^{1/2}, & 2t_0 \lesssim t \ll t_a \\ a f_\alpha \left(\frac{5-\alpha }{2} C_\alpha ^{-1/2} \frac{b^2}{a^2} \frac{t}{t_0}\right)^{2/(5-\alpha)}, & t \gg t_a \end{array}\right.. \end{eqnarray*} $$

(77)

In Fig. 12, we plot the analytical estimates in equations (76) and (77), 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.

The cocoon’s height and width versus time, for b/a = 0.03 and α = 1, 2, and 3. The estimates from equations (76) and (77) are overlaid as dotted lines (t/t0 ≪ a2/b2) and dashed lines (t/t0 ≫ a2/b2).

Figure 12.

The cocoon’s height and width versus time, for b/a = 0.03 and α = 1, 2, and 3. The estimates from equations (76) and (77) are overlaid as dotted lines (t/t₀ ≪ a²/b²) and dashed lines (t/t₀ ≫ a²/b²).

Open in new tab Download slide

We conclude this section with an estimate for the time-scale for the cocoon expansion to become spherical. We define this time-scale, t_sph, as the time when the width equals 90 per cent of the height. Note that t_sph is only defined for |$\alpha \leq 2$|⁠, since for α > 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_α which depends on the density profile. Since |$t_a \simeq \sqrt{2} (a/b)^2 t_0$| (equations 1 and 2), we define Q_α via

$$\begin{eqnarray*} Q_\alpha \equiv (\sqrt{2} a^2/b^2)^{-1} t_{\rm sph}/t_0 \simeq t_{\rm sph}/t_a. \end{eqnarray*} $$

(78)

The numerically computed value of Q_α as a function α 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_α (see below). The value of Q_α grows extremely fast with α, with Q_α = 59, 441, and 4.4 × 10⁶, respectively, for α = 0, 1, and 2. Outflows in steep density profiles therefore take a considerably longer time to become spherical.

$The time it takes for the cocoon’s width to become 60 (yellow), 80 (red), and 90 (blue) per cent of its height, as a function of α, for an initially ellipsoidal cocoon with a/b = 10. All times are scaled to the time for the height to double, $t_a \simeq \sqrt{2} (a/b)^2 t_0$. The solid blue line shows the numerically calculated value of Qα, while the dashed line represents the analytical estimate given by equation (79).$

Figure 13.

The time it takes for the cocoon’s width to become 60 (yellow), 80 (red), and 90 (blue) per cent of its height, as a function of α, for an initially ellipsoidal cocoon with a/b = 10. All times are scaled to the time for the height to double, |$t_a \simeq \sqrt{2} (a/b)^2 t_0$|⁠. The solid blue line shows the numerically calculated value of Q_α, while the dashed line represents the analytical estimate given by equation (79).

Open in new tab Download slide

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

$$\begin{eqnarray*} Q_\alpha \simeq \dfrac{t_{\rm sph}}{t_a} \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{\sqrt{2}}{5-\alpha } (1+10 k_\alpha A_\alpha)^{(5-\alpha)/(2-\alpha)} , & \alpha \lt 2 \\ (\sqrt{2}/3) {\rm e}^{15 \pi ^2/8}, & \alpha = 2 \end{array}\right.. \end{eqnarray*} $$

(79)

This approximation overestimates the true value of Q_α, because equation (67) underestimates the true value of r_c/z_c at late times. For α ≲ 1.5, the analytical estimate of Q_α is roughly twice the numerically determined value. However, the approximation is less accurate when α is close to 2, because the series expansion of r_c(z_c) converges slowly in that case, and we have used only the lowest order term. For α = 2, equation (79) overestimates Q_α by a factor of 9. None the less, Fig. 13 and equation (79) both clearly demonstrate that t_sph is much larger than t_a, especially for α ∼ 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_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_w. Secondly, 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_bo when the shock breaks out. The KA model applies up until t = min (t_eq, t_w, t_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 η. Also, we neglect order-unity constants which depend on the density profile. We first consider t_eq. Pressure equilibrium is attained once the cocoon pressure, P_c, becomes comparable to the ambient pressure, P_a. We suppose that the ambient pressure obeys

$$\begin{eqnarray*} P_{\rm a}= P_{\rm {a0}}(\rho /\rho _0)^{\gamma _a} = P_{\rm {a0}}(R/a)^{-\gamma _a \alpha }, \end{eqnarray*} $$

(80)

where γ_a is the adiabatic index of the ambient medium.⁵ Furthermore, we assume that pressure upon choking satisfies P₀ ≫ 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 ≳ P_c(t_a), pressure equilibrium is reached during the sideways expansion phase, while if P_a0 ≲ P_c(t_a), it happens in the quasi-spherical regime. In the former case, the cocoon’s height is given by z_c ∼ a and its volume is |${V_{\rm c}\sim 4 \pi a r_{\rm c}^2/3}$|⁠. The expansion occurs near the tip of the cocoon, where P_a ∼ P_a0. Pressure equilibrium is attained when the width is |$r_{\rm c}\sim a (E_0/4 \pi P_{\rm {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_{\rm a}(z_{\rm c}) a^3)^{1/3}}$|⁠. Substituting equation (80) for P_c(z_c) then leads to |$z_{\rm c}\sim a (E_0/4 \pi P_{\rm {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 ≪ t ≪ t_a) is given by t ∼ t₀(r_c/b)² ∼ t_a(r_c/a)² via equation (77), and the age in the quasi-spherical phase (t ≫ t_a) is given by t ∼ t_a(z_c/a)^{(5 − α)/2} via equation (76). We then arrive at the pressure equilibrium time-scale

$$\begin{eqnarray*} t_{\rm eq}\sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}t_a \eta _{\rm eq}, & b^2/a^2 \ll \eta _{\rm eq}\lesssim 1 \\ t_a \eta _{\rm eq}^{(5-\alpha)/2(3-\gamma _a \alpha)}, & \eta _{\rm eq}\gtrsim 1 \\ \end{array}\right., \end{eqnarray*} $$

(81)

where

$$\begin{eqnarray*} \eta _{\rm eq}\equiv \dfrac{E_0}{4\pi P_{\rm {a0}}a^3}. \end{eqnarray*} $$

(82)

The condition η_eq ≫ (b/a)² comes from our assumption that P₀ ≫ 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_w exceeding the ambient sound speed. Then, once v ∼ v_w, the shock becomes weak and our solution no longer applies. Let us assume a steady wind with ρ ∝ r⁻² and v_w = const. We also require that v ≫ v_w initially. Setting v ∼ (P_c/ρ)^1/2 ∼ v_w, the calculation then proceeds as before. We ultimately arrive at

$$\begin{eqnarray*} t_{\rm w}\sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}t_a \eta _{\rm w}, & b^2/a^2 \ll \eta _{\rm w}\lesssim 1 \\ t_a \eta _{\rm w}^{3/2}, & \eta _{\rm w}\gtrsim 1 \\ \end{array}\right., \end{eqnarray*} $$

(83)

with η_w given by

$$\begin{eqnarray*} \eta _{\rm w}\equiv \dfrac{E_0}{4\pi \rho _0 v_{\rm w}^2 a^3} = \dfrac{E_0}{\dot{M} v_{\rm w}a}, \end{eqnarray*} $$

(84)

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 ≫ b), so that the breakout takes place in either the sideways expansion regime (if b ≪ R_* − a ≪ a) or in the quasi-spherical regime (if R_* − a ≫ a). Then, taking z_c = R_* in equation (75), the breakout time is

$$\begin{eqnarray*} t_{\rm bo}\sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\eta _{\rm bo}^2 t_a, &b/a \ll \eta _{\rm bo}\lesssim 1 \\ \eta _{\rm bo}^{(5-\alpha)/2} t_a, & \eta _{\rm bo}\gtrsim 1 \\ \end{array}\right., \end{eqnarray*} $$

(85)

where

$$\begin{eqnarray*} \eta _{\rm bo}\equiv \dfrac{R_*-a}{a}. \end{eqnarray*} $$

(86)

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 llGRBs); 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_eq, and t_bo.

In a typical long GRB, it takes the jet ∼15 s to drill through the progenitor WR star (Bromberg et al. 2012). Therefore, we expect jets with a duration of ∼ a few seconds to be choked at an appreciable fraction of R_*, i.e. we have a/R_* ∼ a few tenths. In this case η_bo is order-unity and t_bo ∼ t_a. A WR star can be approximated by an n = 3 polytrope, with the ambient density and pressure scaling as ρ ∝ R^−5/2 and P_a ∝ ρ^4/3 ∝ R^−10/3 in the outer parts of the star. As P_aR³ ∝ R^−1/3 changes slowly with radius, the quantity P_a0a³ appearing in equation (82) is roughly constant, regardless of the choking location. Numerical modelling of pre-supernova stars (e.g. Heger, Woosley & Spruit 2005) suggests a pressure of the order of |$2 \times 10^{18}\, \text{erg cm}^{-3}$| at R = 10¹⁰ cm. We then have |${\eta _{\rm eq}\simeq 80 (E_0/10^{51}\, \text{erg})(P_{\rm {a0}}a^3/10^{48}\, \text{erg})^{-1} \gg 1}$|⁠. Thus, t_a ∼ t_bo ≪ t_eq, and the breakout occurs while the KA model is still valid.

The progenitors of llGRBs, on the other hand, may have an extended optically thick envelope containing a mass |$\sim 10^{-2} \, \mathrm{M}_{\odot }$| within a radius ∼100–1000 R_|$\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 η_eq ≳ 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_* ≪ 1. We therefore have η_bo ≫ 1 and a quasi-spherical breakout occurs at t_bo ≫ t_a. Assuming that the jet successfully escapes the star, a is at least ∼10 R_|$\odot$|, so η_bo is at most ∼100. We then find t_a ≪ t_bo < t_eq, so in this case the KA model also applies up until breakout.

In AGNs, bubbles with a size of ∼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_* ∼ 100 kpc, so we have a/R_* ≲ 0.1 and η_bo ≳ 10. Thus, like in the llGRB case, the breakout time is t_bo ≫ 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₀ ∼ 10⁵⁹ erg and a typical ambient pressure of |$P_{\rm a}(10\, \text{kpc})\sim 0.1$| keV cm⁻³ at R = 10 kpc, we obtain |$\eta _{\rm eq}\gtrsim 1.6 (E_0/10^{59}\, \text{erg})[P_{\rm a}(10\, \mathrm{ kpc})/0.1\, \rm{keV\,cm}^{-3}]^{-1}$| for a ≲ 10 kpc. Assuming a relatively flat (α ≃ 1) and isothermal (γ_a = 1) ambient medium, the lower limits on η_bo and η_eq imply t_bo ≳ 100t_a and t_eq ≳ 2t_a, respectively. Furthermore, in the isothermal case, the ratio t_bo/t_eq is independent of a. Dividing equation (85) by equation (81) and applying |$P_{\rm {a0}}= P(10\, \text{kpc})(a/10\, \text{kpc})^{-1}$|⁠, we obtain

$$\begin{eqnarray*} \dfrac{t_{\rm bo}}{t_{\rm eq}} \sim 60 \left(\dfrac{E_0}{10^{59}\, \text{erg}}\right)^{-1} \left(\dfrac{R_*}{100\, \text{kpc}}\right)^2 \left(\dfrac{P_{\rm a}(10\, \text{kpc})}{0.1\, \rm{keV\,cm}^{-3}}\right). \end{eqnarray*} $$

(87)

Therefore, AGN jets choked in the core of galaxy clusters satisfy t_a ≲ t_eq ≪ t_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. 2019), 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_eq (Churazov et al. 2001; Irwin et al. 2019).

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₀) are related to the parameters of the original jet (the luminosity L_j, the duration t_j, and the initial opening angle θ_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 initial jet properties. We will assume that it is possible to estimate the cocoon’s height (z_c), width (r_c), and energy (E₀) 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 ≪ a, as expected when the cocoon is inflated by a relativistic jet. When the observed aspect ratio is z_c/r_c ≫ 1, the association with a jet is not controversial. However, when z_c/r_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_c/z_c ≳ 0.5, and in the planar or sideways expansion phase when r_c/z_c ≲ 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

$$\begin{eqnarray*} t \sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\left(\dfrac{\rho (z_{\rm c}) z_{\rm c}r_{\rm c}^4}{E_0} \right)^{1/2}, & z_{\rm c}/r_{\rm c}- 1 \gg 1 \\ \left(\dfrac{\rho (z_{\rm c}) z_{\rm c}^5}{E_0}\right)^{1/2}, & z_{\rm c}/r_{\rm c}- 1 \ll 1 \\ \end{array}\right., \end{eqnarray*} $$

(88)

where ρ(z_c) is the observed density at z_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_c/z_c is then directly related to the ratio z_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):

$$\begin{eqnarray*} a \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}z_{\rm c}-r_{\rm c}, & z_{\rm c}/r_{\rm c}- 1 \gg 1 \\ z_{\rm c}\left(1+\dfrac{A_\alpha |k_\alpha |}{1-r_{\rm c}/z_{\rm c}}\right)^{-1/|k_\alpha |}, & z_{\rm c}/r_{\rm c}- 1 \ll 1 \\ \end{array}\right.. \end{eqnarray*} $$

(89)

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 ∼ t₀ and b ∼ r_c.

Statistically speaking, though, it is much more likely to catch the system during the sideways expansion phase, since it lasts roughly (a/b)² 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. On average, we expect to observe evolved systems with t ≫ t₀ more often than recently choked systems with t ∼ t₀. However, the possibility that any particular system was quenched in the recent past cannot be ruled out based on the geometry alone. All we can say for certain is that b < r_c and t₀ < t. For a known choking radius a given by equation (89), the degeneracy between the initial width and the choking time is captured by

$$\begin{eqnarray*} t_0 b^{-2} \sim \left(\dfrac{\rho _0 a}{E_0}\right)^{1/2} = \left(\dfrac{\rho (z_{\rm c}) z_{\rm c}}{E_0}\right)^{1/2} \left(\dfrac{a}{z_{\rm c}}\right)^{(1-\alpha)/2}, \end{eqnarray*} $$

(90)

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

Now that we have established how a and b depend on the observables, let us consider how they relate to the properties of the original jet. Suppose the jet is injected relativistically with luminosity L_j and opening angle θ_op, and operates for a duration t_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₀ when the jet is choked, |$\tilde{L} = \text{max}\left\lbrace \left(\frac{L_{\rm j}}{\rho _0 t_0^2 \theta _{\rm op}^4 c^5}\right)^{2/5}, \frac{L_{\rm j}}{\rho _0 t_0^2 \theta _{\rm op}^2 c^5} \right\rbrace$|⁠. 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 β_h, with |$\beta _{\rm h}\approx \tilde{L}^{1/2}$| for a non-relativistic jet, and |$2 \Gamma _\mathrm{ h}^2=\tilde{L}^{1/2}$| for a relativistic head, where |$\Gamma _\mathrm{ 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_j, how long does it take to be choked? Material continues to flow into the cocoon until t₀, 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_j, and therefore t₀ ≃ t_j. If the head is relativistic instead, choking takes longer, and |$t_0 = t_{\rm j}/(1-\beta _{\rm h}) \approx 2 \Gamma _\mathrm{ h}^2 t_{\rm j}$| (Nakar 2015). Thus we have

$$\begin{eqnarray*} t_0 \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}t_{\rm j}, & \tilde{L} \ll 1 \\ \tilde{L}^{1/2} t_{\rm j}, & \tilde{L} \gg 1 \end{array}\right.. \end{eqnarray*} $$

(91)

The distance the head travels in time t₀ is β_hct_j in the former case, and |$2 \Gamma _\mathrm{ h}^2 c t_{\rm j}$| in the latter; in either regime, we can write

$$\begin{eqnarray*} a \simeq \tilde{L}^{1/2} c t_{\rm j}. \end{eqnarray*} $$

(92)

Then, estimating the cocoon width at the time of choking as b ≃ (β_c/β_h)a, where the sideways expansion speed is |$\beta _{\rm c}\simeq \tilde{L}^{1/2} \theta _{\rm op}$||$\left(\tilde{L} \ll \theta _{\rm op}^{-4/3}\right)$|⁠, or |$\beta _{\rm c}\simeq \tilde{L}^{1/8} \theta _{\rm op}^{1/2}$||$\left(\theta _{\rm op}^{-4/3} \ll \tilde{L} \ll 4 \theta _{\rm op}^{-4}\right)$| (Bromberg et al. 2011), we arrive at

$$\begin{eqnarray*} b \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\tilde{L}^{1/2} \theta _{\rm op}c t_{\rm j}, & \tilde{L} \ll 1 \\ \tilde{L} \theta _{\rm op}c t_{\rm j}, & 1 \ll \tilde{L} \ll \theta _{\rm op}^{-4/3} \\ \tilde{L}^{5/8} \theta _{\rm op}^{1/2} c t_{\rm j}, & \theta _{\rm op}^{-4/3} \ll \tilde{L} \ll 4\theta _{\rm op}^{-4} \end{array}\right.. \end{eqnarray*} $$

(93)

We now have three equations (91)–(93) relating the initial conditions of the cocoon (a, b, and t₀) to three quantities describing the jet (⁠|$\tilde{L}$|⁠, t_j, and θ_op). When the choking radius satisfies a ≪ ct₀, we find that the head was non-relativistic, with |$\tilde{L}$| given by

$$\begin{eqnarray*} \tilde{L} \simeq \left(\dfrac{a}{ct_0}\right)^2 \simeq \dfrac{E_0a}{\rho _0 c^2 b^4} \ll 1. \end{eqnarray*} $$

(94)

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

$$\begin{eqnarray*} \begin{array}{lll}L_{\rm j}& \simeq E_0/t_0 & \simeq E_0^{3/2} \rho _0^{-1/2} a^{-1/2} b^{-2} \\ t_{\rm j}& \simeq t_0 & \simeq E_0^{-1/2} \rho _0^{1/2} a^{1/2} b^{2} \\ \theta _{\rm op}& \simeq b/a & \, \, \\ \end{array}. \end{eqnarray*} $$

(95)

When the jet head is relativistic, however, a degeneracy arises because a and t₀ are related trivially by a = ct₀. 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 θ_op follow a closure relation given by

$$\begin{eqnarray*} b/a \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\tilde{L}^{1/2} \theta _{\rm op}, & 1 \ll \tilde{L} \ll \theta ^{-4/3} \\ \tilde{L}^{1/8} \theta _{\rm op}^{1/2}, & \theta ^{-4/3} \ll \tilde{L} \ll 4\theta ^{-4} \\ \end{array}\right.. \end{eqnarray*} $$

(96)

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_jt_j ≃ E₀. For a collimated jet with a relativistic head (⁠|$1 \ll \tilde{L} \ll \theta ^{-4/3}$|⁠), we obtain the additional constraint

$$\begin{eqnarray*} t_{\rm j}\theta _{\rm op}^{-1} \simeq (a/b) t_0 \simeq E_0^{-1/2} \rho _0^{1/2} a^{3/2} b, \end{eqnarray*} $$

(97)

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

$$\begin{eqnarray*} t_{\rm j}\theta _{\rm op}^{-2} \simeq (a/b)^4 t_0 \simeq E_0^{-1/2} \rho _0^{1/2} a^{9/2} b^{-2}. \end{eqnarray*} $$

(98)

In both GRBs and AGNs, the usual case is that the jet head is Newtonian (see for e.g. 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₀ and a are known, but b is unconstrained, we can only learn about the quantities L_jt_j and |$t_{\rm j}\theta _{\rm op}^{-2}$|⁠, which are independent of b. From equation (95) we have L_jt_j = E₀ (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)

$$\begin{eqnarray*} t_{\rm j}\theta _{\rm op}^{-2} \sim a^2 (t_0 b^{-2}) = \left(\dfrac{\rho (z_{\rm c}) z_{\rm c}^5}{E_0}\right)^{1/2} \left(\dfrac{a}{z_{\rm c}}\right)^{(5-\alpha)/2}. \end{eqnarray*} $$

(99)

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₀ 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

$$\begin{eqnarray*} t_a \simeq (a/b)^2 t_0 \simeq \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\theta _{\rm op}^{-2} t_{\rm j}, & \tilde{L} \ll 1 \\ \tilde{L}^{-1/2} \theta _{\rm op}^{-2} t_{\rm j}, & 1 \ll \tilde{L} \ll \theta _{\rm op}^{-4/3} \\ \tilde{L}^{1/4} \theta _{\rm op}^{-1} t_{\rm j}, & \theta _{\rm op}^{-4/3} \ll \tilde{L} \ll 4\theta _{\rm op}^{-4} \end{array}\right.. \end{eqnarray*} $$

(100)

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_sph ∝ (a/b)²t₀. As |$\tilde{L}$| increases, a/b decreases while t₀ 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⁵¹ erg s⁻¹, a Lorentz factor of 5, and a specific enthalpy of 100 into a nozzle of width 8 × 10⁷ cm for a duration of t_j = 3 s. An adiabatic index of 4/3 is used in all cases. The lower boundary of the simulation is placed at z_inj = 10⁹ cm, and the external density is defined by ρ(R) = 6 × 10⁴(R/z_inj)^−α g cm⁻³. The ambient pressure is set to a negligibly low value.

We compared our results with three simulations, for three values of α relevant to astrophysical jets: α = 1, as is typical for the dense core of a galaxy cluster; α = 2, as expected for a dense circumstellar wind; and α = 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 λ = 1 in equation (16).

The simulation results are shown in Figs 14, 15, and 16 respectively, for α = 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 α = 1 and α = 2 cases, there is good agreement (within 10–15 per cent) between the analytical and numerical results out to about 45 deg 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 λ, with λ ≈ 1.2 giving an improved fit.

Figure 14.

pluto simulations of a choked relativistic jet in a ρ ∝ R⁻¹ density profile. Snapshots are shown at t = 36, 116, and 220 s, when the outflow has reached heights of roughly 2, 3, and 4 × 10¹⁰ cm. The colour scale indicates log pressure in units of dyne cm⁻². The analytical solution obtained from the KA is shown for comparison as a heavy black line.

Open in new tab Download slide

As in Fig. 14, except for ρ ∝ R−2 at t = 13, 37, and 81 s.

Figure 15.

As in Fig. 14, except for ρ ∝ R⁻² at t = 13, 37, and 81 s.

Open in new tab Download slide

Figure 16.

As in Figs 14 and 15, but for ρ ∝ R^−2.5 at t = 16, 30, and 52 s.

Open in new tab Download slide

Beyond 45 deg, 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 α = 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 ∼(P/ρ)^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 α = 1 or 2, but not for α = 2.5. This makes sense considering that, prior to choking, the jet head is accelerating in the α = 2.5 case, whereas in the other cases it is decelerating (α = 1) or has a constant velocity (α = 2). Therefore, for α = 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 α > 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 ρ ∝ 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 ≫ 1, and the ambient pressure is set to a small value (10⁻⁸ 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₀ = 840, when the height is z_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_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.

Figure 17.

The shape of the cocoon for a ρ ∝ R^−2.5 density profile and an initial aspect ratio a/b = 70/3, at a time t/t₀ = 840 when the height is z_c/a = 10.4. The colour scale indicates the ratio of the pressure to the initial cocoon pressure. The black lines show the shape of the shock in the analytical model, scaled in two different ways. The solid line is scaled to the same z_c/a as the simulation, while the dashed line is scaled to the same t/t₀.

Open in new tab Download slide

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 ∝ t^{2/(5 − α)} 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 ≈2.5. For example, whereas the simulations reach a height of z_c/a = 10.4 at a time t/t₀ = 840, the analytical model predicts a height of only z_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. equation 67 and Fig. 10) are rather accurate, while our expressions relating the height and width to time (i.e. equations 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 ρ ∝ R^−α, 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(θ_i) describing the initial shock shape, we obtain exact parametric solutions for the shock coordinates R(θ_i, z_c) and θ(θ_i, z_c) as functions of z_c, the shock height along the axis. The solutions vary in form depending on α. 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 ≫ 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₀, the planar regime persists until the age is ∼2t₀. During this time, dynamical quantities remain roughly constant, and the outflow more or less keeps its original shape. Once t ∼ 2t_b ∼ 2t₀, 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 ∼ t_a ∼ (a/b)²t₀, 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_c) and width (r_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_c/z_c going to zero as |$\propto z_{\rm c}^{(\alpha -2)/2}$| for α < 2, or as ∝ (ln z_c)⁻¹ for α = 2. An approximation for the time it takes for the outflow to become spherical in different density profiles is given in equation (79). We find that in a constant-density medium, the flow spherizes after about ∼10t_a, while in a wind-like medium it takes ∼1000t_a.

In steep density profiles with α > 2, the flow does not become spherical (at least in the idealized case where the KA applies). Instead, the ratio r_c/z_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 α → 2. For 2 < α < 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_c/a ≲ 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 α = 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 α = 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₀. For z_c/a ≫ 1, we find that the size of the simulated outflow is larger than predicted analytically by a factor of ≈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₀) and the parameters characterizing the jet (luminosity L_j, duration t_j, and opening angle θ_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 original jet’s 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₀ 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. 2019).

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_c = R_*. If the choking does not occur too deep inside the star, the width of the outflow upon breakout, r_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⁵¹ 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. Izzo et al. 2019; Piran et al. 2019). Early-time spectroscopy makes it possible to estimate the total mass (M_ej) and energy (E_ej) contained in this high-velocity component, and can even constrain the distribution of mass with velocity. E_ej corresponds to the cocoon energy E₀ in our model, while M_ej serves as a proxy for the width at breakout, with M_ej ∼ M_*(r_c/R_*)². 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, Levin & Ro 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 Council for Higher Education-Israel Science Foundation (CHE-ISF) I-Core Center for Excellence in Astrophysics. TP is supported by an advanced European Research Council (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 Israel Science Foundation (ISF) grant 1114/17.

Footnotes

1

Hereafter, to reduce clutter, we refer to R_i(θ_i) and χ_i(θ_i) as simply R_i and θ_i. These and all other quantities subscripted with ‘i’ should be understood as having an implicit dependence on θ_i.

2

Equation (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 θ is chosen as the independent variable, i.e. when the shock is described by a curve R = R(θ, 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 θ.

3

From here on, we write θ_b0(x) as simply θ_b0.

4

Interestingly, the asymptotic shape often takes a familiar mathematical form when k_α is a rational fraction. For example, when k_α = −1/2, the asymptotic shape is a cardioid.

5

We 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 α < 3/γ_a.

6

The 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{2 a^2}{b^2}$| we do not discuss further, because if a ≫ b, it requires an extreme value of α.

7

If |$\theta _{\rm b}\gt \pi /2$| for some θ_b0, the physical interpretation is that trajectories stemming from θ_b0 never become parallel to the equator.

REFERENCES

Abbott

B. P.

et al. .,

2017

,

Phys. Rev. Lett.

,

119

,

161101

10.1103/PhysRevLett.119.161101

Crossref

Search ADS

PubMed

Aloy

M. A.

,

Müller

E.

,

Ibáñez

J. M.

,

Martí

J. M.

,

MacFadyen

A.

,

2000

,

ApJ

,

531

,

L119

10.1086/312537

Crossref

Search ADS

Ayal

S.

,

Piran

T.

,

2001

,

ApJ

,

555

,

23

10.1086/321431

Crossref

Search ADS

Bannikova

E. Y.

,

Karnaushenko

A. V.

,

Kontorovich

V. M.

,

Shulga

V. M.

,

2012

,

Astron. Rep.

,

56

,

496

10.1134/S1063772912060017

Crossref

Search ADS

Begelman

M. C.

,

Cioffi

D. F.

,

1989

,

ApJ

,

345

,

L21

10.1086/185542

Crossref

Search ADS

Bîrzan

L.

,

Rafferty

D. A.

,

McNamara

B. R.

,

Wise

M. W.

,

Nulsen

P. E. J.

,

2004

,

ApJ

,

607

,

800

10.1086/383519

Crossref

Search ADS

Bisnovatyi-Kogan

G. S.

,

Silich

S. A.

,

1995

,

Rev. Mod. Phys.

,

67

,

661

10.1103/RevModPhys.67.661

Crossref

Search ADS

Bisnovatyj-Kogan

G. S.

,

Blinnikov

S. I.

,

1982

,

Astron. Zh.

,

59

,

876

Blandford

R.

,

Meier

D.

,

Readhead

A.

,

2019

,

ARA&A

,

57

,

467

Crossref

Search ADS

Bromberg

O.

,

Nakar

E.

,

Piran

T.

,

Sari

R.

,

2011

,

ApJ

,

740

,

100

10.1088/0004-637X/740/2/100

Crossref

Search ADS

Bromberg

O.

,

Nakar

E.

,

Piran

T.

,

Sari

R.

,

2012

,

ApJ

,

749

,

110

10.1088/0004-637X/749/2/110

Crossref

Search ADS

Campana

S.

et al. .,

2006

,

Nature

,

442

,

1008

Churazov

E.

,

Forman

W.

,

Jones

C.

,

Böhringer

H.

,

2000

,

A&A

,

356

,

788

Churazov

E.

,

Brüggen

M.

,

Kaiser

C. R.

,

Böhringer

H.

,

Forman

W.

,

2001

,

ApJ

,

554

,

261

10.1086/321357

Crossref

Search ADS

Diehl

S.

,

Li

H.

,

Fryer

C. L.

,

Rafferty

D.

,

2008

,

ApJ

,

687

,

173

10.1086/591310

Crossref

Search ADS

Fabian

A. C.

,

Celotti

A.

,

Blundell

K. M.

,

Kassim

N. E.

,

Perley

R. A.

,

2002

,

MNRAS

,

331

,

369

10.1046/j.1365-8711.2002.05182.x

Crossref

Search ADS

Gottlieb

O.

,

Nakar

E.

,

Piran

T.

,

2018

,

MNRAS

,

473

,

576

10.1093/mnras/stx2357

Crossref

Search ADS

Harrison

R.

,

Gottlieb

O.

,

Nakar

E.

,

2018

,

MNRAS

,

477

,

2128

10.1093/mnras/sty760

Crossref

Search ADS

Heger

A.

,

Woosley

S. E.

,

Spruit

H. C.

,

2005

,

ApJ

,

626

,

350

10.1086/429868

Crossref

Search ADS

Irwin

C. M.

,

Chevalier

R. A.

,

2016

,

MNRAS

,

460

,

1680

10.1093/mnras/stw1058

Crossref

Search ADS

Irwin

C. M.

,

Tang

X.

,

Piran

T.

,

Nakar

E.

,

2019

,

MNRAS

,

488

,

4926

10.1093/mnras/stz2075

Crossref

Search ADS

Ito

H.

,

Matsumoto

J.

,

Nagataki

S.

,

Warren

D. C.

,

Barkov

M. V.

,

2015

,

ApJ

,

814

,

L29

10.1088/2041-8205/814/2/L29

Crossref

Search ADS

Izzo

L.

et al. .,

2019

,

Nature

,

565

,

324

10.1038/s41586-018-0826-3

Crossref

Search ADS

PubMed

Kompaneets

A.

,

1960

,

Sov. Phys. Doklady

,

130

,

46

Koo

B.-C.

,

McKee

C. F.

,

1990

,

ApJ

,

354

,

513

10.1086/168712

Crossref

Search ADS

Korycansky

D. G.

,

1992

,

ApJ

,

398

,

184

10.1086/171847

Crossref

Search ADS

Kulkarni

S. R.

et al. .,

1998

,

Nature

,

395

,

663

10.1038/27139

Crossref

Search ADS

Lazzati

D.

,

Begelman

M. C.

,

2005

,

ApJ

,

629

,

903

10.1086/430877

Crossref

Search ADS

Lazzati

D.

,

Morsony

B. J.

,

Begelman

M. C.

,

2009

,

ApJ

,

700

,

L47

10.1088/0004-637X/700/1/L47

Crossref

Search ADS

Lazzati

D.

,

Morsony

B. J.

,

Blackwell

C. H.

,

Begelman

M. C.

,

2012

,

ApJ

,

750

,

68

10.1088/0004-637X/750/1/68

Crossref

Search ADS

López-Cámara

D.

,

Morsony

B. J.

,

Begelman

M. C.

,

Lazzati

D.

,

2013

,

ApJ

,

767

,

19

10.1088/0004-637X/767/1/19

Crossref

Search ADS

Lyutikov

M.

,

2011

,

MNRAS

,

411

,

2054

10.1111/j.1365-2966.2010.17836.x

Crossref

Search ADS

MacFadyen

A. I.

,

Woosley

S. E.

,

Heger

A.

,

2001

,

ApJ

,

550

,

410

10.1086/319698

Crossref

Search ADS

Martí

J. M.

,

Müller

E.

,

Font

J. A.

,

Ibáñez

J. M. Z.

,

Marquina

A.

,

1997

,

ApJ

,

479

,

151

10.1086/303842

Crossref

Search ADS

Matzner

C. D.

,

2003

,

MNRAS

,

345

,

575

10.1046/j.1365-8711.2003.06969.x

Crossref

Search ADS

Matzner

C. D.

,

Levin

Y.

,

Ro

S.

,

2013

,

ApJ

,

779

,

60

10.1088/0004-637X/779/1/60

Crossref

Search ADS

McNamara

B. R.

et al. .,

2000

,

ApJ

,

534

,

L135

10.1086/312662

Crossref

Search ADS

McNamara

B. R.

,

Nulsen

P. E. J.

,

Wise

M. W.

,

Rafferty

D. A.

,

Carilli

C.

,

Sarazin

C. L.

,

Blanton

E. L.

,

2005

,

Nature

,

433

,

45

Mészáros

P.

,

Waxman

E.

,

2001

,

Phys. Rev. Lett.

,

87

,

171102

10.1103/PhysRevLett.87.171102

Crossref

Search ADS

PubMed

Mignone

A.

,

Bodo

G.

,

Massaglia

S.

,

Matsakos

T.

,

Tesileanu

O.

,

Zanni

C.

,

Ferrari

A.

,

2007

,

ApJS

,

170

,

228

10.1086/513316

Crossref

Search ADS

Mizuta

A.

,

Aloy

M. A.

,

2009

,

ApJ

,

699

,

1261

10.1088/0004-637X/699/2/1261

Crossref

Search ADS

Mizuta

A.

,

Ioka

K.

,

2013

,

ApJ

,

777

,

162

10.1088/0004-637X/777/2/162

Crossref

Search ADS

Mizuta

A.

,

Yamasaki

T.

,

Nagataki

S.

,

Mineshige

S.

,

2006

,

ApJ

,

651

,

960

10.1086/507861

Crossref

Search ADS

Moharana

R.

,

Piran

T.

,

2017

,

MNRAS

,

472

,

L55

10.1093/mnrasl/slx131

Crossref

Search ADS

Morsony

B. J.

,

Lazzati

D.

,

Begelman

M. C.

,

2007

,

ApJ

,

665

,

569

10.1086/519483

Crossref

Search ADS

Nagakura

H.

,

Ito

H.

,

Kiuchi

K.

,

Yamada

S.

,

2011

,

ApJ

,

731

,

80

10.1088/0004-637X/731/2/80

Crossref

Search ADS

Nakar

E.

,

2015

,

ApJ

,

807

,

172

10.1088/0004-637X/807/2/172

Crossref

Search ADS

Nakar

E.

,

Piran

T.

,

2017

,

ApJ

,

834

,

28

10.3847/1538-4357/834/1/28

Crossref

Search ADS

Nakar

E.

,

Sari

R.

,

2012

,

ApJ

,

747

,

88

10.1088/0004-637X/747/2/88

Crossref

Search ADS

Piran

T.

,

Nakar

E.

,

Mazzali

P.

,

Pian

E.

,

2019

,

ApJ

,

871

,

L25

10.3847/2041-8213/aaffce

Crossref

Search ADS

Rees

M. J.

,

1984

,

ARA&A

,

22

,

471

10.1146/annurev.aa.22.090184.002351

Crossref

Search ADS

Rimoldi

A.

,

Rossi

E. M.

,

Piran

T.

,

Portegies Zwart

S.

,

2015

,

MNRAS

,

447

,

3096

10.1093/mnras/stu2630

Crossref

Search ADS

Sedov

L. I.

,

1959

,

Similarity and Dimensional Methods in Mechanics

.

Academic Press

,

New York

Soderberg

A. M.

et al. .,

2006

,

Nature

,

442

,

1014

Tang

X.

,

Churazov

E.

,

2017

,

MNRAS

,

468

,

3516

10.1093/mnras/stx590

Crossref

Search ADS

Taylor

G.

,

1950

,

Proc. R. Soc. A

,

201

,

159

10.1098/rspa.1950.0049

Crossref

Search ADS

Wang

P.

,

Abel

T.

,

Zhang

W.

,

2008

,

ApJS

,

176

,

467

10.1086/529434

Crossref

Search ADS

Waxman

E.

,

Shvarts

D.

,

1993

,

Phys. Fluids A.

,

5

,

1035

10.1063/1.858668

Crossref

Search ADS

Woosley

S. E.

,

Bloom

J. S.

,

2006

,

ARA&A

,

44

,

507

10.1146/annurev.astro.43.072103.150558

Crossref

Search ADS

Zhang

W.

,

Woosley

S. E.

,

Heger

A.

,

2004

,

ApJ

,

608

,

365

10.1086/386300

Crossref

Search ADS

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, θ_i, at the moment the jet was choked (see Section 3), we use θ_i to parametrize the relative position on the shock front. In Table A1, we list the symbols defined throughout the paper, explicitly including any dependence on t or θ_i in parentheses following each symbol.

Table A1.

Open in new tab

Glossary of symbols.

Symbol

Meaning

Introduced

Coordinates

R, θ

Polar coordinates

r, z

Cylindrical coordinates

Initial conditions

E₀

Energy of the cocoon

Section 2

a

Initial height of the cocoon

Section 2

b

Initial width of the cocoon

Section 2

P₀

Initial pressure in the cocoon

Section 2

V₀

Initial volume of the cocoon

Section 2

ρ₀

Ambient density at R = a

Equation (11)

α

Power-law index of the density profile

Equation (11)

k_α

Frequently appearing quantity defined by k_α ≡ 1 − α/2

Equation (29)

R_i(θ_i)

Curve describing the shock’s initial shape

Equation (13)

χ_i(θ_i)

Angle between the axis and the shock’s initial direction of motion

Equations (12) and (14)

ω_i(θ_i)

Angle between the shock’s initial direction of motion and the radial direction

Equation (22)

Time-scales

t₀

Time when the jet is choked

Section 2

t_b

Time it takes for the cocoon’s width to double

Equation (1)

t_a

Time it takes for the cocoon’s height to double

Equation (2)

t_sph

Time for the cocoon’s width to become 90 per cent of its height

Equations (78) and (79)

t_eq

Time to reach pressure equilibrium with the ambient medium

Equations (81) and (82)

t_w

Time when the solution breaks down in an expanding medium

Equations (83) and (84)

t_bo(R_*)

Time when the cocoon breaks out of an ambient medium extending to R_*

Equations (85) and (86)

Cocoon properties

V_c(t)

Volume of the cocoon

Section 2

P_c(t)

Average pressure in the cocoon

Section 2

λ

Dimensionless ratio of post-shock pressure to volume-averaged pressure

Equation (16)

x(t)

Dimensionless time parameter used to solve for the shock’s shape

Equation (18)

z_c(t)

Height of the cocoon

Equation (25)

r_c(t)

Width of the cocoon

Equation (26)

r_eq(t)

Width of the cocoon at the equator

See Fig. 5

ϵ_eff(t)

Effective eccentricity, i.e. |$\epsilon _{\rm eff}\equiv \sqrt{1-r_{\rm c}^2/z_{\rm c}^2}$|

Sections 3.1 and 3.2

R(θ_i, t), θ(θ_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(t), θ_b(t)

Polar coordinates of the bulge

Equation (36) and (37)

r_b(t), z_b(t)

Cylindrical coordinates of the bulge

See Fig. 5

θ_b0(t), R_b0(t), χ_b0(t)

Initial values of θ_i, R_i, and χ_i for a trajectory passing through (R_b, θ_b)

See Fig. 6

Asymptotic quantities

ζ(t)

Quantity important in the expansion of r_c(z_c) at infinity

Equation (43)

A_α

Coefficient appearing in the expansion of r_c(z_c) at infinity

Equations (44) and (63)

f_α, g_α

Ratios of r_c/z_c and r_eq/z_c, respectively, as t → ∞

Equation (39)

C_α

Ratio of V_c to |$4\pi z_{\rm c}^3/3$| as t → ∞

Equations (70) and (71)

θ_b∞

Value of θ_b(t) as t → ∞

Equation (40)

θ_{b0, min}

Value of θ_b0(t) as t → ∞

Section 3.2

Jet properties

L_j

Jet luminosity

Section 4

t_j

Jet duration

Section 4

θ_op

Jet opening angle at injection

Section 4

β_h, Γ_h

Velocity and Lorentz factor of the jet head upon choking

Section 4

β_c

Sideways velocity of the cocoon upon choking

Section 4

|$\tilde{L}$|

Dimensionless parameter governing jet propagation

Section 4

Symbol	Meaning	Introduced
Coordinates
R, θ	Polar coordinates
r, z	Cylindrical coordinates
Initial conditions
E₀	Energy of the cocoon	Section 2
a	Initial height of the cocoon	Section 2
b	Initial width of the cocoon	Section 2
P₀	Initial pressure in the cocoon	Section 2
V₀	Initial volume of the cocoon	Section 2
ρ₀	Ambient density at R = a	Equation (11)
α	Power-law index of the density profile	Equation (11)
k_α	Frequently appearing quantity defined by k_α ≡ 1 − α/2	Equation (29)
R_i(θ_i)	Curve describing the shock’s initial shape	Equation (13)
χ_i(θ_i)	Angle between the axis and the shock’s initial direction of motion	Equations (12) and (14)
ω_i(θ_i)	Angle between the shock’s initial direction of motion and the radial direction	Equation (22)
Time-scales
t₀	Time when the jet is choked	Section 2
t_b	Time it takes for the cocoon’s width to double	Equation (1)
t_a	Time it takes for the cocoon’s height to double	Equation (2)
t_sph	Time for the cocoon’s width to become 90 per cent of its height	Equations (78) and (79)
t_eq	Time to reach pressure equilibrium with the ambient medium	Equations (81) and (82)
t_w	Time when the solution breaks down in an expanding medium	Equations (83) and (84)
t_bo(R_*)	Time when the cocoon breaks out of an ambient medium extending to R_*	Equations (85) and (86)
Cocoon properties
V_c(t)	Volume of the cocoon	Section 2
P_c(t)	Average pressure in the cocoon	Section 2
λ	Dimensionless ratio of post-shock pressure to volume-averaged pressure	Equation (16)
x(t)	Dimensionless time parameter used to solve for the shock’s shape	Equation (18)
z_c(t)	Height of the cocoon	Equation (25)
r_c(t)	Width of the cocoon	Equation (26)
r_eq(t)	Width of the cocoon at the equator	See Fig. 5
ϵ_eff(t)	Effective eccentricity, i.e. \|$\epsilon _{\rm eff}\equiv \sqrt{1-r_{\rm c}^2/z_{\rm c}^2}$\|	Sections 3.1 and 3.2
R(θ_i, t), θ(θ_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(t), θ_b(t)	Polar coordinates of the bulge	Equation (36) and (37)
r_b(t), z_b(t)	Cylindrical coordinates of the bulge	See Fig. 5
θ_b0(t), R_b0(t), χ_b0(t)	Initial values of θ_i, R_i, and χ_i for a trajectory passing through (R_b, θ_b)	See Fig. 6
Asymptotic quantities
ζ(t)	Quantity important in the expansion of r_c(z_c) at infinity	Equation (43)
A_α	Coefficient appearing in the expansion of r_c(z_c) at infinity	Equations (44) and (63)
f_α, g_α	Ratios of r_c/z_c and r_eq/z_c, respectively, as t → ∞	Equation (39)
C_α	Ratio of V_c to \|$4\pi z_{\rm c}^3/3$\| as t → ∞	Equations (70) and (71)
θ_b∞	Value of θ_b(t) as t → ∞	Equation (40)
θ_{b0, min}	Value of θ_b0(t) as t → ∞	Section 3.2
Jet properties
L_j	Jet luminosity	Section 4
t_j	Jet duration	Section 4
θ_op	Jet opening angle at injection	Section 4
β_h, Γ_h	Velocity and Lorentz factor of the jet head upon choking	Section 4
β_c	Sideways velocity of the cocoon upon choking	Section 4
\|$\tilde{L}$\|	Dimensionless parameter governing jet propagation	Section 4

Table A1.

Open in new tab

Glossary of symbols.

Symbol

Meaning

Introduced

Coordinates

R, θ

Polar coordinates

r, z

Cylindrical coordinates

Initial conditions

E₀

Energy of the cocoon

Section 2

a

Initial height of the cocoon

Section 2

b

Initial width of the cocoon

Section 2

P₀

Initial pressure in the cocoon

Section 2

V₀

Initial volume of the cocoon

Section 2

ρ₀

Ambient density at R = a

Equation (11)

α

Power-law index of the density profile

Equation (11)

k_α

Frequently appearing quantity defined by k_α ≡ 1 − α/2

Equation (29)

R_i(θ_i)

Curve describing the shock’s initial shape

Equation (13)

χ_i(θ_i)

Angle between the axis and the shock’s initial direction of motion

Equations (12) and (14)

ω_i(θ_i)

Angle between the shock’s initial direction of motion and the radial direction

Equation (22)

Time-scales

t₀

Time when the jet is choked

Section 2

t_b

Time it takes for the cocoon’s width to double

Equation (1)

t_a

Time it takes for the cocoon’s height to double

Equation (2)

t_sph

Time for the cocoon’s width to become 90 per cent of its height

Equations (78) and (79)

t_eq

Time to reach pressure equilibrium with the ambient medium

Equations (81) and (82)

t_w

Time when the solution breaks down in an expanding medium

Equations (83) and (84)

t_bo(R_*)

Time when the cocoon breaks out of an ambient medium extending to R_*

Equations (85) and (86)

Cocoon properties

V_c(t)

Volume of the cocoon

Section 2

P_c(t)

Average pressure in the cocoon

Section 2

λ

Dimensionless ratio of post-shock pressure to volume-averaged pressure

Equation (16)

x(t)

Dimensionless time parameter used to solve for the shock’s shape

Equation (18)

z_c(t)

Height of the cocoon

Equation (25)

r_c(t)

Width of the cocoon

Equation (26)

r_eq(t)

Width of the cocoon at the equator

See Fig. 5

ϵ_eff(t)

Effective eccentricity, i.e. |$\epsilon _{\rm eff}\equiv \sqrt{1-r_{\rm c}^2/z_{\rm c}^2}$|

Sections 3.1 and 3.2

R(θ_i, t), θ(θ_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(t), θ_b(t)

Polar coordinates of the bulge

Equation (36) and (37)

r_b(t), z_b(t)

Cylindrical coordinates of the bulge

See Fig. 5

θ_b0(t), R_b0(t), χ_b0(t)

Initial values of θ_i, R_i, and χ_i for a trajectory passing through (R_b, θ_b)

See Fig. 6

Asymptotic quantities

ζ(t)

Quantity important in the expansion of r_c(z_c) at infinity

Equation (43)

A_α

Coefficient appearing in the expansion of r_c(z_c) at infinity

Equations (44) and (63)

f_α, g_α

Ratios of r_c/z_c and r_eq/z_c, respectively, as t → ∞

Equation (39)

C_α

Ratio of V_c to |$4\pi z_{\rm c}^3/3$| as t → ∞

Equations (70) and (71)

θ_b∞

Value of θ_b(t) as t → ∞

Equation (40)

θ_{b0, min}

Value of θ_b0(t) as t → ∞

Section 3.2

Jet properties

L_j

Jet luminosity

Section 4

t_j

Jet duration

Section 4

θ_op

Jet opening angle at injection

Section 4

β_h, Γ_h

Velocity and Lorentz factor of the jet head upon choking

Section 4

β_c

Sideways velocity of the cocoon upon choking

Section 4

|$\tilde{L}$|

Dimensionless parameter governing jet propagation

Section 4

Symbol	Meaning	Introduced
Coordinates
R, θ	Polar coordinates
r, z	Cylindrical coordinates
Initial conditions
E₀	Energy of the cocoon	Section 2
a	Initial height of the cocoon	Section 2
b	Initial width of the cocoon	Section 2
P₀	Initial pressure in the cocoon	Section 2
V₀	Initial volume of the cocoon	Section 2
ρ₀	Ambient density at R = a	Equation (11)
α	Power-law index of the density profile	Equation (11)
k_α	Frequently appearing quantity defined by k_α ≡ 1 − α/2	Equation (29)
R_i(θ_i)	Curve describing the shock’s initial shape	Equation (13)
χ_i(θ_i)	Angle between the axis and the shock’s initial direction of motion	Equations (12) and (14)
ω_i(θ_i)	Angle between the shock’s initial direction of motion and the radial direction	Equation (22)
Time-scales
t₀	Time when the jet is choked	Section 2
t_b	Time it takes for the cocoon’s width to double	Equation (1)
t_a	Time it takes for the cocoon’s height to double	Equation (2)
t_sph	Time for the cocoon’s width to become 90 per cent of its height	Equations (78) and (79)
t_eq	Time to reach pressure equilibrium with the ambient medium	Equations (81) and (82)
t_w	Time when the solution breaks down in an expanding medium	Equations (83) and (84)
t_bo(R_*)	Time when the cocoon breaks out of an ambient medium extending to R_*	Equations (85) and (86)
Cocoon properties
V_c(t)	Volume of the cocoon	Section 2
P_c(t)	Average pressure in the cocoon	Section 2
λ	Dimensionless ratio of post-shock pressure to volume-averaged pressure	Equation (16)
x(t)	Dimensionless time parameter used to solve for the shock’s shape	Equation (18)
z_c(t)	Height of the cocoon	Equation (25)
r_c(t)	Width of the cocoon	Equation (26)
r_eq(t)	Width of the cocoon at the equator	See Fig. 5
ϵ_eff(t)	Effective eccentricity, i.e. \|$\epsilon _{\rm eff}\equiv \sqrt{1-r_{\rm c}^2/z_{\rm c}^2}$\|	Sections 3.1 and 3.2
R(θ_i, t), θ(θ_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(t), θ_b(t)	Polar coordinates of the bulge	Equation (36) and (37)
r_b(t), z_b(t)	Cylindrical coordinates of the bulge	See Fig. 5
θ_b0(t), R_b0(t), χ_b0(t)	Initial values of θ_i, R_i, and χ_i for a trajectory passing through (R_b, θ_b)	See Fig. 6
Asymptotic quantities
ζ(t)	Quantity important in the expansion of r_c(z_c) at infinity	Equation (43)
A_α	Coefficient appearing in the expansion of r_c(z_c) at infinity	Equations (44) and (63)
f_α, g_α	Ratios of r_c/z_c and r_eq/z_c, respectively, as t → ∞	Equation (39)
C_α	Ratio of V_c to \|$4\pi z_{\rm c}^3/3$\| as t → ∞	Equations (70) and (71)
θ_b∞	Value of θ_b(t) as t → ∞	Equation (40)
θ_{b0, min}	Value of θ_b0(t) as t → ∞	Section 3.2
Jet properties
L_j	Jet luminosity	Section 4
t_j	Jet duration	Section 4
θ_op	Jet opening angle at injection	Section 4
β_h, Γ_h	Velocity and Lorentz factor of the jet head upon choking	Section 4
β_c	Sideways velocity of the cocoon upon choking	Section 4
\|$\tilde{L}$\|	Dimensionless parameter governing jet propagation	Section 4

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 θ_b0. As discussed in Section 3.2, a trajectory originating from θ_b0 becomes parallel to the equator at the point (R_b, θ_b), with θ_b and R_b given by equations (36) and (37), respectively. How long does it take for a particle at (R_b0, θ_b0) to travel to (R_b, θ_b)? By substituting θ(x(θ_b0), θ_b0) = θ_b into equation (34) or (52), as appropriate, we obtain

$$\begin{eqnarray*} x(\theta _{\rm b0}) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\dfrac{1}{k_\alpha } \left(\dfrac{R_{\rm b0}}{a}\right)^{k_\alpha } \dfrac{\sin (k_\alpha \Delta _{\rm b})}{\sin (\omega _{\rm b0}-k_\alpha \Delta _{\rm b})} & \alpha \ne 2 \\ \dfrac{\Delta _{\rm b}}{\sin \omega _{\rm b0}} & \alpha =2 \end{array}\right. , \end{eqnarray*} $$

(B1)

where we have defined

$$\begin{eqnarray*} \Delta _{\rm b}(\theta _{\rm b0}) \equiv \theta _{\rm b}-\theta _{\rm b0}= \dfrac{2}{\alpha }\left(\pi /2-\chi _{\rm b0}\right) \end{eqnarray*} $$

(B2)

and

$$\begin{eqnarray*} \omega _{\rm b0}(\theta _{\rm b0}) \equiv \chi _{\rm b0}-\theta _{\rm b0} \end{eqnarray*} $$

(B3)

to make the formulae more compact. The values of R_b0 and χ_b0 are given by setting θ_i = θ_b0 in equations (13) and (14), respectively.

If the initial shape is an ellipsoid with R_b0 and χ_b0 given by equations (13) and (14) at θ_i = θ_b0, then equation (B1) only has a solution when x is greater than a critical value,

$$\begin{eqnarray*} x_b = \lim _{\theta _{\rm b0}\rightarrow \pi /2} x(\theta _{\rm b0}) = \left(\dfrac{b}{a}\right)^{k_\alpha } \dfrac{2 b^2/\alpha a^2}{1-2b^2/\alpha a^2}. \end{eqnarray*} $$

(B4)

Physically, this corresponds to the time when the bulge width becomes equal to the equatorial width, i.e. we have r_b = r_eq at the time t = t(x_b). The relation between r_c and z_c is different before and after this transition, with r_c = r_eq up to x = x_b, and r_c = r_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_c = r_eq at all times. However, as long as b/a is small, this is only relevant to very flat density profiles with α ≈ 0. In the case where b ≪ a and α ≫ 2b²/a², 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 α > 2b²/a² 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_b0 = b. We then have r_c = r_eq, and

$$\begin{eqnarray*} z_{\rm c}(x) &=& a(1+k_\alpha x)^{1/k_\alpha } \nonumber\\ r_{\rm c}(x) &=& b\left[1+(a/b)^{k_\alpha } k_\alpha x\right]^{1/k_\alpha }. \end{eqnarray*} $$

(B5)

In this regime, it is possible to express r_c as an explicit function of z_c

$$\begin{eqnarray*} r_{\rm c}(z_{\rm c}) = r_{\rm eq}(z_{\rm c}) = \left\lbrace \begin{array}{ll}b \left\lbrace 1 + \left(\dfrac{a}{b}\right)^{k_\alpha } \left[\left(\dfrac{z_{\rm c}}{a}\right)^{k_\alpha }-1\right]\right\rbrace ^{1/k_\alpha } & \,\,\,\alpha \ne 2 \\ b \left(\dfrac{z_{\rm c}}{a}\right) & \,\,\,\alpha =2 \end{array}\right. . \nonumber\\ \end{eqnarray*} $$

(B6)

For k_α = 1 (constant density), we reproduce equation (27).

When α > 2b²/a² and x > x_b, however, the width is governed by the bulge, and the relation between r_c and z_c becomes more nuanced. It is no longer possible to express r_c as an explicit function of z_c, but the width and height can be parametrized in terms of θ_b0. The height of the cocoon when x = x(θ_b0) is, combining equations (38) or (54) with (B1),

$$\begin{eqnarray*} z_{\rm c}(\theta _{\rm b0}) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}a \left[1 + \left(\dfrac{R_{\rm b0}}{a}\right)^{k_\alpha } \dfrac{\sin (k_\alpha \Delta _{\rm b})}{\cos \theta _{\rm b}} \right]^{1/k_\alpha } & \alpha \ne 2 \\ a {\rm e}^{\Delta _{\rm b}\csc \omega _{\rm b0}} & \alpha =2 \end{array}\right., \end{eqnarray*} $$

(B7)

with θ_b given by equation (36). The width, on the other hand, is r_c = r_b = R_bsin θ_b; from equations (36) and (37), we have

$$\begin{eqnarray*} r_{\rm c}(\theta _{\rm b0}) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}R_{\rm b0}\sin \theta _{\rm b}\left(\dfrac{\sin \omega _{\rm b0}}{\cos \theta _{\rm b}}\right)^{1/k_\alpha } & \alpha \ne 2 \\ R_{\rm b0}\sin \theta _{\rm b}{\rm e}^{\Delta _{\rm b}\cot \omega _{\rm b0}} & \alpha = 2 \end{array}\right.. \end{eqnarray*} $$

(B8)

Similarly, the height of the bulge above the equatorial plane is z_b = R_bcos θ_b, or

$$\begin{eqnarray*} z_{\rm b}(\theta _{\rm b0}) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}R_{\rm b0}\cos \theta _{\rm b}\left(\dfrac{\sin \omega _{\rm b0}}{\cos \theta _{\rm b}}\right)^{1/k_\alpha } & \alpha \ne 2 \\ R_{\rm b0}\cos \theta _{\rm b}{\rm e}^{\Delta _{\rm b}\cot \omega _{\rm b0}} & \alpha = 2 \end{array}\right.. \end{eqnarray*} $$

(B9)

The equatorial width in this case is discussed separately, in Appendix D. In the first part of equations (B7)–(B9), 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 (B2)–(B3).

Since we expect r_c ∼ z_c to hold approximately at late times, it is sometimes more convenient to work with the dimensionless ratio r_c/z_c. Dividing equation (B8) by equation (B7), we obtain

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \!=\! \left\lbrace \!\!\!\begin{array}{ll}\left(\dfrac{R_{\rm b0}}{a}\right) \sin \theta _{\rm b}\left[ \dfrac{\sin \omega _{\rm b0}}{\cos \theta _{\rm b}+ (R_{\rm b0}/a)^{k_\alpha } \sin (k_\alpha \Delta _{\rm b})} \right]^{1/k_\alpha } & \,\,\alpha \ne 2 \\ \left(\dfrac{R_{\rm b0}}{a}\right) \sin \theta _{\rm b}{\rm e}^{\Delta _{\rm b}(\cot \omega _{\rm b0}-\csc \omega _{\rm b0})} & \,\,\alpha =2 \end{array}\right.\! . \nonumber\\ \end{eqnarray*} $$

(B10)

APPENDIX C: EVOLUTION OF THE COCOON WIDTH VERSUS HEIGHT

Approximate formulae for r_c versus z_c can be obtained in certain limiting situations where equations (B7) and (B8) simplify. To understand the interesting limiting cases, let us investigate how the three angular scales appearing in equation (B10) (θ_b, ω_b0, and Δ_b) vary with θ_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 χ_b0, ω_b0 ≡ χ_b0 − θ_b0, θ_b, and Δ_b ≡ θ_b − θ_b0 versus θ_b0 in Fig. C1 for α = 1 (left-hand panel) or 3 (right-hand panel) and a/b = 10. Focusing now on the θ_b and ω_b0 curves, the key features are as follows:

Differentiation of equation (36) implies that θ_b has a minimum if dχ_b0/dθ_b0 = α/2 is satisfied. Applying equation (14), we see that this condition is only fulfilled when |$\frac{2b^2}{a^2} \lt \alpha \lt \frac{2a^2}{b^2}$|⁠. If we assume that b ≪ a, and that |$\frac{2b^2}{a^2} \ll \alpha \ll \frac{2a^2}{b^2}$| holds as well,⁶ then θ_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, θ_b rises to |$\pi /2$| at |$\theta _{\rm b0}=\pi /2$|⁠, and to |$\pi /\alpha$| at θ_b0 = 0.⁷ If α < 2, θ_b becomes equal to |$\pi /2$| at θ_b0 = θ_{b0, min} (see below).
ω_b0 is determined from equation (14), and its behaviour is straightforward since there is no dependence on α: ω_b0 is 0 when θ_b0 = 0, rises to a maximum value of |$\omega _{\rm b0}\approx \pi /2-b/a$| at θ_b0 ≈ b/a, and then goes to zero again as |$\theta _{\rm b0}\rightarrow \pi /2$|⁠.
The θ_b and ω_b0 curves have two crossings, one at θ_b0 ∼ b²/a², and another at |$\theta _{\rm b0}\sim \pi /4$|⁠. Over the range |$b^2/a^2 \lesssim \theta _{\rm b0}\lesssim \pi /4$|⁠, θ_b < ω_b0; otherwise θ_b > ω_b0.

$The angular scales θb0, χb0, θb, ωb0, and Δb as functions of θb0, for an initial ellipsoid with b/a = 0.1 and α = 1 (left-hand panel) or α = 3 (right-hand panel). The four dynamical phases described in the text are labelled with Roman numerals and separated by vertical dashed lines. The critical value θb0, min (which is only defined for α < 2) is indicated by a dotted line. The greyed-out region has $\theta _{\rm b}\gt \pi /2$ and does not affect the evolution of rc/zc.$

Figure C1.

The angular scales θ_b0, χ_b0, θ_b, ω_b0, and Δ_b as functions of θ_b0, for an initial ellipsoid with b/a = 0.1 and α = 1 (left-hand panel) or α = 3 (right-hand panel). The four dynamical phases described in the text are labelled with Roman numerals and separated by vertical dashed lines. The critical value θ_{b0, min} (which is only defined for α < 2) is indicated by a dotted line. The greyed-out region has |$\theta _{\rm b}\gt \pi /2$| and does not affect the evolution of r_c/z_c.

Open in new tab Download slide

As discussed in Section 3.2.1, θ_{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 → ∞. Then, from equation (34), we get an implicit equation,

$$\begin{eqnarray*} \theta _{\rm b0,min}= \dfrac{2}{\alpha }\left(\chi _{\rm b0,min}-k_\alpha \pi /2\right), \end{eqnarray*} $$

(C1)

where χ_{b0, min} is related to θ_{b0, min} by equation (14). Alternatively, formula (C1) may be obtained by taking |$\theta _{\rm b}\rightarrow \pi /2$| in (36). Note that when α = 2b²/a², the solution to equation (C1) is |$\theta _{\rm b0,min}= \chi _{\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 α increases from this value, θ_{b0, min} decreases, and a larger portion of the cocoon becomes subject to collisions. When α → 2, θ_{b0, min} → 0, and all parts of the surface will interact given sufficient time. As long as α ≫ 2b²/a², θ_{b0, min} is small and can be approximated by

$$\begin{eqnarray*} \theta _{\rm b0,min}\approx (b^2/a^2) \tan \left(k_\alpha \pi /2\right). \end{eqnarray*} $$

(C2)

The outflow’s shape evolves in a different way depending on the relative size of the angular parameters. There are four distinct evolutionary phases (labelled by Roman numerals in the figure) as θ_b0 is decreased from |$\pi /2$| to θ_{b0, min} (when α < 2) or to 0 (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:

Phase I: |$\theta _{\rm b0}\approx \pi /2$| and χ_b0 > θ_b0 ≈ θ_b ≫ ω_b0 ≫ Δ_b. This corresponds to the early planar phase, when x ∼ x_b (see also Appendix B), and the bulge is still located close to the equator.
Phase II: |$b/a \ll \theta _{\rm b0}\ll \pi /4$| and χ_b0 > ω_b0 ≫ θ_b0 ≈ θ_b ≫ Δ_b. This applies during the later part of the planar phase, once x ≫ x_b and the bulge is closer to the axis than to the equator.
Phase III: b²/a² ≪ θ_b0 ≪ b/a and χ_b0 ≃ ω_b0 ≫ θ_b ≃ Δ_b ≫ θ_b0. This holds during the sideways expansion phase.
Phase IV: θ_b0 → θ_{b0, min} or θ_b0 → 0 (respectively for α < 2 and |$\alpha \geq 2$|⁠), and θ_b ≈ Δ_b > χ_b0 ≈ ω_b0 ≫ θ_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 δ, we find |$\chi _{\rm b0}\approx \pi /2-(b^2/a^2)\delta$|⁠, ω_b0 ≈ (1 − b²/a²)δ, |$\theta _{\rm b}\approx \pi /2-(1-2 b^2/\alpha a^2)\delta$|⁠, and Δ_b ≈ (2b²/αa²)δ. In this case, the term in brackets in equation (B10) is unity, up to a small correction of order b²/a². Therefore we can write r_c/z_c ≈ (R_b0/a)sin θ_b ≈ b/a.
Phase II (⁠|$b/a \lesssim \theta _{\rm b0}\lesssim \pi /4$|⁠): In this limit |$\delta =\pi /2-\theta _{\rm b0}\geq \pi /4$| is no longer a small quantity. However, (b²/a²)δ is still small. So, to a reasonable approximation θ_b ≈ θ_b0, |$\chi _{\rm b0}\approx \pi /2$|⁠, and |$\omega _{\rm b0}\approx \pi /2-\theta _{\rm b0}$|⁠. Thus the bracketed part of equation (B6) is still order-unity, and r_c/z_c ≈ (R_b0/a)sin θ_b remains valid. Although R_b and θ_b both vary with θ_b0 in this case, one decreases while the other increases. The effects roughly cancel so that r_b turns out to be nearly constant with respect to θ_b0. When phase II ends, θ_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_c/z_c ≃ 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²/a² ≲ θ_b0 ≲ b/a): In phases III and IV, θ_b0 is the smallest angular scale in the problem. In the limit θ_b0 → 0, the shock’s shape is similar to that of an on-axis point explosion at z = a. Since θ_b0 can be neglected, we have θ_b ≈ Δ_b and |$\chi _{\rm b0}\approx \omega _{\rm b0}\approx \pi /2-(\alpha /2)\theta _{\rm b}$|⁠. Equations (B6) and (B7) then simplify to functions of θ_b alone
$$\begin{eqnarray*} z_{\rm c}(\theta _{\rm b}) \approx \left\lbrace \begin{array}{@{}l@{\quad }l@{}}a \left[1+ \left(\dfrac{R_{\rm b0}}{a}\right)^{k_\alpha } \dfrac{\sin (k_\alpha \theta _{\rm b})}{\cos \theta _{\rm b}}\right]^{1/k_\alpha } & \alpha \ne 2 \\ a {\rm e}^{\theta _{\rm b}/\cos (\alpha \theta _{\rm b}/2)} & \alpha =2 \end{array}\right. \end{eqnarray*} $$
(C3)
and
$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \approx \left\lbrace \begin{array}{ll}\left(\dfrac{R_{\rm b0}}{a}\right) \sin \theta _{\rm b}\left[ \dfrac{\cos (\alpha \theta _{\rm b}/2)}{\cos \theta _{\rm b}+ (R_{\rm b0}/a)^{k_\alpha } \sin (k_\alpha \theta _{\rm b})} \right]^{1/k_\alpha } &\,\,\, \alpha \ne 2 \\ \left(\dfrac{R_{\rm b0}}{a}\right) \sin \theta _{\rm b}{\rm e}^{\theta _{\rm b}(\tan \theta _{\rm b}-\sec \theta _{\rm b})} &\,\,\, \alpha =2 \end{array}\right. ; \nonumber\\ \end{eqnarray*} $$
(C4)
these expressions are valid for both phase III and phase IV, since so far we have only assumed χ_b0 ≫ θ_b0 and θ_b ≫ θ_b0. In what follows it will also be useful to know how the parameter ζ varies with θ_b. Combining equation (C3) with the definition of ζ from equation (43), we obtain
$$\begin{eqnarray*} \zeta \approx \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\left(\dfrac{R_{\rm b0}}{a}\right)^{k_\alpha } \dfrac{\sin (k_\alpha \theta _{\rm b})}{k_\alpha \cos \theta _{\rm b}} & \alpha \lt 2 \\ \dfrac{\theta _{\rm b}}{\cos \theta _{\rm b}} & \alpha =2 \\ \dfrac{\sin (k_\alpha \theta _{\rm b})}{k_\alpha [(a/R_{\rm b0})^{k_\alpha } \cos \theta _{\rm b}+\sin (k_\alpha \theta _{\rm b}) ]} & \alpha \gt 2 \end{array}\right.. \end{eqnarray*} $$
(C5)
For phase III specifically, we want θ_b ≪ χ_b0, and since |$\chi _{\rm b0}\leq \pi /2$| this implies θ_b must be small. With the approximation R_b0 ≈ a, expanding (C3) and (C4) to leading order for θ_b ≪ 1 yields simply r_c/z_c ≃ θ_b and z_c ≃ a(1 + θ_b). Therefore,
$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \simeq \dfrac{z_{\rm c}}{a}-1. \end{eqnarray*} $$
(C6)
We caution that this approximation is coarse and is only relevant for sufficiently large a/b. The r_c ∝ z_c − a behaviour can be understood through analogy with a point explosion at z = a. Within a small neighbourhood of size ≪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_c = a is a roughly sphere of radius z_c − a. As phase III ends when θ_b ∼ 1 and z_c ∼ a, it can be identified with the sideways expansion phase discussed in Section 2.
Phase IV (θ_b0 ≲ b²/a²): Finally, we address the case where θ_b > χ_b0 ≫ θ_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 α. Below, we separately discuss the cases α < 2, α = 2, and α > 2. The basic solution method is the same in each case. First, we create some small quantity δ ≡ θ_b∞ − θ_b, with θ_b∞ defined as in equation (40). Next, we replace θ_b by δ in equation (C4), and expand to leading order in δ. We then do the same for equation (C5) to obtain ζ in terms of δ. Finally, we combine the two expressions to determine r_c/z_c to first order in ζ.

For α < 2, we have |$\theta _{\rm b}\rightarrow \pi /2$| as θ_b0 → θ_{b0, min}, so we let |$\delta = \pi /2-\theta _{\rm b}$|⁠. Then, ignoring terms of order δ² or higher, we find cos θ_b ≈ δ; sin θ_b ≈ 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 (C5), we obtain

$$\begin{eqnarray*} \zeta ^{-1} \approx \left[(R_{\rm b0,min}/a)^{k_\alpha } \sin (k_\alpha \pi /2)\right]^{-1} k_\alpha \delta \ll 1, \end{eqnarray*} $$

(C7)

where R_{b0, min} ≡ R_i(θ_i = θ_{b0, min}). Now, expanding equation (C4) and using (C7) to replace δ by ζ, we find

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \approx 1 - \dfrac{1}{k_\alpha ^2}\left[1-\left(\dfrac{R_{\rm b0,min}}{a}\right)^{k_\alpha } \cos \left(\dfrac{k_\alpha \pi }{2}\right) \right] \zeta ^{-1}. \end{eqnarray*} $$

(C8)

If we further assume that θ_{b0, min} is small (i.e. that α ≫ 2b²/a²), then R_{b0, min} ≈ a and equation (C8) reduces to the form in equations (42)–(44).

For α = 2, we can repeat the same procedure using the second expressions in (C4) and (C5), or simply take the limit k_α → 0 in formula (C8). Either way, we recover equation (57).

Lastly, we treat the case of α > 2, for which θ_b0 → 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_b0/a factors. Using k_α = 1 − α/2, the denominator appearing in equations (C4) and (C5) can be rewritten as:

$$\begin{eqnarray*} \cos \theta _{\rm b}+\sin (k_\alpha \theta _{\rm b}) &=& \sin \theta _{\rm b}\cos (\alpha \theta _{\rm b}/2)\nonumber\\ &&+ \,\cos \theta _{\rm b}\left[1- \cos (\delta \alpha /2) \right]. \end{eqnarray*} $$

(C9)

The term in brackets can be neglected, since it vanishes as δ² when δ → 0. Therefore, equations (C4) and (C5) simplify to

$$\begin{eqnarray*} \dfrac{r_{\rm c}}{z_{\rm c}} \approx \left[\sin \theta _{\rm b}\right]^{1-1/k_\alpha } \end{eqnarray*} $$

(C10)

and

$$\begin{eqnarray*} \zeta \approx \dfrac{\sin (k_\alpha \theta _{\rm b})}{k_\alpha \sin \theta _{\rm b}\cos (\alpha \theta _{\rm b}/2)}. \end{eqnarray*} $$

(C11)

When |$\theta _{\rm b}\rightarrow \pi /\alpha$|⁠, we see that r_c/z_c approaches the value f_α given by equation (62).

Expanding once more in terms of δ 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 (αθ_b/2) ≈ (α/2)δ. So, to leading order,

$$\begin{eqnarray*} \zeta ^{-1} \approx (\alpha /2) |k_\alpha | \delta \tan (\pi /\alpha). \end{eqnarray*} $$

(C12)

Expanding formula (C10) up to first order in δ, and then substituting δ for ζ according to equation (C12), 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_eq can also be parametrized in terms of θ_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_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 θ_eq0 is their initial angle. Similar to the way that the width of the bulge can be parametrized in terms of θ_b0 (see Appendix B), the width at the equator can be parametrized in terms of θ_eq0. Equation (34) then leads to the condition

$$\begin{eqnarray*} x(\theta _{\rm eq0}) = \dfrac{1}{k_\alpha } \left(\dfrac{R_{\rm eq0}}{a}\right)^{k_\alpha } \dfrac{\sin (k_\alpha \Delta _{\rm eq})}{\sin (\omega _{\rm eq0}-k_\alpha \Delta _{\rm eq})}, \end{eqnarray*} $$

(D1)

where (in analogy with Appendix B) we defined

$$\begin{eqnarray*} \Delta _{\rm eq}\equiv \pi /2-\theta _{\rm eq0} \end{eqnarray*} $$

(D2)

and

$$\begin{eqnarray*} \omega _{\rm eq0}\equiv \chi _{\rm eq0}-\theta _{\rm eq0}, \end{eqnarray*} $$

(D3)

and R_eq0 and χ_eq0 are now given by taking θ_i = θ_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 θ_eq0 which satisfy x(t) = x(θ_eq0) on the interval |$0 \lt \theta _{\rm eq0}\lt \pi /2$|⁠. The number of solutions for θ_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 θ_eq0, the equatorial radius r_eq is determined by the ejecta which started at |$\theta _i=\pi /2$|⁠, and r_eq is described by equation (B6). Otherwise, the value of r_eq is determined by the intersection point farthest from the axis, which is set by the smallest solution for θ_eq0. We then have r_eq = R(x(θ_eq0), θ_eq0). Inserting x = x(θ_eq0) into equation (33) and simplifying, we obtain

$$\begin{eqnarray*} r_{\rm eq}(\theta _{\rm eq0}) = R_{\rm eq0}\left(\dfrac{\sin \omega _{\rm eq0}}{\sin (\omega _{\rm eq0}-k_\alpha \Delta _{\rm eq})} \right)^{1/k_\alpha }. \end{eqnarray*} $$

(D4)

To compare r_eq to the values of z_c, r_c, and z_b given in Appendix B, it is necessary to determine the relationship between θ_eq0 and θ_b0. This is done by equating the two expressions for x, (B1) and (D1), and solving numerically to obtain θ_eq0 in terms of θ_b0, or vice versa.

Equation (D1) 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 α 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 α 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 < α < 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 (D1) when |$\alpha \leq 2b^2/a^2$|⁠, indicating that no collisions occur in this case. When α > 2b²/a², we find that glancing collisions start to occur at x = x_eq, where

$$\begin{eqnarray*} x_{\mathrm{ eq}} \equiv \lim _{\theta _{\rm eq0}\rightarrow \pi /2} x(\theta _{\rm eq0}) = \left(\dfrac{b}{a}\right)^{k_\alpha } \dfrac{2/\alpha }{1-2b^2/\alpha a^2}=(a^2/b^2)x_b \end{eqnarray*} $$

(D5)

for an initially ellipsoidal cocoon. Just as r_c has a transition at x = x_b (as discussed in Appendix B), r_eq breaks to a different behaviour at x = x_eq. While x < x_eq, there is no solution for θ_eq0 and r_eq is given by equation (B6). Once x > x_eq, ejecta from |$\theta _i \lt \pi /2$| start to arrive at the equator, and there is one solution for θ_eq0. Thereafter, r_eq is determined by equation (D4). Head-on collisions never occur for |$\alpha \leq 2$| (see below).

When α > 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θ = 0 at the point |$(R,\theta)=(r_{\rm eq},\pi /2)$|⁠. From equation (35), we find that ejecta originating from θ_i undergo a head-on collision if the condition

$$\begin{eqnarray*} \theta _i = \dfrac{2}{\alpha }\left[\chi _i-(1+k_\alpha)\pi /2\right], \end{eqnarray*} $$

(D6)

is satisfied. Solving for α yields

$$\begin{eqnarray*} \alpha = \dfrac{2(\pi -\chi _i)}{\pi /2-\theta _i}. \end{eqnarray*} $$

(D7)

For an initially ellipsoidal shape, the right-hand side of equation (D7) has a minimum value. Therefore, equation (D7) only has a solution if |$\alpha \geq \alpha _{\mathrm{ col},\mathrm{ ho}}$|⁠, where

$$\begin{eqnarray*} \alpha _{\mathrm{ col},\mathrm{ ho}} = \min \left(\dfrac{2(\pi -\chi _i)}{\pi /2-\theta _{i}}\right) \approx 2\left(1+\dfrac{4b}{\pi a}\right), \end{eqnarray*} $$

(D8)

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

In the |$\alpha \geq \alpha _{\mathrm{ col},\mathrm{ ho}}$| case, equation (D6) admits two solutions. The relevant one for determining r_eq is approximated by

$$\begin{eqnarray*} \theta _i \approx (b^2/a^2)\tan [\pi (1-\alpha /4)] \end{eqnarray*} $$

(D9)

for b/a ≪ 1. Equation (D9) 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 θ_i = θ_eq0, equation (D6) implies |$\omega _{\rm eq0}\approx \pi (1-\alpha /4)$|⁠, |$\Delta _{\rm eq}\approx \pi /2$|⁠, and R_eq0 ≈ a. Inserting these values into equations (D1) and (D4), we find that the first head-on collision occurs at a time

$$\begin{eqnarray*} x_{\mathrm{ col}} \approx \dfrac{1}{k_\alpha } \sin (k_\alpha \pi /2) \end{eqnarray*} $$

(D10)

and a radius

$$\begin{eqnarray*} r_{\mathrm{ col}} \approx a \left(\sin [(1+k_\alpha)\pi /2]\right)^{1/k_\alpha }. \end{eqnarray*} $$

(D11)

The evolution of r_eq in the 2 < α < 3 case depends on how α compares with α_{col, ho}. If 2 < α < α_{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 \lt \alpha \leq 2$| discussed above. In this case, the equatorial radius r_eq grows continuously with x. However, if |$\alpha _{\mathrm{ col},\mathrm{ ho}} \leq \alpha \lt 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 ≈r_col, so that r_eq has a discontinuity at x_col. The value of r_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. D1, where we plot the critical values of α as functions of b/a. The phase space of α versus b/a can be divided into the four regions labelled A through D in the figure, each with a distinct behaviour:

Region A (⁠|$\alpha \leq 2b^2/a^2$|⁠): There are no collisions.
Region B (⁠|$2b^2/a^2 \lt \alpha \leq 2$|⁠): The part of the surface with θ_i > θ_{b0, min} eventually experiences glancing collisions. The rest of the surface never reaches the equator.
Region C (2 < α < α_{col, ho}): All the ejecta ultimately undergo a glancing collision.
Region D (⁠|$\alpha _{\mathrm{ col},\mathrm{ ho}} \leq \alpha \lt 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.

The critical values of α = 2b2/a2 and α = αcol, ho versus b/a for an initially ellipsoidal cocoon. The dashed line shows the approximation for αcol, ho for b/a ≪ 1, as given by equation (D8). The regions labelled A through D are discussed in the text.

Figure D1.

The critical values of α = 2b²/a² and α = α_{col, ho} versus b/a for an initially ellipsoidal cocoon. The dashed line shows the approximation for α_{col, ho} for b/a ≪ 1, as given by equation (D8). The regions labelled A through D are discussed in the text.

Open in new tab Download slide

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
August 2019	2
September 2019	21
October 2019	23
November 2019	10
December 2019	5
January 2020	11
February 2020	7
March 2020	6
April 2020	3
May 2020	1
June 2020	7
July 2020	5
August 2020	1
September 2020	3
October 2020	1
November 2020	13
December 2020	1
January 2021	2
February 2021	7
March 2021	1
April 2021	2
May 2021	1
June 2021	5
July 2021	5
September 2021	5
October 2021	3
November 2021	10
December 2021	3
April 2022	2
May 2022	3
June 2022	5
August 2022	9
September 2022	4
October 2022	5
November 2022	5
December 2022	5
January 2023	6
February 2023	10
March 2023	9
April 2023	9
June 2023	7
July 2023	3
August 2023	6
September 2023	5
October 2023	4
November 2023	6
December 2023	7
January 2024	9
February 2024	9
March 2024	6
April 2024	6
May 2024	5
June 2024	7
July 2024	25
August 2024	12
September 2024	8
October 2024	9
November 2024	8
December 2024	8
January 2025	15
February 2025	8
March 2025	16
April 2025	7
May 2025	4

Article Contents

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

ABSTRACT

1 INTRODUCTION

2 OVERVIEW

3 ANALYTICAL SOLUTIONS FOR COCOON EVOLUTION

3.1 A constant external density

3.2 A power-law external density

3.2.1 A shallow external density profile (0 < α < 2)

3.2.2 A wind-like external density profile (α = 2)

3.2.3 A steep external density profile (2 < α < 3)

3.3 The volume of the cocoon

3.4 The age of the cocoon

3.5 Eventual breakdown of the solution

4 INFERRING THE JET PARAMETERS FROM THE COCOON PROPERTIES

5 COMPARISON WITH NUMERICAL SIMULATIONS

6 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: GLOSSARY OF SYMBOLS

APPENDIX B: PARAMETRIC EQUATIONS FOR THE COCOON SHAPE

APPENDIX C: EVOLUTION OF THE COCOON WIDTH VERSUS HEIGHT

APPENDIX D: COLLISIONS AND THE EVOLUTION OF THE EQUATORIAL WIDTH VERSUS HEIGHT

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

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

ABSTRACT

1 INTRODUCTION

2 OVERVIEW

3 ANALYTICAL SOLUTIONS FOR COCOON EVOLUTION

3.1 A constant external density

3.2 A power-law external density

3.2.1 A shallow external density profile (0 < α < 2)

3.2.2 A wind-like external density profile (α = 2)

3.2.3 A steep external density profile (2 < α < 3)

3.3 The volume of the cocoon

3.4 The age of the cocoon

3.5 Eventual breakdown of the solution

4 INFERRING THE JET PARAMETERS FROM THE COCOON PROPERTIES

5 COMPARISON WITH NUMERICAL SIMULATIONS

6 DISCUSSION AND CONCLUSIONS

ACKNOWLEDGEMENTS

Footnotes

REFERENCES

APPENDIX A: GLOSSARY OF SYMBOLS

APPENDIX B: PARAMETRIC EQUATIONS FOR THE COCOON SHAPE

APPENDIX C: EVOLUTION OF THE COCOON WIDTH VERSUS HEIGHT

APPENDIX D: COLLISIONS AND THE EVOLUTION OF THE EQUATORIAL WIDTH VERSUS HEIGHT

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

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