Neutrino heating near hyper-accreting black holes

Hyper-accretion discs around black holes emit copious neutrinos and anti-neutrinos. A fraction of the emitted neutrinos convert to electron-positron plasma above the disc through the annihilation reaction $\nu\bar\nu\to e^+e^-$. This process may drive relativistic jets associated with GRB explosions. We calculate the efficiency of energy deposition by neutrinos. Our calculation is fully relativistic and based on a geodesic-tracing method. We find that the efficiency of neutrino heating is a well-defined function of (i) accretion rate and (ii) spin of the black hole. It is practically independent of the details of neutrino transport in the opaque zone of the disc. The results help identify accretion discs whose neutrino emission can power GRBs.


INTRODUCTION
A plausible model for the central engines of gamma-ray bursts (GRBs) pictures a transient, hyper-accreting disc formed around a rotating black hole (see Beloborodov 2008 for a recent review). The disc is a source of copious neutrinos and anti-neutrinos, which partially annihilate above the disc and turn into e ± pairs, ν +ν → e − + e + . This process was proposed as a possible mechanism for creating relativistic, e ± -dominated jets that could power observed GRBs (Eichler et al. 1989). However, its efficiency remained unsettled. Its accurate calculation requires a detailed relativistic model for the neutrino source -the accretion disc -as well as tracing the neutrino propagation in the Kerr spacetime of the black hole. A detailed relativistic model for GRB discs was completed recently (Chen & Beloborodov 2007, hereafter CB07). It describes the disc down to its inner edge and gives accurate energy fluxes carried away by ν andν at all radii r. In the present work we trace the neutrino trajectories and calculate the rate of νν annihilation around the disc. Neutrino annihilation was previously calculated in a number of works (e.g. Popham, Woosley & Fryer 1999.; Asano & Fukuyama 2001;Birkl et. al. 2007). Our work has three motivations: (i) A relativistic calculation of νν annihilation has never been done for a realistic accretion disc around a spinning black hole. Previous works either used a toy model for neutrino source (e.g. Birkl et al. 2007) or replaced neutrino trajec-E-mail: izalamea@phys.columbia.edu † Also at Astro-Space Center of Lebedev Physical Institute, Profsojuznaja 84/32, Moscow 117810, Russia tories by straight lines (Popham et. al. 1999). (ii) The efficiency of νν annihilation strongly depends on the accretion rateṀ and the spin parameter a of the black hole. It is desirable to know this dependence and identify the range ofṀ and a where νν annihilation can provide the observed energy of GRB explosions. (iii) In addition to νν annihilation, neutrinos can create e ± pairs off magnetic field (Kuznetsov & Mikheev 1997;Gvozdev & Ognev 2001). The contribution of this process to the energy deposition rate should be included.
The paper is organized as follows. Section 2 describes the setup and method of our calculations. The results are presented in Section 3 and summarized in Section 4.

MODEL DESCRIPTION
Neutrinos emitted by the disc follow geodesics in Kerr spacetime. The efficiency of their annihilation can be calculated numerically by tracing the geodesics, evaluating the local energy deposition rate [erg s −1 cm −3 ] everywhere around the black hole and then integrating over volume outside the disc to obtain the net energy deposition rateĖνν (energy at infinity per unit time at infinity).
The neutrino emission and annihilation is concentrated near the black hole, where the accretion time-scale is short andṀ may be assumed to be quasi-steady. ThenĖνν depends on four parameters that specify the steady disc model: accretion rateṀ , viscosity parameter α, mass of the black hole M , and spin of the black hole a.

Neutrino source: disc model
As the disc matter spirals into the black hole, it is viscously heated: the gravitational energy is converted to heat. Part of the heat is lost to neutrino emission, and part is stored in the disc and distributed between nuclear matter, radiation, and e ± pairs, in perfect thermodynamic and nuclear statistical equilibrium. The equilibrium microphysics is determined by only three parameters: temperature T , baryon mass density ρ, and electron fraction Ye (equal to the charged nucleon fraction). Other parameters -e.g. the electron chemical potential µe and density of e ± pairs n± -are derived from T , ρ and Ye.
The neutrino emission peaks in the inner region of the disc. The far dominant emission mechanism is the e − /e + capture onto protons/neutrons, The neutrino-cooled discs are nearly in β-equilibrium: the rates of the two reactions in equation (1) are practically equal, which determines the value of Ye(ρ, T ) (Imshennik et al. 1967;Beloborodov 2003). In principle, all three flavors of neutrinos could be emitted from the disc, but only electron neutrinos need to be considered; the emission rates for muon and tau neutrinos are negligible (CB07). A detailed model for neutrino-cooled relativistic discs was developed in CB07, and we use their model in our calculations. The approximate hydrodynamic disc equations are solved with the vertically-integrated α prescription. The equations include radial transport of heat and lepton number. Local microphysics is treated exactly: nuclear composition, electron degeneracy, neutrino emissivity and opacity etc., using the equilibrium distribution functions for all species except neutrinos. Neutrinos are modeled separately in the opaque and transparent zones of the disc, matching at the transition between the two zones.
The disc model provides the vertically averaged T , ρ, Ye, µe (which are approximately equal to their values in the midplane θ = π/2) and the half-thickness of the disc H. The ν andν spectra emitted from the neutrino-transparent (optically thin) zone of the disc are given by for neutrinos (E > 0), and for anti-neutrinos (E > q + 1). Here h = 2πh, λ = h/mec, q = (mn − mp)/me = 2.53, K = 6.5 × 10 −4 s −1 , Θ = kT /mec 2 , µe is electron chemical potential in units of mec 2 and E is neutrino/anti-neutrino energy in units of mec 2 . The spectra emerging from the opaque zone are controlled by neutrino transport through the disc, which cannot be reliably calculated -it depends on the unknown vertical distribution of viscous heating. Fortunately, the νν annihilation rate depends only on the energy fluxes Fν and Fν from the disc surface, which are insensitive to the neutrinotransport details. This fact has a simple analytical explanation (Beloborodov 2008). It follows from the proportionality σνν ∝ Eν Eν , where σνν is the cross-section for annihilation for neutrino and anti-neutrinos with energies Eν and Eν .
To demonstrate that the rate of νν annihilation depends on Fν and Fν but not on the exact shapes of ν andν spectra, we consider below two extreme models A and B, and they give practically the same annihilation rates: Model A: Neutrinos and anti-neutrinos are emitted with the same spectra as found inside the disc. The spectra are normalized so that the emerging emission carries away the known energy fluxes Fν and Fν . In the opaque region, the spectra of ν andν are described by Fermi-Dirac distributions. The temperature and chemical potential for ν and ν inside the disc are obtained from the numerical models of CB07.
Model B: Neutrinos and anti-neutrinos are emitted with a thermal spectrum with zero chemical potential and temperature T = T eff , where T eff is the effective surface temperature defined by (7/8)σT 4 eff = Fν + Fν (σ is the Stefan-Boltzmann constant and the coefficient 7/8 takes into account the difference between statistics of photons and ν,ν).
When the disc is efficiently cooled (neutrino energy flux almost balances viscous heating), T eff is given by the standard thin-disc model of Page & Thorne (1974): T eff ≈ T st eff . This regime occurs in a broad range of accretion rateṡ Mign <Ṁ <Ṁtrap (CB07). IfṀ <Ṁign, the disc temperature is not high enough to ignite the neutrino emitting reactions. IfṀ >Ṁtrap, the emitted neutrinos become trapped in the disc and advected into the black hole. Our third (simplest) model for the neutrino source is defined as follows.
Model C: Neutrinos and anti-neutrinos are emitted with a thermal spectrum that has zero chemical potential and the following temperature, This model does not even require the calculation of the disc structure as T st eff is a known analytical function of r (Page & Thorne 1974). As we show below, this simplest model gives remarkably accurate result forĖνν .
The characteristic accretion ratesṀign andṀtrap depend on the viscosity parameter α. They were calculated in CB07, and the numerical results are well approximated by the following formulae, The coefficients Kign and Ktrap are functions of the blackhole spin a. For a = 0, In all three models A, B and C the neutrino emission is assumed to be isotropic in the local rest frame of the disc (which is in Keplerian rotation around the black hole).

Neutrino transport
To evaluate the νν annihilation rate at a given point we need to know the local ν andν distribution functions. They can be obtained using the known neutrino distribution functions at the surface of the disc. To a first approximation, neutrinos obey the collisionless Boltzmann equation, because most of them do not participate in any interactions (and eventually escape to infinity or get captured into the black hole). The Boltzmann equation in curved spacetime has the same form as in flat spacetime. It states that the phase-space density (or occupation number) of neutrinos remains constant along their trajectories, Here x µ (λ) is a parameterized worldline for a ν orν. The neutrinos emitted by the disc have huge energies compared to their rest mass and we treat them as massless particles propagating along null geodesics in Kerr spacetime. We use Boyer-Lindquist coordinates x α = (t, φ, r, θ), where the Kerr metric has the form, The metric tensor g αβ is specified by two parameters: rg = 2GM/c 2 and the spin parameter |a| < 1; it is given in e.g. Chandrasekhar (1998). Equation (6) is covariant and its solution takes into account Doppler and gravitational redshifts. The null geodesics in Boyer-Lindquist coordinates are described by the known ordinary differential equations of first order (e.g. Chandrasekhar 1998) which we solve numerically. Figure 1 shows the image of the accretion disc observed from three locations near the black hole. Colour represents the redshift (or blueshift) of neutrinos as they propagate from the emission point on the disc to the observation point. The asymmetry of the images is caused by the rotation of the black hole.

Local rates of
The rate of νν annihilation at a given point depends only on the local momentum distribution of ν andν. We use the phase-space occupation number to describe this distribution; it is denoted by f for neutrinos and by f for antineutrinos. The local energy-momentum deposition rate is given by (see e.g. Birkl et. al. 2007), where where p and p are the four-momenta of neutrino and antineutrino, θ is the angle between the spatial (3D) components of p and p , σ0 = 4m 2 e G 2 F /πh 4 ≈ 1.71 × 10 −44 cm 2 , GF is Fermi constant, C1 = 0.78 and C2 = 1.06. We neglect the mass of ν andν, i.e. assume p α pα = 0 and p α p α = 0. Figure 1. Projection of the sky for three different observers receiving neutrinos from a disc of radius 140r h . The colour codes the ratio of the energy at reception to the energy at emission for neutrinos, as measured in the local ZAMO frames (Section 2.4). All three observers are located at 3.1 horizon radii r h around a black hole of spin a = 0.95. Their θ positions are approximately 0, π/4 and π/2. The observer sky is parameterized with two angles α ∈ (0, π) and β ∈ (0, 2π) in the ZAMO frame. The axis α = 0 points into the black hole and α = π points away from the black hole. The x and y coordinates in the figure are (x, y) = ( α π sin β, α π cos β).

Reaction ν → νe + e − in a strong magnetic field
Neutrinos propagating in a strong magnetic field B can create electron-positron plasma through the process ν → νe + e − (Kuznetsov & Mikheev 1997;Gvozdev & Ognev 2001). A neutrino with four-momentum p α and energy E deposits energy-momentum into the e ± plasma with the following rate where ψ is the angle between the magnetic field and the neutrino momentum, cv ≈ 0.96 and ca = 1/2. The local energy-momentum deposition rate via reactions ν → νe + e − and ν → ν e + e − is given by

Integration of the energy deposition rate over volume
The local energy-momentum deposition rates (eqs. 8 and 11) can be calculated in any frame of reference. We chose the frame of Zero Angular Momentum Observer, ZAMO. These observers do not move in r or θ directions. Their φ motion corresponds to zero angular momentum (p φ = 0). We use ZAMO frames because they are well defined everywhere outside the black-hole horizon (no observers can be static in Boyer-Lindquist coordinates inside the ergosphere). The local orthonormal coordinates of ZAMO, dx α , are related to dx α = (dt, dφ, dr, dθ) by where detA = √ −g and the element of four-dimensional volume measured by ZAMO dṼ is related to the Boyer-Lindquist coordinate volume dV = dt dφ dr dθ by Let dP α = Q α dṼ be four-momentum deposited in dṼ as measured by ZAMO. The deposited energy measured by a distant observer is dE = −c dPt. It can be expressed in terms of dP α , where we used the transformation dP α = A α β dP β . This yields The net energy deposition rate outside the horizon dE/dt (measured by a distant observer) is given bẏ where r h is the horizon radius.

Numerical Method
We compute the local rates of four-momentum deposition Q α on a spatial grid. The problem is axially symmetric and the grid is set on the (r, θ) plane. It covers the region r h < r < rmax and 0 < θ < π/2 (we use the symmetry about the equatorial plane θ = π/2). The radius rmax is chosen between 26r h and 38r h , depending on the black hole spin. The grid has 25 × 20 points. Its spacing is logarithmic in the r-direction and uniform in the θ-direction. We calculate Q α according to equations (8) and (11). The distributions f (p) and f (p ) are Lorentz-invariant (scalar) functions of four-momentum. They remain constant along the neutrino trajectories (eq. 6). Therefore, to calculate f (p) measured by ZAMO at a given point, it is sufficient to trace the neutrino with momentum p back to its emission point, find its four-momentum there in the rest-frame of the disc, pem, and use the equality f (p) = fem(pem). The values of fem(pem) and f em (p em ) are provided by the disc models described in § 2.1.
For every point of the grid we trace back the neutrino trajectories for 5000 directions uniformly distributed on the local ZAMO sky. Each direction is followed until the geodesic reaches the disc, goes into the black-hole horizon or reaches a maximum radius that we set equal to 140r h . For trajectories coming from outside 140r h or connecting to the horizon we set f = f = 0. For trajectories connecting our grid point to the disc, we find f and f as described above.
The accuracy of our calculation of Q α can be estimated by doubling the number of sampled geodesics; it is a few per cent. The accuracy of the volume-integrated quantityĖ is controlled mainly by the number of grid points. We checked it by repeating the same calculation with a coarser grid; the estimated error inĖ is smaller than 10 per cent (Table 1).

Comparison with previous works
We compared our calculations with three previous works: Birkl et al. (2007), Asano & Fukuyama (2001) and Popham et al. (1999). Birkl et al. (2007) calculated Q α νν and the corresponding volume-integrated energy deposition rateĖνν outside the ergosphere for several toy models of the neutrino source. They traced exactly the neutrino trajectories in the Kerr metric. Unfortunately, they incorrectly computedĖνν by integrating dP t instead of −dPt. For test purposes we did a similar integration. We computed their models D and REF, which assume that neutrinos are emitted by an isothermal blackbody ring (see Table 1 in Birkl et al. 2007). The results Figure 2. Color-coded and contour plot for log 10 Q t νν around an accretion disc withṀ = 1M /s. The black hole has mass M = 3M and spin parameter a = 0.95. The horizon sphere (black) has the radius r h ≈ 1.3GM/c 2 , and the inner edge of the disc (the marginally stable orbit) is at rms ≈ rg = 2GM/c 2 . Arrows show the projection of Q i νν /Q t νν on the plane of the figure. The white curve is where the radial component of injected momentum Q r νν changes sign. It roughly indicates the region where the deposited energy may be lost into the black hole rather than escape in an outflow. Disc Model A with viscosity parameter α = 0.1 (Section 2.1) was used in the calculation. Practically the same Q t νν is found for Models B and C, and for different α (e.g. α = 0.01).
agreed within the numerical errors of their and our calculations. Asano & Fukuyama (2001) calculatedĖνν on the rotational axis in Kerr metric as a function of r. The neutrino source was modeled as a blackbody disc with the temperature having a power-law dependence on radius. For test purposes, we repeated the calculation for their isothermal disc model and obtainedĖνν (r) on the axis. The functional shape ofĖνν (r) agrees with the shape of function G(r) in Asano & Fukuyama (2001) (see their eq. 18). The normalization of G(r) appears to be arbitrary, so we were unable to compare the numerical values ofĖνν . Popham et al. (1999) made more realistic assumptions about the neutrino source, similar to our Model C (Section 2.1). They evaluated numericallyĖνν for several values of accretion rateṀ and black-hole spin a (see their Table 3). The neutrino geodesics in Kerr metric were replaced by straight lines. When comparing their results with ours (which are presented in the next section) we found a significant disagreement, exceeding factor of 10. Generally, we find thatĖνν of Popham et al. (1999) was overestimated. The dependence ofĖνν on a that is suggested by Table 3 in Popham et al. (1999) is incorrect, apparently steeper than what we find numerically and estimate analytically (see next section).

RESULTS
The two processes of e ± creation considered in this paper (Section 2.3) give the total energy deposition rateĖ =Ėνν + EνB. We first focus onĖνν and discuss the contributionĖνB separately in Section 3.2.
3.1 Energy deposition from νν → e + e − Figures 2 and 3 show Q t νν for an accretion disc withṀ = 1M s −1 around a black hole of mass M = 3M . In Figure 2 the black hole is assumed to be rapidly rotating (a = 0.95), and in Figure 3 it is non-rotating (a = 0). In the case of a = 0.95, the deposition rate Q t νν peaks closer to the black hole and reaches much higher values, because the disc extends to smaller radii and emits a higher neutrino flux.
Volume integration of Q t νν and Q φ νν , as described by equation (16), gives the net energy deposition rate due to νν annihilation outside the black-hole horizon,Ėνν . Figure 4 showsĖνν as a function of the disc accretion rateṀ for the two cases, a = 0 and a = 0.95. For allṀ ,Ėνν is much higher when the black hole is rapidly rotating.
The results are sensitive toṀ and a, but remarkably insensitive to the details of the disc model.
(i) The uncertainty in the vertical structure of the accretion disc leads to a small uncertainty inĖνν as illustrated by two extreme models described in Section 2.1: Model A and Model B. The results of both models are well approximated by simple Model C (eq. 4). Deviations of Model C from Model A arise mainly where Model C does not accurately predict the neutrino flux from the disc, i.e. where the disc is not strongly cooled by neutrino emission.
(ii) The uncertainty in viscosity parameter α ∼ 0.1 has almost no effect onĖνν as long asṀign <Ṁ <Ṁtrap. In particular, Model C in this range ofṀ is explicitly inde- pendent of α. The two characteristic accretion ratesṀign andṀtrap depend on viscosity parameter α (see eq. 5); in Figure 4 we assumed α = 0.1.
The dependence ofĖνν on the black hole spin, for a fixedṀ = 1M /s, is shown in Figure 6. Instead of using the spin parameter a directly, it is more instructive to ploṫ Eνν versus radius of the last (marginally stable) orbit rms. Then one can see the power-law dependence ofĖνν on rms: Eνν ∝ r −4.8 ms . This power-law is accurate for rms > rg which corresponds to a < 0.9. For rms < rg,Ėνν has a somewhat stronger dependence on rms. The standard relation between rms and a (e.g. Page & Thorne 1974) is shown in Figure 5. For non-rotating black holes rms = 6GM/c 2 and for maximally rotating black holes rms = GM/c 2 .
3.2 Energy deposition from ν → νe + e − in a strong magnetic field The four-momentum deposition due to reaction ν → νe + e − is given by equation (11). It depends on the magnetic field. We will assume the strongest field that could be expected, where Pmax is the maximum pressure in the disc; we find this pressure using the numerical disc model of CB07. For illustration, consider a toy model where the magnetic field is perpendicular to the disc, uniform in the region r < 25r h and zero outside this region. The model overestimates a realistic magnetic field around an accretion disc. We will show that even such a strong field gives a modestĖνB. The maximum pressure Pmax and B = (8πPmax) 1/2 depend onṀ and a.
The magnetic field ranges in our models from ∼ 6 × 10 13 G to ∼ 8 × 10 15 G, increasing withṀ and a. Figure 6 showsĖνB as a function of rms for fixedṀ = 1M /s. It approximately follows the same scaling relation as we found forĖνν (rms),ĖνB ∝ r −4.8 ms . Figure 7 shows the ratioĖνB/Ėνν as a function ofṀ . This ratio is small for accretion rates aboveṀign and varies slowly withṀ or a. Thus, we find thatĖνB makes a small contribution to the totalĖ. This conclusion is valid for the most interesting range of accretion ratesṀ >Ṁign. ForṀ <Ṁign, the neutrino flux quickly decreases, which strongly suppresses the νν annihilation above the disc. The reduction in reaction ν → νe + e − is less severe andĖνB ex-ceedsĖνν . As a result,ĖνB dominates the energy deposition rateĖ whenṀ is belowṀign. The dependence of the energy deposition rate onṀ is summarized in Figure 8, a modified version of Figure 4. It approximately represents the totalĖ by showingĖνν aṫ M >Ṁign andĖνB atṀ <Ṁign.

Scaling ofĖνν andĖνB withṀ , M and rms
TheĖνν dependence onṀ , M and a may be estimated using simple analytical arguments (Beloborodov 2008). The total rate of νν annihilation around the discṄνν [s −1 ] is proportional to the volume of the main annihilation region r 3 , which scales as r 3 ms , and to the typical local annihilation rate in this region,ṅνν ∼ σcnν nν . Here nν and nν are the number densities of neutrinos and anti-neutrinos and σ is the annihilation cross-section. The cross-section scales with the energies of the annihilating ν andν approximately as σ ∝ Eν Eν , which leads tȯ The energy fluxes Fν and Fν from an efficiently cooled disc (Ṁign <Ṁ <Ṁtrap) carry the released gravitational energy and do not depend on details of the disc structure. The fluxes scale with M ,Ṁ , and r as Fν ∼ Fν ∝ MṀ /r 3 . This giveṡ   Fig. 4 but now including the contribution of reaction ν → νe + e − . This contributionĖ νB is shown only aṫ M <Ṁ ign , where it is important;Ė νB is small compared witḣ Eνν atṀ >Ṁ ign (cf. Fig. 7).
Here we assumed that the size of the main annihilation region r is proportional to rms = xmsrg, where xms is determined by the black-hole spin parameter a (see Fig. 5).

CONCLUSIONS
We have performed detailed numerical calculations of e ± creation around hyper-accreting spinning black holes. We studied two reactions: νν → e + e − and ν → νe + e − (in a strong magnetic field). The reaction of νν annihilation dominates the energy deposition rate around discs witḣ M >Ṁign, which are strong emitters of ν andν. We found that the net energy deposition rate due to this process,Ėνν , is well approximated by a simple formula (see Fig. 4 whereṁ =Ṁ /M s −1 , xms = rms(a)/rg, rg = 2GM/c 2 , andṀtrap,Ṁign are given in equation (5). The dependence ofĖνν on the black hole spin is huge: x −4.8 ms varies by a factor of 200 for 0 < a < 0.95. Note that α (viscosity parameter of the disc) enters the result only throughṀign andṀtrap. Our numerical simulations in this paper are limited to black holes with mass M = 3M , and theĖνν dependence on M is evaluated analytically.
The efficiency of νν annihilation can be defined as = Eνν /L where L is the total neutrino luminosity of the disc. For example, a = 0.95 (which corresponds to xms = 0.97) gives L ≈ 0.15Ṁ c 2 (CB07) and ≈ 0.05ṁ 5/4 forṀign < M <Ṁtrap. Note thatĖνν is defined in this paper as the total energy deposition rate outside the event horizon. A fraction of the created e ± plasma falls into the black hole and does not contribute to the observed explosion (Fig. 2). The corresponding refinement of depends on the plasma dynamics outside the disc, which is affected by magnetic fields and hard to calculate without additional assumptions.
The obtainedĖνν may be comparable to GRB luminosities L obs . A plausible typical value for L obs is ∼ 10 51 erg/s (it depends on the beaming angle of the observed explosion, which is usually hard to estimate from available data). This power is easily supplied by neutrino heating if the black hole has a large spin, e.g. a = 0.95, for a moderate accretion ratė M > ∼ 0.3M s −1 . For a non-rotating black hole, this mechanism of GRB explosion requires ∼ 10 times higher accretion rates.