Instabilities in the time-dependent neutrino disc in Gamma-Ray
Bursts
A. Janiuk1,2, Y. Yuan3, R. Perna4 & T. Di Matteo5
ABSTRACT
We investigate the properties and evolution of accretion tori formed after the coalescence of
two compact objects. At these extreme densities and temperatures, the accreting torus is cooled
mainly by neutrino emission produced primarily by electron and positron capture on nucleons
(β reactions). We solve for the disc structure and its time evolution by introducing a detailed
treatment of the equation of state which includes photodisintegration of helium, the condition of
β-equilibrium, and neutrino opacities. We self-consistently calculate the chemical equilibrium in
the gas consisting of helium, free protons, neutrons and electron-positron pairs and compute the
chemical potentials of the species, as well as the electron fraction throughout the disc. We find
that, for sufficiently large accretion rates (Ṁ & 10M⊙/s), the inner regions of the disk become
opaque and develop a viscous and thermal instability. The identification of this instability might
be relevant for GRB observations.
Subject headings: accretion, accretion discs – black hole physics – gamma rays: bursts – neutrinos
1. Introduction
Gamma-Ray Bursts are commonly thought to
be produced in relativistic ejecta that dissipate en-
ergy by internal shocks however alternative ideas
based on the Poynting flux dominated jets are also
being proposed (for a review see e.g. Piran 2005;
Meszaros 2006; Zhang 2007).
The enormous power released during the gamma-
ray burst explosion indicates that a relativistic
phenomenon must be involved in creating GRBs
(Narayan et al. 1992). The merger of two neutron
stars or of a neutron star and a black hole (e.g. NS-
NS or NS-BH) has been invoked e.g. by Paczyński
(1986); Eichler et al. (1989); Paczyński (1991);
1Copernicus Astronomical Center, Bartycka 18, 00-716
Warsaw, Poland
2Present address: University of Nevada, Las Vegas, 4505
Maryland Pkwy, NV89154, USA
3 Center for Astrophysics, University of Science and
Technology of China, Hefei, Anhui 230026, P.R. China
4JILA and Department of Astrophysical and Planetary
Sciences, University of Colorado, 440 UCB, Boulder, CO
80309, USA
5Physics Department, Carnegie Mellon University, 5000
Forbes Avenue, Pittsburgh, PA 15232, USA
Narayan, Paczyński & Piran (1992), as well as the
collapse of a massive star, the so called “collap-
sar” scenario (e.g., Woosley 1993 and Paczyński
1998). In both cases a dense, hot accretion disk
is likely to form around a newly born black hole
(Witt et al. 1994). In the collapsar scenario, the
collapsing envelope of the star accretes onto the
newly formed black hole, while a transient debris
disk is formed when the NS-NS or NS-BH binary
merges (see e.g. Ruffert et al. 1997 for numerical
simulations).
The durations of GRBs, which range from mil-
liseconds to over a thousand of seconds, are dis-
tributed in two distinct peaks defining two main
GRB classes: short (≤ 2 sec) and long (& 2 sec)
bursts (Kouvelietou et al. 1993). For some long
bursts, signatures of an accompanying supernova
explosion have been detected in the afterglow spec-
tra (Stanek et al. 2003; Hjorth et al. 2003), which
strongly favors the “collapsar” interpretation for
their origin. Furthermore, the GRB positions in-
ferred from the afterglow observations are consis-
tent with the GRBs being associated with the star
forming regions in their host galaxies. For short
bursts, several pieces of evidence from the analy-
http://arxiv.org/abs/0704.1325v1
sis of the Swift satellite and follow-up observations
(Gehrels et al. 2005; Fox et al. 2005; Villasenor et
al. 2005) argue in favor of a binary merger model
(Hjorth et al. 2005; Berger et al. 2005).
In the merger scenario, the duration of the
burst is comparable to the viscous timescale of the
accretion disc whereas in the collapsar scenario,
the external reservoir of stellar matter can feed
the accretion torus for a much longer timescale.
In general, accretion discs powering GRBs are ex-
pected to have typical densities of the order of
1010−12 g cm−3 and temperatures of 1011 K within
10−20 Schwarzschild radii (RS = 2GM/c2). Thus,
accretion proceeds with rates of a fraction to sev-
eral solar masses per second. In this “hyper-
accreting” regime, photons become trapped and
are not efficient at cooling the disc. Neutrinos,
however, are produced by weak interactions in the
dense and hot plasma, releasing the gravitational
energy of the accretion flow. These discs go under
the name of “neutrino-dominated accretion flows”,
or NDAFs.
Over the last several years a number of stud-
ies have investigated the structure of these discs
(Popham, Woosley & Fryer 1999; Narayan, Piran
& Kumar 2001; Kohri & Mineshige 2002; Di Mat-
teo, Perna & Narayan 2002; Surman & McLaugh-
lin 2004; Kohri, Narayan & Piran 2005; Chen &
Beloborodov 2006; Gu, Liu & Lu 2006; Liu, Gu
& Lu 2007). These models have employed the
customary approximation of one-dimensional hy-
drodynamics (Shakura & Sunyaev 1973), where
the effects of MHD viscous stresses are described
by the dimensionless parameter α, but have been
limited to the steady-state approximation of con-
stant Ṁ . Such an assumption is a good approx-
imation when considering the collapsar scenario,
where the burst duration is much longer than the
viscous time, due to the continuous replenishing of
the disc by the collapsing star envelope. However,
even in this scenario, recent observations suggest
that that engines are “long-lived” (past the torus
feeding phase), requiring a time-dependent com-
putation. In the merger scenario, on the other
hand, a time-dependent calculation is necessary
even for modeling the prompt phase of the burst,
since the duration is set by the viscous timescale
of the disc.
Recently, a fully time-dependent calculation of
the structure of such accretion discs has been pre-
sented by Janiuk et al. (2004). The model was
suitable for the torus being a results of either the
gravitational collapse of a massive stellar core or
the compact binary merger. However the struc-
ture and evolution of the disk (like in many of the
early calculations) was calculated under a num-
ber of simplifying assumptions for the composition
and the equation of state of the accreting matter.
In this paper we improve upon our earlier results
in several ways. Besides the requirement of time-
dependent calculations, the high density and tem-
perature regime in which the accreting gas lies, im-
plies that both multi-dimensional numerical and
semi-analytic calculations for such flows need to
include the detailed microphysics. This includes
photodisintegration of nuclei, the establishment
of statistical equilibrium, neutronization, and the
effects of neutrino opacities in the inner regions.
Here, we introduce a detailed treatment of the
equation of state, and calculate self-consistently
the chemical equilibrium in the gas that consists
of helium, free protons, neutrons and electron-
positron pairs. We compute the chemical poten-
tials of the species, as well as the electron fraction
throughout the disk using the assumption of the
equilibrium between the beta processes. Our EOS
equations include, self-consistently, the contribu-
tion of the neutrino trapping to the beta equilib-
rium. Another important addition compared to
our previous work (Janiuk et al. 2004) is the inclu-
sion of photodisintegration of helium. The pres-
ence of this term can affect the energy balance in
the inner, opaque (to neutrinos as well as pho-
tons) region of the flow and, as it will be shown, it
eventually produces a thermal and viscous insta-
bility in those regions. This is especially relevant
since the GRB phenomenology requires a variable
energy output.
Other time-dependent disc studies of binary
mergers or collapsars have been performed in
2D using hydrodynamical simulations (e.g. Mac-
Fadyen & Woosley 1999; Ruffert & Janka 1999;
Lee et al. 2002; Rosswog et al. 2004; Lee,
Ramirez-Ruiz & Page 2005) and, most recently,
in 3-D simulations (Setiawan et al. 2005). Also,
MHD simulations of the GRB central engine have
been performed, showing that the magnetic field
possibly plays an important role in the genera-
tion of a GRB jet (Proga et al. 2003; Fujimoto
et al. 2006). The advantage of our calculations
is that, whilst including all the relevant physics
to calculate the equation of state, the structure
and stability of the accretion disc, we are able to
study a much larger range of parameter space and
allow our calculations to evolve beyond what can
be reached in higher dimensional calculation and
comparable to at least the short-burst durations.
The paper is organized as follows. In Sec-
tion 2 we describe the basic assumptions of the
model and the method used in the initial sta-
tionary and subsequent time-dependent numerical
simulations. In Section 3 we discuss the structure
of the hyperaccreting disc for various values of the
initial accretion rate, and study the time evolu-
tion of its density and temperature, as well as the
resulting neutrino lightcurve. We also discuss the
physical origin of the instabilities in the disc, and
we compare our model and results with the recent
2D and 3D simulations. We summarize our results
in Section 4.
2. Neutrino-cooled accretion disks
In this Section, we describe how we improve
upon our previous time-dependent calculation (Ja-
niuk et al. 2004) by computing self-consistently
the equation of state of the extremely dense mat-
ter by solving the balance of the β reaction rates.
This allows us to determine the chemical poten-
tials of electrons, protons and neutrons, as well as
the electron fraction, in the initial disc configura-
tion and throughout its evolution.
2.1. Initial disc configuration: 1-D Hydro-
dynamics
We start by considering a steady-state model
of an accretion disc around a Schwarzschild black
hole - formed as a remnant structure either after
a compact binary merger, or in a collapsar after
the birth of a black hole (for a recent calculation
in Kerr spacetime see Chen & Beloborodov 2006).
Throughout our calculations we use the vertically
integrated equations and hence derive a vertically
averaged disc structure. We write the surface den-
sity of the disk as Σ = Hρ, where ρ is the density
and where the disk half thickness (or disk height)
is given by H = cs/ΩK . Here the sound speed
is defined by cs =
P/ρ and ΩK =
GM/r3
is the Keplerian angular velocity with P the to-
tal pressure. We note that, at very high accretion
rates, the disc becomes moderately geometrically
thick (H ∼ 0.5r) in regions where neutrino cooling
becomes inefficient and advection dominates. Our
’slim disk’ approximation neglects terms ∼ (H/r)2
and assumes that the fluid is in Keplerian rota-
tion. For the disc viscous stress we use the stan-
dard α viscosity prescription of Shakura & Sun-
yaev (1973) where the stress tensor is proportional
to the pressure:
τrϕ = −αP. (1)
We adopt a value of α = 0.1
We set the inner radius of the disc at 3 RS, while
the outer radius is at 50 RS. The initial mass of
such a disc is about 0.35M⊙ for an accretion rate
Ṁ = 1 M⊙/s. Throughout the calculations we
adopt a black hole mass of M = 3M⊙.
2.2. The equation of state
We assume that the torus consists of helium,
electron-positron pairs, free neutrons and protons.
The total pressure is contributed by all particle
species in the disc, and the fraction of each species
is determined by self-consistently solving the bal-
ance of the beta reaction rates. In the equation
of state we take into account the pressure due to
the free nuclei and pairs, helium, radiation and the
trapped neutrinos:
P = Pnucl + PHe + Prad + Pν . (2)
The component Pnucl includes free neutrons, pro-
tons, and the electron-positron pair gas in beta
equilibrium:
Pnucl = Pe− + Pe+ + Pn + Pp (3)
(~c)3
F3/2(ηi, βi) +
βiF5/2(ηi, βi)
where Fk are the Fermi-Dirac integrals of the or-
der k, and ηe, ηp and ηn are the reduced chemi-
cal potentials of electrons, protons and neutrons
in units of kT , respectively (where ηi = µi/kT ,
also known as the degeneracy parameter, where µi
the standard chemical potential) calculated from
the chemical equilibrium condition (§ 2.3). The
reduced chemical potential of positrons is ηe+ =
−ηe − 2/βe and the relativity parameters of the
species i are defined as βi = kT/mic
Under the physical conditions in the torus,
helium is generally non-relativistic and non-
degenerate; therefore, its pressure is given by:
PHe = nHekT, (5)
where nHe is the number density of helium. This
is defined as:
nHe =
nb(1 −Xnuc) , (6)
and the fraction of free nucleons is given by
Xnuc = 295.5ρ
11 exp(−0.8209/T11), (7)
with T11 the temperature in unit of 10
11 K (e.g.
Qian & Woosley 1996; Popham et al. 1999).
The radiation pressure is given by:
Prad =
(kT )4
(~c)3
. (8)
When neutrinos become trapped in the disc, the
neutrino pressure is non-zero. Following the treat-
ment of photon transport under the two-stream
approximation (Popham & Narayan 1995; Di Mat-
teo et al. 2002), we have
(kT )4
3(~c)3
i=e,µ,τ
(τa,νi + τs) +
(τa,νi + τs) +
3τa,νi
(kT )4
3(~c)3
b, (9)
where τs is the scattering optical depth due to
the neutrino scattering on free neutrons and pro-
tons and τa,νe and τa,νµ are the absorptive opti-
cal depths for electron and muon neutrinos, re-
spectively (see§2.3). The contribution from tau
neutrinos is the same as that from muon neutri-
nos. These optical depths and neutrino absorption
processes (which are the reverse of the emission
processes) are discussed in more detail in the Ap-
pendix.
In the disc we have to consider both the neu-
trino transparent and opaque regions, as well as
the transition between the two. In the trans-
parent case, the neutrinos are not thermalized
and the chemical potential of neutrinos is negli-
gible. On the other hand, when neutrinos are
totally trapped, the chemical equilibrium condi-
tion yields: µe + µp = µn + µν . The chemi-
cal potential of neutrinos is a parameter depend-
ing on how much neutrinos and anti-neutrinos are
trapped, and assuming that the number densities
of the trapped neutrinos and anti-neutrinos are
the same, µν can be set to zero. In order to de-
termine the distribution function of the partially
trapped neutrinos, in principle one should solve
the Boltzmann equation. To simplify this prob-
lem, we use here a ”gray body” model, and we
introduce a blocking factor b =
i=e,µ,τ bi to de-
scribe the extent to which neutrinos are trapped
(see e.g. Sawyer 2003). In terms of this factor, we
write the distribution function of neutrinos as
f̃νi(p) =
exp(pc/kT ) + 1
= bifνi , (0 ≤ bi ≤ 1).
This simplified assumption is consistent with the
two-stream approximation which we adopt here
(Eq. 9).
2.3. Composition and chemical equilib-
The equilibrium state of the gas in the accret-
ing torus is completely determined by the chem-
ical potentials of neutrons, protons and electrons
(ηn, ηp, ηe), and the trapping factor of neutrinos
(b) which is related to the optical depths of neu-
trinos (cf. Eq. 9). For a given baryon number
density, nb, temperature T , and a value for accre-
tion rate Ṁ and viscous constant α, the chemical
potentials, or equivalently the ratio of free protons
x = np/nb, are determined from the condition of
equilibrium between the transition reactions from
neutrons to protons and from protons to neutrons.
These reactions are:
p + e− → n + νe (11)
p + ν̄e → n + e+ (12)
p + e− + ν̄e → n (13)
n + e+ → p + ν̄e (14)
n → p + e− + ν̄e (15)
n + νe → p + e− (16)
Therefore we have to calculate the ratio of protons
that will satisfy the balance:
np(Γp+e−→n+νe + Γp+ν̄e→n+e+ + Γp+e−+ν̄e→n)
= nn(Γn+e+→p+ν̄e + Γn→p+e−+νe + Γn+νe→p+e−) .(17)
The reaction rates are the sum of forward and
backward rates and are given in the Appendix (see
also Kohri, Narayan and Piran 2005).
These are supplemented by two additional con-
ditions: the conservation of the baryon number,
nn +np = nb×Xnuc, and charge neutrality (Yuan
2005):
ne = ne− − ne+ = np + n0e , (18)
which says that the net number of electrons is
equal to the number of free protons plus the num-
ber of protons in helium:
n0e = 2nHe = (1 −Xnuc)
. (19)
The number density of fermions under arbitrary
degeneracy is determined by the following equa-
tions:
F1/2(ηi, βi) + βiF3/2(ηi, βi)
Finally, the electron fraction is defined as:
ne− − ne+
(Note that this is different from Ye = 1/(1 +
nn/np), which is only valid for free n-p-e gas.)
2.4. Neutrino cooling
The processes that are responsible for the neu-
trino emission in the disc are electron-positron pair
annihilation (e− + e+ → νi + ν̄i), bremsstrahlung
(n + n → n + n + νi + ν̄i), plasmon decay (γ̃ →
νe + ν̄e) and URCA process (reactions 11, 14 and
15). The first two processes produce neutrinos of
all flavors, while the other produce only electron
neutrinos and anti-neutrinos.
The cooling rate due to pair annihilation is ex-
pressed as:
qe+e− = qνe + qνµ + qντ (22)
where the cooling rates for all three neutrino fla-
vors are calculated by means of Fermi-Dirac inte-
grals and are given in the Appendix.
The cooling rate due to nucleon-nucleon bremsstrahlung
(in erg/cm3/s) is given by:
qbrems = 3.35 × 1027ρ210T 5.511 , (23)
where ρ10 is the baryon density in units of 10
g/cm3 and T11 is temperature in units of 10
11 K.
The cooling rate due to the plasmon decay (in
erg/cm3/s) is:
qplasmon = 1.5×1032T 911γ6pe−γp(1+γp)
1 + γp
where γp = 5.565 × 10−2
(π2 + 3(µe/kT )2)/3.
The cooling rate due to the URCA reactions is
given by the three emissivities:
qurca = qp+e−→n+νe + qn+e+→p+ν̄e + qn→p+e−+ν̄e .
The emissivities are given in the Appendix.
Note that the blocking factor of the trapping
neutrinos is used only for the emissivities of the
URCA reactions. For simplicity, we neglect the
blocking effects of neutrinos when calculating the
emissivities for the electron-positron pair annihi-
lation. Two reasons make this approximation rea-
sonable: the emissivities for the electron-positron
pair annihilation is much smaller than those of the
URCA reactions, and the electron-positron pair
annihilation does not change the electron fraction
which sensitively affects the EOS.
Each of the above neutrino emission process has
a reverse process, which leads to neutrino absorp-
tion. These are given by Equations 12, 13 and
16. Therefore we introduce the absorptive optical
depths for neutrinos given by:
τa,νi =
qa,νi (26)
where absorption of the electron neutrinos is de-
termined by:
qa,νe = q
+ qurca + qplasm +
qbrems (27)
and for the muon neutrinos:
qa,νµ = q
qbrems . (28)
In addition, the free escape of neutrinos from
the disc is limited by scattering. The scattering
optical depth is given by:
τs = τs,p + τs,n (29)
= 24.28 × 10−5
H (Cs,pnp + Cs,nnn)
where Cs,p = (4(CV − 1)2 + 5α2)/24, Cs,n = (1 +
5α2)/24, CV = 1/2 + 2 sin
2 θC, with α = 1.25 and
sin2 θC = 0.23.
The neutrino cooling rate is then given by
Q−ν =
i=e,µ
τa,νi+τs
3τa,νi
. (30)
2.5. Energy and momentum Conservation
The hydrodynamic equations we solve to cal-
culate the disc structure are the standard mass,
energy and momentum conservation.
Making use of the standard disk equations, the
vertically integrated viscous heating rate (per unit
area) over a half thickness H is given by:
Ftot =
3GMṀ
f(r) (31)
where the Newtonian boundary condition is as-
sumed: f(r) = 1 −
rmin/r. Note that in the
time-dependent calculations, instead of Eq. 31, we
will solve the viscous diffusion equation (Eq. 44).
Using mass and momentum conservation Ṁ =
4πρRHvr ≈ 6πνρH where vr ≈ (3ν)/(2r) and
ν = (2Pα)/(3ρΩ) is the kinematic viscosity. The
viscous heating rate can be written in terms of α,
Q+visc =
αΩHP. (32)
Cooling in the disc is due to advection, radia-
tion and neutrino emission. The advective cooling
in a stationary disc is determined approximately
= ΣvrT
= qadv
S , (33)
where qadv ∝ d lnS/d ln r ∝ (d lnT/d ln r − (Γ3 −
1)d ln ρ/d ln r) ≈ const and we adopt the value
of 1.0. The entropy density S is the sum of four
components:
S = Snucl + SHe + Srad + Sν . (34)
The entropy density of the gas of free protons,
neutrons and electron-positron pairs is given by:
Snucl = Se− + Se+ + Sp + Sn (35)
where
(ǫi + Pi) − niηi (36)
(~c)3
F3/2(ηi, βi) + βiF5/2(ηi, βi)
is the energy density of electrons, positrons, pro-
tons or neutrons, Pi is their pressure, given by
Equation 4, ni are the number densities and ηi are
the chemical potentials.
The entropy density of helium is given by:
SHe = nHe
log(mHe
− lognHe)
for nHe > 0.
The entropy density of radiation is
Srad = 4
, (39)
while for neutrinos we have
Sν = 4
. (40)
In the initial disc configuration we assume that
qadv is approximately constant and of order of
unity, but in the subsequent time-dependent evo-
lution the advection term is calculated with the
appropriate radial derivatives.
For the case of photon and electron-positron
pairs in the plasma the radiative cooling is equal
3Pradc
11σT 4
where we adopt the Rosseland-mean opacity κ =
0.4 + 0.64 × 1023ρT−3 [cm2g−1].
An important term in the cooling and heating
balance in the disc is due to photodisintegration
of α particles, with rate:
Qphoto = qphotoH (42)
where
qphoto = 6.28 × 1028ρ10vr
dXnuc
and Xnuc is given by Equation 7. Finally, in order
to calculate the initial stationary configuration, we
solve the energy balance: Ftot = Q
+ Q−ν + Qphoto.
2.6. Time evolution
After solving for the initial disc configuration,
we allow the density and temperature to vary with
time. We solve the time-dependent equations of
mass and angular momentum conservation in the
disc:
3r1/2
(r1/2νΣ)
and the energy equation:
4 − 3χ
12 − 9χ
12 − 9χ
(Q+ −Q−).
where χ = (P − Prad)/P . The cooling term Q−
consists of radiative and neutrino cooling, given by
Equations (41) and (30). Advection is included in
the energy equation via the radial derivatives. The
cooling term due to photodisintegration of helium
now must be proportional to the full time deriva-
tive of Xnuc (cf. Eq. 43) :
Qphoto ∝ vr
∂Xnuc
∂Xnuc
. (46)
2.7. Numerical method
The initial configuration of the disc is calcu-
lated by means of the Newton-Raphson method,
iterated with the hydrostatic equilibrium condi-
tion. We interpolate over the matrix of pre-
calculated results for the equation of state (pres-
sure and entropy) and neutrino cooling rate (
the number of points is 1024x1024). The Fermi-
Dirac integrals are calculated using the mixture of
Gauss-Legendre and Gauss-Laguerre quadratures
(Aparicio 1998).
Having determined the initial radial profiles of
density and temperature, as well as the other
quantities at time t = 0, we start the time evolu-
tion of the disc. We solve the set of Equations (44),
(46) and (46) using the convenient change of vari-
ables y = 2r1/2 and Ξ = yΣ, at fixed radial grid,
equally spaced in y (see Janiuk et al. 2002 and ref-
erences therein). The number of radial zones is set
to 200, which we found to be an adequate resolu-
tion. After determining the solutions for the first
100 time steps by the fourth-order Runge-Kutta
method, we use the Adams-Moulton predictor-
corrector method with an adaptive time step. The
code used an explicit communication model that
is implemented with the standardized MPI com-
munication interface and can be run on multipro-
cessor machines.
We choose the no-torque inner boundary con-
dition, Σin = Tin = 0 (see Abramowicz & Kato
1989). The outer boundary of the disc is done by
adding an extra “dead-zone” to the computational
domain, which accounts for the disk expansion and
conservation of angular momentum.
3. Results
We first analyze the pressure, entropy and neu-
trino cooling rate distributions for a given tem-
perature and baryon density in the gas. Then,
we show the disc structure for a converged static
disc model and finally we show examples of time
evolution of the neutrino luminosity, density, tem-
perature and electron fraction for given sets of pa-
rameters.
3.1. EOS solutions for a given tempera-
ture and density
In Figure 1 and 2 we plot the results of the nu-
merically calculated equation of state for the hot
and dense matter. The plots show the dependence
of the electron fraction, pressure, entropy and neu-
trino cooling rate on temperature and density, re-
spectively.
In the upper panels, we show the neutrino cool-
ing rate. At low temperatures, below T = mec
5 × 109K, there are almost no positrons and free
nucleons. Therefore the neutrino emission pro-
cesses switch off, and the cooling of the gas is
either due to advection, or, when the matter be-
comes transparent to photons, radiative cooling
overtakes. For larger temperatures, the neutrino
emission rate increases up to the temperature of
about ∼ 5 × 1011 K. For very high tempera-
tures, the optical depths for neutrinos increase
very rapidly (τ ∝ T 5, see Eq. (7) in Di Matteo
et al. 2002). Therefore the neutrino cooling rate
decreases at high temperatures (Eq. 30). On the
other hand, for a given temperature (e.g. T ∼ 1011
K), the neutrino cooling rate does not sensitively
depend on density. It varies by one order of magni-
tude in the range of 108 ≤ ρ ≤ 1014 g/cm3, where
the optical depth is τ ∼ 100. The middle panels
of Figures 1 and 2, show the entropy and pressure
as a function of temperature and density. At low
temperatures, the entropy of gas is not important.
The highly degenerate electrons do not give contri-
bution to the entropy, while they are a dominant
term in the pressure, which is therefore indepen-
dent of temperature up to T ∼ 5 × 1010 K. When
the temperature increases, helium becomes disin-
tegrated into free nucleons at energy comparable
to the binding energy of helium, and after that the
radiation (including photons and electron-positron
pairs) contributes mainly to the total entropy and
pressure. Therefore, both these quantities rise
with temperature. At high densities, the entropy
is dominated by neutrons. Finally, the electron
fraction is shown in the bottom panel of Figures 1
and 2. At low temperatures, the electron fraction
is equal to 0.5, it then decreases sharply as the
helium nuclei become disintegrated. As the tem-
perature further increases, positrons appear as the
electrons become non-degenerate. The positron
capture again increases the electron fraction (see
Fig.1).
The electron fraction changes significantly as
a function of density for T > 1010K (in Fig.2,
T = 1011K). At low densities, the torus consists
of free neutrons and protons and Ye is close to 0.5
(see also Eq. 7). As density increases, Ye decreases
to satisfy the beta-equilibrium among the free n-p-
e gas. Above some density (when the temperature
is high enough, e.g. for Fig. 2, ρHe ≈ 1013g cm−3)
helium starts forming. Therefore Ye has a kink
and stars rising steeply, asymptotically approach-
ing 0.5 as the torus consists of plenty of ionized
helium and some electrons to keep charge neutral-
3.2. The steady-state disc structure
In Figures 3 and 4 we show the profiles of den-
sity and temperature in the stationary accretion
disk model for three accretion rates: 1 M⊙/s,
10 M⊙/s and 12 M⊙/s. In general, the temper-
ature and density profiles both increase inward.
However, for Ṁ = 12M⊙/s, a distinct branch of
solutions is reached, which appears different than
the so-called “NDAF” branch (see Kohri & Mi-
neshige 2002). The density and temperature pro-
files for this high accretion rate differ also from
what was found in previous work (Di Matteo et
al. 2002; Janiuk et al. 2004). Due to a more
detailed equation of state, in which we allow for
Fig. 1.— The dependence of the electron frac-
tion (bottom), pressure (middle bottom), en-
tropy (middle upper) neutrino cooling rate (up-
per panel) on temperature, for the constant den-
sity ρ = 1012 g/cm3. The accretion rate is Ṁ =
1 M⊙/s. The pressure and neutrino cooling are in
cgs units and the entropy is in units of kB cm
Fig. 2.— The dependence of the electron frac-
tion (bottom), pressure (middle bottom), entropy
(middle upper) and neutrino cooling rate (upper
panel) on density, for the constant temperature
T = 1011 K. The accretion rate is Ṁ = 1 M⊙/s.
The pressure and neutrino cooling are in cgs units
and the entropy is in units of kBcm
a partial degeneracy of nucleons and electrons as
well as neutrino trapping, our solutions reach den-
sities as high as 1012 g/cm3 in the innermost radii
of the disc. The temperature in this inner disc part
is in the range 4×1010−1.25×1011 K, depending
on the accretion rate. For the hottest disk model,
a local peak in the density forms around 7 − 8RS,
while below that radius the density decreases. Be-
tween ∼ 3.5 and 7 RS, the plasma becomes much
hotter and less dense than outside of this region.
This means that the macroscopic state of the sys-
tem is different here due to an abrupt change in
the heat capacity. In order to check what is the
reason for this transition, we investigate the pres-
sure distribution in the disk.
The profile of the pressure is shown in Fig. 5.
The dominant term in the total pressure is due
to the nucleons, while the radiation pressure (in-
cluding electron-positron pairs) is always several
orders of magnitude smaller. The neutrino pres-
sure is large in the inner disc, once it gets opti-
cally thick to neutrinos (i.e. for Ṁ ≥ 10M⊙/s).
A significant contribution to the pressure is due
to helium at densities high enough for helium to
form, albeit at temperatures low enough such that
its nuclei are not fully disintegrated. For the
largest accretion rate shown, in the region of the
temperature excess and inverse density gradient
(3.5 − 7RS), the total pressure distribution flat-
tens. The helium pressure is now vanishingly small
due to the complete photodisintegration, and the
nuclear pressure is slightly decreased due to the
composition change: smaller number density of
neutrons and larger number density of protons.
The substantial contribution to the pressure is now
given by the neutrinos (large optical depths; see
below) and radiation pressure (increased number
of electron-positron pairs). From the comparison
of Figures 3, 4 and the bottom panel of Figure
5, it can be seen that the total pressure becomes
locally correlated with temperature and anticor-
related with density, thus consituting an unstable
phase.
In Figure 6 we show the neutrino optical depths
due to scattering and absorption. The total opti-
cal depth in the outer disc is typically dominated
by scattering processes, while in the inner disc ab-
sorption processes take over for very high accre-
tion rates. For Ṁ = 1M⊙/s only the very in-
ner disk radii have optical depth close to 1. For
Fig. 3.— The baryon density as a function of
the disc radius, calculated in the stationary so-
lution. The accretion rate is Ṁ = 1 M⊙/s (solid
line), Ṁ = 10 M⊙/s (long dashed line) and Ṁ =
12 M⊙/s (short dashed line) .
Fig. 4.— The temperature as a function of the disc
radius, calculated in the stationary solution. The
accretion rate is Ṁ = 1 M⊙/s (solid line), Ṁ =
10 M⊙/s (long dashed line) and Ṁ = 12 M⊙/s
(short dashed line) .
Fig. 5.— The pressure components as a func-
tion of the disc radius, calculated in the station-
ary solution for the three accretion rate values:
Ṁ = 1 M⊙/s (upper panel), Ṁ = 10 M⊙/s (mid-
dle panel) and Ṁ = 12 M⊙/s (bottom panel)
. The total pressure is marked by the solid
line, and its components are: nuclear (gas) pres-
sure (long dashed line), radiation pressure (short
dashed line), helium pressure (dotted line) and
neutrino pressure (dot-dashed line).
Fig. 6.— The neutrino optical depths due to scat-
tering (τs, solid line) and absorption (τa,e for elec-
tron neutrinos, long dashed line, and τa,µ for muon
neutrinos, short dashed line) as a function of ra-
dius for Ṁ = 1M⊙/s (upper panel), Ṁ = 10M⊙/s
(middle panel) and for Ṁ = 12M⊙/s (bottom
panel). The sum of the three quantities is the
total optical depth (τtot, dotted line).
Ṁ = 12M⊙/s, in the radial strip of ∼ 3.5 − 7RS
the disk is optically thick with absorptive optical
depth for electron neutrinos exceeding the scatter-
ing term and reaching values of the order of 100.
3.2.1. Composition and Chemical potentials
In Figure 7 we show the distribution of the
reduced chemical potentials of protons, electrons
and neutrons throughout the disc. Reduced elec-
tron chemical potentials much larger than unity
(indicating strong electron degeneracy) are found
in the inner disc parts for Ṁ = 10M⊙/s and
Ṁ = 12M⊙/s, whereas for 1 M⊙/s electrons are
only slightly degenerate. For the highest accretion
rate, the maximum degeneracies correspond to the
radius of the local peak in the density (cf. Fig. 3)
and the excess of helium number density (cf. Fig.
8). Below this radius, the species become non-
degenerate again, contributing to the increase of
the electron fraction (cf. Fig. 9). In Figure 8 we
plot the mass fraction of free nucleons as a function
of radius for Ṁ = 12M⊙/s, 10 and 1M⊙/s. As the
Figure shows, in the outer regions, Xnuc increases
as the radius decreases, while the temperature and
density increase (Fig. 3 and Fig. 4). Consistent
with the behavior of Ye (Fig.2), Xnuc subsequently
turns around (decreases) at radii where the den-
sity is high enough for significant helium forma-
tion. This trend is reversed sharply for highest
accretion rates, when the temperature in the disk
is high enough (Fig.4) for helium to be fully dis-
sociated. In consequence, the number density of
alpha particles increases at ∼ 7−12RS and sharply
decreases at lower radii. A similar, but far less pro-
nounced fluctuation in Xnuc is seen at smaller radii
for the case of Ṁ = 10M⊙/s. For smallest accre-
tion rate, 1M⊙/s, there is no helium throughout
the disk.
In Figure 9 we show the radial distribution of
the electron fraction throughout the disc for 1,
10 and 12 M⊙/s. For the case of 1 M⊙/s (solid
line), the electron fraction decreases inward in the
disc as the electrons are captured by protons (in
neutronization reactions). Once the electrons be-
come non-degenerate, positrons appear, and the
positron capture by neutrons again increases the
electron fraction. For the hotter plasma (accre-
tion rate of 10 M⊙/s, dashed line), consistently
with the behavior discussed for Xnuc, helium nu-
clei form as the density becomes high enough be-
Fig. 7.— The chemical potentials of neutrons
(ηn, solid line), electrons (ηe, dashed line) and
protons (ηp, dotted line) as a function of the
disc radius, calculated in the stationary solution.
The accretion rate is Ṁ = 1 M⊙/s (triangles),
Ṁ = 10 M⊙/s (squares) and Ṁ = 12 M⊙/s (cir-
cles).
Fig. 8.— The mass fraction of free nucleons
as a function of radius for Ṁ = 1M⊙/s (solid
line), Ṁ = 10M⊙/s (long dashed line) and Ṁ =
12 M⊙/s (short dashed line).
low ∼ 20RS and Ye increases. For the accretion
rate of 12 M⊙/s there is the sharp decrease in Ye,
at ∼ 7−8RS, due to the sudden dissociation of he-
lium. As helium is fully photo-dissociated, there is
an almost equal number of neutrons and protons
due to the balance of the electron and positron
capture. This implies an electron fraction of 0.5.
At the innermost radius, the temperature and den-
sity drop due to the boundary condition, which
affects the behaviour of both Ye and Xnuc.
3.2.2. Cooling and heating rates
In Figure 10 we plot the rates of viscous heat-
ing, advection and cooling due to neutrino emis-
sion and photo-dissociation in the stationary disc.
The accretion rates are Ṁ = 1 M⊙/s (upper
panel), Ṁ = 10 M⊙/s (middle panel), and Ṁ =
12 M⊙/s (lower panel). For the highest accretion
rates, in the innermost disc the neutrino cooling
rate decreases substantially with respect to the
cooling by photodissociacion. This is because the
neutrinos are trapped in the disc due to a large
opacity The smaller the accretion rate, the less
important is the neutrino trapping effect. This
implies that for an accretion rate of ≤ 10 M⊙/s
neutrinos can escape from the innermost disc.
The advective term is a couple orders of magni-
tude smaller than the other terms. The photodis-
sociacion term is negligible for an accretion rate of
1 M⊙/s, since there is no helium in the whole disc,
and Qphoto is equal to zero by definition. For an
accretion rate of 10 M⊙/s there is very little he-
lium down to about 15-20 RS, and therefore Qphoto
is much smaller than other terms. For the accre-
tion rate of 12 M⊙/s , down to 6 − 10RS in the
region of the disc of high density and maximum
degeneracy, helium nuclei form. The nucleosyn-
thesis of alpha particles leads to the plasma heat-
ing instead of cooling, and therefore the relevant
term in the energy balance has a negative value.
Outward, above ∼ 10RS, there is some fraction
of helium which can be photo-dissociated, so the
cooling term due to this reaction is also important
in the total energy balance. In the inner region
helium is fully dissociated and Qphoto is equal to
zero, increasing again only near the inner bound-
ary due to the local density increase and decrease
of temperature.
Fig. 9.— The electron fraction as a function
of the disc radius, calculated in the stationary
solution. The accretion rate is Ṁ = 1 M⊙/s
(solid line), Ṁ = 10 M⊙/s (long dashed line) and
Ṁ = 12 M⊙/s (short dashed line).
Fig. 10.— The heating and cooling rates due
to photodissociation and neutrino emission (solid
lines) as a function of radius for Ṁ = 1M⊙/s (up-
per panel), Ṁ = 10M⊙/s (middle panel) and for
Ṁ = 12M⊙/s (bottom panel). The other terms
are: cooling rate due to advection (long dashed
line) and viscous heating rate (short dashed line).
3.3. Stability analysis: instabilities at
high-Ṁ
The disc is thermally unstable if d logQ+/d logT >
d logQ−/d logT . Then any small increase (de-
crease) in temperature leads to a heating rate
which is more (less) than the cooling rate, and
as a consequence a further increase (decrease) of
the temperature. The viscous instability, which
appears when ∂Ṁ
∂Σ |Q+=Q− < 0, manifests itself
in a faster (slower) evolution of an underdense
(overdense) region. The instabilities can be con-
veniently located in the surface density - temper-
ature diagrams, in which the branch of thermal
equilibrium solutions with a negative slope is not
only unstable to the perturbations in the surface
density, but it is also thermally unstable.
In Figure 11 we show such stability curves for
several radii in the disc. The criterion for a vis-
cously stable disc is generally satisfied through-
out the whole disc for Ṁ ≤ 10M⊙/s. How-
ever, for larger accretion rates, there are unstable
branches at the smallest radii. For Ṁ = 10M⊙/s,
the disc becomes unstable below 5 RS, while for
Ṁ = 12M⊙/s the instability strip is up to ∼ 7RS.
Here helium is almost completely photodisinte-
grated while the electrons and protons become
non-degenerate again. For this high accretion rate,
the electron fraction rises inward in the disc. Un-
der these conditions, the energy balance is affected
leading to the thermal and viscous instability, as
demonstrated by the stability curves. This insta-
bility will be discussed in more detail in Section
3.4. Time dependent solutions
In this Section, we discuss how the tempera-
ture, density, electron fraction and disk luminosity
evolve with time. In Figures 12 and 13 we show the
time evolution of density and temperature, when
the initial accretion rate is 1 M⊙/s. These quan-
tities exponentially decrease with time:
ρ = ρ0(r) exp(−at) , (47)
T = T0(r) exp(−bt) (48)
where a ≈ 1.9 and b ≈ 0.085. The normalization
of these relations depends on the radius, and for
example for r = 6RS it is ρ0 = 2.2 × 1011 and
T0 = 3.5× 1010. The exponential behaviour arises
from the nature of energy equation (45).
In Figure 14 we show the electron fraction as
a function of time for several exemplary radial lo-
cations in the disc, for the disc evolving from a
starting accretion rate of 1M⊙/s. The fraction Ye
is smaller in the inner disc radii, while outward,
the electron fraction is over half an order of mag-
nitude higher. Altogether, during the evolution
of the system, the electron fraction constantly in-
creases with time throughout the disc.
The time-dependent neutrino luminosity of the
disc is given by:
Lν(t) =
∫ Rmax
Q−ν (t)2πrdr (49)
where Q−ν is given by Equation (30).
In Figure 15 we show an example of such a
lightcurve, for our standard model parameters
(M = 3M⊙, α = 0.1), and Rmax = 50RS. The
starting accretion rate was Ṁstart = 1 M⊙/s. At
this accretion rate neutrinos can already escape
from the accretion disc at the beginning of the
evolution. For higher initial accretion rates, e.g
10 − 12M⊙/s, neutrinos are trapped in the inner-
most disc, and, as a consequence, the neutrino
luminosity is lower at the initial stages of disc
evolution until the accretion rate drops to about
∼ 1M⊙/s. This result is qualitatively similar, al-
beit it differs quantitatively, from what was ob-
tained in Di Matteo et al. (2002) and Janiuk et
al. (2004): in those calculations neutrino trap-
ping was far more substantial even for a ’mod-
erate’ accretion rate of Ṁ & 1M⊙/s. The differ-
ence arises from the fact that here we calculate the
neutrino opacities using the β reaction efficiencies,
self-consistently with the equation of state.
For an accretion rate of 1M⊙/s, the solution
does not reach the viscously unstable branch. Ini-
tially, the disc contains almost no alpha particles
(cf. Fig. 8), which appear later on during the evo-
lution and cooling of the plasma. The dynamical
balance between the photodisintegration of helium
and nucleosynthesis leads to an additional non-
zero cooling/heating term in the energy equation
and to only small amplitude flickering at the early
stages of time-evolution.
The situation is much more dramatic when the
starting accretion rate is 12M⊙/s. In this case a
Fig. 11.— The stability curves on the accretion
rate vs. surface density plane, for several chosen
radii in the disc: 3.39RS (solid line), 3.81RS (dot-
ted line), 4.25RS (short dashed line), 5.19RS (long
dashed line) and 8.60RS (dot-dashed line).
Fig. 12.— The density as a function of time, for
several chosen disc radii: 4.01, 6.04, 10.13, 20.7,
35.02, and 45.19 RS. The initial accretion rate is
Ṁ = 1 M⊙/s.
Fig. 13.— The temperature as a function of time,
for several chosen disc radii: 4.01, 6.04, 10.13,
20.7, 35.02, and 45.19 RS. The initial accretion
rate is Ṁ = 1 M⊙/s.
Fig. 14.— The electron fraction as a function
of time, for several chosen disc radii: 4.01, 6.04,
10.13, 20.7, 35.02, and 45.19 RS. The initial ac-
cretion rate is Ṁ = 1 M⊙/s.
large disc strip is viscously and thermally unstable
and the most violent instability takes place around
and below 7 − 12RS.
In Figure 16 we show the behavior of the lo-
cal accretion rate in the unstable disc, at sev-
eral chosen locations within the instability strip.
Near ∼ 12RS, the accretion rate varies due to
the large and rapidly changing photodisintegra-
tion term (locally, it can become larger than the
neutrino cooling rate).
This radius corresponds to the largest local
value of the density of helium (cf. Figure 8 show-
ing its starting model distribution), which is then
being photodissociated. The photodissociacion
process is the cause of the local rapid accretion
rate changes. Then, inside from this highly vari-
able strip, the accretion rate grows too fast to pre-
serve the disc structure. This kind of behaviour
occurs in the locally hotter and less dense region
visible in the starting configuration e.g. in Figures
3 and 4, between 3.5 and 7 RS. In this region the
helium is already totally photodissociated. Due
to the growing accretion rate all the material is
rapidly accreted onto the black hole and the in-
nermost strip of the disc empties.
After the inner strip is destroyed, the outer
parts can still accrete onto the center. As they
approach the black hole, their temperature and
density grow and the above situation can repeat
several times, until the whole disc is completely
broken into rings and destroyed. These later in-
jections of energy, with timescales dictated by the
viscous timescale of each ring, can produce en-
ergy flares following the main GRB activity. Our
results therefore provide another physical mech-
anism1 for the flare model recently proposed by
Perna et al. (2006).
In Figure 17 we show the neutrino lightcurve of
the unstable disc. The instabilities due to photo-
disintegration are reflected in oscillations of vari-
able amplitude and millisecond timescale. This is
of a particular interest if the neutrino annihila-
tion provides the energy input for GRBs, however
it should be pointed out that the oscillations ap-
pearing in the presented lightcurve have a much
1In addition to the gravitational instability in the outer parts
of the disk, which was hinted by the calculations of Di
Matteo et al. (2002) and confirmed by those of Chen &
Beloborodov (2006).
Fig. 15.— The neutrino lightcurve, integrated
over the disc surface. The initial accretion rate
is Ṁ = 1 M⊙/s.
Fig. 16.— The local accretion rate, as a function of
time, for several chosen radial locations in the disc:
6.86, 7.73 and 12.85 RS. The starting accretion
rate is Ṁ = 12 M⊙/s.
smaller amplitude than the observed gamma ray
variability.
4. Discussion
4.1. The unstable neutrino-opaque disc
In our calculations we have shown that, for
large accretion rates, the accreting torus becomes
viscously and thermally unstable. We now discuss
the physical origin of the instability.
The unstable branch appears both in the
steady-state solutions and in the subsequent time-
dependent evolutions. In the steady-state case,
for a chosen value of a constant accretion rate,
this can be seen for instance by plotting the ra-
dial profiles of density and temperature (cf. Figs.
3 and 4, where the distinct branch is found for
the innermost radii), as well as by looking at the
stability curves for a range of accretion rates at
a chosen disc radius (cf. Fig. 11, where the un-
stable inner disc radii exhibit a negative slope in
the curve). In the time-dependent simulations,
the unstable behavior is manifested by the highly
variable accretion rate in certain strips of the disk
and by the subsequent breaking of the disk inside
from these variable strips (cf. Fig. 16). The in-
stability arises from the fact that the accretion
rate rises locally too fast to prevent the disc strip
from emptying, as the material is supplied from
outer strips at much slower rate than it is accreted
inwards. The disc evolves unstably on a viscous
timescale, τvisc = 1/(αΩ) × (r/H)2; for the radii
shown in Figure 16, it is τvisc ≈ 0.05 s (note that
the disc is rather thick, r/H ∼ 2.5, and therefore
the viscous and thermal timescales are close to
each other). Theoretically, in order to find again
a stable solution, the disc would have to increase
the local accretion rate up to about several tens
of M⊙/s within one viscous timescale. However,
this may not be possible if there is not enough
material in the system to support much higher
accretion rates during such violent oscillations.
Therefore the system is unable to be stabilized
and gets broken after a fraction of τvisc . In addi-
tion, the dynamical instability is the source of the
flickering of the local accretion rate at the edge of
the unstable strip.
Let us now discuss in more detail the physical
reason driving this instability. In the inner part of
the disc (below r ∼ 10RS for Ṁ = 12M⊙/s) there
Fig. 17.— The neutrino luminosity, as a func-
tion of time The starting accretion rate is Ṁ =
12 M⊙/s.
Fig. 18.— The total pressure, as a function of
time, for the chosen radial locations in the disc:
7.73 and 6.86 RS. The starting accretion rate is
Ṁ = 12 M⊙/s.
are two important processes, both of which are
incorporated in our equation of state: photodisin-
tegration of helium and neutrino trapping. As it
was already mentioned in Sec. 3.2.1 and can be
seen from Fig. 8, below ∼ 7 − 8RS helium in this
disc is completely photodisintegrated. This part
of the disc is also opaque to neutrinos, as we show
in Fig. 6.
These two mechanisms competitively influence
the electron fraction in the disc (cf. Sec. 3.2, Fig.
6 and Fig. 9). Well outside the unstable strip,
above ∼ 20RS, the electron fraction smoothly de-
creases inwards as positrons appear because of the
neutronization process. Then the scattering opti-
cal depth for neutrinos becomes τs > 1, and the
electron fraction increases again. After photodis-
integration, the electron fraction decreases signif-
icantly from almost 0.3 to much less than 0.1 due
to electron capture. But again, when the disc
becomes optically thick to absorption of electron
neutrinos, the electron fraction gets higher and ap-
proaches almost 0.5.
The total pressure of sub-nuclear matter (cf.
Fig 5) is mainly contributed by electrons, and
therefore it is influenced by the changes in the
electron fraction. In the narrow range of radii
(6.8 < r/RS < 7.8), the pressure decreases due
to photodisintegration. The sudden decrease of
the pressure might drive the dynamical instabil-
ity. (This picture is somewhat similar to that of
the iron core collapse in the core collapse super-
nova explosions: electron capture consumes most
of the electrons and makes the EOS softer, conse-
quently, it triggers the collapse of the iron core.)
However, the transition from neutrino transpar-
ent to opaque disc, and the increase of the elec-
tron fraction due to the beta equilibrium (see also
Yuan & Heyl 2005), are the reason for a steeper
increase of the total pressure of the system.
The same effect can also be observed in the
time-dependent plot (Fig. 18), in which we show
the pressure changes in the characteristic radii
of the unstable part of the disc (cf. Fig. 16).
The pressure decreases with time up to a radius
R = 6.8RS, since the temperature and density
gradually drop, as well as the neutrino opacities,
so the electron fraction gets smaller. Then, in a
strip between ∼ 6.8 and 12 RS, the pressure rises
with time: in fact, when α particles appear, the
electron fraction rises and matter locally piles up,
thus increasing the pressure.
At the border of these radii the disc breaks up,
when the thermal-viscous instability induces an
avalanche-like increase of the local accretion rate
below ∼ 8RS. This happens because the increase
in the pressure causes an excess in the local energy
dissipation rate and the disc heats up, while at ad-
jacent radii the pressure decreases and heating is
insufficient. The system tries to compensate these
temperature gradients by decreasing/ increasing
the temperature in the outer/inner radius, respec-
tively. But since in the unstable mode of the ther-
mal balance this causes a further increase/decrease
of density, the pressure drops further and the disc
heats up in the outer radius, while cools down in
the inner one. As long as it cannot find any sta-
ble track of evolution, the emptying of the inner
strip continues and finally the whole material is
accreted towards the black hole or blown out.
The radial extent of the unstable part depends
on the initial accretion rate and in our model for
Ṁ = 12M⊙/s it is up to ∼ 8RS, for Ṁ = 10M⊙/s
it is up to ∼ 5RS, while for Ṁ = 1M⊙/s it is below
∼ 3.5RS. In the latter case, since the inner radius
is located at ∼ 3RS, the instability hardly affects
the disc. The extension of the instability strip de-
pends also on the mass of the accreting compact
object, and since for lower mass black holes the ac-
cretion disc is generally denser, it reaches a density
∼ 1012 g/cm3 around ∼ 15RS.
Above 25-30 RS, where the plasma is already
optically thin and the evolution is stable, both
the pressure and the accretion rate smoothly drop
with time. The dominant source of cooling of the
disc in this region is the neutrino emission (ad-
vective cooling decreases as the disc transits from
neutrino opaque to transparent). The photodis-
integration term (if non-zero), is usually by 1-2
orders of magnitudes smaller than neutrino cool-
ing, and in the inner disc, up to about 6.8 RS,
there are no helium nuclei and the photodisinte-
gration term is negligible, while at 7.5 RS it has
a value of about Qphoto ∼ 1038 erg/s/cm2 with
rapid fluctuations. These fluctuations induce the
local accretion rate flickering (cf. Fig. 16), on
a timescale and amplitude much smaller than for
the viscous instability.
4.2. Comparison with previous work
The neutrino dominated accretion flow has al-
ready been studied in a number of papers, includ-
ing both 1-D models and multi-D simulations. The
steady-state 1-D models (e.g Popham et al. 1999;
Kohri & Mineshige 2002) assumed the disk opti-
cally thin to neutrinos, and neglected photodis-
integration cooling. Di Matteo et al. (2002) took
these two effects into account, and showed that the
trapped neutrinos dominate the pressure in the
inner region of the hyperaccreting disc, however
their equation of state did not include the numeri-
cal calculation of chemical equilibrium and did not
incorporate the opacities directly in the EOS itera-
tions (see also the time-dependent model of Janiuk
et al. 2004). Kohri, Narayan & Piran (2005) con-
sidered the neutrino opaque disk and the equilib-
rium between neutrons and protons and calculated
the number densities of species by numerically in-
tegrating their distribution functions. However,
these authors calculated the gas pressure from the
ideal gas approximation, and neglected the contri-
bution of helium to the pressure. In all of these
papers the disc occurred to be stable against any
kind of instability.
On the other hand, in their recent work, Chen &
Beloborodov (2006) find that the outskirts of the
disk are gravitationally unstable. The approach
used by these authors provides a detailed treat-
ment of the microphysics which is very similar to
ours; however, some differences between our work
and theirs must be crucial to the development of
viscous and thermal instabilities. One difference
deals with the approximation made for the treat-
ment of transition region between the neutrino-
opaque and the transparent matter. In our work,
we adopted a gray body model, i.e., we introduced
the b factor to describe the distribution function
(c.f. Eq. 10). This assumption is consistent with
the two fluids approximation we have made, which
has recently been studied numerically by Sawyer
(2003), and shown to be appropriate for the con-
ditions of these disks. On the other hand, Chen
& Beloborodov (2006) smoothly connect the opti-
cally thin and thick regimes by means of interpo-
lation. A further difference lies in the description
of the mass fraction of free nucleons. In this work
we use an expression for Xnuc developed by Qian
et al (1996), while in Chen & Beloborodov (2006)
Xnuc is a function of Ye, which couples the nucle-
osynthesis to the electron fraction.
In our calculations we reach the range of densi-
ties and temperatures where the nucleons start to
become partially degenerate. This is accompanied
by the neutrinos being more and more trapped in
the gas and helium being destroyed by photodisso-
ciation. As a result of these calculations, we found
an additional, unstable branch of solutions for the
disc thermal balance.
This supports the recent results of 2-D simu-
lations by Lee, Ramirez-Ruiz & Page (2005), who
found the disc opaque to neutrinos to be thermally
unstable. Their simulations showed that large cir-
culations develop in the accretion flow. Setiawan,
Ruffert and Janka (2005) found small fluctuations
of the accretion rate and neutrino luminosity on
the dynamical timescale, after the 10-20 msec of
relaxation period (note, that in our calculations we
start from the steady-state disc model at a given
accretion rate, thus having no need for a relaxation
to the quasi-steady configuration). The equation
of state used in their work (see also Janka et al.
1999) is based on the work of Lattimer & Swesty
(1991). Given the electron fraction, this EOS as-
sumes the condition of nuclear statistical equilib-
rium without neutrino trapping, but the evolution
of the electron fraction is affected by the asym-
metric neutrino emission from the hot and dense
matter, which is called ‘neutrino leakage scheme’.
The neutrino leakage scheme focuses on the ef-
fects of the neutrino trapping on the net neutrino
emissivities, not on the nuclear statistical equilib-
rium. The equation of state used in the work of
Rosswog et al. (2004) is temperature and compo-
sition dependent, based on the relativistic mean
field theory (Shen et al. 1998a,b), and the neu-
trino cooling is accounted for by the multiflavor
scheme (Rosswog & Liebendoerfer 2003).
In our work, we use an equation of state based
on the β equilibrium, including the contribution
from the trapped neutrino, and neutrino trapping
effects are accounted for by the appropriate opac-
ities. It should be emphasized that most previous
multi-D simulations neglected the effects of neu-
trino trapping on the β equilibrium, as well as the
contribution of the trapped neutrinos to the ther-
modynamical properties of the dense matter. An-
other difference between our treatment of the EOS
and the previous numerical simulations is that we
include the cooling of the photodisintegration of
helium. Even though the original EOS of Lattimer
& Swesty (1991) can provide detailed information
about the composition of the dense matter, this
information was not considered in order to keep
the table of the EOS as small as possible (see e.g.
Ruffert et al. 1996) just for numerical reasons. In
this way, the disintegration cooling had not been
investigated without the information on the com-
position. Our results indicate that photodisinte-
gration significantly affects the energy balance.
4.3. Limitations of our model
We find the thermal-viscous instability to be an
intrinsic property of the disc for extremely large
torus densities (about 1012g cm−3) and high ac-
cretion rates (Ṁ ≥ 10M⊙/s). This is seen both
in the steady-state results (radial profiles of den-
sity and temperature) and in the subsequent time
evolution.
Thermal and viscous instabilities have been
studied in the case of standard accretion discs
around compact objects (Lightman & Eardley
1974; Pringle 1977; Shakura & Sunyaev 1976).
Two main physical processes that lead to disc
instabilities were invoked to explain the time-
dependent behavior of various objects: partial
ionization of hydrogen in the discs of Dwarf No-
vae (e.g. Meyer & Meyer-Hofmeister 1981; Smak
1984) and domination of radiation pressure in the
X-ray transients (e.g Taam & Lin 1984). Such in-
stabilities do not have to lead to a total disc break-
down, but rather to a limit-cycle behavior, if only
an additional (i.e. upper) stable branch of solu-
tions can be found. This might be a hot state with
a temperature above 104 K, or a slim disc, domi-
nated by advection (Abramowicz et al. 1988). In
our 1-D calculation the disc in the GRB central en-
gine is not stabilized but rather breaks down into
rings, as no stable solutions are reached (possibly,
for even higher accretion rates again a stable part
near the black hole could be formed - but these
extremely high accretion rates would not be pro-
duced by any compact merger scenario). There-
fore, instead of a limit-cycle activity, what we
find here are several dramatic accretion episodes
on the viscous timescale. The remaining parts of
the torus will subsequently accrete and, while ap-
proaching the central black hole, will get hotter
and denser, breaking at ∼ 7RS.
Of course, it would be interesting to study
whether such a violent instability would occur also
in the 2D or 3D simulations. This is indeed likely
to be the case, since as the multi-dimensional sim-
ulations of accretion discs show, the instabilities
derived first in 1D are still present in the hydro-
dynamical simulations of flows with non-Keplerian
velocity fields (e.g. Agol et al. 2001; Turner
2004; Ohsuga 2006). Possibly, the instability re-
gion would be located at other (larger) radii if
the calculations included the vertical structure of
the disc: this is dependent on temperature and
density, which above the meridional plane may be
larger than the mean value considered in the ver-
tically averaged model.
We need to note that our 1D calculations do
not take into account the possible effects of non-
radial velocity components in the fluid. For exam-
ple, the inverse composition gradient that leads
to the disk instability, might be stabilized by ro-
tation (e.g. Begelman & Meier 1982; Quataert &
Gruzinov 2000). In the 2-D simulation of Lee et al
(2005) the neutrino opaque disk exhibits circula-
tions in the r-z direction. Such meridional circula-
tions are known to be present in the Keplerian ac-
cretion disks (e.g. Siemiginowska 1988), however
it is unclear if they could always provide a stabiliz-
ing mechanism for the thermal-viscous instability.
Possibly, if the nonradial motions of the flow pro-
vided a stronger stabilizing effect, the disk would
exhibit oscillations in the viscous timescale, with-
out breaking, similarly to the outbursts of Dwarf
Nova disks.
The assumption of the β equilibrium (justified,
as the mixture of protons, electrons, neutrons and
positrons is able to achieve the equilibrium con-
ditions) might also have an effect on this result,
as the equilibrium conditions reduce the heating
and entropy in the gas. In fact, the β equilibrium
condition which is satisfied in the innermost part
of a hyperaccreting disc that is optically thick to
neutrinos, is µn = µp +µe. Once the disc becomes
transparent in its outer part, this condition is no
longer valid. Analytically, it has been derived by
Yuan (2005) that the condition for β equilibrium
in this case is µn = µp + 2µe.
4.4. Observational consequences
Our findings might be relevant for interpret-
ing some recent observations. The flickering due
to the photodisintegration of alpha particles may
lead to a variable energy output on small (millisec-
ond) timescales. The consequence of this may be
variability in the gamma ray luminosity, although
the changes in the local accretion rate may be
spread by viscous effects (in the lightcurve Lν(t)
integrated over the whole surface of the disc, the
millisecond variability is somewhat smeared, and
the amplitudes are not very large). Therefore the
mass accreted by the black hole may not be vary-
ing substantially, while some irregularity in the
overall outflow could help produce internal shocks.
The thermal-viscous instability, if accompa-
nied by the disc breaking, may lead to the sev-
eral episodic accretion events and several re-
brightenings of the central engine on longer
timescales, possibly detected in the later stages
of the evolution. A similar kind of a long-term
activity is possible also if the disk was not com-
pletely broken, but exhibited some large accretion
rate fluctuations on the viscous timescale.
We thank Bożena Czerny, Pawe l Haensel and
Daniel Proga for helpful discussions. We also
thank the anonymous referee for detailed reports
which helped us to improve our model and its
presentation. This work was supported in part
by grant 1P03D 00829 of the Polish State Com-
mittee for Scientific Research and by NASA un-
der grant NNG06GA80G. Y.-F. Y. is partially
supported by Program for New Century Excel-
lent Talents in University, and the National Natu-
ral Science Foundation (10233030,10573016). RP
acknowledges support from NASA under grant
NNG05GH55G, and from the NSF under grant
AST 0507571.
A. Appendix
The neutrino absorption and production rates in the beta processes for all participating particles at
arbitrary degeneracy have been obtained in the previous works (Reddy, Prakash & Lattimer 1998; Yuan
2005). In the subnuclear dense matter with high temperatures, the nucleons are generally nondegenerate,
therefore, the transition reaction rates from neutrons to protons and from protons to neutrons can be
simplified as follows:
Γp+e−→n+νe =
|M |2
dEeEepe(Ee −Q)2fe(1 − befνe), (A1)
Γp+e−←n+νe =
|M |2
dEeEepe(Ee −Q)2(1 − fe)befνe , (A2)
Γn+e+→p+ν̄e =
|M |2
dEeEepe(Ee + Q)
2fe+(1 − befν̄e), (A3)
Γn+e+←p+ν̄e =
|M |2
dEeEepe(Ee + Q)
2(1 − fe+)befν̄e , (A4)
Γn→p+e−+ν̄e =
|M |2
dEeEepe(Q− Ee)2(1 − fe)(1 − befν̄e), (A5)
Γn←p+e−+ν̄e =
|M |2
dEeEepe(Q− Ee)2febefν̄e . (A6)
Here Q = (mn −mp)c2, |M |2 is the averaged transition rate which depends on the initial and final states
of all participating particles, for nonrelativistic noninteracting nucleons, |M |2 = G2F cos2 θC(1 + 3g2A), here
GF ≃ 1.436×10−49 erg cm3 is the Fermi weak interaction constant, θC (sin θC = 0.231) is the Cabibbo angle,
and gA = 1.26 is the axial-vector coupling constant. fe,νe are the distribution functions for electrons and
neutrinos, respectively. The “chemical potential” of neutrinos is generally assumed to be zero. The factor
be reflects the percentage of the partially trapped neutrinos. When neutrinos completely trapped, be = 1.
The corresponding neutrino emissivities for the URCA reactions are given by:
qp+e−→n+νe =
|M |2
dEeEepe(Ee −Q)3fe(1 − befνe), (A7)
qn+e+→p+ν̄e =
|M |2
dEeEepe(Ee + Q)
3fe+(1 − befν̄e), (A8)
qn→p+e−+ν̄e =
|M |2
dEeEepe(Q− Ee)3(1 − fe)(1 − befν̄e). (A9)
The emissivities due to the electron-positron pair annihilation, following the notation of Yakovlev et al
(2001), is written as:
qe−+e+→νi+ν̄i =
C2+νi [8(Φ1U2 + Φ2U1) − 2(Φ−1U2 + Φ2U−1) + 7(Φ0U1 + Φ1U0)
+ 5(Φ0U−1 + Φ−1U0)] + 9C
−νi [Φ0(U1 + U−1) + (Φ−1 + Φ1)U0]
, (A10)
where
= 1.023 × 1023 erg cm−3 s−1, (A11)
C+νi = C
+ C2Ai and C−νi = C
− C2Ai , here CVi and CAi are the vector and axial-vector constants for
neutrinos (CV e = 2 sin
2 θC + 0.5, CAe = 0.5, CV µ = CV τ = 2 sin
2 θC − 0.5 and CAµ = CAτ = −0.5). The
dimensionless functions Uk and Φk (k= −1, 0, 1, 2) in the above equation can be expressed in terms of the
Fermi-Dirac functions:
U−1 =
β3/2F1/2(ηe, βe) (A12)
β3/2[F1/2(ηe, βe) + βeF3/2(ηe, βe)] (A13)
β3/2[F1/2(ηe, βe) + 2βeF3/2(ηe, βe) + β
eF5/2(ηe, βe)] (A14)
β3/2[F1/2(ηe, βe) + 3βeF3/2(ηe, βe) + 3β
eF5/2(ηe, βe) + β
eF7/2(ηe, βe)]. (A15)
Replacing ηe with ηe+ in Uk, we get the corresponding expressions for Φk.
REFERENCES
Aparicio J., 1998, ApJS, 117, 627
Abramowicz M.A., Czerny, B., Lasota J.P.,
Szuszkiewicz E., 1988, ApJ, 332, 646
Abramowicz M.A., Kato S., 1989, ApJ, 336, 304
Agol E., Krolik J., Turner N.J., Stone J.M., 2001,
ApJ, 558, 543
Begelman M.C., Meier D.L., 1982, ApJ, 253, 873
Berger E., et al., 2005, submitted to Nature
(astro-ph/0508115)
Chen, W.-X. & Beloborodov, A. preprint
astro-ph/0607145
Di Matteo T., Perna R., Narayan R., 2002, ApJ,
579, 706
Eichler D., Livio M., Piran T., Schramm D. N.,
1989, Nature, 340, 126
Falcone, A. et al. 2006, ApJ, 641, 1010
Fox D.B., et al., 2005, Nature, 437, 845
Fujimoto S., Kotake K., Yamada S., Hashimoto
M., Sato K., 2006, ApJ, 644, 1040
Gehrels N., et al., 2005, Nature, 437, 851
Gu, W.-M., Liu, T., & Lu, J.-F. 2006, ApJ, 643,
Hjorth, J. et al. 2003, Nature, 423, 847
Hjorth J., et al., 2005, Nature, 437, 859
Janiuk A., Perna R., Di Matteo T., Czerny, B.
2004, MNRAS, 355, 950
Janka H.-T., Eberl T., Ruffert M., Fryer C. L.
1999, ApJ, 527L, 39
Kohri K., Mineshige S., 2002, ApJ, 577, 311
Kohri K., Narayan R., Piran T., 2005, ApJ, 629,
Lattimer J.M., Swesty F.D., 1991, NuPhA, 535,
Lee W.H., Ramirez-Ruiz E., 2002, ApJ, 577, 893
Lee W.H., Ramirez-Ruiz E., Page D., 2005, ApJ,
632, 421
Lightman A.P., Eardley D.M., 1974, ApJL, 187,
Liu, T., Gu, W.-M., Xue, L., & Lu, J.-
F. 2007, ArXiv Astrophysics e-prints,
arXiv:astro-ph/0702186
MacFadyen A. I., Woosley S. E. 1999, ApJ, 524,
Meszaros P., 2006, Rept.Prog.Phys., 69, 2259
Meyer F., Meyer-Hofmeister E., 1981, A&A, 104,
Narayan R., Paczynski B., Piran T., 1992, ApJ,
395, L83
Narayan R., Piran T., Kumar P., 2001, ApJ, 557,
Ohsuga K., 2006, ApJ, 640, 923
Qian Y.Z., Woosley S.E., 1996, ApJ, 471, 331
Quataert E., Gruzinov A., 2000, ApJ, 539, 809
Paczynski B., 1986, ApJ, 308, L43
Paczynski B., 1991, AcA, 41, 257
Paczynski B., 1998, ApJ, 494, L45
Perna R., Armitage P., Zhang B. 2006, ApJ, 636L,
Piran T., 2005, Rev.Mod.Phys., 76, 1143
Popham R., Narayan R., 1995, ApJ, 442, 337
Popham R., Woosley S.E., Fryer C., 1999, ApJ,
518, 356
Pringle J.E., 1977, MNRAS, 177, 65
Proga D., MacFadyen A.I., Armitage P.J., Begel-
mann M.C., 2003, ApJ, 599, L5
Reddy S., Prakash M., Lattimer J.M., 1998, Phys.
Rev. D58, 300
Rosswog S., Liebendoerfer M. 2003, MNRAS, 342,
Rosswog S., Speith, R., Wynn G. A. 2004, MN-
RAS, 351, 1121
http://arxiv.org/abs/astro-ph/0508115
http://arxiv.org/abs/astro-ph/0607145
http://arxiv.org/abs/astro-ph/0702186
Rosswog S. 2005, ApJ, 634, 1202
Ruffert M., Janka H.-T., 1999, A&A, 344, 573
Ruffert M., Janka H.-T., 2001, A&A, 380, 544
Ruffert M., Janka H.-T., Schaefer G., 1996, A&A,
311, 532
Sawyer, R.F. 2003, Physical Review D, vol. 68,Is-
sue 6, id. 06300
Setiawan S., Ruffert M., Janka H.-Th. 2004, MN-
RAS, 352, 753
Setiawan S., Ruffert M., Janka H.-Th. 2005, A&A,
in press (astro-ph/0509300)
Shakura N.I, Sunyaev R., 1976, MNRAS, 175, 613
Shen H., Toki H., Oyamatsu K., Sumiyoshi K.,
1998a, PThPh, 100, 1013
Shen H., Toki H., Oyamatsu K., Sumiyoshi K.,
1998b, NuPhA, 637, 435
Surman R., McLaughlin G.C., 2004, ApJ, 603, 611
Smak J., 1984, AcA, 34, 161
Taam R.E., Lin D.N.C., 1984, ApJ, 287, 761
Siemiginowska A., 1988, AcA, 38, 21
Stanek K. Z., et al., 2003, ApJ, 591, L17
Szuszkiewicz E., 1990, MNRAS, 244, 377
Turner N.J., 2004, ApJ, 605, L45
Villasenor J.S., et al., 2005, Nature, 437, 855
Watarai K., Mineshige S., 2001, PASJ, 53, 915
Witt H. J., Jaroszynski M., Haensel P., Paczynski
B., Wambsganss J., 1994, ApJ, 422, 219
Yakovlev D.G., Kaminker A.D., Gnedin O.Y.,
Haensel P., 2001, Phys. Rep., 354, 1
Yuan Y., 2005, Phys. Rev. D, 72, 013007
Yuan Y., Heyl J.S., 2005, MNRAS, 360, 1493
Zhang B., 2007, astro-ph/0701520
This 2-column preprint was prepared with the AAS LATEX
macros v5.2.
http://arxiv.org/abs/astro-ph/0509300
http://arxiv.org/abs/astro-ph/0701520
	Introduction
	Neutrino-cooled accretion disks
	Initial disc configuration: 1-D Hydrodynamics
	The equation of state
	Composition and chemical equilibrium
	Neutrino cooling
	Energy and momentum Conservation
	Time evolution
	Numerical method
	Results
	EOS solutions for a given temperature and density
	The steady-state disc structure
	Composition and Chemical potentials
	Cooling and heating rates
	Stability analysis: instabilities at high-"705FM
	Time dependent solutions
	Discussion
	The unstable neutrino-opaque disc
	Comparison with previous work
	Limitations of our model
	Observational consequences
	Appendix

  We investigate the properties and evolution of accretion tori formed after
the coalescence of two compact objects. At these extreme densities and
temperatures, the accreting torus is cooled mainly by neutrino emission
produced primarily by electron and positron capture on nucleons (beta
reactions). We solve for the disc structure and its time evolution by
introducing a detailed treatment of the equation of state which includes
photodisintegration of helium, the condition of beta-equilibrium, and neutrino
opacities. We self-consistently calculate the chemical equilibrium in the gas
consisting of helium, free protons, neutrons and electron-positron pairs and
compute the chemical potentials of the species, as well as the electron
fraction throughout the disc. We find that, for sufficiently large accretion
rates (> 10 solar masses per second), the inner regions of the disk become
opaque and develop a viscous and thermal instability. The identification of
this instability might be relevant for GRB observations.

Introduction
Gamma-Ray Bursts are commonly thought to
be produced in relativistic ejecta that dissipate en-
ergy by internal shocks however alternative ideas
based on the Poynting flux dominated jets are also
being proposed (for a review see e.g. Piran 2005;
Meszaros 2006; Zhang 2007).
The enormous power released during the gamma-
ray burst explosion indicates that a relativistic
phenomenon must be involved in creating GRBs
(Narayan et al. 1992). The merger of two neutron
stars or of a neutron star and a black hole (e.g. NS-
NS or NS-BH) has been invoked e.g. by Paczyński
(1986); Eichler et al. (1989); Paczyński (1991);
1Copernicus Astronomical Center, Bartycka 18, 00-716
Warsaw, Poland
2Present address: University of Nevada, Las Vegas, 4505
Maryland Pkwy, NV89154, USA
3 Center for Astrophysics, University of Science and
Technology of China, Hefei, Anhui 230026, P.R. China
4JILA and Department of Astrophysical and Planetary
Sciences, University of Colorado, 440 UCB, Boulder, CO
80309, USA
5Physics Department, Carnegie Mellon University, 5000
Forbes Avenue, Pittsburgh, PA 15232, USA
Narayan, Paczyński & Piran (1992), as well as the
collapse of a massive star, the so called “collap-
sar” scenario (e.g., Woosley 1993 and Paczyński
1998). In both cases a dense, hot accretion disk
is likely to form around a newly born black hole
(Witt et al. 1994). In the collapsar scenario, the
collapsing envelope of the star accretes onto the
newly formed black hole, while a transient debris
disk is formed when the NS-NS or NS-BH binary
merges (see e.g. Ruffert et al. 1997 for numerical
simulations).
The durations of GRBs, which range from mil-
liseconds to over a thousand of seconds, are dis-
tributed in two distinct peaks defining two main
GRB classes: short (≤ 2 sec) and long (& 2 sec)
bursts (Kouvelietou et al. 1993). For some long
bursts, signatures of an accompanying supernova
explosion have been detected in the afterglow spec-
tra (Stanek et al. 2003; Hjorth et al. 2003), which
strongly favors the “collapsar” interpretation for
their origin. Furthermore, the GRB positions in-
ferred from the afterglow observations are consis-
tent with the GRBs being associated with the star
forming regions in their host galaxies. For short
bursts, several pieces of evidence from the analy-
http://arxiv.org/abs/0704.1325v1
sis of the Swift satellite and follow-up observations
(Gehrels et al. 2005; Fox et al. 2005; Villasenor et
al. 2005) argue in favor of a binary merger model
(Hjorth et al. 2005; Berger et al. 2005).
In the merger scenario, the duration of the
burst is comparable to the viscous timescale of the
accretion disc whereas in the collapsar scenario,
the external reservoir of stellar matter can feed
the accretion torus for a much longer timescale.
In general, accretion discs powering GRBs are ex-
pected to have typical densities of the order of
1010−12 g cm−3 and temperatures of 1011 K within
10−20 Schwarzschild radii (RS = 2GM/c2). Thus,
accretion proceeds with rates of a fraction to sev-
eral solar masses per second. In this “hyper-
accreting” regime, photons become trapped and
are not efficient at cooling the disc. Neutrinos,
however, are produced by weak interactions in the
dense and hot plasma, releasing the gravitational
energy of the accretion flow. These discs go under
the name of “neutrino-dominated accretion flows”,
or NDAFs.
Over the last several years a number of stud-
ies have investigated the structure of these discs
(Popham, Woosley & Fryer 1999; Narayan, Piran
& Kumar 2001; Kohri & Mineshige 2002; Di Mat-
teo, Perna & Narayan 2002; Surman & McLaugh-
lin 2004; Kohri, Narayan & Piran 2005; Chen &
Beloborodov 2006; Gu, Liu & Lu 2006; Liu, Gu
& Lu 2007). These models have employed the
customary approximation of one-dimensional hy-
drodynamics (Shakura & Sunyaev 1973), where
the effects of MHD viscous stresses are described
by the dimensionless parameter α, but have been
limited to the steady-state approximation of con-
stant Ṁ . Such an assumption is a good approx-
imation when considering the collapsar scenario,
where the burst duration is much longer than the
viscous time, due to the continuous replenishing of
the disc by the collapsing star envelope. However,
even in this scenario, recent observations suggest
that that engines are “long-lived” (past the torus
feeding phase), requiring a time-dependent com-
putation. In the merger scenario, on the other
hand, a time-dependent calculation is necessary
even for modeling the prompt phase of the burst,
since the duration is set by the viscous timescale
of the disc.
Recently, a fully time-dependent calculation of
the structure of such accretion discs has been pre-
sented by Janiuk et al. (2004). The model was
suitable for the torus being a results of either the
gravitational collapse of a massive stellar core or
the compact binary merger. However the struc-
ture and evolution of the disk (like in many of the
early calculations) was calculated under a num-
ber of simplifying assumptions for the composition
and the equation of state of the accreting matter.
In this paper we improve upon our earlier results
in several ways. Besides the requirement of time-
dependent calculations, the high density and tem-
perature regime in which the accreting gas lies, im-
plies that both multi-dimensional numerical and
semi-analytic calculations for such flows need to
include the detailed microphysics. This includes
photodisintegration of nuclei, the establishment
of statistical equilibrium, neutronization, and the
effects of neutrino opacities in the inner regions.
Here, we introduce a detailed treatment of the
equation of state, and calculate self-consistently
the chemical equilibrium in the gas that consists
of helium, free protons, neutrons and electron-
positron pairs. We compute the chemical poten-
tials of the species, as well as the electron fraction
throughout the disk using the assumption of the
equilibrium between the beta processes. Our EOS
equations include, self-consistently, the contribu-
tion of the neutrino trapping to the beta equilib-
rium. Another important addition compared to
our previous work (Janiuk et al. 2004) is the inclu-
sion of photodisintegration of helium. The pres-
ence of this term can affect the energy balance in
the inner, opaque (to neutrinos as well as pho-
tons) region of the flow and, as it will be shown, it
eventually produces a thermal and viscous insta-
bility in those regions. This is especially relevant
since the GRB phenomenology requires a variable
energy output.
Other time-dependent disc studies of binary
mergers or collapsars have been performed in
2D using hydrodynamical simulations (e.g. Mac-
Fadyen & Woosley 1999; Ruffert & Janka 1999;
Lee et al. 2002; Rosswog et al. 2004; Lee,
Ramirez-Ruiz & Page 2005) and, most recently,
in 3-D simulations (Setiawan et al. 2005). Also,
MHD simulations of the GRB central engine have
been performed, showing that the magnetic field
possibly plays an important role in the genera-
tion of a GRB jet (Proga et al. 2003; Fujimoto
et al. 2006). The advantage of our calculations
is that, whilst including all the relevant physics
to calculate the equation of state, the structure
and stability of the accretion disc, we are able to
study a much larger range of parameter space and
allow our calculations to evolve beyond what can
be reached in higher dimensional calculation and
comparable to at least the short-burst durations.
The paper is organized as follows. In Sec-
tion 2 we describe the basic assumptions of the
model and the method used in the initial sta-
tionary and subsequent time-dependent numerical
simulations. In Section 3 we discuss the structure
of the hyperaccreting disc for various values of the
initial accretion rate, and study the time evolu-
tion of its density and temperature, as well as the
resulting neutrino lightcurve. We also discuss the
physical origin of the instabilities in the disc, and
we compare our model and results with the recent
2D and 3D simulations. We summarize our results
in Section 4.
2. Neutrino-cooled accretion disks
In this Section, we describe how we improve
upon our previous time-dependent calculation (Ja-
niuk et al. 2004) by computing self-consistently
the equation of state of the extremely dense mat-
ter by solving the balance of the β reaction rates.
This allows us to determine the chemical poten-
tials of electrons, protons and neutrons, as well as
the electron fraction, in the initial disc configura-
tion and throughout its evolution.
2.1. Initial disc configuration: 1-D Hydro-
dynamics
We start by considering a steady-state model
of an accretion disc around a Schwarzschild black
hole - formed as a remnant structure either after
a compact binary merger, or in a collapsar after
the birth of a black hole (for a recent calculation
in Kerr spacetime see Chen & Beloborodov 2006).
Throughout our calculations we use the vertically
integrated equations and hence derive a vertically
averaged disc structure. We write the surface den-
sity of the disk as Σ = Hρ, where ρ is the density
and where the disk half thickness (or disk height)
is given by H = cs/ΩK . Here the sound speed
is defined by cs =
P/ρ and ΩK =
GM/r3
is the Keplerian angular velocity with P the to-
tal pressure. We note that, at very high accretion
rates, the disc becomes moderately geometrically
thick (H ∼ 0.5r) in regions where neutrino cooling
becomes inefficient and advection dominates. Our
’slim disk’ approximation neglects terms ∼ (H/r)2
and assumes that the fluid is in Keplerian rota-
tion. For the disc viscous stress we use the stan-
dard α viscosity prescription of Shakura & Sun-
yaev (1973) where the stress tensor is proportional
to the pressure:
τrϕ = −αP. (1)
We adopt a value of α = 0.1
We set the inner radius of the disc at 3 RS, while
the outer radius is at 50 RS. The initial mass of
such a disc is about 0.35M⊙ for an accretion rate
Ṁ = 1 M⊙/s. Throughout the calculations we
adopt a black hole mass of M = 3M⊙.
2.2. The equation of state
We assume that the torus consists of helium,
electron-positron pairs, free neutrons and protons.
The total pressure is contributed by all particle
species in the disc, and the fraction of each species
is determined by self-consistently solving the bal-
ance of the beta reaction rates. In the equation
of state we take into account the pressure due to
the free nuclei and pairs, helium, radiation and the
trapped neutrinos:
P = Pnucl + PHe + Prad + Pν . (2)
The component Pnucl includes free neutrons, pro-
tons, and the electron-positron pair gas in beta
equilibrium:
Pnucl = Pe− + Pe+ + Pn + Pp (3)
(~c)3
F3/2(ηi, βi) +
βiF5/2(ηi, βi)
where Fk are the Fermi-Dirac integrals of the or-
der k, and ηe, ηp and ηn are the reduced chemi-
cal potentials of electrons, protons and neutrons
in units of kT , respectively (where ηi = µi/kT ,
also known as the degeneracy parameter, where µi
the standard chemical potential) calculated from
the chemical equilibrium condition (§ 2.3). The
reduced chemical potential of positrons is ηe+ =
−ηe − 2/βe and the relativity parameters of the
species i are defined as βi = kT/mic
Under the physical conditions in the torus,
helium is generally non-relativistic and non-
degenerate; therefore, its pressure is given by:
PHe = nHekT, (5)
where nHe is the number density of helium. This
is defined as:
nHe =
nb(1 −Xnuc) , (6)
and the fraction of free nucleons is given by
Xnuc = 295.5ρ
11 exp(−0.8209/T11), (7)
with T11 the temperature in unit of 10
11 K (e.g.
Qian & Woosley 1996; Popham et al. 1999).
The radiation pressure is given by:
Prad =
(kT )4
(~c)3
. (8)
When neutrinos become trapped in the disc, the
neutrino pressure is non-zero. Following the treat-
ment of photon transport under the two-stream
approximation (Popham & Narayan 1995; Di Mat-
teo et al. 2002), we have
(kT )4
3(~c)3
i=e,µ,τ
(τa,νi + τs) +
(τa,νi + τs) +
3τa,νi
(kT )4
3(~c)3
b, (9)
where τs is the scattering optical depth due to
the neutrino scattering on free neutrons and pro-
tons and τa,νe and τa,νµ are the absorptive opti-
cal depths for electron and muon neutrinos, re-
spectively (see§2.3). The contribution from tau
neutrinos is the same as that from muon neutri-
nos. These optical depths and neutrino absorption
processes (which are the reverse of the emission
processes) are discussed in more detail in the Ap-
pendix.
In the disc we have to consider both the neu-
trino transparent and opaque regions, as well as
the transition between the two. In the trans-
parent case, the neutrinos are not thermalized
and the chemical potential of neutrinos is negli-
gible. On the other hand, when neutrinos are
totally trapped, the chemical equilibrium condi-
tion yields: µe + µp = µn + µν . The chemi-
cal potential of neutrinos is a parameter depend-
ing on how much neutrinos and anti-neutrinos are
trapped, and assuming that the number densities
of the trapped neutrinos and anti-neutrinos are
the same, µν can be set to zero. In order to de-
termine the distribution function of the partially
trapped neutrinos, in principle one should solve
the Boltzmann equation. To simplify this prob-
lem, we use here a ”gray body” model, and we
introduce a blocking factor b =
i=e,µ,τ bi to de-
scribe the extent to which neutrinos are trapped
(see e.g. Sawyer 2003). In terms of this factor, we
write the distribution function of neutrinos as
f̃νi(p) =
exp(pc/kT ) + 1
= bifνi , (0 ≤ bi ≤ 1).
This simplified assumption is consistent with the
two-stream approximation which we adopt here
(Eq. 9).
2.3. Composition and chemical equilib-
The equilibrium state of the gas in the accret-
ing torus is completely determined by the chem-
ical potentials of neutrons, protons and electrons
(ηn, ηp, ηe), and the trapping factor of neutrinos
(b) which is related to the optical depths of neu-
trinos (cf. Eq. 9). For a given baryon number
density, nb, temperature T , and a value for accre-
tion rate Ṁ and viscous constant α, the chemical
potentials, or equivalently the ratio of free protons
x = np/nb, are determined from the condition of
equilibrium between the transition reactions from
neutrons to protons and from protons to neutrons.
These reactions are:
p + e− → n + νe (11)
p + ν̄e → n + e+ (12)
p + e− + ν̄e → n (13)
n + e+ → p + ν̄e (14)
n → p + e− + ν̄e (15)
n + νe → p + e− (16)
Therefore we have to calculate the ratio of protons
that will satisfy the balance:
np(Γp+e−→n+νe + Γp+ν̄e→n+e+ + Γp+e−+ν̄e→n)
= nn(Γn+e+→p+ν̄e + Γn→p+e−+νe + Γn+νe→p+e−) .(17)
The reaction rates are the sum of forward and
backward rates and are given in the Appendix (see
also Kohri, Narayan and Piran 2005).
These are supplemented by two additional con-
ditions: the conservation of the baryon number,
nn +np = nb×Xnuc, and charge neutrality (Yuan
2005):
ne = ne− − ne+ = np + n0e , (18)
which says that the net number of electrons is
equal to the number of free protons plus the num-
ber of protons in helium:
n0e = 2nHe = (1 −Xnuc)
. (19)
The number density of fermions under arbitrary
degeneracy is determined by the following equa-
tions:
F1/2(ηi, βi) + βiF3/2(ηi, βi)
Finally, the electron fraction is defined as:
ne− − ne+
(Note that this is different from Ye = 1/(1 +
nn/np), which is only valid for free n-p-e gas.)
2.4. Neutrino cooling
The processes that are responsible for the neu-
trino emission in the disc are electron-positron pair
annihilation (e− + e+ → νi + ν̄i), bremsstrahlung
(n + n → n + n + νi + ν̄i), plasmon decay (γ̃ →
νe + ν̄e) and URCA process (reactions 11, 14 and
15). The first two processes produce neutrinos of
all flavors, while the other produce only electron
neutrinos and anti-neutrinos.
The cooling rate due to pair annihilation is ex-
pressed as:
qe+e− = qνe + qνµ + qντ (22)
where the cooling rates for all three neutrino fla-
vors are calculated by means of Fermi-Dirac inte-
grals and are given in the Appendix.
The cooling rate due to nucleon-nucleon bremsstrahlung
(in erg/cm3/s) is given by:
qbrems = 3.35 × 1027ρ210T 5.511 , (23)
where ρ10 is the baryon density in units of 10
g/cm3 and T11 is temperature in units of 10
11 K.
The cooling rate due to the plasmon decay (in
erg/cm3/s) is:
qplasmon = 1.5×1032T 911γ6pe−γp(1+γp)
1 + γp
where γp = 5.565 × 10−2
(π2 + 3(µe/kT )2)/3.
The cooling rate due to the URCA reactions is
given by the three emissivities:
qurca = qp+e−→n+νe + qn+e+→p+ν̄e + qn→p+e−+ν̄e .
The emissivities are given in the Appendix.
Note that the blocking factor of the trapping
neutrinos is used only for the emissivities of the
URCA reactions. For simplicity, we neglect the
blocking effects of neutrinos when calculating the
emissivities for the electron-positron pair annihi-
lation. Two reasons make this approximation rea-
sonable: the emissivities for the electron-positron
pair annihilation is much smaller than those of the
URCA reactions, and the electron-positron pair
annihilation does not change the electron fraction
which sensitively affects the EOS.
Each of the above neutrino emission process has
a reverse process, which leads to neutrino absorp-
tion. These are given by Equations 12, 13 and
16. Therefore we introduce the absorptive optical
depths for neutrinos given by:
τa,νi =
qa,νi (26)
where absorption of the electron neutrinos is de-
termined by:
qa,νe = q
+ qurca + qplasm +
qbrems (27)
and for the muon neutrinos:
qa,νµ = q
qbrems . (28)
In addition, the free escape of neutrinos from
the disc is limited by scattering. The scattering
optical depth is given by:
τs = τs,p + τs,n (29)
= 24.28 × 10−5
H (Cs,pnp + Cs,nnn)
where Cs,p = (4(CV − 1)2 + 5α2)/24, Cs,n = (1 +
5α2)/24, CV = 1/2 + 2 sin
2 θC, with α = 1.25 and
sin2 θC = 0.23.
The neutrino cooling rate is then given by
Q−ν =
i=e,µ
τa,νi+τs
3τa,νi
. (30)
2.5. Energy and momentum Conservation
The hydrodynamic equations we solve to cal-
culate the disc structure are the standard mass,
energy and momentum conservation.
Making use of the standard disk equations, the
vertically integrated viscous heating rate (per unit
area) over a half thickness H is given by:
Ftot =
3GMṀ
f(r) (31)
where the Newtonian boundary condition is as-
sumed: f(r) = 1 −
rmin/r. Note that in the
time-dependent calculations, instead of Eq. 31, we
will solve the viscous diffusion equation (Eq. 44).
Using mass and momentum conservation Ṁ =
4πρRHvr ≈ 6πνρH where vr ≈ (3ν)/(2r) and
ν = (2Pα)/(3ρΩ) is the kinematic viscosity. The
viscous heating rate can be written in terms of α,
Q+visc =
αΩHP. (32)
Cooling in the disc is due to advection, radia-
tion and neutrino emission. The advective cooling
in a stationary disc is determined approximately
= ΣvrT
= qadv
S , (33)
where qadv ∝ d lnS/d ln r ∝ (d lnT/d ln r − (Γ3 −
1)d ln ρ/d ln r) ≈ const and we adopt the value
of 1.0. The entropy density S is the sum of four
components:
S = Snucl + SHe + Srad + Sν . (34)
The entropy density of the gas of free protons,
neutrons and electron-positron pairs is given by:
Snucl = Se− + Se+ + Sp + Sn (35)
where
(ǫi + Pi) − niηi (36)
(~c)3
F3/2(ηi, βi) + βiF5/2(ηi, βi)
is the energy density of electrons, positrons, pro-
tons or neutrons, Pi is their pressure, given by
Equation 4, ni are the number densities and ηi are
the chemical potentials.
The entropy density of helium is given by:
SHe = nHe
log(mHe
− lognHe)
for nHe > 0.
The entropy density of radiation is
Srad = 4
, (39)
while for neutrinos we have
Sν = 4
. (40)
In the initial disc configuration we assume that
qadv is approximately constant and of order of
unity, but in the subsequent time-dependent evo-
lution the advection term is calculated with the
appropriate radial derivatives.
For the case of photon and electron-positron
pairs in the plasma the radiative cooling is equal
3Pradc
11σT 4
where we adopt the Rosseland-mean opacity κ =
0.4 + 0.64 × 1023ρT−3 [cm2g−1].
An important term in the cooling and heating
balance in the disc is due to photodisintegration
of α particles, with rate:
Qphoto = qphotoH (42)
where
qphoto = 6.28 × 1028ρ10vr
dXnuc
and Xnuc is given by Equation 7. Finally, in order
to calculate the initial stationary configuration, we
solve the energy balance: Ftot = Q
+ Q−ν + Qphoto.
2.6. Time evolution
After solving for the initial disc configuration,
we allow the density and temperature to vary with
time. We solve the time-dependent equations of
mass and angular momentum conservation in the
disc:
3r1/2
(r1/2νΣ)
and the energy equation:
4 − 3χ
12 − 9χ
12 − 9χ
(Q+ −Q−).
where χ = (P − Prad)/P . The cooling term Q−
consists of radiative and neutrino cooling, given by
Equations (41) and (30). Advection is included in
the energy equation via the radial derivatives. The
cooling term due to photodisintegration of helium
now must be proportional to the full time deriva-
tive of Xnuc (cf. Eq. 43) :
Qphoto ∝ vr
∂Xnuc
∂Xnuc
. (46)
2.7. Numerical method
The initial configuration of the disc is calcu-
lated by means of the Newton-Raphson method,
iterated with the hydrostatic equilibrium condi-
tion. We interpolate over the matrix of pre-
calculated results for the equation of state (pres-
sure and entropy) and neutrino cooling rate (
the number of points is 1024x1024). The Fermi-
Dirac integrals are calculated using the mixture of
Gauss-Legendre and Gauss-Laguerre quadratures
(Aparicio 1998).
Having determined the initial radial profiles of
density and temperature, as well as the other
quantities at time t = 0, we start the time evolu-
tion of the disc. We solve the set of Equations (44),
(46) and (46) using the convenient change of vari-
ables y = 2r1/2 and Ξ = yΣ, at fixed radial grid,
equally spaced in y (see Janiuk et al. 2002 and ref-
erences therein). The number of radial zones is set
to 200, which we found to be an adequate resolu-
tion. After determining the solutions for the first
100 time steps by the fourth-order Runge-Kutta
method, we use the Adams-Moulton predictor-
corrector method with an adaptive time step. The
code used an explicit communication model that
is implemented with the standardized MPI com-
munication interface and can be run on multipro-
cessor machines.
We choose the no-torque inner boundary con-
dition, Σin = Tin = 0 (see Abramowicz & Kato
1989). The outer boundary of the disc is done by
adding an extra “dead-zone” to the computational
domain, which accounts for the disk expansion and
conservation of angular momentum.
3. Results
We first analyze the pressure, entropy and neu-
trino cooling rate distributions for a given tem-
perature and baryon density in the gas. Then,
we show the disc structure for a converged static
disc model and finally we show examples of time
evolution of the neutrino luminosity, density, tem-
perature and electron fraction for given sets of pa-
rameters.
3.1. EOS solutions for a given tempera-
ture and density
In Figure 1 and 2 we plot the results of the nu-
merically calculated equation of state for the hot
and dense matter. The plots show the dependence
of the electron fraction, pressure, entropy and neu-
trino cooling rate on temperature and density, re-
spectively.
In the upper panels, we show the neutrino cool-
ing rate. At low temperatures, below T = mec
5 × 109K, there are almost no positrons and free
nucleons. Therefore the neutrino emission pro-
cesses switch off, and the cooling of the gas is
either due to advection, or, when the matter be-
comes transparent to photons, radiative cooling
overtakes. For larger temperatures, the neutrino
emission rate increases up to the temperature of
about ∼ 5 × 1011 K. For very high tempera-
tures, the optical depths for neutrinos increase
very rapidly (τ ∝ T 5, see Eq. (7) in Di Matteo
et al. 2002). Therefore the neutrino cooling rate
decreases at high temperatures (Eq. 30). On the
other hand, for a given temperature (e.g. T ∼ 1011
K), the neutrino cooling rate does not sensitively
depend on density. It varies by one order of magni-
tude in the range of 108 ≤ ρ ≤ 1014 g/cm3, where
the optical depth is τ ∼ 100. The middle panels
of Figures 1 and 2, show the entropy and pressure
as a function of temperature and density. At low
temperatures, the entropy of gas is not important.
The highly degenerate electrons do not give contri-
bution to the entropy, while they are a dominant
term in the pressure, which is therefore indepen-
dent of temperature up to T ∼ 5 × 1010 K. When
the temperature increases, helium becomes disin-
tegrated into free nucleons at energy comparable
to the binding energy of helium, and after that the
radiation (including photons and electron-positron
pairs) contributes mainly to the total entropy and
pressure. Therefore, both these quantities rise
with temperature. At high densities, the entropy
is dominated by neutrons. Finally, the electron
fraction is shown in the bottom panel of Figures 1
and 2. At low temperatures, the electron fraction
is equal to 0.5, it then decreases sharply as the
helium nuclei become disintegrated. As the tem-
perature further increases, positrons appear as the
electrons become non-degenerate. The positron
capture again increases the electron fraction (see
Fig.1).
The electron fraction changes significantly as
a function of density for T > 1010K (in Fig.2,
T = 1011K). At low densities, the torus consists
of free neutrons and protons and Ye is close to 0.5
(see also Eq. 7). As density increases, Ye decreases
to satisfy the beta-equilibrium among the free n-p-
e gas. Above some density (when the temperature
is high enough, e.g. for Fig. 2, ρHe ≈ 1013g cm−3)
helium starts forming. Therefore Ye has a kink
and stars rising steeply, asymptotically approach-
ing 0.5 as the torus consists of plenty of ionized
helium and some electrons to keep charge neutral-
3.2. The steady-state disc structure
In Figures 3 and 4 we show the profiles of den-
sity and temperature in the stationary accretion
disk model for three accretion rates: 1 M⊙/s,
10 M⊙/s and 12 M⊙/s. In general, the temper-
ature and density profiles both increase inward.
However, for Ṁ = 12M⊙/s, a distinct branch of
solutions is reached, which appears different than
the so-called “NDAF” branch (see Kohri & Mi-
neshige 2002). The density and temperature pro-
files for this high accretion rate differ also from
what was found in previous work (Di Matteo et
al. 2002; Janiuk et al. 2004). Due to a more
detailed equation of state, in which we allow for
Fig. 1.— The dependence of the electron frac-
tion (bottom), pressure (middle bottom), en-
tropy (middle upper) neutrino cooling rate (up-
per panel) on temperature, for the constant den-
sity ρ = 1012 g/cm3. The accretion rate is Ṁ =
1 M⊙/s. The pressure and neutrino cooling are in
cgs units and the entropy is in units of kB cm
Fig. 2.— The dependence of the electron frac-
tion (bottom), pressure (middle bottom), entropy
(middle upper) and neutrino cooling rate (upper
panel) on density, for the constant temperature
T = 1011 K. The accretion rate is Ṁ = 1 M⊙/s.
The pressure and neutrino cooling are in cgs units
and the entropy is in units of kBcm
a partial degeneracy of nucleons and electrons as
well as neutrino trapping, our solutions reach den-
sities as high as 1012 g/cm3 in the innermost radii
of the disc. The temperature in this inner disc part
is in the range 4×1010−1.25×1011 K, depending
on the accretion rate. For the hottest disk model,
a local peak in the density forms around 7 − 8RS,
while below that radius the density decreases. Be-
tween ∼ 3.5 and 7 RS, the plasma becomes much
hotter and less dense than outside of this region.
This means that the macroscopic state of the sys-
tem is different here due to an abrupt change in
the heat capacity. In order to check what is the
reason for this transition, we investigate the pres-
sure distribution in the disk.
The profile of the pressure is shown in Fig. 5.
The dominant term in the total pressure is due
to the nucleons, while the radiation pressure (in-
cluding electron-positron pairs) is always several
orders of magnitude smaller. The neutrino pres-
sure is large in the inner disc, once it gets opti-
cally thick to neutrinos (i.e. for Ṁ ≥ 10M⊙/s).
A significant contribution to the pressure is due
to helium at densities high enough for helium to
form, albeit at temperatures low enough such that
its nuclei are not fully disintegrated. For the
largest accretion rate shown, in the region of the
temperature excess and inverse density gradient
(3.5 − 7RS), the total pressure distribution flat-
tens. The helium pressure is now vanishingly small
due to the complete photodisintegration, and the
nuclear pressure is slightly decreased due to the
composition change: smaller number density of
neutrons and larger number density of protons.
The substantial contribution to the pressure is now
given by the neutrinos (large optical depths; see
below) and radiation pressure (increased number
of electron-positron pairs). From the comparison
of Figures 3, 4 and the bottom panel of Figure
5, it can be seen that the total pressure becomes
locally correlated with temperature and anticor-
related with density, thus consituting an unstable
phase.
In Figure 6 we show the neutrino optical depths
due to scattering and absorption. The total opti-
cal depth in the outer disc is typically dominated
by scattering processes, while in the inner disc ab-
sorption processes take over for very high accre-
tion rates. For Ṁ = 1M⊙/s only the very in-
ner disk radii have optical depth close to 1. For
Fig. 3.— The baryon density as a function of
the disc radius, calculated in the stationary so-
lution. The accretion rate is Ṁ = 1 M⊙/s (solid
line), Ṁ = 10 M⊙/s (long dashed line) and Ṁ =
12 M⊙/s (short dashed line) .
Fig. 4.— The temperature as a function of the disc
radius, calculated in the stationary solution. The
accretion rate is Ṁ = 1 M⊙/s (solid line), Ṁ =
10 M⊙/s (long dashed line) and Ṁ = 12 M⊙/s
(short dashed line) .
Fig. 5.— The pressure components as a func-
tion of the disc radius, calculated in the station-
ary solution for the three accretion rate values:
Ṁ = 1 M⊙/s (upper panel), Ṁ = 10 M⊙/s (mid-
dle panel) and Ṁ = 12 M⊙/s (bottom panel)
. The total pressure is marked by the solid
line, and its components are: nuclear (gas) pres-
sure (long dashed line), radiation pressure (short
dashed line), helium pressure (dotted line) and
neutrino pressure (dot-dashed line).
Fig. 6.— The neutrino optical depths due to scat-
tering (τs, solid line) and absorption (τa,e for elec-
tron neutrinos, long dashed line, and τa,µ for muon
neutrinos, short dashed line) as a function of ra-
dius for Ṁ = 1M⊙/s (upper panel), Ṁ = 10M⊙/s
(middle panel) and for Ṁ = 12M⊙/s (bottom
panel). The sum of the three quantities is the
total optical depth (τtot, dotted line).
Ṁ = 12M⊙/s, in the radial strip of ∼ 3.5 − 7RS
the disk is optically thick with absorptive optical
depth for electron neutrinos exceeding the scatter-
ing term and reaching values of the order of 100.
3.2.1. Composition and Chemical potentials
In Figure 7 we show the distribution of the
reduced chemical potentials of protons, electrons
and neutrons throughout the disc. Reduced elec-
tron chemical potentials much larger than unity
(indicating strong electron degeneracy) are found
in the inner disc parts for Ṁ = 10M⊙/s and
Ṁ = 12M⊙/s, whereas for 1 M⊙/s electrons are
only slightly degenerate. For the highest accretion
rate, the maximum degeneracies correspond to the
radius of the local peak in the density (cf. Fig. 3)
and the excess of helium number density (cf. Fig.
8). Below this radius, the species become non-
degenerate again, contributing to the increase of
the electron fraction (cf. Fig. 9). In Figure 8 we
plot the mass fraction of free nucleons as a function
of radius for Ṁ = 12M⊙/s, 10 and 1M⊙/s. As the
Figure shows, in the outer regions, Xnuc increases
as the radius decreases, while the temperature and
density increase (Fig. 3 and Fig. 4). Consistent
with the behavior of Ye (Fig.2), Xnuc subsequently
turns around (decreases) at radii where the den-
sity is high enough for significant helium forma-
tion. This trend is reversed sharply for highest
accretion rates, when the temperature in the disk
is high enough (Fig.4) for helium to be fully dis-
sociated. In consequence, the number density of
alpha particles increases at ∼ 7−12RS and sharply
decreases at lower radii. A similar, but far less pro-
nounced fluctuation in Xnuc is seen at smaller radii
for the case of Ṁ = 10M⊙/s. For smallest accre-
tion rate, 1M⊙/s, there is no helium throughout
the disk.
In Figure 9 we show the radial distribution of
the electron fraction throughout the disc for 1,
10 and 12 M⊙/s. For the case of 1 M⊙/s (solid
line), the electron fraction decreases inward in the
disc as the electrons are captured by protons (in
neutronization reactions). Once the electrons be-
come non-degenerate, positrons appear, and the
positron capture by neutrons again increases the
electron fraction. For the hotter plasma (accre-
tion rate of 10 M⊙/s, dashed line), consistently
with the behavior discussed for Xnuc, helium nu-
clei form as the density becomes high enough be-
Fig. 7.— The chemical potentials of neutrons
(ηn, solid line), electrons (ηe, dashed line) and
protons (ηp, dotted line) as a function of the
disc radius, calculated in the stationary solution.
The accretion rate is Ṁ = 1 M⊙/s (triangles),
Ṁ = 10 M⊙/s (squares) and Ṁ = 12 M⊙/s (cir-
cles).
Fig. 8.— The mass fraction of free nucleons
as a function of radius for Ṁ = 1M⊙/s (solid
line), Ṁ = 10M⊙/s (long dashed line) and Ṁ =
12 M⊙/s (short dashed line).
low ∼ 20RS and Ye increases. For the accretion
rate of 12 M⊙/s there is the sharp decrease in Ye,
at ∼ 7−8RS, due to the sudden dissociation of he-
lium. As helium is fully photo-dissociated, there is
an almost equal number of neutrons and protons
due to the balance of the electron and positron
capture. This implies an electron fraction of 0.5.
At the innermost radius, the temperature and den-
sity drop due to the boundary condition, which
affects the behaviour of both Ye and Xnuc.
3.2.2. Cooling and heating rates
In Figure 10 we plot the rates of viscous heat-
ing, advection and cooling due to neutrino emis-
sion and photo-dissociation in the stationary disc.
The accretion rates are Ṁ = 1 M⊙/s (upper
panel), Ṁ = 10 M⊙/s (middle panel), and Ṁ =
12 M⊙/s (lower panel). For the highest accretion
rates, in the innermost disc the neutrino cooling
rate decreases substantially with respect to the
cooling by photodissociacion. This is because the
neutrinos are trapped in the disc due to a large
opacity The smaller the accretion rate, the less
important is the neutrino trapping effect. This
implies that for an accretion rate of ≤ 10 M⊙/s
neutrinos can escape from the innermost disc.
The advective term is a couple orders of magni-
tude smaller than the other terms. The photodis-
sociacion term is negligible for an accretion rate of
1 M⊙/s, since there is no helium in the whole disc,
and Qphoto is equal to zero by definition. For an
accretion rate of 10 M⊙/s there is very little he-
lium down to about 15-20 RS, and therefore Qphoto
is much smaller than other terms. For the accre-
tion rate of 12 M⊙/s , down to 6 − 10RS in the
region of the disc of high density and maximum
degeneracy, helium nuclei form. The nucleosyn-
thesis of alpha particles leads to the plasma heat-
ing instead of cooling, and therefore the relevant
term in the energy balance has a negative value.
Outward, above ∼ 10RS, there is some fraction
of helium which can be photo-dissociated, so the
cooling term due to this reaction is also important
in the total energy balance. In the inner region
helium is fully dissociated and Qphoto is equal to
zero, increasing again only near the inner bound-
ary due to the local density increase and decrease
of temperature.
Fig. 9.— The electron fraction as a function
of the disc radius, calculated in the stationary
solution. The accretion rate is Ṁ = 1 M⊙/s
(solid line), Ṁ = 10 M⊙/s (long dashed line) and
Ṁ = 12 M⊙/s (short dashed line).
Fig. 10.— The heating and cooling rates due
to photodissociation and neutrino emission (solid
lines) as a function of radius for Ṁ = 1M⊙/s (up-
per panel), Ṁ = 10M⊙/s (middle panel) and for
Ṁ = 12M⊙/s (bottom panel). The other terms
are: cooling rate due to advection (long dashed
line) and viscous heating rate (short dashed line).
3.3. Stability analysis: instabilities at
high-Ṁ
The disc is thermally unstable if d logQ+/d logT >
d logQ−/d logT . Then any small increase (de-
crease) in temperature leads to a heating rate
which is more (less) than the cooling rate, and
as a consequence a further increase (decrease) of
the temperature. The viscous instability, which
appears when ∂Ṁ
∂Σ |Q+=Q− < 0, manifests itself
in a faster (slower) evolution of an underdense
(overdense) region. The instabilities can be con-
veniently located in the surface density - temper-
ature diagrams, in which the branch of thermal
equilibrium solutions with a negative slope is not
only unstable to the perturbations in the surface
density, but it is also thermally unstable.
In Figure 11 we show such stability curves for
several radii in the disc. The criterion for a vis-
cously stable disc is generally satisfied through-
out the whole disc for Ṁ ≤ 10M⊙/s. How-
ever, for larger accretion rates, there are unstable
branches at the smallest radii. For Ṁ = 10M⊙/s,
the disc becomes unstable below 5 RS, while for
Ṁ = 12M⊙/s the instability strip is up to ∼ 7RS.
Here helium is almost completely photodisinte-
grated while the electrons and protons become
non-degenerate again. For this high accretion rate,
the electron fraction rises inward in the disc. Un-
der these conditions, the energy balance is affected
leading to the thermal and viscous instability, as
demonstrated by the stability curves. This insta-
bility will be discussed in more detail in Section
3.4. Time dependent solutions
In this Section, we discuss how the tempera-
ture, density, electron fraction and disk luminosity
evolve with time. In Figures 12 and 13 we show the
time evolution of density and temperature, when
the initial accretion rate is 1 M⊙/s. These quan-
tities exponentially decrease with time:
ρ = ρ0(r) exp(−at) , (47)
T = T0(r) exp(−bt) (48)
where a ≈ 1.9 and b ≈ 0.085. The normalization
of these relations depends on the radius, and for
example for r = 6RS it is ρ0 = 2.2 × 1011 and
T0 = 3.5× 1010. The exponential behaviour arises
from the nature of energy equation (45).
In Figure 14 we show the electron fraction as
a function of time for several exemplary radial lo-
cations in the disc, for the disc evolving from a
starting accretion rate of 1M⊙/s. The fraction Ye
is smaller in the inner disc radii, while outward,
the electron fraction is over half an order of mag-
nitude higher. Altogether, during the evolution
of the system, the electron fraction constantly in-
creases with time throughout the disc.
The time-dependent neutrino luminosity of the
disc is given by:
Lν(t) =
∫ Rmax
Q−ν (t)2πrdr (49)
where Q−ν is given by Equation (30).
In Figure 15 we show an example of such a
lightcurve, for our standard model parameters
(M = 3M⊙, α = 0.1), and Rmax = 50RS. The
starting accretion rate was Ṁstart = 1 M⊙/s. At
this accretion rate neutrinos can already escape
from the accretion disc at the beginning of the
evolution. For higher initial accretion rates, e.g
10 − 12M⊙/s, neutrinos are trapped in the inner-
most disc, and, as a consequence, the neutrino
luminosity is lower at the initial stages of disc
evolution until the accretion rate drops to about
∼ 1M⊙/s. This result is qualitatively similar, al-
beit it differs quantitatively, from what was ob-
tained in Di Matteo et al. (2002) and Janiuk et
al. (2004): in those calculations neutrino trap-
ping was far more substantial even for a ’mod-
erate’ accretion rate of Ṁ & 1M⊙/s. The differ-
ence arises from the fact that here we calculate the
neutrino opacities using the β reaction efficiencies,
self-consistently with the equation of state.
For an accretion rate of 1M⊙/s, the solution
does not reach the viscously unstable branch. Ini-
tially, the disc contains almost no alpha particles
(cf. Fig. 8), which appear later on during the evo-
lution and cooling of the plasma. The dynamical
balance between the photodisintegration of helium
and nucleosynthesis leads to an additional non-
zero cooling/heating term in the energy equation
and to only small amplitude flickering at the early
stages of time-evolution.
The situation is much more dramatic when the
starting accretion rate is 12M⊙/s. In this case a
Fig. 11.— The stability curves on the accretion
rate vs. surface density plane, for several chosen
radii in the disc: 3.39RS (solid line), 3.81RS (dot-
ted line), 4.25RS (short dashed line), 5.19RS (long
dashed line) and 8.60RS (dot-dashed line).
Fig. 12.— The density as a function of time, for
several chosen disc radii: 4.01, 6.04, 10.13, 20.7,
35.02, and 45.19 RS. The initial accretion rate is
Ṁ = 1 M⊙/s.
Fig. 13.— The temperature as a function of time,
for several chosen disc radii: 4.01, 6.04, 10.13,
20.7, 35.02, and 45.19 RS. The initial accretion
rate is Ṁ = 1 M⊙/s.
Fig. 14.— The electron fraction as a function
of time, for several chosen disc radii: 4.01, 6.04,
10.13, 20.7, 35.02, and 45.19 RS. The initial ac-
cretion rate is Ṁ = 1 M⊙/s.
large disc strip is viscously and thermally unstable
and the most violent instability takes place around
and below 7 − 12RS.
In Figure 16 we show the behavior of the lo-
cal accretion rate in the unstable disc, at sev-
eral chosen locations within the instability strip.
Near ∼ 12RS, the accretion rate varies due to
the large and rapidly changing photodisintegra-
tion term (locally, it can become larger than the
neutrino cooling rate).
This radius corresponds to the largest local
value of the density of helium (cf. Figure 8 show-
ing its starting model distribution), which is then
being photodissociated. The photodissociacion
process is the cause of the local rapid accretion
rate changes. Then, inside from this highly vari-
able strip, the accretion rate grows too fast to pre-
serve the disc structure. This kind of behaviour
occurs in the locally hotter and less dense region
visible in the starting configuration e.g. in Figures
3 and 4, between 3.5 and 7 RS. In this region the
helium is already totally photodissociated. Due
to the growing accretion rate all the material is
rapidly accreted onto the black hole and the in-
nermost strip of the disc empties.
After the inner strip is destroyed, the outer
parts can still accrete onto the center. As they
approach the black hole, their temperature and
density grow and the above situation can repeat
several times, until the whole disc is completely
broken into rings and destroyed. These later in-
jections of energy, with timescales dictated by the
viscous timescale of each ring, can produce en-
ergy flares following the main GRB activity. Our
results therefore provide another physical mech-
anism1 for the flare model recently proposed by
Perna et al. (2006).
In Figure 17 we show the neutrino lightcurve of
the unstable disc. The instabilities due to photo-
disintegration are reflected in oscillations of vari-
able amplitude and millisecond timescale. This is
of a particular interest if the neutrino annihila-
tion provides the energy input for GRBs, however
it should be pointed out that the oscillations ap-
pearing in the presented lightcurve have a much
1In addition to the gravitational instability in the outer parts
of the disk, which was hinted by the calculations of Di
Matteo et al. (2002) and confirmed by those of Chen &
Beloborodov (2006).
Fig. 15.— The neutrino lightcurve, integrated
over the disc surface. The initial accretion rate
is Ṁ = 1 M⊙/s.
Fig. 16.— The local accretion rate, as a function of
time, for several chosen radial locations in the disc:
6.86, 7.73 and 12.85 RS. The starting accretion
rate is Ṁ = 12 M⊙/s.
smaller amplitude than the observed gamma ray
variability.
4. Discussion
4.1. The unstable neutrino-opaque disc
In our calculations we have shown that, for
large accretion rates, the accreting torus becomes
viscously and thermally unstable. We now discuss
the physical origin of the instability.
The unstable branch appears both in the
steady-state solutions and in the subsequent time-
dependent evolutions. In the steady-state case,
for a chosen value of a constant accretion rate,
this can be seen for instance by plotting the ra-
dial profiles of density and temperature (cf. Figs.
3 and 4, where the distinct branch is found for
the innermost radii), as well as by looking at the
stability curves for a range of accretion rates at
a chosen disc radius (cf. Fig. 11, where the un-
stable inner disc radii exhibit a negative slope in
the curve). In the time-dependent simulations,
the unstable behavior is manifested by the highly
variable accretion rate in certain strips of the disk
and by the subsequent breaking of the disk inside
from these variable strips (cf. Fig. 16). The in-
stability arises from the fact that the accretion
rate rises locally too fast to prevent the disc strip
from emptying, as the material is supplied from
outer strips at much slower rate than it is accreted
inwards. The disc evolves unstably on a viscous
timescale, τvisc = 1/(αΩ) × (r/H)2; for the radii
shown in Figure 16, it is τvisc ≈ 0.05 s (note that
the disc is rather thick, r/H ∼ 2.5, and therefore
the viscous and thermal timescales are close to
each other). Theoretically, in order to find again
a stable solution, the disc would have to increase
the local accretion rate up to about several tens
of M⊙/s within one viscous timescale. However,
this may not be possible if there is not enough
material in the system to support much higher
accretion rates during such violent oscillations.
Therefore the system is unable to be stabilized
and gets broken after a fraction of τvisc . In addi-
tion, the dynamical instability is the source of the
flickering of the local accretion rate at the edge of
the unstable strip.
Let us now discuss in more detail the physical
reason driving this instability. In the inner part of
the disc (below r ∼ 10RS for Ṁ = 12M⊙/s) there
Fig. 17.— The neutrino luminosity, as a func-
tion of time The starting accretion rate is Ṁ =
12 M⊙/s.
Fig. 18.— The total pressure, as a function of
time, for the chosen radial locations in the disc:
7.73 and 6.86 RS. The starting accretion rate is
Ṁ = 12 M⊙/s.
are two important processes, both of which are
incorporated in our equation of state: photodisin-
tegration of helium and neutrino trapping. As it
was already mentioned in Sec. 3.2.1 and can be
seen from Fig. 8, below ∼ 7 − 8RS helium in this
disc is completely photodisintegrated. This part
of the disc is also opaque to neutrinos, as we show
in Fig. 6.
These two mechanisms competitively influence
the electron fraction in the disc (cf. Sec. 3.2, Fig.
6 and Fig. 9). Well outside the unstable strip,
above ∼ 20RS, the electron fraction smoothly de-
creases inwards as positrons appear because of the
neutronization process. Then the scattering opti-
cal depth for neutrinos becomes τs > 1, and the
electron fraction increases again. After photodis-
integration, the electron fraction decreases signif-
icantly from almost 0.3 to much less than 0.1 due
to electron capture. But again, when the disc
becomes optically thick to absorption of electron
neutrinos, the electron fraction gets higher and ap-
proaches almost 0.5.
The total pressure of sub-nuclear matter (cf.
Fig 5) is mainly contributed by electrons, and
therefore it is influenced by the changes in the
electron fraction. In the narrow range of radii
(6.8 < r/RS < 7.8), the pressure decreases due
to photodisintegration. The sudden decrease of
the pressure might drive the dynamical instabil-
ity. (This picture is somewhat similar to that of
the iron core collapse in the core collapse super-
nova explosions: electron capture consumes most
of the electrons and makes the EOS softer, conse-
quently, it triggers the collapse of the iron core.)
However, the transition from neutrino transpar-
ent to opaque disc, and the increase of the elec-
tron fraction due to the beta equilibrium (see also
Yuan & Heyl 2005), are the reason for a steeper
increase of the total pressure of the system.
The same effect can also be observed in the
time-dependent plot (Fig. 18), in which we show
the pressure changes in the characteristic radii
of the unstable part of the disc (cf. Fig. 16).
The pressure decreases with time up to a radius
R = 6.8RS, since the temperature and density
gradually drop, as well as the neutrino opacities,
so the electron fraction gets smaller. Then, in a
strip between ∼ 6.8 and 12 RS, the pressure rises
with time: in fact, when α particles appear, the
electron fraction rises and matter locally piles up,
thus increasing the pressure.
At the border of these radii the disc breaks up,
when the thermal-viscous instability induces an
avalanche-like increase of the local accretion rate
below ∼ 8RS. This happens because the increase
in the pressure causes an excess in the local energy
dissipation rate and the disc heats up, while at ad-
jacent radii the pressure decreases and heating is
insufficient. The system tries to compensate these
temperature gradients by decreasing/ increasing
the temperature in the outer/inner radius, respec-
tively. But since in the unstable mode of the ther-
mal balance this causes a further increase/decrease
of density, the pressure drops further and the disc
heats up in the outer radius, while cools down in
the inner one. As long as it cannot find any sta-
ble track of evolution, the emptying of the inner
strip continues and finally the whole material is
accreted towards the black hole or blown out.
The radial extent of the unstable part depends
on the initial accretion rate and in our model for
Ṁ = 12M⊙/s it is up to ∼ 8RS, for Ṁ = 10M⊙/s
it is up to ∼ 5RS, while for Ṁ = 1M⊙/s it is below
∼ 3.5RS. In the latter case, since the inner radius
is located at ∼ 3RS, the instability hardly affects
the disc. The extension of the instability strip de-
pends also on the mass of the accreting compact
object, and since for lower mass black holes the ac-
cretion disc is generally denser, it reaches a density
∼ 1012 g/cm3 around ∼ 15RS.
Above 25-30 RS, where the plasma is already
optically thin and the evolution is stable, both
the pressure and the accretion rate smoothly drop
with time. The dominant source of cooling of the
disc in this region is the neutrino emission (ad-
vective cooling decreases as the disc transits from
neutrino opaque to transparent). The photodis-
integration term (if non-zero), is usually by 1-2
orders of magnitudes smaller than neutrino cool-
ing, and in the inner disc, up to about 6.8 RS,
there are no helium nuclei and the photodisinte-
gration term is negligible, while at 7.5 RS it has
a value of about Qphoto ∼ 1038 erg/s/cm2 with
rapid fluctuations. These fluctuations induce the
local accretion rate flickering (cf. Fig. 16), on
a timescale and amplitude much smaller than for
the viscous instability.
4.2. Comparison with previous work
The neutrino dominated accretion flow has al-
ready been studied in a number of papers, includ-
ing both 1-D models and multi-D simulations. The
steady-state 1-D models (e.g Popham et al. 1999;
Kohri & Mineshige 2002) assumed the disk opti-
cally thin to neutrinos, and neglected photodis-
integration cooling. Di Matteo et al. (2002) took
these two effects into account, and showed that the
trapped neutrinos dominate the pressure in the
inner region of the hyperaccreting disc, however
their equation of state did not include the numeri-
cal calculation of chemical equilibrium and did not
incorporate the opacities directly in the EOS itera-
tions (see also the time-dependent model of Janiuk
et al. 2004). Kohri, Narayan & Piran (2005) con-
sidered the neutrino opaque disk and the equilib-
rium between neutrons and protons and calculated
the number densities of species by numerically in-
tegrating their distribution functions. However,
these authors calculated the gas pressure from the
ideal gas approximation, and neglected the contri-
bution of helium to the pressure. In all of these
papers the disc occurred to be stable against any
kind of instability.
On the other hand, in their recent work, Chen &
Beloborodov (2006) find that the outskirts of the
disk are gravitationally unstable. The approach
used by these authors provides a detailed treat-
ment of the microphysics which is very similar to
ours; however, some differences between our work
and theirs must be crucial to the development of
viscous and thermal instabilities. One difference
deals with the approximation made for the treat-
ment of transition region between the neutrino-
opaque and the transparent matter. In our work,
we adopted a gray body model, i.e., we introduced
the b factor to describe the distribution function
(c.f. Eq. 10). This assumption is consistent with
the two fluids approximation we have made, which
has recently been studied numerically by Sawyer
(2003), and shown to be appropriate for the con-
ditions of these disks. On the other hand, Chen
& Beloborodov (2006) smoothly connect the opti-
cally thin and thick regimes by means of interpo-
lation. A further difference lies in the description
of the mass fraction of free nucleons. In this work
we use an expression for Xnuc developed by Qian
et al (1996), while in Chen & Beloborodov (2006)
Xnuc is a function of Ye, which couples the nucle-
osynthesis to the electron fraction.
In our calculations we reach the range of densi-
ties and temperatures where the nucleons start to
become partially degenerate. This is accompanied
by the neutrinos being more and more trapped in
the gas and helium being destroyed by photodisso-
ciation. As a result of these calculations, we found
an additional, unstable branch of solutions for the
disc thermal balance.
This supports the recent results of 2-D simu-
lations by Lee, Ramirez-Ruiz & Page (2005), who
found the disc opaque to neutrinos to be thermally
unstable. Their simulations showed that large cir-
culations develop in the accretion flow. Setiawan,
Ruffert and Janka (2005) found small fluctuations
of the accretion rate and neutrino luminosity on
the dynamical timescale, after the 10-20 msec of
relaxation period (note, that in our calculations we
start from the steady-state disc model at a given
accretion rate, thus having no need for a relaxation
to the quasi-steady configuration). The equation
of state used in their work (see also Janka et al.
1999) is based on the work of Lattimer & Swesty
(1991). Given the electron fraction, this EOS as-
sumes the condition of nuclear statistical equilib-
rium without neutrino trapping, but the evolution
of the electron fraction is affected by the asym-
metric neutrino emission from the hot and dense
matter, which is called ‘neutrino leakage scheme’.
The neutrino leakage scheme focuses on the ef-
fects of the neutrino trapping on the net neutrino
emissivities, not on the nuclear statistical equilib-
rium. The equation of state used in the work of
Rosswog et al. (2004) is temperature and compo-
sition dependent, based on the relativistic mean
field theory (Shen et al. 1998a,b), and the neu-
trino cooling is accounted for by the multiflavor
scheme (Rosswog & Liebendoerfer 2003).
In our work, we use an equation of state based
on the β equilibrium, including the contribution
from the trapped neutrino, and neutrino trapping
effects are accounted for by the appropriate opac-
ities. It should be emphasized that most previous
multi-D simulations neglected the effects of neu-
trino trapping on the β equilibrium, as well as the
contribution of the trapped neutrinos to the ther-
modynamical properties of the dense matter. An-
other difference between our treatment of the EOS
and the previous numerical simulations is that we
include the cooling of the photodisintegration of
helium. Even though the original EOS of Lattimer
& Swesty (1991) can provide detailed information
about the composition of the dense matter, this
information was not considered in order to keep
the table of the EOS as small as possible (see e.g.
Ruffert et al. 1996) just for numerical reasons. In
this way, the disintegration cooling had not been
investigated without the information on the com-
position. Our results indicate that photodisinte-
gration significantly affects the energy balance.
4.3. Limitations of our model
We find the thermal-viscous instability to be an
intrinsic property of the disc for extremely large
torus densities (about 1012g cm−3) and high ac-
cretion rates (Ṁ ≥ 10M⊙/s). This is seen both
in the steady-state results (radial profiles of den-
sity and temperature) and in the subsequent time
evolution.
Thermal and viscous instabilities have been
studied in the case of standard accretion discs
around compact objects (Lightman & Eardley
1974; Pringle 1977; Shakura & Sunyaev 1976).
Two main physical processes that lead to disc
instabilities were invoked to explain the time-
dependent behavior of various objects: partial
ionization of hydrogen in the discs of Dwarf No-
vae (e.g. Meyer & Meyer-Hofmeister 1981; Smak
1984) and domination of radiation pressure in the
X-ray transients (e.g Taam & Lin 1984). Such in-
stabilities do not have to lead to a total disc break-
down, but rather to a limit-cycle behavior, if only
an additional (i.e. upper) stable branch of solu-
tions can be found. This might be a hot state with
a temperature above 104 K, or a slim disc, domi-
nated by advection (Abramowicz et al. 1988). In
our 1-D calculation the disc in the GRB central en-
gine is not stabilized but rather breaks down into
rings, as no stable solutions are reached (possibly,
for even higher accretion rates again a stable part
near the black hole could be formed - but these
extremely high accretion rates would not be pro-
duced by any compact merger scenario). There-
fore, instead of a limit-cycle activity, what we
find here are several dramatic accretion episodes
on the viscous timescale. The remaining parts of
the torus will subsequently accrete and, while ap-
proaching the central black hole, will get hotter
and denser, breaking at ∼ 7RS.
Of course, it would be interesting to study
whether such a violent instability would occur also
in the 2D or 3D simulations. This is indeed likely
to be the case, since as the multi-dimensional sim-
ulations of accretion discs show, the instabilities
derived first in 1D are still present in the hydro-
dynamical simulations of flows with non-Keplerian
velocity fields (e.g. Agol et al. 2001; Turner
2004; Ohsuga 2006). Possibly, the instability re-
gion would be located at other (larger) radii if
the calculations included the vertical structure of
the disc: this is dependent on temperature and
density, which above the meridional plane may be
larger than the mean value considered in the ver-
tically averaged model.
We need to note that our 1D calculations do
not take into account the possible effects of non-
radial velocity components in the fluid. For exam-
ple, the inverse composition gradient that leads
to the disk instability, might be stabilized by ro-
tation (e.g. Begelman & Meier 1982; Quataert &
Gruzinov 2000). In the 2-D simulation of Lee et al
(2005) the neutrino opaque disk exhibits circula-
tions in the r-z direction. Such meridional circula-
tions are known to be present in the Keplerian ac-
cretion disks (e.g. Siemiginowska 1988), however
it is unclear if they could always provide a stabiliz-
ing mechanism for the thermal-viscous instability.
Possibly, if the nonradial motions of the flow pro-
vided a stronger stabilizing effect, the disk would
exhibit oscillations in the viscous timescale, with-
out breaking, similarly to the outbursts of Dwarf
Nova disks.
The assumption of the β equilibrium (justified,
as the mixture of protons, electrons, neutrons and
positrons is able to achieve the equilibrium con-
ditions) might also have an effect on this result,
as the equilibrium conditions reduce the heating
and entropy in the gas. In fact, the β equilibrium
condition which is satisfied in the innermost part
of a hyperaccreting disc that is optically thick to
neutrinos, is µn = µp +µe. Once the disc becomes
transparent in its outer part, this condition is no
longer valid. Analytically, it has been derived by
Yuan (2005) that the condition for β equilibrium
in this case is µn = µp + 2µe.
4.4. Observational consequences
Our findings might be relevant for interpret-
ing some recent observations. The flickering due
to the photodisintegration of alpha particles may
lead to a variable energy output on small (millisec-
ond) timescales. The consequence of this may be
variability in the gamma ray luminosity, although
the changes in the local accretion rate may be
spread by viscous effects (in the lightcurve Lν(t)
integrated over the whole surface of the disc, the
millisecond variability is somewhat smeared, and
the amplitudes are not very large). Therefore the
mass accreted by the black hole may not be vary-
ing substantially, while some irregularity in the
overall outflow could help produce internal shocks.
The thermal-viscous instability, if accompa-
nied by the disc breaking, may lead to the sev-
eral episodic accretion events and several re-
brightenings of the central engine on longer
timescales, possibly detected in the later stages
of the evolution. A similar kind of a long-term
activity is possible also if the disk was not com-
pletely broken, but exhibited some large accretion
rate fluctuations on the viscous timescale.
We thank Bożena Czerny, Pawe l Haensel and
Daniel Proga for helpful discussions. We also
thank the anonymous referee for detailed reports
which helped us to improve our model and its
presentation. This work was supported in part
by grant 1P03D 00829 of the Polish State Com-
mittee for Scientific Research and by NASA un-
der grant NNG06GA80G. Y.-F. Y. is partially
supported by Program for New Century Excel-
lent Talents in University, and the National Natu-
ral Science Foundation (10233030,10573016). RP
acknowledges support from NASA under grant
NNG05GH55G, and from the NSF under grant
AST 0507571.
A. Appendix
The neutrino absorption and production rates in the beta processes for all participating particles at
arbitrary degeneracy have been obtained in the previous works (Reddy, Prakash & Lattimer 1998; Yuan
2005). In the subnuclear dense matter with high temperatures, the nucleons are generally nondegenerate,
therefore, the transition reaction rates from neutrons to protons and from protons to neutrons can be
simplified as follows:
Γp+e−→n+νe =
|M |2
dEeEepe(Ee −Q)2fe(1 − befνe), (A1)
Γp+e−←n+νe =
|M |2
dEeEepe(Ee −Q)2(1 − fe)befνe , (A2)
Γn+e+→p+ν̄e =
|M |2
dEeEepe(Ee + Q)
2fe+(1 − befν̄e), (A3)
Γn+e+←p+ν̄e =
|M |2
dEeEepe(Ee + Q)
2(1 − fe+)befν̄e , (A4)
Γn→p+e−+ν̄e =
|M |2
dEeEepe(Q− Ee)2(1 − fe)(1 − befν̄e), (A5)
Γn←p+e−+ν̄e =
|M |2
dEeEepe(Q− Ee)2febefν̄e . (A6)
Here Q = (mn −mp)c2, |M |2 is the averaged transition rate which depends on the initial and final states
of all participating particles, for nonrelativistic noninteracting nucleons, |M |2 = G2F cos2 θC(1 + 3g2A), here
GF ≃ 1.436×10−49 erg cm3 is the Fermi weak interaction constant, θC (sin θC = 0.231) is the Cabibbo angle,
and gA = 1.26 is the axial-vector coupling constant. fe,νe are the distribution functions for electrons and
neutrinos, respectively. The “chemical potential” of neutrinos is generally assumed to be zero. The factor
be reflects the percentage of the partially trapped neutrinos. When neutrinos completely trapped, be = 1.
The corresponding neutrino emissivities for the URCA reactions are given by:
qp+e−→n+νe =
|M |2
dEeEepe(Ee −Q)3fe(1 − befνe), (A7)
qn+e+→p+ν̄e =
|M |2
dEeEepe(Ee + Q)
3fe+(1 − befν̄e), (A8)
qn→p+e−+ν̄e =
|M |2
dEeEepe(Q− Ee)3(1 − fe)(1 − befν̄e). (A9)
The emissivities due to the electron-positron pair annihilation, following the notation of Yakovlev et al
(2001), is written as:
qe−+e+→νi+ν̄i =
C2+νi [8(Φ1U2 + Φ2U1) − 2(Φ−1U2 + Φ2U−1) + 7(Φ0U1 + Φ1U0)
+ 5(Φ0U−1 + Φ−1U0)] + 9C
−νi [Φ0(U1 + U−1) + (Φ−1 + Φ1)U0]
, (A10)
where
= 1.023 × 1023 erg cm−3 s−1, (A11)
C+νi = C
+ C2Ai and C−νi = C
− C2Ai , here CVi and CAi are the vector and axial-vector constants for
neutrinos (CV e = 2 sin
2 θC + 0.5, CAe = 0.5, CV µ = CV τ = 2 sin
2 θC − 0.5 and CAµ = CAτ = −0.5). The
dimensionless functions Uk and Φk (k= −1, 0, 1, 2) in the above equation can be expressed in terms of the
Fermi-Dirac functions:
U−1 =
β3/2F1/2(ηe, βe) (A12)
β3/2[F1/2(ηe, βe) + βeF3/2(ηe, βe)] (A13)
β3/2[F1/2(ηe, βe) + 2βeF3/2(ηe, βe) + β
eF5/2(ηe, βe)] (A14)
β3/2[F1/2(ηe, βe) + 3βeF3/2(ηe, βe) + 3β
eF5/2(ηe, βe) + β
eF7/2(ηe, βe)]. (A15)
Replacing ηe with ηe+ in Uk, we get the corresponding expressions for Φk.
REFERENCES
Aparicio J., 1998, ApJS, 117, 627
Abramowicz M.A., Czerny, B., Lasota J.P.,
Szuszkiewicz E., 1988, ApJ, 332, 646
Abramowicz M.A., Kato S., 1989, ApJ, 336, 304
Agol E., Krolik J., Turner N.J., Stone J.M., 2001,
ApJ, 558, 543
Begelman M.C., Meier D.L., 1982, ApJ, 253, 873
Berger E., et al., 2005, submitted to Nature
(astro-ph/0508115)
Chen, W.-X. & Beloborodov, A. preprint
astro-ph/0607145
Di Matteo T., Perna R., Narayan R., 2002, ApJ,
579, 706
Eichler D., Livio M., Piran T., Schramm D. N.,
1989, Nature, 340, 126
Falcone, A. et al. 2006, ApJ, 641, 1010
Fox D.B., et al., 2005, Nature, 437, 845
Fujimoto S., Kotake K., Yamada S., Hashimoto
M., Sato K., 2006, ApJ, 644, 1040
Gehrels N., et al., 2005, Nature, 437, 851
Gu, W.-M., Liu, T., & Lu, J.-F. 2006, ApJ, 643,
Hjorth, J. et al. 2003, Nature, 423, 847
Hjorth J., et al., 2005, Nature, 437, 859
Janiuk A., Perna R., Di Matteo T., Czerny, B.
2004, MNRAS, 355, 950
Janka H.-T., Eberl T., Ruffert M., Fryer C. L.
1999, ApJ, 527L, 39
Kohri K., Mineshige S., 2002, ApJ, 577, 311
Kohri K., Narayan R., Piran T., 2005, ApJ, 629,
Lattimer J.M., Swesty F.D., 1991, NuPhA, 535,
Lee W.H., Ramirez-Ruiz E., 2002, ApJ, 577, 893
Lee W.H., Ramirez-Ruiz E., Page D., 2005, ApJ,
632, 421
Lightman A.P., Eardley D.M., 1974, ApJL, 187,
Liu, T., Gu, W.-M., Xue, L., & Lu, J.-
F. 2007, ArXiv Astrophysics e-prints,
arXiv:astro-ph/0702186
MacFadyen A. I., Woosley S. E. 1999, ApJ, 524,
Meszaros P., 2006, Rept.Prog.Phys., 69, 2259
Meyer F., Meyer-Hofmeister E., 1981, A&A, 104,
Narayan R., Paczynski B., Piran T., 1992, ApJ,
395, L83
Narayan R., Piran T., Kumar P., 2001, ApJ, 557,
Ohsuga K., 2006, ApJ, 640, 923
Qian Y.Z., Woosley S.E., 1996, ApJ, 471, 331
Quataert E., Gruzinov A., 2000, ApJ, 539, 809
Paczynski B., 1986, ApJ, 308, L43
Paczynski B., 1991, AcA, 41, 257
Paczynski B., 1998, ApJ, 494, L45
Perna R., Armitage P., Zhang B. 2006, ApJ, 636L,
Piran T., 2005, Rev.Mod.Phys., 76, 1143
Popham R., Narayan R., 1995, ApJ, 442, 337
Popham R., Woosley S.E., Fryer C., 1999, ApJ,
518, 356
Pringle J.E., 1977, MNRAS, 177, 65
Proga D., MacFadyen A.I., Armitage P.J., Begel-
mann M.C., 2003, ApJ, 599, L5
Reddy S., Prakash M., Lattimer J.M., 1998, Phys.
Rev. D58, 300
Rosswog S., Liebendoerfer M. 2003, MNRAS, 342,
Rosswog S., Speith, R., Wynn G. A. 2004, MN-
RAS, 351, 1121
http://arxiv.org/abs/astro-ph/0508115
http://arxiv.org/abs/astro-ph/0607145
http://arxiv.org/abs/astro-ph/0702186
Rosswog S. 2005, ApJ, 634, 1202
Ruffert M., Janka H.-T., 1999, A&A, 344, 573
Ruffert M., Janka H.-T., 2001, A&A, 380, 544
Ruffert M., Janka H.-T., Schaefer G., 1996, A&A,
311, 532
Sawyer, R.F. 2003, Physical Review D, vol. 68,Is-
sue 6, id. 06300
Setiawan S., Ruffert M., Janka H.-Th. 2004, MN-
RAS, 352, 753
Setiawan S., Ruffert M., Janka H.-Th. 2005, A&A,
in press (astro-ph/0509300)
Shakura N.I, Sunyaev R., 1976, MNRAS, 175, 613
Shen H., Toki H., Oyamatsu K., Sumiyoshi K.,
1998a, PThPh, 100, 1013
Shen H., Toki H., Oyamatsu K., Sumiyoshi K.,
1998b, NuPhA, 637, 435
Surman R., McLaughlin G.C., 2004, ApJ, 603, 611
Smak J., 1984, AcA, 34, 161
Taam R.E., Lin D.N.C., 1984, ApJ, 287, 761
Siemiginowska A., 1988, AcA, 38, 21
Stanek K. Z., et al., 2003, ApJ, 591, L17
Szuszkiewicz E., 1990, MNRAS, 244, 377
Turner N.J., 2004, ApJ, 605, L45
Villasenor J.S., et al., 2005, Nature, 437, 855
Watarai K., Mineshige S., 2001, PASJ, 53, 915
Witt H. J., Jaroszynski M., Haensel P., Paczynski
B., Wambsganss J., 1994, ApJ, 422, 219
Yakovlev D.G., Kaminker A.D., Gnedin O.Y.,
Haensel P., 2001, Phys. Rep., 354, 1
Yuan Y., 2005, Phys. Rev. D, 72, 013007
Yuan Y., Heyl J.S., 2005, MNRAS, 360, 1493
Zhang B., 2007, astro-ph/0701520
This 2-column preprint was prepared with the AAS LATEX
macros v5.2.
http://arxiv.org/abs/astro-ph/0509300
http://arxiv.org/abs/astro-ph/0701520
	Introduction
	Neutrino-cooled accretion disks
	Initial disc configuration: 1-D Hydrodynamics
	The equation of state
	Composition and chemical equilibrium
	Neutrino cooling
	Energy and momentum Conservation
	Time evolution
	Numerical method
	Results
	EOS solutions for a given temperature and density
	The steady-state disc structure
	Composition and Chemical potentials
	Cooling and heating rates
	Stability analysis: instabilities at high-"705FM
	Time dependent solutions
	Discussion
	The unstable neutrino-opaque disc
	Comparison with previous work
	Limitations of our model
	Observational consequences
	Appendix
