Neutrino-dominated accretion flows as the central engine of gamma-ray bursts

TL;DR

This paper reviews the theoretical understanding of neutrino-dominated accretion flows around black holes and their role as the central engine in gamma-ray bursts, covering structure, jet production, and observational signatures.

Contribution

It provides a comprehensive review of NDAF properties, their application to GRB phenomenology, and the potential for multi-messenger observations to probe these accretion flows.

Findings

NDAFs can produce the high luminosities observed in GRBs.

Jet power depends on black hole parameters and accretion rates.

Multi-messenger signals can help identify NDAF activity in GRBs.

Abstract

Neutrino-dominated accretion flows (NDAFs) around rotating stellar-mass black holes (BHs) are plausible candidates for the central engines of gamma-ray bursts (GRBs). NDAFs are hyperaccretion disks with accretion rates in the range of around 0.001-10 $M_\odot ~\rm s^{-1}$, which have high density and temperature and therefore are extremely optically thick and geometrically slim or even thick. We review the theoretical progresses in studying the properties of NDAFs as well as their applications to the GRB phenomenology. The topics include: the steady radial and vertical structure of NDAFs and the implications for calculating neutrino luminosity and annihilation luminosity, jet power due to neutrino-antineutrino annihilation and Blandford-Znajek mechanism and their dependences on parameters such as BH mass, spin, and accretion rate, time evolution of NDAFs, effect of magnetic fields, applications of NDAF theories to the GRB phenomenology such as lightcurve variability, extended emission, X-ray flares, kilonovae, etc., as well as probing NDAFs using multi-messenger signals such as MeV neutrinos and gravitational waves.

Equations

$$\dot{M}=-4\pi \rho H{\Delta }^{1/2}\frac{v_r}{\sqrt{1-{v_r}^2}},$$

$$\frac{v_r}{1-{v_r}^2}\frac{d{v}_r}{dr}=\frac{\mathcal{A}}{r}-(1-{v_r}^2)\frac{1}{\lambda \rho }\frac{dp}{dr},$$

$$\mathcal{A}\equiv -\frac{A}{r^3\Delta {\Omega }_K^{+}{\Omega }_K^{-}}\frac{(\Omega -{\Omega }_K^{+})(\Omega -{\Omega }_K^{-})}{1-{\tilde{\Omega }}^2{\tilde{R}}^2},$$

$$\dot{M}(\mathcal{L}-{\mathcal{L}}_{\mathrm{in}})=\frac{4\pi pH{A}^{1/2}{\Delta }^{1/2}\gamma }{r},$$

$$\gamma =\sqrt{\frac{1}{1-{v_r}^2}+\frac{{\mathcal{L}}^2{r}^2}{A}}.$$

$$\frac{p}{\lambda \rho H^2}=\frac{{\mathcal{L}}^2-a^2({ϵ}^2-1)}{r^4},$$

$$ϵ=-\gamma \frac{r{\Delta }^{1/2}}{A^{1/2}}-\frac{2ar}{A}\mathcal{L}.$$

$$-\frac{\dot{M}}{2\pi r^2}\left(\frac{u}{\rho }\frac{d\ln u}{d\ln r}-\frac{p}{\rho }\frac{d\ln \rho }{d\ln r}\right)=-\frac{2\alpha pHA{\gamma }^2}{r^3}\frac{d\Omega }{dr}-2Q^{-},$$

$${\tau }_{{\nu }_i}={\tau }_{s,{\nu }_i}+{\tau }_{a,{\nu }_i},$$

$${\tau }_{s,{\nu }_i}\approx H({\sigma }_{\mathrm{e},{\nu }_i}{n}_{\mathrm{e}}+\sum _j{\sigma }_{j,{\nu }_i}{n}_j),$$

$${\sigma }_{\mathrm{e},{\nu }_i}\approx \frac{3k_{\mathrm{B}}T{\sigma }_0{e}_{{\nu }_i}}{8m_{\mathrm{e}}{c}^2}(1+\frac{{\eta }_{\mathrm{e}}}{4})[(C_{V,{\nu }_i}+C_{A,{\nu }_i}{)}^2$$

$$+\frac{1}{3}(C_{V,{\nu }_i}-C_{A,{\nu }_i}{)}^2],$$

$${\sigma }_{{\mathrm{n}}_1,{\nu }_i}\approx \frac{{\sigma }_0{e}_{{\nu }_i}^2}{4}[(C_{V,{\nu }_i}-1{)}^2+3g_A^2(C_{A,{\nu }_i}-1{)}^2],$$

$${\sigma }_{{\mathrm{n}}_2,{\nu }_i}\approx \frac{{\sigma }_0{e}_{{\nu }_i}^2}{4}\frac{1+3g_A^2}{4},$$

$${\sigma }_{n_j,{\nu }_i}\approx \frac{{\sigma }_0}{16}{e}_{{\nu }_i}^2(Z_j+N_j)[1-\frac{2Z_j}{Z_j+N_j}(1-2{\sin }^2{\theta }_W){]}^2,$$

$$n_{{\mathrm{e}}^{\mp }}=\frac{1}{{\hslash }^3{\pi }^2}{\int }_0^{\infty }dp{p}^2\frac{1}{e^{(\sqrt{p^2{c}^2+{m_{\mathrm{e}}}^2{c}^4}\mp {\mu }_{\mathrm{e}})/k_{\mathrm{B}}T}+1},$$

$${\tau }_{a,{\nu }_i}=\frac{q_{{\nu }_i}H}{4(7/8)\sigma T^4},$$

$$q_{{\nu }_i}=q_{\mathrm{Urca}}+q_{{\mathrm{e}}^{-}+{\mathrm{e}}^{+}\to {\nu }_i+{\overline{\nu }}_i}+q_{\mathrm{n}+\mathrm{n}\to \mathrm{n}+\mathrm{n}+{\nu }_i+{\overline{\nu }}_i}+q_{\tilde{\gamma }\to {\nu }_i+{\overline{\nu }}_i}.$$

$$q_{\mathrm{Urca}}=q_{\mathrm{p}+{\mathrm{e}}^{-}\to \mathrm{n}+{\nu }_{\mathrm{e}}}+q_{\mathrm{n}+{\mathrm{e}}^{+}\to \mathrm{p}+{\overline{\nu }}_{\mathrm{e}}}+q_{\mathrm{n}\to \mathrm{p}+{\mathrm{e}}^{-}+{\overline{\nu }}_{\mathrm{e}}}$$

$$+q_{X_j+{\mathrm{e}}^{-}\to X_j^{\prime }+{\nu }_{\mathrm{e}}},$$

$$q_{\mathrm{p}+{\mathrm{e}}^{-}\to \mathrm{n}+{\nu }_{\mathrm{e}}}=\frac{G_F^2{\cos }^2{\theta }_c}{2{\pi }^2{\hslash }^3{c}^2}(1+3g_A^2)n_1$$

$$\times {\int }_Q^{\infty }d{E}_{\mathrm{e}}{E}_{\mathrm{e}}\sqrt{{E_{\mathrm{e}}}^2-{m_{\mathrm{e}}}^2{c}^4}(E_{\mathrm{e}}-Q{)}^3{f}_{{\mathrm{e}}^{-}},$$

$$q_{\mathrm{n}+{\mathrm{e}}^{+}\to \mathrm{p}+{\overline{\nu }}_{\mathrm{e}}}=\frac{G_F^2{\cos }^2{\theta }_c}{2{\pi }^2{\hslash }^3{c}^2}(1+3g_A^2)n_2$$

$$\times {\int }_{m_{\mathrm{e}}{c}^2}^{\infty }d{E}_{\mathrm{e}}{E}_{\mathrm{e}}\sqrt{{E_{\mathrm{e}}}^2-{m_{\mathrm{e}}}^2{c}^4}(E_{\mathrm{e}}+Q{)}^3{f}_{{\mathrm{e}}^{+}},$$

$$q_{\mathrm{n}\to \mathrm{p}+{\mathrm{e}}^{-}+{\overline{\nu }}_{\mathrm{e}}}=\frac{G_F^2{\cos }^2{\theta }_c}{2{\pi }^2{\hslash }^3{c}^2}(1+3g_A^2)n_2$$

$$\times {\int }_{m_{\mathrm{e}}{c}^2}^{Q}d{E}_{\mathrm{e}}{E}_{\mathrm{e}}\sqrt{{E_{\mathrm{e}}}^2-{m_{\mathrm{e}}}^2{c}^4}(Q-E_{\mathrm{e}}{)}^3(1-f_{{\mathrm{e}}^{-}}),$$

$$q_{X_j+{\mathrm{e}}^{-}\to X_j^{\prime }+{\nu }_{\mathrm{e}}}=\frac{G_F^2{\cos }^2{\theta }_c}{2{\pi }^2{\hslash }^3{c}^2}{g}_A^2\frac{2}{7}{N}_p(Z_j)N_h(N_j)n_j$$

$$\times {\int }_{Q^{\prime }}^{\infty }d{E}_{\mathrm{e}}{E}_{\mathrm{e}}\sqrt{{E_{\mathrm{e}}}^2-{m_{\mathrm{e}}}^2{c}^4}(E_{\mathrm{e}}-Q^{\prime }{)}^3{f}_{{\mathrm{e}}^{-}},$$

$$N_p(Z_j)=\{\begin{array}{l}0,Z_j<20,\\ Z_j-20,20<Z_j<28,\\ 8,Z_j>28,\end{array}$$

$$N_h(N_j)=\{\begin{array}{l}6,N_j<34,\\ 40-N_j,34<N_j<40,\\ 0,N_j>40.\end{array}$$

$$q_{{\mathrm{e}}^{-}+{\mathrm{e}}^{+}\to {\nu }_i+{\overline{\nu }}_i}=\frac{Q_c}{36\pi }\{(C_{V,{\nu }_i}^2+C_{A,{\nu }_i}^2{)}^2[8({\Phi }_1{U}_2+{\Phi }_2{U}_1)$$

$$-2({\Phi }_{-1}{U}_2+{\Phi }_2{U}_{-1})+7({\Phi }_0{U}_1+{\Phi }_1{U}_0)]\}$$

$$+\{[5({\Phi }_0{U}_{-1}+{\Phi }_{-1}{U}_0)]$$

$$+9(C_{V,{\nu }_i}^2-C_{A,{\nu }_i}^2{)}^2[{\Phi }_0(U_1+U_{-1})+({\Phi }_{-1}+{\Phi }_1)U_0]\},$$

$$q_{\mathrm{n}+\mathrm{n}\to \mathrm{n}+\mathrm{n}+{\nu }_i+{\overline{\nu }}_i}\approx 1.5\times 10^{27}{\rho }_{10}^2{T}_{11}^{5.5}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3}{\mathrm{s}}^{-1},$$

Full text

††journal: New Astronomy Reviews

1 Introduction

In astrophysics, accretion is a process that matter falls to a central object, which prompts part of gravitational binding energy of the infalling matter converted into heat and radiation due to viscous dissipation. The conservation of angular momentum, however, forces the accreted matter to form an accretion disk around the central object. Accretion disks widely exist in astrophysical systems, such as cataclysmic variable stars, X-ray binaries, protoplanetary disks, active galactic nuclei (AGNs), and gamma-ray bursts (GRBs).

As a fundamental physical model, black hole (BH) accretion disks have been widely studied (see reviews by Frank et al. 2002; Kato et al. 2008; Abramowicz & Fragile 2013; Blaes 2014; Yuan & Narayan 2014). Three classic accretion disk models, namely the Shakura-Sunyaev disk (SSD, Shakura & Sunyaev 1973), the slim disk (Abramowicz et al. 1988), and the advection-dominated accretion disk (ADAF, Narayan & Yi 1994; Abramowicz et al. 1995), have been successfully applied to different systems. The SSD model is geometrically thin, optically thick, and Keplerianly rotating, where the viscous heating is balanced by the radiative cooling. Such a model is very successful in interpreting the high/soft state of X-ray binaries and can be even used to measure the spin of the BH (e.g., Zhang et al. 1997). It is also widely applied to high luminosity AGNs such as quasars. The slim disk model was introduced mainly for systems with super-Eddington accretion, where the disk is geometrically slim and optically thick. The fundamental difference from the SSD model is that the large amount of photons generated by the viscosity cannot escape from the disk. Most of the photons are carried by the accretion flow and finally fall into the BH. In other words, the main cooling mechanism is advection rather than radiation. The slim disk model is often applied to super-Eddington systems such as ultraluminous X-ray sources and narrow line Seyfert 1 galaxies, which is also applied to study cosmology (Wang et al. 2013). Different from the above two models, an ADAF is optically thin and has extremely high temperature. The main cooling mechanism is also advection rather than radiation. The difference from the slim disk is that, the energy advection in an ADAF is related to the internal energy of the flow instead of the photons. The ADAF model has been successfully applied to the low/hard and quiescent states of X-ray binaries and low luminosity AGNs.

Apart from the above three classic accretion models, there are some significant progresses in this field. The advection-dominated inflow outflow solution (ADIOS, Blandford & Begelman 1999) shows that the outflows may have an essential effect on the structure and radiation of the flow. The convection-dominated accretion flow (CDAF, Narayan et al. 2000) includes energy and angular momentum transfers by radially convective process. A luminous hot accretion flow (LHAF, Yuan 2001) can provide both high luminosity and hard photons.

A well-known unified description of different accretion models is in the $\dot{M}$-$\Sigma$ parameter space (Figure 1), where $\dot{M}$ is the mass accretion rate and $\Sigma$ is the mass surface density. By including the thermally unstable Shapiro-Lightman-Eardley disk (SLE, Shapiro et al. 1976), Abramowicz et al. (1995) first presented a unified picture of the thermal equilibrium solutions containing the above three classic models. Such a picture was improved by Yuan et al. (2003), where the LHAF solution is included.

All the above mentioned accretion flows are photon radiation dominated and neutrino radiation is negligible. A main ingredient of this review, on the other hand, we will focus on the neutrino radiation-dominated disk which may be related to GRBs.

GRBs are extremely energetic transient events in the Universe and isotropically distributed over the sky (see e.g., Mészáros 2002; Zhang & Mészáros 2004; Zhang 2007; Gehrels et al. 2009; Kumar & Zhang 2015; Wang et al. 2015 for reviews), which are sorted into two categories, i.e., short- and long-duration GRBs (SGRBs and LGRBs, see Kouveliotou et al. 1993) or Type I and II GRBs (Zhang 2006; Zhang et al. 2007). Milliseconds variability requires that GRBs are the stellar scale events, and they may be related to the stellar evolution in the external galaxies as well as cosmology. Moreover, their typically isotropic energy is about $10^{50}-10^{52}$ $\rm erg$, which asks for efficiently released energy approach. That is, a compact object should exist in the center of GRBs. SGRBs and LGRBs are generally considered to be mergers of two compact objects, i.e., two neutron stars (NSs) or a BH and an NS (for reviews, see, e.g., Nakar 2007; Berger 2014), and collapses of massive stars (e.g., Woosley & Bloom 2006 for reviews), respectively. The popular models on the central engines of GRBs are either a rotating stellar BH surrounded by a hyperaccretion disk (e.g., Paczyński 1991; Narayan et al. 1992; MacFadyen & Woosley 1999) or a quickly rotating magnetar (or protomagnetar, e.g., Usov 1992; Dai & Lu 1998; Zhang & Mészáros 2001; Dai et al. 2006; Metzger et al. 2011; Lü & Zhang 2014; Lü et al. 2015; Li et al. 2016).

In the inner region of a hyperaccretion disk the density and temperature are so high ($\rho\sim$ $10^{10}-10^{13}$ g cm*-3*, $T\sim 10^{10}-10^{11}$ K) that photons are completely trapped. A large amount of energetic neutrinos are emitted from the surface of the disk, carrying away the viscous dissipation energy of accreted gas via the reactions that neutrinos participate. Annihilation of a fraction of neutrinos and antineutrinos produces a relativistic electron-positron pair dominated outflow, which can power a GRB. The hyperaccretion disk is referred as the neutrino-dominated accretion disk (NDAF). The properties of NDAFs were first investigated in details by Popham et al. (1999). They proposed that NDAFs around stellar-mass rotating BHs are a plausible candidate for the central engine of GRBs.

Figure 1 shows the unified description of the thermal equilibrium solutions of different accretion models, where the NDAF model is also included at extremely high accretion rates. This figure can be regarded as an extension of the unified picture of Abramowicz et al. (1995), where the SLE disk is included whereas the LHAF is not considered. As the figure shows, slim disks and NDAFs have no separation boundary, and NDAFs can be considered as the naturally extended branch of slim disks for very high accretion rates. Certainly, NDAFs are extremely optically thick due to the high accretion rates and the corresponding high surface density. Figure 2 shows the very different components between the classic accretion disks and NDAFs, which result from the significant difference in the mass accretion rates. The inner region of NDAFs is dominated by free baryons, which will be discussed in detail when in the results of the radial and vertical directions of NDAFs in Subsections 2.1 and 4.5.

According to the observations on GRBs and their afterglows, NDAF as a plausible and central engine model of GRBs must satisfy the following requirements: First, the released energy from NDAFs should be adequate for GRBs, and NDAFs should be able to launch an ultra-relativistic outflow with an opening angle of $\sim$ 0.1 and very low baryon contamination. NDAFs should be able to produce GRB emission with diverse temporal behaviors from single pulse to multi-pulses. These have been reviewed in Sections 2-5. Second, NDAFs should be able to reactivate to produce X-ray flares and extended X-ray emission after the prompt emission phase. These are reviewed in Section 6. Third, NDAFs should be produce diverse GRBs with different jet compositions, including both matter dominated fireball (neutrino-anti-neutrino annihilations) and Poynting flux dominated outflow [the Blandford-Znajek (BZ) mechanism (Blandford & Znajek 1977)]. These are discussed in Sections 3 and 6.

In this review, we summarize the theoretical progresses and applications of NDAFs, and their possible detectable effects. In Sections 2 and 4, we describe the radial and vertical dynamics and radiation of NDAFs. In Section 3, magnetized NDAF models are discussed. Some simulations of NDAFs are briefly introduced in Section 5. In Section 6, numerical applications to GRB observations are presented. We review the validation of the existence of NDAFs in Section 7. Summary is given in Section 8.

2 Radial dynamics and radiation of NDAFs

2.1 Radial dynamics of NDAFs

The steady radial structure and neutrino luminosity of BH NDAFs are most widely studied (see e.g, Popham et al. 1999; Narayan et al. 2001; Di Matteo et al. 2002; Kohri & Mineshige 2002; Kohri et al. 2005; Lee et al. 2005; Gu et al. 2006; Chen & Beloborodov 2007; Liu et al. 2007; Janiuk et al. 2007; Rossi et al. 2008; Janiuk & Yuan 2010; Kawanaka & Mineshige 2007; Kawanaka et al. 2013b; Li & Liu 2013; Luo & Yuan 2013; Xue et al. 2013). After the first work of Popham et al. (1999), Di Matteo et al. (2002) recalculated the NDAF radial structure by taking into account neutrino opacity with Newtonian potential. They found that neutrinos are sufficiently trapped in the inner region of the disk, which may suggest that NDAFs may be inadequate to power bright GRBs. Gu et al. (2006) argued that when the general relativistic effects are considered, NDAFs can be powerful enough to be the central engine for GRBs. In another discussion, Kohri & Mineshige (2002) pointed out that electron degeneracy can suppress neutrino emission. Kohri et al. (2005) started to introduce detailed neutrino physics into NDAF model calculations, which brings significant improvement to the NDAF model.

As most detailed study of NDAFs, Xue et al. (2013) investigated the relativistic one-dimensional (1D) steady global solutions of NDAFs by taking into account detailed neutrino physics, balance of chemical potentials, photodisintegration, and nuclear statistical equilibrium. This is mainly introduced below.

2.1.1 Relativistic hydrodynamics

The hydrodynamical model of disks is based on the ADAF model of Abramowicz et al. (1996), the slim disk model of Sa̧dowski (2009) and Abramowicz & Fragile (2013), and the NDAF model of Popham et al. (1999) and Xue et al. (2013), which are all research on the 1D global solutions of accretion disks in Kerr metric. We describe in the units of $G=c=M=1$ ($M$ is the BH mass) just in this subsection.

The basic equations in Kerr geometry include:

(I) The equation of mass conservation (or the continuity equation) is

[TABLE]

where $\dot{M}$ is the rest-mass accretion rate, $\rho$ is the rest-mass density, $H$ is the half thickness of the disk, $v_{r}$ is the radial velocity measured in the corotating frame, $\Delta\equiv r^{2}-2r+a^{2}$ is a function of the Boyer-Lindquist radial coordinate $r$, and $a$ is the total specific angular momentum of the BH.

(II) The equation of radial momentum conservation is

[TABLE]

where $p$ is pressure, $\lambda\equiv(\rho+p+u)/\rho$ is relativistic enthalpy, $u$ is specific internal energy, and the $\mathcal{A}$ term combines the effects of gravity and rotation,

[TABLE]

where $\Omega$ is the angular velocity with respect to the stationary observer, $\Omega^{\pm}_{K}\equiv\pm(r^{3/2}\pm a)^{-1}$ are the angular velocities of the corotating and counter rotating Keplerian orbits, $\tilde{\Omega}\equiv\Omega-2ar/A$ is the angular velocity with respect to the local inertial observer, $\tilde{R}\equiv A/(r^{2}\Delta^{1/2})$ is the radius of gyration, and $A\equiv r^{4}+r^{2}a^{2}+2ra^{2}$.

(III) The equation of angular momentum conservation is

[TABLE]

where $\mathcal{L}$ is the specific angular momentum of the disk, $\mathcal{L}_{\rm{in}}$ is the specific angular momentum at the inner edge of the disk, and the Lorentz factor $\gamma$ is written as

[TABLE]

(IV) The equation of vertical mechanical equilibrium (Abramowicz et al. 1997) is

[TABLE]

where $\epsilon$ is the energy at infinity, which is expressed as

[TABLE]

(V) The equation of energy conservation is

[TABLE]

where $\alpha$ is the viscosity parameter, and $Q^{-}$ is the total cooling rate (per unit area of a half-disk above or below the equator). The term on the left hand side is the advective cooling rate $Q_{\rm adv}$ (per unit area of a whole disk), and the first term on the right represents the viscous heating rate $Q_{\rm vis}$ (per unit area of a whole disk).

2.1.2 Neutrino physics

The cooling mechanism is the main difference between NDAFs and typical accretion disks. Neutrino radiation is dominantly in NDAFs, so microphysics, especially neutrino physics should be considered. In great detail, the microphysics in NDAFs is extended from that in NSs (e.g., Shapiro & Teukolsky 1986) because of close resemblance.

The cooling process of NDAFs is much more complicated than that of classic accretion disks. In most previous works, neutrino physics is simplified through parameterizing the electron fraction $Y_{\rm e}=0.5$ or considering the effect of 4He. Below we present a detailed description of neutrino physics in NDAFs.

We first introduce the total optical depth for neutrinos

[TABLE]

where $\tau_{s,\nu_{i}}$ and $\tau_{a,\nu_{i}}$ are the neutrino optical depth from scattering and absorption, respectively, the subscript “$i$” runs for the three species of neutrinos $\nu_{\rm e}$ , $\nu_{\rm\mu}$ , and $\nu_{\rm\tau}$.

Specifically, the optical depth for neutrinos through scattering off electrons and nucleons $\tau_{s,\nu_{i}}$ is given by

[TABLE]

where $\sigma_{{\rm e},\nu_{i}}$, $\sigma_{j,\nu_{i}}$, $n_{\rm e}$ and $n_{j}$ ($j=1$, 2,…) are the cross sections of electron and nucleons ($n_{1}$ and $n_{2}$ are the number density of free protons and free neutrons), and the number density of electrons and nucleons ($j\geq 3$), respectively (e.g., Kohri et al. 2005; Chen & Beloborodov 2007; Kawanaka & Mineshige 2007; Liu et al. 2007, 2012a; Xue et al. 2013). The major cross sections from scattering off electrons, free protons, free neutrons and other elements particles are given by (Burrows & Thompson 2004; Chen & Beloborodov 2007)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $k_{\rm B}$ is the the Boltzmann constant and $\eta_{\rm e}$ is electron degeneracy, $\sigma_{0}=4G^{2}_{F}(m_{\rm e}c^{2})^{2}/\pi(\hbar c)^{4}\approx 1.71\times 10^{-44}{\rm cm^{2}}$, $G_{F}\approx 1.436\times 10^{-49}{\rm erg~{}cm^{3}}$, $e_{\nu_{i}}$ is the mean energy of neutrinos in units of $(m_{\rm e}c^{2})$, $g_{A}\approx 1.26$, $\sin^{2}\theta_{W}\approx 0.23$, $Z_{j}$ and $N_{j}$ are defined as the number of protons and neutrons of a nucleus $X_{j}$, and $C_{V,\nu_{e}}=1/2+2\sin^{2}\theta_{W}$ , $C_{V,\nu_{\mu}}=C_{V,\nu_{\tau}}=-1/2+2\sin^{2}\theta_{W}$ , $C_{A,\nu_{e}}=C_{A,\overline{\nu}_{\mu}}=C_{A,\overline{\nu}_{\tau}}=1/2$, $C_{A,\overline{\nu}_{e}}=C_{A,\nu_{\mu}}=C_{A,\nu_{\tau}}=-1/2$. The electron degeneracy is an important physical parameter that affects electron fraction, degeneracy pressure, and neutrino cooling (Kohri et al. 2005; Chen & Beloborodov 2007; Liu et al. 2007; Kawanaka & Mineshige 2007).

The number density of electrons and positrons can be calculated by the Fermi-Dirac integration (see, e.g., Kohri et al. 2005; Kawanaka & Mineshige 2007; Liu et al. 2007; Xue et al. 2013),

[TABLE]

where $\mu_{\rm e}=\eta_{\rm e}k_{\rm B}T$ is the chemical potential of electrons.

The absorption depth for neutrinos $\tau_{a,\nu_{i}}$ is defined by

[TABLE]

where $q_{\nu_{i}}$ is the total neutrino cooling rate (per unit volume) and is the sum of four terms,

[TABLE]

Urca processes have been included in the proton-rich nuclear statistical equilibrium (NSE, Seitenzahl et al. 2008; Liu et al. 2013; Xue et al. 2013). The neutrino cooling rate due to the Urca processes $q_{\rm Urca}$ relates only to electron neutrino and antineutrino. There are four major terms by electrons, positrons, free protons, free neutrons and nucleons (Liu et al. 2007; Kawanaka & Mineshige 2007), which is expressed by

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\cos^{2}\theta_{c}\approx 0.947$, $Q=(m_{\rm n}-m_{\rm p})c^{2}$, $Q^{\prime}\approx\mu^{\prime}_{\rm n}-\mu^{\prime}_{\rm p}+\Delta$, $\mu^{\prime}_{\rm n}$ and $\mu^{\prime}_{\rm p}$ are the chemical potential of protons and neutrons in their own nuclei, $\Delta\approx 3\rm MeV$ is the energy of the neutron 1$f_{5/2}$ state above the ground state, and $f_{{\rm e}^{\mp}}=\{\exp[(E_{\rm e}\mp\mu_{\rm e})/k_{\rm B}T]+1\}^{-1}$ is the Fermi-Dirac function. $N_{p}(Z_{j})$ and $N_{h}(N_{j})$ are satisfied with

[TABLE]

[TABLE]

The electron-positron pair annihilation rate into neutrinos $q_{{\rm e^{-}+e^{+}}\rightarrow\nu_{i}+\overline{\nu}_{i}}$ is (e.g., Itoh et al. 1989; Yakovlev et al. 2001; Janiuk et al. 2007)

[TABLE]

where $Q_{c}=({m_{\rm e}c}/{\hbar})^{9}{G^{2}_{F}}/{\hbar}\approx 1.023\times 10^{23}\rm erg~{}cm^{-3}~{}s^{-1}$, and the dimensionless functions $U_{k}$ and $\Phi_{k}$ ($k=-1,0,1,2$) in the above equation can be expressed in terms of the Fermi-Dirac functions (Kawanaka & Mineshige 2007). Once electrons are in degenerate state, this process can be ignored.

The nucleon-nucleon bremsstrahlung rate $q_{{\rm n+n\rightarrow n+n+}\nu_{i}+\overline{\nu}_{i}}$ is the same for the three species of neutrinos (e.g., Itoh et al. 1996; Di Matteo et al. 2002; Liu et al. 2007), which can be simplified as

[TABLE]

where ${\rho}_{10}\equiv{\rho}/10^{10}{\rm g~{}cm^{-3}}$ and $T_{11}\equiv T/10^{11}{\rm K}$.

The plasmon decay rate for the three species of neutrinos $q_{\tilde{\gamma}\rightarrow\nu_{i}+\overline{\nu}_{i}}$ also needs to be considered, where plasmons $\tilde{\gamma}$ are photons interacting with electrons (e.g., Kawanaka & Mineshige 2007; Xue et al. 2013),

[TABLE]

[TABLE]

where $\alpha^{*}\approx 1/137$ is the fine-structure constant, and $\gamma\approx 5.565\times 10^{-2}{[(\pi^{2}+3\eta_{e}^{2})/3]}^{1/2}$. Here $q_{{\rm n+n\rightarrow n+n+}\nu_{i}+\overline{\nu}_{i}}$ and $q_{\tilde{\gamma}\rightarrow\nu_{i}+\overline{\nu}_{i}}$ may become important only for very high electron degeneracy state.

Second, the electron fraction can be written as (Liu et al. 2013; Xue et al. 2013)

[TABLE]

Moreover, the condition of electrical neutrality naturally constrains $Y_{\rm e}$ (Liu et al. 2007, 2013), i.e.,

[TABLE]

where $m_{u}$ is the mean mass of nucleus, and the mass fraction is considered to approximately equal the number density.

For NDAFs, in order to allow for a transition from the optically thin [$\mu_{\rm n}=\mu_{\rm p}+2\mu_{\rm e}$ (Yuan 2005), where $\mu_{\rm n}$, $\mu_{\rm p}$, and $\mu_{\rm e}$ are the chemical potential of free neutrons, free protons, and electrons, respectively] to optically thick ($\mu_{\rm n}=\mu_{\rm p}+\mu_{\rm e}$) regimes, the bridging formula of free protons and neutrons can be established, which is given by (Liu et al. 2007; Xue et al. 2013)

[TABLE]

where $f(\tau_{\nu})=\rm exp(-\tau_{\nu_{\rm e}})$ is a weight factor. In addition, the bridging formula can be generally used in NDAFs if when nucleosynthesis is taken (see Subsection 2.1.3 below), because the outer region of the NDAF is optically thin.

If we assume that the heaviest nuclei is ${}^{4}\rm He$, by combing with the above equations, the bridging formula can be simplified as (Liu et al. 2007)

[TABLE]

where $X_{\rm nuc}$ is the mass fraction of free nucleons (e.g., Di Matteo et al. 2002; Kohri et al. 2005; Chen & Beloborodov 2007; Liu et al. 2007), which can be reduced from the equation of NSE (e.g., Meyer 1994).

2.1.3 Nucleosynthesis

As mentioned above, NSE is established by all nuclear reactions in the chemical equilibrium. Seitenzahl et al. (2008) developed a method of NSE especially for proton-rich material, which is applicable to almost the entire range of the electron fraction. The complicated and detailed balance has been included under the condition of chemical potential equilibrium. Although they focus on the proton-rich state of matter, it can be used in the description of all the state of matter. The conclusion is that, $\rm{}^{56}Ni$ is favored under proton-rich conditions ($Y_{\rm e}\simeq 0.5$), which is very different from the case of Fe-peak nuclei domination with the largest binding energy per nucleon that have a proton to nucleon ratio close to the prescribed electron fraction (e.g., Lattimer & Swesty 1991). This method is suitable for the calculations of NDAFs whatever origins from compact objects mergers or collapsars.

One in particular is that the lower limit of the temperature in the NSE calculation is identified at about $2\times 10^{9}$ K. Once the temperature is lower than this limit, their NSE solution is unreliable. Furthermore, Kawanaka & Mineshige (2007) assumed that the inflowing gas is composed primarily of neutron-rich iron group nuclei as the outer boundary conditions, and the maximum electron fraction is less than 0.42. Anyway, it is predictably that the cooling process of the matter in the outflow of NDAFs can reveal the history of nucleosynthesis in very short timescale.

2.1.4 Thermodynamics

Distinguished from the classic accretion disks, the pressure from degenerate electrons and neutrinos are included in the equation of state (EoS), i.e.,

[TABLE]

The gas pressure from free nucleons $p_{\rm gas}$ can be estimated by

[TABLE]

The disk is definitely optically thick for the photons, so photon radiation pressure $p_{\rm rad}$ can be given by

[TABLE]

The electron pressure $p_{\rm e}$ is the sum of the electron and positron pressure, which is described by the exact Fermi-Dirac distribution (e.g., Chen & Beloborodov 2007; Liu et al. 2007; Xue et al. 2013). In NDAFs, electrons are neither nondegenerate nor strongly degenerate, and they are not ultrarelativistic at all radii. One has

[TABLE]

where

[TABLE]

The neutrino pressure $p_{\rm\nu}$ is

[TABLE]

where ${u}_{\nu}$ is the energy density of the neutrinos, which is adopted by a bridging formula connecting the optically thin region with the optically thick region (e.g., Di Matteo et al. 2002; Kohri et al. 2005; Liu et al. 2007; Xue et al. 2013),

[TABLE]

In fact, this equation is derived from ADAF models via the approximation of the radiation transfer equations (e.g., Hubeny 1990; Narayan & Yi 1995b; Popham & Narayan 1995; Artemova et al. 1996).

Li & Liu (2013) revisited various properties of the hot nuclear matter possible in the inner regions of GRBs and supernovae (SNe). They employed the microscopic Brueckner-Hartree-Fock approach to account for the strong interaction between nucleons, and calculated the nucleonic chemical potentials and the nucleonic EoS. A parameterized chemical potential equilibrium bridging between neutrino optically thin and thick regions was introduced. This method can update the above descriptions of EoS and neutrino physics in the inner regions of NDAFs.

A fraction of the viscous heating energy following inflows has been advected into the BH, and the rest part equals the cooling energy. The cooling rate $Q^{-}$ mentioned in Equation (8) is composed of the cooling rates of photodisintegration, neutrino emission, and photon radiation

[TABLE]

However, cooling due to photon radiation can be ignored because NDAFs are optically thick.

The cooling rate by photodisintegration can be given by the NSE. For the case that the heaviest nucleus is $\alpha$-particles, $Q_{\rm ph}$ can be written as

[TABLE]

where $v_{r,10}=v_{r}/(10^{10}~{}\rm cm~{}s^{-1})$. The cooling rate by disintegration of other heavy nuclei can be ignored because of the lower number density of these nuclei and the absolutely dominant advective cooling rate in the outer region.

The cooling rate due to neutrino loss $Q_{\nu}$ is expressed in accordance to the above bridging formula of ${u}_{\nu}$ (e.g., Di Matteo et al. 2002; Kohri et al. 2005; Liu et al. 2007; Xue et al. 2013),

[TABLE]

2.1.5 Characteristic radii in NDAFs

One can define three characteristic radii in NDAFs (Chen & Beloborodov 2007) as shown in Figure 3.

(I) Neutrino trapping radius $r_{\rm trap}$

Photon trapping may occur in the supercritical accretion disk if the time of photon diffusion from the equatorial plane of the disk to the surface is longer than the accretion timescale (e.g., Katz 1977; Begelman 1978; Ohsuga et al. 2002, 2005). Photons generated near the equatorial plane diffuse toward the disk surface at a speed of $\sim c/3\tau$ (Mihalas & Mihalas 1984), where $\tau$ is the total optical depth. The timescale of photon diffusion is about $t_{\rm diff}\sim H/(c/3\tau)$. Similarly, neutrinos may also be trapped in NDAFs (Chen & Beloborodov 2007; Liu & Xue 2012; Xue et al. 2013). The neutrino velocity $v_{n}$ replaces $c$ in the above equation of $t_{\rm diff}$, which can be estimated by $\sim(3.7k_{\rm B}Tc^{2}/0.07\rm eV)^{1/2}$, where $\sim 3.7k_{\rm B}T$ and $0.07\rm eV$ roughly equals to the neutrino energy and the lower limit of neutrino rest-mass energy, respectively. The accretion timescale can be estimated by $t_{\rm acc}\sim-r/{v_{r}}$. One can then use $\tilde{t}=t_{\rm diff}/t_{\rm acc}=1$ to define the trapping radius,

[TABLE]

Once neutrinos are trapped in the inner region, the neutrinos launched from the surface would decrease. As a result, neutrino trapping would greatly affect the neutrino luminosity and annihilation luminosity for high accretion rates ($\dot{M}\gtrsim 2~{}M_{\odot}~{}\rm s^{-1}$) and large spin ($a_{*}\sim 0.95$, where $a_{*}$ is the dimensionless spin parameter of the BH) (Liu et al. 2012a; Xue et al. 2013). Moreover, one can set a critical accretion rate $\dot{M}_{\rm trap}$ when $r_{\rm trap}$ appears in NDAFs.

(II) Ignition radius $r_{\rm ign}$

The region where $r<r_{\rm ign}$ meets that the neutrino emission switches on and neutrino cooling dominates. Accurately, $r_{\rm ign}$ is defined as the radius satisfied with $Q^{-}/Q_{\rm vis}=0.5$ (Chen & Beloborodov 2007; Zalamea & Beloborodov 2011; Liu & Xue 2012). Once $r_{\rm ign}$ emerges, NDAFs are ignited, the corresponding $\dot{M}$ is defined as $\dot{M}_{\rm ign}$. For a rapidly rotating BH with $3~{}M_{\odot}$ surrounded by a NDAF with low viscosity, $\dot{M}_{\rm ign}$ is about 0.001 $M_{\odot}~{}\rm s^{-1}$. $\dot{M}_{\rm trap}$ and $\dot{M}_{\rm ign}$ mainly depend on the spin of BHs and the viscous parameter of disks in BH-NDAF systems.

(III) Neutrino opaque radii $r_{\nu}$ and $r_{\bar{\nu}}$

Radii $r_{\nu}$ and $r_{\bar{\nu}}$ locate where the disks become opaque for neutrinos and antineutrino (Chen & Beloborodov 2007), and the radiation becomes thermal. $r_{\nu}$ is usually larger than $r_{\bar{\nu}}$ due to the neutrino optical depth as shown in Subsection 2.1.2, and they are definitely less than $r_{\rm ign}$.

2.1.6 Results

Figure 4 shows the basic structure of NDAFs with $M=3M_{\odot}$ and $\alpha=0.1$. The six panels correspond with the variations of density $\rho$, temperature $T$, radial velocity $v_{r}$, electron degeneracy $\eta_{\rm e}$, optical depth of electron neutrino $\tau_{\nu_{\rm e}}$, and electron fraction $Y_{\rm e}$ with $r/r_{\rm g}$, where $r_{\rm g}=2GM/c^{2}$. From the outer to inner region of the disk, $\rho$, $T$, and $\tau_{\nu_{\rm e}}$ increase by about 6, 1, and 5 orders of magnitude. In the innermost region, they reach about $10^{13}~{}\rm g~{}cm^{-3}$, $10^{12}~{}\rm K$, and $10^{3}$, respectively. The value of $Y_{\rm e}$ at the outer boundary of the disk in sixteen solutions all tend to about 0.46. There is a little different value $\sim 0.42$ in Kawanaka & Mineshige (2007), but this difference leads to hugely different nucleosynthesis products.

Figure 5 shows the radial distributions of the mass fractions of seven major nucleons cover almost $99\%$ mass of the flow, including $\rm{}^{1}n$, $\rm{}^{1}H$, $\rm{}^{4}He$, $\rm{}^{52}Cr$, $\rm{}^{54}Cr$, $\rm{}^{56}Fe$ and $\rm{}^{58}Fe$. The mass fraction of $\rm{}^{56}Fe$, $\rm{}^{4}He$, and free neutrons and protons dominate in the outer, middle, and inner regions of the disk for all accretion rates. It is an implication for the origin of heavy nuclei in GRBs accounting for the suspected detection of Fe K$\alpha$ X-ray lines and other emission lines (e.g., Lazzati et al. 1999; Kallman et al. 2003; Gou et al. 2005; Butler 2007), which can play an important role in understanding the nature of GRBs, especially its central engine. For a comparison, Janiuk (2014) obtained that 56Ni dominates out of about 500 $r_{\rm g}$ of the disk with $a_{*}=0.9$ by using simplified neutrino physics, which is similar to the vertical results in Liu et al. (2013). We will detailed discuss nucleosynthesis of NDAFs in Subsection 6.4.

2.2 Neutrino luminosity and annihilation luminosity

Once the neutrino cooling rate $Q_{\nu}$ is obtained, the neutrino radiation luminosity before annihilation, $L_{\rm\nu}$, can be calculated by

[TABLE]

where $r_{\rm in}$ and $r_{\rm out}$ are the inner and outer edge of the disk. Here we sum two jets in the opposite direction. For comparison with the observations, the coefficient “4$\pi$” should be replaced by “2$\pi$” in the above equation. The same applies to the equation of the neutrino annihilation luminosity below.

In the neutrino annihilation calculations, a Newtonian approach is introduced in Ruffert et al. (1997), Popham et al. (1999), and Rosswog et al. (2003). Strictly, the general relativistic effects on the neutrino trajectory near the BH need to be considered (e.g., Birkl et al. 2007; Zalamea & Beloborodov 2011). The disk is modeled as a grid of cells in the equatorial plane. As shown in Figure 6, a cell $k$ has its mean neutrino energy $\varepsilon_{\nu_{i}}^{k}$, neutrino radiation luminosity $l_{\nu_{i}}^{k}$, and distance to a spatial point above (or below) the disk $d_{k}$. At any spatial point, the angle at which a neutrino from cell $k$ encounter an antineutrino from another cell $k^{\prime}$ is denoted as $\theta_{kk^{\prime}}$. Then the neutrino annihilation luminosity at that point is given by the summation over all pairs of cells,

[TABLE]

where $A_{1,i}=(1/12\pi^{2})[\sigma_{0}/c{(m_{\rm e}c^{2})}^{2}][{(C_{V,\nu_{i}}-C_{A,\nu_{i}})}^{2}+{(C_{V,\nu_{i}}+C_{A,\nu_{i}})}^{2}]$, and $A_{2,i}=(1/6\pi^{2})(\sigma_{0}/c)$ $(2C_{V,\nu_{i}}^{2}-C_{A,\nu_{i}}^{2})$, with $C_{V,\nu_{i}}$ and $C_{A,\nu_{i}}$ given below Equation (14). The total neutrino annihilation luminosity is the integration over the entire space outside the disk,

[TABLE]

The neutrino radiation luminosity and annihilation luminosity depend on with the BH spin and accretion rate. Two analytic formulae can be derived by fitting results in Xue et al. (2013),

[TABLE]

[TABLE]

which are applicable for the accretion rate in the range of $0.01\lesssim\dot{m}\lesssim 10$, where $\dot{m}=\dot{M}/(M_{\odot}~{}\rm s^{-1})$.

Including the effect of BH mass, the annihilation luminosity is approximated as (Liu et al. 2016b)

[TABLE]

where $m=M/M_{\odot}$. It should be emphasized that this formula is applicable for $0.01\lesssim\dot{m}\lesssim 0.5$.

Fryer et al. (1999) also displayed the approximate fit to the annihilation luminosity results of Popham et al. (1999), i.e.,

[TABLE]

which is satisfied with $0.01\lesssim\dot{m}\lesssim 0.1$. The comparison of the above three annihilation formulae is shown in Figure 7.

Another well-known analytic formula of neutrino annihilation luminosity is given by Zalamea & Beloborodov (2011), which is expressed as

[TABLE]

where $x_{\rm ms}=r_{\rm ms}/r_{\rm g}$ is the dimensionless marginally stable orbit radius of the disk, $r_{\rm ms}$ is radius of the last marginally stable orbit, and $\dot{m}_{\rm ign}=\dot{M}_{\rm ign}/(M_{\odot}~{}\rm s^{-1})$, $\dot{m}_{\rm trap}=\dot{M}_{\rm trap}/(M_{\odot}~{}\rm s^{-1})$. Here $x_{\rm ms}=3+Z_{2}-\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})}$, where $Z_{1}=1+(1-a_{*}^{2})^{1/3}[(1+a_{*})^{1/3}+(1-a_{*})^{1/3}]$ and $Z_{2}=\sqrt{3a_{*}^{2}+Z_{1}^{2}}$ for $0<a_{*}<1$ (e.g., Bardeen et al. 1972; Kato et al. 2008; Zalamea & Beloborodov 2011; Liu et al. 2015c).

The reasons of the different forms include: (1) different levels of neutrino physics is considered, which cause the slight effects on the annihilation luminosity; (2) a different applying range of $\dot{M}$ is adopted; (3) we have fitted the numerical results directly for BH mass instead of introducing some analytical results as did in Zalamea & Beloborodov (2011); (4) the difference between Newtonian approach or relativistic method also leads to the slight effects (Yang et al. 2017).

In order to reveal the spatial distribution of neutrino annihilation luminosity, we plot in Figure 8 contours of ($2\pi rl_{\nu\overline{\nu}}$) in units of (${\rm erg~{}cm^{-2}~{}s^{-1}}$) with cylindrical coordinates $r$ and $z$. It is easy to find that most of the annihilation events occur in the region of $r\lesssim 20~{}r_{\rm g}$, and it varies more rapidly along the $z$ coordinate than along the $r$ coordinate (also see, Popham et al. 1999). The results indicate that the annihilation luminosity is anisotropic and most of the annihilation energy escapes outward along the angular momentum axis of the disk.

3 Magnetized NDAFs

3.1 Magnetized BH-NDAFs

Blandford & Znajek (1977) suggested that the rotational energy of a BH can be efficiently extracted to power a Poynting jet via a large-scale poloidal magnetic field threading the horizon of the BH. Besides neutrino annihilation above the NDAF, magnetohydrodynamics (MHD) processes, such as the BZ mechanism (Blandford & Znajek 1977) can also power a relativistic jet launched from the central engine of GRBs, even for accretion rates much lower than $\dot{M}_{\rm ign}$ (e.g., Lee et al. 2000a, b; Globus & Levinson 2014; Pan & Yu 2015).

The BZ luminosity can be estimated as (e.g., Lei et al. 2005; Krolik & Piran 2011; Kawanaka et al. 2013b; Liu et al. 2015b)

[TABLE]

where $f(a_{*})$ is a factor depending on the specific configuration of the magnetic field including the information of the field configuration (e.g., Tchekhovskoy et al. 2008; Kawanaka et al. 2013b), and $B_{\rm in}$ is the poloidal magnetic field strength near the horizon. The analytical expression of $f(a_{*})$ was attempted (e.g., Blandford & Znajek 1977; Tchekhovskoy et al. 2008), which is an increasing function of $a_{*}$, and in the range from a small number to $\sim 1$ (Hawley & Krolik 2006). Many 2- or 3-dimensional (2D or 3D) MHD simulations have investigated on how a large-scale vertical magnetic field evolves with an accretion disk and which configurations may power a jet (e.g., McKinney & Gammie 2004; Beckwith et al. 2008, 2009; McKinney & Blandford 2009), yet the analytic form of $f(a_{*})$ is still unclear. As simplification, $f(a_{*})=1$ is adopted (Kawanaka et al. 2013b; Liu et al. 2015b), which is suitable for the fast-spinning BH in the center of GRBs.

The magnetic field energy can be estimated by the disk pressure near the horizon $p_{\rm in}$, which is presented as

[TABLE]

where $\beta_{h}\sim 1$ is the ratio of the midplane pressure near the horizon of the BH to the magnetic pressure in the stretched horizon.

Cao et al. (2014) calculated the global solutions of the NDAFs considering the radial advection and diffusion (e.g., Balus & Hawley 1998; Rossi et al. 2008) of the large-scale magnetic field. Although the configuration of the magnetic field and the structure of NDAFs affect each other, they found that only the structure of the inner NDAF changes significantly once magnetic fields are considered. For the magnetic field with $10^{14}~{}\rm G$, the BZ jet luminosity can reach $\sim 10^{53}-10^{54}~{}\rm erg~{}s^{-1}$ for an extreme Kerr BH. Furthermore, the thermal instability111The instabilities in NDAFs will be discussed in Subsection 6.2. in magnetized NDAFs has been studied to explain the GRB variability (Janiuk & Yuan 2010; Xie et al. 2017).

Lei et al. (2013) studied the baryon loading problem of a GRB jet launched by a NDAF under either the neutrino annihilation process and BZ mechanism. They argued that in no matter mechanism, more luminous jets tends to be more baryon poor. Comparing with jet produced by annihilation, the magnetically dominated jet is much cleaner. Thus they suggested that at least a good fraction of GRBs should have a magnetically dominated central engine. Furthermore, the relation between the Lorentz factor $\Gamma$ and isotropic luminosity $L_{\rm iso}$ of the jet can be interpreted for both two models.

Yuan & Zhang (2012) investigated the closed magnetic field lines that continuously emerge out of the hyperaccretion flow. Since the differential rotation of the accretion flow is shear and turbulent, the line may form flux ropes. Once the system loses its equilibrium, the flux rope is thrust outward and then an episodic jet launches. This mechanism can provide enormous amount of energy to trigger GRBs.

Alternatively, the magnetic coupling (MC) from the BH horizon to the inner region of the disk (Li 2000), which can effectively transfer the angular momentum and rotational energy of the BH to heat the inner region of the disk, then radiate a larger number of neutrinos from the disk than those from non-magnetized NDAFs to produce the primordial fireball (e.g., Lei et al. 2005, 2009; Luo et al. 2013). The MC effect affects the EoS and angular momentum of NDAFs, and effectively heats the inner region of the disk. In addition, the inner region of MC-NDAF becomes thermally and viscously unstable, which may be invoked to interpret the GRB variability. Xie et al. (2016) studied magnetized NDAFs with non-zero boundary stresses. They also found that the boundary torque has strong effects on the properties of inner disk. The neutrino annihilation luminosity of such NDAFs may be powerful enough to account for most of GRBs including the ultra-LGRBs (ULGRBs), and viscously unstable emerges in the inner region of the disk, which may be interpreted the origin of GRBs variabilities.

Wu et al. (2013) presented that a BZ jet launched from the BH hyperaccretion system in the fall-back framework (Kumar et al. 2008a, b) can explain the giant X-ray bump in GRB 121027A. Gao et al. (2016) investigated that the BZ jet and Blandford-Payne (BP, Blandford & Payne 1982) outflow can interpret the ULGRB GRB 111209A and its associated SN 2011kl.

3.2 NS-NDAF system

After merger of a compact object binary or collapse of a massive star, a stellar-mass BH around a hyperaccretion disk may form to produce GRBs. Another possible product is an accreting NS or magnetar instead of a BH (e.g., Duncan & Thompson 1992; Usov 1992; Dai & Lu 1998; Zhang & Mészáros 2001; Gao & Fan 2006). Chevalier (1996) first analyzed a system consisting of a neutrino-cooled accretion disk and a NS. In order to maintain the balance of neutrino cooling and viscous heating, the viscous parameter should be small ($\lesssim 2\times 10^{-6}$). Alternatively, a possible option to relax neutrino cooling rate is advection.

Zhang & Dai (2008, 2009) studied the radial structure of a NDAF around a NS based on a two-region scenario, and calculated the corresponding neutrino luminosity and annihilation luminosity. As a result, they found in the weak field case ($\lesssim 10^{15}~{}\rm G$), the structure of the outer disk is very similar to that of a BH NDAF. A NS has a solid surface, which is different from a BH. Most of the advection energy has to be released in the inner region near the NS surface, so that an outflow may be produced. Moreover, the NS NDAF has a brighter neutrino luminosity compared with the BH NDAF caused by the additional neutrino emission from the NS surface boundary layer, and the neutrino annihilation efficiency of a NS-NDAF system should be higher than that of BH-NDAF system. The luminosity $L_{s}$ from the boundary layer at the NS surface is (Frank et al. 2002; Zhang & Dai 2009)

[TABLE]

where $r_{\rm NS}$ and $\Omega_{\rm NS}$ are the radius and angular velocity of the NS, respectively, $\Omega_{\rm in}$ is the angular velocity of the disk inner boundary.

For the strong field case ($\gtrsim 10^{15}~{}\rm G$) corresponding to a magnetar, the NDAF properties may be significantly changed. Zhang & Dai (2010) found that the quantum effects (Landau levels) and the MC processes play two competitive roles in changing the disk properties, with the latter being the main factor to increase pressure, density and neutrino luminosity with increasing magnetic field strength. The strange star-NDAF model has also been studied (Hao & Dai 2013).

4 Vertical structure and luminosity of NDAFs

The simple well-known relationship “$H=c_{\rm s}/\Omega_{\rm K}$” (or “$H\Omega_{\rm K}/c_{\rm s}={\rm constant}$”) was widely adopted in the description of the vertical structure of accretion disks. Such a relationship can be derived from the simply vertical hydrostatic equilibrium

[TABLE]

with two additional assumptions in cylindrical coordinates ($r$, $\phi$, $z$) adopted, i.e., the approximated Hōshi form of gravitational potential (e.g., Hōshi 1977)

[TABLE]

and a one-zone approximation or a polytropic relation in the vertical direction (e.g., Hōshi 1977)

[TABLE]

Obviously, the above assumptions work well for geometrically thin disks, but may be inaccurate when the mass accretion rate $\dot{M}$ approaches the Eddington rate $\dot{M}_{\rm Edd}$, for which the disk is likely not thin as shown in Figure 1. Consequently, the relationship “$H=c_{\rm s}/\Omega_{\rm K}$” may be invalid for $\dot{M}\gtrsim\dot{M}_{\rm Edd}$ (Gu & Lu 2007; Gu et al. 2009; Liu et al. 2010a; Gu 2012; Gu et al. 2016), as well as $\dot{M}\gg\dot{M}_{\rm Edd}$ for NDAFs. The vertical structure of NDAFs should be revisited by using a more accurate description of the gravitational potential and vertically radiative transfer. The results would also affect the value of the neutrino luminosity.

4.1 Gravitational potential

As discussion in Subsection 2.1.1, the hydrodynamics of NDAFs are expected to be similar to those of slim disks. If the angular velocity $\Omega=\Omega_{\rm K}$, the continuity and angular momentum equations of slim disks and NDAFs in cylindrical coordinates can be simplified as (Liu et al. 2008)

[TABLE]

[TABLE]

where $\Omega_{\rm K}=(GM/r)^{1/2}/(r-r_{\rm g})$ is the Keplerian angular velocity, $j=1.8cr_{\rm g}$ is an integration constant representing the specific angular momentum accreted by the BH, and $\Sigma$ and $\Pi$ are the surface density and vertically integrated pressure, respectively, which can be defined as

[TABLE]

[TABLE]

The vertical hydrostatic equilibrium equation is provided by Equation (55), and here we adopt Paczyński-Wiita potential (Paczyński & Wiita 1980)

[TABLE]

to replace the Hōshi form. This modification reverts the real effects of gravity in the vertical direction of the disk, which can bring significant change to the structure of the geometrically thick disks (Liu et al. 2008).

The equation of energy conservation reads

[TABLE]

Here the viscous heating rate is given by

[TABLE]

where $f=1-j/\Omega_{\rm K}r^{2}$, and $g=-{\rm dln}\Omega_{\rm K}/{\rm dln}r$. The advective cooling rate is

[TABLE]

with $\xi=3/2$ being a dimensionless quantity of the order of unity (e.g., Kato et al. 2008; Liu et al. 2008). The sound speed is further defined as $c_{\rm s}=(\Pi/\Sigma)^{1/2}$. The half thickness of the disk can be estimated via $H=\Sigma/2\rho_{0}$, where $\rho_{0}$ is the mass density on the equatorial plane. Note that we can roughly use $Q_{\nu}$ instead of $Q^{-}$ here (Liu et al. 2008).

Figure 9 shows thermal equilibria of SSDs, slim disks, and NDAFs at each radius $r$ with the corresponding accretion rate $\dot{M}$, which is plotted in units of the Eddington accretion rate ${\dot{M}}_{\rm Edd}=64\pi GM/c\kappa_{\rm es}$, where $\kappa_{\rm es}=0.34~{}{\rm cm}^{2}~{}{\rm g}^{-1}$ is the electron scattering opacity, and $M_{\odot}$ s*-1* , respectively. The $\dot{M}$-$r$ plane can be divided into five regions: (1) The region below line $a$ corresponds to the stable, photon radiation-cooled, and gas pressure-supported SSDs; (2) the region between lines $a$ and $b$ corresponds to the unstable, photon radiation pressure-supported SSDs; (3) the region between lines $b$ and $c$ corresponds to the stable, advective cooling-dominated, and photon radiation pressure-supported slim disks (Gu & Lu 2007); (4) line $c$ extends upward till line $d$, and the region between lines $d$ and $e$ corresponds to NDAFs; (5) the region on the right of lines $c$, $d$, and $e$ represents the ‘no solution’ region.

Line $c$ represents a maximal possible accretion rate of slim disks for each radius, and lines $d$ and $e$ represent the lower and upper limits of $\dot{M}$ needed for NDAFs. They are defined by that no thermal equilibrium solutions exist because of viscous heating being always larger than total cooling. The filled circle at the connection between lines $d$ and $e$ at r\approx 185$$r_{\rm g} defines the maximal possible outer boundary of an NDAF. The physical reason is that the correct BH’s gravitational force in the vertical direction can only gather a limited amount of accreted gas. Once the pressure force is larger than the gravitational force, the balance is broken and outflows are produced. Outflows fall back to the BH, restart accretion process and may launch X-ray flares. We also notice that there is no boundary separating slim disks and NDAFs. It is easy to understand because both are very optically thick for photons. Along with the increase of the accretion rate, the neutrino emission processes operate and gradually become important. Thus the accretion flow changes from the slim disk form to the NDAF form. This picture is consistent with the results in the $\dot{M}$-$\Sigma$ plane in Figure 1. It should be noticed that the constraint on the accretion rate of NDAFs in this framework may be too tight because $H\sim r$ in cylindrical coordinates is considered.

4.2 Self-gravity

If the mass density of a disk becomes comparable to $M/r^{3}$, self-gravity becomes important and its resulting local instabilities may develop (see, e.g., Paczyński 1978a, b; Abramowicz et al. 1984; Goodman & Narayan 1988). The effects of self-gravity are generally important in some astrophysical processes, such as AGNs (e.g., King et al. 2008; Hopkins & Quataert 2010) and protostars and protostellar disks (e.g., Goodman 2003; McKee & Ostriker 2007; Rice et al. 2010). Due to the high density of NDAFs, Liu et al. (2014b) argued that the self-gravity effect may be important to the structure of NDAF, and the neutrino luminosity.

First, we define the surface density of the disk in the range from the equatorial plane to a certain height $z$ ($z\leq H$),

[TABLE]

We assume that it varies slowly with radius. A new term $4\pi G\Sigma_{z}$ should be added in the vertical equilibrium equation to represents the effect of self-gravity. Thus, Equation (55) can be rewritten as

[TABLE]

For comparison, we define cases I and II, which correspond to the vertical equilibrium forms with and without self-gravity, respectively.

The Toomre parameter is introduced to measure the local gravitational stability of the accretion disks, which is expressed as

[TABLE]

where $Q_{\rm T}<1$ implies instability. If the effects of the self-gravity are considered in the vertical structure of NDAFs, the gravitational instability should be also reviewed in the framework.

As shown in the left panel of Figure 10, the effects of the self-gravity are mainly relevant in the outer region of the disk, especially for higher accretion rates. The right panel displays the equal $Q$ contours against the disk structure. One can see that $Q=1$ appears at the disk surface at $r/r_{\rm g}\sim 250$ for an accretion rate $1M_{\odot}~{}\rm s^{-1}$. Since the self-gravity has limited effects on the structure of NDAF, there is also no significant influence on the neutrino luminosity (Liu et al. 2014b). The criterion, $Q<1$, indicates that the disk is gravitationally unstable, which may cause two classes of possible behaviors (e.g., Perna et al. 2006): (1) If the local cooling of the disk is rapid, the disk may fragment into two or more parts (e.g., Nelson 2000). Since the fall back timescale is long enough, the accretion processes would restart, which may be related to the origin of late-time X-ray flares in GRBs. (2) The disk may evolve to a quasi-steady spiral structure transferring angular momentum outward and mass inward. This mode may drive long-duration, violent explosions if the disk mass is large enough (e.g., Lodato & Rice 2005), which may give rise to SGRBs with extended emission (e.g., Liu et al. 2012b; Cao et al. 2014) or X-ray flares in LGRBs.

4.3 Vertical hydrostatic equilibrium

The above work is based on the simple vertical hydrostatic equilibrium equation as shown in Equation (55), instead of the general form (Abramowicz et al. 1997; Gu et al. 2009; Liu et al. 2010a),

[TABLE]

where $v_{z}$ is the vertical velocity. Since $v_{z}$ is not negligible for geometrically thick disks, the solutions in Gu & Lu (2007) and Liu et al. (2008) are still not self-consistent. To resolve this issue, we revisit the vertical structure of NDAFs in spherical coordinates ($r$, $\theta$, $\phi$) substituting cylindrical coordinates as shown in Figure 11.

We adopt the radial self-similar assumptions (Narayan & Yi 1995a) to simplify the basic equations of continuity and momentum (see, e.g., Xue & Wang 2005; Gu et al. 2009), then obtain

[TABLE]

[TABLE]

[TABLE]

where $v_{r}$ and $v_{\phi}$ are the radial and azimuthal components of the velocity ($v_{\theta}=0$), and the sound speed $c_{s}$ is defined as $c_{s}^{2}=p/\rho$, the Keplerian angular velocity is $\Omega_{\rm K}=(GM/r^{3})^{1/2}$ here. The continuity equation can be adapted:

[TABLE]

where $\theta_{0}$ is the polar angle of the surface.

The EoS follows Equation (33), and the energy equation is written as

[TABLE]

where the cooling of photodisintegration of $\alpha$-particles is ignored (Liu et al. 2010a, 2012a). The viscous heating rate per unit volume $q_{\rm vis}=\nu\rho r^{2}[\partial(v_{\phi}/r)/\partial r]^{2}$ ($\nu$ is kinematic coefficient of viscosity) and the advective cooling rate per unit volume $q_{\rm adv}=\rho v_{r}(\partial e/\partial r-(p/\rho^{2})\partial\rho/\partial r)$ ($e$ is the internal energy per unit volume) are expressed, after self-similar simplification, we obtain

[TABLE]

[TABLE]

where the entropy of degenerate particles is neglected. Thus the vertical integration of $Q_{\rm vis}$ and $Q_{\rm adv}$ are the following (Liu et al. 2010a, 2012a):

[TABLE]

[TABLE]

The cooling rate due to neutrino radiation $Q_{\nu}$ can be defined as (Lee et al. 2005; Liu et al. 2012a)

[TABLE]

where $k$ represents different types of neutrinos and antineutrinos, and $q_{\nu_{k}}$ is the sum of the cooling rates per unit volume due to the Urca processes, electron-positron pair annihilation, nucleon-nucleon bremsstrahlung, and plasmon decay, as introduced in Subsection 2.1.2. The neutrino optical depth is calculated by using integral to $\theta$ instead of the forms by timing the half thickness of the disk $H$ in Equation (10).

Once gravity cannot hold the radiation pressure gradient, outflows may occur, so a boundary condition is essential. Analogous to the principle of Eddington luminosity, the mechanical equilibrium can be written as (Liu et al. 2012a)

[TABLE]

combined with the photon radiation pressure, the surface temperature can be derived as (Liu et al. 2012a)

[TABLE]

where $\sigma_{\rm T}$ is the Thompson scattering cross section.

Figure 12 shows the variations of density $\rho$, temperature $T$ and electron fraction $Y_{\rm e}$ with $\pi/2-\theta$. The lower limit of the temperature is $\sim 2\times 10^{9}~{}\rm K$ in NSE equations, which is marked in Figure 12(b). The values $\dot{M}=0.05~{}M_{\odot}~{}\rm s^{-1}$ and $\dot{M}=1~{}M_{\odot}~{}\rm s^{-1}$ correspond to $Y_{\rm e}$ around 0.49 and 0.47 near the disk surface, respectively. The half-opening angles of the disks have the positive correlation with accretion rate and radius.

4.4 Outflow

The outflow from the critical accretion disk must be strong, which has been repeatedly demonstrated by theories (e.g., Liu et al. 2008; Gu 2015), simulations (e.g., Jiang et al. 2014; Sa̧dowski & Narayan 2015), and observations (e.g., Wang et al. 2013; Cheung et al. 2016; Parker et al. 2017). In order to discuss the possible outflow from NDAFs, we introduce the Bernoulli parameter of the accreted matter, which is expressed as (e.g., Narayan & Yi 1995a)

[TABLE]

This equation reflects the balance of energy, including the kinetic energy, potential energy and enthalpy of accreted matter. The region where $B>0$ is satisfied implies a possible outflow. Conversely, $B<0$ everywhere is a sufficient condition for the absence of an outflow (Abramowicz et al. 2000).

Figure 13 shows the variation of the half-opening angles of the disk surface $(\pi/2-\theta_{0})$, and the angle at which $B=0$ is satisfied, with radii. Figure 13(a) shows that for $\alpha$ = 0.1 and $\dot{M}$ = 1 $M_{\odot}~{}\rm s^{-1}$ the region with $B>0$ appears at $r\sim 10~{}r_{\rm g}$, and expands continuously with increasing radius. There is no region of neutrino trapping. In Figure 13(b), for $\alpha$ = 0.01 and $\dot{M}$ = 10 $M_{\odot}~{}\rm s^{-1}$ the region with $B>0$ appears at $\sim 9~{}r_{\rm g}$, and increases to near the surface of the disk at $\sim 500~{}r_{\rm g}$.

In the research on the vertical structure of NDAFs, we can further discuss the effects of neutrino trapping. Different from the trapping radius in radial structure of NDAFs, by equating the neutrino diffusing time and the accretion time, we can obtain a region of trapped neutrinos. Figure 13(b) shows the region of neutrino trapping occurs from inner area to $\sim 46~{}r_{\rm g}$, and its opening angle is from $\sim 0.4$ to [math]. We consider that the radius and the open angle of $B=0$ and $\tilde{t}=1$ are closely determined by the accretion rate and the viscous parameter.

4.5 Vertical composition

In the NSE (Seitenzahl et al. 2008; Liu et al. 2013; Xue et al. 2013), the independent variables are the density $\rho$, temperature $T$ and electron fraction $Y_{\rm e}$, which are essential in the vertical NDAF model according to the description above. We thus obtain the vertical distribution of the mass fraction (also approximately equals the number density) of the free neutron and proton, and main elements (include $\rm{}^{4}He$, $\rm{}^{52}Cr$, $\rm{}^{54}Cr$, $\rm{}^{54}Fe$, $\rm{}^{56}Fe$, $\rm{}^{56}Ni$, and $\rm{}^{58}Ni$), as shown in Figure 14. $\rm{}^{56}Ni$ dominates at the disk surface for $\dot{M}=0.05~{}M_{\odot}~{}\rm s^{-1}$, and $\rm{}^{56}Fe$ dominates for $\dot{M}=1~{}M_{\odot}~{}\rm s^{-1}$, corresponding to $Y_{\rm e}$ around 0.49 and 0.47 as shown in Figure 12(c), respectively. The solutions show that the proportion of the nuclear matter increases with radius for the same accretion rate. The mass fraction of $\rm{}^{56}Ni$ or $\rm{}^{56}Fe$ near the surface increases with radius. In the middle region, $\rm{}^{4}He$ is dominant for all the accretion rates. The free neutrons and protons are dominant near the equatorial plane of the disk in the hot and dense state. Most of the free protons turn into free neutrons due to the Urca process (see, e.g., Liu et al. 2007; Xue et al. 2013), which causes the dominance of free neutrons and the decrease of electron fraction. In simple terms, the change of the electron fraction is inversely associated with the accretion rate and radius when the free baryons dominate.

4.6 Convection

In the vertical direction of the disk, the convective motion (or vertical advection) plays an important role, which can carry photons (or neutrinos) from the equatorial plane to the surface of the disk by magnetic buoyancy. For slim disk, the key issue is whether the vertical radiation transfer due to magnetic buoyancy is faster than the radial advection process, thus photons can escape from the surface before being advected into the BH. The results of Jiang et al. (2014) indicated that the radiation can be widely increased because of the vertical advection of radiation caused by magnetic buoyancy. We would like to examine whether the same mechanism exists in NDAFs to increase the efficiency of neutrino emission.

First of all, we estimate two typical timescales, the vertically convective timescale, $t_{\rm c}$, and the radially advective timescale, $t_{\rm adv}$. The vertical convective speed along the vertical direction can be approximatively expressed as (Kawanaka & Kohri 2012)

[TABLE]

This is the vertically “infalling” speed onto the equatorial plane due to the vertical component of the BH’s gravity. Thus the convective timescale is

[TABLE]

where $r_{0}$ is the radius at the equatorial plane.

Meanwhile, the advection timescale (or accretion timescale $t_{\rm acc}$) can be expressed as

[TABLE]

The vertical convection should exist if $t_{\rm c}<t_{\rm adv}$, i.e., the polar angle of the convective region should satisfy $\theta_{0}<\theta<\theta_{\rm c}$. Otherwise, convection would be destroyed by gravity. The critical polar angle $\theta_{\rm c}$ is defined by

[TABLE]

When vertical convection is included in the NDAF model and for $\theta_{0}<\theta_{\rm c}$, Equation (78) should be rewritten as

[TABLE]

which means that the advection has been suppressed.

Figure 15 shows the variations of the half-opening angle $(\pi/2-\theta)$ with the dimensionless radius $r$ for different cases. The half-opening angle of the disk in the case excluding vertical convection is similar to the solutions in Liu et al. (2012a, 2013), which has a slight difference compared with the case including convection. As the cooling modes, advection and neutrino cooling dominate in the outer and inner region of the disk, respectively. Thus the effects of vertical convection mainly operate in the inner region until very near the BH. The vertically convective energy transfer can be effective to suppress radial advection in NDAFs, but the region dominated by vertical convection is deviated from the equatorial plane of the disk, thus it can be expected that the neutrino emission rate should increase slightly.

In the research of NDAFs in the vertical direction, the luminosity of neutrino annihilation $L_{\nu\bar{\nu}}$ can be roughly evaluated. We define the annihilation efficiency as $\eta\equiv L_{\nu\bar{\nu}}/L_{\nu}$, which satisfies $\eta\propto V_{\rm ann}^{-1}$ (see, e.g., Mochkovitch et al. 1993; Liu et al. 2010a, 2015a), where $V_{\rm ann}$ is the volume above the NDAF. We can estimate $V_{\rm ann}$ by integrating the region of $\theta<\theta_{0}$ and $r<r_{\rm out}$. Thus we can obtain the proportionality constant of the efficiency to the volume above the disk as a function of accretion rate from the vertically-integrated NDAF model where $\theta_{0}$ corresponds to $\pi/2$ (Liu et al. 2007, 2015a).

Figure 16 displays the neutrino luminosity $L_{\nu}$ and annihilation luminosity $L_{\nu\bar{\nu}}$ as a function of $\dot{m}$. The solid and dashed lines correspond to the cases including (Liu et al. 2015a) and excluding (Liu et al. 2010a, 2012a) the vertical convection, respectively. The thick and thin lines correspond to the neutrino luminosity and annihilation luminosity. The certain increase of luminosity comes from the suppressed advection and correspondingly increased neutrino cooling. We also notice that for $\dot{M}\sim 5~{}M_{\odot}~{}{\rm s}^{-1}$, the density of radiated neutrino is so large that the annihilation efficiency is close to 1. The extremely thick disk restrains the ejection in a very narrow empty funnel along the rotation axis. Similar to the discussion in Liu et al. (2010a, 2015a), thanks to the large opening angle of the disk and the sufficient output of energy, the neutrino annihilable ejection required by GRBs may be naturally realized by NDAFs.

4.7 Radiative transfer

At mentioned above, the polytropic relation should be substituted by the equations of neutrino radiative transfer. Sawyer (2003) firstly proposed the two-stream approximation to solve the neutrino transport problem in the analysis of some accretion disks. Rossi et al. (2007) derived a very simple form of neutrino radiative transfer, which can be used in the solutions of NDAFs.

Pan & Yuan (2012a) calculated the 1D Boltzmann equation of neutrino transfer in the accretion disk and obtained the neutrino spectra by given the distribution of density, temperature and chemical components of the disk. They also verified that the approaches of Sawyer (2003) is a good approximation. Pan & Yuan (2012b) further investigated the vertical structure, neutrino luminosity, and annihilation luminosity of NDAFs by considering neutrino radiative transfer, and they considered that the effects of the BH spin and magnetic field might be introduced in NDAF models to power GRBs.

5 Numerical simulations of NDAFs

The time-dependent 2D or 3D hydrodynamical simulations of hyperaccretion disks in BH-NS or NS-NS mergers and in collapsars have been widely studied (e.g., Ruffert & Janka 1999; Lee et al. 2004; Setiawan et al. 2006; Lee et al. 2009; Carballido & Lee 2011; Sekiguchi & Shibata 2011; Caballero et al. 2012; Fernández & Metzger 2013; Janiuk et al. 2013; Just et al. 2015; Fernández et al. 2015; Foucart et al. 2015; Batta & Lee 2016; Just et al. 2016).

Janiuk et al. (2013) calculated the structure and violent evolution of a turbulent torus accreting onto a BH. In the 2D simulations, neutrino cooling makes the disk much denser, geometrically thinner and less magnetized. For the accretion rate $\dot{m}\sim$ 0.03-0.1, the neutrino luminosity reaches $10^{53}-10^{54}$ erg s*-1*, which is 1-2 orders of magnitude larger than the BZ jet luminosity. The conclusion is similar to the results on the vertical structure and luminosity of NDAFs by Liu et al. (2010a, 2015a). Moreover, the neutrino cooled torus launches a fast, rarefied wind that is responsible for a powerful mass outflow, the power of which is correlated with the net accretion rate. The neutrino cooling rates are similar for the inner $\sim$ 20-30 $r_{\rm g}$ in the 1D and 2D calculations. Janiuk (2017) studied the general relativistic, MHD NDAF models including self-consistent nuclear EoS. The structure of inflows and outflows, neutrino radiation, and nucleosynthesis are studied. They argued that the central BH characteristics may be estimated by the observations on the decay of the radioactive elements in the inflows and outflows (e.g., Surman et al. 2011).

Different from the traditional simulation method in accretion disks, the 2D or 3D Lagrangian smooth particle hydrodynamics methods are also introduced in NDAFs (e.g., Lee et al. 2009; Batta & Lee 2016). Batta & Lee (2016) investigated the collapse and accretion onto BHs of spherically rotating envelopes by considering the angular momentum distribution. Contrary to the results obtained in previous 2D hydrodynamical simulations, they found that the collapsing gas with angular momentum between supercritical and subcritical causes the production of large quiescent times originating from the absence of an accretion disk near the BH. Moreover, the collapse of extremely subcritical materials on to the disk would result in a shutdown of the inner engine, if the critical angular momentum increases beyond the angular momentum of the materials in the disk before supercritical materials fill in.

6 Applications to GRBs

6.1 Jet luminosity

The most important motivation in study of NDAFs is to explain the source of GRB power. Two mechanisms, i.e., the neutrino annihilation process and BZ mechanism, have been widely discussed in the literatures. In the following, we confront both mechanisms against the GRB data.

6.1.1 Neutrino annihilation for SGRBs

As mentioned above, Popham et al. (1999) and Liu et al. (2007) investigated the spatial distribution of neutrino annihilation rate and found that most of the annihilation luminosity is ejected from the region $r\lesssim 20~{}r_{\rm g}$. In the studies on the vertical structure of NDAF model, Liu et al. (2010a, 2012a, 2013), found that the very large half-opening angle of the disk for a typical accretion rate can naturally constrain the neutrino annihilable ejection to produce the primary fireball of a GRB.

The observed fireball mean power outputting $\dot{E}$ from a central engine is a fraction of $L_{\nu\bar{\nu}}$, i.e.,

[TABLE]

where $\varepsilon$ is the conversion factor (e.g., Aloy et al. 2005; Fan & Wei 2011; Liu et al. 2012b). The output power can be estimated from the observation data, i.e.,

[TABLE]

where $z$ is the redshift, $E_{\rm\gamma,iso}$ is the isotropic radiated energy in the prompt emission phase, $E_{\rm k,iso}$ is the isotropic kinetic energy of the afterglow, $T_{\rm 90}$ is the duration of GRBs, which can be roughly considered as the duration of the the central engine activity, and $\theta_{\rm jet}$ is the half opening angle of the ejecta.

Hence, for the cases of $\dot{m}_{\rm ign}<\dot{m}<\dot{m}_{\rm trap}$ in the analytic formula of Zalamea & Beloborodov (2011), one has the mean accretion rate (Fan & Wei 2011; Liu et al. 2015b)

[TABLE]

where $E_{\rm k,iso,51}=E_{\rm k,iso}/(10^{51}~{}\rm erg)$, $E_{\rm\gamma,iso,51}=E_{\rm\gamma,iso}/(10^{51}~{}\rm erg)$, and $T_{90,\rm s}=T_{90}/(1~{}\rm s)$. The dimensionless mean disk mass derived from $m_{\rm disk}=\dot{m}T_{\rm 90,s}/(1+z)$, which reads (Liu et al. 2015b)

[TABLE]

Notice that $E_{\rm\gamma,iso}$ can be derived from the prompt emission data, and $E_{\rm k,iso}$ and $\theta_{\rm jet}$ can be deduced from afterglow modeling (e.g., Sari et al. 1998, 1999; Frail et al. 2001; Panaitescu & Kumar 2002; Lloyd-Ronning & Zhang 2004; Zhang et al. 2007; Fong et al. 2012; Wang X et al. 2015).

According to the above equations, one can estimate the disk mass by using the observational data. Liu et al. (2015b) collected 31 SGRBs with the authentic X-ray detections and known redshifts discovered by Swift and HETE-2. Figure 17 displays the distributions of the disk masses $m_{\rm disk}$ with $\varepsilon=0.3$ for different typical BH masses and spins. One can see that the spin parameters are more effective than BH mass on the values of disk mass. The disk masses of most SGRBs are below $0.2-0.4~{}M_{\odot}$, and occasionally reach the limit of $0.5~{}M_{\odot}$. Even for the extreme case of $m=3$, $a_{*}=0.9$, there still exists one SGRB, whose disk mass is larger than $0.45~{}M_{\odot}$. These cases with massive disks may point towards an origin that invoke a massive star (Type II) rather than a compact stars merger (Type I) (e.g., Zhang et al. 2009; Lazzati et al. 2010).

6.1.2 Neutrino annihilation for LGRBs

Song et al. (2016) collected the isotropic gamma-ray radiation energy and jet kinetic energy of 48 LGRBs with known redshifts. Similar to Liu et al. (2015b), they estimated the mean accreted masses $m_{\rm acc}$ (instead of $m_{\rm disk}$ in SGRBs) of LGRBs in the sample to investigate whether NDAFs can power LGRBs with reasonable BH parameters and conversion factor $\varepsilon$. Figure 18 shows the distributions of the accreted masses for different typical mean BH masses and spins and conversion factors. Since the BHs in the centre of collapsars should be rotating very rapidly, the spin parameter is set to be larger than 0.9 (e.g., MacFadyen & Woosley 1999; Popham et al. 1999; Zhang et al. 2003; Woosley & Bloom 2006). If one considers $a_{*}$ = 0.998, $m=3$, and $\varepsilon$ = 0.5, all the values of the accreted masses are less than $7~{}M_{\odot}$. Obviously, if the accreted mass is larger than $7~{}M_{\odot}$, the mean BH mass must be much larger than $3~{}M_{\odot}$. As a result, most of the values of the accreted masses are less than that from the collapsar simulation, $\sim 5~{}M_{\odot}$, for extreme Kerr BHs and a high conversion factor. It suggests that the NDAFs may be suitable for most LGRBs except for some extremely energetic sources.

In the above studies, the effects of the outflow are ignored (Liu et al. 2012a; Janiuk et al. 2013). Including them may significant influence the disk mass of LGRBs even not so much for SGRBs. Furthermore, X-ray flares (e.g., Burrows et al. 2005; Chincarini et al. 2007; Falcone et al. 2007) and the shallow decay phase (e.g., Zhang et al. 2006) as the common feature in GRB afterglows have not been considered in this framework, which may bring more challenges to NDAFs.

6.1.3 Evolution of NDAFs for GRBs

High accretion rate of NDAFs should trigger violent evolution of BH characteristics, which further leads to the evolution of the neutrino annihilation luminosity (Janiuk et al. 2004; Song et al. 2015). The evolution equations of a Kerr BH, based on the conservation of energy and angular momentum, can be expressed by (e.g., Liu et al. 2012b; Song et al. 2015)

[TABLE]

[TABLE]

where $J=a_{*}GM^{2}/c$ is the angular momentum of the BH, $e_{\rm ms}$ and $l_{\rm ms}$ are the specific energy and angular momentum corresponding to the marginally stable orbit radius $r_{\rm ms}$ of the disk, i.e., (e.g., Novikov & Thorne 1973; Wu et al. 2013; Hou et al. 2014b)

[TABLE]

[TABLE]

Therefore the evolution of the BH spin is expressed by

[TABLE]

According to the above equations, one can obtain the characteristics of the BH if the initial mass $M_{0}$ and spin $a_{0}$ of the BH are given. Similar to the above discussion, one can calculate along with the evolution of the BH mass and spin, the relations between the neutrino annihilation energy $E_{\nu\bar{\nu}}$ and $T_{90}$ by using the analytic formula of Zalamea & Beloborodov (2011) with $\dot{m}_{\rm ign}<\dot{m}<\dot{m}_{\rm trap}$, i.e.,

[TABLE]

where $t_{\rm s}=t/(1~{}\rm s)$.

Figure 19(a) displays the individual neutrino annihilation energy of the SGRB data in Liu et al. (2015a) compared with the typical theoretical lines. All curves correspond to the initial value of BH mass $m_{0}=2.3$, conversion factor $\varepsilon=0.3$, and typical redshift $z=0.5$, and all the lines are truncated at the disk mass $m_{\rm disk}=0.5$. We notice that all the points are under our predicted lines, which means that NDAFs can interpret the profiles in SGRB sample even considering BH evolution.

Figure 19(b) shows a comparison between the predictions of NDAFs and the data of LGRBs in Nemmen et al. (2012) with parameters $m_{0}=3$, $\varepsilon=0.3$, and $z=2$. The duration of LGRBs shorter than 2 s or longer than 300 s are not considered (Leng & Giannios 2014; Liu et al. 2015c). All lines are truncated when $m_{\rm disk}=5$. In this figure, more than half of the LGRB data are below our predicted lines, which means that the neutrino annihilation processes can power these LGRBs with reasonable initial conditions. Other LGRB data still exceed the range of the lines, suggesting that alternative MHD processes may be required.

6.1.4 Neutrino annihilation versus BZ mechanism

Since the neutrino annihilation process cannot explain all the observed GRBs, BZ mechanism should be further studied. Kawanaka et al. (2013b) investigated the BZ jet luminosity and efficiency expected from BH-NDAF systems, and obtain the analytic descriptions of BZ luminosity and compared it with the neutrino annihilation luminosity.

Liu et al. (2015b) used the analytic formula of $p_{\rm in}$ from Xue et al. (2013),

[TABLE]

and combined with the BZ luminosity described by Equations (52) and (53) by ignoring the effects of the magnetic field configuration to estimate the BZ and neutrino annihilation luminosities as the functions of the disk masses and BH spin parameters, and contrasted the observational GRB jet luminosities. As the results, the BZ mechanism is more effective than the neutrino annihilation processes, especially for LGRBs. Actually, if the energy of afterglows and flares is included, the distinction between these two mechanisms is more significant. Mounting evidence suggests that at least for some GRBs outflow carries significant Poynting flux: missing or weak thermal components (Zhang & Pe’er 2009), strong linear polarization in gamma-rays (Fan et al. 2005; Lai 2015; Yonetoku et al. 2011) and early afterglows (Steele et al. 2009; Mundell et al. 2013), stringent upper limits of neutrino flux (Zhang & Kumar 2013), and bulk acceleration in GRBs and X-ray flares (Uhm & Zhang 2016a, b). Future GRB polarization observations by the POLAR detector may give more evidence and further distinguish these two mechanisms.

As the summary and comparison, Figure 20 shows the applications of three central engine mechanisms, i.e., neutrino annihilation, BZ jet, and magnetar, to GRBs and flares. Most of SGRBs, about half of LGRBs, and non ULGRBs (e.g., Virgili et al. 2013; Levan et al. 2014; Zhang et al. 2014) can be explained by the neutrino annihilation (e.g., Nathanail & Contopoulos 2015). If we consider that some flares originate from central engine with low accretion rates and low disk masses (e.g., Bernardini et al. 2011; Margutti et al. 2011; Mu et al. 2016), most of flares can be covered by the power of the neutrino annihilation, however, ultra-long flares (UL flares) lasting about tens of thousand seconds may come from magnetars (Mu et al. 2016).

6.2 Variability

Rapid variability has been observed in GRB prompt emission lightcurve. Variability may originate from internal shock model (e.g., Kobayashi et al. 1997), relativistic mini-jets (e.g., Lyutikov & Blandford 2003; Yamazaki et al. 2004; Zhang & Zhang 2014), or relativistic turbulence (Narayan & Kumar 2009; Kumar & Narayan 2009; Lazar et al. 2009; Lin et al. 2013). More fundamentally variability arises from the central engine, which may be modified as the jet interacts with the stellar envelope (Morsony et al. 2010).

Variability can naturally arise from the hyperaccreting central BH engine. The thermal and viscously instabilities in NDAFs have been widely studied (e.g., Janiuk et al. 2004, 2007; Lee et al. 2005; Lei et al. 2009; Kawanaka et al. 2013a; Kimura et al. 2015; Xie et al. 2016) to explain the variability. Carballido & Lee (2011) characterized the time variability of energy release at small scales in NDAFs through shearing box MHD simulations. Masada et al. (2007) investigated the magnetorotational instability in the hyperaccretion disk. Once the dead zone gains a large amount of mass and becomes gravitationally unstable, the episodically intense mass accretion can cause short-term variabilities. Kawanaka & Kohri (2012) found that convective motion in the vertical direction can also trigger sporadic mass accretion, which causes variability in GRBs. Differently, Lin et al. (2016) showed that the propagating fluctuations mechanism (e.g., Lyubarskii 1997; King et al. 2004; Lin et al. 2012) that drives variabilities in BH binaries and AGNs may be responsible for the observed variability in prompt emission.

6.3 Jet precession

Another area of research in NDAF-GRB connection is the possibility of jet precession (e.g., Blackman et al. 1996; Portegies Zwart et al. 1999; Reynoso et al. 2006; Lei et al. 2007; Liu et al. 2010b; Stone et al. 2013). Blackman et al. (1996) first investigated a relativistic precessing blob-emitting NS jet in a binary pulsar. They considered Newtonian tidal torque and gravitomagnetic interaction between pulsar binary to account for precession and nutation. Portegies Zwart et al. (1999) proposed that the BH forces the accretion disk and jet to precess due to the Newtonian tidal torque, and fitted the observational data. Since the gravitomagnetic interaction between the BH and disk is, however, much stronger than the tidal force (Thorne et al. 1986), Reynoso et al. (2006) suggested that the gravitomagnetic interaction can precess jet from the central engine of GRBs. Lei et al. (2007) considered the effects of the half-opening angle of a precessing jet on shaping GRB lightcurves. In their model, since the whole disk precesses around the BH (i.e., the angular momentum of the disk is less than that of the BH), the required size of the disk needs to be small enough.

Naturally, the angular momentum of the disk is much larger than that of the BH in BH-NDAF systems, thus a new jet precession model should be built (Liu et al. 2010b; Sun et al. 2012). The basic picture (see Figure 21) is that a Kerr BH surrounded by a tilted NDAF whose initial orbital axis is misaligned with the BH spinning axis. In this framework, the outer disk whose angular momentum is sufficiently larger than that of the BH can maintain its orientation and force the BH to precess. Meanwhile, the inner disk whose angular momentum is significantly smaller than that of the BH, should be aligned with the BH spin axis (Bardeen & Petterson 1975).

The angular momentum per ring at radius $r$ with width $dr$ is $dJ_{\rm disk}=2\pi r^{2}\Sigma v_{\varphi}dr$, so the typical angular momentum of the disk is (e.g., Sarazin et al. 1980)

[TABLE]

Due to the Lense-Thirring effect (Lense & Thirring 1918), there exists a critical radius $r_{\rm{p}}$ where the typical angular momentum $J_{\rm disk}|_{r=r_{\rm p}}$ is equal to the BH angular momentum $J$, i.e.,

[TABLE]

The outer disk ($r>r_{\rm{p}}$) will maintain its orientation and therefore force the BH and the inner disk ($r<r_{\rm{p}}$) to be a whole precessing system. The precession rate is expressed as (e.g., Sarazin et al. 1980; Lu 1990)

[TABLE]

According to the mass conservation equation, the precession period can be derived as

[TABLE]

Jet precession model can also interpret some lightcurves that show $E_{p}$ tracking behavior. Usually a GRB pulse shows a fast-rise-exponential-decay shape with the $\nu F_{\nu}$ spectrum $E_{p}$ evolves from hard to soft. This can be naturally explained by synchrotron radiation of an expanding shell (Uhm & Zhang 2014, 2016b). However, some well-separated GRB pulses show a symmetric structure, and their $E_{p}$ traces the lightcurve behavior (e.g., Liang & Kargatis 1996; Liang & Nishimura 2004; Lu et al. 2010, 2012; Liang et al. 2015). Both the temporal and spectral properties of these symmetric pulse are difficult to be explained with internal shocks. Figure 22 shows the flux $F$ and $E_{p}$ in jet precession model as compared with the observations (Liu et al. 2010b).

Using the above equations and the analytic formulae of the structure of NDAFs derived from Popham et al. (1999), the analytic expression of the precession period $P$ can be estimated by

[TABLE]

Assuming $\alpha$ as a constant, the time derivative of $P$ is expressed as

[TABLE]

These equations suggest a quasi-periodic variability and its evolution of GRBs. Hou et al. (2014a) exhibited that a $86^{+5.9}_{-9.4}~{}\rm s$ periodic oscillation may exist in the data from about 5300 s to about 6100 s in the bump of GRB 121027A using the stepwise filter correlation method (Gao et al. 2012) and the Lomb-Scargle method (Scargle 1982), which can be interpreted by jet precession model (Liu et al. 2010b). Hou et al. (2014b) discovered that there is a remarkable time evolution in the data of the flares in GRB 130925A. It can be also explained by Equation (104) with reasonable initial BH mass and spin.

6.4 Nucleosynthesis

The central engine of GRBs is an ideal location to supply an extremely hot and dense environment to produce heavy elements. The radial and vertical components of NDAFs have been studied in detail in Liu et al. (2013) and Xue et al. (2013), which can be connected to the origin of strong Fe K$\alpha$ emission lines and the SN bump observed in long GRB afterglow lightcurves (see Figure 23).

The SN bump may be driven by the decay of ${}^{56}\rm Ni$ (e.g., Galama et al. 1999; Woosley & Bloom 2006), but the product of massive ${}^{56}\rm Ni$ in SNe accompanied GRBs remains unsolved. The possibilities including ${}^{56}\rm Ni$ being originated from the central engine transported by the outflow (e.g., MacFadyen & Woosley 1999; MacFadyen 2003; Surman et al. 2011) or due to explosive burning (e.g., Maeda & Nomoto 2003; Maeda & Tominaga 2009). Reeves et al. (2002) reported on an XMM-Newton observation of the X-ray afterglow of GRB 011211, which spectrum reveals evidence for emission lines of Magnesium, Silicon, Sulphur, Argon, Calcium, and possibly Nickel, arising in enriched material with an outflow velocity of order $\rm 0.1c$. However, there are no reported metal lines in GRBs from Swift XRT observatory.

The production processes of the heavy nuclei from the central engine of GRBs are widely studied recently. As shown in Figures 14 (c) and (d), for the low accretion rate of NDAFs exactly corresponding to LGRBs that are associated with SNe, ${}^{56}\rm Ni$ dominates near the disk surface (Liu et al. 2013). The decay of ${}^{56}\rm Ni$ in the outflows can produce the optical lightcurve SN bump. Surman et al. (2011) focused on nucleosynthesis, particularly for the progenitor of ${}^{56}\rm Ni$ in the hot outflows from GRB accretion disks. Metzger et al. (2011) suggested that the composition of ultrahigh energy cosmic rays becomes dominated by heavy nuclei at high energies forming GRB jets or outflows (also see e.g., Sigl et al. 1995; Horiuchi et al. 2012; Hu 2015).

6.5 Extended emission

Significant extended emission (EE) after prompt emission of SGRBs has been observed up to $\sim$ 100 s by CGRO/BATSE (Norris et al. 1995; Lazzati et al. 2001; Connaughton 2002), BeppoSAX (Montanari et al. 2005), and Swift/BAT (Norris et al. 2010). The most prominent case is GRB 060614 lasting $\sim 110$ s, whose lightcurve is composed of some initial hard spikes and a long soft gamma-ray tail (Gehrels et al. 2006). There is no accompanied SNe detected for this nearby GRB (Della Valle et al. 2006; Fynbo et al. 2006; Gal-Yam et al. 2006), disfavoring the collapse of a massive star as the progenitor (Zhang et al. 2007).

A magnetar model is proposed to explain the origin of SGRBs with EE (e.g., Bucciantini et al. 2012; Gompertz et al. 2013, 2014). Barkov & Pozanenko (2011) suggested that BZ jet and neutrino annihilation may coexist in BH-NDAF system, and GRBs with and without EE may depend on observer’s viewing angle from the jet axis.

In the framework of the BH-NS merger models, Liu et al. (2012b) studied the roles of radial angular momentum transfer in the disk and time evolution of the BH mass and spin as shown in Equations (92) and (93), then calculated the magnetic barrier around the BH, as Equation (53), to estimated the timescale to separate prompt emission and EE. The basic picture is shown in Figure 24, i.e., the radial angular momentum transfer may significantly prolong the lifetime of the accretion process and multiple episodes may be switched by the magnetic barrier. Our numerical calculations suggest that the model can fit most data of SGRBs with EE with a reasonable disk mass and BH characteristics.

6.6 X-ray flares

X-ray flares (e.g., Burrows et al. 2005; Chincarini et al. 2007; Falcone et al. 2007; Margutti et al. 2011) are common in GRBs and occur at times well after the initial prompt emission, with a time lag of the order of hundreds or thousands of seconds, which might be related to the activities of the central engine (e.g., Fan et al. 2005; Lazzati & Perna 2007; Maxham & Zhang 2009; Luo et al. 2013; Yi et al. 2016). Several mechanisms have been proposed to generally explain the episodic X-ray flares, including gravitational instability in the hyperaccretion disk (Perna et al. 2006; Liu et al. 2014b), fragmentation of a rapidly rotating core (King et al. 2005), a magnetic switch of the accretion flow (Proga & Zhang 2006; Liu et al. 2012b; Cao et al. 2014), differential rotation in a post-merger millisecond pulsar (Dai et al. 2006), transition from a thin to a thick disk (Lazzati et al. 2008), He-synthesis-driven wind (Lee et al. 2009), instability in the jet (Lazzati et al. 2011), outflow caused by the maximal and minimal possible mass accretion rates at each radius of NDAFs (Liu et al. 2008), episodic jets produced by the magnetohydrodynamic mechanism from the disk (Yuan & Zhang 2012), and MC-NDAFs (Luo et al. 2013).

Mu et al. (2016) investigated UL flares with duration $\gtrsim 10^{4}$ s, and argued that the corresponding central engine may not be associated with BHs, but a fast rotating NS with strong dipolar magnetic fields as shown in Figure 20.

6.7 Constraint on the mass and spin of central engine BHs in GRBs

BHs are mysterious and fascinating compact objects, which are generally related to multiband electromagnetic radiation, gravitational waves (GWs), neutrino emission, and cosmic rays. Two essential properties of BHs, i.e., mass and spin, are not easy to measure. Some dynamical or statistical methods have been introduced to constrain these parameters of supermassive BHs (e.g., Natarajan & Pringle 1998; Brenneman & Reynolds 2006; Tchekhovskoy et al. 2010; Lei & Zhang 2011; Kormendy & Ho 2013; Wang et al. 2013) and stellar-mass BHs (e.g, Bahcall 1978; Zhang et al. 1997; McClintock et al. 2014).

In NDAFs, neutrinos can tap the thermal energy gathered by the viscous dissipation and liberate a tremendous amount of binding energy, and neutrino-antineutrino annihilation above the disk would launch a hot fireball. A GRB powered by this mechanism is therefore thermally dominated. The launch site of the fireball should be above the typical neutrino annihilation radius, which can be constrained by the observed thermally spectral component (Pe’er et al. 2007).

From the analytic solutions in Xue et al. (2013), we obtain that the dimensionless annihilation height $h=H_{\rm ann}/r_{\rm g}$ is

[TABLE]

Combing with Equation (49), using the GRB luminosity and fireball launch radius, one can constrain the mass and spin of a GRB central engine BH. As shown in Figure 25, Liu et al. (2016b) estimate the following constraints on the central engine BH of GRB 101219B: mass $M_{\rm BH}\sim 5-9~{}M_{\odot}$, spin parameter $a_{*}\gtrsim 0.6$, and disk mass $3~{}M_{\odot}\lesssim M_{\rm disk}\lesssim 4~{}M_{\odot}$. The results also suggest that the NDAF model is a competitive candidate for the central engine of GRBs with a strong thermal component.

6.8 Kilonovae

BH-BH, BH-NS and NS-NS mergers are the strong GWs transient sources, which may be associated with short GRBs (e.g., Paczyński 1986; Eichler et al. 1989; Narayan et al. 1992), short-GRB-less X-ray transients (Zhang 2013; Sun et al. 2017), kilonovae/mergernovae (e.g., Li & Paczyński 1998; Metzger et al. 2010; Yu et al. 2013), even fast radio bursts (e.g., Liu et al. 2016a; Zhang 2016; Wang et al. 2016).

(Li & Paczyński 1998) first predicts a type of optical transient lasting a few days, which is powered by radioactive decay of the neutron-rich ejecta from an NS-NS or BH-NS merger system. The kilonova signature has been investigated in detail recently (e.g., Metzger et al. 2010; Barnes & Kasen 2013; Metzger & Fernández 2014; Kasen et al. 2015; Metzger 2017), and the enhanced merger nova signature due to energy injection of a post-merger magnetar engine has been studied in detail (e.g., Yu et al. 2013, 2015; Metzger & Piro 2014; Li & Yu 2016; Gao et al. 2017). Observationally, several kilonova/mergernova candidates have been discovered in association with some short GRBs, including typical kilonovae in GRB 130603B (Tanvir et al. 2013; Berger et al. 2013), GRB 060614 (Jin et al. 2015; Yang et al. 2015), and GRB 050704 (Jin et al. 2016), and a few magnetar-powered merger nova (Gao et al. 2016, 2017). The applicability of NDAFs for some of these GRBs have been studied (Shen et al. 2017).

For NDAF models, Yi et al. (2017) studied the effects of vertical advection (Jiang et al. 2014) on the structure and luminosity of NDAFs. Due to the strong cooling through the vertical advection, the neutrino luminosity from the disk is lower than that of NDAFs without vertical advection. Interestingly, the gamma-ray emission from the disk surface can be extremely super-Eddington, which may have significant contribution to the kilonovae as shown in Figure 26. Following radioactive decay from the outflow of the disk, this process may power a transient even in the direction off the jet. Song & Liu (2017) proposed that a BH hyperaccretion disk with a strong outflow (about 99$\%$ accretion matter) via radioactive decay may trigger a kilonova-like transient as powerful as the one powered by a magnetar.

6.9 Progenitor stars

For the same settings of the accretion rate and BH characteristics, the BH hyperaccretion inflow-outflow model including BZ mechanism is more powerful than the NDAF model, which can be used to constrain the progenitor stars of GRBs (Liu et al. 2017b; Song & Liu 2017).

Considering outflows from an NDAF, Liu et al. (2017b) showed tight constraints on the properties of progenitor stars from long GRB observations. Specifically, only the solar-metallicity, massive stars or parts of zero-metallicity stars can be as the progenitors of LGRBs lasting from several seconds to tens of seconds in rest frame. The fraction of bursts in the LGRB population is more than 40$\%$, which cannot be accounted for by these rare required progenitors. It implies that the activity timescale of central engine may be much longer than the duration of prompt gamma-ray emission, as indicated by the extended X-ray activities observed from the Swift data (Zhang et al. 2014; Lü et al. 2014). Alternatively, LGRBs may be powered by the magnetars rather than NDAF systems.

7 Multi-messenger signals from NDAFs

There may be two observational messenger signals to probe the invisible central engine of GRBs and test the existence of NDAFs, i.e., MeV neutrinos and GWs.

7.1 Detectable MeV neutrinos from NDAFs

Whereas the leading candidates of MeV neutrino sources are SNe similar to SN 1987A. NDAFs around rotating BHs or NSs can be another source of cosmic MeV neutrinos. Assuming a distance of 10 kpc, the detectability of a nominal NDAF with a large accretion rate by Super-Kamiokande was discussed (e.g., Nagataki & Kohri 2002; Caballero et al. 2012, 2015).

Based on the solutions of Xue et al. (2013), Liu et al. (2016c) derived the fitting formulae for the mean cooling rate due to electron neutrino and antineutrino losses, $Q_{\nu_{\rm e}}$ and $Q_{\bar{\nu}_{\rm e}}$, and the temperature of the disk $T$, as the functions of the mean BH spin parameter, mean accretion rate, and radius ($m=3$ adopted), i.e.,

[TABLE]

[TABLE]

[TABLE]

According to the above equations, Liu et al. (2016c) calculated the electron neutrino and antineutrino spectra of NDAFs by fully taking into account the general relativistic effects by using the null geodesic equation (Carter 1968), and studied the effects of viewing angle, BH spin, and accretion rate on the spectra. Figure 27 shows that even though a typical NDAF has a neutrino luminosity lower than that of a typical SN, it can reach $10^{50}-10^{51}~{}{\rm erg~{}s^{-1}}$ peaking at $\sim 10-20$ MeV, making them potentially detectable with the upcoming sensitive MeV neutrino detectors if they are close enough to Earth. Based on the detailed star formation rate and metallicity in Local Group, we estimate a detection rate up to $\sim$ (0.10-0.25) per century for GRB-related NDAFs by the Hyper-Kamiokande (Hyper-K) detector (Abe et al. 2011; Mirizzi et al. 2016), if one neglects neutrino oscillation and demands at least 3 neutrinos being detected. If one assumes that all Type Ib/c SNe have NDAFs, the Hyper-K detection rate would be $\sim$ (1-3) per century. By considering neutrino oscillations, the detection rate may decrease by a factor of 2-3. By the way, it can be expected that neutrinos from NDAFs may slightly effect on the universal neutrino background.

7.2 GWs from NDAFs

One of the methods to infer the existence of NDAFs may be through detections of the GW signals from NDAFs due to precession (Romero et al. 2010; Sun et al. 2012) or anisotropic neutrino emission (Epstein 1978; Sago et al. 2004; Hiramatsu et al. 2005; Suwa & Murase 2009; Kotake et al. 2012) in the disk. Liu et al. (2017a) compared the detectabilities of the GWs from NDAFs, BH hyperaccretion disks with BZ jets, and magnetars.

The jet precession mechanism introduced by Liu et al. (2012b) leads to a non-axisymmetric mass distribution, which can cause gravitational radiation. Sun et al. (2012) estimated that the typical quadrupole power of NDAFs is about $10^{44}~{}\rm erg~{}s^{-1}$ at 10 Hz, which is much less than the typical GRB luminosity. According to Figure 28, the GWs from NDAFs caused by jet precession may be detectable by the future GW detectors if the sources are in the Local Group. GW power and frequency from NDAFs may be different from these from compact objects mergers or collapsars (e.g., Nakar & Piran 2011; Ott et al. 2011; Gao et al. 2013; Messenger et al. 2014) and relativistic jet (e.g., Akiba et al. 2013; Birnholtz & Piran 2013), but detailed studies are needed to model the wave forms of these events.

Two confirmed GW events and another candidate have been discovered by LIGO (Abbott et al. 2016a, b). The Fermi/GBM recorded a possible gamma-ray transient 0.4 s after GW 150914 (Connaughton et al. 2016). Several models have been proposed to explain the possible electromagnetic counterpart of GW 150914 (e.g., Li X et al. 2016; Zhang 2016; Loeb 2016; Liu et al. 2016a; Perna et al. 2016; Woosley 2016; Zhang et al. 2016; Janiuk et al. 2017). If such an association is true, the combined GW and electromagnetic signals can be used to test Einstein’s Equivalence Principle (Wu et al. 2016). Detecting NDAF-related GWs and the electromagnetic counterparts can also serve the purpose.

8 Summary and discussion

As a theoretical concept, NDAFs are studies within the context of accretion systems with extremely high accretion rates such as cosmological GRBs. NDAFs are the naturally extension of slim disks for extreme accretion rates. As the main cooling mechanism, neutrino cooling process leads to a copious neutrinos that may power GRBs via neutrino annihilation. Due to the extreme density, temperature involved in such a system with extreme accretion rate, rich physics (gravitational, thermal dynamical, nuclear, and particle physics) is involved. We have reviewed the recent progress in studying the physical processes and properties of NDAFs, including the steady radial and vertical structure of NDAFs and the implications for calculating neutrino luminosity and annihilation luminosity, jet power due to neutrino-antineutrino annihilation and BZ mechanism and their dependences on parameters such as BH mass, spin, and accretion rate, time evolution of NDAFs, and effect of magnetic fields. The applications of NDAF theories to the GRB phenomenology such as lightcurve variability, precession, extended emission, X-ray flares, etc., are reviewed, and possible probes NDAFs using multi-messenger signals such as MeV neutrinos and gravitational waves are discussed.

Since GRB observations are directly related to the relativistic jets rather than the central engine, directly inferring the existence of NDAFs as well as their properties are not easy. The uncertainties lie in the unknown composition, energy dissipation mechanism, and particle acceleration and radiation mechanism of the jet (Zhang 2011; Kumar & Zhang 2015, for recent reviews). Furthermore, the identity of the central engine, be it a hyperaccreting BH or a spinning down millisecond magnetars, is still uncertain. It is possible and even likely that both types of engines are operating in different GRBs. In any case, NDAFs must exist in BH systems, and may also exist and operate in magnetar systems as well. For BH engines, the BZ-dominated systems tend to produce Poynting-flux-dominated outflows whereas the neutrino annihilation dominated systems tend to produce hot fireballs. GRB spectral observations seem to suggest that both types of jet composition may exist in different GRBs (e.g., Abdo et al. 2009a, b; Ryde et al. 2010; Zhang 2011; Guiriec et al. 2011, 2013; Axelsson et al. 2012; Zhang & Pe’er 2009; Gao & Zhang 2015). For thermally-dominated fireballs, the connection between the central engine and observational properties is closer, and interesting constraints can be placed to the properties of the BH and NDAFs (Liu et al. 2016b). Future more sensitive GRB spectral detectors and hopefully multi-messenger observations of some nearby NDAF sources may shed light into more detailed information about NDAFs and may even directly infer its existence.

Acknowledgements

We thank Zi-Gao Dai, En-Wei Liang, Yi-Zhong Fan, Xue-Feng Wu, Da-Bin Lin and Li Xue for beneficial suggestions and discussion. This work was supported by the National Basic Research Program of China (973 Program) under grant 2014CB845800, the National Natural Science Foundation of China under grants 11473022, 11333004, and 11573023.