ULTRA-RELATIVISTIC, NEUTRINO-DRIVEN FLOWS IN GAMMA-RAY BURSTS: A DOUBLE TRANSONIC FLOW SOLUTION IN A SCHWARZSCHILD SPACETIME

The relativistic outflows producing the prompt and afterglow emissions in gamma-ray bursts (GRBs) are commonly thought to be powered by hyperaccreting black holes. The high Lorentz factors inferred from energy considerations, compactness arguments, and afterglow models, Γ ∼ 10²–10³, require extremely low baryon load at the outflow injection point, which seems difficult to achieve in disk outflows (Levinson 2006; Barzilay & Levinson 2008; Metzger et al. 2008). A natural way to avoid baryonic contamination is to accelerate the outflow in the region above the horizon of the black hole (Levinson & Eichler 1993). Two popular jet production mechanisms, which are widely discussed in the literature, are magnetic extraction of the rotational energy of a Kerr black hole (Blandford & Znajek 1977) and pair creation on horizon threading field lines by annihilation of neutrinos that emanate from the accretion disk surrounding the black hole (Eichler et al. 1989; Jaroszynski 1996; Popham et al. 1999; Asano & Fukuyama 2001; Birkl et al. 2007; Zalamea & Beloborodov 2011, hereafter ZB11). The former mechanism requires nearly maximal rotation of the black hole and sufficiently high magnetization. The latter mechanism can operate also in a Schwarzschild spacetime, although extremely high accretion rates (≳ 3 M_☉ s⁻¹) are required in this case to account for the characteristic GRB luminosities (ZB11). It is worth noting that the presence of a rapidly spinning Kerr black hole may not be sufficient to guarantee effective production of a Blandford–Znajek (hereafter BZ) outflow; overloading of magnetic field lines by relativistic pairs may lead to a complete shutdown of this mechanism, in which case one has to rely on the neutrino-driven outflow to power a GRB. However, the accretion rate required to power a typical GRB by the latter mechanism can be considerably reduced when the specific angular momentum of the black hole, a, is high enough. The reason is that the neutrino luminosity emitted by the disk and, hence, the rate at which energy is deposited above the horizon by $\nu \bar{\nu }$ annihilation, increases sharply with increasing a (ZB11).

A question of interest is what fraction of the total energy deposited above the horizon via $\nu \bar{\nu }$ annihilation emerges at infinity. To address this question, we constructed a model for the double-transonic flow established in the polar region, which incorporates energy injection by the external neutrino source in a self-consistent manner.

Neutrino annihilation in the baryon-free region above the horizon generates a relativistically hot fluid consisting of a mixture of e^± pairs and radiation in equilibrium. Part of this fluid is accelerated outward by pressure gradient forces, and the other part falls into the black hole, owing to the strong gravitational force exerted on it. In what follows, we consider the structure of the double transonic flow thereby produced. Schematic illustration of the model is shown in Figure 1. We restrict our analysis to the Schwarzschild spacetime and suppose that the flow is stationary, radial, unmagnetized, and non-rotating. For clarity, we consider baryon-free plasma, although our analysis applies also when baryons are present, provided the enthalpy per baryon is sufficiently large (w/n_b ≫ 1). The idea is to find a self-consistent solution that starts from a stagnation point, where the flow velocity vanishes, and crosses two sonic points; an inner one located in the inflow section below the stagnation point and an outer one located in the outflow section above the stagnation point. As shown below, for a given choice of an energy deposition profile the stagnation radius and stagnation pressure are uniquely determined by the requirement that the solution pass smoothly through both sonic points. A similar analysis has been performed by Jaroszynski (1996) for specific conditions. Here, we investigate the behavior of the solution for a range of conditions and also compute the entropy profile in the outflow. Henceforth, we use geometrical units (c = G = 1), unless otherwise specified.

Zoom In Zoom Out Reset image size

2.1. Energy Deposition Rate

The annihilation of neutrinos emitted from a hyperaccretion disk was calculated in a number of recent works, under different simplifying assumptions (e.g., Popham et al. 1999; Asano & Fukuyama 2001; Birkl et al. 2007; ZB11). The most detailed analysis is presented in ZB11 who adopted the relativistic disk model of Chen & Beloborodov (2007) as the neutrino source and employed a geodesic-tracing method in the Kerr spacetime to evaluate the local energy–momentum deposition rates by the reaction $\nu \bar{\nu }\rightarrow$ e⁺e⁻ in the vicinity of the black hole.

The local energy–momentum deposition rates, $Q_{\nu \bar{\nu }}^\alpha$, were computed in ZB11 in the frame of a zero angular momentum observer (ZAMO). Quite generally, these rates are functions of r and θ, when expressed in Boyer–Lindquist coordinates, x^u = (t, r, θ, ϕ). Since we restrict our analysis to a conical flow, we shall adopt a local energy deposition rate of the form:

$$\begin{equation} Q^t_{\nu \bar{\nu }}(r)=\dot{Q}_0f(r/r_H), \end{equation} \tag{ 1 }$$

where $Q^t_{\nu \bar{\nu }}(r)$ represents the angle-averaged rate at radius r, r_H is the horizon scale, and f(1) = 1. From Figures 2 and 3 in ZB11, we estimate f(x) ≃ x^−4.5 for a = 0.95 and f(x) ≃ x^−3.5 for a = 0, in the range r_H ⩽ r ⩽ 30m delineated in the figures, where m is the black hole mass in geometrical units and a its specific angular momentum. The profile should steepen as the radius increases, approaching f(x)∝x⁻⁸ at radii much larger than the size R_ν of the neutrino source (e.g., Goodman et al. 1987; Qian & Woosley 1996). For the disk model outlined in Chen & Beloborodov (2007), R_ν can be identified with the ignition radius of the disk, r_ign, within which neutrino emission is efficient.

The net energy deposition rate, as measured by a distant observer, is given in terms of the ZAMO rates $Q_{\nu \bar{\nu }}^\alpha$ by ZB11,

$$\begin{equation} \dot{E}_{\nu \bar{\nu }}=\int _{r\ge r_H}{\left(\frac{Q^t_{\nu \bar{\nu }}}{\sqrt{-g^{tt}}}-\frac{g_{t\phi }}{\sqrt{g_{\phi \phi }}}Q_{\nu \bar{\nu }}^\phi \right)\sqrt{-g}dr d\theta d\phi }, \end{equation} \tag{ 2 }$$

here g_μν are the metric components of the Kerr spacetime in Boyer–Lindquist coordinates. The results exhibited in Figure (4) of ZB11 (see also their Equation (22) for a fitting formula) indicate that $\dot{E}_{\nu \bar{\nu }}$ depends sensitively on the accretion rate $\dot{m}$ (given in units of M_☉ s⁻¹) and the spin of the black hole a, with $\dot{E}_{\nu \bar{\nu }}\simeq 5\times 10^{49}\dot{m}^{9/4}$ erg s⁻¹at a = 0 and $\dot{E}_{\nu \bar{\nu }}\simeq 10^{52}\dot{m}^{9/4}$ erg s⁻¹ at a = 0.95, for a black hole mass M_BH = 3 M_☉, and accretion rates in the range $\dot{M}_{{\rm ign}}<\dot{M}<\dot{M}_{{\rm trap}}$. At $\dot{M}<\dot{M}_{{\rm ign}}$ neutrino emission is severely suppressed by virtue of the low disk temperature. At $\dot{M}_{{\rm trap}}<\dot{M}$, the neutrinos are trapped in the disk and advected into the black hole, whereby the neutrino luminosity saturates.

2.2. Flow Equations

The stress–energy tensor of a purely hydrodynamic flow takes the form

$$\begin{equation} T^{\alpha \beta } = wu^{\alpha }u^{\beta } + pg^{\alpha \beta }, \end{equation} \tag{ 3 }$$

here u^α is the four-velocity measured in units of c, w is the specific enthalpy, p is the pressure, and g_μν = diag(− α², α⁻², r², r²sin ²θ), α² = 1–2m/r, is the metric tensor of the Schwarzschild spacetime.

The dynamics of the flow is governed by the energy–momentum equations:

$$\begin{equation} \frac{1}{\sqrt{-g}}(\sqrt{-g}T^{\alpha \beta })_{,\alpha }+ \Gamma ^{\beta }_{\ \mu \nu }T^{\mu \nu } = q^\beta, \end{equation} \tag{ 4 }$$

where q^β denotes the source terms associated with energy–momentum transfer by the external agent and $\Gamma ^\beta _{\mu \nu }$ denotes the affine connection. By contracting g_βγ with Equation (4), using the relation $(\sqrt{-g}g^{\alpha \beta })_{,\alpha }+\sqrt{-g} \Gamma ^\beta _{\mu \nu }g^{\mu \nu }=0$, taking the t component, and noting that Γ_μtνu^μu^ν = 0 for a stationary flow, one obtains

$$\begin{equation} \frac{1}{\sqrt{-g}}(\sqrt{-g}w u^{\alpha } u_t)_{,\alpha }=q_t. \end{equation} \tag{ 5 }$$

Likewise, contracting u_β with Equation (4), using the identity $\Gamma ^{\beta }_{\ \mu \nu }u_\beta u^\mu u^\nu =-u_\beta u^\nu (u^\beta)_{,\nu }$, yields

$$\begin{equation} \frac{1}{\sqrt{-g}}(\sqrt{-g}w u^{\alpha })_{,\alpha }-u^\alpha p_{,\alpha }=-u_\nu q^\nu. \end{equation} \tag{ 6 }$$

Combining Equation (6) with the thermodynamic identity d(wV) − Vdp = kTd(Vs), where V is the volume of a fluid element, T is its temperature, and and s is its specific entropy, yields the change in s:

$$\begin{equation} \frac{kT}{\sqrt{-g}}(\sqrt{-g}s u^{\alpha })_{,\alpha }=-u_\nu q^\nu. \end{equation} \tag{ 7 }$$

Henceforth, we consider a conical flow, for which u^μ = (u^t, u^r, 0, 0), u^μ∂_μ = u^r∂_r, and denote by u_p = u^r/α the poloidal velocity, and by $\gamma =(1+u_p^2)^{1/2}$ and v = u_p/γ the corresponding Lorentz factor and three-velocity. The convention here is that u_p < 0 on inflow lines and u_p > 0 on outflow lines. Since the created pairs are relativistic, we adopt the equation of state w = 4p. Furthermore, we neglect momentum transfer to the flow by the created pairs so that q^β = (q^t, 0, 0, 0). The source term q^t is measured in the frame of a distant observer and can be expressed in terms of the local energy deposition rate adopted in Equation (1) as $cq^t=\sqrt{-g^{tt}}Q^t_{\nu \bar{\nu }}=\dot{Q}_0f(x)/\alpha$. Under the above simplifications, Equations (5) and (6) reduce to

$$\begin{eqnarray} &&(3v^2-1)\partial _x\ln u_p=\frac{(3-4\gamma ^2)}{\gamma u_p 4\tilde{p}}\frac{f(x)}{(1-x^{-1})^{1/2}}+\frac{2}{x}-\frac{1}{x(x-1)}\nonumber\\ \end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray} \partial _x\ln \tilde{p}&=&\frac{f(x)}{\gamma u_p 4\tilde{p}(1-x^{-1})^{1/2}}-\frac{2}{x}-\frac{1}{x(x-1)}\nonumber\\ &&-(1+v^2)\partial _x\ln u_p, \end{eqnarray} \tag{ 9 }$$

in terms of the dimensionless radius x = r/r_s, r_s = 2m, and the normalized pressure $\tilde{p}=p/(\dot{Q}_0t_d)$, with t_d = r_s/c being the dynamical time of the black hole. Equation (8) has critical points at v = ±c_s, where $c_s=1/\sqrt{3}$ is the sound speed.

At the stagnation point x = x_st, where u_p = 0, the above equations yield

$$\begin{eqnarray} &&\partial _x u_p=\frac{f(x_{{\rm st}})}{4 \tilde{p}_{{\rm st}}\big(1-x_{{\rm st}}^{-1}\big)^{1/2}}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &&\partial _x\ln (\tilde{p})=-\frac{2}{x_{{\rm st}}(x_{{\rm st}}-1)}, \end{eqnarray} \tag{ 11 }$$

where $\tilde{p}_{{\rm st}}=\tilde{p}(x_{{\rm st}})$ is the stagnation pressure. Thus, for a given choice of f(x) the solution is fully determined once x_st and p_st are known.

The regularity conditions at the sonic points, obtained from Equation (8), read

$$\begin{equation} \sqrt{3}x^2_{c1}\big(1-x^{-1}_{c1}\big)^{1/2}f(x_{c1})=-2(2x_{c1}-3)\tilde{p}_{c1} \end{equation} \tag{ 12 }$$

and

$$\begin{equation} \sqrt{3}x^2_{c2}\big(1-x^{-1}_{c2}\big)^{1/2}f(x_{c2})=2(2x_{c2}-3)\tilde{p}_{c2}, \end{equation} \tag{ 13 }$$

denoting the sonic point of the inflow (outflow) by x_c1(x_c2) and noting that $u_p= -1/\sqrt{2}$ at x_c1 and $u_p= 1/\sqrt{2}$ at x_c2. Evidently, x_c1 < 3/2 and x_c2 > 3/2. The existence of two sonic points is a consequence of energy injection. When f(x) = 0, the system has only one critical point at x_c = 3/2.

2.3. Asymptotic Power and Entropy of the Outflow

Well above the sonic radius, at x ≫ x_c2, the power of the outflow is given by $L_{j\infty }=2\pi c\int _0^{\theta _0}{T^{0r}r^2\sin \theta d\theta }=2\pi (1-\cos \theta _0)w\gamma u_p$, here θ₀ is the opening angle of the flow. Integrating Equation (5) from the stagnation point x_st to infinity, using the boundary condition u_p(x_st) = 0, and noting that q_t = g_ttq^t = −α²q^t, one obtains

$$\begin{eqnarray} L_{j\infty }&=&2\pi (1-\cos \theta _0)\int _{r_{{\rm st}}}^\infty {cq_t r^2dr}\nonumber\\ &=&2\pi (1-\cos \theta _0)\dot{Q}_0r_s^3\int _{x_{{\rm st}}}^\infty {(1-x^{-1})^{1/2}f(x) x^{2}dx}.\nonumber\\ \end{eqnarray} \tag{ 14 }$$

Let us define the outflow production efficiency as the fraction of $\dot{E}_{\nu \bar{\nu }}$ that emerges at infinity. From Equation (2) with g_tϕ = 0, $\sqrt{-g^{tt}}=\alpha ^{-1}$, and Equation (14), one finds

$$\begin{equation} \epsilon \equiv \frac{L_{j\infty }}{\dot{E}_{\nu \bar{\nu }}}=\frac{\int _{x_{{\rm st}}}^\infty {(1-x^{-1})^{1/2}f(x) x^{2}dx}}{{\int _{1}^\infty {(1-x^{-1})^{1/2}f(x) x^{2}dx}}}. \end{equation} \tag{ 15 }$$

The rate at which entropy is carried by the flow at some radius r is obtained by integrating Equation (7) from the stagnation radius r_st to r > r_st:

$$\begin{eqnarray} \frac{dS(r)}{dt} &\equiv& 2\pi (1-\cos \theta _0) r^2cu^r s\nonumber\\ &=& 2\pi (1-\cos \theta _0)\frac{\dot{Q}_0r_s^3}{kT_0}\int _{x_{{\rm st}}}^x \gamma \tilde{p}^{-1/4}f(x^\prime)x^{\prime 2} dx^\prime,\nonumber\\ \end{eqnarray} \tag{ 16 }$$

where T₀ is a fiducial temperature defined by¹ $T_0=(12\dot{Q}_0t_d/11a)^{1/4}$. If the outflow also carries baryons at a rate $\dot{N}_b=2\pi (1-\cos \theta _0)n_bu^rr^2$, for a baryon density n_b, then the dimensionless entropy per baryon can be expressed as

$$\begin{equation} \sigma (x)=\frac{1}{\dot{N}_b}\frac{dS}{dt}=\frac{m_pc^2\gamma _\infty }{kT_0}F(x), \end{equation} \tag{ 17 }$$

with $\gamma _{\infty }=L_{j\infty }/(\dot{N}_bm_pc^2)$ being the maximum Lorentz factor of the outflow, and

$$\begin{equation} F(x)=\frac{\int _{x_{{\rm st}}}^x \gamma \tilde{p}^{-1/4}f(x^\prime)x^{\prime 2} dx^{\prime }}{\int _{x_{{\rm st}}}^\infty {(1-x^{\prime -1})^{1/2}f(x^\prime) x^{\prime 2}dx^\prime }}. \end{equation} \tag{ 18 }$$

Equations (8) and (9) have been integrated numerically using the energy deposition profile of the form f(x) = x^−b. The pressure at the inner sonic point, $\tilde{p}_{c1}$, is used as a free parameter that we adjust to converge to the desired solution. The integration starts at the inner sonic point x_c1, where $u_p(x_{c1})=-1/\sqrt{2}$, and is repeated iteratively by changing the value of $\tilde{p}_{c1}$ until a smooth transition across the outer sonic point is obtained. The value of x_c1 is computed, in each run, from the regularity condition (12). The stagnation radius x_st is then found numerically from the condition u_p(x_st) = 0.

Solutions are exhibited in Figure 2 for different values of b. The horizontal dashed lines in the left panel mark the sonic velocities of the inflow ($u_p=-1/\sqrt{2}$) and the outflow ($u_p=1/\sqrt{2}$). As is seen, in all cases the sonic point of the outflow is located well above that of an adiabatic flow (x_c = 3/2). This is a consequence of the injection of energy over extended scales, which loads the flow and delays its acceleration. The values of the stagnation radii found from the integration are x_st = 1.57, 1.70, 1.81, and 1.95 for b = 5, 4, 3.5, and 3, respectively. For comparison, note that the radius at which the escape velocity, $v_{{\rm esc}}=\sqrt{r_s/r}=x^{-1/2}$, equals the sound speed, $c_s=1/\sqrt{3}$, is x = 3, slightly above the stagnation point. The corresponding efficiencies, computed using Equation (15), are = 0.58, 0.73, and 0.83 for b = 5, 4, and 3.5, respectively. For b = 3 the energy deposition rate diverges logarithmically and needs to be cut off at some radius x_ν. We arbitrarily invoked x_ν = 100, whereby = 0.93.

Zoom In Zoom Out Reset image size

By employing Equations (14) and (17), we can express the entropy per baryon in the form:

$$\begin{eqnarray} \sigma (x)&=&1.3\times 10^5\eta _b\left(\frac{\gamma _\infty }{300}\right)\left(\frac{L_{{\rm iso}}}{{1}0^{51}\, {\rm erg\,s^{-1}}}\right)^{-1/4}\nonumber\\ &&\times \left(\frac{M_{{\rm BH}}}{3\, M_\odot }\right)^{1/2}F(x), \end{eqnarray} \tag{ 19 }$$

where L_iso = 2L_j∞/(1 − cos θ₀) is the isotropic equivalent luminosity, M_BH is the mass of the black hole, and η_b = 1.06, 0.84, and 0.62 for b = 3.5, 4, and 5, respectively. A plot of F(x) is shown in Figure 3. The asymptotic value of F(x) represents enhancement of the entropy relative to that produced in an adiabatic outflow injected from a radius R₀ ∼ r_s.

Zoom In Zoom Out Reset image size

3.1. Asymptotic Behavior

We have computed the structure of a double transonic flow generated above the horizon of a Schwarzschild black hole by annihilation of neutrinos emitted from a hyperaccretion disk. We have shown that, for a given choice of the energy deposition profile, there exists a unique solution that passes smoothly through the inner and outer sonic points. The stagnation point was found to be located slightly below the escape radius of sonic material—in the range 1.5–2 Schwarzschild radii, depending on the energy deposition profile. The asymptotic power of the outflow is given as the integral of the energy deposition rate above the stagnation radius. The outflow production efficiency was found to be typically large, with over one half of the injected energy emerging at infinity. Our analysis ignores collimation of the outflow, which seems to be required by the small opening angles inferred in many sources. This will, most likely, lead to the formation of recollimation shocks, as seen in numerical simulations (Aloy et al. 2000, 2005), that may alter the dynamics of the outflow. The outflow production efficiency should not be significantly affected by these details.

The continuous injection of energy pushes out the outer sonic point and delays the acceleration of the outflow above the sonic point. The linear acceleration phase commences once the energy deposition rate becomes small and the outflow enters the adiabatic regime. The radius at which this happens depends on the energy deposition profile and may exceed 100m at large accretion rates. The deposition of energy over extended scales leads to generation of a relatively large specific entropy. As a consequence, in models wherein the prompt emission originates from the photosphere (e.g., Ryde & Pe'er 2009; Lazzati et al. 2009; Pe'er et al. 2012; Levinson 2012; Beloborodov 2013), the thermal peak may be located at an energy lower than previously estimated. For instance, if the prompt emission is produced by sub-photospheric, radiation-mediated shocks (Bromberg et al. 2011; Levinson 2012), then the observed temperature behind the shock is

$$\begin{eqnarray*} kT_{{\rm obs}}&\simeq& 50 \left(\frac{F_\infty }{10}\right)^{-1}\left(\frac{L_{{\rm iso}}}{{1}0^{51}\,{\rm erg\,s^{-1}}}\right)^{1/4}\nonumber\\ &&\times \left(\frac{M_{{\rm BH}}}{M_\odot }\right)^{-1/2}u_s\quad {\rm keV}, \end{eqnarray*}$$

where u_s is the shock four-velocity (i.e., the four-velocity of the upstream fluid, as measured in the shock frame) and F_∞ is the asymptotic value of the function F(x) exhibited in Figure 3. Consequently, the peak of the spectrum produced by mildly relativistic shocks (u_s ≳ 1) can be located at an energy E_peak ∼ 3kT_obs ∼ 200 keV, significantly lower than those exhibited in Levinson (2012). Additional dissipation and entropy generation may result from the interaction of the outflow with the stellar envelope, via the formation of re-collimation shocks (Bromberg & Levinson 2007).

In situations where the central black hole is rapidly rotating, the outflow can be produced by the BZ mechanism. The resultant Poynting power is expected to exceed the net power deposited above the horizon by the reaction $\nu \bar{\nu }\rightarrow e^+ e^-$ (e.g., Kawanaka et al. 2013). However, the activation of the BZ mechanism requires sufficiently high magnetization near the horizon, and it is therefore anticipated that overloading of magnetic field lines by relativistic pairs may lead to a complete shutdown of this process. In the latter case, the outflow will be driven by the pressure of the injected plasma, as described in Section 2, rather than by magnetic extraction of the spin energy of the black hole. Interestingly, relatively low accretion rates may be favorable for producing powerful outflows. The critical load above which the activation of the BZ process is prevented is a key issue, currently under investigation (N. Globus & A. Levinson, in preparation).

We thank Miguel Aloy for useful comments.