Neutrino-Dominated Accretion Flows with Magnetic Prandtl Number-Dependent MRI-driven Turbulence

TL;DR

This paper explores how magnetic Prandtl number-dependent turbulence affects the stability of neutrino-dominated accretion flows, revealing potential causes for gamma-ray burst variability and implications for neutron star merger emissions.

Contribution

It introduces a variable-$\\alpha$ model based on magnetic Prandtl number dependence, providing new insights into NDAF stability and variability mechanisms.

Findings

NDAFs can become highly unstable with variable-$\alpha$ due to Pm dependence.

Unstable NDAFs exhibit significant time variability in accretion rates.

The instability mechanism may explain rapid GRB prompt emission variability.

Abstract

We investigate the stability of a neutrino-dominated accretion flow (NDAF), which is expected to be formed in the gravitational collapse of a massive star or the merger of a neutron star binary, based on the variable-$\alpha$ prescription. Recent magnetohydrodynamic (MHD) simulations shows that the viscosity parameter $\alpha$ is proportional to the power of the magnetic Prandtl number ${\rm Pm}=\nu/\eta$, where $\nu$ and $\eta$ are the kinematic viscosity and electric resistivity of the fluid, respectively. In the inner region of a hyperaccretion flow, the viscosity and resistivity are carried by mildly, relativistically degenerated electrons. We fit the dependence of the magnetic Prandtl number on density and temperature by a simple analytic form, and derive the condition for an NDAF to be dynamically unstable. As demonstrations we perform simple one-dimensional simulations of NDAFs with the variable-$\alpha$ and show that the mass accretion becomes highly time-variable in the unstable branch. This mechanism may account for the rapid variability observed in the prompt emission of gamma-ray bursts (GRBs). The mass ejection from a hyperaccretion flow due to viscous heating, which makes a kilonova/macronova emission in the merger of a neutron star binary, is also briefly discussed.

Equations

$$\dot{M}=-2\pi r{\text{varv}}_r\Sigma ,$$

$$\Omega =\sqrt{\frac{G{M}_{\mathrm{BH}}}{r^3}},$$

$$\frac{p}{\rho }={\Omega }^2{H}^2,$$

$$-{\nu }_t\Sigma \frac{d\Omega }{d\ln r}=\frac{\dot{M}\Omega }{2\pi }\left(1-\sqrt{\frac{R_{\mathrm{in}}}{r}}\right),$$

$$Q^{+}\equiv {\nu }_t\Sigma {\left(\frac{d\Omega }{d\ln r}\right)}^2=Q^{-},$$

$$-{\nu }_t\Sigma \frac{d\Omega }{d\ln r}=\alpha pH.$$

$$Q^{-}=Q_{\mathrm{adv}}^{-}+Q_{\nu }^{-},$$

$$q_{eN}$$

$$q_{e^{-}{e}^{+}\to {\nu }_x{\bar{\nu }}_x}=1.4\times 10^{33}\mathrm{erg}{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}{T}_{11}^9,$$

$$Q_{\nu }^{-}=2\sum _i\frac{(7/8)\sigma T^4}{(3/4)({\tau }_{{\nu }_i}/2+1/\sqrt{3}+1/(3{\tau }_{a,{\nu }_i}))},$$

$${\tau }_{a,{\nu }_i}$$

$${\tau }_{s,{\nu }_i}\simeq 2.7\times 10^{-7}{T}_{11}^2{\rho }_{10}H,$$

$$p\simeq \frac{11}{12}a{T}^4+\frac{\rho k_{B}T}{m_p}+\frac{2\pi hc}{3}{\left(\frac{3}{8\pi m_p}\right)}^{4/3}{\left(Y_e\rho \right)}^{4/3}+\frac{u_{\nu }}{3},$$

$$p$$

$$Q^{-}$$

$$Q^{-}\approx 7.4\times 10^{46}\mathrm{erg}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}{T}_{11}^2{\rho }_{10}^{-1}{H}^{-1},$$

$$\frac{\partial \left(T_{r\varphi }H\right)}{\partial \Sigma }<0,$$

$$\alpha \propto T^5,$$

$$T\propto H^2,$$

$$\frac{\partial \left({\alpha }^{6/5}\Sigma \right)}{\partial \Sigma }<0,$$

$$\frac{\partial \left({\alpha }^2{\Sigma }^3\right)}{\partial \Sigma }<0,$$

$$\mathrm{Pm}\simeq 4.5\times 10^2{T}_{11}^2{\rho }_{10}^{-1}.$$

$$\alpha \propto {\rho }^{-\beta }{T}^{2\beta }.$$

$$\frac{10}{17}<\beta <2,$$

$$\beta >\frac{2}{5}.$$

$$\frac{\partial \Sigma }{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[r^{1/2}\frac{\partial }{\partial r}\left({\nu }_t\Sigma r^{1/2}\right)\right],$$

$$\frac{\partial \ln T}{\partial t}+v_r\frac{\partial \ln T}{\partial r}$$

$$v_r=-\frac{3}{\Sigma r^{1/2}}\frac{\partial }{\partial r}\left[{\nu }_t\Sigma r^{1/2}\right].$$

$$\alpha (\mathrm{Pm})={\alpha }_{\min }+({\alpha }_{\max }-{\alpha }_{\min })\left(\frac{{\mathrm{Pm}}^{\beta }}{{\mathrm{Pm}}^{\beta }+{\mathrm{Pm}}_c}\right).$$

$$\mathcal{Q}\equiv \frac{c_s\Omega }{\pi G\Sigma }<{\mathcal{Q}}_{\mathrm{crit}},$$

$$\frac{\partial {\alpha }_{\mathrm{grav}}}{\partial t}=-\sigma \frac{{\alpha }_{\mathrm{grav}}^2-{\alpha }_0^2}{{\alpha }_{\mathrm{grav}}^2+{\alpha }_0^2}+{\nu }_g\frac{{\partial }^2{\alpha }_{\mathrm{grav}}}{\partial r^2},$$

Full text

1. Introduction

It is widely believed that gamma-ray bursts (GRBs) are powered by hyperaccreting black holes (or neutron stars, possibly) which are formed after the gravitational collapse of a massive star or the compact binary (neutron star - neutron star, or neutron star - black hole) merger (Narayan et al., 1992; Popham et al., 1999; Narayan et al., 2001). The accretion rate would be as high as $0.01-1M_{\odot}~{}{\rm s}^{-1}$, and due to its high density ($\sim 10^{9-10}~{}{\rm g}~{}{\rm cm}^{-3}$) and temperature ($\sim 10^{9-11}~{}{\rm K}$), the accretion flow is mainly cooled via neutrino emission. It is often called a neutrino-dominated accretion flow (NDAF), and its structure and behavior have been investigated by many researchers (Kohri & Mineshige, 2002; Di Matteo et al., 2002; Kohri et al., 2005; Gu et al., 2006; Chen & Beloborodov, 2007; Kawanaka & Mineshige, 2007; Liu et al., 2007; Komissarov & Barkov, 2009; Lei et al., 2013; Kawanaka et al., 2013; Xue et al., 2013; Liu et al., 2015; Lei et al., 2017; Janiuk, 2017; for the recent comprehensive review see Liu et al., 2017). In this model, it is considered that an ultrarelativistic jet is launched from this accretion system due to the energy deposition process such as neutrino pair annihilation (Eichler et al., 1989; Asano & Fukuyama, 2000; Zalamea & Beloborodov, 2011; Suwa, 2013) or magnetohydrodynamical mechanisms such as Blandford-Znajek process (Blandford & Znajek, 1977; McKinney & Gammie, 2004; Hawley & Krolik, 2006; Kawanaka et al., 2013).

Observations show that the prompt emission of GRBs is highly variable, and its typical timescale is as short as $10^{-3}-10^{-1}~{}{\rm s}$ (e.g. MacLachlan et al., 2013). Although the origin of such violent behavior is still uncertain, it is often considered to be driven by a certain instability of the NDAF (Kawanaka & Mineshige, 2007; Janiuk et al., 2007; Masada et al., 2007; Kawanaka & Kohri, 2012; Kawanaka et al., 2013; Kimura et al., 2015). When a disk becomes unstable, violent and sporadic mass accretion is expected, which would launch a highly variable jet from the vicinity of a central black hole. Most of the works about the disk instability in NDAFs (except for Masada et al., 2007) are based on the $\alpha$-viscosity model proposed by Shakura & Sunyaev (1973), where the ratio of the viscous stress to the isotropic pressure in an accretion disk is constant. This model has been widely accepted and used to explain the liberation of gravitational energy and behavior such as outbursts occurring in accretion disk systems such as dwarf novae, X-ray binaries, and active galactic nuclei (see Kato et al., 2008 for the detail). Balbus & Hawley (1991) pointed out the possibility that the viscosity in accretion disks comes from the magnetohydrodynamic (MHD) turbulence operated by the magnetorotational instability (MRI). Since their indication, the non-linear growth of MRI and the turbulent angular momentum transport have been investigated in many studies (Hawley & Balbus, 1992; Hawley et al., 1995; Sano et al., 2004; Fromang et al., 2007). In this context, the shear viscosity in an accretion disk is determined by the saturation level of the growth of MRI, which is assumed to be simply determined by the value of isotropic pressure in the $\alpha$-viscosity model.

According to the recent numerical studies, the saturated value of the viscosity depends on the magnetic Prandtl number, ${\rm Pm}=\nu/\eta$, where $\nu$ and $\eta$ are the kinetic viscosity and magnetic diffusivity, respectively (Fromang et al., 2007; Lesur & Longaretti, 2007; Simon & Hawley, 2009; Longaretti & Lesur, 2010; Simon et al., 2011). Especially, when the net magnetic flux is non-zero, the effective value of $\alpha$-viscosity is found to be proportional to the power of the magnetic Prandtl number, $\alpha\propto{\rm Pm}^{\beta}$, where $\beta$ is a positive number (Fromang et al., 2007; Simon & Hawley, 2009). Taking into account the effect of ${\rm Pm}$-dependence of $\alpha$-viscosity, the stability criteria for an accretion disk would be modified. For example, Takahashi & Masada (2011) constructed the steady state model for a geometrically thin, optically thick disk using ${\rm Pm}$-dependent $\alpha$-viscosity and found that the disk becomes gravitationally, thermally, or secularly unstable for a sufficiently large power-law index $\beta$. More recently, Potter & Balbus (2014) investigated the spectrum and time-dependent behavior of such an unstable disk and found that they are analogous to what are observed in flaring X-ray binary systems (see also Potter & Balbus, 2017).

In this paper we firstly discuss the stability of an NDAF with Pm-dependent $\alpha$-viscosity, and investigate its time-dependent behavior by one-dimensional simulations of an NDAF. In evaluating the magnetic Prandtl number in NDAFs, we use the results presented in Rossi et al. (2008), which discussed the viscosity and magnetic diffusivity in an NDAF taking into account the degeneracy of electrons. In Section 2 we construct an analytic model of an NDAF using variable $\alpha$-viscosity prescription and evaluate the condition that an NDAF becomes unstable. In Section 3 we present the time evolution of an NDAF with variable $\alpha$-viscosity and show that the mass accretion is highly time-variable in the unstable branch. Section 4 is devoted to the interpretation of our results and application to the prompt emission of GRBs, which show highly violent behavior. We summarize this work in Section 5.

2. Analysis

2.1. Fundamental Equations of an NDAF

In this section we construct the model of a neutrino-cooled accretion disk taking into account the Pm-dependence of $\alpha$, and investigate its stability. The vertically-integrated equations governing the structure of a neutrino-cooled accretion disk are the following (see Kawanaka & Mineshige, 2007; Kawanaka et al., 2013 and references therein):

1. Mass conservation equation:

[TABLE]

2. Kepler rotation (Newtonian gravity):

[TABLE]

3. Hydrostatic balance:

[TABLE]

4. Angular momentum conservation:

[TABLE]

5. Energy balance:

[TABLE]

6. $\alpha$ prescription:

[TABLE]

Here $\varv_{r}$, $\Sigma$, $\nu_{t}$, $R_{\rm in}$, $Q^{-}$, $p$, $\rho$ and $H$ are the radial velocity, the surface density, the kinematic viscosity due to the turbulence driven by MRI, the inner radius of the disk, the cooling rate per unit surface area, the disk pressure, the mass density of the disk, and the disk scale height, respectively. The surface density is defined as $\Sigma=2\rho H$.

The cooling rate $Q^{-}$ is in general described by the sum of two components:

[TABLE]

where $Q_{\rm adv}^{-}=T\Sigma\varv_{r}ds/dr$ is the advective cooling rate, where $s$ is the specific entropy and $T$ is the temperature in the disk midplane, and $Q_{\nu}^{-}$ is the cooling due to neutrino emission from the disk. The most effective emission process for electron-type neutrinos in an NDAF is the URCA process ($e^{-}/e^{+}$-capture onto a proton/neutron), whose emissivity (i.e. energy loss rate per unit volume) can be approximated as

[TABLE]

where $\rho_{10}=\rho/10^{10}~{}{\rm g}~{}{\rm cm}^{-3}$ and $T_{11}=T/10^{11}~{}{\rm K}$ (Di Matteo et al., 2002). As for the other types of neutrinos (i.e., mu- and tau-type neutrinos), electron-positron pair annihilation ($e^{-}+e^{+}\rightarrow\nu_{x}+\bar{\nu}_{x}$) is dominant, and its emissivity is approximately described as

[TABLE]

which is valid in the nondegenerate limit, and when electrons are strongly degenerated we can neglect these processes.

The neutrino cooling rate can be described by the two-stream approximation (Di Matteo et al., 2002; Kohri et al., 2005):

[TABLE]

where $\tau_{a,\nu_{i}}$ and $\tau_{s,\nu_{i}}$ represent absorptive and scattering optical depths, respectively, and $\tau_{\nu_{i}}=\tau_{a,\nu_{i}}+\tau_{s,\nu_{i}}$ is the sum of absorptive and scattering optical depth with respect to $\nu_{i}$ ($=\nu_{e},\bar{\nu}_{e}$). From the Kirchhoff’s law, the absorptive optical depth is given by

[TABLE]

and the scattering optical depth, which is dominated by nucleons, is given by

[TABLE]

both of which have the same temperature and density dependence (Di Matteo et al., 2002).

The pressure in a disk, $p$, is described as the sum of the contributions from radiation, baryonic gas, degenerate electrons, and neutrinos (if they are trapped):

[TABLE]

where $Y_{e}=n_{p}/(n_{p}+n_{n})$ is the electron fraction, $a\simeq 7.56\times 10^{-15}~{}{\rm erg}~{}{\rm cm}^{-3}~{}{\rm deg}^{-4}$ is the radiation constant, and $k_{B}$ is the Boltzmann constant.

When the mass accretion rate is the order of $\sim 0.01M_{\odot}~{}{\rm s}^{-1}$, the cooling rate $Q^{-}$ is dominated by the neutrino emission via URCA process, the optical depth of the disk with respect to neutrinos is thin, and the disk pressure is dominated by baryonic gas in the inner disk region (see Kawanaka & Mineshige, 2007 and references therein). Then we can approximate the disk pressure and the neutrino energy flux as

[TABLE]

and when the mass accretion rate is the order of $\sim 0.1M_{\odot}~{}{\rm s}^{-1}$, the disk pressure is still dominated by baryonic gas but the disk becomes optically-thick with respect to neutrinos and we can approximate the neutrino energy flux as

[TABLE]

where $\rho_{11}\equiv\rho/10^{11}~{}{\rm g}~{}{\rm cm}^{-3}$ and $T_{10}\equiv T/10^{10}~{}{\rm K}$ (Popham et al., 1999; Kawanaka et al., 2013).

2.2. Instability Criteria

Let us derive the stability criterion for this hyperaccretion flow. The condition for the secular instability can be described as

[TABLE]

where $T_{r\phi}=\rho\nu_{t}(d\Omega/d\ln r)$ is the shear viscous stress. The physical meaning of this condition is the following. When this inequality is fulfilled, a local increase of the surface density gives rise to the local decrease of the shear stress, which would suppress the local mass accretion rate. Then mass would accumulate locally and the surface density grows further. In this sense the initial perturbation given in an accretion flow would be amplified when Eq. (17) is satisfied.

Let us evaluate the condition for the secular instability in a disk with optically-thin neutrino cooling. From Eqs. (5), (6), (14) and (15) we obtain

[TABLE]

and from Eqs. (3) and (14) we obtain

[TABLE]

at a fixed radius. Then we can rewrite the instability condition (17) as

[TABLE]

at the fixed radius. On the other hand, the instability criterion for a disk with optically-thick neutrino cooling is obtained from Eqs. (5), (6), and (16) as

[TABLE]

at a fixed radius. Note that if $\alpha$ is constant and independent of density or temperature, these conditions are never fulfilled, which means that the accretion flow would be stable as long as its cooling is dominated by neutrino emission due to the URCA process.

Here we evaluate the magnetic Prandtl number, ${\rm Pm}=\nu/\eta$, in an NDAF to find the disk instability criterion. In Balbus & Henri (2008) they consider the microscopic kinematic viscosity $\nu$ is dominated by the Coulomb scattering or radiation, and the resistivity $\eta$ by the Coulomb scattering (Spitzer, 1962; Mihalas & Weibel Mihalas, 1984). However, these values are not applicable to NDAFs because they assume that electrons are non-relativistic and non-degenerated, which is generally not the case in NDAFs. Rossi et al. (2008) investigated the microscopic physical properties in NDAFs taking into account the relativistic degeneracy of electrons and found that the viscosity and resistivity are modified from the values adopted Balbus & Henri (2008). We fit the viscosity and resistivity presented in Fig. 3 of Rossi et al. (2008) by simple analytic functions and obtain the fitting formula for the magnetic Prandtl number as

[TABLE]

Then we can rewrite our assumption on the dependence of the viscosity parameter on the magnetic Prandtl number, $\alpha\propto{\rm Pm}^{\beta}$, as

[TABLE]

Using Eqs. (18), (19) and (20), we can evaluate the criterion for $\beta$ that the accretion disk cooled dominantly by neutrino emissions would suffer from the secular instability. In the case when the disk is optically-thin, the criterion is

[TABLE]

while in the case when the disk is optically-thick, the criterion reads as

[TABLE]

2.3. Disk Evolution Scenario

As shown in the previous section, a hyperaccretion flow can be secularly unstable when the index $\beta$ is within a certain range. When the instability sets in locally, the mass would be accumulated in the unstable zone and the surface density of a disk grows locally. As a result, the accretion flow would be gravitationally unstable and the gravitational torque would become effective. This causes the intense mass accretion onto a black hole from that zone. Just after that the local surface density becomes lower but as the mass accretion from the outer part continues that zone would be filled again and become secularly unstable again. In this way the intense mass accretion onto a black hole occurs repeatedly, which may be the origin of the short-term variability observed in the prompt emission of GRBs.

The suppression of MRI-driven turbulence and following evolution of an accretion flow have been discussed in the context of protoplanetary disks (Nakano, 1991; Jin, 1996; Gammie, 1996; Sano & Miyama, 1999; Sano et al., 2000; Armitage et al., 2001), in which the Ohmic dissipation of the electric current in a disk can suppress the growth of MRI. In the context of a hyperaccretion flow as a central engine of GRBs, Masada et al. (2007) discuss the suppression of MRI-driven turbulence due to neutrino diffusion in a neutrino-thick hyperaccretion flow whose mass accretion rate is as high as $\sim 1M_{\odot}~{}{\rm s}^{-1}$, followed by intermittent mass accretion driven by gravitational instability. However, in the current model it is not necessary for a hyperaccretion flow to be optically-thick with respect to neutrinos, as shown above. As long as an accretion flow is dominantly cooled via neutrino emissions the secular instability and subsequent intermittent accretion can appear for a certain range of $\beta$. In the next section, we present the results of one-dimensional simulations of a hyperaccretion flow according to this evolution scenario.

3. 1D Simulation of Disk Evolution

3.1. Basic Equation and Simulation Setup

We perform one-dimensional numerical simulations as a demonstration of the dynamical evolution of the NDAF with Pm-dependent $\alpha$-viscosity. The basic equations are the mass and angular momentum conservation in the disk,

[TABLE]

and the energy conservation,

[TABLE]

where $\chi$ is the ratio of gas pressure to total pressure, $Q^{+}\equiv\nu_{t}\Sigma({\rm d}\Omega/{\rm d}\ln r)^{2}$ is the viscous heating rate, and $Q^{-}$ is the neutrino cooling rate described by eq.(10). The radial velocity of the accretion flow is given by

[TABLE]

The advection terms in equation (27) are neglected in our actual simulation for the simplicity since it can be expected that they have little impact on the energy balance in the relatively small mass accretion regime ($\dot{M}=\mathcal{O}(10^{-2})\ M_{\odot}/{\rm sec}$) of our interest.

We take into account the Pm-dependence of $\alpha$ in the similar way as Potter & Balbus (2014):

[TABLE]

While the ceiling value of $\alpha$ is fixed to $\alpha_{\rm max}=0.1$, the floor value $\alpha_{\rm min}$ is remained as a control parameter of our simulation. This functional form bridges the small and large Pm regimes smoothly, and behaves approximately as being proportional to ${\rm Pm}^{\beta}$ when ${\rm Pm}\lesssim{\rm Pm}_{c}$. Note that ${\rm Pm}$ is given, as a simple function of $\rho$ and $T$, by eq.(22) in this study. The critical value ${\rm Pm}_{c}$ is varied depending on the simulation model.

As stated in § 2.2, when a part of an accretion disk becomes secularly unstable, mass would accumulate in that part and the surface density locally grows with time, resulting in the formation of the gravitationally unstable region at some evolutionary stage.

In our simulations, the condition for the onset of the gravitational instability is given, with Toomre’s $\mathcal{Q}$ parameter, by

[TABLE]

where $c_{s}=(p/\rho)^{1/2}$ is the speed of sound and $\mathcal{Q}_{\rm crit}$ is the critical value for the instability which is evaluated to be in the range $1\lesssim\mathcal{Q}_{\rm crit}\lesssim 5$ from the numerical simulation of the accretion disk. Once an accretion flow becomes gravitationally unstable, non-axisymmetric structure and hence the gravitational torque arises, leading to the outward angular momentum transport in the disk. To model this process, we simultaneously solve the evolution equation of $\alpha_{\rm grav}$ ($\alpha$-parameter due to the gravitational instability), similar to the one proposed by Zhu et al. (2010),

[TABLE]

where $\alpha_{0}=0.1$ is the reference value, $\sigma$ is a control parameter for the evolution of $\alpha_{\rm grav}$ and chosen so as to match the typical growth rate of the gravitational instability, i.e., $\sigma=\Omega\sqrt{(\mathcal{Q}_{\rm crit}/\mathcal{Q})^{2}-1}$. The last term of the RHS represents the turbulent diffusion of $\alpha_{\rm grav}$ and the coefficient $\nu_{g}$ is assumed to be the product of the typical size and velocity of the turbulence, i.e., $\nu_{g}=c_{s}^{2}/\Omega$ based on the mixing-length theory.

Computationally, we solve equations (26), (27), (28), and (31) with explicit finite difference method, and then include the self-gravity as an additional sweep over the grid using $\nu_{\rm grav}\equiv\alpha_{\rm grav}c_{s}^{2}/\Omega$ and the surface density $\Sigma$, and as an additional heating term in the energy equation when the condition (30) is satisfied. The computational grid of $256$ radial mesh-points is uniform in a scaled radial variable $x\propto r^{1/2}$.

For the numerical setup, we adopt the mass of the central black hole $M_{\rm BH}=3M_{\odot}$, an inner disk radius $r_{\rm in}=3r_{g}$, and an outer disk radius $r_{\rm out}=50r_{g}$, where $r_{g}=GM_{\rm BH}/c^{2}=4.5\times 10^{5}$ cm is the Schwarzschild radius. A zero-torque boundary condition is imposed at $r_{\rm in}$, while at $r_{\rm out}$ we impose two types of boundary conditions: (1) the constant mass accretion rate, i.e., $\dot{M}_{\rm out}={\rm const.}$ (§ 3.2), and (2) the time-decreasing mass accretion rate, i.e., $\dot{M}_{\rm out}\propto t^{-5/3}$ (§ 3.3), which is expected in the fallback accretion just after the failed supernova explosion or the merger of a neutron star binary. For the initial profiles of $\Sigma$ and $T$, the analytic solutions of equations (1)–(6), and (10) with constant mass accretion rate $\dot{M}$ and constant viscous parameter $\alpha$ are adopted (Kawanaka et al., 2013). After the relaxation to the equilibrium state with the constant $\alpha$, we take into account the Pm-dependence of $\alpha$ and restart the time-dependent calculation.

3.2. Models with Constant Mass Injection

Shown in Figure 1 is the temporal evolution of the mass accretion rate in unit of $M_{\odot}{\rm/sec}$ for the models with $\beta=0.5$ (blue), $\beta=0.6$ (green) and $\beta=0.75$ (red). The mass accretion rate is measured at the innermost radius. The mass injection rate from the outer boundary is constant and common among these models at $\dot{M}_{\rm out}\simeq 0.025M_{\odot}{\rm/sec}$. The parameters $\alpha_{\rm min}=10^{-4}$, $\rm{Pm}_{\rm c}=10$, and $\mathcal{Q}_{\rm crit}=5.0$ are also adopted commonly to these models.

In the case with lower value of $\beta$ (see models with $\beta=0.5$ and $0.6$), the disk gradually evolves to the equilibrium state with the Pm-dependent $\alpha$ and then converges into the time-steady accretion stage compatible with the mass injection from the outer boundary. In contrast to them, the model with $\beta=0.75$ shows a time-dependent episodic mass accretion while the average value of $\dot{M}$ is almost same with other models. In this model, as predicted from the linear theory, the inner region becomes secularly unstable because of the Pm-dependent $\alpha$ and as a result the gravitational instability is developed and makes the mass accretion sporadic.

Figure 2 shows the radial profiles of physical variables during the early evolution of the disk with $\beta=0.75$ from the initial state to the quasi-steady state (around $t\simeq 4$ sec just before the initial spike), where the linear growth of the physical variables are terminated and the mass accretion becomes almost steady when taking long-term average of it. Panels (a)–(f) are for (a) $\Sigma$, (b) $T$, (c) $\alpha$, (d) $\tau$ (neutrino depth), and (e) Toomre’s $\mathcal{Q}$. In each panel, the blue line corresponds to the initial equilibrium state with constant $\alpha$, the red line is the quasi-steady state and the gray lines are in-between them, where the lines with darker tones depict the profiles at later times.

In the initial state the disk is geometrically thin and neutrino-thin, and the Toomre’s $\mathcal{Q}$ is $\mathcal{O}(100)$ over the whole disk. When applying the Pm-dependent $\alpha$-viscosity to it, $\alpha$ becomes lower in the inner part of the disk and then the secular instability starts to grow there. Since the secular instability leads to the mass accumulation, the surface density and thus the neutrino depth also start to grow. As a result, Toomre’s $\mathcal{Q}$ decreases with time and finally reaches $\mathcal{Q}_{\rm crit}$ which causes the gravitational instability, providing a drastic mass accretion. Once the gravitational instability grows, the inner part of the disk repeats the mass accumulation stage and the gravitationally unstable stage, alternately. This would be a reason why the mass accretion falls into the episodic mode in the secularly unstable disk with $\beta=0.75$.

Figure 3 depicts the time profiles of the $\dot{M}$ (black) and the maximum value of $\alpha_{\rm grav}$ (red) focusing on $6~{}{\rm sec}\lesssim t\lesssim 20~{}{\rm sec}$, showing that the mass accretion rate and the torque driven by gravitational instability are synchronized. We can see that the cycle of mass accretion coincides well with the modulation cycle of the maximum value of $\alpha_{\rm grav}$ with relatively long period, $t_{\rm cycle}\sim 3~{}{\rm sec}$.

This cycle period can be estimated from the redistribution time of the baryonic matter through the secularly unstable region with the suppressed value of $\alpha$. Since $\alpha$ is reduced to $10^{-4}$ in the inner secularly unstable region, the viscous redistribution time $t_{\rm vis}$ of the matter through there is evaluated as

[TABLE]

equivalent to the modulation period, where $r_{\rm SI}\sim\mathcal{O}(1)r_{g}$ is the radius of the secularly unstable region and $\nu_{\rm SI}$ is the magnitude of the viscosity there, which is evaluated from the simulation result as $\nu_{\rm SI}\simeq 10^{12}\ {\rm cm^{2}/sec}$. This implies that the global change of the disk structure is regulated by the floor value of Pm-dependent $\alpha$, which is one of the control parameters in our simulation model.

In addition to the longer period mass accretion cycle, we can find the shorter time-variability with a duration $\lesssim 1\ {\rm sec}$ in $\dot{M}$. Since the shorter wavelength mode of the secular instability has a faster growth rate, we can expect that such an additional component should be generated due to the development of smaller-scale structures by the localized secular instability. Another possibility is that this component comes from the thermal instability occurring in the disk, whose timescale is generally much shorter than the viscous timescale (see Sec.4).

In Figure 3, the temporal evolution of the neutrino luminosity integrated over the whole disk, $L_{\nu}$, is also shown by a blue line. As seen from this, $L_{\nu}$ also changes following the modulation of $\alpha_{\rm grav}$, suggesting again the change of the global disk structure with this period.

3.3. Models with Power-law Decay of Mass Injection

The mass accretion rate from the outer boundary of the disk is expected to decrease with time in the realistic situation of the collapse of a massive star, which is prerequisite for the NDAF formation. As the demonstration of such situations, under the condition of decaying mass injection, we simulate the models with variety of parameter sets which give the secularly and gravitationally unstable disks. In these models, the outer boundary condition is time-dependent and is controlled as

[TABLE]

where $\Sigma_{{\rm out},0}$ is the initial surface density at the outer boundary, and $t_{\rm ds}$ is the start time of the decay. We choose here $t_{\rm ds}=5\ {\rm sec}$, which corresponds to the timing of the initial spike of the mass accretion.

Shown in Figure 4 is $\dot{M}$ (red) and $L_{\nu}$ (blue) for the model with the same parameter sets as that in Figure 3 but with decaying mass injection. Due to the decrease of $\Sigma_{\rm out}$, the duration of both $\dot{M}$ and $L_{\nu}$ becomes shorter and converges into $\sim 10$ sec. Such a temporal behavior is reminiscent of the light curves of GRBs. Depending on the parameter sets, the temporal evolution of $\dot{M}$ exhibits various behaviors.

The temporal evolution of $\dot{M}$ for the models with different parameter sets are demonstrated in Figure 5. Panels are for the models with (a) $\beta=0.8,\ \alpha_{\rm min}=10^{-4},\ \mathcal{Q}_{\rm crit}=5.0,\ {\rm Pm}_{\rm c}=10.0$, (b) $\beta=0.8,\ \alpha_{\rm min}=10^{-5},\ \mathcal{Q}_{\rm crit}=3.0,\ {\rm Pm}_{\rm c}=10.0$, (c) $\beta=0.9,\ \alpha_{\rm min}=2\times 10^{-4},\ \mathcal{Q}_{\rm crit}=3.0,\ {\rm Pm}_{\rm c}=15.0$, and (d) $\beta=0.7,\ \alpha_{\rm min}=2\times 10^{-4},\ \mathcal{Q}_{\rm crit}=3.0,\ {\rm Pm}_{\rm c}=10.0$. The inset of each panels focus on the shorter-time variability in the first spike of $\dot{M}$.

We can see that with the smaller value of $\alpha_{\rm min}$, the modulation cycle becomes longer (see, panels (a)–(c)). While the time-decaying behavior is common among them, in the case of model (d) with $\beta=0.7$ which is located near the critical value, the cyclic behavior almost disappears, showing one longer humped structure. Of course, we understand that the time-scale of the variability should be studied by more realistic MHD models of the NDAF because the various important effects lacks in our 1D simulation, such as advection effects, pressure gradient effect, and non-linear interaction between them. Even though our model is far from the realistic one, such a variety of temporal behavior of the mass accretion rate is reminiscent of the diversity of the light curve of GRBs.

4. Discussion

4.1. Disk Evolution on $\dot{M}$–$\Sigma$ plane

The episodic mass accretion process seen in our simulations naturally provokes a question “how an accretion disk changes its state on $\dot{M}-\Sigma$ or $T-\Sigma$ plane in this model ?”. To answer it, we discuss the evolution path of a disk on these planes here. In the following, we cast a spotlight on the model demonstrated in Figure 5(b) as a typical one.

Figure 6(a) shows the time variation of $\dot{M}$ focusing on the period around the secondary burst ($t\sim 10$ sec) for that model. The thick green line, that denotes the low resolution time-series, is superimposed on the gray curve that denotes highly time-resolved data. Shown in the panels (b) and (c) in Figure 6 are the disk evolution on $\dot{M}$–$\Sigma$ and $T$–$\Sigma$ planes at the disk inner edge, $r=3r_{g}$. The reference points $1$–$12$ in the panels (b) and (c) correspond to the states at the reference times $1$–$12$ in the panel (a), respectively.

In the earlier stage when a disk is optically thin with respect to neutrinos, both the $\dot{M}$ and $T$ are lower because of the low $\alpha$-value suppressed at the inner region due to its $Pm$-dependence (states $1$–$2$). Since the $\dot{M}$ is higher in the outer region, the baryonic matter is accumulated in the inner part during this stage and thus the $\Sigma$ gradually increases with time (states $3\rightarrow 5$). When $\Sigma$ gets larger, $T$ becomes higher because viscous heating becomes more effective, resulting in the increase of $\alpha$ and thus $\dot{M}$. The growths of the $\Sigma$ and $T$ are terminated when the Toomre’s $\mathcal{Q}$ reaches the critical value at the inner region. Then, the disk enters the bursting mass accretion stage with high turbulent viscosity ($\because\alpha_{\rm grav}>\alpha$) (state $6$), and finally returns to the thinner and cooler state (states $6\rightarrow 12$) similar to the starting point.

Although one typical mode is focused here, the other models also show the similar behavior qualitatively on the $\dot{M}$–$\Sigma$ plane. To summarize, our results thus show that the evolution of an NDAF in our model has the limit-cycle property which is reminiscent of that in the other astrophysical disk models, such as dwarf novae (e.g. Honma et al., 1991).

4.2. Application to Astrophysical Objects

The origin of short-term variability in the prompt emission of GRBs is still unclear. One of the most promising mechanisms is the accretion disk instability, which has been discussed to interpret the flaring behavior of X-ray binaries and dwarf novae, etc. Several models of the disk instability in a hyperaccretion flow have been proposed (Janiuk et al., 2007; Masada et al., 2007; Kawanaka & Kohri, 2012; Kawanaka et al., 2013; Kimura et al., 2015). However, most of these models rely on the prescription with a constant $\alpha$-viscosity, and since the mass accretion process is governed by the angular momentum transportation driven by the growth of MRI inside the disk, it is necessary to follow the fundamental MHD process in the neutrino-emitting gas in order to investigate the variability occurring in the hyperaccretion flow. Only Masada et al. (2007) discussed the possibility that the angular momentum transport is suppressed because MRI in a neutrino-thick accretion flow would be suppressed, which would cause the gravitational instability in the inner disk. Such a mechanism works only if the accretion rate is so high that the accretion flow becomes dense and the optical depth with respect to neutrinos is larger than unity. We have applied the prescription in which the $\alpha$-viscosity is not constant but depends on the magnetic Prandtl number, ${\rm Pm}=\nu/\eta$, to a hyperaccretion flow for the first time. The $\alpha$ parameter represents the saturation level of MRI in the disk, which would be impacted by magnetic reconnection rate of the MRI-amplified magnetic field. The higher viscosity suppresses efficiently the small scale turbulent structure, sustaining the larger size of the turbulent vortices and thus the larger size of the frozen-in magnetic structure. Since magnetic reconnection rate is expected to be smaller with the larger size of the turbulent magnetic structure and the smaller magnetic resistivity, the saturation level of the MRI-driven turbulence would become higher with the increase of the Pm. We make use of the results of Rossi et al. (2008) to express the magnetic Prandtl number in a dense gas, where electrons are relativistically degenerate, as an analytic function of temperature and density. By deriving the thermal equilibrium solution of a hyperaccretion flow which is cooled efficiently via neutrino emission, assuming $\alpha$ depends on the magnetic Prandtl number as $\alpha\propto{\rm Pm}^{\beta}$, we show that when the exponent $\beta$ takes a value in a certain range a neutrino-cooled accretion flow becomes secularly unstable. To investigate what will occur in a hyperaccretion flow when this condition is fulfilled we perform one-dimensional simulations of a hyperaccretion flow with Pm-dependent $\alpha$-viscosity and show that if the index $\beta$ is within this range the mass would be accumulated in the inner part of the accretion flow because the angular momentum transportation due to the turbulence driven by MRI would be suppressed as the surface density grows. This mass accumulation makes an accretion flow gravitationally unstable at some stage and intermittent mass accretion onto a black hole due to the gravitational torque would occur. This may cause the ejection of inhomogeneous jets, which can be the origin of the short-term variability of GRBs.

As can be seen in Fig. 1 and 3, the variability of mass accretion can be decomposed into long-term cycles whose timescale is $\sim\mathcal{O}(1)~{}{\rm sec}$ and short-term oscillations whose timescale is $\lesssim 0.1-1~{}{\rm sec}$. Actually it has been suggested that a GRB lightcurve contains two components whose timescales are different from each other (e.g. Vetere et al., 2006), which may be produced by the mechanism shown in our simulations. The long-term oscillation corresponds to the viscous redistribution of disk material which is regulated by the floor value of Pm-dependent $\alpha$, while it is not obvious what makes the latter oscillation. One of the possibilities is that it comes from the shorter wavelength modes of the secular instability, which have a faster growing rate. Another possibility is that the thermal instability, whose timescale is generally much shorter than the viscous timescale. Let us discuss the conditions for thermal instability to occur in a hyperaccretion flow. We can check if the disk is thermally stable or not by comparing the temperature dependence of $Q^{+}$ and $Q^{-}$ in the thermal equilibrium, keeping the surface density $\Sigma$ constant (Pringle, 1976). In the optically-thin regime, from Eqs.(3), (5), (6), and (23), we have the temperature dependence of the heating rate as $Q^{+}\propto T^{5\beta/2+1}$, and from Eq. (15) we have $Q^{-}\propto T^{6}$. Therefore, an optically-thin disk becomes thermally unstable when $\beta>2$, and we can see that this condition contradicts with Eq.(24). In other words, an optically-thin disk cannot become secularly unstable and thermally unstable at the same time. In the optically-thick disk, however, the situation changes completely. Using Eq.(16), we have $Q^{-}\propto T^{2}$ and the condition that an optically-thick disk becomes thermally unstable can be described as $\beta>2/5$. This is identical with the condition that an optically-thick disk becomes secularly unstable, Eq. (25). Therefore, in the optically-thick regime, a disk that is secularly unstable is always thermally unstable. As can be seen in Fig.2 the innermost region of the disk becomes optically-thick as the secular instability develops, and so it should be thermally unstable at the same time. The thermal timescale in an accretion disk is described as $t_{\rm th}\sim 1/(\alpha\Omega)$, while the viscous timescale is described as $t_{\rm vis}\sim 1/(\alpha\Omega)(r/H)^{2}$. Since the disk is geometrically thin (i.e. $H\ll r$), the thermal timescale is much shorter than the viscous timescale. This means that when a disk is in the optically-thick regime, the thermal instability develops faster than the secular instability. As shown in Fig.3 the time profiles of mass accretion rate have the shorter-time variability with duration $\lesssim 1~{}{\rm sec}$. Since the disk simulated here is in the optically-thick regime at the evolved state, it is not only secularly unstable but also thermally unstable. Thus, the short-time variability may originate essentially from the thermal instability whose timescale is much shorter than the viscous timescale (e.g., Janiuk et al., 2007). It is beyond the scope of this paper to specify the origin of this short-term variability but is one of the main targets of our future work with “more realistic” multi-dimensional time-dependent model.

The viscosity in a hyperaccretion flow is important not only in the context of the central engine of GRBs but also in the context of the electromagnetic counterpart of gravitational waves from binary neutron star mergers or black hole - neutron star mergers. From the remnant accretion disks of compact binary mergers, one can expect the disk outflow driven by viscous heating (Metzger et al., 2008, 2009; Fernández & Metzger, 2013; Metzger & Fernández, 2014; Martin et al., 2015; Siegel & Metzger, 2017; Lippuner et al., 2017; Siegel & Metzger, 2018; Fujibayashi et al., 2017, 2018; Fernández et al., 2018; Radice et al., 2018). For example, Fujibayashi et al. (2018) investigated the mass ejection from a torus surrounding a massive neutron star, which is formed after the merger of binary neutron stars by general relativistic neutrino-radiation hydrodynamics simulations. They discussed the importance of mass ejection due to the viscous heating in the torus, using $\alpha$-prescription to model the shear viscosity. However, they assume that $\alpha$ is a constant ($\sim 0.01-0.04$) in a whole torus. As we have shown above, due to the degeneracy of electrons the $\alpha$-parameter may no longer be constant in the innermost region of a hyperaccretion flow. In fact, we can see in Fig.2 (c) that the $\alpha$-parameter largely depends on radius. Therefore, the behavior of mass ejection due to the viscosity can change if taking into account this effect, which is beyond the scope of this paper.

5. Summary

We investigate the stability of a hyperaccretion flow, which is considered to be the central engine of GRBs. Especially we take into account the variation of the viscosity parameter $\alpha$ depending on the magnetic Prandtl number, which has recently been implied by several MHD simulations, in the context of a hyperaccretion flow for the first time. We have shown that a hyperaccretion flow becomes secularly unstable when $\alpha$ is proportional to the magnetic Prandtl number with an exponent $\beta$ within a certain range, which is consistent with the value inferred from recent MHD simulations. Once a hyperaccretion flow becomes secularly unstable, the surface density of the inner region increases due to the mass accumulation. As a result, the inner region of the accretion flow becomes gravitationally unstable, and intense mass accretion would occur induced by non-axisymmetric pattern in the accretion flow. Since such an accretion process can occur repeatedly as long as matter is sufficiently supplied from the outer region, we can expect sporadic mass accretion in this model, which may account for the short-term variability observed in the prompt emission of GRBs. We also perform one-dimensional simulations of a hyperaccretion flow assuming the Pm-dependence of the $\alpha$-parameter, and we can reproduce the intermittent behavior of an accretion flow which is secularly unstable.

This work was supported by JSPS KAKENHI Grants No. 18K03700, No. 18H04444, and No. 18H01212. NK acknowledges the support from the Hakubi project at Kyoto University. Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.