text:
# Accretion Models of Gamma-Ray Bursts

## 1  Introduction

The fireball model (see Piran, 1999, 2000, for reviews) provides a good understanding of conditions within the $`\gamma `$-ray-emitting and afterglow-emitting regions of gamma-ray bursts (GRBs). According to this model, GRBs are produced when relativistic ejecta from a “central engine” are slowed down by interactions, either with an external medium (the external shock model) or among different layers within the ejecta themselves (the internal shock model). In the interactions the kinetic energy in the ejecta is converted to relativistic electrons which produce the observed radiation. Among the many successes of the model we note the observational confirmation of relativistic motion in the afterglow (Frail et al., 1997; Katz & Piran, 1997).

Despite the successes of the fireball model, the nature of the central engine remains a mystery. The problem is that the central engine is hidden from view; no radiation (apart from gravitational radiation and neutrinos that may possibly be detected in the distant future) reaches the observer directly from the engine. For a number of so-called “long bursts” accurate positions have been determined through observations of their afterglows. Based on this, there is circumstantial evidence that these bursts are associated with star-forming regions (e.g. Bloom et al. 2000). There is no information at present on the other class of GRBs, the so-called “short bursts.”

Although we lack direct evidence on the nature of the central engine, it is nevertheless widely accepted that GRBs are the result of cataclysmic events involving either neutron stars or stellar-mass black holes. The arguments in support of this hypothesis are straightforward. (i) Since bursts radiate the bulk of their energy in the $`\gamma `$-ray band, it seems likely that a relativistic object is behind their production. (ii) The energy budget ($`10^{51}`$ erg) is comparable to the kinetic energy of ejecta in a supernova explosion. (iii) Most long bursts are highly variable in gamma-rays, and so are many short bursts (Nakar & Piran, 2001a). In particular, the ratio of the total duration of the burst to the variability time scale is large, from which one concludes that the gamma-rays must be produced in internal shocks (Sari & Piran,, 1997). A key feature of internal shocks is that the observed gamma-ray variability reflects the variability in the activity of the central engine (Sari & Piran,, 1997). Since variability time scales as short as a millisecond are observed, the engine must contain a compact object of no more than a few solar masses (otherwise the light-crossing time would exceed the variability time).

An interesting clue to the nature of GRBs is provided by the durations of bursts. While the fastest variability time scale is under a millisecond, burst durations are usually very much longer. Long bursts have durations ranging from 10–1000 seconds, and even short bursts have a median duration of about 0.3 seconds. Clearly, whatever is the physical mechanism behind GRB production, it acts on a much longer time scale than the fastest dynamical time of the central engine.

Narayan, Paczyński & Piran (1992) suggested that the central engine in GRBs involves the accretion of matter onto a compact star, and that the energy in the burst is provided by the gravitational energy released by the accreting gas. In such a model, the duration of the burst is set by the viscous time scale of the accreting gas. In most accretion flows, the viscous time is significantly longer than the dynamical time, and so the accretion model naturally explains the large difference between the durations of bursts and their fastest variability time.

The formation of an accretion disk is a natural outcome of most popular models of GRBs, e.g. the mergers of double neutron star binaries (Eichler et al., 1989; Narayan, Paczyński & Piran 1992), neutron star black hole binaries (Paczyński, 1991; Narayan Paczyński & Piran 1992), black hole white dwarf binaries (Fryer et al., 1999), black hole helium star binaries (Fryer & Woosley, 1998), and models based on “failed supernovae” or “collapsars” (Woosley, 1993; Paczynski, 1998; MacFadyen & Woosley, 1999). An important exception is Usov’s (1992) model in which the GRB energy is provided by the magnetic and rotational energy of a newly formed rapidly rotating neutron star.

In this paper, we consider a generic accretion model of a GRB in which a certain amount of mass, $`m_dM_{}`$, goes into orbit around a relativistic star of mass $`3m_3M_{}`$. We assume that the orbiting mass is initially inserted into a torus at a radius $`r_{\mathrm{out}}R_S`$, where $`R_S`$ is the Schwarzschild radius of the central star: $`R_S=2GM/c^2=8.85\times 10^5m_3\mathrm{cm}`$. Starting from the initial toroidal configuration, the mass spreads out by viscosity and becomes an accretion flow extending from $`r=1`$ (the horizon of the central black hole) to $`rr_{\mathrm{out}}`$. (There might also be an outflow at radii $`>r_{\mathrm{out}}`$ as we discuss below.) We work out in this paper the time scale, $`t_{\mathrm{acc}}`$, of the accretion flow, the average mass accretion rate, $`\dot{M}=\dot{m}M_{}\mathrm{s}^1`$, onto the central star, the amount of mass, $`m_{\mathrm{acc}}=t_{\mathrm{acc}}\dot{m}`$, accreted by the star, and the accretion efficiency, $`\xi =m_{\mathrm{acc}}/m_d`$.

These parameters are constrained by observations. In binary merger models, the durations of bursts should be comparable to the accretion time $`t_{\mathrm{acc}}`$. Durations are of order a second or less for the class of “short duration GRBs” and in the range 10–1000 seconds for the class of “long duration GRBs.” The energy in long duration bursts is estimated to be $`5\times 10^{50}`$ erg (Panaitescu & Kumar, 2001; Frail et al. 2001). With a reasonable efficiency of converting accretion energy to relativistic flow ($`\stackrel{<}{}0.01`$) this corresponds to $`m_{\mathrm{acc}}\stackrel{>}{}0.1M_{}`$ and to a peak accretion rate $`\dot{M}\stackrel{>}{}10^2M_{}\mathrm{s}^1`$. The values of these parameters for short bursts are less certain, as the distance scale to this population is uncertain. However, it is unlikely that they are smaller by more than an order of magnitude.

The mass accretion rate tends to be extremely high for a typical GRB model — it is of order a fraction of a solar mass per second. At such accretion rates, the optical depth of the accreting gas is enormous and radiation is trapped inside the gas. In the normal course, the accretion would proceed via a radiatively inefficient flow, such as an advection-dominated accretion flow or the related convection-dominated accretion flow. We consider such flows in §2. If the accreting gas has a sufficiently high temperature and density, however, it can cool via neutrino emission, leading to a neutrino-dominated accretion flow, as discussed by Popham, Woosley & Fryer (1999) (see also Ruffert & Janka, 1999, for a numerical simulation of a binary merger which included neutrino losses). We discuss this case in §3. In §4 we go beyond the analytical results of §§2,3 and present numerical results, delineating the regions of $`(r_{\mathrm{out}},m_d)`$ space where different types of accretion occur. We discuss in §5 the implications of the results for various models of GRBs.

## 2 Radiatively Inefficient Accretion: CDAF

At the mass accretion rates of interest to us, the optical depth of the gas is extremely high, and the radiation is very effectively trapped (we demonstrate this below). The accretion flow then corresponds to an advection-dominated accretion flow, or ADAF (Narayan & Yi 1994, 1995; Abramowicz et al. 1995; Chen et al. 1995; see Narayan, Mahadevan & Quataert 1998 and Kato, Fukue & Mineshige 1998 for reviews). Radiation-trapped ADAFs were initially discussed by Katz (1977) and Begelman (1978), and later analyzed via height-integrated “slim disk equations” by Abramowicz et al. (1988).

Igumenshchev and co-workers have carried out hydrodynamic simulations of ADAFs and have discovered several interesting properties of these flows (Igumenshchev, Chen & Abramowicz 1996; Igumenshchev & Abramowicz 1999, 2000). They find that, when the dimensionless viscosity parameter $`\alpha `$ is large, say $`\alpha \stackrel{>}{}0.3`$, the accretion flow has a strong bipolar outflow, as anticipated in some previous papers (Narayan & Yi 1994, 1995; Blandford & Begelman 1999). On the other hand, when the viscosity is relatively weak, say $`\alpha \stackrel{<}{}0.1`$, the flow has well-developed convection (which was again anticipated in earlier work, cf. Narayan & Yi 1994, 1995). A convective ADAF has quite unusual properties (Stone, Pringle & Begelman 1999; Narayan, Igumenshchev & Abramowicz 2000; Quataert & Gruzinov 2000; Igumenshchev, Abramowicz & Narayan 2000), which arise because convection moves angular momentum inward rather than outward (Narayan et al. 2000; Quataert & Gruzinov 2000). Convective ADAFs have been given the name of convection-dominated accretion flows, or CDAFs.

As in other accretion flows, “viscosity” in radiation-trapped flows is believed to arise from magnetic stresses resulting from the Balbus-Hawley instability (Balbus & Hawley,, 1991). Numerical simulations of shearing MHD flows give values of $`\alpha `$ in the range $`\alpha \stackrel{<}{}0.1`$. We therefore assume that radiation-trapped accretion flows in GRBs also have $`\alpha `$ in this range (no 3D MHD simulations of these flows have been done so far, but there has been some 2D work by Stone & Pringle 2001). Given the relatively low value of $`\alpha `$, we expect the accretion flow to take the form of a CDAF.

Ball, Narayan & Quataert (2001) have given approximate scalings for various fluid variables in a CDAF. The density and velocity scale as

$`\rho (r)`$ $``$ $`\rho _{\mathrm{out}}\left({\displaystyle \frac{r}{r_{\mathrm{out}}}}\right)^{1/2}=2.97\times 10^{14}m_3^3m_dr_{\mathrm{out}}^{5/2}r^{1/2}\mathrm{g}\mathrm{cm}^3,`$
$`v(r)`$ $``$ $`cr^{3/2}=3\times 10^{10}r^{3/2},`$ (1)
where, in the first equation, we have used $`4\pi r_{\mathrm{out}}^3R_S^3(H/R)\rho _{\mathrm{out}}=m_dM_{}`$ to relate $`\rho _{\mathrm{out}}`$ in the CDAF to the mass in the accretion flow:

$$\rho _{\mathrm{out}}=2.97\times 10^{14}m_3^3m_dr_{\mathrm{out}}^3\mathrm{g}\mathrm{cm}^3.$$
(2)
Note the use of $`4\pi `$ rather than $`4\pi /3`$ in the equation for the mass. This has been done to obtain a better match with the results for an NDAF (see §3).

We assume that the accreting gas consists of photo-disintegrated nuclei with roughly equal numbers of neutrons and protons (and electrons). The optical depth through the flow is

$$\tau _{\mathrm{out}}\frac{\rho _{\mathrm{out}}}{2m_p}\sigma _Tr_{\mathrm{out}}R_S=5.21\times 10^{19}m_3^2m_dr_{\mathrm{out}}^2.$$
(3)
The diffusion time for radiation to leak out of the flow is then

$$t_{\mathrm{diff}}\tau _{\mathrm{out}}\frac{r_{\mathrm{out}}R_S}{c}=1.54\times 10^{15}m_3^1m_dr_{\mathrm{out}}^1\mathrm{s}.$$
(4)
This time is very much longer than the accretion time for the parameter ranges of interest to us, namely $`m_31`$, $`m_d0.11`$, $`r_{\mathrm{out}}1010^4`$. Thus we expect radiation to be very effectively trapped within the accretion flow.

Ball et al. (2001) also estimate the isothermal sound speed $`c_s`$ and the scale height $`H`$ of the accretion flow. We use slightly different coefficients here which are more appropriate for a radiation-pressure dominated $`\gamma =4/3`$ gas (as opposed to the $`\gamma =5/3`$ gas that Ball et al. 2001, considered):

$`c_s^2`$ $``$ $`0.3c^2r^1=2.70\times 10^{20}r^1\mathrm{cm}^2\mathrm{s}^2,`$
$`{\displaystyle \frac{H}{R}}`$ $`=`$ $`{\displaystyle \frac{c_s}{\mathrm{\Omega }_KR}}0.77.`$ (5)
These scalings are approximate, but they are likely to be accurate enough for the purposes of this paper. The mass accretion rate in the CDAF is estimated to be $`\dot{M}=4\pi r^2R_S^2(H/R)\rho (r)v(r)`$, which gives, using the above relation for $`H/R`$,

$$\dot{m}=3.39\times 10^4m_3^1m_dr_{\mathrm{out}}^{5/2}.$$
(6)
(Recall that $`\dot{m}`$ is defined as the mass accretion rate in units of solar masses per second.)

The accretion time scale is not simply equal to $`m_d/\dot{m}`$. The reason is that much of the mass in a CDAF actually flows out of the system rather than into the central black hole (Stone et al. 1999; Narayan et al. 2000; Quataert & Gruzinov 2000Igumenshchev & Abramowicz 2000). We therefore proceed as follows.

The random velocities of convective blobs in a CDAF are typically less than the local Keplerian velocity by a factor which is proportional to the viscosity parameter $`\alpha `$ (see Narayan et al. 2000). We thus write approximately

$$v_{\mathrm{turb}}\alpha v_K=2.12\times 10^9\alpha _1r^{1/2}\mathrm{cm}\mathrm{s}^1.$$
(7)
The residence time of a convective blob at radius $`r`$ is then

$$t_{\mathrm{res}}\frac{rR_S}{v_{\mathrm{turb}}}=4.17\times 10^4\alpha _1^1m_3r^{3/2}\mathrm{s}.$$
(8)
We make the reasonable assumption that the accretion time is of the same order as $`t_{\mathrm{res}}`$ at $`r=r_{\mathrm{out}}`$:

$$t_{\mathrm{acc}}t_{\mathrm{res}}(r_{\mathrm{out}})=4.17\times 10^4\alpha _1^1m_3r_{\mathrm{out}}^{3/2}\mathrm{s}.$$
(9)
We then find that the amount of mass accreted by the black hole is $`m_{\mathrm{acc}}=t_{\mathrm{acc}}\dot{m}`$, but since this mass cannot exceed the total available mass $`m_d`$ we write:

$`m_{\mathrm{acc}}`$ $`=`$ $`m_d,r_{\mathrm{out}}14.1\alpha _1^1m_d,`$ (10)
$`=`$ $`14.1\alpha _1^1m_dr_{\mathrm{out}}^1,r_{\mathrm{out}}>14.1\alpha _1^1m_d.`$
Note that, when $`r_{\mathrm{out}}`$ is large, the accreted mass is much less than $`m_d`$. The reason is that the bulk of the mass is ejected from the system, flowing out at $`rr_{\mathrm{out}}`$. The energy for the ejection is provided by convective energy flux from the interior of the flow.

The coefficient 14 in equation (10) is somewhat uncertain since we do not know the exact relation between $`t_{\mathrm{acc}}`$ and $`t_{\mathrm{res}}`$; there could well be a numerical factor other than unity relating the two. A different (and more more detailed) derivation of equation (10) is given in the Appendix, where we again find that the coefficient is uncertain. In the rest of the paper we use equation (10) as written, but we should keep in mind that it could be in error by a factor of a few.

Let us now calculate the temperature of the CDAF. The pressure in the accreting gas is given by

$$p=\rho c_s^2=8.02\times 10^{34}m_3^3m_dr_{\mathrm{out}}^{5/2}r^{3/2}\mathrm{erg}\mathrm{cm}^3.$$
(11)
The pressure has three contributions: radiation pressure, gas pressure, and degeneracy pressure (Popham et al. 1999):

$$p=\frac{11}{12}aT^4+\frac{\rho kT}{m_p}+\frac{2\pi hc}{3}\left(\frac{3}{8\pi m_p}\right)^{4/3}\left(\frac{\rho }{\mu _e}\right)^{4/3}.$$
(12)
The quantity $`a`$ is the radiation constant, and the factor 11/12 includes the contribution of relativistic electron-positron pairs (asssuming that the temperature is sufficiently above the pair threshold limit). The gas pressure term includes the contributions from non-relativistic particles. We assume that we have an equal mix of protons and neutrons (i.e. we assume that all complex nuclei have been photo-disintegrated). In the degeneracy pressure term we use an electron molecular weight $`\mu _e=2`$, assuming an equal mix of protons and neutrons; the ratio of neutrons to protons in the interior of a NS is about 0.2. The contribution of $`e^{}`$$`e^+`$ pairs has been ignored in the degeneracy pressure term. For the calculations presented in §4, we solve for $`T`$ by using the full expression for $`p`$ given in equation (12). But here, in order to obtain simple analytic estimates, we simplify the equation and assume that radiation pressure dominates. We then find that

$$T=1.84\times 10^{12}m_3^{3/4}m_d^{1/4}r_{\mathrm{out}}^{5/8}r^{3/8}K.$$
(13)
This estimate of $`T`$ is valid only if gas pressure is smaller than radiation pressure. The gas pressure is equal to

$$p_{\mathrm{gas}}=\frac{\rho kT}{m_p}=4.52\times 10^{34}m_3^{15/4}m_d^{5/4}r_{\mathrm{out}}^{25/8}r^{7/8}\mathrm{erg}\mathrm{cm}^3.$$
(14)
For this to be smaller than radiation pressure we require

$$\frac{r}{r_{\mathrm{out}}}<1.5m_3^{6/5}m_d^{2/5}.$$
(15)
For our fiducial black hole of mass $`3M_{}`$, gas pressure dominates only if the disk is quite massive, $`m_d\stackrel{>}{}5`$, and if $`r`$ is close to $`r_{\mathrm{out}}`$. (The detailed numerical calculations in §4 show that gas pressure dominates for $`m_d\stackrel{>}{}1`$ rather than 5, but even for such values of $`m_d`$, if $`r`$ is much less than $`r_{\mathrm{out}}`$, radiation pressure takes over.)

Let us next calculate the rate of cooling of the accreting gas. We showed earlier that cooling via radiative diffusion is negligible. However, at the temperatures found in these flows, cooling via neutrino losses may be important (cf. Popham et al. 1999). The cooling rate per unit volume due to neutrinos, $`q^{}`$, takes the form

$$q^{}=q_{\nu \overline{\nu }}^{}+q_{eN}^{}5\times 10^{33}T_{11}^9+9.0\times 10^{23}\rho T_{11}^6\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1,$$
(16)
where the first term on the right describes cooling via pair annihilation (the so-called URCA process) and the second term describes cooling via pair capture on nuclei (estimated for $`X_{\mathrm{nuc}}=1`$ as appropriate for our fully photo-disintegrated nuclear gas, cf. Popham et al. 1999). For the range of parameters of interest to us, $`q_{eN}^{}`$ invariably dominates over $`q_{\nu \overline{\nu }}^{}`$. We therefore neglect $`q_{\nu \overline{\nu }}^{}`$ hereafter.

We then obtain the following estimate for the cooling rate in the CDAF:

$$q_{eN}^{}=1.04\times 10^{46}m_3^{15/2}m_d^{5/2}r_{\mathrm{out}}^{25/4}r^{11/4}\mathrm{erg}\mathrm{cm}^3\mathrm{s}^1.$$
(17)
Since the energy density of the radiation-dominated gas is equal to $`3p_{\mathrm{rad}}=(11/4)aT^4`$, we estimate the cooling time of the gas at radius $`r`$ to be

$$t_{\mathrm{cool}}=\frac{3p_{\mathrm{rad}}}{q_{eN}^{}}=2.31\times 10^{11}m_3^{9/2}m_d^{3/2}r_{\mathrm{out}}^{15/4}r^{5/4}\mathrm{s}.$$
(18)
One of the important requirements for the accretion flow to behave like a CDAF is that it should be radiatively inefficient. We, therefore, require the cooling time $`t_{\mathrm{cool}}`$ at each radius to be longer than the residence time $`t_{\mathrm{res}}`$ of a convecting blob at that radius. The ratio of cooling to residence time is

$$\frac{t_{\mathrm{cool}}}{t_{\mathrm{res}}}=5.54\times 10^8\alpha _1m_3^{7/2}m_d^{3/2}r_{\mathrm{out}}^{15/4}r^{1/4}.$$
(19)
Since the ratio decreases with increasing $`r`$, we set $`r=r_{\mathrm{out}}`$. Then we see that a CDAF is possible only if the following condition is satisfied:

$$r_{\mathrm{out}}>118\alpha _1^{2/7}m_3^1m_d^{3/7}.$$
(20)
If the mass in the accretion flow is initially at a radius greater than the above limit, then the flow will become a CDAF at $`r=r_{\mathrm{out}}`$, and indeed for all $`r<r_{\mathrm{out}}`$. The scaling relations derived in this section would then be valid. If, however, the initial radius of the gas is below the above limit, then neutrino-cooling will prevail, and we will have a cooling-dominated accretion flow at all radii from $`r_{\mathrm{out}}`$ down to the black hole.

For completeness, we consider also the case when gas pressure dominates (which requires a large value of $`m_d`$, and $`r`$ close to $`r_{\mathrm{out}}`$, as shown above). The temperature is then obtained by setting $`p=p_{\mathrm{gas}}`$. This gives

$$T=\frac{m_pc_s^2}{k}=3.27\times 10^{12}r^1K.$$
(21)
The analysis presented so far assumes that we are given the initial mass in the disk $`m_d`$ and the initial radius of the gas $`r_{\mathrm{out}}`$. Another possibility is that we have steady injection of mass at a rate $`\dot{m}_{\mathrm{inj}}M_{}\mathrm{s}^1`$ at a circularization radius $`r_{\mathrm{out}}`$. This is the case, for instance, in the collapsar model, where the material is supplied by fallback from a supernova explosion. We may assume that the accretion flow achieves an approximate steady state in which the rate of injection of mass equals the rate of mass loss from the flow (both inflow and outflow): $`\dot{m}_{\mathrm{inj}}=m_d/t_{\mathrm{acc}}`$. Therefore, we obtain

$$m_d=4.17\times 10^4\alpha _1^1m_3\dot{m}_{\mathrm{inj}}r_{\mathrm{out}}^{3/2}.$$
(22)
This expression, which gives the mapping between $`m_d`$ and $`\dot{m}_{\mathrm{inj}}`$, may be substituted into the various relations derived in this section to obtain the corresponding results for the case of steady mass injection.

## 3 Radiatively Efficient Accretion: NDAF

We saw in the previous section that when the mass is initially injected at an outer radius smaller than the limit given in equation (20), cooling via neutrino emission becomes significant and the flow is no longer radiatively inefficient. We then have a regime of accretion in which the viscous energy dissipation is balanced by neutrino cooling. Popham et al. (1999) named this a neutrino-dominated accretion flow (NDAF) and worked out its properties. We extend their results in this section.

Because the gas cools efficiently, an NDAF behaves like a thin accretion disk and we are entitled to use the basic theory of thin disks (Shakura & Sunyaev 1973; Frank, King & Raine 1992). As before, the pressure of the gas has three contributions and is described by equation (12). The isothermal sound speed $`c_s`$ and the vertical scale height $`H`$ are given by $`c_s^2=p/\rho `$ and $`H=c_s/\mathrm{\Omega }_K`$, where $`\mathrm{\Omega }_K=(GM/R^3)^{1/2}=2.40\times 10^4m_3^1r^{3/2}\mathrm{s}^1`$ is the Keplerian angular velocity. We write the coefficient of kinematic viscosity in the usual form as

$$\nu =\alpha \frac{c_s^2}{\mathrm{\Omega }_K},$$
(23)
where $`\alpha `$ is a dimensionless parameter which we expect to have a value $`0.1`$. We denote, therefore, $`\alpha _1=\alpha /0.1`$, and write down our scalings in terms of $`\alpha _1`$.

A thin accretion disk satisfies two equations. Angular momentum balance gives

$$\dot{M}=3\pi \nu \mathrm{\Sigma }\left[1\left(\frac{R}{R_{}}\right)^{1/2}\right]6\pi \nu \rho H,$$
(24)
where $`\mathrm{\Sigma }=2\rho H`$ is the surface density and $`R_{}`$ is the radius of the inner edge of the disk ($`3R_S`$ for a Schwarzschild black hole). The approximation on the right is valid for $`RR_{}`$. The radial velocity of the gas is given by $`v=3\nu /2R`$ and so the accretion time is

$$t_{\mathrm{acc}}=\frac{R_{\mathrm{out}}}{v(R_{\mathrm{out}})}=\frac{2R_{\mathrm{out}}^2}{3\nu }.$$
(25)
If we write the disk mass as $`m_d=t_{\mathrm{acc}}\dot{M}`$, then by combining the previous two equations we find that

$$m_d=4\pi R_{\mathrm{out}}^3\left(\frac{H}{R}\right)_{\mathrm{out}}\rho _{\mathrm{out}}.$$
(26)
Note that we used this relation, with the same coefficient $`4\pi `$, also for a CDAF.

The condition of energy balance (viscous heating equals radiative losses) gives:

$$\frac{3GM\dot{M}}{8\pi R^3}=q^{}H=(q_{eN}^{}+q_{\nu \overline{\nu }}^{})H.$$
(27)
This relation closes our set of equations and allows us to solve for $`\rho `$, $`T`$ and other quantities. We show numerical solutions of the full equations in §4.

In the rest of this section we derive analytical scalings by making some approximations. First, the numerical calculations show that the $`q_{eN}^{}`$ term in the cooling law dominates over $`q_{\nu \overline{\nu }}^{}`$ for all parameters of interest. We will therefore assume this. The calculations also show that for most parameters, either gas pressure or degeneracy pressure dominates. (There is a small region of parameter space where radiation pressure dominates; we ignore this region in the analytical work.)

Let us first assume that gas pressure dominates. Expressing all results in terms of the scaled temperature, $`T_{11}=T/10^{11}K`$, and substituting in the angular momentum equation, we obtain a relation between $`\rho `$ and $`T_{11}`$

$$\rho =2.56\times 10^{13}\alpha _1^1m_3^2\dot{m}r^3T_{11}^{3/2}\mathrm{g}\mathrm{cm}^3.$$
(28)
Substituting this in the energy equation, we can solve for $`T_{11}`$, and thereby obtain the various other quantities (Popham et al. have derived similar relations):

$`T_{11}`$ $`=`$ $`0.548\alpha _1^{1/5}m_3^{1/5}r^{3/10},`$
$`t_{\mathrm{acc}}`$ $`=`$ $`2.76\times 10^2\alpha _1^{6/5}m_3^{6/5}r_{\mathrm{out}}^{4/5}\mathrm{s},`$
$`\dot{m}`$ $`=`$ $`36.2\alpha _1^{6/5}m_3^{6/5}m_dr_{\mathrm{out}}^{4/5},`$ (29)
$`\rho _{\mathrm{out}}`$ $`=`$ $`2.28\times 10^{15}\alpha _1^{1/10}m_3^{29/10}m_dr_{\mathrm{out}}^{67/20}\mathrm{g}\mathrm{cm}^3.`$
In a cooling-dominated thin disk, very little mass is expected to be lost to outflows, so we expect nearly all the mass in the disk to be accreted by the star, i.e.

$$m_{\mathrm{acc}}m_d.$$
(30)
The above results are valid provided gas pressure dominates. By comparing gas pressure to degeneracy pressure we can determine the condition for this to be true. We find that we require

$$r_{\mathrm{out}}>26.2\alpha _1^{2/7}m_3^{46/49}m_d^{20/49}.$$
(31)
When the outer radius is below this limit, degeneracy pressure takes over from gas pressure and we obtain a different set of analytical scalings:

$`\dot{m}`$ $`=`$ $`355\alpha _1m_3^{13/7}m_d^{9/7}r_{\mathrm{out}}^{3/2},`$
$`\rho _{\mathrm{out}}`$ $`=`$ $`7.32\times 10^{14}m_3^{18/7}\dot{m}^{6/7}r_{\mathrm{out}}^3\mathrm{g}\mathrm{cm}^3.`$
$`T_{11}`$ $`=`$ $`0.800\alpha _1^{1/6}m_3^{13/42}m_d^{1/21}r_{\mathrm{out}}^{5/12},`$ (32)
$`t_{\mathrm{acc}}`$ $`=`$ $`2.82\times 10^3\alpha _1^1m_3^{13/7}m_d^{2/7}r_{\mathrm{out}}^{3/2}\mathrm{s}.`$
In our analysis of both the gas-pressure-dominated and degeneracy-pressure-dominated regimes, we assumed that the accreting gas is optically thin to its own neutrino emission. This assumption breaks down at sufficiently small radii. We may estimate the neutrino optical depth as

$$\tau _\nu =\frac{q_{eN}^{}H}{4\sigma T^4}=5.20\times 10^3\alpha _1^{2/3}m_3^{55/21}m_d^{23/21}r_{\mathrm{out}}^{17/6}.$$
(33)
The radius at which the optical depth goes to unity is

$$r_{\tau =1}=20.5\alpha _1^{4/17}m_3^{110/119}m_d^{46/119}.$$
(34)
This radius lies within the degeneracy-pressure-dominated zone. Inside this radius, we need to consider neutrino transport in more detail. This is beyond the scope of the paper.

It is interesting to examine whether the NDAF solution is stable. Following Piran (1978) we use the general condition for thermal stability:

$$\left(\frac{d\mathrm{ln}Q^+}{d\mathrm{ln}H}\right)_{|\mathrm{\Sigma }}<\left(\frac{d\mathrm{ln}Q^{}}{d\mathrm{ln}H}\right)_{|\mathrm{\Sigma }},$$
(35)
where $`Q^\pm `$ are the integrated (over the height of the disk) heating ($`+`$) and cooling ($``$) rates. For an NDAF with viscosity described by the $`\alpha `$-prescription, the heating rate goes as $`Q^+H^2`$, and the vertically integrated cooling rate (for pair capture on nuclei) goes as $`Q^{}\mathrm{\Sigma }T^6`$. For an NDAF in which gas pressure dominates, $`TH^2`$ and $`Q^{}\mathrm{\Sigma }H^{12}`$. The criterion (35) is satisfied easily and we conclude that such an NDAF is thermally stable.

If radiation pressure dominates, $`T^4p\mathrm{\Sigma }H`$ and $`Q^{}\mathrm{\Sigma }^{5/2}H^{3/2}`$. The condition (35) is not satisfied and the NDAF is unstable. This resembles the situation in “conventional” accretion disks whose inner regions become thermally unstable when radiation pressure dominates (Frank et al. 1992).

If degeneracy pressure dominates, the temperature is independent of either $`H`$ or $`\mathrm{\Sigma }`$. The temperature of the disk can then freely adjust such that $`Q^{}`$ balances any $`Q^+`$. The situation is clearly thermally stable.

The condition for viscous stability is:

$$\frac{d\dot{M}}{d\mathrm{\Sigma }}>0.$$
(36)
We have $`\dot{M}\mathrm{\Sigma }`$ for the gas pressure case, $`\dot{M}\mathrm{\Sigma }^7`$ for the radiation pressure case, and $`\dot{M}\mathrm{\Sigma }^{9/7}`$ for the degeneracy pressure case. All three cases are viscously stable.

Thus, combining all these results on stability, we find that an optically thin NADF is unstable only if it is radiation pressure dominated. There is a very narrow region of parameter space near the boundary between NDAFs and CDAFs in $`(r_{\mathrm{out}},\dot{m})`$ space where radiation pressure does dominate. This unstable region could conceivably play a role in determining the temporal behavior of some bursts. The stability properties of the optically thick regions of the NDAF remain to be worked out.

The disks that we are considering are massive and we should also consider gravitational instabilities. The CDAF zone is always gravitationally stable in our models; the Toomre $`Q`$ parameter is invariably much greater than unity. For gas pressure dominated NDAFs, we find that the Toomre $`Q`$ parameter is given by

$$Q=8.5\alpha _1^{1/10}m_3^{9/10}m_d^1r^{9/20}r_{out}^{4/5}.$$
(37)
We see that $`Q`$ decreases with increasing $`r`$, so that the flow is most unstable on the outside. However, even for $`r=r_{out}`$, $`Q`$ is sufficiently larger than unity (for all reasonable values of parameters) that we are guaranteed stability. A similar result applies to radiation pressure dominated NDAFs.

The case of degeneracy pressure dominated NDAFs is more complicated. For most of the parameter space $`(m_d,r_{out})`$ we find $`Q`$ to be greater than unity, but it is only marginally so ($`Q3`$). When $`r_{out}<10`$, we find that $`Q`$ might become less than unity, signifying gravitational instability. But this is also the region in which the disk becomes optically thick to neutrinos and the analysis we have carried out breaks down.

## 4 Numerical Results

We discuss in this section numerical results which we have obtained by solving the full equations. We assumed that we are given the initial mass of the accretion disk $`m_d`$ and the initial radius $`r_{\mathrm{out}}`$. For each choice of $`r_{\mathrm{out}}`$ and $`m_d`$, we first assumed that the flow consists of a CDAF and used the relations described in §2 to calculate the properties of the accreting gas. In particular, we used equations (11) and (12), with $`\rho `$ given by equation (1), to solve for the temperature $`T`$ as a function of $`r`$. We checked the flow at all radii from $`r=r_{\mathrm{out}}`$ down to $`r=1`$ to make sure that the gas is radiatively inefficient at all radii. Specifically, we estimated the cooling time $`t_{\mathrm{cool}}`$, using the cooling formula given in equation (16), and checked that $`t_{\mathrm{cool}}`$ is longer than the residence time of convective blobs $`t_{\mathrm{res}}`$ at that radius. If $`t_{\mathrm{cool}}>t_{\mathrm{res}}`$ at all $`r`$, then we identified the flow as a pure CDAF and estimated the accretion time $`t_{\mathrm{acc}}`$ and the total mass accreted $`m_{\mathrm{acc}}`$, using the formulae given in §2.

For a small region of parameter space near the transition between a pure CDAF and a pure NDAF, we found that the flow starts off as a CDAF at $`r=r_{\mathrm{out}}`$ and switches to an NDAF at a smaller radius. We estimated $`t_{\mathrm{acc}}`$ and $`m_{\mathrm{acc}}`$ for these cases by using a reasonable matching formula at the transition radius. (We do not provide the details here, since this case is seen for only a small range of parameters).

For $`r_{\mathrm{out}}`$ less than a critical value (whose value depends on $`m_d`$), we found that the flow cools too rapidly at $`r=r_{\mathrm{out}}`$ to be a CDAF. In such cases, the flow becomes an NDAF at $`r=r_{\mathrm{out}}`$ and remains an NDAF all the way down to $`r=1`$. We solved the corresponding set of equations (see §3), and estimated $`t_{\mathrm{acc}}`$ and $`m_{\mathrm{acc}}`$ appropriately.

Figures 1 & 2 show contours of $`t_{\mathrm{acc}}`$ and the accretion “efficiency” $`\xi m_{\mathrm{acc}}/m_d`$, plotted in the space of the two principal parameters of the problem, $`r_{\mathrm{out}}`$ and $`m_d`$. These models have the “canonical values” $`m_3=\alpha _1=1`$. The boundary between the CDAF and NDAF zones is clearly seen in both figures, but especially in Fig. 1. The numerically determined location of the CDAF-NDAF boundary is fairly close to the analytical approximation given in equation (20) and shown as the right-most dotted line in Fig. 1.

As equation (9) shows, the accretion time scale $`t_{\mathrm{acc}}`$ for a CDAF is independent of $`m_d`$ and depends only on $`r_{\mathrm{out}}`$. In contrast, $`t_{\mathrm{acc}}`$ has a more complicated dependence on $`r_{\mathrm{out}}`$ and $`m_d`$ for an NDAF (eqs. 29 and 32). These behaviors are seen clearly in Fig. 1.

For a given initial mass $`m_d`$, the accreted mass $`m_{\mathrm{acc}}`$ is equal to $`m_d`$, independent of $`r_{\mathrm{out}}`$, in the case of an NDAF, and so $`\xi =1`$ (see Fig. 2). However, if the flow becomes a CDAF (which happens for $`r_{\mathrm{out}}`$ greater than a critical value), $`m_{\mathrm{acc}}`$ is significantly less than $`m_d`$, and the rest of the mass is ejected from the system; $`m_{\mathrm{acc}}`$ scales as $`r_{\mathrm{out}}^1`$ in this case (cf eq 10). Therefore, $`\xi `$ becomes significantly less than unity. Figure 2 illustrates this dependence.

In some GRB models, such as those involving a BH-NS merger or the collapse of a very massive star, the mass of the black hole could be larger than $`3M_{}`$. We have therefore computed results for an $`M=30M_{}`$ black hole, i.e. $`m_3=10`$. Equations (20), (31) and (34) show that the various critical radii vary roughly inversely as $`m_3`$ (which is equivalent to saying that the physical radii $`R=rR_S`$ at which the corresponding transitions occur are independent of the black hole mass). We have confirmed this result in the detailed numerical calculations. As a result, for $`m_3=10`$, the NDAF zone shrinks by a factor of 10 in $`r_{\mathrm{out}}`$, and the various variants of the NDAF (degeneracy pressure dominated zone, optically thick zone) practically disappear. Thus, large black hole masses are not conducive to the formation of an interesting NDAF zone.

In addition, we investigated models with a smaller value of the viscosity coefficient: $`\alpha =0.01`$, i.e. $`\alpha _1=0.1`$. In this case, the various critical radii become larger by roughly a factor of 2, as expected from equations (20), (31) and (34). More importantly, $`t_{\mathrm{acc}}`$ increases by a factor of $`10`$ (see equation 25).

We also considered the case when the accretion is fed at a constant rate $`\dot{m}_{\mathrm{inj}}`$ (rather than being initiated with an instantaneous addition of mass $`m_d`$ as we have assumed so far). This case could be relevant for collapsar models. Figure 3 shows contours of accretion efficiency defined in this case as $`\xi =\dot{m}/\dot{m}_{\mathrm{inj}}`$, plotted in the space of $`r_{\mathrm{out}}`$ and $`\dot{m}_{\mathrm{inj}}`$ for the “canonical values” $`m_3=\alpha _1=1`$. The overall behavior is consistent with the scaling relations derived in §§2,3.

## 5 Discussion

The starting point for this work is the fact that the accretion flow in a putative GRB central engine can have two very different forms: we expect the accretion to occur as a radiatively inefficient CDAF (convection-dominated accretion flow, cf Narayan et al. 2000, Quataert & Gruzinov 2000) if the mass $`m_d`$ is introduced at a somewhat large outer radius, namely $`r_{\mathrm{out}}`$ greater than the limit given in equation (20), while we expect the accretion to proceed via a radiatively efficient NDAF (neutrino-dominated accretion flow, cf. Popham et al. 1999) if $`r_{\mathrm{out}}`$ is smaller than this limit. For a narrow zone in parameter space close to the CDAF/NDAF transition, it is possible for the flow to be a CDAF on the outside ($`r\stackrel{<}{}r_{\mathrm{out}}`$) and to switch to an NDAF on the inside. But this is rare. By and large, for most choices of $`r_{\mathrm{out}}`$ and $`m_d`$, the flow is either a CDAF at all radii or an NDAF at all radii. Since the two kinds of flow are very different from each other, there is a rather large difference in what an observer would see in the two cases.

In most of the region of $`(r_{\mathrm{out}},m_d)`$ space where the flow is a CDAF, the mass accretion rate $`\dot{m}`$ and the amount of mass accreted $`m_{\mathrm{acc}}`$ are both very small. This is because in a CDAF, especially when $`r_{\mathrm{out}}`$ is large, much more mass flows out of the system than into the black hole. If, as seems reasonable, a GRB engine requires a relatively large $`m_{\mathrm{acc}}`$ in order to produce a viable burst, then our work suggests that systems which form CDAFs are less likely to produce observable bursts.

In contrast, an NDAF has substantially larger values of $`\dot{m}`$ and $`m_{\mathrm{acc}}`$ for a given $`m_d`$. (In fact, we assume that $`m_{\mathrm{acc}}=m_d`$ for an NDAF; there is probably some mass loss in a wind, but we expect it to be a small fraction of the total mass.) Therefore, a source with an NDAF is a much more plausible model of a GRB engine. By this argument, neutrino cooling (the key ingredient of an NDAF) is important for an efficient GRB. Note that we are not assuming anything about the actual mechanism of a burst. The fireball may be produced through neutrino-antineutrino annihilation (Eichler et al.,, 1989), or it could be the result of some other mechanism. Our suggestion is that whatever the mechanism may be, if it is based on accretion, then it probably requires a large $`m_{\mathrm{acc}}`$ to operate efficiently; such an $`m_{\mathrm{acc}}`$ can be achieved only through neutrino cooling of the accreting gas in an NDAF model.

We note a possible caveat to the above conclusion. The energy in long duration GRBs is estimated to be about $`5\times 10^{50}`$ erg (Panaitescu & Kumar, 2001; Frail et al. 2001). This energy could in principle be generated in a CDAF model with $`m_{\mathrm{acc}}\stackrel{<}{}0.1M_{}`$, provided the gravitational energy of the accreted matter is used very efficiently to launch a high Lorentz factor wind. We feel, however, that it is reasonable to hypothesize that the efficiency for converting the gravitational energy of $`m_{\mathrm{acc}}`$ to a GRB-producing relativistic wind is low, say not larger than 1%, and that bursts arise only from systems that form NDAFs and have large $`m_{\mathrm{acc}}`$.

Figures 1 and 2 show results corresponding to a simple model in which we assume that the accretion disk is formed instantaneously, e.g. by the disruption of a companion star. We specify the initial state of the disk by giving its mass $`m_d`$ and initial radius $`r_{\mathrm{out}}`$. For such a model, we see that the region of parameter space where NDAFs form corresponds to short accretion times $`t_{\mathrm{acc}}`$ of order a few tenths of a second. This suggests that a binary-disruption-based accretion model is capable of producing short GRBs. It is, however, very hard to see how such an accretion system could produce a long burst.

One way to make a long burst from an NDAF is to decrease the value of the viscosity parameter $`\alpha `$ (see eqs. 29 and 32). Could $`\alpha `$ be significantly lower, say 0.01? MHD simulations of thin accretion disks generally give $`\alpha `$ in the range 0.01 upwards, and it is widely agreed that the values obtained are lower limits since the simulations have limited spatial resolution. Empirical estimates of $`\alpha `$ in cataclysmic variables (obtained by comparing observations of dwarf nova outbursts with model predictions) give $`\alpha 0.1`$ in the high state (which is most relevant for our models). There is clear evidence that $`\alpha `$ is smaller, perhaps $`0.01`$, in the low state of CVs. This is probably the result of the cold gas becoming neutral and losing its coupling to the magnetic field (Gammie & Menou 1998), which is clearly not relevant for our ultra-hot plasma. We feel that $`\alpha 0.1`$ is a reasonable estimate for the viscosity parameter.

Sakimoto & Coroniti (1981) suggested that $`\alpha `$ may be lower in radiation-dominated gases because the magnetic pressure may be in equipartition with only the gas pressure rather than the total pressure. Their proposal requires that the gas and the radiation be able to slip past each other (cf. Blaes & Socrates 2001). At the extraordinarily large radiation optical depths found in our models (both CDAFs and NDAFs), the gas and the radiation are extremely tightly coupled. It is therefore very unlikely that the value of $`\alpha `$ would be modified.

The models we consider have large disk masses. The accreting gas may therefore develop gravitational instabilities (see the discussion near the end of §3) and lose angular momentum via gravitational waves (see Bonnell & Pringle 1995 and references therein). This might lead to an increase in the effective value of $`\alpha `$. Another potential uncertainty is that, for some choices of the parameters, electron degeneracy becomes important. It is not understood how the Balbus-Hawley instability (which is thought to produce the shear stresses behind $`\alpha `$) behaves in such a gas.

While the time scale of the accretion $`t_{\mathrm{acc}}`$ depends fairly sensitively on $`\alpha `$, the size of the NDAF zone is insensitive to $`\alpha `$ (cf eq. 20). Thus, Figs. 1 and 2 give fairly reliable limits on the radius inside which the mass needs to be introduced if we wish to have an NDAF. We may use this information to deduce a few interesting results on GRB models.

NS-NS and BH-NS merger models, with $`(r_{\mathrm{out}},m_d)=`$ (10,0.1) and (10,0.5) (see Popham et al. 1999), are well inside the NDAF zone and, according to our calculations, are capable of producing GRBs. However, this is only if the black hole is small (few $`M_{}`$). If the black hole is larger than $`10M_{}`$ its Schwarszchild radius becomes too large and there is not enough “room” for an NDAF solution around it. Moreover, the neutron star in this case is swallowed whole by the BH and it is not tidally disrupted to create an accretion disk. On the other hand, as already noted, unless the viscosity is much smaller than what we have assumed (which we consider unlikely), such disks cannot produce long bursts lasting hundreds or even tens of seconds. This suggests that NS-NS mergers and BH-NS mergers with smallish BH masses produce the class of short duration GRBs, but not the long duration GRBs.

Other merger models, specifically the BH-WD and the BH-He star merger models, would appear not to be viable GRB engines. As the secondaries in these systems are not compact, they would form accretion flows with large values of $`r_{\mathrm{out}}`$. For instance, Popham et al. (1999) estimate $`r_{\mathrm{out}}3000`$ for a BH-WD binary and $`r_{\mathrm{out}}5000`$ for a BH-He star binary. At these radii, the accretion flow will be a very exteneded CDAF and hardly any mass will be accreted. Although the time scales of these models are consistent with long bursts, the extremely small value of $`m_{\mathrm{acc}}`$ suggests that these models do not produce GRBs of any kind. It is interesting to speculate what kind of observable events these binaries might produce (as undoubtedly they do merge in nature).

All of the discussion so far is concerned with binary mergers, where we have imagined that a certain fixed amount of mass is instantaneously input into the accretion flow. The popular collapsar model (MacFadyen & Woosley 1999) corresponds to a different scenario in which mass is steadily fed over a period of time by fallback from the supernova explosion. MacFadyen and Woosley (1999) show that the time scale of the GRB is set by the physics of fallback rather than by accretion. Further, the time scales they obtain are consistent with observations of long GRBs.

While the time scale may be set by fallback, the efficiency of the burst still depends on the nature of the post-fallback accretion. Let us ssume that fallback supplies mass at a certain injection rate $`\dot{m}_{\mathrm{inj}}`$ at a characteristic radius $`r_{\mathrm{out}}`$; the latter depends on the specific angular momentum of the material. Figure 3 shows the numerical results. As expected, efficient accretion, where most of the fallback material reaches the black hole, is possible only if $`r_{\mathrm{out}}`$ is small and falls within the NDAF zone. If collapsars have a distribution of $`r_{\mathrm{out}}`$, then our calculations suggest that only those systems that have $`r_{\mathrm{out}}\stackrel{<}{}100\alpha _1^{2/7}m_3^1`$ will make bursts. Systems with larger specific angular momentum, and hence larger $`r_{\mathrm{out}}`$, will form CDAFs and will eject most of the mass. Such systems may make very interesting supernova explosions, but if at all they make GRBs the bursts are likely to be very weak.

We have assumed in this paper that the energy of a GRB is proportional to the total mass accreted on the central compact object. For NS-NS and NS-BH binaries, which give rise to short duration bursts, the disk mass is expected to be nearly constant and hence the total energy might be roughly the same for all bursts. However, for long bursts, which we assume result from the collapse of a massive star, the energy release depends on the angular momentum of the stellar core (which determines $`r_{\mathrm{out}}`$) and the mass of the stellar envelope that falls back on the collapsed core. Observations of long duration bursts indicate that the energy does not vary much from one burst to another (Panaitescu and Kumar 2001, Frail et al. 2001). This means that, for some reason, the total accreted mass on the central object depends only weakly on the properties of the progenitor star.

In our calculations we neglected photo-disintegration of nuclei, which Popham et al. (1999) and MacFadyen & Woosley (1999) show to be an important coolant of the accreting gas. If we include this effect, the boundary between CDAFs and NDAFs will move to somewhat higher values of $`r_{\mathrm{out}}`$, perhaps by a factor of a few. However, it will not change our key conclusions.

We are grateful to Eliot Quataert and Andrew MacFadyen for useful discussions, and Martin Rees for comments on the manuscript. We thank the organizers of the Jerusalem Winter School in Physics for hospitality while part of this research was carried out. This work was supported in part by NSF grant AST 9820686 and by a US-Israel BSF grant 98-00225.

## Appendix A Mass Accretion versus Ejection in CDAF Models

We provide an alternate “derivation” of equation (10), and discuss the relevant uncertainties in the result.

In a CDAF, there is a flow of energy from the inside to the outside as a result of convection. Let us write the convective luminosity as

$$L_c=ϵ_c\dot{M}_{\mathrm{in}}c^2,$$

where $`\dot{M}_{\mathrm{in}}`$ is the mass accretion rate onto the black hole. Since the convective luminosity is proportional to $`\alpha `$ (Quataert & Gruzinov 2000; Narayan et al. 2000), we may write $`ϵ_c=\alpha \eta `$, with $`\eta `$ roughly a constant. Ball et al. (2001) computed a global model of a CDAF with $`\alpha =0.03`$, and found that $`ϵ_c=0.0045`$. This suggests that $`\eta 0.15`$.

In our GRB model, there is practically no radiative emission from the gas. Therefore, the entire convective luminosity $`L_c`$ must be converted into thermal and mechanical energy of the gas on the outside. This energy will cause the gas at $`rr_{\mathrm{out}}`$ to move to larger radii.

How much energy do we need to evaporate gas from $`r=r_{\mathrm{out}}`$? Since the gas in a CDAF is hardly bound to the black hole — the Bernoulli parameter is close to zero or even positive (cf Quataert & Gruzinov 2000; Narayan et al. 2000) — very little energy is needed to drive a mass outflow to infinity. Let us write the energy carried away per unit mass of escaping gas as $`\zeta GM/r_{\mathrm{out}}R_S`$. The rate at which mass is ejected is then given by

$$\dot{M}_{\mathrm{out}}=\frac{L_cr_{\mathrm{out}}R_S}{\zeta GM}=\frac{0.2\eta \alpha _1}{\zeta }r_{\mathrm{out}}\dot{M}_{\mathrm{in}}.$$

Thus we have

$$\frac{m_{\mathrm{acc}}}{m_d}=\frac{\dot{M}_{\mathrm{in}}}{\dot{M}_{\mathrm{out}}}=\frac{5\zeta }{\eta }\alpha _1^1r_{\mathrm{out}}^1.$$

Equation (10) of the main paper (derived by a different approach) has the same scaling, but with a coefficient equal to 14.

By the Bernoulli parameter argument mentioned earlier, we expect $`\zeta `$ to be significantly less than unity; we would arbitrarily guess that $`\zeta 0.10.3`$. Thus we expect $`5\zeta /\eta 10`$ (compared to 14 in the main text). However, the coefficient is uncertain since $`\zeta `$ could be much smaller than our assumed value, or conceivably much larger. In the latter case, the convective luminosity must be carried away by a small amount of gas that is accelerated to a speed much greater than the escape velocity. This possibility is not supported by the numerical simulations carried out to date. Our estimate of $`\eta `$ is also uncertain since it is based on a single global model calculated by Ball et al. (2001). A better estimate could be obtained from full-scale numerical simulations.

The simulations reported by Igumenshchev & Abramowicz (2000) may be used to obtain a direct estimate of $`\dot{M}_{\mathrm{in}}/\dot{M}_{\mathrm{out}}`$. For their model L ($`\gamma =4/3`$, $`\alpha =0.03`$, $`r_{\mathrm{out}}7000`$), they estimate $`\dot{M}_{\mathrm{in}}/\dot{M}_{\mathrm{out}}0.003`$, whereas equation (10) of the present paper predicts a value of 0.007. The two estimates agree to within a factor of about 2, which is reassuring. On the other hand, for $`\alpha =0.1`$ (Model I), they obtain $`\dot{M}_{\mathrm{in}}/\dot{M}_{\mathrm{out}}0.02`$ which is much larger than our formula would predict. A possible explanation is that Model I does not strictly correspond to a CDAF. The flow appears to be a transition case where turbulent convection is replaced by a large scale circulation. In the case of MHD simulations (Stone & Pringle 2000) the mass outflow rate is again found to be much larger then the mass accretion rate on the central object. Therefore, the main results discused here should still be valid.


format: text

text:
# Giant anharmonicity and non-linear electron-phonon coupling in MgB2; Combined first-principles calculations and neutron scattering study

\[

## Abstract

We report first-principles calculations of the electronic band structure and lattice dynamics for the new superconductor MgB2. The excellent agreement between theory and our inelastic neutron scattering measurements of the phonon density of states gives confidence that the calculations provide a sound description of the physical properties of the system. The numerical results reveal that the in-plane boron phonons (with E2g symmetry) near the zone-center are very anharmonic, and are strongly coupled to the partially occupied planar B $`\sigma `$ bands near the Fermi level. This giant anharmonicity and non-linear electron-phonon coupling is key to quantitatively explaining the observed high Tc and boron isotope effect in MgB2

preprint: TO BE SUBMITTED TO SCIENCE!preprint: Submitted to Phys. Rev. Lett.

\]

The recent discovery of superconductivity at 40 K in the MgB2 binary alloy system has triggered enormous interest in the structural and electronic properties of this class of materials. The system has a very simple crystal structure, where the boron atoms form graphite-like sheets separated by hexagonal layers of Mg atoms (see inset to Fig. 1). Our pseudopotential plane wave band structure calculations show that the bands near the Fermi level arise mainly from the $`p_{x,y}`$ $`\sigma `$ bonding orbitals of boron, while the Mg does not contribute appreciably to the conductivity, in good agreement with initial reports from other groups. In the case of graphite these $`\sigma `$ bands are full, but for MgB2 they are partially unoccupied, creating a hole-type conduction band like the high-Tc cuprates. In contrast to the cuprates, however, the normal-state conductivity is three-dimensional in nature instead of being highly anisotropic, thus eliminating the “weak-link” problem that has plagued widespread commercialization of the cuprates. The normal-state conductivity is also one to two orders-of-magnitude higher than either the Nb-based alloys or Bi-based cuprates used in present day wires, and this feature combined with low cost and easy fabrication could make this class of materials quite attractive for applications.

From a fundamental point of view the central question is whether the high Tc in this new system can be understood within the framework of a conventional electron-phonon mechanism, or a more exotic mechanism is responsible for the superconducting pairing. The observed boron isotope effect argues for an electron-phonon mechanism, while the positive Hall coefficient suggests similarities with the cuprates. To answer this question, we have carried out inelastic neutron scattering measurements of the phonon density of states, and compare these results with detailed first-principles calculations of the lattice dynamics (and electronic) calculations for MgB2. Excellent agreement is found between theory and experiment. More importantly, the numerical results demonstrate that the in-plane boron phonons near the zone-center (with E2g symmetry at $`\mathrm{\Gamma }`$) are very anharmonic and strongly coupled to the partially occupied conduction band (planar B $`\sigma `$ bands) near the Fermi level, providing the large electron-phonon interaction in this system. Due to these strong anharmonic phonon modes, the electron-phonon coupling is non-linear, providing the essential ingredient to explain the high Tc and boron isotope effect in MgB2.

The results of the neutron measurements and calculations for the phonon density of states are summarized in Fig. 1. For the experimental data, the 5g polycrystalline sample was prepared in the usual way with the 11B isotope to avoid neutron absorption problems, and was characterized by magnetization, x-ray, and neutron diffraction. Inelastic neutron scattering measurements were made to determine the generalized phonon density of states (GPDOS), which is the phonon density of states weighted by the cross section divided by the mass of each atom. The data were collected from 7 K to 325 K on the Filter Analyzer Spectrometer in the range 5-130 meV, and on the Fermi chopper instrument for energies of 0.5-30 meV. Additional details of the experiment and analysis can be found elsewhere. The experimental data indicate two bands of phonons, one below 40 meV corresponding primarily to the acoustic phonon modes, and one above 50 meV that mostly involves the boron motions. The high-energy portion of the spectrum shows a sharp cutoff at about 100 meV. The data are in basic agreement with the results of Osborn et al., although the present data were obtained with better overall energy resolution. Our results do not agree with Sato et al., who reported a feature at 17 meV in the GPDOS. We did not observe this feature at any temperature and conclude that it is not intrinsic to the GPDOS. We also have measured the temperature dependence of the spectrum and find no substantial changes in the basic features from 7 K up to 200 K, while a very modest softening of some of the modes was observed in going to 325 K.

Fig. 1 (a) Generalized phonon density of states for MgB2 (top curve) determined from inelastic neutron scattering measurements, and the calculated phonon density of states (DOS) (bottom curve). The inset shows the unit cell of MgB2 along with the Brillouin zone and the high-symmetry directions. (b) Calculated dispersion curves of all phonons and the corresponding DOS (right panel). At the zone center ($`\mathrm{\Gamma }`$), the symmetries of the modes are also indicated. The modes indicated by the red-line dominate the electron-phonon coupling.

In order to understand the origin of the features observed in the neutron data, we have carried out detailed first-principles calculations of the electronic band structure and lattice dynamics. The calculations were performed using the pseudopotential plane wave method, within the generalized gradient approximation. We employed plane waves with an energy cutoff of 500 eV, and the ultra-soft pseudopotentials for Mg and B. The total energy and forces converged within 0.5 meV/atom and 0.01 eV/Å, respectively. Brillouin-zone integrations were carried out using $`dk=0.02`$ Å-1, generating $`15\times 15\times 11`$ and $`9\times 9\times 6`$ k-points for $`1\times 1\times 1`$ and $`2\times 2\times 2`$ supercells, respectively. The lattice parameters a and c were relaxed before the lattice dynamics calculations, and the force constant matrix was obtained by the direct-force method using periodically repeated supercells. The optimized lattice parameters are in good agreement with our experimental values within 1%. The calculated density of states obtained from these calculations is shown in Fig. 1a, and we see that the agreement between experiment and theory for the energies of the modes is excellent. The agreement for the intensities is also good considering that the observed spectrum is a weighted density of states rather than the actual DOS, and the wave vector averaging may include some coherency effects. The calculated phonon dispersion curves along the high-symmetry directions of the Brillouin-zone (see inset to Fig. 1) are shown in Fig. 1b. All the modes are found to be quite dispersive as expected for the small size of the unit cell, but divide nearly completely into acoustic (below $`40`$ meV) and optical (50-100 meV) bands, with the only exception being along the A-L direction. All the optical modes are weakly dispersive along the $`\mathrm{\Gamma }`$-A direction, reflecting the layered nature of the MgB2 crystal structure.

There are four distinct phonon modes at the zone center $`\mathrm{\Gamma }`$, as indicated in Fig. 1b. The A2u and B1g singly-degenerate modes involve only vibrations along the c axis; for B1g the boron atoms move in opposite directions while the Mg ion is stationary, while for the A2u mode the Mg and B planes move in opposite directions along the z-axis. The other two modes are doubly degenerate and involve only in-plane motions (along the x or y axes). For the E1u mode the Mg and B planes vibrate in opposite directions along the x (or y) axis. The calculated energies of all zone-center phonons are listed in Fig. 2, and for the three types of modes just discussed there is good agreement between calculations reported by other groups. For the E2g mode, however, there are large discrepancies between the results reported for various calculations. The boron ions for this mode vibrate in opposite directions along the x (or y) axis, with the Mg ions stationary. The discrepancies noted above are resolved when the strong anharmonicity associated with the in-plane displacements of the boron atoms is taken into account.

To further investigate this behavior, we plot in Fig. 2 the energy as the atoms are displaced by an amount $`u`$ according to each of the four eigenmodes. For the B1g, A1u, and E1u modes the distortion energy can be very well approximated by the harmonic expression $`E(u)=A_2u^2`$ up to $`u/a=0.065`$. It is unusual that the phonons remain harmonic at these large distortions, but this is consistent with the fact that we did not observe any substantial changes in the measured GPDOS with temperature.

Fig. 2 Energy curves as a function of boron displacement for each zone center frozen-phonon mode. The top panel shows that the B1g, A2u, and E1u modes are harmonic up to $`u/a0.065`$. The bottom panel shows that the E2g boron in-plane mode is very anharmonic. The inset indicates one of the E2g modes. Animations of these modes and their coupling with the electronic structure can be found at http://www.ncnr.nist.gov/staff/taner/mgb2.

The most interesting observation, though, is the totally opposite behavior of the E2g in-plane boron mode, as shown in Fig. 2b. The potential well is very shallow at small displacements and increases rapidly at large displacements, and can be fit very well to $`E(u)=A_2u^2+A_4u^4`$ with a large ratio of$`A_4/A_2^28`$. This indicates that this mode is unusually anharmonic and that we are in a non-perturbative regime. Hence the estimated harmonic energy, $`\omega _H`$(E2g) = 60.3 meV will be lower than the actual energy. Using Hartree-Fock decoupling, one gets $`E(u)=(A^2+3A_4<u^2>)u^2`$, where $`<u^2>=/2M\omega _{sch}`$ and $`\omega _{sch}`$ is the Self Consistent Harmonic (SCH) solution of $`\omega _{sch}^2=[(A_2+3A_4<u^2>)/M]^{1/2}`$. This yields a better estimate of $`\omega _{sch}`$ $`70`$ meV, a 17% enhancement of the harmonic phonon energy $`\omega _H(E_{2g})`$. We can calculate the exact energy levels of the potential by numerically solving the Schrodinger equation and obtain $`\omega (E_{2g})=74.5`$ meV, a 25% enhancement of the harmonic value and also the TC, which then matches with the peak in the GPDOS.

Since the in-plane motions of the boron will change the boron orbital overlap, significant electron-phonon coupling can be expected for the planer $`\sigma `$ boron conduction band at the Fermi level and this plays an important role in the superconducting pairing. This can be most easily seen by comparing the band structure of the undistorted lattice to that of distorted one by a zone center phonon. For the harmonic phonons we did not see any significant changes near the Fermi level. However, for E2g modes, there is a significant splitting of the bands as well as shifts near the Fermi level as shown in Fig. 3, indicating strong electron-phonon coupling.

Fig. 3 Band structure of the undistorted (left) and distorted structures (right) by E2g phonons ($`u_B0.06`$Å). Band structure of the other three frozen-phonon structures do not show any significant changes.

To determine the electron-phonon (EP) coupling quantitatively, we evaluated the Fermi-surface averaged deformation potential,

$$\mathrm{\Delta }=<[\delta ϵ(𝐤)\delta \mu ]^2>,$$
(1)
for each zone-center frozen phonon. Here $`\delta ϵ(𝐤)`$ is the change in the one-particle energy with momentum $`𝐤`$ due to the frozen phonon, $`\delta \mu `$ the corresponding change in the chemical potential. $`<>`$ denotes an average of $`𝐤`$ over the Fermi surface, which we have carried out numerically.

For the harmonic B1g, A2u, and E1u phonons we calculated an insignificant coupling, and conclude that the electron-phonon coupling is negligible for these modes. For the E2g modes, Fig. 4a shows that the coupling $`\mathrm{\Delta }(u)`$ is large and has both quadratic and quartic terms in frozen-phonon amplitude, $`u`$, indicating significant non-linear electron-phonon coupling. Neglecting non-linear cross terms, the electron-phonon coupling constant $`\lambda `$ for our anharmonic phonon is approximately given by

$$\lambda =N(E_F)\left(B_2^{}\underset{n}{}\frac{|<n|Q|0>|^2}{E_nE_0}+B_4^{}\underset{n}{}\frac{|<n|Q^2|0>|^2}{E_nE_0}\right),$$
(2)
where $`N(E_F)`$ is the total density of states at the Fermi energy and equal to 0.69 states/eV/cell. $`E_n`$ and $`|n>`$ are the eigenvalues and eigenfunctions of the oscillator E2g in its adiabatic potential shown in Fig. 2b. $`B_2^{}`$ and $`B_4^{}`$ are the first and second-order EP coupling, respectively, and obtained from expansion of $`\mathrm{\Delta }(u)`$. $`Q`$ is the normal coordinate which is related to boron displacement $`u`$ via a normalization constant. We calculate a total electron-phonon coupling $`\lambda =0.907(0.922)`$ for MB=10 (11). Using the McMillan expression for TC and taking a typical value for $`\mu ^{}=0.15`$, we obtain a TC of 39.4 K and 38.6 K for MB=10 and MB=11, respectively, yielding an isotope effect $`\alpha =0.21`$ (see Fig. 4b). These estimates are in excellent agreement with experiments. Within the above approximation, the non-linear EP coupling increases TC by about 10%. Our preliminary calculations of $`\lambda `$ including other non zone-center phonons indicates that the numbers will not change significantly. Therefore we conclude that MgB2 is a conventional electron-phonon superconductor in the strong-coupling regime.

Fig. 4 (a) Calculated Fermi-surface averaged deformation potential $`\mathrm{\Delta }(u)`$. The large quartic term indicates significant non-linear electron-phonon coupling. (b) Tc versus boron mass, indicating a boron isotope effect $`\alpha =0.21`$, significantly reduced from harmonic BCS value of 0.5 due to giant anharmonicity. Within our theory, we expect zero isotope effect for Mg, in good agreement with recent report of $`\alpha _{Mg}`$=0.02.

It is interesting to note that there seems to be a close correlation between the anharmonicity of the E2g in-plane modes and the B-B bond length (dBB). For MgB2, dBB is 1.764 Å, significantly stretched from its optimal value of 1.65 Å in elemental planar boron, probably due to repulsive interactions between the Mg and B ions. This explains the unusual anharmonicity and observed high TC in MgB2. For hypothetical CaB2, the calculated dBB is quite large (1.84 Å) and we calculate very large anharmonicity, which may explain the instability of CaB2. It then appears likely that MgB2 is fortuitously just at the phase boundary. It will be interesting to see if it is possible to expand the dBB in MgB2 by Ca substitution on the Mg site to increase TC further.

NOTE ADDED IN PROOF: After submission of our work, a recent Raman study has reported 77 meV for the E2g mode with a very large width, confirming our predictions.

ACKNOWLEDGEMENTS: We acknowledge many useful discussion with R. L. Cappelletti and D. A. Neumann.


text:
# Magneto-transport Properties near the Superconductor-Insulator Transition in 2D

From the earliest experiments on superconducting thin metal alloy films, a simple picture emerged for the magnetic-field tuned transition to the insulating state. Namely, the transition is continuous with a well-defined critical field, $`H_c`$, below which superconductivity occurs and above which an insulator obtains. At the critical field, the resistivity is independent of temperature. Within the dirty boson model in an applied magnetic field, these results follow naturally as a direct consequence of the inherent duality between Cooper pairs and vortices: above $`H_c`$, vortices condense while Cooper pairs remain localized but below $`H_c`$ it is the Cooper pair condensation and vortex localization that affords the zero-resistance state. At criticality, Cooper pairs on the brink of forming a phase coherent state diffuse with the quantum of conductance $`\sigma _Q=4e^2/h`$ for charge $`2e`$ bosons. However, more recent experiments have shown that this simple picture breaks down dramatically at and below $`H_c`$. In particular, a metallic phase emerges below $`H_c`$ with the zero-resistance state occurring only at a much smaller field $`H_MH_c`$. On phenomenological grounds, it has been suggested that the metallic phase arises from the quantum tunneling of field-induced vortices in the presence of dissipation.

As the metallic phase has been observed at zero-field in the disorder-tuned insulator-superconductor transition as well, we reinvestigated the conductivity of the standard quantum-disordered regime in zero field. We found that quite generally Cooper pairs lacking phase coherence form a Bose metal phase. This result arises from the fact that in the quantum-disordered regime the elementary excitations obey Boltzmann statistics. Hence, their concentration $`n`$ is exponentially suppressed. As a consequence, so is the inverse scattering time $`1/\tau `$ that is determined by the probability that two quasiparticles collide. But because the conductivity is a product of the density and the scattering time, the exponentials cancel one another, thereby leading to a finite conductivity as $`T0`$. In the current study, we focus on the fate of the Bose metal phase in a magnetic field. Near the quantum critical point where vortex dynamics are subdominant, the relevant physics can be accessed by the standard $`O(M)`$ Ginzburg-Landau (GL) theory. In this work, we utilize this approach for a system near the IST in the presence of a perpendicular magnetic field. We establish that at low temperature $`T`$ and close to the quantum critical point, $`\mathrm{\Delta }1`$, the application of a field $`H`$ drives the system to the high-field magnetically activated (MA) regime once $`\sqrt{H}\mathrm{\Delta },T`$, thereby compounding the analysis of the transport properties in the quantum disordered regime. In this regime, the energy levels of the bosonic quasiparticle excitations are quantized and have a gap $`m\sqrt{H}`$ that is inversely proportional to the correlation length $`\xi `$. Working within the large $`M`$ approximation, we discuss the crossovers and limiting forms of the scaling functions for the correlation length $`\xi =\xi (H,T,\mathrm{\Delta })`$ and the static conductivity $`\sigma `$ in the vicinity of the critical point. We present the qualitative form for the conductivity in the presence of dissipation based on the standard Kubo formula neglecting the effects of frustration. The dc conductivity is regularized by inclusion of the weak coupling to an Ohmic heat bath $`\eta `$ as well by quasiparticle scattering $`1/\tau `$ due to the quartic non-linearity in the GL action. Collectively, internal and external dissipation can be treated jointly with an effective coupling constant $`\stackrel{~}{\eta }=\eta +1/\tau `$. In the absence of the bath and at a small constant magnetic field such that $`\sqrt{H}\mathrm{\Delta }`$, the magnetic energy scale is irrelevant. Consequently, the Bose metal obtained previously in the quantum disordered regime is recovered. We show that in the MA regime, the conductivity $`\sigma `$ is directly, rather than inversely, proportional to the total dissipation $`\stackrel{~}{\eta }`$. Further, in contrast to the low-field limit, the conductivity in the MA regime is found to be a decreasing function of temperature and inversely related to the field.

Our starting point is the $`M`$-component Ginzburg-Landau free energy,

$`F[\psi ]`$ $`=`$ $`{\displaystyle }d^2r{\displaystyle }d\tau \{\left[(+{\displaystyle \frac{ie^{}}{\mathrm{}}}𝐀_0(𝐫))\psi _a^{}(𝐫,\tau )\right]`$ (3)
$`\left[\left({\displaystyle \frac{ie^{}}{\mathrm{}}}𝐀_0(𝐫)\right)\psi _a(𝐫,\tau )\right]+\left|_\tau \psi _a(𝐫,\tau )\right|^2`$
$`+\delta |\psi _a(𝐫,\tau )|^2+{\displaystyle \frac{U}{2M}}|\psi _a(𝐫,\tau )|^4\},`$
in the presence of a magnetic field, where $`𝐀_0(𝐫)=\{0,Hx,0\}`$ is the vector potential due to the applied field, $`e^{}=2e`$, and $`\delta `$ is the $`T=0`$ unrenormalized distance to the quantum critical point. For a time-independent magnetic field, it is convenient to expand the order parameter $`\psi _a(𝐫,\tau )`$ in the corresponding eigenfunctions that diagonalize the Gaussian part of the free-energy density

$`\psi _a(𝐫,\tau )={\displaystyle \underset{\omega _m,p_y,n}{}}e^{ip_yyi\omega _m\tau }\mathrm{\Psi }_n(x{\displaystyle \frac{\mathrm{}p_y}{e^{}H}})c_a(n,\omega _m,p_y)`$ (4)
where $`\mathrm{\Psi }_n(x\mathrm{}p_y/e^{}H)`$ is the solution of the corresponding eigenequation, expressible in Hermite polynomials, and $`c_a(n,\omega _m,p_y)`$ are the expansion coefficients. The related eigenvalues are given by $`\omega _H(n+1/2)`$ with $`\omega _H=2e^{}Ha/\mathrm{}=4\pi f/a`$, where $`a`$ is the lattice constant and $`f=\mathrm{\Phi }/\mathrm{\Phi }_0`$ is the magnetic frustration. Henceforth, we rescale $`H`$ such that $`\omega _HH`$. In the large $`M`$ (mean-field) treatment the bare $`T=0`$ distance to the critical point $`\delta `$ is renormalized, and the equation for the quasiparticle excitation gap $`m`$,

$`m^2=\delta +{\displaystyle \frac{H}{2}}+{\displaystyle \frac{UTH}{4\pi }}{\displaystyle \underset{\omega _m,n}{}}{\displaystyle \frac{1}{m^2+Hn+\omega _n^2}},`$ (5)
must be solved to obtain the effective Gaussian action. Performing the summation over the Matsubara frequencies exactly and employing the Poisson formula for the sum over ‘Landau’ levels $`n`$, we obtain (assuming $`m1`$)

$`\mathrm{\Delta }`$ $`=`$ $`2T\mathrm{ln}\left(2\mathrm{sinh}{\displaystyle \frac{m}{2T}}\right){\displaystyle \frac{H}{4m}}\mathrm{coth}{\displaystyle \frac{m}{2T}}{\displaystyle \frac{H}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dx}{e^{2\pi x}1}}`$ (7)
$`{\displaystyle \frac{a_{}\mathrm{sinh}(a_+/T)+a_+\mathrm{sin}(a_{}/T)}{\sqrt{m^4+H^2x^2}(\mathrm{sinh}^2(a_+/2T)+\mathrm{sin}^2(a_{}/2T))}},`$
where

$`a_\pm ={\displaystyle \frac{1}{\sqrt{2}}}\left[\sqrt{m^4+H^2x^2}\pm m^2\right]^{1/2}`$ (8)
and $`\mathrm{\Delta }=(4\pi /U)\delta +\mathrm{\Lambda }`$ ($`\mathrm{\Lambda }`$ is the upper momentum cutoff) is the renormalized distance to the zero-temperature critical point in the absence of the field.

Three limiting regimes follow from Eq. (7). 1) If $`T,\sqrt{H}\mathrm{\Delta }`$, then $`m=\mathrm{\Delta }`$ and we are in the quantum disordered regime (QD). This regime has been studied extensively recently, the key surprise being the onset of a Bose metal phase once the Cooper pairs loose their phase coherence. 2) The regime $`\sqrt{H},\mathrm{\Delta }T`$ is quantum critical (QC) in which $`m=\mathrm{\Theta }T`$ ($`\mathrm{\Theta }=2\mathrm{ln}((\sqrt{5}+1)/2)`$). Extensive studies of this regime exist in the literature. In both of these regimes, the magnetic field effects are subdominant, leaving the conclusions of the previous analyses in tact. 3) Finally, when $`\sqrt{H}T,\mathrm{\Delta }`$, the gap $`m=𝒜\sqrt{H}`$ emerges, where $`𝒜=0.418945`$ is the solution to the equation

$`𝒜={\displaystyle \frac{1}{4𝒜}}+\sqrt{{\displaystyle \frac{1}{2}}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dx}{e^{2\pi x}1}}\left[{\displaystyle \frac{\sqrt{𝒜^4+x^2}𝒜^2}{𝒜^4+x^2}}\right]^{1/2}.`$ (9)
In this regime, which we refer to as magnetically activated (MA), the correlation length, $`m^1`$, exhibits universal behavior and is determined (in the leading approximation) solely by the magnetic field $`H`$. The energies of bosonic quasiparticles, which represent the low-energy excitations close to the transition point, are quantized and occupy the degenerate ‘Landau levels’ with energies $`ϵ_n=\sqrt{Hn+m^2}`$. The $`1/M`$ correction is expected to modify the universal value of $`𝒜`$. Its calculation is prohibitively difficult because the coefficient in the quartic term of the GL action is a complicated function of the energy levels $`ϵ_n`$. All of these regimes imply that the elementary excitations are spin-wave like. This takes place only if $`\mathrm{\Delta }>0`$. In the opposite case, close enough to the finite $`T`$, finite $`H`$ transition, transport is determined by the motion of vortices as illustrated schematically in Fig. (1). The finite-size scaling considerations suggest that the correlation length near the transition is a universal function of two arguments: $`\xi =T^1f(\mathrm{\Delta }/T,\sqrt{H}/T)`$ and assumes the appropriate limiting forms if one of the parameters tends to zero. The $`M=\mathrm{}`$ form of $`m=\xi ^1`$ as a function of $`\sqrt{H/T}`$ obtained from the numerical solution of Eq. (7) for different values of $`\mathrm{\Delta }/T`$ is presented on Fig. (2). As is evident, a crossover occurs to $`\sqrt{H}`$ dependence as the field increases.

The linear conductivity (per flavor) near the IST can be obtained from the GL action with the help of the Kubo formula

$`\sigma _{\alpha \beta }(i\omega _n)={\displaystyle \frac{\mathrm{}}{\omega _n}}{\displaystyle d^2r𝑑\tau \frac{\delta ^2\mathrm{ln}Z}{\delta 𝐀_\alpha (\tau ,𝐫)\delta 𝐀_\beta (\mathrm{𝟎})}e^{i\omega _n\tau }}.`$ (10)
While calculating the averages in the partition function we expand $`\psi `$ using Eq. (4) which yields for the longitudinal conductivity:

$`\sigma (i\omega _\nu )`$ $`=`$ $`{\displaystyle \frac{(e^{})^2}{2h\omega _\nu }}TH^2{\displaystyle \underset{\omega _n}{}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(n+1)[2G(\omega _m,n)`$ (13)
$`G(\omega _m,n+1)G(\omega _m,n)G(\omega _m+\omega _\nu ,n+1)`$
$`G(\omega _m,n+1)G(\omega _m+\omega _\nu ,n)]`$
The experimentally observable finite temperature conductivity must be calculated once dissipative mechanisms are taken into account . Weak dissipation, $`\stackrel{~}{\eta }=\eta +1/\tau `$, is accounted for by introducing the term $`\stackrel{~}{\eta }|\omega _m||\psi (𝐫,\omega _m)|^2`$ into the free-energy density. Its inclusion changes the analytical properties of the Green functions in Eq. (13). Performing the analytical continuation, one finds that the static conductivity is given by

$`\sigma `$ $`=`$ $`{\displaystyle \frac{(e^{})^2H^2}{4\pi h}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(n+1)`$ (16)
$`\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{\mathrm{sinh}^2x}}{\displaystyle \frac{8\stackrel{~}{\eta }^2T^2x^2}{(ϵ_n4T^2x^2)^2+4\stackrel{~}{\eta }^2T^2x^2}}`$
$`\times {\displaystyle \frac{1}{(ϵ_{n+1}4T^2x^2)^2+4\stackrel{~}{\eta }^2T^2x^2}}.`$
The analysis of Eq. (16) simplifies if the dissipation is weak compared to the quasiparticle energies, implying $`\stackrel{~}{\eta }m`$. Asymptotic analytical expressions are readily obtained for $`mT`$, which define QD and MA regimes. In these cases, there are two competitive contributions to the integral over $`x`$. One comes from the minima of the denominator at $`x_n^2(ϵ_n^2\stackrel{~}{\eta }^2/2)/4T^2`$ and $`x_{n+1}^2(ϵ_{n+1}^2\stackrel{~}{\eta }^2/2)/4T^2`$ and the second from the region near $`x=0`$. Performing the approximate integration over $`x`$, and analyzing subsequently the sums over $`n`$, we arrive at the following limiting cases.

In the QD regime, $`m=\mathrm{\Delta }\sqrt{H}`$, and the conductivity is given by $`\sigma =\sigma ^{(1)}+\sigma ^{(2)}`$, where

$`\sigma ^{(1)}={\displaystyle \frac{4e^2}{\pi h}}\left({\displaystyle \frac{\pi \stackrel{~}{\eta }T}{3\mathrm{\Delta }^2}}\right)^2`$ (17)
and

$`\sigma ^{(2)}={\displaystyle \frac{4e^2}{\pi h}}\{\begin{array}{cc}{\displaystyle \frac{2\pi T}{\stackrel{~}{\eta }}}e^{\mathrm{\Delta }/T}\hfill & H\eta \mathrm{\Delta }\hfill \\ {\displaystyle \frac{8\eta T\mathrm{\Delta }^2}{H^2}}e^{\mathrm{\Delta }/T}\hfill & \stackrel{~}{\eta }\mathrm{\Delta }H\mathrm{\Delta }T\hfill \\ {\displaystyle \frac{\pi \stackrel{~}{\eta }}{T}}e^{\mathrm{\Delta }/T}\hfill & \mathrm{\Delta }TH\mathrm{\Delta }^2\hfill \end{array}.`$ (21)
If we are in the MA regime with $`\sqrt{H}m`$,

$`\sigma ={\displaystyle \frac{4e^2}{h}}\left[{\displaystyle \frac{\stackrel{~}{\eta }}{T}}e^{𝒜\sqrt{H}/T}+{\displaystyle \frac{2\pi }{3}}J\left({\displaystyle \frac{\stackrel{~}{\eta }T}{H}}\right)^2\right].`$ (22)
In the formula above,

$`J={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{n+1}{(n+1+𝒜^2)^2(n+𝒜^2)^2}}=23.913332`$ (23)
is a rather large numerical prefactor, arising from the $`M=\mathrm{}`$ limit. These expressions describe qualitatively the crossover from the QD to MA regime, provided the quasiparticles are well defined– that is, dissipation is sufficiently weak. As the field increases, the conductivity becomes directly proportional to dissipation, as dictated by Eq. (21). This can be attributed to the splitting of the minima in the denominator of the general formula, Eq. (16).

There are two sources of dissipation for which we must account. Internal dissipation arises from mutual scattering of quasiparticles, leading to a finite quasiparticle scattering rate $`1/\tau `$. External Ohmic dissipation $`\eta `$ may come from the resistive shunting of superconducting grains in a system of fabricated Josephson-junction arrays (JJA). The coefficient $`\eta `$ is then inversely proportional to the resistivity of a shunting resistor, $`R_g`$, which may be itself a function of temperature and magnetic field. In addition to Ohmic dissipation, disorder can also serve as a dissipative channel as has been shown for static disorder in the presence of d-wave pairing . In the QC regime, interactions alone lead to a large inverse scattering time, $`1/\tau Tm`$. Hence, quasiparticles are poorly defined, and small $`\epsilon =3d`$ or $`1/M`$ expansions are necessary to determine the conductivity in a controlled way. The implementation of this procedure leads to the universal value for the conductivity close to $`\sigma _Q=4e^2/h`$. The weak external dissipation $`\eta mT`$ is clearly subdominant in this regime, except at low temperatures. On the contrary, in the QD regime, $`1/\tau Te^{\mathrm{\Delta }/T}`$, (up to logarithmic corrections), as a result of the Boltzmann statistics of the quasiparticles. In the absence of Ohmic dissipation, as follows from Eq. (17)- (21), metallic, rather than insulating behavior of the conductivity obtains provided that $`\sqrt{H}\mathrm{\Delta }`$. The persistence of the Bose metal phase for a finite range of magnetic field is possibly relevant to the recent experiments of Mason and Kapitulnik.

In the MA regime, quasiparticles also obey Boltzmann statistics. Based on dimensionality arguments similar to those valid for the QD regime, we suggest that in the MA regime, $`1/\tau Te^{𝒜\sqrt{H}/T}`$, where the constant of proportionality is determined by a numerical coefficient that is of the first order in $`1/M`$. The calculation of the corresponding lifetime is complicated because of the discrete character of the quasiparticle spectrum. However, as follows from Eq. (21), the conductivity becomes directly, rather than inversely, proportional to dissipation for sufficiently large $`H`$. This means that inclusion of external dissipation $`\eta `$ and determines only how fast the conductivity decreases with increasing field. To illustrate this, we resort to a semi-phenomenological calculation in which we set $`1/\tau =2.3300502Te^{m/T}`$ and $`\eta =0.1T`$. The interpolation formula for $`1/\tau `$ is chosen to satisfy the behavior in all regimes near the critical point. The numerical prefactor in $`1/\tau `$ is taken in such a way, that for $`\eta =0`$, $`\sigma \sigma _Q`$ in the QC regime. The conductivity was then calculated numerically as a function of $`\sqrt{H}/T`$ for different values of $`\mathrm{\Delta }/T`$, using Eqs. (7) and (16). The results that are depicted in Fig. (2) support rather rapid decrease of $`\sigma `$ with increasing $`\sqrt{H}/T`$ but their dependence on $`\eta /T`$ is relatively weak. We conclude, then, based on the mean-field treatment, that the conductivity decreases when the crossover to the MA regime occurs even if $`\eta =0`$. This behavior is in agreement with experiments on the magnetic field-tuned IST in disordered indium-oxide thin films. However, these results are in disagreement with the suggestion of Kim and Weichman that the conductivity saturates with increasing field. This discrepancy can be explained by the fact that for $`\eta =0`$ in the general universal scaling function,

$`\sigma =(4e^2/h)\mathrm{\Sigma }(\omega /T,\mathrm{\Delta }^z/T,\sqrt{H}/\mathrm{\Delta }),`$ (24)
applicable for the conductivity near the $`d=2`$ IST, the limits $`T=0`$, $`\omega 0`$ and $`\omega =0`$, $`T0`$ do not commute. The former is the collisionless limit, apparently considered in Ref. (14), in which the conductivity is likely to experience a crossover between two universal values depending on the ratio $`\sqrt{H}/\mathrm{\Delta }`$. However, it is the collision-dominated limit that is experimentally relevant and hence the main subject of our treatment. It is straightforward to verify that Eq.(16) satisfies the general scaling relation (24) with $`\omega =0`$ and $`z=1`$, provided that $`\eta =0`$ and the inverse scattering time $`1/\tau `$ obeys itself the scaling form $`1/\tau =Tf(\mathrm{\Delta }^z/T,\sqrt{H}/\mathrm{\Delta })`$ in the critical region.

The results above were obtained for low enough magnetic fields, $`\sqrt{H}f1`$ and hence, magnetic frustration can be neglected. In translationally periodic systems such as JJA, magnetic frustration plays an important role giving rise to peaks in the conductivity for integer values of frustration $`f`$. Moreover, the peaks corresponding to the rational values $`f=\frac{1}{2},\frac{1}{3},\frac{3}{8}\mathrm{}`$ may also be observed. However, disorder disrupts the lattice periodicity smearing the peaks in the conductivity giving rise to a disappearance of frustration effects. This situation takes place in the disorder-tuned IST in thin films. Hence, we expect, that as long as Cooper pairs remain intact and the transition is described by Eq.(3), our calculations provide a qualitative picture for the behavior of the magneto-conductivity near the point of nominal the IST in thin films including the observation of a Bose metal phase.

This work was funded by the DMR of the NSF grant No. DMR98-96134.