# RELATIVISTIC JETS FROM COLLAPSARS

## 1 INTRODUCTION

Catastrophic collapse events have been proposed to explain the energies released in a gamma-ray burst (GRB) including mergers of compact binaries (Pacyński 1986; Goodman 1986; Eichler et al. 1989; Mochkovitch et al. 1993), collapsars (Woosley 1993) and hypernovae (Pacyński 1998). According to the current view these models require a stellar mass black hole (BH) which accretes up to several solar masses of matter powering a pair fireball. If the baryon load of the fireball is not too large, baryons are accelerated together with e$`^+`$e<sup>-</sup> pairs to Lorentz factors $`>10^2`$ (Cavallo & Rees 1978). Such relativistic flows are supported by radio observations of GRB 980425 (Kulkarni et al. 1998b). Spherically symmetric fireballs have been studied by several authors by means of 1D Lagrangian hydrodynamic simulations (e.g., Panaitescu et al. 1997; Panaitescu & Mészáros 1998; Kobayashi, Piran & Sari 1999). Recently, it has been argued that the rapid temporal decay of several GRB afterglows is more consistent with the evolution of a relativistic jet after it slows down and spreads laterally than with a spherical blast wave (Sari, Piran & Halpern 1999; Halpern et al. 1999; Kulkarni et al. 1999; Rhoads 1999). The lack of a radio afterglow in GRB 990123 provides independent evidence for jet-like geometry (Kulkarni et al. 1999b).

## 2 INITIAL MODEL AND NUMERICAL SETUP

MacFadyen & Woosley (1999; MW99) have explored the evolution of rotating helium stars ($`M_\alpha >\text{ }\mathrm{\hspace{0.17em}10}M_{}`$) whose iron core collapse does not produce a successful outgoing shock, but instead forms a BH surrounded by a compact accretion torus. Assuming an enhanced efficiency of energy deposition in polar regions MW99 obtain relativistic jets along the rotation axis, which are highly focused and seem to be capable of penetrating the star. However, as their simulations are Newtonian, they obtain flow speeds which are superluminal.

We have performed axisymmetric relativistic simulations using a $`14M_{}`$ collapsar model from MW99. When the central BH has acquired a mass of $`3.762M_{}`$ we map the model to our computational grid. In a consistent collapsar model a jet will be launched by any process which gives rise to a local deposition of energy and/or momentum, as e.g.,$`\nu \overline{\nu }`$–annihilation, or magneto-hydrodynamic processes. We mimic such a process by depositing energy at a prescribed rate homogeneously within a $`30^{}`$ cone around the rotation axis. In the radial direction the deposition region extends from the inner grid boundary located at 200 km to a radius of 600 km. We have investigated constant energy deposition rates $`\dot{E}=10^{50}`$erg s<sup>-1</sup>, and $`\dot{E}=10^{51}`$erg s<sup>-1</sup>, and a varying deposition rate with a mean value of $`10^{50}`$erg s<sup>-1</sup>. The constant rates roughly bracket the expected $`\dot{E}`$ of collapsar models, while the varying rate mimics, e.g., time-dependent mass accretion rates resulting in time-dependent $`\nu \overline{\nu }`$–annihilation (MW99).

The simulations were performed with the multidimensional relativistic hydrodynamic code GENESIS (Aloy et al. 1999) using a 2D spherical grid with 200 radial zones spaced logarithmically between the inner boundary and the surface of the helium star at $`R_{}=2.98\times 10^{10}`$cm. Assuming equatorial symmetry we performed simulations with $`2^{},1^{}`$ and $`0.5^{}`$ angular resolution. In the latter case the grid consists of 60 uniform zones covering the polar region ($`0^{}\theta 30^{}`$) and 40 nonuniform zones logarithmically distributed between $`30^{}\theta 90^{}`$. The gravitational field of the BH is described by the Schwarzschild metric. Effects due to the self-gravity of the star are neglected, i.e., we consider only the gravitational potential of the BH. The equation of state includes non-relativistic nucleons treated as a mixture of Boltzmann gases, radiation, and an approximate correction due to $`e^+e^-`$–pairs as described in Witti, Janka & Takahashi (1994). Complete ionization is assumed, and the effects due to degeneracy are neglected. We advect nine non-reacting nuclear species which are present in the initial model: C<sup>12</sup>, O<sup>16</sup>, Ne<sup>20</sup>, Mg<sup>24</sup>, Si<sup>28</sup>, Ni<sup>56</sup>, He<sup>4</sup>, neutrons and protons.

## 3 RESULTS

### 3.1 Constant small energy deposition rate (Model C50)

For a constant $`\dot{E}=10^{50}`$erg s<sup>-1</sup> a relativistic jet forms within a fraction of a second and starts to propagate along the rotation axis with a mean speed of $`7.8\times 10^9`$ cm s<sup>-1</sup> (Fig. 1). The jet exhibits all the morphological elements of the Blandford & Rees (1974) jet model in the context of classical double radio sources: a terminal bow shock, a narrow cocoon, a contact discontinuity separating stellar and jet matter, and a hot spot. Fig. 1 shows that the density structure of the star does not change noticeably during the whole evolution. This, a posteriori, justifies our neglect of the self-gravity of the star.

The propagation of the jet is unsteady, because of density inhomogeneities in the star. The Lorentz factor of the jet, $`\mathrm{\Gamma }`$, increases non-monotonically with time, while the density drops to $`10^6`$g cm<sup>-3</sup> (Fig. 2). The density profile shows large variations (up to a factor of 100) due to internal shock waves. The mean density in the jet is $`10^1`$ g cm<sup>-3</sup>. Some of the internal biconical shocks, which develop during the jet’s propagation, recollimate the beam. They may provide the “internal shocks” proposed to explain the observed gamma-ray emission. A particularly strong recollimation shock wave (hardly evident at low resolution) forms early in the evolution. A very strong rarefaction wave behind this recollimation shock causes the largest local acceleration of the beam material giving rise to a maximum in the Lorentz factor. When the jet encounters a region along the axis where the density gradient is positive (at $`\mathrm{log}r8.1`$ and $`\mathrm{log}r8.6`$) the jet’s head is decelerated, while a central channel in the beam is cleaned by outflow into the cocoon through the head, which accelerates the beam. The combination of both effects (deceleration of the head and beam acceleration) increases the strength of the internal shocks. Within the jet the mean value of the specific internal energy, $`10^{20}-10^{21}`$ erg g<sup>-1</sup>, or $`O(c^2)`$. The mean temperature $`5\times 10^8`$ K (well below the pair creation threshold) implying that the pressure is radiation dominated in accordance with our simplified EOS.

The relativistic treatment of the hydrodynamics leads to a qualitatively similar (formation of a jet), but quantitatively very different evolution than in MW99. According to their Fig. 27 the jet propagates 7 000 km within the first 0.82 s. Furthermore, MW99 infer an asymptotic $`\mathrm{\Gamma }10`$, and find a half opening angle, $`\mathrm{\Omega }`$, for their jet of $`10^{}`$. In our simulation, at the same time for the same angular resolution ($`2^{}`$) and $`\dot{E}`$ the head reaches a radius of 30 000 km, but the maximum Lorentz factor ($`\mathrm{\Gamma }_{\mathrm{max}}`$) is only 4.62 at $`\mathrm{12\hspace{0.17em}200}`$ km. Such quantitative difference is expected, in part, due to different mapping time and inner boundary radius between the two calculations. Initially, in our simulations $`\mathrm{\Omega }`$ is between $`6^{}`$ to $`8^{}`$ depending on angular resolution. At $`1.5`$s the strong recollimation shock reduces $`\mathrm{\Omega }<\text{ }1^{}`$.

We find that some results strongly depend on angular resolution, the minimum acceptable one being $`0.5^{}`$ (at least near the axis). The morphology of the jet is richer at higher resolution. At $`0.5^{}`$ angular resolution $`\mathrm{\Gamma }_{\mathrm{max}}15-20`$ at a radius $`8\times 10^9`$ cm at jet breakout. Within the uncertainties of the jet mass determination due to finite zoning and the lack of a precise numerical criterion to identify jet matter, the baryon load ($`\eta Mc^2/E_{\mathrm{depos}}`$ with $`E_{\mathrm{depos}}=\dot{E}𝑑t`$) decreases with increasing resolution. In the highest resolution run we find an average baryon load of $`\eta 1.3`$ at jet breakout (see also Sect. 4).

### 3.2 Constant large energy deposition rate (Model C51)

Enhancing $`\dot{E}`$ by a factor of ten (to $`10^{51}`$ erg s<sup>-1</sup>), the jet flow reaches larger Lorentz factors. We observe transients during which the Lorentz factor becomes as large as 40. After 1.2 s the Lorentz factor steadily increases from 22 to 33. The jet propagates faster than in model C50. The time required to reach the surface of the star is 2.27 s instead of 3.35 s. At breakout the jet is less collimated ($`\mathrm{\Omega }10^{}`$). The strong recollimation shock present in model C50 is not so evident here. Instead, several biconical shocks are observed within a very knotty beam and the Lorentz factor near the head of the jet is larger ($`22`$ in the final model) because, due to the larger $`\dot{E}`$, the central funnel is evacuated faster, and because the mean density of the jet is 5 times smaller than in model C50 ($`\eta `$ being twice as large).

### 3.3 Varying energy deposition rate (Model V50)

We have computed a model where the mean energy deposition rate ($`10^{50}`$ erg s<sup>-1</sup>) randomly varies on time scales of a few milliseconds and the amplitude by a factor of ten. Compared to model C50 the jet structure is more knotty and also richer in shocks, particularly inside the first $`10^9`$ cm (where the radial resolution is large enough to capture the finest structures imprinted on the flow by the time variability of the deposition rate). At breakout $`\mathrm{\Gamma }_{\mathrm{max}}=26.81`$, which is almost twice as large as the one found in model C50. Thus, a variable $`\dot{E}`$ is more efficient in converting internal energy into kinetic energy, and in this case, the internal shocks are stronger and more numerous. The mean propagation speed is similar in both models, although the instantaneous velocity of the jet’s head is clearly different. Behind the strongest recollimation shock $`\mathrm{\Omega }<1^{}`$ in both models.

### 3.4 Evolution after jet breakout

The structure of the circumstellar medium will influence the characteristics of the GRB and of the subsequent afterglow. Thus, a continuation of the simulations beyond jet breakout is necessary. In order to satisfy the conditions for accelerating shocks (Shapiro 1979) we endowed the star with a Gaussian atmosphere, which at $`R_\mathrm{a}=1.8R_{}`$ passes over into an external uniform medium with a density $`10^5`$ g cm<sup>-3</sup> and a pressure $`10^8p(R_{})`$. The computational domain is extended to $`R_t=2.54R_{}`$ with 70 additional zones. The evolution after jet breakout has been computed for models C50 and C51. In both cases the jet reaches $`R_t`$ after $`1.8`$ s (measured from breakout). Its mean propagation velocity is $`0.85`$ c, which is almost three times faster than the velocity of the head inside the star (0.30 c in model C50; 0.44 c in model C51).

The evolution after jet breakout consists of three distinct epochs (Fig. 3). The first one lasting 0.35 s is characterized by a head velocity of 0.48 c and a small sideways expansion. During the second phase (of 0.3 s) the jet head accelerates to 0.91 c, because of the steep external density gradient, and because the flux of axial momentum is still important compared to pressure. The sideways expansion is still sub-relativistic ($`0.008`$ c), and $`\mathrm{\Omega }`$ of the beam increases to $`10^{}`$. During the final 1.2 s the bow-shock propagates within the uniform part of the ambient medium leading to a rapid ($`\mathrm{\Gamma }5`$) lateral spreading (Fig. 3).

The shape of the expanding bubble is prolate (Figs. 1 and 3) during the post-breakout evolution. However, when the jet reaches the uniform part of the circumstellar environment, the bubble widens due to the faster sideways expansion. We expect a more isotropic expansion when most of the bubble is inside the uniform medium, and when it is pressure driven (in particular if the energy deposition is switched off). The Lorentz factor near the boundary of the cavity blown by the jet grows from $`1`$ (at jet breakout) to $`3`$ in both models decreasing with latitude. At the end of the simulation $`\mathrm{\Gamma }_{\mathrm{max}}`$ is 29.35 (44.17) for model C50 (C51), which is still smaller than the ones required for the fireball model (Piran 1999).

## 4 DISCUSSION AND CONCLUSIONS

Energy deposition in the polar regions of a collapsar model gives rise to both the formation and propagation of a relativistic jet through the mantle and envelope of the star, and to a supernova explosion. The jet has a small opening angle ($`8^{}`$) and possesses a highly collimated ($`1^{}`$), ultra-relativistic core in which the Lorentz factor reaches a value of $`\mathrm{\Gamma }_{\mathrm{max}}=44`$ (model C51) at the end of the simulation about 2 s after shock breakout. The equivalent isotropic kinetic energy (see MacFadyen, Woosley & Heger 1999) slightly exceeds $`10^{54}`$ erg for model C51 ($`10^{53}`$ erg for model C50) within $`2^{}`$ ($`5^{}`$) of the rotation axis dropping by a factor of 10 within $`17^{}`$ ($`10^{}`$). The inner region contains $`8\times 10^4M_{}`$ with $`<\mathrm{\Gamma }>4`$ (Fig. 1). For a larger $`\dot{E}`$, the jet and in particular the cocoon are less collimated, because a harder driven jet also expands stronger laterally.

The rest-mass density and the internal energy strongly vary in space and time within the jet giving rise to a very inhomogeneous baryon load, i.e., the concept of $`\eta `$ as a global parameter is useless. Instead it is more appropriate to discuss the efficiency of energy conversion in terms of incremental baryon loads considering only matter within a given range of $`\mathrm{\Gamma }`$–values. Although we find an average baryon load of the jet of $`\overline{\eta }1`$, some parts of the flow have a baryon load as low as $`10^5`$ or even less. After jet breakout $`\overline{\eta }`$ decreases by a factor 4 in less than 1.8 s. If this trend continues even $`\overline{\eta }10^3`$ within 9 s.

In model C51, at the end of the simulation, $`2M_{}`$ have a Lorentz factor of less than three, $`3\times 10^4M_{}`$ move with $`3\mathrm{\Gamma }<10`$, and for $`2\times 10^6M_{}`$ the Lorentz factor $`\mathrm{\Gamma }10`$ (Fig. 4). The latter two masses reduce to $`2\times 10^5M_{}`$ and $`2\times 10^7M_{}`$ for model C50. Except for the very early evolution ($`t<1`$ s) the amount of matter moving at moderate ($`3\mathrm{\Gamma }<10`$) and highly ($`\mathrm{\Gamma }10`$) relativistic velocities increases by a factor $`3`$ every second, i.e., if the central engine is active for another 5 s, at the assumed energy deposition rate, $`10^4M_{}`$ will move with $`\mathrm{\Gamma }10`$. As $`\mathrm{\Gamma }_{\mathrm{max}}`$ is also rapidly increasing, Lorentz factors of several hundreds might be reached before the central engine is switched off. For models where the released total energy is equal (C50 and V50) $`\mathrm{\Gamma }_{\mathrm{max}}`$ is higher (by a factor of two) for a time-dependent $`\dot{E}`$. Determining the efficiency of energy conversion $`E_k/E_{\mathrm{depos}}`$ is hampered by the fact, that the kinetic energy $`E_k`$ of a relativistic fluid is not a well defined quantity. If $`E_k\rho \mathrm{\Gamma }(\mathrm{\Gamma }-1)𝑑V`$, we find $`=1.6`$ for model C50 (and $`=2.8`$ for model C51) at the end of the simulation. These efficiencies are obtained considering only matter with radial velocities $`>0.3c`$ and specific internal energy densities $`>5\times 10^{19}`$ ergg<sup>-1</sup>, i.e., matter in the jet. Note that efficiencies larger than one can arise, because (i) there are also other large sources of energy (e.g., gravitational, internal) available, and (ii) because matter is entrained into the jet which does not originate from the deposition region. Thus, efficiencies larger than one suggest that the local energy deposition is efficient triggering conversion of other forms of energy into kinetic energy.

In our simulations the jet reaches the stellar surface intact (propagating over three decades in radius). This result may also hold for other less specific initial conditions. A more spherical density stratification might decrease the collimation of the jet, but the outflow might also be initiated mostly by momentum deposition instead of pure energy deposition (e.g., by MHD effects). The propagation of the jet after breakout will depend on the density stratification of the circumstellar medium. Thus, further simulations with different environments are planned. We note in this respect that the post–breakout propagation is similar in models C50 and C51 suggesting that a lower value of $`\rho _{\mathrm{ext}}`$ will not change the dynamics quantitatively. In our models the jet has only reached a radius of $`7.5\times 10^{10}`$ cm at the end of our simulations, which is $`10^2`$ to $`10^4`$ times smaller than the distance at which the fireball becomes optically thin. Determining whether a GRB will eventually be produced requires to compute the further evolution of the jet. As the jet has stayed collimated in the star it might remain focused over the next three decades, too.

This work has been supported in part by the Spanish DGES (grant PB97-1432) and the CSIC. MAA expresses his gratitude to the Conselleria d’Educació i Ciència de la Generalitat Valenciana for a fellowship. We would like to thank Stan Woosley for his enthusiasm in promoting this work, and Thomas Janka for many helpful remarks. The calculations were performed on two SGI Origin 2000 at the CEPBA and at the SIUV.