PROPAGATION OF RELATIVISTIC, HYDRODYNAMIC, INTERMITTENT JETS IN A ROTATING, COLLAPSING GRB PROGENITOR STAR

Gamma-ray bursts (GRBs) are extremely energetic explosions with enormous gamma-ray radiation. The observational properties of GRBs demand that they originate from relativistic jets beaming toward Earth (see Piran 2004; Mészáros 2006; Kumar & Zhang 2015, for reviews). Observations of GRB afterglows further reveal that these jets are narrowly beamed with a typical opening angle around 0.1 rad (Frail et al. 2001; Liang et al. 2008; Wang et al. 2015). The so-called long GRBs (with durations longer than 2 s) are widely believed to be associated with the death of massive stars (Woosley 1993; Galama et al. 1998; Hjorth et al. 2003; Stanek et al. 2003; Woosley & Bloom 2006). The possible mechanisms to launch a relativistic jet from the central engine (likely a hyperaccreting black hole [BH] or a millisecond magnetar) include magnetohydrodynamical processes (e.g., Blandford & Znajek 1977; Narayan & McClintock 2012) or neutrino annihilation processes (e.g., Popham et al. 1999; Kohri & Mineshige 2002; Gu et al. 2006; Liu et al. 2015). After launching, the jet propagates through and breaks out of the stellar envelope of the progenitor star before emitting γ-ray photons at a large radius.

This paper focuses on the GRB jet propagation process inside the progenitor star. Previous studies have revealed some main characteristics of a propagating, constant-luminosity jet. MacFadyen & Woosley (1999) first used hydrodynamical simulations to show that a light jet can penetrate through a star and remains highly beamed. Later, jet propagation was investigated by many authors using different codes that implement special relativity and improved numerical resolutions (e.g., Aloy et al. 2000; Zhang et al. 2003; Mizuta et al. 2006; Morsony et al. 2007; Mizuta & Aloy 2009; Mizuta & Ioka 2013). These works showed detailed temporal and angular properties of the jets emerging from massive stars. In the meantime, analytical studies have been carried out on the topic (Matzner 2003; Waxman & Mészáros 2003; Bromberg et al. 2011). According to these analytical and numerical approaches, the basic picture for jet propagation may be summarized in the evolution of two phases defined by the time domains before and after the jet breaks out of the surface of the star, respectively. During the first jet-penetrating phase, one may spatially define five structures: beam, envelope, jet head, cocoon, and multiple oblique collimation shocks. When the beam is injected from the base of the star, a collimation shock (Komissarov & Falle 1997, 1998; Bromberg & Levinson 2009; Mizuta & Ioka 2013) would emerge to generate the pressure needed to counterbalance the pressure from the surrounding materials (the cocoon). The jet begins to propagate by pushing aside the matter in front of it. As a consequence, a jet head consisting of a forward shock (FS) and a reverse shock (RS) would form. The beam material shocked by the RS would expand due to high pressure, which forms the jet cocoon. The stellar envelope materials shocked by the FS would also flow sideways to form the jet cavity. The jet cocoon and jet cavity are called cocoon as a whole in this paper for simplicity. During the jet's expedition to the stellar surface, some further collimating oblique shocks occur due to a cycle of expansion and collimation of the jet (Mizuta & Ioka 2013). Meanwhile, the shear motion between the shocked beam and the shocked envelope materials would produce some turbulent structures in the cocoon. After the jet breaks out of the star, the jet would be accelerated to a higher Lorentz factor (∼5 times of the initial Lorentz factor), and the opening angle gets relatively smaller than the initial opening angle (Morsony et al. 2007; Mizuta & Ioka 2013). Finally, the jet becomes internally free expanding, bounded by a shear layer.

In most of these previous simulations, a steady, constant luminosity has been assumed. In reality, GRB jets are likely intermittent in nature, as revealed by their erratic light curves (Fishman & Meegan 1995). Even though the emission sites and processes that give rise to the observed GRB light curves are still subject to debate (e.g., photospheric emission, Paczynski 1986; Rees & Mészáros 2005; Lazzati et al. 2013; Pe'er & Ryde 2016; internal shocks, Rees & Mészáros 1994; Kobayashi et al. 1997; internal collision-induced magnetic reconnection and turbulence, Zhang & Yan 2011; Zhang & Zhang 2014; Deng et al. 2015), all these models attribute all or part of the variability timescales to the intrinsic intermittent activities (e.g., Fenimore et al. 1999; Yuan & Zhang 2012) or the precession activities (e.g., Ito et al. 2015) at the central engine.

Numerical studies of propagation of intermittent jets in GRB progenitor stars have been carried out by Morsony et al. (2010) and López-Cámara et al. (2014, 2016). They showed that the fast engine variability in the injected jets could be preserved after the jet emerges from the star. López-Cámara et al. (2016) showed that the interaction between the episodic jet and the progenitor stellar envelope can lead to an asymmetric behavior in the light curves of GRBs. Their results suggested a correlation between the duty cycle period of the central engine variability and the breakout time. Furthermore, they pointed out that the general trend and behavior in the three-dimensional (3D) simulations are basically similar to those in the two-dimensional (2D) simulations. However, the envelope model used in their simulations is simplified, with the progenitor star as a background without rotation or infall due to core collapse.

In this paper, we investigate the propagation of relativistic, hydrodynamic, intermittent GRB jets in a star with more realistic physical conditions through numerical simulations. The progenitor star used in our simulations is rotating and collapsing due to gravity. In order to find out how intermittent jets drill out the star, we perform a series of 2.5-dimensional (2.5D, 2D plus rotation) simulations in which jets are injected intermittently from the central engine. For simplicity, the intermittent jets are characterized as a periodic step function with the same active and quiescent durations, T, which is half of the duty cycle period. From these simulations, we intend to study how the parameter T affects the ability of the jet to break out of the star.

This paper is organized as follows. Numerical simulation setup is presented in Section 2. Detailed information for the numerical methods is described in Section 3. In Section 4, we show the numerical results, which are followed by a discussion of the analytical understanding of these results in Section 5. The conclusions are summarized in Section 6.

2.1. The Envelope of a Rotating Star

Our first step is to obtain the structure of a rotating massive star envelope in dynamical equilibrium. According to Woosley & Heger (2006), the 16TI model is one of the promising GRB progenitor models. In this subsection, we try to construct our envelope model to mimic the 16TI model.

Since the envelope of the 16TI model is nearly radiation dominated, a polytropic equation of state p = K ρ^4/3 is adopted, where p is pressure and ρ is density. We consider a simple barotropic star, i.e., K is a constant. After assuming that the star is axisymmetric about its rotational axis, the main task left is to solve the Euler equations with rotation (see Fujisawa 2015)

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial r}=-\displaystyle \frac{\partial \phi }{\partial r}+r{\sin }^{2}\theta {{\rm{\Omega }}}^{2},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial p}{\partial \theta }=-\displaystyle \frac{\partial \phi }{\partial \theta }+{r}^{2}\sin \theta \cos \theta {{\rm{\Omega }}}^{2},\end{eqnarray} \tag{ 2 }$$

and the integral form of the Poisson equation

$$\begin{eqnarray}&&\phi =-G\int \displaystyle \frac{\rho ({{\boldsymbol{r}}}^{\prime })}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}^{\prime }| }{d}^{3}{{\boldsymbol{r}}}^{\prime }.\end{eqnarray} \tag{ 3 }$$

Here p, ρ, the gravitational potential ϕ, and angular velocity Ω are all functions of r and θ, and G is the gravitational constant. Spherical polar coordinates (r, θ) are adopted in the formulation.

Before solving Equations (1)–(3), a one-dimensional (1D) spherical equilibrium nonrotating star is used as the initial guess for the 2D rotating star. We first solve the simple equation

$$\begin{eqnarray}&&\displaystyle \frac{{dp}}{{dr}}=-\displaystyle \frac{{Gm}\rho }{{r}^{2}},\end{eqnarray} \tag{ 4 }$$

where m is the enclosed mass. We integrate this equation from a starting radius r_s (10⁶ cm), using the Runge–Kutta method. The mass inside r_s is set as 1.7 M_⊙, and the density at r_s is set as 6.3 × 10⁹ g cm⁻³, so that the radius of this 1D star is ∼4 × 10¹⁰ cm, and the total mass of the star is ∼14 M_⊙. Using this 1D solution as an initial guess, we then apply Hachisu's Self-Consistent Field method (see Hachisu 1986a, 1986b; Fujisawa 2015 for details) to get the 2D configuration of a rotating star. A realistic star may have different layers obeying different rotation laws due to differential rotation (Kiuchi et al. 2010). Detailed simulations have shown that the distribution of the angular velocity as a function of the mass coordinate may be approximated as a step function (Maeder & Meynet 2012). The angular velocity may be treated as a constant in the core region that encloses a mass M ≤ 7 M_⊙ (Hirschi et al. 2004; see also Figure 16 of Maeder & Meynet 2012). Since our problem mostly concerns the core region of the star where the jet penetrates, for simplicity, we assume rigid rotation, i.e., Ω = constant, throughout the paper to ease the numerical simulation problem. Moreover, the axis ratio q = R_pol/R_eq (R_pol and R_eq are the polar radius and the equatorial radius, respectively) is taken as 0.7, which corresponds to a critical (or breakup) rotation configuration (Maeder & Meynet 2012). Our results give Ω = 5.3 × 10⁻³ rad s⁻¹, which is consistent with the value in Maeder & Meynet (2012) and Nagakura et al. (2011).⁵ One may note that the surface velocity of this star is larger than observed (e.g., Dufton et al. 2011), which is mainly due to the assumption of rigid rotation. However, this would not strongly affect the results of the jet propagation study below since the jet is moving along the polar direction. The upper panel of Figure 1 shows the density distribution of the rotating star, which is used as an input in the next step.

Zoom In Zoom Out Reset image size

2.2. Collapse of the Star

According to the popular GRB central engine models, the GRB jet is launched through the accretion of envelope materials into the central BH. Therefore, as the jet is launched, the star should undergo collapsing after the core collapses to a BH and loses radiation support. In the following, we try to set up a collapsing star into which the GRB jet is injected.

In our simulations, we treat the inner region (r ≤ 4 × 10⁸ cm) of the star as a newly formed BH, of which the mass is ∼2.0 M_⊙. The collapsing of the envelope can be modeled by setting the boundary condition to be a type of outflow (the gradients of all quantities are zero; see also Nagakura et al. 2011). The rarefaction wave triggered at the inner boundary would cause infall of the matter at the radius it reaches. Since the self-gravity module in the PLUTO code (see Section 3 for details) is not fully implemented yet,⁶ we adopt an approximation method to treat the effects of gravity. For a point at radius r, its acceleration in the radial direction is calculated as

$$\begin{eqnarray}&&g=-\displaystyle \frac{{{GM}}_{r}}{{r}^{2}},\end{eqnarray} \tag{ 5 }$$

where ${M}_{r}=4\pi {\int }_{0}^{\pi /2}{\int }_{{r}_{0}}^{r}\rho ({r}^{\prime }){r}^{\prime 2}\sin \theta {{dr}}^{\prime }d\theta +{M}_{\mathrm{BH}}$, r₀ = 4 × 10⁸ cm is the inner boundary, and M_BH is the mass of the central BH. Notice that M_BH is increasing with time due to the infall of matter in our calculation. The justification of this treatment can be seen in the lower panel of Figure 1, where we plot the specific angular momentum j at the innermost stable circular orbit (ISCO) for both a Schwarzschild BH (green) and a Kerr BH (red) as a function of the BH mass, together with the equatorial specific angular momentum j of the star (blue) as a function of the mass enclosed. During the jet-launching phase, the BH in our simulation grows from 2 to ∼6 M_☉. One can see that at the equator, the stellar j is larger than j at ISCO, so that infall is impossible. This is consistent with the formation of an accretion torus that powers the GRB (Woosley & Heger 2006). On the other hand, the j values of the star decrease with decreasing polar angle θ, so that below a critical angle, vertical infall into the BH becomes possible. At the extreme case with θ = 0, j is essentially zero for the stellar materials. In our simulation, when gravity is turned on, steady growth of the BH is observed.

Using the envelope solution obtained in Section 2.1, we simulate the collapsing process using the hydrodynamic module provided in PLUTO and get the structure of the star at a specific time. We then move on to launch a relativistic jet from r₀ starting from this epoch. When the jet is launched, we neglect self-gravity of the envelope, but consider the point-mass gravity of the central BH.

2.3. Jet Propagation

Since our objective is to study jet propagation, we simply inject a relativistic jet via boundary conditions without simulating the jet-launching processes (e.g., Nagataki et al. 2007; Nagataki 2009) from the first principles. The jet is characterized by three parameters, i.e., luminosity (L_j), initial Lorentz factor (Γ₀, corresponding velocity β₀c), and half-opening angle (θ₀). The boundary conditions can be set up correspondingly. We denote density, pressure, and internal energy density of the jet in the fluid frame as ρ_j, p_j, and e_j, respectively, and then the energy density in the lab frame can be expressed as

$$\begin{eqnarray}&&{e}_{\mathrm{lab}}={{\rm{\Gamma }}}_{0}^{2}({\rho }_{j}{c}^{2}+{e}_{j}+{\beta }_{0}^{2}{p}_{j})\simeq {{\rm{\Gamma }}}_{0}^{2}\left(1+\displaystyle \frac{4{p}_{j}}{{\rho }_{j}{c}^{2}}\right){\rho }_{j}{c}^{2},\end{eqnarray} \tag{ 6 }$$

where we have assumed the adiabatic index as 4/3 (e_j = 3 p_j), and β₀ ≈ 1. On the other hand, e_lab is related to L_j by

$$\begin{eqnarray}&&{e}_{\mathrm{lab}}=\displaystyle \frac{{L}_{j}}{2\pi {r}_{0}^{2}(1-\cos {\theta }_{0}){\beta }_{0}c}.\end{eqnarray} \tag{ 7 }$$

When the internal energy of the material is fully converted to kinetic energy, it would reach a terminal Lorentz factor ${{\rm{\Gamma }}}_{\infty }$, which is calculated as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\infty }=\left(1+\displaystyle \frac{4{p}_{j}}{{\rho }_{j}{c}^{2}}\right){{\rm{\Gamma }}}_{0}.\end{eqnarray} \tag{ 8 }$$

Equations (6)–(8) are essential to infer ρ_j, p_j (initial conditions of the jet) given L_j, θ₀, Γ₀, and ${{\rm{\Gamma }}}_{\infty }$.

In this paper, we fix L_j = 5 × 10⁵⁰ erg s⁻¹, Γ₀ = 5, and ${{\rm{\Gamma }}}_{\infty }=400$, unless noted otherwise. The initial opening angle θ₀ is set to vary to investigate its effect on jet propagation.

All the simulations presented here were performed using the PLUTO code, version 4.2. The PLUTO code is built on the Godunov-type shock-capturing schemes, aiming to the solution of Newtonian, relativistic flows (see Mignone et al. 2007, for a full description).

In this work, spherical coordinates (r, θ) are employed and axisymmetry is assumed for all simulations. The collapsing process of the star is simulated using the hydrodynamic (HD) module provided in PLUTO. The star is immersed in a medium with a constant density of 10⁻¹⁰ g cm⁻³. The computational domain covers a region of 4 × 10⁸ cm ≤ r ≤ 1.0 × 10¹¹ cm and 0° ≤ θ ≤ 90°. The radial grid consists of 2048 points and is logarithmically distributed, while the angular grid is uniform with 256 points. A Riemann solver, called HLL solver (Harten-Lax-Van Leer approximate Riemann solver; see Harten et al. 1983), a linear-type spatial reconstruction, and a second-order Runge–Kutta time integration were chosen in the simulations.

Jet propagation is simulated using the relativistic hydrodynamic (RHD) module in PLUTO. Specifically, we choose an extended Harten–Lax–Van Leer Riemann solver (called HLLC solver; see Harten et al. 1983; Mignone & Bodo 2005), a so-called "OSPRE_LIM" flux limiter using a piecewise linear interpolation (Mignone et al. 2005), and a second-order Runge–Kutta time integration. As a result, we achieve second-order accuracy in both space and time. The computational domain covers a region of 4 × 10⁸ cm ≤ r ≤ 1.0 × 10¹¹ cm and 0° ≤ θ ≤ 90°. Adaptive mesh refinement (Mignone et al. 2012) is adopted to improve the computational efficiency. There are altogether six grid levels in our setup, of which the base grid consists of 128 × 64 uniform grid cells, and the refinement ratio is 2. With this setup, the highest resolution of the grid is (Δr, Δθ) = (2.4 × 10⁷ cm, 004), which is comparable to previous 2D studies (Morsony et al. 2007; Mizuta & Ioka 2013).

In order to ensure that our simulation results are reliable, we first perform one baseline simulation to directly compare with the results of Morsony et al. (2007). In this simulation, we use a simple power-law model for a nonrotating star with a mass of 15 M_⊙ and a surface radius of 10¹¹ cm, without considering gravity. All the jet parameters are the same as those of Morsony et al. (2007), i.e., L_j = 5.32 × 10⁵⁰ erg s⁻¹, ${{\rm{\Gamma }}}_{\infty }=400$, and θ₀ = 10°. The inner boundary is set to r = 10⁹ cm, and the grid setup is the same as that in Section 3. The upper panels of Figure 2 show the density distributions at three epochs, t = 9.0, 13.0, and 16.0 s, respectively. It can be noticed that the jet breakout time is ∼13 s, which is almost the same as the value of Morsony et al. (2007). Moreover, the structures of the collimation shock, the cocoon, and the bow shock after the jet breakout time are also in good agreement with previous works. On the other hand, low-resolution spherical coordinates would lead to numerical baryon loading (Mizuta & Ioka 2013). We then diagnosed the resolution of our setup by checking whether Bernoulli's constant hΓ ($h=1+{e}_{j}/({\rho }_{j}{c}^{2})+{p}_{j}/({\rho }_{j}{c}^{2})$ is the specific enthalpy) along the jet axis is conserved. Figure 3 shows that hΓ along the jet axis is conserved up to the jet head, except for some fluctuations due to internal shocks. Also in this figure, the evolution of Γ shows that the core of the unshocked jet is free streaming up to the converging position of the collimation shock. Therefore, using PLUTO, our grid setup is reliable for the purpose of this research.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The star is already collapsing before jet launching. However, the exact jet-launching delay time is not known. We choose a moderate value for this delay time, i.e., 24 s, in our simulations. Figure 4 shows the density profiles in the equatorial and axial directions for t = 0 and 24 s after the collapse, respectively. A shock wave produced by centrifugal bounce is predicted to emerge during the collapse (Nagakura et al. 2011). This shock wave does not appear in Figure 4. This may be caused by the approximate treatment used here on self-gravity (Equation (5)). However, this would not affect much our following results, since the primary density profile in Figure 4 is still similar to that in Nagakura et al. (2011). Using this envelope as the initial condition, we perform two series of simulations for intermittent jets. The jets are characterized by periodic injections with top-hat constant luminosity in each episode and with the same duration during the on and off states. We simulate four intermittent jets with half period T = 0.2, 1.0, 2.0, and 4.0 s, respectively. For comparison, we also simulate a baseline constant-luminosity jet (see Figure 5). For each case, we perform two simulations with the jet opening angle being 5^o and 10^o, respectively. The simulation models are named in the form of "mN1TN2," with N1 denoting the value of θ₀ and N2 denoting the value of T. For constant-luminosity jets, N2 is replaced by "con." The detailed parameters for all these models are listed in Table 1. The density distributions of the m10Tcon model at three epochs, t = 3.0, 4.4, and 6.2 s, respectively, are presented in the lower panels of Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Initial Conditions of the Episodic Jet Models

Model	T	θ0	Γ0	Luminosity (Lj)	Inner Boundary Radius	${{\rm{\Gamma }}}_{\infty }$	tbo
	(s)			(erg s−1)	(r0)		(s)
m5Tcon	0.0a						4.2
m5T0.2	0.2						4.86
m5T1.0	1.0	5°	5.0	5 × 1050	4 × 108 cm	400	5.8
m5T2.0	2.0						6.1
m5T4.0	4.0						4.45
m10Tcon	0.0						4.4
m10T0.2	0.2						6.47
m10T1.0	1.0	10°	5.0	5 × 1050	4 × 108 cm	400	6.0
m10T2.0	2.0						6.16
m10T4.0	4.0						4.88

Note.

^aT = 0.0 means that the jet is continuous.

Download table as:  ASCIITypeset image

From these simulations, we investigate the relationship between the breakout time (t_bo) and the half period (T), which is presented in Figure 6. It is interesting to note that at T = 0.2 s, t_bo is significantly different for the two different θ₀ values, while for other values of T, t_bo depends weakly on θ₀. In comparison with the results of López-Cámara et al. (2016), under the same initial condition (θ₀ = 10°), the trend of our t_bo–T relation is similar to the 3D result in López-Cámara et al. (2016), but is different from their 2D result. We should note that the oblate configuration of the envelope is the main reason for the difference of the absolute values of t_bo between our results and the results in López-Cámara et al. (2016).

Zoom In Zoom Out Reset image size

Figures 7 and 8 show the snapshot Lorentz factor maps (color) and density distributions (gray) for two sets of intermittent jet simulations. The cases of constant-luminosity jets are presented in Figure 5 for comparison. The epoch for the snapshot for each model roughly corresponds to the time when the first episode of the jet reaches the simulation boundary. To identify fast-moving ejecta, only regions with Γ > 2 are shown in color. From Figures 7 and 8, one can see spatially separated high-Γ regions that manifest the intermittent injection of the central engine.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One most striking finding from our simulations is that intermittent jets with a relatively short half period do not have every subpulse emerging from the star. After the first pulses escape the star at t_bo, some subpulses vanish during their propagation inside the star. Here we define "jet vanishing" as the Lorentz factor dropping below 2. For example, for model m10T0.2, after t_bo (6.47 s), the next jet pulse that can reach the outer boundary breaks out of the star surface at ∼11.6 s. For model m10T1.0, after t_bo (6.0 s), only the seventh pulse (ejected during 12–13 s) and later ones are able to reach the outer boundary. Figure 9 gives some details of the m10T1.0 simulation. The top panel (t = 7.2 s) shows that the first three pulses (p1-3) merge and break out of the star at Y ≥ 150, and that the fourth pulse (p4) was just ejected from the engine (head at Y ≥ 40). In the next two panels (t = 8.2 and 9.2 s), one can see that p4 gradually vanishes. The fifth pulse (p5) seen in the t = 9.2 s panel and the sixth pulse (p6) seen in the t = 11.2 s panel all gradually vanish before penetrating through the star. This vanishing continues until pulse 7 (p7) is ejected. One can see that p7 manages to breakout of the star again in the t = 15.6 s panel.

Zoom In Zoom Out Reset image size

For relatively large half-period intermittent jets (T = 2, 4 s), the subpulse vanishing events do not happen. Every subpulse can manage to break out of the star separately.

The numerical results can be understood using some analytical treatments.

Let us first consider the propagation of a continuous jet. Our numerical results (m5Tcon and m10Tcon) can be directly compared with the results of the analytical model (e.g., Bromberg et al. 2011). The jet head velocity is given by

$$\begin{eqnarray}&&{\beta }_{h}=\displaystyle \frac{{\beta }_{j}}{1+{\tilde{L}}^{-1/2}},\end{eqnarray} \tag{ 9 }$$

where $\tilde{L}={L}_{j}/{{\rm{\Sigma }}}_{j}{\rho }_{a}{c}^{3}$, ${\rho }_{a}(r,\theta =0)$ is the ambient medium density in the polar direction, and Σ_j is the jet cross section and can be further derived by considering the pressure balance between the jet and the cocoon. In order to get the analytical formulae for Σ_j, four coefficients, ξ_a, ξ_h, ξ_c, and η, are introduced to simplify the derivation (see Bromberg et al. 2011; Mizuta & Ioka 2013, for details). In short, η is a parameter to correct the approximation of the cylindrical cocoon shape, ξ_a represents the correction to the mean density of the medium in the cocoon, and ξ_h and ξ_c are the corrections to the approximations of jet head position and cocoon width, respectively. When the density profile in the poloidal direction is a power law ρ_a ∝ z^−α, ξ_a, ξ_h, ξ_c can be expressed as ${\xi }_{a}=\tfrac{3}{3-\alpha }$, ${\xi }_{h}={\xi }_{c}=\tfrac{5-\alpha }{3}$. After introducing these coefficients, the jet head velocity and jet head position can be finally calculated as

$$\begin{eqnarray}&&{\beta }_{h}\simeq {\tilde{L}}^{1/2}\simeq {\left(\displaystyle \frac{{L}_{j}}{{c}^{5}{t}^{2}{\rho }_{a}{\theta }_{0}^{4}}\right)}^{1/5}{\left(\displaystyle \frac{16}{3\pi }\displaystyle \frac{\eta {\xi }_{a}}{{\xi }_{h}{\xi }_{c}^{2}}\right)}^{1/5},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}{z}_{h} & = & {\xi }_{h}{\beta }_{h}{ct}\simeq {\left(\displaystyle \frac{{t}^{3}{L}_{j}}{{\rho }_{a}{\theta }_{0}^{4}}\right)}^{1/5}{\left(\displaystyle \frac{16}{3\pi }\displaystyle \frac{\eta {\xi }_{a}{\xi }_{h}^{4}}{{\xi }_{c}^{2}}\right)}^{1/5}\\ & \simeq & 1.4\times {10}^{10}\,\mathrm{cm}{\left(\displaystyle \frac{t}{1{\rm{s}}}\right)}^{3/5}{\left(\displaystyle \frac{{L}_{j}}{{10}^{51}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{1/5}\\ & & \times {\left(\displaystyle \frac{{\rho }_{a}}{{10}^{3}{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{-1/5}{\left(\displaystyle \frac{{\theta }_{0}}{0.1}\right)}^{-4/5}{\left(\displaystyle \frac{\eta }{0.01}\right)}^{1/5},\end{eqnarray} \tag{ 11 }$$

where α = 2 is assumed (which gives ξ_a = 3, ξ_h = ξ_c = 1). By adjusting the parameter η, we can fit the numerical results with the analytical solution in Equation (11). Figure 10 shows the comparison between the analytical solution (η = 0.02 for θ₀ = 5° and η = 0.006 for θ₀ = 10°) and the numerical results. The numerical results are in good agreement with the analytical results. In this figure, the analytical solutions are obtained using Equation (11) with the density slope α = −dlnρ_a/dlnz at z = z_h.

Zoom In Zoom Out Reset image size

For T = 0.2 s intermittent jets, t_bo in m5T0.2 and m10T0.2 are significantly different from each other. The narrow jet (m5T0.2) travels much faster (t_bo = 4.86 s) than the wide jet (m10T0.2) (t_bo = 6.47 s). This is caused by two reasons. First, for a smaller θ₀, the energy is deposited on a smaller area so that the jet is more penetrating (jet pulse is narrower and moves faster) in the polar direction. Second, the breakout in these short-T models is triggered by the expanding cocoon. In these models, stellar materials are shocked into a hot cocoon matter by the vanishing jet pulses. When a bow shock forms at the stellar surface, there is only a small amount of jet matter (with high Γ) near the surface, so that the outgoing materials are mostly from the cocoon. The larger the θ₀, the wider the cocoon (see Figure 11), and the smaller the expanding velocity in the polar direction.

Zoom In Zoom Out Reset image size

For longer duty cycle jets, i.e., T ≥ 1.0 s, t_bo is not sensitive to θ₀. This is because the breakout of these cases is triggered by the jet pulses rather than the cocoon. As a result, given the same T, the jet breakout times for different θ₀ values are defined by the same ejection time from the central engine. For example, in m5T2.0 and m10T2.0, the second jet pulse ejected during 4–6 s is powerful enough to lead to the breakout in both cases.

The jet pulse vanishing in m10T0.2 and m10T1.0 but not in m10T2.0 and m10T4.0 can be understood. Like a continuous jet, a jet head is also formed to evacuate the envelope matter in front of the jet pulse. The velocity of the jet head is smaller than that of the jet itself, so that the jet materials would flow across the RS into the jet head, and finally flow into the cocoon. If the jet pulse does not last for a long time, all the jet materials would go into the cocoon before the breakout of the jet head. For a jet pulse of an active duration of T, we assume that the average velocity of the jet head is ${\bar{\beta }}_{h}$. It would vanish after a duration

$$\begin{eqnarray}&&\tau \simeq \displaystyle \frac{{\bar{\beta }}_{h}T}{1-{\bar{\beta }}_{h}}=\displaystyle \frac{T}{1/{\bar{\beta }}_{h}-1}.\end{eqnarray} \tag{ 12 }$$

By equaling ${\bar{\beta }}_{h}c\tau$ and R_pol, one can obtain a critical duration

$$\begin{eqnarray}&&{T}_{c}=\displaystyle \frac{1}{{\bar{\beta }}_{h}}\left(\displaystyle \frac{1}{{\bar{\beta }}_{h}}-1\right)\displaystyle \frac{{R}_{\mathrm{pol}}}{c},\end{eqnarray} \tag{ 13 }$$

below which the jet crosses the jet head before penetrating through the star, so that it would vanish. In our calculations, R_pol ≃ 2.8 × 10¹⁰ cm is adopted. Our simulations show that ${\bar{\beta }}_{h}\simeq 0.6$ could be reached for each pulse. One then gets T_c ∼ 1.0 s. As a result, a jet pulse would be able to drill out the stellar surface and reach the outer boundary if T > T_c ∼ 1.0 s is satisfied. This is consistent with the fact that jet vanishing does not appear in m5T2.0 or m10T2.0. Another condition for a successful jet pulse is that the next jet would not catch up with the prior jet within the duration τ, so a quiescent duration should not exceed $(1-{\bar{\beta }}_{h})\tau$. Both conditions are satisfied for the T = 2.0 and 4.0 s cases, but not satisfied for the T = 0.2 and 1 s cases (see Figure 12). This explains the jet pulse vanishing phenomenon in our simulations with small T.

Zoom In Zoom Out Reset image size

The collapsing envelope also plays a key role in destroying jet pulses. In comparison with m10T0.2, we run a simulation in which the star is not rotating or collapsing and the jet is launched at r₀ = 10⁹ cm (a common value used by other authors). It is found that the jet vanishing does not exist in this simulation (see Figure 13). This difference also holds if we move the injection radius in m10T0.2 from 4 × 10⁸ cm to 10⁹ cm. In Figure 13, we can see that the cocoon structure without gravity is obviously wider than that in m10T0.2. This in turn gives a larger pressure in the hot cocoon, which makes narrower jet pulses and stronger internal shocks to dissipation jet energy (Figure 14). Even though gravity does not affect jet dynamics directly, our results suggest that gravity in addition to rotation may influence the cocoon structure and hence the jet dynamics.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The jet pulse vanishing phenomenon in m10T0.2 or m10T1.0 indicates that there may be a time delay between the first breakout and the breakouts of subsequent jet pulses. In the case of m10T1.0, the time delay can reach ∼7 s. This value could be even larger, if the engine quiescent time duration is even longer than the active time (in our simulation we assumed the same duration). This motivates us to relate our simulation results to the light-curve quiescent time between the GRB precursor and its main emission. It has been found that ∼15% of GRBs have a precursor (Burlon et al. 2008, 2009). The mean duration of the quiescence is ∼50 s, of which the origin is still uncertain (Mészáros & Rees 2000; Ramirez-Ruiz et al. 2002; Lazzati & Begelman 2005; Wang & Mészáros 2007; Hu et al. 2014). Our results can provide a reasonable explanation for such a quiescent time. According to our picture, the GRB central engine may continuously eject episodic jets with short active timescales. Initially, the several leading subpulses merge and break out of the star with the mixed cocoon and jet matter, which produces the precursor. Later, for a certain duration, the episodic jet pulses injected from the central engine vanish since the vanishing time τ is shorter than the jet-penetrating time. The GRB therefore displays a long gap of quiescence. Finally, the funnel near the jet axis becomes more and more hollow; some episodic jets of long active periods can eventually escape the star to make radiation, making the main episode of the burst. According to our simulations, pulse vanishing happens only if T is around or shorter than 1 s. Observationally, the shortest variability timescale for long GRBs is of the similar order. This also lends support to our model. Recently, Lloyd-Ronning et al. (2016) showed that within the context of a magnetically arrested disk around a BH, the characteristic timescale for the variability of the Poynting flux jet is ∼1 s. This also provides a physical standpoint for our scenario, i.e., the jet pulses are characterized by a duty cycle period of the order of a second. However, we have not considered the role of magnetic fields on jet propagation in this work.

There are two facts that support our proposal. The precursors and main events have very similar spectral properties (Burlon et al. 2009; Hu et al. 2014), which implies that they share similar radiation physics. In our explanation, the precursor and the later main bursts are all from episodic jet pulses ejected from the same central engine, so that their emission properties should be similar. Second, multiple precursors have been observed in many cases (e.g., 73 cases in the Burlon et al. [2009] sample). This fact is also a natural outcome in our scenario. After the initial breakout to produce the first precursor, in some GRBs, the episodic jets may have relatively long active durations, so that they can subsequently break out of the star to produce more precursors before the main burst comes. Some others may have relatively short active time durations, so that they can only produce one precursor with a long quiescent gap before the main burst comes after the funnel is cleared.

Another interesting finding in our simulations is that later-injected pulses are likely to catch up with the pulses injected earlier. GRB internal shock models require that late shells catch up with the early ones. In our simulations, even if all the pulses have the same initial energy and internal energy, later-injected jet pulses are accelerated more quickly when the episodic jets propagate in the collapsing envelope. This may be understood as follows. First, the group of pulses that break out of the star at t_bo and make the precursor would evacuate a funnel along the jet axis direction, which cannot be refilled quickly. As a consequence, pulses injected at later epochs would travel in a more stratified medium and therefore can be accelerated more quickly. In m5T0.2, later pulses can even catch up with the early ones within the simulation box (see Figure 15).

Zoom In Zoom Out Reset image size

In this paper, we numerically modeled the propagation of a relativistic, hydrodynamic, intermittent jet in the envelope of a massive, rotating, and collapsing star for the first time. Using the PLUTO code with specific setups, we performed a set of 2.5D simulations. Because of the introduction of a more realistic rotating star and central gravity during the jet propagation phase, the following interesting results (some of which were not obtained previously) are obtained.

For intermittent jets, the shock breakout time t_bo depends on the jet opening angle when the half period T is short enough (e.g., T = 0.2 s). This is because for short T, the RS to each pulse crosses the shell before it penetrates through the star, so that the shock breakout is defined by that of a hot cocoon, the speed of which depends on the opening angle of the jet. When the half period is large enough (T ≥ 1 s), the RS crossing time becomes comparable to or longer than the jet-penetrating time, so that t_bo no longer sensitively depends on the jet opening angle. The t_bo–T relation trend obtained from our 2.5D results is similar to that from 3D results in López-Cámara et al. (2016).

For short half-period jets (T = 0.2, 1 s), our simulations revealed some vanishing pulses after the first group of pulses break out of the star. The vanishing of subpulses is due to two reasons: that the RS crosses the pulse before it penetrates the star, and that a later pulse may have caught up with a preceding jet. Such a vanishing effect gives rise to a quiescent time after the initial shock breakout even if the engine is continuously ejecting intermittent pulses, and can be a mechanism to interpret the quiescent gap between the precursor and the main burst. Such a behavior is only observed when a rotating, collapsing star is simulated.

Growing evidence suggests that GRB prompt emission jets may carry a significant Poynting flux (e.g., Zhang & Kumar 2013; Uhm & Zhang 2014, 2016). In the future, we plan to introduce strong magnetization in relativistic jets and numerically investigate the effect of strong magnetization on jet propagation. A more self-consistent setup of the collapsar including the magnetic fields may be also considered (Mösta et al. 2014, 2015).

We thank the anonymous referee for valuable suggestions and Kotaro Fujisawa, Davide Lazzati, Diego López-Cámara, Andrea Mignone, Brian Morsony, Philipp Mösta, and Shigehiro Nagataki for helpful comments on the paper. This work is partially supported by NASA under grants NNX15AK85G and NNX14AF85G, by the National Basic Research Program of China with grant no. 2014CB845800, and by the National Natural Science Foundation of China (grant nos. 11473012 and 11303013). J.-J.G. acknowledges the China Scholarship Program to conduct research at UNLV. R.K. acknowledges financial support within the Emmy Noether research group on "Accretion Flows and Feedback in Realistic Models of Massive Star Formation" funded by the German Research Foundation under grant no. KU 2849/3-1. The visualization of the simulation data in this article is conducted with software VisIt (Childs et al. 2005). This work made use of the Cherry Creek cluster at UNLV's National Supercomputing Institute.