text:
# Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?

## 1 Introduction

Research on gamma-ray bursts (GRBs) has made rapid progress over the past decade, particularly after the discovery of afterglow by Beppo-SAX (see recent reviews by Mészáros, Rees, & Wijers 1999; Piran 2000; Mészáros 2001). The current common belief is that GRBs are the most energetic explosions which ever occurred in the Universe, releasing energy of $`10^{51}10^{53}`$ erg only in a few to few tens of seconds. It is also widely argued that GRBs result from the conversion of the kinetic energy of ultra-relativistic particles created within a fire ball. A relativistic fireball shock model was proposed to resolve the compactness problem and has had some success in explaining how the GRB and afterglow radiation arise (Rees & Mészáros 1992, 1994; Wijers & Galama 1999).

However, the central engines of GRBs creating initial hot plasma or extremely energetic particles are not well understood yet due mainly to the fact that they are hidden from our view. It is usually argued that relativistic phenomena should be somehow involved, since otherwise it is difficult to explain huge fluence and rapid burst-like bust profiles. (Narayan, Paczyński & Piran 1992). Along this line, many interesting possibilities have been proposed so far, such as (1) mergers of double neutron star (NS) binaries or NS and BH (black hole) (Paczyński 1986, 1991; Goodman 1986; Eichler et al. 1989), (2) BH-WD (white dwarf) merger (Fryer et al. 1999); (3) BH-helium star merger (Fryer et al. 1999); (4) failed supernova (or collapsars, Woosley 1993; Paczyński 1998; MacFadyen & Woosley 1999), and (5) magnetar, rapidly spinning neutron star with extremely large magnetic fields (Usov 1992, 1994; Thompson 1994; Spruit 1999). It is of great importance to note that almost all the models (except for the magnetar hypothesis) predicts a similar configuration as an end result; namely, formation of a few solar mass black hole surrounded by a temporary debris torus (or disk) with mass of $`0.011.0M_{}`$ whose accretion can provide a sudden release of gravitational energy (Mészáros, Rees, & Wijers 1999). Then, the key issue will be to understand the properties of such a compact and massive disk associated with a huge mass accretion rate.

Accretion models for gamma-ray bursts were first considered by Narayan, Paczyński & Piran (1992), and recently discussed in more details by Popham, Woosley, & Fryer (1999) and Narayan, Piran, & Kumar, (2001, hereafter NPK). According to them, neutrino cooling dominates over radiative loss at very huge accretion rates (see, however, Chevalier 1996; see also Ruffert & Janka 1999 for a numerical simulation of a binary merger including neutrino losses. Therefore, such disks (or flows) are called as neutrino-dominated accretion flow (NDAF). In the present study, we elucidate the theory of NDAF based on the $`\alpha `$ viscosity prescription, paying special attention to their thermal structure and the implications of expected huge magnetic fields. We present the resultant equilibrium solutions and properties of neutrino-cooled disks in section 2. We then discuss the effects of huge magnetic fields and implications on neutrino emission in section 3. Final section is devoted to conclusions.

## 2 Structure of Neutrino-Cooled Disk

### 2.1 $`\mathrm{\Sigma }T`$ diagram

In this study we adopt the following basic equations based on the Newtonian dynamics, for simplicity (see, e.g., Kato, Fukue, & Mineshige 1998). We use the cylindrical coordinates, $`(r,\phi ,z)`$. Then, the mass conservation in the vertical and radial direction is, respectively, represented by

$`\rho ={\displaystyle \frac{\mathrm{\Sigma }}{2H}}\mathrm{and}\dot{M}=2\pi rv_r\mathrm{\Sigma },`$ (1)
where $`\rho `$ is the matter density, $`\mathrm{\Sigma }`$ is the surface density of the disk, $`H`$ is the disk half-thickness, $`\dot{M}`$ is the mass accretion rate, and $`v_r`$ is the radial velocity (negative for inflow). The disk half-thickness is, from hydrostatic balance in the vertical direction, given by

$`H=c_\mathrm{s}/\mathrm{\Omega },`$ (2)
where the (isothermal) sound velocity $`c_\mathrm{s}`$ is defined by

$`c_\mathrm{s}\sqrt{p/\rho },`$ (3)
with the pressure $`p`$. The angular velocity, $`\mathrm{\Omega }`$, is in the Newtonian approximation given by

$`\mathrm{\Omega }=\sqrt{GM_{\mathrm{BH}}/r^3},`$ (4)
with $`G`$ being the gravitational constant and $`M_{\mathrm{BH}}`$ being the black-hole mass. The pressure is composed of three terms:

$`p=p_{\mathrm{rad}}+p_{\mathrm{gas}}+p_\mathrm{d},`$ (5)
where $`p_{\mathrm{rad}}`$ is the radiation pressure, $`p_{\mathrm{gas}}`$ is the gas pressure of the nonrelativistic and non-degenerate particles, and $`p_\mathrm{d}`$ is the pressure of the degenerate particles. The radiation pressure is

$`p_{\mathrm{rad}}={\displaystyle \frac{a}{6}}g_{}T^4,`$ (6)
where $`g_{}`$ is statistical degree of freedom in the radiation (=2 for $`\gamma `$, or 11/2 for $`\gamma `$ and $`e^+e^{}`$ plasma) and $`a`$ is the radiation constant. The gas pressure is

$`p_{\mathrm{gas}}={\displaystyle \underset{i}{}}n_ik_\mathrm{B}T,`$ (7)
where $`k_\mathrm{B}`$ is the Boltzmann constant, $`n_i`$ is the number density of particle “$`i`$”, and the suffix $`i`$ runs the nonrelativistic and non-degenerate particles at the temperature (e.g., $`e^{}`$ and nucleons at $`T\stackrel{<}{}\text{ }m_\mathrm{e}c^2/k_\mathrm{B}`$). The pressure of the degenerate particles is written as

$`p_\mathrm{d}=p_{\mathrm{d},\mathrm{nonrel}}+p_{\mathrm{d},\mathrm{rel}}.`$ (8)
Then,

$`p_{\mathrm{d},\mathrm{nonrel}}={\displaystyle \underset{i=\mathrm{nonrel}}{}}{\displaystyle \frac{1}{5}}\left({\displaystyle \frac{6\pi ^2}{g_i}}\right)^{2/3}{\displaystyle \frac{\mathrm{}^2n_i^{5/3}}{m_i}},`$ (9)
is for nonrelativistic particles which satisfy the condition of the degeneracy,

$`n_ig_i(m_ik_\mathrm{B}T/2\pi \mathrm{}^2)^{3/2},`$ (10)
for $`Tm_ic^2/k_\mathrm{B}`$ and $`\mu _im_ic^2`$ with their mass $`m_i`$, statistical degree of freedom $`g_i`$ and their chemical potential $`\mu _i`$, and

$`p_{\mathrm{d},\mathrm{rel}}={\displaystyle \underset{i=\mathrm{rel}}{}}{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{6\pi ^2}{g_i}}\right)^{1/3}\mathrm{}^2cn_i^{4/3},`$ (11)
is for relativistic particles which satisfy the condition of the degeneracy,

$`n_ig_i(k_\mathrm{B}T/\mathrm{}c)^3/\pi ^2,`$ (12)
for $`Tm_ic^2/k_\mathrm{B}`$ or $`\mu _im_ic^2`$ . In this study, for simplicity we assume the complete degeneracy whenever the above condition of the degeneracy is satisfied.

In Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?, we plot the contours of the matter density on the ($`\mathrm{\Sigma }`$, $`T`$) plane. From the figure we can read off the matter density. It also discriminates the regions where each pressure component is dominant and also the regions where electrons and/or nucleons are degenerate. We understand that degeneracy pressure is important at high density and low temperature regimes. Such regimes inevitably appear in very massive disks, and one can never neglect the contributions by degeneracy pressure in GRBs unlike the cases of binary systems or galactic centers. We also notice that one should distinguish the regime where only electrons are degenerate from the one where both electrons and nucleons are degenerate.

### 2.2 chemical potential of electrons

When electrons are degenerate, there emerges an important consequence on the process of neutrino emission, that is, e.g., since the pair creation of electron and positron in the electromagnetic thermal bath is suppressed by the charge neutrality (see below), the neutrino emission is also suppressed (This is not discussed in NPK). The electron chemical potential $`\mu _e`$ is determined by the condition of charge neutrality,

$`n_p=n_e^{}n_{e^+},`$ (13)
with the number density of electron and positron,

$`n_e^{}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑pp^2{\displaystyle \frac{1}{e^{(\sqrt{p^2c^2+m_e^2c^4}\mu _e)/k_\mathrm{B}T}+1}},`$ (14)
$`n_{e^+}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑pp^2{\displaystyle \frac{1}{e^{(\sqrt{p^2c^2+m_e^2c^4}+\mu _e)/k_\mathrm{B}T}+1}}.`$ (15)
In the following sections, it is convenient to introduce a dimensionless chemical potential of electrons which is defined by

$`\eta _e=\mu _e/k_\mathrm{B}T.`$ (16)
Then, the complete degeneracy of electrons corresponds to $`\eta _e1`$, and the charge neutrality in Eq. (13) becomes

$`n_pn_e^{}{\displaystyle \frac{1}{3\pi ^2}}\left({\displaystyle \frac{\mu _e}{\mathrm{}c}}\right)^3.`$ (17)
The chemical potentials of protons and neutrons are defined through the relation,

$`n_p`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑pp^2{\displaystyle \frac{1}{e^{(\sqrt{p^2c^2+m_p^2c^4}\mu _p)/k_\mathrm{B}T}+1}},`$ (18)
$`n_n`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}^3\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑pp^2{\displaystyle \frac{1}{e^{(\sqrt{p^2c^2+m_n^2c^4}\mu _n)/k_\mathrm{B}T}+1}}.`$ (19)
If the nucleons are not degenerate, we have

$`n_i`$ $`=`$ $`2\left({\displaystyle \frac{m_ik_\mathrm{B}T}{2\pi \mathrm{}^2}}\right)^{3/2}e^{(m_ic^2\mu _i)/k_\mathrm{B}T},`$ (20)
for $`i=p`$, and $`n`$. On the other hand, if nucleons are degenerate, we have

$`n_i={\displaystyle \frac{1}{3\pi ^2}}\left({\displaystyle \frac{2\stackrel{~}{\mu }_im_i}{\mathrm{}^2}}\right)^{3/2},`$ (21)
for the complete degenerate limit with

$`\stackrel{~}{\mu }_i=\mu _im_ic^2.`$ (22)
In this study, we did not directly solve the exact equation in Eq. (13) to derive the chemical potential of electrons. Instead, we performed the following approximated computations in the complete-degeneracy regime ($`\eta _e1`$). As will be shown later, nucleons are in $`\beta `$-equilibrium at a high temperature $`T5.0\times 10^{10}\mathrm{K}(\alpha /0.1)^{1/5}(r/4r_g)^{3/10}(M_{\mathrm{BH}}/3M_{})^{1/5}`$ where the neutrino emissions are effective and mainly contribute to the cooling process. Then, the chemical equilibrium is realized,

$`\mu _p+\mu _e=\mu _n,`$ (23)
or

$`\stackrel{~}{\mu }_p+\mu _eQ=\stackrel{~}{\mu }_n,`$ (24)
where $`Q=(m_nm_p)c^21.29`$ MeV. In the regime where only electrons are degenerate, $`\beta `$-equilibrium means

$`{\displaystyle \frac{n_p}{n_n}}=\mathrm{exp}[(Q\mu _e)/(k_\mathrm{B}T)].`$ (25)
The electron fraction $`Y_en_p/(n_p+n_n)`$ is related with the above ratio by $`Y_e=1/(1+(n_p/n_n)^1)`$. Combined with the condition of the charge neutrality, i.e., Eq. (17), we may only solve the equation,

$`{\displaystyle \frac{\rho }{m_p}}\mathrm{exp}(\mu _e/(k_\mathrm{B}T))={\displaystyle \frac{1}{3\pi ^2}}\left({\displaystyle \frac{\mu _e}{\mathrm{}c}}\right)^3,`$ (26)
for $`\mu _ek_\mathrm{B}TQ`$. On the other hand, it is relatively easy for us to calculate $`\mu _e`$ in the regime where both electrons and nucleons are degenerate. From the charge neutrality (Eq. (17)) and chemical equilibrium (Eq.(24)), we find that $`\mu _e\stackrel{~}{\mu }_n\stackrel{~}{\mu }_p`$ for $`k_\mathrm{B}T\mu _e<𝒪(10^2)\mathrm{MeV}`$. Then, from Eq. (21) we can approximately estimate the chemical potential of electrons as

$`\mu _e`$ $``$ $`\stackrel{~}{\mu _n}`$ (27)
$``$ $`6.628\mathrm{MeV}(\rho /10^{13}\mathrm{g}\mathrm{cm}^3)^{2/3},`$ (28)
in the complete-degeneracy limit of both electrons and nucleons.

We also plot the contours of the chemical potential of electrons in units of MeV ($`\mu _{e,\mathrm{MeV}}`$ = 10, and 20) in Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?. The chemical potential of electrons play important roles to estimate the emission rates of neutrinos which will be discussed in the next subsection.

### 2.3 Heating and cooling rates

The most important relation in the accretion disk theory is the energy balance between the heating and cooling processes. Distinct branches (such as a standard disk or a slim disk, etc.) appear because of different heating and/or cooling sources being dominant. According to the standard $`\alpha `$ viscosity model, we express the vertically integrated heating rate (over a half thickness, $`H`$) as

$`Q^+=Q_{\mathrm{vis}}^+={\displaystyle \frac{9}{8}}\nu \mathrm{\Sigma }\mathrm{\Omega }^2,`$ (29)
where $`Q_{\mathrm{vis}}^+`$ denotes the viscous heating rate per unit surface area, the kinetic viscosity $`\nu `$ is related to the viscosity parameter $`\alpha `$ by

$`\nu ={\displaystyle \frac{2}{3}}\alpha c_\mathrm{s}H.`$ (30)
If we assume the angular momentum conservation, we obtain

$`\nu \mathrm{\Sigma }={\displaystyle \frac{\dot{M}}{3\pi }}\left(1\sqrt{{\displaystyle \frac{r_{\mathrm{in}}}{r}}}\right),`$ (31)
where $`r_{\mathrm{in}}`$ is the radius of the inner edge of the disk (= 3 $`r_\mathrm{g}`$). From Eqs. (29), (30) and (31), we can relate the heating rate with the mass accretion rate;

$`\dot{M}={\displaystyle \frac{8\pi }{3}}\mathrm{\Omega }^2\left(1\sqrt{{\displaystyle \frac{r_{\mathrm{in}}}{r}}}\right)^1Q^+.`$ (32)
The cooling rate is, on the other hand, summation of three major contributions;

$`Q^{}=Q_{\mathrm{rad}}^{}+Q_{\mathrm{adv}}^{}+Q_\nu ^{},`$ (33)
where $`Q_{\mathrm{rad}}^{}`$ is the radiative cooling rate, $`Q_{\mathrm{adv}}^{}`$ is the advective energy transport (Abramowicz et al. 1988), and $`Q_\nu ^{}`$ is the cooling rate due to neutrino loss. Note that instead of including advective energy transport in the energy equation NPK adopted the view of CDAF (convection-dominated accretion flow; Igumenshchev, Abramowicz, & Narayan 2000, Narayan, Igumenshchev, & Abramowicz, 2000, Quataert & Gruzinov 2000). We, in the present study, retain the classical picture based on the (vertically) one-zone treatment (e.g., Kato, Fukue, & Mineshige 1998), since it is not yet clear if CDAF provides precise description to the flow structure not only in the low-luminosity regimes but also in the hyper-accretion regimes. However, such a distinction is not essential here, since we are concerned with the regimes of even higher accretion rates, in which neutrino emission is substantial (discussed later).

The radiative cooling rate is

$`Q_{\mathrm{rad}}^{}={\displaystyle \frac{g_{}\sigma _\mathrm{s}T^4}{2\tau _{\mathrm{tot}}}},`$ (34)
where $`\sigma _\mathrm{s}=\pi ^2k_\mathrm{B}^4/(60\mathrm{}^3c^2)`$ is the Stefan-Boltzmann constant and the optical depth, $`\tau _{\mathrm{tot}}`$, is given by

$`\tau _{\mathrm{tot}}=\kappa _\mathrm{R}\rho H={\displaystyle \frac{\kappa _\mathrm{R}\mathrm{\Sigma }}{2}},`$ (35)
with $`\kappa _\mathrm{R}`$ being the Rosseland-mean opacity,

$`\kappa _\mathrm{R}=0.40+0.64\times 10^{23}\left({\displaystyle \frac{\rho }{\mathrm{g}\mathrm{cm}^3}}\right)\left({\displaystyle \frac{T}{\mathrm{K}}}\right)^3\mathrm{g}^1\mathrm{cm}^2.`$ (36)
The advective cooling rate is given by (Kato, Fukue, & Mineshige 1998),

$`Q_{\mathrm{adv}}^{}=\mathrm{\Sigma }Tv_r{\displaystyle \frac{ds}{dr}},`$ (37)
where the radial velocity is $`v_r=\dot{M}/(2\pi r\mathrm{\Sigma })`$ and $`s`$ denotes entropy per particle,

$`s=\left(s_{\mathrm{rad}}+s_{\mathrm{gas}}\right)/\rho .`$ (38)
Here, the entropy density of the radiation is

$`s_{\mathrm{rad}}={\displaystyle \frac{2}{3}}ag_{}T^3,`$ (39)
and entropy density of the gas (i.e., nonrelativistic particles) is

$`s_{\mathrm{gas}}={\displaystyle \underset{i}{}}n_i\left({\displaystyle \frac{5}{2}}+\mathrm{ln}\left[{\displaystyle \frac{\mathrm{g}_i}{n_i}}\left({\displaystyle \frac{m_iT}{2\pi }}\right)^{3/2}\right]\right),`$ (40)
where the suffix $`i`$ runs over nonrelativistic nucleon and electron, $`g_i`$ is the statistical degree of freedom of the particle. In the following, we approximated $`ds/dr`$ as $`s/r`$. Note that since the entropy of degenerate particles is small in the complete-degeneracy limit, we neglect it here.

The neutrino cooling rate is composed of four terms;

$`Q_\nu ^{}=(\dot{q}_{Ne}+\dot{q}_{e^+e^{}}+\dot{q}_{\mathrm{brems}}+\dot{q}_{\mathrm{plasmon}})H,`$ (41)
where $`\dot{q}_{Ne}`$ is the electron-positron capture rate by a nucleon “N”, $`\dot{q}_{e^+e^{}}`$ is the electron-positron pair annihilation rate, $`\dot{q}_{\mathrm{brems}}`$ is the nucleon-nucleon bremsstrahlung rate, and $`\dot{q}_{\mathrm{plasmon}}`$ is the rate of plasmon decays. Note that NPK considered only the case of the non-degenerate electrons in first two terms on the right-hand side of Eq. (41). However, the case of the degenerate electrons in the first two terms and the other two terms could be also important.

The electron-positron capture rate is represented by two terms:

$`\dot{q}_{Ne}=\dot{q}_{p+e^{}n+\nu _e}+\dot{q}_{n+e^+p+\overline{\nu }_e},`$ (42)
with

$`\dot{q}_{p+e^{}n+\nu _e}`$ $`=`$ $`{\displaystyle \frac{G_\mathrm{F}^2}{2\pi ^3\mathrm{}^3c^2}}(1+3g_A)n_p{\displaystyle _Q^{\mathrm{}}}𝑑E_eE_e\sqrt{E_e^2m_e^2c^4}\left(E_eQ\right)^3{\displaystyle \frac{1}{e^{\left(E_e\mu _e\right)/k_\mathrm{B}T}+1}},`$ (43)
$`\dot{q}_{n+e^+p+\overline{\nu }_e}`$ $`=`$ $`{\displaystyle \frac{G_\mathrm{F}^2}{2\pi ^3\mathrm{}^3c^2}}(1+3g_A)n_n{\displaystyle _{m_ec^2}^{\mathrm{}}}𝑑E_eE_e\sqrt{E_e^2m_e^2c^4}\left(E_e+Q\right)^3{\displaystyle \frac{1}{e^{\left(E_e+\mu _e\right)/k_\mathrm{B}T}+1}},`$ (44)
where $`G_\mathrm{F}`$ is the Fermi coupling constant ($`=2.302\times 10^{22}\mathrm{cm}\mathrm{MeV}^1`$), and $`g_A`$ is the axial-vector coupling constant of nucleon ($``$ 1.39) which is normalized by the experimental value of neutron lifetime $`\tau _n886.7`$ s. In the non-degeneracy limit ($`\mu _ek_\mathrm{B}T`$), it is easily estimated by

$`\dot{q}_{Ne}=9.2\times 10^{33}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^6\left({\displaystyle \frac{\rho }{10^{10}\mathrm{g}\mathrm{cm}^3}}\right).`$ (45)
On the other hand, it is a little complicated to estimate the electron capture rate in the electron-degeneracy regime. In the only electron-degeneracy regime, we have

$`\dot{q}_{Ne}=1.1\times 10^{31}\eta _e^9\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^9,`$ (46)
for complete degeneracy limit of electrons. Note that in this limit, independently we can also obtain the chemical potential of electrons by the balance between the two dominant weak interaction rates, i.e., $`\dot{q}_{p+e^{}n+\nu _e}=5.3\times 10^{30}\eta _e^9\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1(T/10^{11}\mathrm{K})^9`$ and $`\dot{q}_{n+e^+p+\overline{\nu }_e}=2.6\times 10^{37}\mathrm{exp}(m_e/(k_\mathrm{B}T))\mathrm{exp}(\eta _e)\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1(\rho /10^{13}\mathrm{g}\mathrm{cm}^3)(T/10^{11}\mathrm{K})^6`$. Compared with the chemical potential which was obtained in Eq (26) where we assumed the complete $`\beta `$–equilibrium, they agree with each other just within a few percent, e.g., $`\eta _e3.5`$ at $`\rho =10^{13}\mathrm{g}\mathrm{cm}^3`$ and $`T`$ = 10 MeV/$`k_\mathrm{B}`$. Therefore, our assumptions of the complete-degeneracy limits would be reasonable.

When both electrons and nucleons are degenerate, the reactions $`p+e^{}n+\nu _e`$ are suppressed by the Fermi blocking of degenerate nucleons in the final states. Then, for example the $`dE_e`$ integration in Eq. (43) is limited as \[$`\mu _e`$, $`\mathrm{}`$\], and we approximately obtain the relatively small cooling rate, $`\dot{q}_{Ne}=5.0\times 10^{32}\eta _e^7\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1(T/10^{11}\mathrm{K})^9`$. Compared with Eq. (46), of course, we find that this rate is smaller by a factor of $`1/\eta _e^2`$ in the complete degeneracy regime, i.e., $`\eta _e1`$.

The electron-positron pair annihilation rate through $`e^++e^{}\nu +\overline{\nu }`$ is

$`\dot{q}_{e^+e^{}}=4.8\times 10^{33}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^9,`$ (47)
in the non-degeneracy regime, e.g., see Itoh et al. (1989, 1990). When the electrons are degenerate, the electron-positron pair annihilation rate is too small to compete the other cooling process, and we can neglect it.

The nucleon-nucleon bremsstrahlung rate through $`n+nn+n+\nu +\overline{\nu }`$ is represented by

$`\dot{q}_{\mathrm{brems}}=3.4\times 10^{33}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^8\left({\displaystyle \frac{\rho }{10^{13}g\mathrm{cm}^3}}\right)^{1/3},`$ (48)
in the degeneracy regime of nucleons (Hannestad & Raffelt 1998), and

$`\dot{q}_{\mathrm{brems}}=1.5\times 10^{33}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^{5.5}\left({\displaystyle \frac{\rho }{10^{13}g\mathrm{cm}^3}}\right)^2,`$ (49)
in the non-degeneracy regime of nucleons (Hannestad & Raffelt 1998; Burrows et al. 2000).

It is also known that the plasmon decay is effective in high densities and high electron degeneracy region (Schinder et al. 1987). The decay rate of the transverse plasmons which are normal photons interacting with the electron gas through $`\stackrel{~}{\gamma }\nu _e+\overline{\nu }_e`$ is estimated by Ruffert, Janka & Schäfer (1996) as

$`\dot{q}_{\mathrm{plasmon}}=1.5\times 10^{32}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^9\gamma _p^6e^{\gamma _p}(1+\gamma _p)\left(2+{\displaystyle \frac{\gamma _p^2}{1+\gamma _p}}\right),`$ (50)
where $`\gamma _p=5.565\times 10^2\sqrt{(\pi ^2+3\eta _e^2)/3}`$. In particular, $`\stackrel{~}{\gamma }\nu _e+\overline{\nu }_e`$ is a dominant process by a factor of $`163`$ compared with the other flavor ($`\nu _\mu \overline{\nu }_\mu `$ or $`\nu _\tau \overline{\nu }_\tau `$ ).

Note that in NPK only Eqs. (45) and (47) cases were considered in $`\dot{q}_\nu `$. To see what component of the cooling rates mainly contributes to the process, in Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts? we plot the most dominant one on the ($`\mathrm{\Sigma }`$, $`T`$) plane. It is interesting to note that the region where advective energy transport dominates covers a rather wide, large $`\mathrm{\Sigma }`$ and high $`T`$ part of the $`\mathrm{\Sigma }T`$ diagram. This is because $`Q_{\mathrm{adv}}`$ ($`T^{16}`$) has a stronger $`T`$ dependence than $`Q_{\mathrm{rad}}`$ ($`T^4`$) (shown later). Also notice that neutrino cooling is essential only at very high temperatures ($`T10^{11}\mathrm{K}`$) and very high matter density $`\rho 10^{10}`$$`10^{13}\mathrm{g}\mathrm{cm}^3`$. The striped region represents the place where heavy elements are produced which are expected from the point of view in the equation of state of nuclear matter (Shen et al, 1998; Ishizuka, Ohnishi & Sumiyoshi 2002). Here, we adopt the case of the electron fraction $`Y_e0.1`$ which the numerical simulations of coalescing neutron stars predict, e.g., see Ruffert et al (1997). Note that the left side of the boundary of the striped region is just virtual one. Of course, the tail realistically continues to the left further. As a result, there exist just a small amount of free nucleons. Namely, the neutrino processes which we consider in this study no longer become effective there.

Next, in Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?, we plot the heating and cooling rates as a function of the temperature for some representative values of surface density: (a) $`\mathrm{\Sigma }=10^{3.5}\mathrm{g}\mathrm{cm}^2`$, (b) $`\mathrm{\Sigma }=10^8\mathrm{g}\mathrm{cm}^2`$, and (c) $`\mathrm{\Sigma }=10^{20}\mathrm{g}\mathrm{cm}^2`$. In Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?(a), we find three intersection points. They contain two points in which the energy balance is realized due to the radiative or advective cooling processes. As for Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?(b), there exists only one intersection point in which the energy balance is realized because of the advective cooling processes. It is interesting that in Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?(c), a new branch of the energy balance appears at $`T\stackrel{>}{}\text{ }m_\mathrm{e}c^2/k_\mathrm{B}`$, since the neutrino-cooling process becomes effective. By plotting only the values at those intersecting points on the ($`\mathrm{\Sigma }`$, $`T`$) or ($`\mathrm{\Sigma }`$, $`\dot{M}`$) plane for various values of $`\mathrm{\Sigma }`$, we obtain a sequence of the thermal equilibrium solutions; that is the thermal equilibrium curve.

### 2.4 Various time scales

Here it is necessary for us to check the time scale of the physical processes of the neutrino cooling. The dynamical timescale is represented by the accretion time,

$`\tau _{\mathrm{acc}}`$ $``$ $`{\displaystyle \frac{1}{\alpha }}\sqrt{{\displaystyle \frac{r^3}{GM_{\mathrm{BH}}}}}\left({\displaystyle \frac{r}{H}}\right)^2,`$ (51)
$``$ $`1.3\times 10^2\mathrm{s}\left({\displaystyle \frac{\alpha }{0.1}}\right)^1\left({\displaystyle \frac{r}{4r_g}}\right)^{3/2}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right),`$ (52)
where we assumed that the disk half-thickness is approximately $`Hr/2`$. On the other hand, the time scale of weak interaction of a nucleon off the background electrons and positrons is estimated by

$`\tau _\beta `$ $``$ $`\left(\sigma _{Ne}n_ec\right)^1`$ (53)
$``$ $`1.0\times 10^4\mathrm{s}\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^5,`$ (54)
where we assume that the cross section is $`\sigma _{Ne}G_\mathrm{F}^2(k_\mathrm{B}T)^2`$, and the electron number density is $`n_e(k_\mathrm{B}T/\mathrm{}c)^3`$. Then, the condition $`\tau _\beta \tau _{\mathrm{acc}}`$ is realized when the following relation is satisfied,

$`T3.8\times 10^{10}\mathrm{K}\left({\displaystyle \frac{\alpha }{0.1}}\right)^{1/5}\left({\displaystyle \frac{r}{4r_g}}\right)^{3/10}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^{1/5}.`$ (55)
The interaction rate of a neutrino off the background nucleons is roughly given by

$`\mathrm{\Gamma }_{N\nu }G_\mathrm{F}^2(k_\mathrm{B}T)^2n_Nc,`$ (56)
with the nucleon density $`n_N`$. Since $`Hr/2`$, the ratio of the interaction time ($`1/\mathrm{\Gamma }_{N\nu }`$) and the neutrino crossing time ($`H/c`$) is estimated by

$`\mathrm{\Gamma }_{N\nu }\times H/c42.6\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^2\left({\displaystyle \frac{\rho }{10^{13}\mathrm{g}\mathrm{cm}^3}}\right)\left({\displaystyle \frac{r}{4r_g}}\right)\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right).`$ (57)
Since the nucleon to electron-positron ratio is $`n_N/n_e𝒪(10)(\rho /10^{13}\mathrm{g}\mathrm{cm}^3)/(T/10^{11}\mathrm{K})^3`$, we find that the neutrino interaction is rapid enough to realize $`\beta `$–equilibrium which we have assumed in the previous subsections. When $`\mathrm{\Gamma }_{N\nu }H/c`$ is greater than unity, we should modify the neutrino cooling rate in Eq. (41). To correctly estimate it in that case, we should divide the right-hand side of Eq. (41) by $`\mathrm{\Gamma }_{N\nu }H/c`$. However, since the neutrino cooling rate becomes $`𝒪(10^2)`$ times larger than the heating rate suddenly at intersection points between the cooling and the heating rates, to estimate the equilibrium curve, our simple treatments in the previous section are not so wrong and the results are not changed.

Next, we estimate the vertical diffusion time of neutrinos. It is given by

$`\tau _{\mathrm{diff}}`$ $``$ $`{\displaystyle \frac{(H/c)^2}{\mathrm{\Gamma }_{N\nu }^1}},`$ (58)
$``$ $`2.6\times 10^3\mathrm{s}\left({\displaystyle \frac{T}{10^{11}\mathrm{K}}}\right)^2\left({\displaystyle \frac{\rho }{10^{13}\mathrm{g}\mathrm{cm}^3}}\right)\left({\displaystyle \frac{r}{4r_g}}\right)^2\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^2.`$ (59)
Compared it with Eq. (51), the condition of $`\tau _{\mathrm{diff}}\tau _{\mathrm{acc}}`$ is translated into

$`T2.3\times 10^{11}\mathrm{K}\left({\displaystyle \frac{\alpha }{0.1}}\right)^{1/2}\left({\displaystyle \frac{\rho }{10^{13}\mathrm{g}\mathrm{cm}^3}}\right)^{1/2}\left({\displaystyle \frac{r}{4r_g}}\right)^{1/4}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^{1/2}.`$ (60)
Therefore, the neutrino cooling process is effective and realistic in the equilibrium solutions of the accretion disk at a high temperature $`T10^{11}\mathrm{K}`$ and a high density $`\rho 10^{13}\mathrm{g}\mathrm{cm}^3`$ which will be studied in the next subsection.

### 2.5 Thermal Equilibrium Solutions

According to the procedure discussed in Sec. 2.3, we can find the thermal equilibrium solutions; that is, intersecting points of the heating and cooling curves. In Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts? we plot them on the ($`\mathrm{\Sigma }`$, $`T`$) plane for $`r=4r_g`$ (left panel) and for $`r=40r_g`$ (right panel), respectively. In the parameter region where neutrino cooling is effective, the dominant processes are the electron capture rate ($`N+eN^{}+\nu `$).

Using Eq. (32), we can also express the mass accretion rate $`\dot{M}`$ at the intersection points as a function of $`\mathrm{\Sigma }`$. The thermal equilibrium curves on the ($`\mathrm{\Sigma }`$, $`\dot{M}`$) plane are plotted in Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?. Here, mass accretion rate is normalized by the critical mass accretion rate,

$`\dot{M}_{\mathrm{crit}}=16L_{\mathrm{Edd}}/c^2,`$ (61)
where $`L_{\mathrm{Edd}}`$is the Eddington luminosity,

$`L_{\mathrm{Edd}}=4\pi GM_{\mathrm{BH}}m_\mathrm{p}c/\sigma ,`$ (62)
with proton mass $`m_\mathrm{p}`$ and the cross section $`\sigma `$ of the matter and the radiation through the electromagnetic scattering, i.e., $`L_{\mathrm{Edd}}/c^27.3\times 10^{17}M_{}\mathrm{s}^1(\sigma /\sigma _\mathrm{T})^1(M_{\mathrm{BH}}/M_{})`$ with the Thomson cross section $`\sigma _\mathrm{T}0.6652\mathrm{mb}`$.

There are 5 distinct branches seen in these plots (see Table 1). In the lower-left parts of both figures, the equilibrium sequence has an S shape (see Abramowicz et al. 1988), which arises because of changes in sources of pressure (gas and radiation pressure) and in sources of cooling (radiation and advection). In the upper branch of the S shape, advective energy transport takes over radiative cooling. Then, generated photons inside the disk take long time to go out from the disk surface so that photons are advected inward and are finally swallowed by a central black hole with accreting matter. Then, disk luminosity is largely reduced. Such trapping effects should occur in branches III and IV with even higher $`\dot{M}`$. Neutrino cooling is dominant only in the high temperature and density regimes (branch V), in which neutrinos, instead of radiation, can carry away the generated energy inside the disk. In branch V of Table 1, we wrote $`Q_\nu ^{}`$ in the case that only electrons are degenerate. On the other hand, in the case that both electrons and nucleons are degenerate, $`Q_\nu ^{}`$ becomes smaller by a factor of $`1/\eta _e^2`$. Then, the thermal equilibrium solution is just realized on the line of the condition of nucleon degeneracy. Namely, the $`T`$$`\mathrm{\Sigma }`$ relation of the branch V coincides with that of nucleon degeneracy. In the intermediate branch (IV) still advective cooling is dominant, however, note that entropy which is transported by advection is contributed by non-degenerate nucleons, since entropy of degenerate particles is little.

It is interesting to examine the stability of each branch. The criterion for the thermal instability is

$$\left(\frac{dQ^+}{dT}\right)_\mathrm{\Sigma }>\left(\frac{dQ^{}}{dT}\right)_\mathrm{\Sigma },$$
(63)
whereas that for the viscous (secular) instability is

$$\left(\frac{d\dot{M}}{d\mathrm{\Sigma }}\right)_{Q^+=Q^{}}<0,$$
(64)
(see, e.g., Kato, Fukue, & Matsumoto 1998). We understand immediately from Table 1 that only branch II (radiation pressure-dominant standard disk branch) is unstable for both modes and all other branches, including the NDAF branch, are stable for any modes, as already noted by NPK.

### 2.6 Properties of Neutrino-Cooled Disks

¿From the equilibrium solutions found in Sec. 2.5, we derive the following expressions for the temperature and density on the equatorial plane and for the corresponding mass-accretion rates as functions of radius for a given value of $`\mathrm{\Sigma }`$:

$`T=1.1\times 10^{11}\mathrm{K}\left({\displaystyle \frac{\mathrm{\Sigma }}{10^{20}\mathrm{g}\mathrm{cm}^2}}\right)^{4/7}\left({\displaystyle \frac{r}{4r_g}}\right)^{6/7}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^{2/7},`$ (65)
$`\rho =2.3\times 10^{13}\mathrm{g}\mathrm{cm}^3\left({\displaystyle \frac{\mathrm{\Sigma }}{10^{20}\mathrm{g}\mathrm{cm}^2}}\right)^{6/7}\left({\displaystyle \frac{r}{4r_g}}\right)^{9/7}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^{3/7},`$ (66)
and

$`\dot{M}=1.1\times 10^{17}\dot{M}_{\mathrm{crit}}\left({\displaystyle \frac{\alpha }{0.1}}\right)\left({\displaystyle \frac{\mathrm{\Sigma }}{10^{20}\mathrm{g}\mathrm{cm}^2}}\right)^{9/7}\left({\displaystyle \frac{r}{4r_g}}\right)^{15/14}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^{5/14}.`$ (67)
Note that NPK derived similar expressions but in terms of the disk mass, $`M_{\mathrm{disk}}`$, and the size of the initial torus, $`R_{\mathrm{out}}`$, which are roughly related to surface density as

$$\mathrm{\Sigma }\frac{2M_{\mathrm{disk}}}{\pi bR_{\mathrm{out}}^2}\frac{10^{20}}{b}\mathrm{g}\mathrm{cm}^2\left(\frac{M_{\mathrm{disk}}}{M_{}}\right)\left(\frac{R_{\mathrm{out}}}{5\times 10^6\mathrm{cm}}\right)^2,$$
(68)
where $`b`$ is a constant of the order of unity, which varies, depending on the geometry of the initial torus and we assumed uniform surface density distribution. We thus understand that neutrino cooling dominates only if a solar-mass material is contained within a disk with a size of $`5\times 10^6`$ cm. In other words, if the disk mass is less or if the disk size is larger, neutrino cooling never becomes dominant. This constraint is more severe than that obtained by NPK.

By comparing the left and right panels in Fig. Can neutrino-cooled accretion disk be an origin of gamma-ray bursts?, we notice that even at such a high surface density as $`\mathrm{\Sigma }10^{20}\mathrm{g}\mathrm{cm}^2`$, the neutrino cooling does not become dominant at a somewhat larger radius (i.e., $`r40r_g`$). This also supports the claim by NPK that NDAF solution only appears in a rather compact region around the center.

We can get a crude approximation for the strength of magnetic fields, by using the above results. The equi-partition field strength, for which magnetic energy is equal to gas energy, is

$`B_{\mathrm{eq}}`$ $`=`$ $`\left({\displaystyle \frac{8\pi \rho k_BT}{m_\mathrm{p}}}\right)^{1/2}`$ (69)
$``$ $`7.3\times 10^{16}\mathrm{G}\left({\displaystyle \frac{\alpha }{0.1}}\right)^{1/12}\left({\displaystyle \frac{\mathrm{\Sigma }}{10^{20}\mathrm{g}\mathrm{cm}^2}}\right)^{5/7}\left({\displaystyle \frac{r}{4r_g}}\right)^{15/14}\left({\displaystyle \frac{M_{\mathrm{BH}}}{3M_{}}}\right)^{5/14}.`$
Therefore, the strength of the magnetic field is very likely to exceed the critical value (e.g. Mészáros 1992),

$$B_{\mathrm{crit}}\frac{m_\mathrm{e}^2c^3}{e\mathrm{}}4.4\times 10^{13}\mathrm{G},$$
(70)
even if the fields strength is only 0.1 % of the equi-partition value (for which magnetic energy is only $`10^4`$ % of the particle energy). It is interesting to note that equilibrium $`B`$ value is about the critical field for the proton,

$$B_{\mathrm{crit},\mathrm{p}}=(m_\mathrm{p}/m_\mathrm{e})B_{\mathrm{crit}}8\times 10^{16}\mathrm{G}.$$
(71)
Then various QED (quantum-electrodynamics) effects should manifest themselves (see the next section).

## 3 Discussion

### 3.1 CDAF and slim disk

The present study is similar to that by NPK, although there are some differences. We newly take into account the effects of the electron and nucleon degeneracy in neutrino cooling processes. The resultant disk density and temperature differ in the neutrino-dominated regimes. In the present study, we also take the slim-disk picture, instead of the CDAF picture adopted by NPK, at moderately large $`\mathrm{\Sigma }`$ (and thus $`\dot{M}`$) ranges. As mentioned above, the properties of the hyper-critical accretion flow is not yet clear, since full multi-dimensional simulations coupled with radiation hydrodynamics (RHD) and possibly with magnetohydrodynamics (MHD) are finally necessary to settle on this issue. Probably convection plays a key role there (Agol et al. 2001), but RHD simulations so far made in the hyper-critical regimes are based on the approximation of the flux-limited diffusion (e.g. Eggum, Coroniti, & Katz 1988; Fujita & Okuda 1998), which may not hold within highly turbulent media. Also, magnetic reconnection, which is not easy to simulate in realistic situations with extremely high magnetic Reynolds number, seems to play an important role (see below).

Our conclusion basically agrees with that by NPK in that NDAF appears at small radii only for high $`M_{\mathrm{disk}}/r_{\mathrm{out}}^3`$ cases. NPK concluded that models involving the mergers of BH-NS binaries and of NS-NS binaries are favored to explain short GRBs because then the flow is likely to be NDAF. However, our condition for the NDAF we found is more stringent that they found and may not be generally satisfied.

NPK also argue that the mergers of BH-WD binaries and BH-helium star binaries are difficult to explain GRBs for the following reason. Such mergers will create disks with large dimensions, thus being unable to reach the NDAF regimes but leading to the formation of CDAF. Then, convective motion mostly transports energy outward, thereby liberating little energy at small radii. It is thus unlikely to cause violent explosions like GRBs. However, we would like to point out that the flow pattern is rather uncertain in the presence of magnetic fields.

If there are strong magnetic fields with poloidal (vertical) components, for example, formation of outflow (or jet) is unavoidable (e.g. Blandford & Payne 1982; Kudoh, Matsumoto, & Shibata 1998), and such outflowing matter can produce enormous amount of synchrotron photons. Magnetic reconnections, giving rise to intense flare emission, may also occur. That is, mergers of BH and WD or helium star cannot be rejected as possible origins of GRBs. More careful discussion is needed regarding this issue.

### 3.2 Proof of $`p_{\mathrm{rad}}>p_{\mathrm{gas}}`$ in electron non-degeneracy

There are some discrepancies in the results between the present study and NPK. Those mainly come from a difference in the method. In NPK, for example, they analytically studied the scaling law of physical parameters in NDAF in two regimes, i.e., i) the gas pressure-dominated regime, and ii) the degeneracy pressure-dominated regime. Here, we show that gas pressure never becomes dominant whenever electrons are not degenerate at $`T\stackrel{>}{}\text{ }m_\mathrm{e}c^2/k_\mathrm{B}`$. Namely only when electrons are degenerate, either gas pressure or degeneracy pressure can become dominant, i.e., the case i) is not realized in electron non-degeneracy.

Using Eqs. (6) and (7) with the condition $`n_pn_n`$, we obtain the relation,

$`p_{\mathrm{rad}}/p_{\mathrm{gas}}13.02\left({\displaystyle \frac{\rho }{10^7\mathrm{g}\mathrm{cm}^3}}\right)^1\left({\displaystyle \frac{T}{\mathrm{MeV}}}\right)^3.`$ (72)
On the other hand, the non-degeneracy condition for relativistic electrons, see Eq. (12), is represented by

$`\left({\displaystyle \frac{\rho }{10^7\mathrm{g}\mathrm{cm}^3}}\right)^1\left({\displaystyle \frac{T}{\mathrm{MeV}}}\right)^30.1142.`$ (73)
From Eqs. (72) and (73), we see that the following relation is always realized,

$`p_{\mathrm{rad}}p_{\mathrm{gas}},`$ (74)
in the non-degenerate electron regime at $`T\stackrel{>}{}\text{ }m_\mathrm{e}c^2/k_\mathrm{B}`$.

### 3.3 Expected QED effects arising due to super-critical magnetic fields

When the field strength exceeds the critical value ($`𝒪(10^{13})`$ G), quantum electrodynamical processes become important. Then, as it were, magnetic fields themselves behave as “particles” in the strongly magnetized vacuum, thus producing several unique features. For example, strong magnetic fields induce large energy splitting of Landau levels, large refractive indices of photons, photon splitting effect, and so on (see, e.g., Adler 1971; Shabad 1975; Melrose & Stoneham 1976; Chistyakov et al. 1998; Kohri & Yamada 2002a). In particular, the photon splitting effect is interesting, when we consider the present model of the hyper-critical accretion disk to a central engine of GRBs. That is because if the rate of photon splitting ($`\gamma 2\gamma `$) is larger than that of $`e^\pm `$ pair production ($`\gamma +\gamma e^++e^{}`$), high energy photons with their energy $`E_\gamma `$MeV can lose kinetic energy rapidly and produce a lot of soft photons ($`E_\gamma \stackrel{<}{}\text{ }0.511`$MeV) without producing copious electron-positron pairs. Then, it may be possible to produce a fire ball as a source of GRB without requiring a large Lorentz factor.

Unfortunately, however, any concrete values of the photon splitting rate near the threshold energy of the pair production has not been calculated precisely except for some formal formulations (Adler 1971; Shabad 1975; Melrose & Stoneham 1976). The numerical study above the threshold is currently underway (Kohri & Yamada 2002b) and needs further attention.

In addition, radiative processes in a ultra-magnetized plasma are inherently complicated and we are still far away from having the basic physics under control. At the enormous densities we are dealing with the optical depth is huge, so radiative transfer must be solved. Therefore, the above predictions may be premature at present. Only a detailed radiative transfer calculation, which is outside the present capabilities of numerical codes, could address this point.

### 3.4 Detectability of neutrinos from NDAFs

As we have shown in the previous sections, we expect that a lot of neutrinos would be emitted from NDAF. It is thus of great importance to check whether the signals of the neutrinos are detectable or not in GRB events. Moreover, it is also informative to find some correlations between neutrino emission and photon emission from GRBs in terms of the observed times and/or directions to the source of the neutrinos.

Such investigation has already been initiated by Nagataki and Kohri (2001), who computed the neutrino luminosity and its detectability, properly considering the time evolution of the central BH. Their calculations are based not on the merger models but on the collapsar model, and they adopt the simple analytical fitting formula to the hyper-critical (slim) disk model taken from Fujimoto et al. (2001). (Basic numbers are in good agreement with those obtained by the numerical evaluations by Popham et al 1999.) Their conclusion is that we will be able to marginally detect the neutrinos from collapsars occurring with “TITAND,” a next-generation multi-megaton water Cherenkov detector (Suzuki 2001). The first phase of TITAND is planed to be with 2 Mt water inside. In addition, we will be also able to know the direction to the source of neutrinos in TITAND. For the total number of emitted $`\overline{\nu }_e`$, we roughly evaluate

$`N_{\overline{\nu }_e}`$ $``$ $`{\displaystyle \frac{\dot{q}_{Ne}V\mathrm{\Delta }t}{\overline{E}_{\overline{\nu }_e}}}`$ (75)
$``$ $`10^{57}\left({\displaystyle \frac{\dot{q}_{Ne}}{10^{35}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1}}\right)\left({\displaystyle \frac{\overline{E}_{\overline{\nu }_e}}{10\mathrm{M}\mathrm{e}\mathrm{V}}}\right)^1\left({\displaystyle \frac{V}{10^{20}\mathrm{cm}^3}}\right)\left({\displaystyle \frac{\mathrm{\Delta }t}{10^3\mathrm{s}}}\right),`$ (76)
with their mean energy $`\overline{E}_{\overline{\nu }_e}`$ to be emitted during one GRB event, where $`\dot{q}_{Ne}`$ is estimated at $`T10\mathrm{M}\mathrm{e}\mathrm{V}/k_\mathrm{B}`$ and $`\rho =10^{13}\mathrm{g}\mathrm{cm}^3`$, $`V`$ is the volume of the emitting region, and $`\mathrm{\Delta }t`$ is the duration. For the event number of $`\overline{\nu }_e`$’s,

$`R_{\mathrm{event}}10\left({\displaystyle \frac{d}{3\mathrm{Mpc}}}\right)^2\left({\displaystyle \frac{N_{\overline{\nu }_e}}{10^{57}}}\right)\left({\displaystyle \frac{\overline{E}_{\overline{\nu }_e}}{10\mathrm{MeV}}}\right)^2\left({\displaystyle \frac{V_{\mathrm{H}_2\mathrm{O}}}{2\mathrm{Mt}}}\right)`$ (77)
will be detected by TITAND, where $`d`$ is the distance from the source and $`V_{\mathrm{H}_2\mathrm{O}}`$ is the volume of water. Hence, we optimistically expect to detect neutrinos from NDAFs, if GRB occurs within a distance of 3 Mpc. If detected, neutrino signals and correlations with GRB will be able to make clear the unknown explosion mechanism of GRBs. In order to derive the reliable estimations on neutrino flux, however, we should perform detailed computations as was done by Nagataki and Kohri (2001). This is left as future work.

## 4 Conclusion

We study the properties of hyper-critical accretion flow and for surface density exceeding about $`10^{20}`$ g cm<sup>-2</sup>, which realizes when about a solar mass material is contained within a disk with a size of $`5\times 10^6`$ cm, we find the following unique features:

1. Radiation luminosity of such flow is practically zero due to significant photon trapping, although mass accretion rate enormously exceeds the critical rate, $`\dot{M}L_{\mathrm{Edd}}/c^2`$.
2. Neutrino-cooling dominates over advective cooling. Thus the flow can cool via neutrino emission. We expect to detect neutrinos from GRBs with the next-generation, multi-megaton water Cherenkov detector in future.
3. Electron degeneracy pressure dominates over gas and radiation pressure, and the degenerate electron definitely influences the processes of neutrino emission. This feature is distinct from that of the solution found by NPK.
4. The disks are stable both against thermal and viscous instabilities, even when we take into account electron degeneracy.
5. Magnetic field strength exceeds the critical value, even if we take only 0.1% of the equi-partition value. Then, photon splitting may occur, producing significant number of photons of energy below $`511`$ keV, thereby possibly suppressing e<sup>±</sup> pair creation. However, quantitative discussion is left as future work.

We are grateful to R. Narayan, Shoichi Yamada and S. Nagataki for useful discussions and comments. The comments by the anonymous referee are also suggestive and helpful in making the final draft. We also thank the Yukawa Institute for Theoretical Physics at Kyoto University for the YITP workshops YITP-W-01-13? on “Gamma-Ray Bursts,” and YITP-W-01-17 on “Black Holes, Gravitational Lens, and Gamma-Ray Bursts,” in which this work was initiated and completed. This work was supported in part by the Grants-in Aid of the Ministry of Education, Science, Sports, and Culture of Japan (13640238, SM). Numerical computation in this work was carried out at the Yukawa Institute Computer Facility.


format: text

text:
# Time series analysis for minority game simulations of financial markets

## 1 Introduction

The pricing of contingent claims contracts in financial economics is often based on very restrictive assumptions about the time evolution of the underlying instrument . In recent years researchers have endeavored to remove some of these restrictions by proposing more realistic models which would incorporate features found in real markets . There is however a trade-off between analytic tractability and adherence to stylized facts observed from empirical financial time series. The most attractive features of the usual Black-Scholes type of models is the possibility of obtaining closed, exact formulas for the premium of derivative securities, and to build a risk free replicating strategy. Such qualities are amply used in financial institutions which require fast calculation and tools to hedge risky assets.

However, the geometric Brownian motion assumption of the Black-Scholes models ignores several empirically observed features of the real markets such as, volatility clusters, fat tails, scaling, occurrence of crashes, etc. It is as yet unknown which stochastic process is responsible for the motion of risky assets, but physicists have taken some important steps in the right direction. In this work we implement a microscopic agents-based model of market dynamics which gives rise to a quite complex and rich behaviour and whose output are macroscopic quantities such as price returns. By varying its parameters the model exhibits market crashes, Gaussian statistics and short ranged correlation, fat tailed returns and long range correlation. This model retains the nontrivial opinion formation structure of the grand canonical minority game because it incorporates two new features. The first one is to allow two categories of agents, producers (who use the market for exchanging goods), and speculators (whose aim is to profit from price fluctuations). The second feature is that speculators might choose not to trade, and in this sense the model is similar to the grand canonical ensemble of statistical physics since the number of active traders is not constant.

We perform an analysis of the time series generated by this model in order to classify its dynamical behaviour. We first test the data for unit roots and remove a simple kind of nonstationarity by taking differences wherever necessary. The BDS statistic , originated in the chaos literature, uses the correlation integral  as the basis of a test for the hypothesis that the data is independent and identically distributed (i.i.d.). We apply this statistic as a model specification test by applying it to the residue of an ARIMA, autoregressive integrated moving average, process (although this might not remove all kinds of linearity) . If the null hypothesis is rejected then this is indication that the data is nonlinear. Another test based on surrogate data  is used to confirm the nonlinearity of the model. Since the alternative hypothesis for the BDS procedure is not specified other tests have to be applied in order to determine whether the nonlinearity comes from a stochastic or deterministic mechanism. Such distinction is subtle and here we approach this question by computing more parameters, the correlation dimension and two other procedures which do not require embeddings: recurrence plots , and the Lempel-Ziv complexity (LZC) . In complement to the above procedures we implement a Bayesian approach called cluster weighted modelling (CWM)  in order to find further indication of determinism.

The paper is organized as follows. The market model and its trajectories are analysed in Section 2. The time series analysis, including all statistical tests is discussed in Section 3. A succinct presentation of the CWM and its results are found in Section 4. Analysis of the results obtained and a classification of the model evolution are given in Section 5 while in the last section we summarize the main results and comment on future work.

## 2 Market Model

Inspired by the El Farol Problem proposed by Arthur , the so called Minority Game (MG) model introduced by Challet and Zhang  represents a fascinating toy-model for financial market. Now it is becoming a paradigm for complex adaptive systems, in which individual members or traders repeatedly compete to be in the winning group. The game consists of $`N`$ agents that participate in the market buying or selling some kind of asset, e.g. stocks. At any given time agent $`i`$ can take two possible actions $`a_i=\pm 1`$, meaning buy or sell. Those players whose bets fall in the minority group are the winners, i.e., the sellers win if there is an excess of buyers, and vice versa. To determine the minority group we just consider the sign of the global action $`A(t)=_ia_i`$, so that if $`A(t)`$ is positive the minority is the group of sellers; in this case the majority of players expect asset prices to go up.

In other words, this dynamics follows the law of demand and supply. In the end of each turn, the output is a single label, $`1`$ or $`1`$, denoting the winning group at each time step. This output is made available to all traders, and it is the only information they can use to make decisions in subsequent turns. Indeed, they store the $`m`$ most recent output of winners set $`\mu (t)`$. In this way a limited memory of length $`m`$ is assigned to the traders corresponding to the most recent history bit-string that traders use to make decisions for the next step. In order to decide what action to take, agents use strategies. A strategy is an object that processes the outcomes of the winning sets in the last $`m`$ bets and from this information it determines whether a given agent should buy or sell for the next turn. When a tie results from the strategies, buying or selling is decided by coin tossing. The memory defines $`D=2^m`$ possible past histories so the strategy of one agent can be viewed as a D-dimensional vector, whose elements can be 1 or -1. The space $`\mathrm{\Gamma }`$ of strategies is an hypercube of dimension $`D`$, and the total number of strategies in this space is $`2^D`$. At the beginning of the game each agent draws randomly a number $`S`$ of strategies from the space $`\mathrm{\Gamma }`$ and keeps them forever. After each turn, the traders assigns one (virtual) point to each of the strategies which would have predicted the correct outcome. Along the game the traders always choose the strategy with the highest score. The MG is a very simple model, capable of exhibiting complex behaviour, like phase transitions between an information-efficient phase and information-inefficient phase, by just varying a control parameter $`\alpha =\frac{2.2^m}{NS}`$, e.g., the ratio between information complexity and number of strategies present in the game.

More realistic features are reached with the grand canonical minority game. In this model we define the price process in terms of excess demand $`A(t)`$,

$$\mathrm{log}p(t+1)=\mathrm{log}p(t)+A(t),$$
(1)
and introduce two kinds of agents.

* The first kind is called producers, who go to the market only for the purpose of exchanging goods; they only have one strategy and in this sense they behave in a deterministic way with respect to $`\mu (t)`$. The number of producers is $`N_p`$.
* The other kind of agents are the speculators, who go to the market to make profit from price fluctuations. Since they are endowed with at least two strategies during the game they need to use the best strategy; in this sense speculators are adaptative agents with bounded rationality. The number of speculators is $`N_s`$.

In this version of the game, the number of speculators can change anytime, because the agent may decide not to trade in which case $`a_i=0`$. Since the strategy chosen by speculators is the one with highest score, it is very important to update the scores of each strategy at each time step. The updating of the scoring $`U_{i,s}`$ of each strategy $`s`$ belonging to agent $`i`$ is given by the following equations

$`U_{i,s}(t+1)`$ $`=`$ $`U_{i,s}(t)a_{i,s}^{\mu (t)}A(t)\text{for}s>0`$
$`U_{i,0}(t+1)`$ $`=`$ $`U_{i,0}(t)+ϵ`$ (2)
where $`ϵ`$ is a threshold parameter. A sample of the trajectory generated by this model is shown in Fig. 1(a). An important piece of information to understand the generating mechanism of the system can be obtained from the periods where no coin is tossed. The game is necessarily deterministic during this mode of evolution and we selected a rather long period of 1114 points to perform a time series analysis in this regime.

This work is based upon the grand canonical MG model, henceforth called MG model.

## 3 Time Series Analysis

In this section we discuss the tools employed in time series analysis starting with the BDS statistics. We generate $`10000`$ points for each of the benchmark data described in the Appendix. Another benchmark is the set comprising $`9998`$ closing prices for the index SP500 from January 01, 1965 to January 01, 1995, shown in Fig. 1(b). All time series used here are unit root stationary with the exception of the financial index: this nonstationary is removed by taking first difference.

The BDS test uses the correlation integral  as the basis of a statistic to test whether a series of data is i.i.d. In the chaos literature the correlation integral is part of an efficient tool to compute the fractal dimension of objects called attractors (for a formal definition of attractor see, e.g., Ref. ). Given a sample of empirical data $`\{x_t\}_{t=1}^N`$, the theory of state-space reconstruction  requires that the $`d`$-histories of the data be constructed, where $`d`$ is called the embedding dimension. Under certain conditions it is possible to reproduce in this space the dynamics of the system for a correct choice of $`d`$. The correlation integral is a function defined on the trajectories in this space and from it one can compute the correlation dimension. A simple test for determinism consists of increasing the embedding dimension and observing the occurrence of a corresponding increase in the correlation dimension. Some conditions have to be met in order to apply this method, mainly stationarity and sufficient number of data points . Our results show that the increase of the correlation dimension for the MG model with respect to the embedding dimension is practically identical to the stochastic benchmark series, including the index SP500. In this work we resort to better ways of analysing the occurrence of stochastic behaviour in complex time series evolution.

The BDS statistic will be applied to the residue of ARIMA processes in order to detect nonlinearity in the data. In this sense the statistic is used as a specification test. The asymptotic distribution of the statistic under the null of pure whiteness, is the standard normal distribution. The alternative hypothesis is not specified . The code implemented here is taken from Ref. . From Table 1, the null hypothesis for the MG model and for the SP500 index are rejected more strongly than for the nonlinear stochastic model $`NLMA(2)`$.

In the surrogate data analysis, a null hypothesis is tested under a measure $`M`$, usually some nonlinear statistic. Surrogates are copies of the original time series preserving all its linear stochastic structure. Let $`M_{\text{orig}}`$ and $`M_{\text{surr}}`$ be the test measure for the original and the surrogate data sets, respectively. The null hypothesis is rejected only if $`M_{\text{orig}}<M_{\text{surr}}`$ or $`M_{\text{orig}}>M_{\text{surr}}`$ for all surrogates. Some care should be taken to test the null hypothesis, otherwise false-positive rejection of null hypothesis could result. A detailed analysis on this is found in .

Here surrogate data sets  are used in complement to the BDS statistic to test for nonlinearity. The null hypothesis is that the data is described by a stationary linear stochastic process with Gaussian inputs. The test statistic $`M`$ is a simple measure of predictability as in  and the procedure is to generate $`19`$ samples of rescaled surrogates and to compare with the MG model. If the test statistic falls outside the interval defined by the surrogates, then the null is rejected at the $`95`$% level. We found that, for embedding dimensions from 2 to 5, the MG is nonlinear at this level of significance, corroborating the previous result using the BDS.

In the remaining of this Section we discuss methods which do not require embeddings.

Another measure of randomness that provides further insight into time series dynamics is the Lempel-Ziv complexity . No embedding is necessary and the data is interpreted as a binary signal generated by some kind of source. This idea is ever present in communication theory where one wishes to determine the minimum alphabet required to code a source whose signal is to be sent through a noisy channel. Let us consider the length $`L(N)`$ of the minimal program that reproduces a sequence with $`N`$ symbols. The Lempel-Ziv algorithm is constructed by a special parsing which splits the sequence into words of the shortest length and that has not appeared previously. For example, the sequence $`0011101001011011`$ is parsed as $`\mathrm{0.01.1.10.100.101.1011}`$. One can show that $`L(N)N_w(\mathrm{log}N_w(N)+1)`$ where $`N_w`$ is the number of distinct words in a parsing and $`N`$ the size of the sequence. From this one can see that $`L(N)`$ contains a measure of randomness where a source that produces a greater number of new words is more random than a source producing a more repetitive pattern. In analogy with dynamical evolution, those systems that are composed of well defined cycles are predictable while chaotic motion and stochastic processes are always producing new kinds of trajectories that never repeat themselves. A comparison between chaos generated by differential equations and stochasticity can now be obtained. In the former case the Lempel-Ziv complexity is well below 1 while in the latter it is close to this value. More specifically, if we consider an oscillatory system such as the well known van der Pol oscillator, then $`C=0.049`$. For the Lorenz attractor with 2000 points $`C=0.181`$. In Section 5 we comment on the discrepancy between this value and that in Table 2, and discuss the use of this complexity measure in dynamical systems generated by maps. The complexity for the MG dynamics is found to be $`0.82`$ while the financial index has $`C1`$.

The implementation of recurrence plots  used here is taken from  since it provides a clear distinction of the systems we intend to classify. The idea is to detect regions of “close returns” in a data set. The construction of the plots is simple: just compute the absolute values $`x_ix_n`$ for all the points in the data base. If the horizontal axis is designated by $`i`$, corresponding to $`x_i`$, and the vertical axis by $`i+n`$, corresponding to $`x_{i+n}`$, then plot a black dot at the site $`(i,n)`$, whenever the absolute value difference is lower than $`\delta `$; otherwise plot it white. Actually the time series is normalized to $`(0,1)`$ and then $`\delta `$ is taken as $`26`$% of the average distance between successive points. The black/white pattern can be used to detect determinism in the data. There is a clear difference amongst plots generated by differential equations, maps, random data and stochastic processes as shown in Fig. 2. Patterns of horizontal segments in the recurrence plot indicate the presence of unstable periodic orbits in maps or differential equations. In this sense, these plots detect low dimensional chaos in relatively small and even noisy data sets.

## 4 Cluster-Weighted Modelling

An interesting probability density estimation approach to characterize and forecast time series developed by Gershenfeld, Schoner and Metois  is the so called cluster-weighted modelling. This seems to be a powerful technique as it characterizes extremely well the time series of nonlinear, nonstationary, non-Gaussian and discontinuous systems using probabilistic dependence of local models. The cluster-weighted modelling technique estimates the functional dependence of time series in terms of delay coordinates. The main task of this approach is to find the conditional forecast by estimating the joint probability density.

Let $`\{y_n,\stackrel{}{x}_n\}_{n=1}^N`$ be the $`N`$ observations in which $`\stackrel{}{x}_n`$ are known inputs and $`y_n`$ are the corresponding outputs. By knowing the joint probability density $`p(y,\stackrel{}{x})`$, we can derive the conditional forecast, the expectation value of $`y`$ given $`\stackrel{}{x}`$, $`y|\stackrel{}{x}`$. We can also deduce other quantities such as the variance of the above estimation. Actually, the joint density $`p(y,\stackrel{}{x})`$ is expanded terms of clusters which describe the local models. Each cluster contains three terms namely, the weight $`p(c_m)`$, the domain of influence in the input space $`p(\stackrel{}{x}|c_m)`$, and finally the dependence in the output space $`p(y|\stackrel{}{x},c_m)`$. Thus the joint density can be written as ,

$`p(y,\stackrel{}{x})={\displaystyle \underset{m=1}{\overset{M}{}}}p(y|\stackrel{}{x},c_m)p(\stackrel{}{x}|c_m)p(c_m)`$ (3)
Once the joint density is known the other quantities can be derived from $`p(y,\stackrel{}{x})`$. For example, the conditional forecast is given by,

$`y|\stackrel{}{x}={\displaystyle yp(y|\stackrel{}{x})𝑑y}={\displaystyle \frac{_{m=1}^Mf(\stackrel{}{x},\beta _m)p(\stackrel{}{x}|c_m)p(c_m)}{_{m=1}^Mp(\stackrel{}{x}|c_m)p(c_m)}}.`$ (4)
Here $`f(\stackrel{}{x},\beta _m)`$ describes the local relationship between $`\stackrel{}{x}`$ and $`y`$. The parameters $`\beta _m`$ are found by maximizing the cluster-weighted log-likelihood. The simplest approximation for the local model is with linear coefficients of the form,

$`f(\stackrel{}{x},\beta _m)={\displaystyle \underset{i=1}{\overset{I}{}}}\beta _{m,i}f_i(\stackrel{}{x}).`$ (5)
The method just described is capable of modelling a wide range of deterministic time series. Here we use cluster weighted modelling to distinguish between deterministic and stochastic time series. In deterministic systems we observe that the variances of the different clusters converge to values lower than the variance of the original time series and one can verify that this property is robust under changes of the number of clusters. On the other hand stochastic systems do not have this property. Fig. 3 illustrates the comparison between the fitting of a deterministic system (Lorenz) and the minority game data using cluster weighted modelling.

## 5 Analysis and Classification

The main objective is to understand the minority game mode of evolution and other similar time series behaviour. Although this system is not generated by any kind of differentiable dynamics or even stochastic differential equations, we use in its analysis methods from dynamical systems, stochastic processes and complex systems theory. Tables 1 & 2 and Fig. 2 summarize the main results. We will make frequent reference to the SP500 index since the minority game model is supposed to reproduce the dynamical evolution of financial markets and this index is used as a benchmark for comparison.

The BDS statistic and surrogate data analysis provide clear evidence that the minority game is nonlinear. The benchmark systems described in the Appendix reproduce the expected results for the BDS in Table 2 and we can clearly see that the null hypothesis is strongly rejected for known deterministic nonlinear systems. As for the nonlinear stochastic system, the index and the MG model, the rejection occurs at all dimensions and the probability of type I error is practically zero. In particular, the BDS was used in other instances of the MG model, using different parameters and time intervals with the same result. These findings support the conclusion that there is a nonlinear mechanism in operation which drives the MG dynamics and that this property is robust to the extent tested herein.

The Lempel-Ziv complexity is an important parameter that can be used in the analysis of complex systems. Its advantage is that it does not require embeddings and can be easily employed in conjunction with other methods. The results in Table 1 supports the idea that there is a stochastic mechanism in operation driving the MG model. There is a higher degree of indeterminacy in the SP500 and this is perhaps due to the fact that in this index there is a certain amount of measurement noise. The surprisingly low complexity of the Lorenz system is comparable to that of limit cycles, e.g. Van der Pol oscillator. The explanation for this comes about when we compute its complexity for shorter time series. For example at 2000 points the complexity is about $`2.7`$ times higher than the complexity for a 10000 length series as reported in Table 2. This phenomenon does not occur for the discrete system like Henon attractor. Due to the fact that the Lorenz attractor contains a dense set of unstable periodic orbits, long time evolution affects the computation of the complexity and reveals some resemblance with periodic systems. Effects of this magnitude did not appear in the other time series analysed. In the deterministic intervals of the MG model, the complexity has a value comparable to that of the chaotic Henon map.

The Lempel-Ziv complexity, as any other test or statistic, should always be used in conjunction with other diagnostic tools. In all simulations performed so far we have never found a stochastic process with complexity less than 0.8. However there are chaotic systems with complexity beyond this value, for example the family of maps $`x_{n+1}=px_n\text{mod}\mathrm{\hspace{0.17em}1}`$, $`p`$ a prime number. Another result is that known chaotic systems described by differential equations do not have high complexity. To conclude that a system is stochastic we employ recurrence plots and the cluster weighted modelling approach. The complexity is then used as a confirmation and, more importantly, to associate a level of stochasticity. This is done in the same way as using Lyapunov exponents to quantify chaos once determinism has been found.

The recurrence plot is a visual method which helps in the identification of similarities and differences amongst diverse modes of evolution. Several tests with differential equations and maps, represented here by the Lorenz and Henon systems, show that it is unlikely that low dimensional chaos can describe the kind of evolution found in the markets and in the MG model. In particular, recurrences are clearly identified in the Lorenz system and nothing of the kind will ever appear in stochastic models. In this sense stochastic processes are better suited to describe the financial index and the MG model. However, in the deterministic regime one can identify close return patterns similar to those of low dimensional chaotic systems.

Cluster weighted analysis reveals another aspect of the minority game behaviour. Its use in the modelling of general deterministic systems produces clusters whose variances are always smaller than the variance of the data. In contrast to this, when applied to a stochastic system the variances of the clusters are comparable to the variances of the data. Such distinction is preserved when we vary the number of clusters.

## 6 Conclusions and Further Work

The issue of nonstationarity is a subtle one. We limited ourselves in this study to unit root stationarity, but more sophisticated methods need to be used for the several brands of MG models and financial indices. Using the complexity parameter we found that long time evolution reveals some intrinsic features of chaotic attractors described by differential equations. Also, this parameter associated a higher complexity to the SP500 index as compared to the MG model and a possible explanation for this was given above. The recurrence plots confirm that, in general, the MG model cannot be described by low dimensional chaotic systems. The nonlinear character of the model and the index are clearly indicated by the BDS test and surrogate data analysis.

The MG model parameters employed in our simulations were chosen in the information efficient region. In the inefficient case our simulations show that the system is a non random process with complexity close to zero. The recurrence plot in this case provides a strong evidence of recurring orbits with several traces of horizontal segments. Thus, randomness is directly related to efficiency.

An important finding in our analysis of the MG game is that low dimensional chaotic regimes are possible during the evolution. This occurs for intervals of about 1000 iterations, while shorter intervals of deterministic evolution are difficult to classify. A complete study of the probabilistic structure of the deterministic phases, and their statistical significance, requires a more extensive investigation . For the full motion, comprised of transitions between deterministic modes induced by random perturbations, the time series analysis has indicated the operation of a nonlinear stochastic process which is similar, but not identical, to the SP500 index. In real markets it is possible to find periods of deterministic behaviour in exchange rates for certain instruments and specific periods of time . In the SP500 such occurrence is unlikely but an extensive search in the data base for several lengths and starting points is feasible. Even if deterministic modes can be found in historical stock prices, indices or rates, turning this into potential profit is a remote, but not an entirely discarded, possibility worthwhile pursuing.

Another interest in this investigation is that it provides a test amongst several models. In this sense one selects the models whose time series properties reproduce better real market evolution. In derivatives pricing a widely used test, and parameter estimation procedure, is to minimize the hedging error. This can be implemented only if a model is available, e.g. Black-Scholes. In a non-gaussian context there is no consensus on which model is appropriate. An extension of this work, which will be pursued in a future work, is to use time series methods in the several brands of agents based models and to compare them with market data. In addition to the stylized facts we should include nonlinearity, at least for the SP500, and this is an additional reason to discard simple models based on the geometric Brownian motion.

Additional investigations in this study is to model the trajectories of the Minority Game model having in view the pricing of derivatives instruments, and to use a more extensive set of statistics to examine other financial series in order to confront the similarities and differences of this model with real market data.

The authors acknowledge Damien Challet for useful e-mail on the implementation of the grand canonical MG model. We acknowledge Fundação de Amparo à Pesquisa do Estado de São Paulo of Brazil for fully supporting this research (F.F.F. and P.M.) and Coordenadoria de Aperfeiçoamento de Pessoal de Ensino Superior (B.S.M.)

## Appendix A Description of the data sets

The time series used as benchmarks were chosen to represent the kind of behaviour we intend to identify in the evolution of the minority game model. The Lorenz system and the Henon mapping are prototypes of deterministic behavior generated by differential equations and differentiable mappings. ARMA models and the NLMA are examples of linear and nonlinear stochastic processes. In the following we describe briefly the models used in the present study.

### A.1 Lorenz system:

Nonlinear differential equations

$`{\displaystyle \frac{dx}{dt}}=\sigma (yx);{\displaystyle \frac{dy}{dt}}=xz+rxy;{\displaystyle \frac{dz}{dt}}=xy\beta z`$ (6)
The parameter values are chosen as $`\sigma =16`$, $`\beta =4`$ and $`r=40`$. The data series is obtained by solving Eqs. (6) numerically using fourth order Runge-Kutta method.

### A.2 Henon map:

Nonlinear differentiable mapping

$`x_t=1ax_{t1}^2+y_{t1};y_t=bx_{t1}`$ (7)
where $`a=1.4`$ and $`b=0.3`$.

### A.3 Random (Noise):

The data set is taken by drawing random points from a uniform distribution in the unit interval $`x_i\text{rand}(0,1)`$.

### A.4 ARMA(2,1):

Linear stochastic process

$`x_t=0.8x_{t1}0.5x_{t2}+0.4ϵ_{t1}+ϵ_t`$ (8)
where $`ϵ_ti.i.d.`$, a normal distribution with zero mean and unit variance $`N(0,1)`$.

### A.5 NLMA(2):

Nonlinear stochastic process

$`y_t=ϵ_t0.3ϵ_{t1}+0.2ϵ_{t2}+0.4ϵ_{t1}ϵ_{t2}0.25ϵ_{t2}^2`$ (9)
where $`ϵ_ti.i.d.`$ $`N(0,1)`$.


text:
# Fluxes of cosmic rays: A delicately balanced stationary state

## Abstract

The analysis of cosmic rays fluxes as a function of energy reveals a knee slightly below $`10^{16}`$ eV and an ankle close to $`10^{19}`$ eV. Their physical origins remain up to now quite enigmatic; in particular, no elementary process is known which occurs at energies close to $`10^{16}`$ eV. We propose a phenomenological approach along the lines of nonextensive statistical mechanics, a formalism which contains Boltzmann-Gibbs statistical mechanics as a particular case. The knee then appears as a crossover between two fractal-like thermal regimes, the crossover being caused by process occurring at energies ten million times lower than that of the knee, in the region of the quark hadron transition ($`10^9`$ eV). This opens the door to an unexpected standpoint for further clarifying the phenomenon.

PACS number: 95.30.-k, 95.35.+d, 95.85.Ry, 05.90.+m

Cosmic rays fascinate since long. They provide galactic, extragalactic and cosmological information, related to recent or very old events concerning various sources, going back to the early times of the universe . They reflect all types of elementary process and interactions, and are associated with phenomena of very different space and time scales. The complete physical scenario still remains quite enigmatic, although a variety of specific mechanisms for acceleration and propagation have been advanced along the years for various energy regions. The most known of these mechanisms is the Fermi one , which addresses acceleration in magnetized turbulent plasma, but many others have been advanced  in relation with the knee and energies below it; for the energies beyond these, for the ankle, as well as for general reviews, see .

Through various types of detectors, the flux of cosmic rays at the top of the Earth atmosphere has been measured  and varies from $`10^4`$ down to $`10^{29}[m^2srs\text{GeV}]^1`$ for energies increasing from $`10^8`$ up to near $`10^{21}`$ eV: See Fig. 1. This distribution (which spans 13 decades in energy and 33 decades in flux!) is not exponential, hence it does not correspond to Boltzmann-Gibbs (BG) statistical mechanics thermal equilibrium. Consistently, even at a phenomenological level, i.e., without specifying any concrete model or mechanism, this problem represents a challenge. This is the one we address here. We shall use a point of view based on a current generalization of Boltzmann-Gibbs statistical mechanics, referred to as nonextensive statistical mechanics, we shall briefly describe later on. The first step will be to remark that the fluxes of cosmic rays in general, and the studies of the“knee” and the “ankle” in particular, involve phenomena such as turbulence (see, for instance, ), anomalous diffusion and fractality (see, for instance, ), self-organized criticality (see, for instance, ), long-range interactions (classical and quantum gravitation), among other complex phenomena (such as, for example, possible nonmarkovianity ). It is precisely such phenomena that constitute the scope of nonextensive statistical mechanics; for turbulence and related matters see, for instance, ; for anomalous diffusion and fractality see, for instance, ; for self-organized criticality see, for instance, ; for long-range interactions see, for instance, .

Let us now first briefly review the usual, BG thermostatistics. If we optimize under appropriate constraints the BG $`S=k_ip_i\mathrm{ln}p_i`$ ($`k`$ Boltzmann constant; $`\{p_i\}`$ microscopic probabilities) we obtain the celebrated equilibrium distribution $`p_i=\frac{e^{\beta E_i}}{Z}e^{\beta E_i}`$ ($`\beta 1/kT`$, $`E_i`$ energy of the i-th state; $`Z_je^{\beta E_j}`$ partition function). Excepting for the trivial normalizing factor $`1/Z`$, this distribution can alternatively be obtained as the solution of the linear differential equation

$$dp_i/dE_i=\beta p_i.$$
(1)
In order to deal with a variety of thermodynamically anomalous systems, a more general formalism, nonextensive statistical mechanics, was introduced in 1988 . It is based on the generalized entropic form $`S_q=k(1_ip_i^q)/(q1)`$ ($`q`$ and $`S_1=S`$). Its optimization under appropriate constraints yields  a power-law, $`p_i[1(1q)\beta _qE_i]^{\frac{1}{1q}}e_q^{\beta _qE_i}`$ (definition), which recovers the BG weight for $`q=1`$ ($`\beta _1\beta `$). As usual, $`kT_q1/\beta _q`$ characterizes the conveniently averaged energy. This anomalous equilibrium-like distribution can be alternatively obtained (excepting for the normalizing factor) by solving the nonlinear differential equation

$$dp_i/dE_i=\beta _qp_i^q.$$
(2)
This generalized weight naturally emerges in ubiquitous problems such as fully developed turbulence  (which is relevant for the Fermi mechanism, and most probably for others as well), electron-positron annihilation , motion of Hydra viridissima , long-range many-body Hamiltonians , among many others.

We now use the above differential equation path in order to further generalize the anomalous equilibrium distribution, in such a way as to have a crossover from anomalous ($`q1`$) to normal ($`q=1`$) thermostatistics, while increasing the energy. We consider then the differential equation

$$dp_i/dE_i=\beta _1p_i(\beta _q\beta _1)p_i^q,$$
(3)
whose solution is $`p_i[1\frac{\beta _q}{\beta _1}+\frac{\beta _q}{\beta _1}e^{(q1)\beta _1E_i}]^{\frac{1}{q1}}`$. The crossover typically occurs for $`q>1`$ and $`\beta _1\beta _q`$, the distribution being anomalous at low energies and BG at high energies. It is undoubtedly interesting to notice that this differential equation precisely coincides, for $`q=2`$, with the heuristic one that in 1900 led Planck to the discovery of the black-body radiation law and ultimately to quantum mechanics .

Finally, by doing one more step along the same direction, we can further generalize the differential equation, now becoming

$$dp_i/dE_i=\beta _q^{}p_i^q^{}(\beta _q\beta _q^{})p_i^q.$$
(4)
This manner of writing the coefficient of the $`p_i^q`$ term has the advantage of recovering the simplest generalization of the BG distribution by considering $`\beta _q^{}=\beta _q`$ or $`\beta _q^{}=0`$ or even $`q^{}=q`$. For $`1<q^{}<q`$ and $`\beta _q^{}\beta _q`$, a crossover occurs at

$$E_{crossover}=[(q1)\beta _q]^{\frac{q^{}1}{qq^{}}}/[(q^{}1)\beta _q^{}]^{\frac{q1}{qq^{}}}.$$
(5)
For $`EE_{crossover}`$ we have an anomalous distribution characterized by $`(q,\beta _q)`$ (namely $`p_ie_q^{\beta _qE_i}`$), whereas for $`EE_{crossover}`$ we have a different anomalous distribution characterized by $`(q^{},\beta _q^{})`$ (namely $`p_ie_q^{}^{\beta _q^{}E_i}`$). The exact solution of the above differential equation (the most general one considered here) is given by $`p_if(E_i)`$ where $`f^1(x)`$ is an explicit monotonic function of $`x`$ involving hypergeometric functions (see Ref.  for details). Interestingly enough, this precise solution arrives in the discussion of the re-association of $`CO`$ molecules in Myoglobin , where time plays a role very analogous to the one played by energy in our cosmic rays problem. This time-energy analogy is not surprising after all if we take into account that, in the history of the universe after the big-bang, the time scale reflects the energy scale, as discussed in detail in Ref. .

The flux $`\mathrm{\Phi }(E)`$ can be obtained straightforwardly from $`p_if(E_i)`$ by calculating the density of states $`\omega (E)`$. In the ultrarelativistic limit $`E|𝐩|`$ ($`𝐩`$ being the momentum), which we adopt here for simplicity given the high values of the involved energies, the density of states of an ideal gas in three dimensions is given by $`\omega (E)E^2`$, hence $`\mathrm{\Phi }(E)=AE^2f(E)`$, where $`A`$ is a normalizing factor (and where red shift effects have been neglected). With this expression we fit the observational data and obtain the results displayed in Fig. 1. As we can see, the agreement is quite remarkable.

Our summarizing comments are:

(i) The high quality agreement over so many decades, including crossovers between different regimes, suggests that the phenomenological approach is correct, and specific models clarifying the various physical mechanisms that are involved should essentially satisfy it;

(ii) The deep explanation of the knee might well be found at energies extremely lower (ten million times lower, in fact), basically at energies related to the characteristic temperatures obtained from the fitting, namely $`\mathrm{9.615\hspace{0.33em}10}^7`$ eV (energy comparable to the pion mass) and $`\mathrm{1.562\hspace{0.33em}10}^9`$ eV (energy corresponding to the quark-hadron transition , as well as to the proton mass), and given the fact that $`q`$ and $`q^{}`$ differ by only $`3\%`$ (which in log-log representation asymptotically determines two almost parallel straight lines which very slowly approach to each other and intersect at $`E=E_{crossover}`$). The existence of two thermostatistical regimes, respectively related to ($`q,\beta _q`$) and ($`q^{},\beta _q^{}`$), could correspond to two different mechanisms of acceleration/propagation, for instance related to galactic and extra-galactic contributions. It is worthy stressing at this point that the location of the knee emerges here through Eq. (5), i.e., using only the basic four phenomenological parameters $`(q,\beta _q,q^{},\beta _q^{})`$, and not by introducing an extra parameter on top of the previous ones (whose role would be to fix the knee at the right value);

(iii) Since the entropic index $`q`$ is known to reflect (multi) fractality , the present results strongly suggest that either the generation or the transport (or both) of cosmic rays occur in scale invariant media, which is consistent with Ref. ;

(iv) With the better observational statistics (i.e., with refined precision) at these very high energies, expected from the Pierre Auger Observatory (or from similar projects), it might happen that the ankle disappears; if it does not, then it is probably associated with strongly nonstationary phenomena (perhaps related to memory effects from early cosmic stages), certainly out from the present thermodynamical description — such study would eventually clarify the physical meaning, if any, of the GZK feature  (see also ).

(v) The Kaskade collaboration  has shown the relevance of the various constituents which compose cosmic rays. In order to take this into account, the present approach could be further improved by using the grand-canonical instead of the canonical ensemble used here.

(vi) The average energy $`E_0^{\mathrm{}}𝑑EE\mathrm{\Phi }(E)/_0^{\mathrm{}}𝑑E\mathrm{\Phi }(E)`$ has been calculated, within the present approach, to be $`E2.489\text{GeV}`$. Any connection of this value with other cosmological or astrophysical quantities is of course very welcome.

Stimulating remarks from I. Bediaga, L. Galgani, R.C. Shellard and L. Masperi are acknowledged. Partial support from CNPq, PRONEX, CAPES and FAPERJ (Brazilian agencies) is also acknowledged.