# Relativistic Jets from Collapsars

## 1 Motivation and numerical setup

Various catastrophic collapse events have been proposed to explain the energies released in a gamma–ray burst (GRB) including compact binary system mergers , collapsars  and hypernovae . These models all rely on a common engine, namely a stellar mass black hole (BH) which accretes several solar masses of matter from a disk (formed during a merger or by a non–spherical collapse). A fraction of the gravitational binding energy released by accretion is converted into a pair fireball. Provided the baryon load of the fireball is not too large, the baryons are accelerated together with the e$`^+`$e<sup>-</sup> pairs to ultra–relativistic speeds (Lorentz factors $`>10^2`$; ). The existence of such relativistic flows is supported by radio observations of GRB 980425 .

The dynamics of spherically symmetric relativistic fireballs has been studied by several authors by means of 1D Lagrangian hydrodynamic simulations (e.g.,). 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 . The lack of a significant radio afterglow in GRB 990123 provides independent evidence for jet–like geometry . Motivated by these observations and by the collapsar model of , we have simulated the propagation of jets from collapsars using relativistic hydrodynamics.

In  the continued evolution of rotating helium stars, whose iron core collapse does not produce a successful outgoing shock but instead forms a BH surrounded by a compact accretion disk, has been explored. Assuming that the efficiency of energy deposition by $`\nu \overline{\nu }`$–annihilation or, e.g., magneto-hydrodynamic processes is higher in the polar regions,  obtained relativistic jets along the rotation axis, which remained highly focused, and capable of penetrating the star. However, as these simulations were performed with a Newtonian hydrodynamic code, appreciably superluminal speeds in the jet flow were obtained.

We have performed axisymmetric relativistic simulations of jets from collapsars starting from Model 14A of . The simulations have been performed with GENESIS a multidimensional relativistic hydrodynamic code (based on Godunov-type schemes) developed by  using 2D spherical coordinates ($`r,\theta `$). GENESIS employs a 3th order explicit Runge–Kutta method  to advance in the time the relativistic Euler equations written in conservation form. High spatial order is provided by a PPM reconstruction  that sets up the values of the physical variables in order to solve linearized Riemann problems at every cell interface (using the Marquina’s flux formula ).

The innermost $`2.03M_{}`$ representing the iron core were removed from the helium star model by introducing an inner boundary at a radius of $`200`$km. When the central BH has acquired a mass of $`3.762M_{}`$, we mapped the model to our computational grid. In the $`r`$–direction the computational grid consists of 200 zones spaced logarithmically between the inner boundary and the surface of the helium star at $`R_{}=2.98\times 10^{10}`$cm. Assuming equatorial–plane symmetry we use four different zonings in the angular direction: 44, 90 and 180 uniform zones (i.e.,$`2^{},1^{}`$ and $`0.5^{}`$ angular resolution), and 100 nonuniform zones covering the polar region $`0^{}\theta 30^{}`$ with 60 equidistant zones ($`0.5^{}`$ resolution) and the remaining 40 zones being logarithmically distributed between $`30^{}\theta 90^{}`$.

The gravitational field of the BH is described by the static Schwarzschild metric, neglecting the effects due to self–gravity of the star. We used the EOS of  which includes the contribution of non–relativistic nucleons treated as a mixture of Boltzmann gases, and radiation, as well as an approximate correction due to pairs $`e^+e^{}`$. Full ionization and non-degeneracy of the electrons is assumed. We advect (i.e., we do not solve additional Riemann problems for each component) nine non-reacting nuclear species which are present in the initial model.

In a consistent collapsar model the jet will be launched by any physical process which gives rise to a local deposition of energy and/or momentum. We mimic this process by depositing energy at a constant rate, $`\dot{E}`$, within a $`30^{}`$ cone around the rotation axis of the progenitor star. In radial direction the deposition region extends from the inner boundary to a radius of $`6\times 10^7`$cm. We consider two cases that bracket the expected $`\dot{E}`$ of the collapsar models: $`10^{50}`$erg/s, and $`10^{51}`$erg/s.

## 2 Results

Low energy deposition rate (Model A). Using a constant $`\dot{E}=10^{50}`$erg/s a relativistic jet forms within a fraction of a second and starts to propagate along the rotation axis (Fig. 1). The jet exhibits all the typical morphological elements : a terminal bow shock, a narrow cocoon, a contact discontinuity separating ste-

llar and jet matter, and a hot spot. The propagation of the jet is unsteady, because of density inhomogeneities in the star. The Lorentz factor of the jet, $`W`$, increases non–monotonically with time, while the density drops to $`10^6`$ gr/cm<sup>3</sup>. The density profile shows large variations (up to a factor of 100) due to internal shocks. The mean density in the jet is $`10^21`$ g/cm<sup>3</sup>.

Some of the internal shocks are biconical and recollimate the beam. These shocks develop during the jet’s propagation and may provide the “internal shocks” proposed to explain the observed gamma–ray emission . A particularly strong recollimation shock forms during the early stages of the evolution, followed by a strong rarefaction that causes the largest acceleration of the beam material giving rise to a maximum in $`W`$. When the jet encounters a region along the axis where the density gradient is positive the jet’s head is decelerated, while the a central channel in the beam is cleaned by outflow into the cocoon through the head. This leads to an acceleration of the beam. The combination of both effects (deceleration of the head and beam acceleration) increases the strength of the internal shocks.

The relativistic treatment of the hydrodynamics leads to an overall qualitatively similar evolution than in  (formation of a jet), being, however, a quantitatively very different. We find that the results strongly depend on the angular resolution, and the minimum acceptable one is $`0.5^{}`$ (at least near the axis). At this resolution we find $`W_{\mathrm{max}}1520`$ (at shock break–out) at a radius $`8\times 10^9`$ cm. Within the uncertainties of the jet mass determinations due to finite zoning and the lack of a precise numerical criterium to identify jet matter, the baryon load, $`\eta `$, seems to decrease with increasing resolution. In the highest resolution run we find $`\eta 1.3\pm 1.2`$ at shock break-out (see also Sect. 4).

High energy deposition rate (Model B). Enhancing $`\dot{E}`$ by a tenfold ($`\dot{E}=10^{51}`$ erg/s), the jet flow reaches larger values of $`W_{\mathrm{max}}`$. We observe transients during which $`W_{\mathrm{max}}`$ becomes as large as 40 ($`W_{\mathrm{max}}=33.3`$ at shock breakout). The jet propagates faster than in model A. The time required to reach the surface of the star is 2.27 s instead of 3.35 s. The opening angle of the jet at shock breakout is $`10^{}`$, i.e., the jets is less collimated than model A. The strong recollimation shock present in the model A is not so evident here. Instead, several biconical shocks are observed, and $`W`$ 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 in the jet is 5 times smaller than in model A ($`\eta `$ being twice as large).

Evolution after shock breakout. After reaching the stellar surface the relativistic jet propagates through a medium of decreasing density continuously releasing energy into a medium whose pressure is negligible compared to that in the jet cavity, and whose density is (initially) of the same order as that of the jet. These are jump conditions that generate a strong blast wave. The external density gradient determines whether the shock will accelerate or decelerate with time (). In order to satisfy the conditions for accelerating shocks (), we have generated a Gaussian atmosphere matching an external uniform medium. We use models A and B to simulate the evolution after shock breakout. The computational domain is extended for this purpose to a radius of $`R_t=7.6\times 10^{10}`$ cm. The jet reaches $`R_t`$ (from the stellar surface) after 1.8 s in both models, i.e., the mean propagation velocity is $`0.85c`$ (almost three times larger than that inside the star).

The evolution after shock breakout can be distinguished into three epochs (see Figs. 1 and 2), which are related with (i) the external thermodynamical gradients and (ii) the importance of the axial momentum flux relative to the pressure into the jet cavity. Both effects determine the shape of the expanding bubble –prolate– (see Figs. 1 and 2) during the post–breakout evolution. However, when the jet reaches the uniform part of the circumstellar environment, the shape changes appreciably, because the sideways expansion is faster. We have not followed the evolution long enough to see what happens when most of the bubble has reached the uniform part of the environment. Nevertheless, we can infer from Fig. 2 than the widening rate reduces with time in a way similar to what has happened to the axial expansion. At latter times most of the bubble is inside the uniform medium, and the bubble will eventually be pressure driven. Hence a isotropic expansion is expected.

After shock breakout there are transients in which $`W_{\mathrm{max}}`$ becomes almost 50 in some parts of the beam, $`W_{\mathrm{max}}`$ is again obtained behind the strongest recollimation shock. The Lorentz factor near the boundary of the cavity blown by the jet grows from $`1`$ (at shock breakout) to $`3`$ in both models decreasing with latitude. At the end of the simulation $`W_{\mathrm{max}}`$ is 29.35 (44.17) for model A (B), which is still smaller than the ones required for the fireball model (). However, our simulations have not been pushed far enough in time yet and, therfore, they can (at the present stage) neither account for the observational properties of GRBs nor of their afterglows. Instead, our set of numerical models can be regarded as simulations of a proto-GRB, because the scales treated in the simulations are still by more than 100 times smaller than the typical distances at which the fireball eventually becomes optically thin ($`10^{13}`$ cm).