Neutrino Pair Annihilation in Collapsars

# Neutrino Pair Annihilation in Collapsars; Ray-Tracing Method in Special Relativity

Seiji Harikae11affiliationmark: 22affiliationmark: , Kei Kotake 22affiliationmark: 33affiliationmark: , and Tomoya Takiwaki22affiliationmark: 11affiliationmark: Department of Astronomy, The Graduate School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,Tokyo, 113-0033, Japan 22affiliationmark: Division of Theoretical Astronomy, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan 33affiliationmark: Center for Computational Astrophysics, National Astronomical Observatory of Japan, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan
###### Abstract

We develop a numerical scheme and code for estimating the energy and momentum transfer via neutrino pair annihilation (), bearing in mind the application to the collapsar models of gamma-ray bursts (GRBs). To calculate the neutrino flux illuminated from the accretion disk, we perform a ray-tracing calculation in the framework of special relativity. The numerical accuracy of the developed code is certificated by several tests, in which we show comparisons with the corresponding analytical solutions. Using hydrodynamical data in our collapsar simulation, we estimate the annihilation rates in a post-processing manner. We show that the neutrino energy deposition and momentum transfers are strongest near the inner edge of the accretion disk. The beaming effects of special relativity are found to change the annihilation rates by several factors in the polar funnel region. After the accretion disk settles into a stationary state (typically later than s from the onset of gravitational collapse), we find that the neutrino-heating timescale in the vicinity of the polar funnel ( km) can become shorter than the hydrodynamical timescale, indicating that the neutrino-heated outflows can be launched there. We point out that the momentum transfer can play as important role as the energy deposition for the efficient acceleration of neutrino-driven outflows. Our results suggest that the neutrino pair annihilation has a potential importance equal to the conventional magnetohydrodynamic mechanism for igniting the GRB fireballs.

accretion, accretion disks — gamma-ray burst: general — methods: numerical — magnetohydrodynamics — neutrinos — supernovae: general
slugcomment: Accepted to the ApJ

## 1 Introduction

Accumulating observational evidences have been reported so far, identifying a massive stellar collapse as the origins of long-duration gamma-ray bursts (GRBs), such as the preponderance of short-lived massive star formation for their host galaxies, as well as the identification of SN Ib/c light curves in the afterglows (Paczynski, 1998; Galama et al., 1998; Hjorth et al., 2003; Stanek et al., 2003). The duration of the long bursts may correspond to the accretion of debris falling into the central black hole (BH)(Piro et al., 1998). Pushed by those observations, the so-called collapsar has received quite some interest for their central engines (Woosley, 1993; Paczynski, 1998; MacFadyen & Woosley, 1999).

In the collapsar model, the central cores with significant angular momentum collapse into a BH. Neutrinos emitted from the accretion disk heat the matter of the polar funnel region to launch the GRB outflows. Paczynski (1990); Meszaros & Rees (1992) pioneerlingly proposed that the energy deposition proceeds predominantly via neutrino and antineutrino annihilation into electron and positron (e.g., , hereafter “neutrino pair annihilation”). The relativistic outflows are expected to ultimately form a fireball, which is good for explaining the observed afterglow (e.g., Piran (1999)). In addition, it is suggested that the strong magnetic fields in the cores of order of play also an active role both for driving the magneto-driven jets and for extracting a significant amount of energy from the central engine (e.g., Wheeler et al. (2000); Thompson et al. (2004); Uzdensky & MacFadyen (2007) and see references therein).

However, it is still controversial whether the generation of the relativistic outflows proceeds predominantly via magnetohydrodynamic (MHD) or neutrino-heating processes. So far, much attention has been paid to the MHD processes (e.g., Proga (2003); Proga et al. (2003); Mizuno et al. (2006); Fujimoto et al. (2006, 2007); Nagataki et al. (2007); Nagataki (2009); Harikae et al. (2009)). A general outcome of these extensive MHD simulations is that the magneto-driven shock waves, produced mainly by the field-wrapping of the strong precollapse magnetic fields or by the field amplification due to magnetorotational instability (MRI), can blow up massive stars along the rotational axis. Although those primary jet-like explosions are at most mildly relativistic due to too much baryons, they may be related to the X-ray flashes, which is a low energy analogue of the GRBs (Soderberg et al., 2006; Ghisellini et al., 2007), and those jets, which make a low density funnel along the rotational axis, could provide a good birthplace for the subsequent relativistic outflows (Harikae et al., 2009).

In contrast to such blossoms in MHD studies, there are only a few studies pursuing the possibility of generating jets by the energy deposition via neutrino pair annihilation. This is mainly because the neutrino emission from the accretion disk generally becomes highly aspherical, thus demanding one in principle, to solve the multidimensional neutrino transfer problem. For the first time in the numerical simulations, MacFadyen & Woosley (1999) pointed out the importance of the energy deposition via neutrino pair annihilation in their simulations, however the energy deposition rates to the polar funnel region were adjusted by hand to produce jets. Without those artificial energy injections, neutrino-driven outflows have been yet to be realized so far in the numerical simulations, to our best knowledge. It should be noted in the MHD-driven collapsars that too much precollapse magnetic fields with rapid rotation can lead to the MHD-driven explosion, which can prohibit the formation of a BH (e.g., Dessart et al. (2008); Bucciantini et al. (2009); Burrows et al. (2007)), thus not good for the collapsar scenario. Therefore it is still important to address the possibility of the neutrino-driven mechanism in collapsars.

Thus far, several methods for the implementation of the neutrino pair annihilation, have been reported. One straightforward method is to solve the full radiation field by the Monte-Carlo method (Tubbs, 1978; Janka & Hillebrandt, 1989). Still at present, this seems too expensive to combine with the hydrodynamic simulations. By estimating the local fluxes and spectra of the neutrino emission via the so-called neutrino leakage scheme, Ruffert et al. (1997) proposed to estimate the heating rate by summing up the local contributions of the neutrino and antineutrino radiation incident from all directions. However this scheme requires the double angular integration for neutrino and antineutrino for all the grid points in three dimension, which is still highly time-consuming and prevents its application to hydrodynamical simulation. Ruffert & Janka (1998) have upgraded this scheme to be more feasible, by defining the position of the neutrinosphere. Instead of the three-dimensional stellar volume as neutrino and antineutrino sources like in Ruffert et al. (1997), the energy deposition rates can be estimated by summing only over the two-dimensional neutrinospheres, which reduces the computational cost significantly.

Along this prescription, Nagataki et al. (2007) have estimated the neutrino heating rates, and included them to the hydrodynamical simulation. For reducing the computational time, they added one more assumption of the optically thinness of the accretion disk to the prescription by Ruffert & Janka (1998). Even with this potential overestimation of the heating rates, the neutrino-driven outflows were not observed in their simulations. As for the methodology of Ruffert & Janka (1998), the neutrino fluxes are estimated by summing up the local rates only along the direction (:perpendicular to the equatorial plane), which simplifies the directional dependence of the incoming neutrinos. The directional dependence can be more appropriately handled by doing the so-called ray-tracing analysis. Moreover, even when the spacial structure of the accretion disk becomes highly inhomogeneous, which is the case for our collapsar simulations, the ray-tracing analysis is good for omitting the surface of the accretion disk, from which the neutrino rays cannot travel, such as from the regions opposite to the central black hole.

More recently Dessart et al. (2009) have developed a new scheme to estimate the energy deposition rates via neutrino pair annihilation using the mutliangle neutrino-transport solver (Ott et al., 2008a). They discussed the possibility of the neutrino-driven outflow production in the postmerger phase of binary neutron-star coalescence. Without the general relativistic effects, this is the state-of-the-art method to estimate the neutrino pair annihilation. As for the general relativistic effects, Birkl et al. (2007) have performed a series of ray-tracing calculations in a Kerr spacetime. They found that the total energy deposition rates measured by an observer at infinity can increase by a factor of two by the general relativistic effect, which is consistent with the result of previous studies (e.g., Asano & Fukuyama (2001, 2000)).

Joining in the extensive efforts mentioned above, here we present a numerical scheme and code for calculating the neutrino heating via pair neutrino annihilation, which is designed to be incorporable to the hydrodynamical simulations of collapsars. Relying on the neutrino leakage scheme as in Ruffert et al. (1997); Ruffert & Janka (1998), we estimate the incident neutrino flux by doing a ray-tracing calculation from the neutrino sphere, on which the neutrino distribution function is assumed to take a Fermi-Dirac distribution. The ray-tracing calculation is done in the Minkowskian spacetime (e.g., Kotake et al. (2009a)), neglecting the general relativistic (GR) effects (e.g., Birkl et al. (2007)) for simplicity. Thus our scheme is limited to capture the special relativistic effects such as beaming effects. One of the main characteristics of our method is that the neutrino pair annihilation can be estimated only with quantities defined on the numerical grids of the hydrodynamical simulation. This can be done by carrying out a Lorentz transformation of radiation variables in the comoving frame back to the Eulerian frame. By doing so, the computational cost of the double angular integration that is required in estimating the annihilation rates, can be reduced significantly. We derive the formulation required for this, and describe it elaborately. We check the numerical accuracy of the developed code by showing the comparisons with analytical solutions, some of which we newly give in this paper. Based on the results of our long-term collapsar simulation (Harikae et al., 2009), we run our new code to estimate the heating rate in a post-processing manner. Although the presented scheme is not the state-of-the-art in comparison with Birkl et al. (2007); Dessart et al. (2009), we elaborately study the possibility of the neutrino-driven mechanism in collapsars and will point out that the neutrino pair annihilation has a potential importance equal to the conventional MHD mechanism for igniting the GRB fireballs.

This paper is organized as follows. We summarize the special relativistic formulation of the neutrino energy deposition and momentum transfer in Section 2. Section 3 is devoted to the main results. The numerical tests and its comparison with the analytical solutions are shown in section 3.1. In section 3.2, we estimate the annihilation rates in a post-processing manner using hydrodynamical data in our collapsar simulation. We summarize our results and discuss their implications in Section 4.

## 2 Neutrino pair annihilation in special relativity

Before going into details, we first show an illustrative figure, which presents our calculation concept for the neutrino pair annihilation in collapsars (Figure 1). For estimating the neutrino pair annihilation at a given target, we trace each neutrino trajectory (white lines) backwards till it hits to the surface of the accretion disk (colored by red). It should be noted that the surface closely coincides with the surface of the neutrino sphere (footpoints of the rays on the tori), from which neutrinos and antineutrinos emerge freely out. As indicated in this figure, some fraction of neutrinos comes out from the the inner edge of the disk, whose rotational speed is close to the speed of light. This suggests that the special relativistic effects become important towards the precise determination of the heating rates.

For later convenience, we define the two frames namely of and in the spherical coordinate systems (Figure 2). Here is the coordinate system that we use for describing the variables obtained in the hydrodynamical calculations. On the other hand, is the one for describing the target points where the pair annihilation occurs. Note that the coordinates between and system, differ only in the position of their basis. in the figure represents the position of neutrino/antineutrino source. is usually described by the system. In the following, variables measured in the system, are described with the subscript or , indicating that they are related to the neutrino sources. With these definitions, we move on to describe the special relativistic formulation from the next section.

### 2.1 Formulation of energy deposition and momentum transfer rates

The energy deposition rate via neutrino pair annihilation (e.g., Goodman et al. (1987); Asano & Fukuyama (2000)) is written as,

 dq+ν¯ν(\boldmathr)dtdV = ∬fν(pν,\boldmathr)f¯ν(p¯ν,\boldmathr)σ|\boldmathvν−\boldmathv¯ν|(ϵν+ϵ¯ν)d3\boldmathpνd3\boldmathp¯ν, (1)

where is the number density of neutrinos in phase space, , , , is the momentum, energy, and velocity of neutrino, respectively. Those definitions are the same for antineutrino by changing the notation to . is the cross section given by

 σ=2c2KG2F(\boldmathpν⋅\boldmathp¯ν), (2)

where the dimensionless parameter is written as

 K(νe,¯νe)=1+4sin2θW+8sin4θW6π, (3) K(νμ,¯νμ)=K(ντ,¯ντ)=1−4sin2θW+8sin4θW6π. (4)

Here, the Fermi constant , and the Weinberg angle is . Note in our setting of the coordinate system (e.g., Figure 2) that corresponds to the origin of the (targets) system measured in the system. Variables with subscripts of () are quantities upon the neutrinosphere, which is labeled by in Figure 2. They are also measured in the system.

Aiming at a straightforward incorporation of the ray-tracing calculation into the hydrodynamical simulations, we calculate the energy deposition(momentum transfer) rates in the laboratory (lab) frame. On the other hand, the quantities related to radiation, such as the distribution function and emissivity/absorptivity are locally defined in the rest frame of fluid, in which the radiation isotropy is maintained. Designating the subscript 0 to the variables in the rest frame, the Lorentz transformations between the two frames are given in Mihalas & Mihalas (1984) as,

 dtdV = dt0dV0 (5) dϵ = ϵϵ0dϵ0 (6) dΩ = ϵ20ϵ2dΩ0 (7) ϵ = γ(1+βμ0)ϵ0 (8) = ϵ0/[γ(1−βμ)],

where , and with and , being the velocity of fluid and the unit vector describing a given ray of neutrino, within a solid angle of .

Substituting from Equations (5) to (8) in Equation (1) gives the heating rate in the lab frame as,

 dq+ν¯ν(\boldmathr)dtdV = 2cKG2F∫dθνdϕνdθ¯νdϕ¯ν (9) ×[ξ5ν(\boldmathr,Ων)ξ4¯ν(\boldmathr,Ω¯ν)Eν,0(\boldmathr,Ων)N¯ν,0(\boldmathr,Ω¯ν) +ξ4ν(\boldmathr,Ων)ξ5¯ν(\boldmathr,Ω¯ν)Nν,0(\boldmathr,Ων)E¯ν,0(\boldmathr,Ω¯ν)] ×[1−sinθνsinθ¯νcos(φν−φ¯ν)−cosθνcosθ¯ν]2sinθνsinθ¯ν.

Here reflects the special relativistic correction to the neutrino energy as

 ξν(\boldmathr,Ων) = ϵν/ϵν,0 (10) = 1/[γν(1−μνβν)].

Here it should be noted again that the variables with the subscript 0 and denote the ones defined in the rest frame of the neutrino source and on the neutrino sphere respectively, such that

 βν=|\boldmathVν|c, (11)

is the velocity of the fluid on the neutrino sphere with its Lorentz factor of

 γν=1√1−β2ν, (12)

and the directional cosine of

 μν=\boldmathnν⋅\boldmathVν|\boldmathVν|, (13)

where is the direction of the rays of neutrino from the neutrino sphere,

 \boldmathnν = \boldmathpν|\boldmathpν|=(sinθνcosϕν, sinθνsinϕν, cosθν). (14)

The following two quantities are the energy-weighted integration of the neutrino distribution function on the neutrino sphere (namely )) as,

 Eν,0(\boldmathr,Ων) = ∫ϵ4ν,0fν,0(\boldmathrν,0,\boldmathpν,0)dϵν,0, (15) Nν,0(\boldmathr,Ων) = ∫ϵ3ν,0fν,0(\boldmathrν,0,\boldmathpν,0)dϵν,0,. (16)

As in Ruffert et al. (1997); Ruffert & Janka (1998); Nagataki et al. (2007), the neutrino number flux is simply assumed to be conserved as along its trajectory to the target. The neutrino distribution function on the neutrino sphere is assumed to take a Fermi-Dirac shape with a vanishing chemical potential as,

 fν(\boldmathrν,0,\boldmathpν,0) = 1(hc)3dn0dϵ0dΩ0dt0dV0 (17) = 1(hc)31exp(ϵν,0/kTν,0)+1,

where is taken to be the same with the matter temperature on the neutrino sphere of . Now the energy integrals in Equations (15) and (16) can be expressed by the Fermi integrals as

 Fk(y) ≡ ∞∫0xkexp(x−y)+1dx, (18) Eν,0(\boldmathr,Ω) = (kT(\boldmathrν))5(hc)3F4(0), (19) Nν,0(\boldmathr,Ω) = (kT(\boldmathrν))4(hc)3F3(0). (20)

The surface of the neutrino spheres is defined per each lateral angular grid in the hydrodynamic calculation () as

 ∫∞Rνi(θ)1λi dr=2/3, (21)

where is the radius of the neutrinosphere for the neutrino species of , namely for , and with the corresponding mean free path calculated by a multi-flavor leakage scheme (Epstein & Pethick (1981); Rosswog & Liebendörfer (2003); Kotake et al. (2003); Takiwaki et al. (2009)) with special relativistic corrections (Harikae et al., 2009). As for the neutrino opacities, electron capture on proton and free nuclei, positron capture on neutron, neutrino scattering with nucleon and nuclei, photo-pair, plasma processes are included (e.g., Rosswog & Liebendörfer (2003)).

Following the same procedure above, the momentum transfer rate in special relativity can be also derived as,

 d\boldmathm+ν¯ν(\boldmathr)dtdV = 2KG2F∫dθνdϕνdθ¯νdϕ¯ν (22) ×[ξ5ν(\boldmathr,Ων)ξ4¯ν(\boldmathr,Ω¯ν)Eν,0(\boldmathr,Ων)N¯ν,0(\boldmathr,Ω¯ν)% \boldmathnν +ξ5¯ν(\boldmathr,Ω¯ν)ξ4ν(\boldmathr,Ων)E¯ν,0(\boldmathr% ,Ω¯ν)Nν,0(\boldmathr,Ων)% \boldmathn¯ν] ×[1−sinθνsinθ¯νcos(φν−φ¯ν)−cosθνcosθ¯ν]2sinθνsinθ¯ν.

Finally the remaining procedure is to change the coordinates of the angle integrations in Equations (9) and (22) from the to coordinate system. Albeit sacrificing the formulae to be messy (e.g., appendix A), this procedure, by which the neutrino source and the target can be connected one by one, merits for reducing the numerical costs significantly. The coordinate transformation with respect to the solid angle is written as,

 dΩν = Jrμ(μν,ϕν,rs,ν,μs,ν)|ϕs,ν=ϕsph,νdrs,νdμs,ν (23) + Jμϕ(μν,ϕν,μs,ν,ϕs,ν)|rs,ν=rsph,νdμs,νdϕs,ν + Jϕr(μν,ϕν,ϕs,ν,rs,ν)|μs,ν=μsph,νdϕs,νdrs,ν,

where is the position of measured in the system. Using the expressions of Jacobian (, see appendix A for details), we can calculate the annihilation rates based only on the quantities defined on the hydrodynamical grid-points.

### 2.2 Implementation of ray-tracing calculation

We have shown how to calculate the annihilation rates when neutrinos from the neutrino sphere can reach directly to the target region. However in reality, neutrinos can be absorbed in some situations before hitting to the target, for example, by encountering with the optically thick region. To correctly count the real contribution, we set the following three criteria by utilizing the ray-tracing calculation.

The first task is to define an outward neutrino emission from the neutrino sphere. This can be satisfied by , where the former is the normalized neutrino momentum from the neutrino sphere , and the latter is the normal direction to the differential surface of the neutrino sphere. As mentioned, we define the surface of the neutrino sphere to be dependent only on the lateral direction (e.g., Equation (21)). Thus the analytic form of can be simply given using the valuables such as and that denote the surface position of the neutrino sphere. Furthermore the criterion of would provide more realistic regulation of the neutrino flux than in Ruffert et al. (1997); Ruffert & Janka (1998) whose criterion is determined only by the density gradient at the neutrino sphere.

The second criterion is the positivity of () along each trajectory from the surface of the neutrino sphere to the target region. Here denotes the (locally-defined) opacity for each neutrino species (e.g., the integrand of the left-hand-side of Equation (21)). This is necessary to omit the neutrino flux which comes from the optically thick region, such as inside the neutrino sphere. For not counting radiation absorbing to the black hole (BH), the third criterion is the positivity of , where with being the gravitational constant and the mass of the BH, respectively. For simplicity, we will set to be the Schwartzshild radius when we apply the ray-tracing calculation to the collapsar simulations in the later section. However in reality, the BH is rotating so that the accurate estimation of the geodesics is necessary, which can be an another possible extension of this study.

## 3 Result

### 3.1 Numerical Tests

Before applying the newly developed code to collapsars, we shall check the accuracy of our code. For the purpose, we derive some analytical solutions for test problems and present comparison with the numerical results in this section.

Amongst a list of tests, the first requirement we set is to check whether the newly derived Jacobians (Equation (23)) and their numerical implementations are correct. This will be checked by reproducing radiation fields shedding from a light-bulb in spherical and non-spherical geometry (section 3.1.1 and 3.1.2). The second one is to check the accuracy of our ray-tracing calculation. Note that this can be checked in each test, because the ray-tracing is done in whichever tests. In addition, we present an additional test, which is specially designed to test the validity for the application to the collapsar simulations in section 3.1.3.

In the following tests, our numerical domain has the spherical coordinates with equally spaced 100 radial meshes covering from 0 to 500 km in radius. To check the numerical convergence, we vary the numerical resolution for and direction, respectively. For simplicity, the axial and equatorial symmetry is assumed for the hydrodynamic quantities (such as for density, electron fraction and so on) as a background. The density and temperature on the neutrino sphere is set to be and , respectively.

To measure the deviation from the corresponding analytical solutions, we estimate the error function of , which is defined by the norm as,

 E ≡ L1,errL1,ana, (24) L1,err ≡ Σ|xnum−xana|, (25) L1,ana ≡ Σxana, (26)

where represents the numerical values such as the energy deposition rate () and the absolute value of the momentum transfer rate (), while represents the corresponding analytical solutions. indicates the summation, which is taken over the whole computational domain.

As for the momentum transfer rate, the error function of cannot always be a useful measure because the analytical solution becomes zero in several situations (see sections 3.1.2 and 3.1.3). In such a case, we define another error function of as,

 E2 ≡ MnumMmax, (27) Mnum ≡ Σ|\boldmathm+ν¯ν|, (28) Mmax ≡ Σm+max, (29)

where is the analytical solution, which is defined from the energy deposition rate as , thus meaning the maximum momentum transfer. With this , we will evaluate how precisely the cancellation can be maintained.

#### 3.1.1 Radiation from Spherical Neutrino Sphere

As a simplest test, we first present how accurately we can reproduce the radiation field from a spherical source, which emits radiation isotropically from its surface. In this test, we set a spherical inner boundary of 50 km in radius, inside which material is taken to be completely optically thick for neutrinos, while outside is taken to be completely optically thin. Albeit very crude, such a situation is akin to the neutrino radiation from the neutrino sphere in the postbounce phase of core-collapse supernovae. In addition to this static source, we consider another case in which the spherical neutrino sphere rotates relativistically. By this test, we hope to see and test the special relativistic effects, which is one of the main focus in this paper.

\subsubsubsection

Non Relativistic Case

In the case of the static neutrino sphere, the analytic formulae of the heating and momentum transfer rate have been explicitly given in Goodman et al. (1987); Cooperstein et al. (1987); Salmonson & Wilson (1999); Asano & Fukuyama (2000), as

 dq+ν¯ν(\boldmathr)dtdV=4cKG2F3(kT)9(hc)6F4(0)F3(0)F(r), (30)

and

 dm+ν¯ν(\boldmathr)dtdV=4KG2F3(kT)9(hc)6F4(0)F3(0)G(r)\boldmathnr, (31)

where geometrical factor and is given by

 F(r)=2π23(1−X)4(5+4X+X2), (32)

and

 G(r)=π26(1−X)4(1+X)(8+9X+3X2), (33)

with

 X=√1−(Rνr), (34)

where is the radius of the neutrino sphere.

Figure 3 shows the energy (left panel) and momentum (right panel) rates numerically calculated with our code (left-half) and the corresponding analytical solutions from Equation (30) (right-half). Note that the label such as ”20X16” at the top right edge indicates the resolution of our ray-tracing calculation, here , the equally spaced mesh numbers in the lateral () and azimuthal direction (), respectively. Comparing the two panel, we find no visible differences in each panel. In Figure 4, the deviation from the analytical solution for the different numerical resolutions (e.g., Equation (24)) is shown. As shown, the relative error decreases with as low as 0.1 %. Compared to other studies treating the radiation problems (e.g., Stone et al. (1992); Turner & Stone (2001)), the obtained accuracy here seems high enough. These results assure the validity both of the newly derived expression of (from Equation (A1) to (A16)) and its implementation.

\subsubsubsection

Relativistic Case

For the neutrino sources with relativistic motion, it is not generally possible to write down the special relativistic factor of (Equation (10)) in an analytic form. This is because it depends locally on the angle between the velocity field and the flight direction of neutrinos. To see special relativistic effects, we examine a special case, in which the analytic solution can be found. The annihilation rates along the polar axis, emitted from the relativistically rotating neutrino source, can be analytically given as,

 dq+ν¯ν(\boldmathr)dtdV = 4cKG2F3(kT)9γ9(hc)6F4(0)F3(0)F(r), (35) dm+ν¯ν(\boldmathr)dtdV = 4KG2F3(kT)9γ9(hc)6F4(0)F3(0)G(r)\boldmathnr, (36)

where the geometrical factors and are just same as Equations (32) and (33), respectively. The sharp suppression with a factor of is the outcome of the relativistic beaming, simply because the rotational velocity is perpendicular to the polar direction. We will see it soon in the following.

Figure 5 shows the energy deposition (left) and momentum transfer rates (right) from the relativistically rotating neutrino sphere. Here we assume with vanishing and . Note in this figure that the numerical results are only shown, owing to the inaccessibility of the analytical solution as mentioned above. It is shown that the annihilation rates both for energy and momentum is greatly reduced near on the rotational axis as expected from Equations (35) and (36), while they are enhanced around the equatorial plane. Figure 6 shows the comparison of the numerical solution (solid line) with the analytical solution (dashed line) along the rotational axis for energy (left) and momentum (right) transfer rates. In support of the validity of our special relativistic implementation, no visible differences are seen between the analytic and numerical solutions. Figure 7 shows the deviation of (e.g., Equation (24)) for the energy (left) and momentum (right) from the analytic solution. From this test, it is shown that the lateral grid points () should be at least more than 20 to maintain the percent levels of agreement with the analytic solution.

#### 3.1.2 Spherical Cavity

In this test, we set a spherical cavity of km in radius, outside which is the uniform opaque neutrino source. Although this is not a realistic astrophysical situation, it provides one of the most simplest null test of the momentum transfer because the neutrino flux is isotropic everywhere inside the cavity. For the heating rate, this test merits because the geometrical factor (Equation (30)) is analytically given by

 F(r)=64π23. (37)

Comparing the numerical (left-half) and analytical solution (right-half) in each panel of Figure 8, the discrepancy is shown to become smaller for the better numerical resolution (from left to right). From Figure 9 showing the and errors for the energy (left) and momentum transfer rate (right), it can be seen that the errors mainly depend on . By setting , we can get the agreement with an accuracy of 1 % for and 0.1 % for , which would satisfy the requirement of the vanishing momentum transfer sufficiently.

#### 3.1.3 Non-Spherical Cavity

We now move on to the non-spherical cavity test, by which we can check if the radiation contribution from the face element of is accurately estimated. For simplicity, we just modify the shape of the spherical cavity in the previous section. We set the position of the neutrino sphere to change with the polar angle . For and , the neutrino sphere is situated at 300 km, while for , at 100 km. Region inside the neutrino sphere is taken to be completely thin, while completely thick outside. Even with the geometrical difference, it should be noted that the items to be checked are the same for the case of the spherical cavity, because the net flux is isotropic and does not change.

From Figure 10, it can be seen that some inhomogenities seen closely at the edge of the boundary (left-half: numerical) become smaller for the higher resolutions (to the right panel). From Figure 11, the decreasing rate of the relative errors for and is shown to be smaller than the ones from the previous tests due to the relatively complicated geometry. By taking the resolution of (, , our code can reproduce the analytic solutions with percent levels of accuracy here.

As evident from the numerical convergence with increasing the numerical resolution (e.g., see Figures 4, 7, 9, and 11), the previous tests have supported the validity of the newly derived expression of the Jacobians (from eq.(A1) to (A16)), its implementation, and the ray-tracing calculation. From those test calculations, it is suggested that the grid numbers of are required for obtaining the percent levels of agreement with the analytical solution. Considering the computational cost, we choose to employ the grid numbers of in the actual implementation for collapsar simulations in the following section.

### 3.2 Application to Collapsar

Having demonstrated the accuracy of our code in the previous sections, we now proceed to show an application to collapsars. In this section, we estimate the annihilation rates in a post-processing manner based on our simulation results of the long-term evolution of collapsar and then discuss the possible impacts on the collapsar dynamics.

#### 3.2.1 Recipes for Collapsar Simulation

For obtaining the hydrodynamic profiles for the ray-tracing calculation, we perform the special relativistic simulations of collapsar (Harikae et al., 2009). Without repeating the numerical procedures again, we will shortly summarize major modifications added in this work, namely about the initial model and the position of the inner-boundary.

Since our focus is not on the MHD-driven collapsars, we construct initial models with no magnetic fields. Taking the precollapse density/electron-fraction/entropy profiles of progenitor of 35OC in Woosley & Heger (2006), we embed the angular velocity analytically as

 Ω(r,θ)=Ω0X40+αΩlso(M(X))X4X40+X4, (38)

where , . is the mass coordinate at , and and is the specific angular momentum and angular velocity in the last stable orbit. Model parameters are and . This distribution provides co-rotation with at connected to constant specific angular momentum at . would correspond to the size of co-rotating region in pre-collapse phase, such as Fe core. To be closer to the original rotational profile in Woosley & Heger (2006), we set , , and , where is the size of Fe core ( for 35OC).

Another major modification is the position of the absorbing inner boundary of the computational domain. For reducing the computational cost, the radius of the boundary, which mimics the surface of black hole, was set to be km (Harikae et al., 2009). This manipulation may lead to the underestimation of the neutrino luminosity by excising the radiation from more inner edge of the accretion disk. This is obviously disadvantageous for producing the neutrino-driven outflow. In this paper, we take more compact inner-boundary as , extending to the outer boundary of 30000 km.

#### 3.2.2 When can the pair neutrino annihilation work ?

For energizing the jets via pair neutrino annihilation, the higher neutrino luminosity from the accretion disk and the lower density along the polar region, are favorable. Our collapsar model starts with a rapid mass accretion into the central objects. Within 2.0 s after the onset of gravitational collapse, this model accretes more than 2 in the center. This indicates the black hole formation in the center because the maximum mass of the neutron star, which the Shen equation of state employed here can sustain, is in the same mass range (Kiuchi & Kotake, 2008). Simultaneous with the rapid infall, the neutrino luminosities also show a drastic increase, and then shift to a gradual increase reflecting the mass accretion to the newly formed accretion disk. At 8 s from the onset of gravitational collapse, the total neutrino luminosities gets as high as , which consists of for , for , for and . Afterwards, the total luminosity keeps almost constant with time, reflecting that the accretion disk comes to a (quasi)stationary state.

Figure 12 shows the distributions of density, temperature, electron fraction, and Lorentz factor of our model at 9.1 s after the onset of gravitational collapse. From the left panel, the structure of the accretion disk with a polar funnel along the rotation axis is shown. The disk temperature reaches typically MeV, which provides the high neutrino luminosity of . From the right panel, the Lorentz factor at the inner edge of the disk becomes , indicating the importance of the special relativistic effects. Taking these hydrodynamical data, we estimate the neutrino pair annihilation by the ray-tracing calculation in a post-processing manner.

#### 3.2.3 Energy deposition and momentum transfer rates in collapsar

The left- and right- hand sides of the left panel of Figure 13 shows the energy deposition and the momentum transfer rates (its absolute value). The energy and momentum transport from the accretion disk (regions colored by black) to the polar funnel regions (regions colored by red) can be apparently seen. Interestingly, the transfer rates become highest near the inner edge of the accretion disk. Inside 80 km in radius from the center, the typical energy deposition rates along the polar axis are , where the momentum transfer rates are the order of . Thus the momentum transfer rates multiplied by the speed of light becomes comparable to the energy deposition rate there.

For our model, the density along the polar funnel regions is as low as (see Figure 12(left)). Assuming that the momentum transfer there is maintained for 10 ms with conserving its momentum, we can estimate that the material along the funnel can be accelerated up to the Lorentz factor of . This suggests that the momentum transfer plays as important role as the energy deposition for the efficient acceleration of the neutrino-driven outflows. The right panel of Figure 13 is the annihilation rates without the special relativistic correction ( in Equation (10)). Comparing the two panels in Figure 13, we find that the special relativistic effect decreases the transfer rates by several factors inside the polar funnel. As already discussed in section 3.1.1, this is the outcome of the relativistic beaming by the rapidly rotating accretion disk. On the other hand, the reduction in the total deposition rates, given by the volume integral of the local rates in Figure 13, can stay relatively smaller, namely of and with and without the correction, respectively. These values are comparable to what expected by Nagataki et al. (2003), since the mass accretion rate in this epoch is about . The resulting conversion efficiency of the heating, which is defined by the ratio of the total deposition rates to the total neutrino luminosity, is 0.514% and 0.630 % with and without the correction. Those numbers are comparable to the ones in the previous studies (e.g., Ruffert et al. (1997); Ruffert & Janka (1998, 1999); Narayan et al. (2001); Setiawan et al. (2004, 2006)), and this negative effect of special relativity on the energy deposition rate is also mentioned by Birkl et al. (2007).

#### 3.2.4 Comparison of timescales

Based on the annihilation rates in the last section, we compare various timescales in this section, such as the neutrino-heating timescale and the dynamical timescale. Then we anticipate if the neutrino-heating outflows could be produced or not in the polar funnel regions.

To trigger the neutrino-heating explosion, the neutrino-heating timescale should be smaller than the advection timescale, which is characterised by the free-fall timescale in the polar funnel regions. This condition is akin to the condition of the successful neutrino-driven explosion in the case of core-collapse supernovae (e.g., Bethe (1990) and see collective references in Janka et al. (2007)). The heating timescale is the timescale for a fluid to absorb the comparable energy to its gravitational binding energy to be gravitationally unbound, which is defined as , where is the average density at a certain radius and we take , is the local gravitational potential. The dynamical timescale is defined as . Figure 14 depicts the ratio of the dynamical to the heating timescales , showing that the ratio becomes greater than unity inside 80 km in the vicinity of the rotational axis. This indicates the possible formation of the neutrino-driven outflows, if coupled to the collapsar’s hydrodynamics.

Figure 15 shows various timescales for the material along the rotational axis in Figure 14. characterizes the timescale for matter to become relativistic by the neutrino heating. As mentioned, the first criterion of , is satisfied for the regions inside 80 km (compare green and blue lines). More interestingly, inside 50 km, gets slightly shorter than (compare red and blue lines), indicating that matter could attain relativistic motions there. To see these outcomes needs the implementation of the ray-tracing calculation in the hydrodynamic calculation, which we pose as a next task of this study. It is noted that with our choice of the grid numbers (), it generally takes several minutes (in the CPU time) for each ray-tracing calculation, using 256 processors in the Cray XT4 system at the National Astronomical Observatory in Japan. To follow the dynamics typically later than 10 s after the onset of core-collapse for satisfying the condition, 10000000 hydro-steps are generally needed because the hydrodynamical timescale is an order of s in our simulation. Therefore it is still computationally expensive to turn on the ray-tracing calculation throughout the whole simulation. In the actual implementation, we plan to switch on the ray-tracing calculation by monitoring intermittently during the simulations in the computational domain. After the ray-tracing calculation is turned on, seeing , the neutrino-driven outflows will be launched within 0.1 s (Harikae et al. in preparation). Thus the computational time using the above supercomputer can be reduced to several months, which are expected to be quite feasible.

## 4 Summary and discussion

Bearing in mind the application to the collapsar models of gamma-ray bursts, we have developed a numerical scheme and code for calculating the energy and momentum transfer via neutrino pair annihilation in the framework of special relativity. To estimate the transfer rates, we perform a ray-tracing calculation of the neutrino flux from the accretion disk of collapsars. We have tested the numerical accuracy of the developed code by several numerical tests, in which the comparisons with the corresponding analytical solutions were shown. Using hydrodynamical data obtained in our collapsar simulation, we have run our code to estimate the annihilation rates in a post-processing manner. We have obtained the following results.

For the progenitor chosen in this study, the accretion disk settles into a stationary state typically s after the onset of gravitational collapse. The accretion disk later on can be luminous source of neutrinos as high as . At this epoch, the polar funnel along the rotation axis is formed, which is heated intensively from the disk via neutrino pair annihilation. The beaming effects of special relativity are found to decrease the transfer rates by several factors in the polar funnel regions. Inside 80 km in the vicinity of the polar funnel, we point out that the neutrino-heating timescale can become shorter than the hydrodynamical timescale. This indicates the possible formation of the neutrino-heated outflows there. Our results also suggest that the momentum transfer via neutrino pair annihilation plays as important role as the energy deposition for the efficient acceleration of the neutrino-driven outflows.

Finally we shall make comparison with previous work and also mention some imperfections of this study. As mentioned, general relativistic effects ignored in this study have been considered as an important ingredient of the neutrino pair annihilation in collapsars (Jaroszynski, 1993, 1996; Salmonson & Wilson, 1999; Asano & Fukuyama, 2001, 2000; Birkl et al., 2007). It is noted that we have calculated the energy and momentum deposition rates for the thin disk models in Birkl et al. (2007) to assess the validity of our code. Our code can reproduce their results well both for the Newtonian and special relativistic case (see Figure 16 for the 2D configurations). In fact, the total heating rates are and for the Newtonian (model DN) and the special relativistic case (model DN4L), which are in Birkl et al. (2007), and , respectively (see Table 1 in Birkl et al. (2007)). When we boldly extrapolate results in Birkl et al. (2007) to ours, the local heating rates can be enhanced significantly (by several factors) in the vicinity of the central BH due to the GR effects. To update the present scheme to the one in general relativity requires not only the ray-tracing in the curved space, but also the general relativistic formulation of radiative transport, which we are to investigate as extension of this study. At the same time, we should improve the ray-tracing calculation to take into account the charged-current absorption occurring along a given ray, which is ignored in this study. Without the neutrino absorption, the heating rates estimated in this paper should give a upper bound. We also plan to include the effects of magnetic fields on this study. By changing the strength of magnetic fields systematically, we hope to clearly understand how the outflow production in collapsars, depending on the precollapse rotation, change from the neutrino-driven mechanism to the MHD-driven one. The magnetic effects on the annihilation rates and its impacts on the dynamics are also remained to be studied (e.g., Zhang & Dai (2009)).

Furthermore, the neutrino oscillation by Mikheyev-Smirnov-Wolfenstein (MSW) effect (see collective references in Kotake et al. (2006); Kawagoe et al. (2009)) could be important, albeit in much later phase than we considered in this paper. It should be remembered that the typical density in the polar funnel was in the present model. When the density there drops as low as later, the neutrino oscillation could operate for neutrinos traveling from the accretion disk to the polar funnel. If this is the case, the incoming neutrino spectra to the polar funnel regions and the pair annihilation rates there could be affected significantly. Together with the effects of neutrino self-interaction (e.g., Duan et al. (2006)), there seem a lot of issues to be clarified about the neutrino oscillation in the collapsar environments. As in the case of core-collapse supernovae (e.g., Kotake et al. (2009b); Ott et al. (2008b) and see references therein), studies of gravitational-wave emissions from collapsars might provide us a new window to probe into the central engine (e.g., Hiramatsu et al. (2005); Suwa & Murase (2009)). As a sequel of this work, we are planning to implement the ray-tracing calculation in the hydrodynamic simulation and clarify those aspects. We just stand at a starting line to study those interesting topics, which will be presented elsewhere one by one soon.

S.H. is grateful to T. Kajino for fruitful discussions. T.T. and K.K. express thanks to K. Sato, S. Yamada, and S. Nagataki for continuing encouragements. Numerical computations were in part carried on XT4 and general common use computer system at the center for Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. This study was supported in part by the Grants-in-Aid for the Scientific Research from the Ministry of Education, Science and Culture of Japan (Nos. S19104006, 19540309 and 20740150).

## Appendix Appendix A Calculations of Jacobian

In this section, we give an exact expression of the Jacobians () namely,

 Jrμ(μν,ϕν,rs,ν,μs,ν) ≡ ∣∣ ∣ ∣∣∂μν∂rs,ν∂μν∂μs,ν∂ϕν∂rs,ν∂ϕν∂μs,ν∣∣ ∣ ∣∣, (A1) Jμϕ(μν,ϕν,μs,ν,ϕs,ν) ≡ ∣∣ ∣ ∣∣∂μν∂μs,ν∂μν∂ϕs,ν∂ϕν∂μs,ν∂ϕν∂ϕs,ν∣∣ ∣ ∣∣, (A2) Jϕr(μν,ϕν,ϕs,ν,rs,ν) ≡ ∣∣ ∣ ∣∣∂μν∂ϕs,ν∂μν∂rs,ν∂ϕν∂ϕs,ν∂ϕν∂rs,ν∣∣ ∣ ∣∣. (A3)

by which we can estimate the heating/momentum-transfer rates only by the quantities defined on the numerical grids of the hydrodynamical simulations.

From a geometrical calculation, albeit a little bit messy, the following partial derivatives can be found straightforward as,

 ∂μν∂rs,ν = 1rs,ν(1+\boldmathns,ν⋅\boldmathnr)(\boldmathnr)z+(|\boldmathn% r|2−1)cosθs,ν(1−2\boldmathns,ν⋅\boldmathnr+|\boldmathnr|2)32, (A4) ∂μν∂μs,ν = sinθs,ν(1−2\boldmathns,ν⋅\boldmathnr+|\boldmathnr|2)−(cosθs,ν−(\boldmathnr)z)\boldmathnr⋅\boldmathks,νsinθs,ν(1−2\boldmathns,ν⋅\boldmathnr+|\boldmathnr|2)32, (A5) ∂μν∂ϕs,ν = (cosθs,ν−(\boldmathnr)z)(%\boldmath$n$r)⊥z⋅\boldmathls,ν(1−2\boldmathns,ν⋅\boldmathnr+|\boldmathnr|2)32, (A6) ∂ϕν∂rs,ν = 1rs,ν(|(\boldmathnr)⊥z|2−\boldmathns,ν⋅(\boldmathnr)⊥z)(\boldmathns,ν)x+(sin2θs,ν−%\boldmath$n$s,ν⋅(\boldmathnr)⊥z)(%\boldmath$n$r)x((\boldmathnr)y−(% \boldmathns,ν)y)(sin2θs,ν−2% \boldmathns,ν⋅(\boldmathnr)⊥z+|(% \boldmathnr)⊥z|2), (A7) ∂ϕν∂μs,ν = −1sinθs,ν(\boldmathnr)⊥z⋅\boldmathms,νsin2θs,ν−2% \boldmathns,ν⋅(\boldmathnr)⊥z+|(% \boldmathnr)⊥z|2, (A8) ∂ϕν∂ϕs,ν = sin2θs,ν−\boldmathns,ν⋅(\boldmathnr)⊥zsin2θs,ν−2\boldmathns,ν⋅(\boldmathnr)⊥z+|(\boldmathnr)⊥z|2, (A9)

where the following valuables are dependent only on the position of the neutrino sphere () and the direction of a specified direction () as,

 \boldmathns,ν = \boldmathrs,ν|\boldmathrs,ν|=(sinθs,νcosϕs,ν,sinθs,νsinϕs,ν,cosθs,ν), (A10) \boldmathnr = \boldmathr|\boldmathrs,ν|=|\boldmathr||\boldmathrs,ν|(sinθcosϕ,sinθsinϕ,cosθ), (A11) (\boldmathnr)z = z|\boldmathrs,ν|=|\boldmath% r||\boldmathrs,ν|cosθ, (A12) (\boldmathnr)⊥z = \boldmathr⊥z|\boldmathrs,ν|=|\boldmathr||\boldmathrs,ν|(cosϕ,sinϕ,0), (A13) \boldmathks,ν = (cosθs,νcosϕs,ν,cosθs,νsinϕs,ν,−sinθs,ν), (A14) \boldmathls,ν = (−sinθs,νsinϕs,ν,sinθs,νcosϕs,ν,cos