NEUTRINO PAIR ANNIHILATION IN COLLAPSARS: A RAY-TRACING METHOD IN SPECIAL RELATIVITY

Accumulating observational evidence has been reported, identifying a massive stellar collapse as the origin 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 Type Ib/c supernova (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 collapsars have received much 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) and Meszaros & Rees (1992) pioneeringly proposed that the energy deposition proceeds predominantly via neutrino and antineutrino annihilation into electron and positron (e.g., $\nu + {\bar{\nu }} \rightarrow e^{-}+ e^{+}$, 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 10¹⁵ G 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. 2004; 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 many baryons, they may be related to the X-ray flashes, which are a low energy analog 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 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 many 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), which is 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). But 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 double angular integration for neutrino and antineutrino for all the grid points in three dimensions, which is still highly time consuming and prevents its application to hydrodynamical simulation. Ruffert & Janka (1998) have made this scheme 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.

Using this method, Nagataki et al. (2007) have estimated the neutrino heating rates and included them in the hydrodynamical simulation. For reducing the computational time, they added one more assumption of the optical 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 z-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 spatial 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 BH.

More recently, Dessart et al. (2009) have developed a new scheme to estimate the energy deposition rates via neutrino pair annihilation using the mutli-angle neutrino-transport solver (Ott et al. 2008). They discussed the possibility of the neutrino-driven outflow production in the postmerger phase of binary neutron-star coalescence. Without the general relativistic (GR) effects, this is the state-of-the-art method to estimate the neutrino pair annihilation. As for the GR 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 2 by the GR 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 incorporated into the hydrodynamical simulations of collapsars. Relying on the neutrino leakage scheme as in Ruffert et al. (1997) and 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 GR effects (e.g., Birkl et al. 2007) for simplicity. Thus our scheme is limited to capturing 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 give for the first time 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 start-of-the-art in comparison with Birkl et al. (2007) and 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 their 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.

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) backward till it hits 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 freely emerge. As indicated in this figure, some fraction of neutrinos comes out from 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 toward the precise determination of the heating rates.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray} &&\frac{dq_{\nu {\bar{\nu }}}^+({\bm r})}{dt dV} =\! \int\!\!\int f_{\nu }(p_{\nu },{\bm r}) f_{\bar{\nu }}(p_{\bar{\nu }}, {\bm r}) \sigma |{\bm v}_{\nu } - {\bm v}_{\bar{\nu }}| (\epsilon _{\nu } + \epsilon _{\bar{\nu }}) d^3 {\bm p}_{\nu } d^3 {\bm p}_{\bar{\nu }},\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where f_ν is the number density of neutrinos in phase space, and p_ν, _ν, v_ν, are the momentum, energy, and velocity of the neutrino, respectively. These definitions work for antineutrinos by changing the notation ν to ${\bar{\nu }}$. σ is the cross section given by

$$\begin{equation} \sigma = 2c^2KG_{\rm F}^2({\bm p}_\nu\, {\bm \cdot}\, {\bm p}_{\bar{\nu }}), \end{equation} \tag{ 2 }$$

where the dimensionless parameter K is written as

$$\begin{equation} K(\nu _e, {\bar{\nu }}_e) = \frac{1+4\sin ^2\theta _{\rm W}+8\sin ^4\theta _{\rm W}}{6 \pi }, \end{equation} \tag{ 3 }$$

$$\begin{equation} K(\nu _\mu, {\bar{\nu }}_\mu) = K(\nu _\tau, {\bar{\nu }}_\tau) = \frac{1-4\sin ^2\theta _{\rm W}+8\sin ^4\theta _{\rm W}}{6 \pi }. \end{equation} \tag{ 4 }$$

Here, the Fermi constant G²_F = 5.29 × 10⁻⁴⁴ cm² MeV⁻², and the Weinberg angle is sin²θ_W = 0.23. Note in our setting of the coordinate system (e.g., Figure 2) that r corresponds to the origin of the O' (targets) system measured in the O system. Variables with subscripts of ν $(\bar{\nu })$ are quantities upon the neutrinosphere, which is labeled by $S_{\nu (\bar{\nu })}$ in Figure 2. They are also measured in the O' 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

$$\begin{equation} dtdV = dt_0dV_0 \end{equation} \tag{ 5 }$$

$$\begin{equation} d\epsilon = \frac{\epsilon }{\epsilon _0}d\epsilon _0 \end{equation} \tag{ 6 }$$

$$\begin{equation} d\Omega = \frac{\epsilon _0^2}{\epsilon ^2} d\Omega _0 \end{equation} \tag{ 7 }$$

$$\begin{eqnarray} \epsilon &=& \gamma (1+\beta \mu _0) \epsilon _0 \nonumber \\ &=& \epsilon _0 / [\gamma (1-\beta \mu)], \end{eqnarray} \tag{ 8 }$$

where β = |V|/c, γ = (1 − β²)^−1/2 and μ = V · n/|V| with V and n being the velocity of fluid and the unit vector describing a given ray of neutrino, within a solid angle of dΩ.

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

$$\begin{eqnarray} \frac{dq_{\nu {\bar{\nu }}}^+({\bm r})}{dt dV} &=& 2 c K G_{\rm F}^2 \int d\theta _{\nu } d \phi _{\nu } d\theta _{\bar{\nu }} d \phi _{\bar{\nu }} \nonumber \\ && \times [ \xi _\nu ^5 ({\bm r},\Omega _\nu) \xi _{\bar{\nu }}^4 ({\bm r},\Omega _{\bar{\nu }}) E_{\nu,0}({\bm r},\Omega _\nu) N_{\bar{\nu },0}({\bm r},\Omega _{\bar{\nu }}) \nonumber \\ && + \xi _\nu ^4 ({\bm r},\Omega _\nu) \xi _{\bar{\nu }}^5 ({\bm r},\Omega _{\bar{\nu }}) N_{\nu,0}({\bm r},\Omega _\nu) E_{\bar{\nu },0}({\bm r},\Omega _{\bar{\nu }}) ] \nonumber \\ && \times \left[ 1-\sin {\theta _{\nu }} \sin {\theta _{\bar{\nu }}} \cos {(\varphi _{\nu }-\varphi _{\bar{\nu }})} -\cos {\theta _{\nu }} \cos {\theta _{\bar{\nu }}} \right]^2\nonumber\\ &&\times \sin \theta _{\nu }\sin \theta _{\bar{\nu }}. \end{eqnarray} \tag{ 9 }$$

Here ξ_ν reflects the special relativistic correction to the neutrino energy as

$$\begin{eqnarray} \xi _\nu ({\bm r},\Omega _\nu) &=& \epsilon _\nu / \epsilon _{\nu,0} \nonumber \\ &=& 1/[\gamma _\nu (1-\mu _{\nu }\beta _{\nu })]. \end{eqnarray} \tag{ 10 }$$

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

$$\begin{equation} \beta _\nu = \frac{|{\bm V}_\nu |}{c}, \end{equation} \tag{ 11 }$$

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

$$\begin{equation} \gamma _\nu = \frac{1}{\sqrt{1-\beta _{\nu }^2}}, \end{equation} \tag{ 12 }$$

and the directional cosine of

$$\begin{equation} \mu _{\nu } = \frac{{\bm n}_\nu \,{\bm \cdot}\, {\bm V}_\nu }{|{\bm V}_\nu |}, \end{equation} \tag{ 13 }$$

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

$$\begin{equation} {\bm n}_\nu = \frac{{\bm p}_\nu }{|{\bm p}_\nu |} = (\sin \theta _\nu \cos \phi _\nu,\ \sin \theta _\nu \sin \phi _\nu, \ \cos \theta _\nu). \end{equation} \tag{ 14 }$$

The following two quantities are the energy-weighted integration of the neutrino distribution function on the neutrino sphere (namely f_ν,0(r_ν,0, p_ν,0)) as follows:

$$\begin{equation} E_{\nu,0}({\bm r},\Omega _\nu) = \int \epsilon _{\nu, 0}^4 f_{\nu,0}({\bm r}_{\nu, 0}, {\bm p}_{\nu,0}) d\epsilon _{\nu,0}, \end{equation} \tag{ 15 }$$

$$\begin{equation} N_{\nu,0}({\bm r},\Omega _\nu) = \int \epsilon _{\nu, 0}^3 f_{\nu,0}({\bm r}_{\nu, 0}, {\bm p}_{\nu,0}) d\epsilon _{\nu, 0},. \end{equation} \tag{ 16 }$$

As in Ruffert et al. (1997), Ruffert & Janka (1998), and Nagataki et al. (2007), the neutrino number flux is simply assumed to be conserved as f_ν(r, p_ν) = f_ν(r_ν,0, p_ν,0) 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 follows:

$$\begin{eqnarray} f_\nu ({\bm r}_{\nu, 0}, {\bm p}_{\nu, 0}) &=& \frac{1}{(hc)^3} \frac{dn_0}{d\epsilon _0 d\Omega _0 dt_0 dV_0} \nonumber \\ &=&\frac{1}{(hc)^3}\frac{1}{\exp (\epsilon _{\nu,0} / kT_{\nu,0}) + 1}, \end{eqnarray} \tag{ 17 }$$

where T_ν,0 is taken to be the same with the matter temperature on the neutrino sphere of T(r_ν). Now the energy integrals in Equations (15) and (16) can be expressed by the Fermi integrals ${\cal F}_k$ as

$$\begin{equation} {\cal F}_k(y) \equiv \int _0^\infty \frac{x^k}{\exp (x-y) + 1} dx, \end{equation} \tag{ 18 }$$

$$\begin{equation} E_{\nu,0}({\bm r},\Omega) = \frac{(kT({\bm r}_\nu))^5}{(hc)^3}{\cal F}_4(0), \end{equation} \tag{ 19 }$$

$$\begin{equation} N_{\nu,0}({\bm r},\Omega) = \frac{(kT({\bm r}_\nu))^4}{(hc)^3}{\cal F}_3(0). \end{equation} \tag{ 20 }$$

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

$$\begin{equation} \int _{R_{\nu _{i}}(\theta)}^{\infty } \frac{1}{\lambda _{i}} dr = 2/3, \end{equation} \tag{ 21 }$$

where $R_{\nu _{i}}(\theta)$ is the radius of the neutrinosphere for the neutrino species of i, namely for $\nu _e,\bar{\nu }_e$, and $\nu _{\mu,\tau }(\bar{\nu }_{\mu,\nu })$ 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, and 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

$$\begin{eqnarray} \frac{d {{\bm m}}_{\nu \bar{\nu }}^+({\bm r})}{dt dV}&=& 2 K G_{\rm F}^2 \int d\theta _{\nu } d \phi _{\nu } d\theta _{\bar{\nu }} d \phi _{\bar{\nu }} \nonumber \\ && \times [ \xi _\nu ^5 ({\bm r},\Omega _\nu) \xi _{\bar{\nu }}^4 ({\bm r},\Omega _{\bar{\nu }}) E_{\nu,0}({\bm r},\Omega _\nu) N_{\bar{\nu },0}({\bm r},\Omega _{\bar{\nu }}) {\bm n}_\nu \nonumber \\ && + \xi _{\bar{\nu }}^5 ({\bm r},\Omega _{\bar{\nu }}) \xi _{\nu }^4 ({\bm r},\Omega _{\nu }) E_{\bar{\nu },0}({\bm r},\Omega _{\bar{\nu }}) N_{\nu,0}({\bm r},\Omega _{\nu }) {\bm n}_{\bar{\nu }} ] \nonumber \\ && \times [ 1-\sin {\theta _{\nu }} \sin {\theta _{\bar{\nu }}} \cos {(\varphi _{\nu }-\varphi _{\bar{\nu }})} -\cos {\theta _{\nu }} \cos {\theta _{\bar{\nu }}} ]^2\nonumber\\ &&\times \sin \theta _{\nu }\sin \theta _{\bar{\nu }}. \end{eqnarray} \tag{ 22 }$$

Finally, the remaining procedure is to change the coordinates of the angle integrations in Equations (9) and (22) from the O' to O coordinate system. Although the formulae are messy (see the Appendix), this procedure, by which the neutrino source and the target can be connected one by one, reduces the numerical costs significantly. The coordinate transformation with respect to the solid angle is written as

$$\begin{eqnarray} d\Omega _{\nu } &=& J_{r \mu }(\mu _{\nu }, \phi _{\nu }, r_{s,\nu }, \mu _{s,\nu }) |_{\phi _{s,\nu } = \phi _{{\rm sph},\nu }} dr_{s,\nu } d\mu _{s,\nu } \nonumber \\ &&+ J_{\mu \phi }(\mu _{\nu }, \phi _{\nu }, \mu _{s,\nu }, \phi _{s,\nu }) |_{r_{s,\nu }=r_{{\rm sph},\nu }} d\mu _{s,\nu } d\phi _{s,\nu } \nonumber \\ &&+ J_{\phi r}(\mu _{\nu }, \phi _{\nu }, \phi _{s,\nu }, r_{s,\nu }) |_{\mu _{s,\nu } = \mu _{{\rm sph},\nu }} d\phi _{s,\nu } dr_{s,\nu }, \end{eqnarray} \tag{ 23 }$$

where r_s,ν(r_s,ν, θ_s,ν, ϕ_s,ν) is the position of S_ν measured in the O system. Using the expressions of Jacobian (J_rμ, J_μϕ, J_ϕr, see the Appendix 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

3.1. Numerical Tests

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

In our 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 shed from a light bulb in spherical and non-spherical geometry (Sections 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 of the application to the collapsar simulations in Section 3.1.3.

In the following tests, our numerical domain has the spherical coordinates with 100 equally spaced radial meshes covering from 0 to 500 km in radius. To check the numerical convergence, we vary the numerical resolution for the θ and ϕ directions, respectively. For simplicity, the axial and equatorial symmetry is assumed for the hydrodynamic quantities (such as density, electron fraction, and so on) as a background. The density and temperature on the neutrino sphere are set to be ρ = 10¹² g cm⁻³ and T = 5 MeV, respectively.

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

$$\begin{equation} E \equiv \frac{L_{1,{\rm err}}}{L_{1,{\rm ana}}}, \end{equation} \tag{ 24 }$$

$$\begin{equation} L_{1,{\rm err}} \equiv \Sigma | {x_{\rm num} - x_{\rm ana}} |, \end{equation} \tag{ 25 }$$

$$\begin{equation} L_{1,{\rm ana}} \equiv \Sigma {x_{\rm ana}}, \end{equation} \tag{ 26 }$$

where x_num represents the numerical values such as the energy deposition rate $(q^+_{\nu {\bar{\nu }}})$ and the absolute value of the momentum transfer rate $(|{\bm m}^+_{\nu {\bar{\nu }}}|)$, while x_ana 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 E 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 E2 as

$$\begin{equation} E2 \equiv \frac{M_{\rm num}}{M_{\rm max}}, \end{equation} \tag{ 27 }$$

$$\begin{equation} M_{\rm num} \equiv \Sigma \big| {\bm m}^+_{\nu {\bar{\nu }}} \big|, \end{equation} \tag{ 28 }$$

$$\begin{equation} M_{\rm max} \equiv \Sigma m^+_{\rm max}, \end{equation} \tag{ 29 }$$

where m⁺_max is the analytical solution, which is defined from the energy deposition rate as m⁺_max ≡ q⁺_ana/c, thus meaning the maximum momentum transfer. With this E2, we will evaluate how precisely the cancellation can be maintained.

3.1.1. Radiation from Spherical Neutrino Sphere

In the 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 foci of this paper.

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

$$\begin{equation} \frac{dq_{\nu {\bar{\nu }}}^+({\bm r})}{dt dV} = \frac{ 4 c K G_{\rm F}^2 }{3} \frac{(kT)^9}{(hc)^6}{\cal F}_4(0){\cal F}_3(0) F(r), \end{equation} \tag{ 30 }$$

and

$$\begin{equation} \frac{d{\bf m}_{\nu {\bar{\nu }}}^+({\bm r})}{dt dV} = \frac{ 4 K G_{\rm F}^2 }{3} \frac{(kT)^9}{(hc)^6}{\cal F}_4(0){\cal F}_3(0) G(r) {{\bm n}}_r, \end{equation} \tag{ 31 }$$

where geometrical factor F(r) and G(r) are given by

$$\begin{equation} F(r) = \frac{2\pi ^2}{3} (1-X)^4 (5 + 4X + X^2) \end{equation} \tag{ 32 }$$

and

$$\begin{equation} G(r) = \frac{\pi ^2}{6} (1-X)^4 (1+X)(8 + 9X + 3X^2), \end{equation} \tag{ 33 }$$

with

$$\begin{equation} X = \sqrt{1-\left(\frac{R_\nu }{r}\right)}, \end{equation} \tag{ 34 }$$

where R_ν 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 labels such as "20 × 16" at the top right edge indicates the resolution of our ray-tracing calculation, here n_θ = 60, n_ϕ = 64, the equally spaced mesh numbers in the lateral (θ) and azimuthal directions (ϕ), respectively. Comparing the two panels, 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 n_θ, n_ϕ 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 J_μϕ (from Equations (A1) to (A16)) and its implementation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{equation} \frac{dq_{\nu {\bar{\nu }}}^+({\bm r})}{dt dV} = \frac{ 4 c K G_{\rm F}^2 }{3} \frac{(kT)^9}{\gamma ^9 (hc)^6}{\cal F}_4(0){\cal F}_3(0) F(r), \end{equation} \tag{ 35 }$$

$$\begin{equation} \frac{d{\bf m}_{\nu {\bar{\nu }}}^+({\bm r})}{dt dV} = \frac{ 4 K G_{\rm F}^2 }{3} \frac{(kT)^9}{\gamma ^9 (hc)^6}{\cal F}_4(0){\cal F}_3(0) G(r) {{\bm n}}_r, \end{equation} \tag{ 36 }$$

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

Figure 5 shows the energy deposition (left) and momentum transfer rates (right) from the relativistically rotating neutrino sphere. Here we assume γ = 2 with vanishing V_r and V_θ. Note in this figure that only the numerical results are shown, owing to the inaccessibility of the analytical solution as mentioned above. It is shown that the annihilation rates both for energy and momentum are greatly reduced near 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 (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 (n_θ) should be at least more than 20 to maintain the percent levels of agreement with the analytic solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1.2. Spherical Cavity

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

$$\begin{equation} F(r) = \frac{64 \pi ^2}{3}. \end{equation} \tag{ 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 E1 and E2 errors for the energy (left) and momentum transfer rate (right), it can be seen that the errors mainly depend on n_ϕ. By setting n_ϕ = 32, we can get the agreement with an accuracy of 1% for E and 0.1% for E2, which would satisfy the requirement of the vanishing momentum transfer sufficiently.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 drdϕ 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 θ < π/4 and θ>3π/4, the neutrino sphere is situated at 300 km, while for π/4 ⩽ θ ⩽ 3π/4, at 100 km. The 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 as 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 E and E₂ is shown to be smaller than the ones from the previous tests due to the relatively complicated geometry. By taking the resolution of (n_θ, n_ϕ) ≳ (32, 32), our code can reproduce the analytic solutions with percent levels of accuracy here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 Equations (A1) to (A16)), its implementation, and the ray-tracing calculation. From those test calculations, it is suggested that the grid numbers of (n_θ, n_ϕ) ≳ (30, 30) 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 (n_r, n_θ, n_ϕ) = (300, 40, 32) 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 collapsars 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 collapsars (Harikae et al. 2009). Without repeating the numerical procedures again, we will briefly 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 Ω(r, θ) analytically as

$$\begin{equation} \Omega (r,\theta)=\frac{\Omega _0 X_0^4 + \alpha \Omega _\mathrm{lso}(M(X))X^4}{X_0^4+X^4}, \end{equation} \tag{ 38 }$$

where X = rsin θ, Ω_lso(M(X)) = j_lso(M(X))/X². M(X) is the mass coordinate at X, and j_lso and Ω_lso are the specific angular momentum and angular velocity in the last stable orbit. Model parameters are α, Ω₀, and X₀. This distribution provides co-rotation with Ω = Ω₀ at X < X₀ connected to constant specific angular momentum j = αj_lso at X > X₀. X₀ would correspond to the size of a co-rotating region in the pre-collapse phase, such as an Fe core. To be closer to the original rotational profile in Woosley & Heger (2006), we set α = 0.8, Ω₀ = 1.2 rad s⁻¹, and X₀ = R_Fe, where R_Fe is the size of the Fe core (≈3000 km 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 the BH, was set at 50 km (Harikae et al. 2009). This manipulation may lead to the underestimation of the neutrino luminosity by excising the radiation from the more inner edge of the accretion disk. This is obviously disadvantageous for producing the neutrino-driven outflow. In this paper, we take a more compact inner boundary as max(10 km, 2r_g), extending to the outer boundary of 30,000 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–3 M_☉ in the center. This indicates the BH 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 10⁵³ erg s⁻¹, which consists of 2.2 × 10⁵² erg s⁻¹ for ν_e, 9.4 × 10⁵² erg s⁻¹ for ${\bar{\nu }}_{\rm e}$, 5.1 × 10⁵¹ erg s⁻¹ for ν_μ and ν_τ. Afterward, 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 typically reaches ∼10 MeV, which provides the high neutrino luminosity of 10⁵³ erg s⁻¹. From the right panel, the Lorentz factor at the inner edge of the disk becomes γ ≈ 1.3, 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.

Zoom In Zoom Out Reset image size

3.2.3. Energy Deposition and Momentum Transfer Rates in Collapsar

The left- and right-hand sides of the left panel of Figure 13 show the energy deposition and the momentum transfer rates (their absolute value). The energy and momentum transport from the accretion disk (regions colored black) to the polar funnel regions (regions colored red) can apparently be 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 10³⁰ erg s⁻¹ cm⁻³, where the momentum transfer rates are the order of 10¹⁹ erg cm⁻⁴. Thus, the momentum transfer rates multiplied by the speed of light c become comparable to the energy deposition rate there.

Zoom In Zoom Out Reset image size

For our model, the density along the polar funnel regions is as low as 10^5–6 g cm⁻³ (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 10^1–2. This suggests that the momentum transfer plays as important a role as the energy deposition for the efficient acceleration of the neutrino-driven outflows. The right panel of Figure 13 is the annihilation rate without the special relativistic correction (ξ_ν = 1 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 6.22 and 7.62/(10⁵⁰ erg s⁻¹) with and without the correction, respectively. These values are comparable to what was expected by Nagataki et al. (2003), since the mass accretion rate in this epoch is about 0.1 M_☉ s⁻¹. 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 characterized 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 τ_heat ≡ ρΦ/q⁺, where Φ is the local gravitational potential. The dynamical timescale is defined as $\tau _{\rm {dyn}} \equiv \sqrt{3\pi /16G\bar{\rho }}$, where $\bar{\rho }$ is the average density at a certain radius and we take $\bar{\rho }(r) \equiv 3 M(r) /4 \pi r^3$. Figure 14 depicts the ratio of the dynamical τ_dyn to the heating timescales τ_heat, 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.

Zoom In Zoom Out Reset image size

Figure 15 shows various timescales for the material along the rotational axis in Figure 14. τ_rel ≡ ρc²/q⁺ characterizes the timescale for matter to become relativistic by the neutrino heating. As mentioned, the first criterion of τ_dyn > τ_heat is satisfied for the regions inside ∼80 km (compare green and blue lines). More interestingly, inside ∼50 km, τ_rel gets slightly shorter than τ_dyn (compare red and blue lines), indicating that matter could attain relativistic motions there. To see these outcomes requires the implementation of the ray-tracing calculation in the hydrodynamic calculation, which we pose as the next task of this study. It is noted that with our choice of the grid numbers ((n_r, n_θ, n_ϕ) = (300, 40, 32)), it generally takes several minutes (in 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 τ_dyn > τ_heat condition, 10,000,000 hydro-steps are generally needed because the hydrodynamical timescale is an order of 10⁻⁶ 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 τ_dyn/τ_heat during the simulations in the computational domain. After the ray-tracing calculation is turned on, seeing τ_dyn/τ_heat ≳ 1, the neutrino-driven outflows will be launched within 0.1 s (S. Harikae et al. 2010 in preparation). Thus, the computational time using the above supercomputer can be reduced to several months, which is expected to be quite feasible.

Zoom In Zoom Out Reset image size

Bearing in mind the application to the collapsar models of GRBs, 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 ∼9 s after the onset of gravitational collapse. The accretion disk later on can be a luminous source of neutrinos as high as ∼10⁵³ erg s⁻¹. 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 a role as the energy deposition for the efficient acceleration of the neutrino-driven outflows.

Finally, we shall make comparisons with previous work and also mention some imperfections of this study. As mentioned, GR 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 2000, 2001; 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 two-dimensional configurations). In fact, the total heating rates are Q^tot = 2.0 × 10⁴⁸ and 1.5 × 10⁴⁸(erg s⁻¹) for the Newtonian (model DN) and the special relativistic case (model DN4L), which are in Birkl et al. (2007), 1.5 × 10⁴⁸ and 1.3 × 10⁴⁸(erg s⁻¹), 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 GR formulation of radiative transport, which we will investigate as an 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 an 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, changes 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).

Zoom In Zoom Out Reset image size

Furthermore, the neutrino oscillation by Mikheyev–Smirnov–Wolfenstein (MSW) effect (see collective references in Kotake et al. 2006 and Kawagoe et al. 2009) could be important, albeit in a much later phase than we considered in this paper. It should be remembered that the typical density in the polar funnel was ρ ∼ 10⁵ g cm⁻³ in the present model. When the density there drops as low as ρ ≲ 10³ g cm⁻³ 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 remain 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, 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 these 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 the 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).

In this section, we give an exact expression of the Jacobians (J_rμ, J_μϕ, J_ϕr), namely,

$$\begin{equation} J_{r \mu }(\mu _{\nu }, \phi _{\nu }, r_{s,\nu }, \mu _{s,\nu }) \equiv \left| \begin{array}{@{}cc@{}} \dsty \frac{\partial {\mu _{\nu } }}{\partial {r_{s,\nu }}} & \dsty \frac{\partial {\mu _{\nu } }}{\partial {\mu _{s,\nu }}}\\\\[-6pt] \dsty \frac{\partial {\phi _{\nu }}}{\partial {r_{s,\nu }}} & \dsty \frac{\partial {\phi _{\nu }}}{\partial {\mu _{s,\nu }}}\\ \end{array} \right|, \end{equation} \tag{ A1 }$$

$$\begin{equation} J_{\mu \phi }(\mu _{\nu }, \phi _{\nu }, \mu _{s,\nu }, \phi _{s,\nu }) \equiv \left| \begin{array}{@{}cc@{}} \dsty \frac{\partial {\mu _{\nu } }}{\partial {\mu _{s,\nu }}} &\dsty \frac{\partial {\mu _{\nu } }}{\partial {\phi _{s,\nu }}}\\\\[-6pt] \dsty \frac{\partial {\phi _{\nu }}}{\partial {\mu _{s,\nu }}} &\dsty \frac{\partial {\phi _{\nu }}}{\partial {\phi _{s,\nu }}}\\ \end{array} \right|, \end{equation} \tag{ A2 }$$

$$\begin{equation} J_{\phi r}(\mu _{\nu }, \phi _{\nu }, \phi _{s,\nu }, r_{s,\nu }) \equiv \left| \begin{array}{@{}cc@{}} \dsty \frac{\partial {\mu _{\nu } }}{\partial {\phi _{s,\nu }}} &\dsty \frac{\partial {\mu _{\nu } }}{\partial {r_{s,\nu }}}\\\\[-6pt] \dsty \frac{\partial {\phi _{\nu }}}{\partial {\phi _{s,\nu }}} &\dsty \frac{\partial {\phi _{\nu }}}{\partial {r_{s,\nu }}}\\ \end{array} \right| \end{equation} \tag{ 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 straightforwardly as

$$\begin{equation} \frac{\partial {\mu _{\nu }}}{\partial {r_{s,\nu }}} = \frac{1}{r_{s,\nu }} \frac{ (1 + {{\bm n}}_{s,\nu } \cdot {{\bm n}}_r)({{\bm n}}_r)_z + (|{{\bm n}}_r|^2 - 1)\cos \theta _{s,\nu } }{(1-2{{\bm n}}_{s,\nu } \cdot {{\bm n}}_r + |{{\bm n}}_r|^2)^{\frac{3}{2}}}, \end{equation} \tag{ A4 }$$

$$\begin{eqnarray} &&\frac{\partial {\mu _{\nu }}}{\partial {\mu _{s,\nu }}} \nonumber\\ &&= \frac{\sin \theta _{s,\nu } (1-2{{\bm n}}_{s,\nu } \cdot {{\bm n}}_r + |{{\bm n}}_r|^2) - (\cos \theta _{s,\nu } - ({{\bm n}}_r)_z) {{\bm n}}_r \cdot {{\bm k}}_{s,\nu } }{ \sin \theta _{s,\nu } (1-2{{\bm n}}_{s,\nu } \cdot {{\bm n}}_r + |{{\bm n}}_r|^2)^{\frac{3}{2}}},\nonumber\\ \end{eqnarray} \tag{ A5 }$$

$$\begin{equation} \frac{\partial {\mu _{\nu }}}{\partial {\phi _{s,\nu }}} = \frac{ (\cos \theta _{s,\nu } - ({{\bm n}}_r)_z) ({{\bm n}}_r)_{\bot z} \cdot {{\bm l}}_{s,\nu } }{(1-2{{\bm n}}_{s,\nu } \cdot {{\bm n}}_r + |{{\bm n}}_r|^2)^{\frac{3}{2}}}, \end{equation} \tag{ A6 }$$

$$\begin{eqnarray} &&\frac{\partial {\phi _{\nu }}}{\partial {r_{s,\nu }}} = \frac{1}{r_{s,\nu }}\nonumber\\ &&\times \frac{ (|({{\bm n}}_r\!)_{\bot z}|^2 \,{-}\, {{\bm n}}_{s,\nu } \cdot ({{\bm n}}_r\!)_{\bot z}) ({{\bm n}}_{s,\nu })_x \,{+}\, (\sin ^2 \theta _{s,\nu } \,{-}\, {{\bm n}}_{s,\nu } \cdot ({{\bm n}}_r\!)_{\bot z}) ({{\bm n}}_r\!)_x }{ (({{\bm n}}_r)_y \,{-}\, ({{\bm n}}_{s,\nu })_y) (\sin ^2 \theta _{s,\nu } \,{-}\, 2{{\bm n}}_{s,\nu } \cdot ({{\bm n}}_r)_{\bot z} \,{+}\, |({{\bm n}}_r)_{\bot z}|^2)}\!,\nonumber\\ \end{eqnarray} \tag{ A7 }$$

$$\begin{equation} \frac{\partial {\phi _{\nu }}}{\partial {\mu _{s,\nu }}} = \frac{-1}{\sin \theta _{s,\nu }} \frac{ ({{\bm n}}_r)_{\bot z} \cdot {{\bm m}}_{s,\nu }}{\sin ^2 \theta _{s,\nu } - 2{{\bm n}}_{s,\nu } \cdot ({{\bm n}}_r)_{\bot z} + |({{\bm n}}_r)_{\bot z}|^2 }, \end{equation} \tag{ A8 }$$

$$\begin{equation} \frac{\partial {\phi _{\nu }}}{\partial {\phi _{s,\nu }}} = \frac{ \sin ^2 \theta _{s,\nu } - {{\bm n}}_{s,\nu } \cdot ({{\bm n}}_r)_{\bot z} }{\sin ^2 \theta _{s,\nu } - 2{{\bm n}}_{s,\nu } \cdot ({{\bm n}}_r)_{\bot z} + |({{\bm n}}_r)_{\bot z}|^2 }, \end{equation} \tag{ A9 }$$

where the following values are dependent only on the position of the neutrino sphere (θ_s,ν, ϕ_s,ν) and the direction of a specified direction (θ, ϕ) as

$$\begin{equation} {{\bm n}}_{s,\nu } = \frac{{{\bm r}}_{s,\nu }}{|{{\bm r}}_{s,\nu }|} = (\sin \theta _{s,\nu } \cos \phi _{s,\nu }, \sin \theta _{s,\nu } \sin \phi _{s,\nu }, \cos \theta _{s,\nu }), \end{equation} \tag{ A10 }$$

$$\begin{equation} {{\bm n}}_r = \frac{{{\bm r}}}{|{{\bm r}}_{s,\nu }|} = \frac{|{{\bm r}}|}{|{{\bm r}}_{s,\nu }|}(\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta), \end{equation} \tag{ A11 }$$

$$\begin{equation} ({{\bm n}}_r)_{z}= \frac{z}{|{{\bm r}}_{s,\nu }|} = \frac{|{{\bm r}}|}{|{{\bm r}}_{s,\nu }|} \cos \theta, \end{equation} \tag{ A12 }$$

$$\begin{equation} ({{\bm n}}_r)_{\bot z}= \frac{{{\bm r}}_{\bot z}}{|{{\bm r}}_{s,\nu }|} = \frac{|{{\bm r}}|}{|{{\bm r}}_{s,\nu }|}(\cos \phi, \sin \phi, 0), \end{equation} \tag{ A13 }$$

$$\begin{equation} {{\bm k}}_{s,\nu } = (\cos \theta _{s,\nu } \cos \phi _{s,\nu }, \cos \theta _{s,\nu } \sin \phi _{s,\nu }, {-}{\rm sin} \theta _{s,\nu }), \end{equation} \tag{ A14 }$$

$$\begin{equation} {{\bm l}}_{s,\nu } = ({-}{\rm sin} \theta _{s,\nu } \sin \phi _{s,\nu }, \sin \theta _{s,\nu } \cos \phi _{s,\nu }, \cos \theta _{s,\nu }), \end{equation} \tag{ A15 }$$

$$\begin{equation} {{\bm m}}_{s,\nu } = ({-}{\rm cos} \theta _{s,\nu } \sin \phi _{s,\nu }, \cos \theta _{s,\nu } \cos \phi _{s,\nu }, {-} {\rm sin} \theta _{s,\nu }). \end{equation} \tag{ A16 }$$

Finally, it should be noted that above equations can be very simple along the polar axis as,

$$\begin{equation} \frac{\partial {\mu _{\nu }}}{\partial {r_{s,\nu }}} = \frac{1}{r_{s,\nu }} \frac{ \big(\frac{Z}{r_{s,\nu }} - \cos \theta _{s,\nu } \big) \frac{Z}{r_{s,\nu }} \cos \theta _{s,\nu }}{\big(1 - 2 \frac{Z}{r_{s,\nu }} \cos \theta _{s,\nu } + \big|\frac{Z}{r_{s,\nu }}\big|^2 \big)^{\frac{3}{2}}}, \end{equation} \tag{ A17 }$$

$$\begin{equation} \frac{\partial {\mu _{\nu }}}{\partial {\mu _{s,\nu }}} = \frac{1 - \frac{Z}{r_{s,\nu }} \cos \theta _{s,\nu }}{\big(1 - 2 \frac{Z}{r_{s,\nu }} \cos \theta _{s,\nu } + \big|\frac{Z}{r_{s,\nu }}\big|^2 \big)^{\frac{3}{2}}}, \end{equation} \tag{ A18 }$$

$$\begin{equation} \frac{\partial {\mu _{\nu }}}{\partial {\phi _{s,\nu }}} = 0, \end{equation} \tag{ A19 }$$

$$\begin{equation} \frac{\partial {\phi _{\nu }}}{\partial {r_{s,\nu }}} = 0, \end{equation} \tag{ A20 }$$

$$\begin{equation} \frac{\partial {\phi _{\nu }}}{\partial {\mu _{s,\nu }}} = 0, \end{equation} \tag{ A21 }$$

$$\begin{equation} \frac{\partial {\phi _{\nu }}}{\partial {\phi _{s,\nu }}} = 1, \end{equation} \tag{ A22 }$$

where Z is the distance from the equatorial plane.