format: text

text:
# References

Rearrangement of the Fermi Surface of Dense Neutron Matter
and Direct Urca Cooling of Neutron Stars

D. N. Voskresensky1

Moscow Institute for Physics and Engineering, Moscow 115409, Russia

V. A. Khodel2 and M. V. Zverev

Russian Research Center Kurchatov Institute, Moscow 123182, Russia

and

J. W. Clark

McDonnell Center for the Space Sciences and Department of Physics, Washington University,
St. Louis, MO 63130 USA

1Gesellschaft für Schwerionenforschung GSI, D-64220 Darmstadt, Germany

2McDonnell Center for the Space Sciences and Department of Physics, Washington University,
St. Louis, MO 63130 USA

## Abstract

It is proposed that a rearrangement of single-particle degrees of freedom may occur in a portion of the quantum fluid interior of a neutron star. Such a rearrangement is associated with the pronounced softening of the spin-isospin collective mode which, under increasing density, leads to pion condensation. Arguments and estimates based on fundamental relations of many-body theory show that one realization of this phenomenon could produce very rapid cooling of the star via a direct nucleon Urca process displaying a $`T^5`$ dependence on temperature.

1. INTRODUCTION

The EINSTEIN, EXOSAT, and ROSAT orbiting X-ray observatories have measured surface temperatures of certain neutron stars and set upper limits on surface temperatures of others (). The data for the supernova remnants in 3C58, the Crab, and RCW103 indicate relatively slow cooling, while that for Vela, PSR2334+61, PSR0656+14, and Geminga point to substantially more rapid cooling. In the so-called standard scenario for neutron-star cooling, the primary role is played by the modified Urca process $`(nnnpe^{}\overline{\nu }_e;npe^{}nn\nu _e)`$, first considered by Bahcall & Wolf (1965) and later reexamined by Friman & Maxwell (1979) in terms of an in-vacuum one-pion exchange model. Cooling simulations based on the results of these works play out the slow scenario of thermal evolution and fail to explain the rapid cooling of some stars. This picture is profoundly altered when in-medium effects are taken into account. Modification of $`NN`$, $`\pi N`$, and $`KN`$ interactions with increasing density may be so strong that pion () and kaon (\[Brown 1994\]) condensates form in the interior region of a high-mass neutron star, leading to a dramatic increase of the neutrino luminosity ().

Here we shall demonstrate that softening of the spin-isospin (pion) collective mode in dense neutron matter as predicted by Migdal (1978) could give rise to a rearrangement of single-particle degrees of freedom prior to the onset of pion condensation and open a new channel of neutrino cooling of neutron stars, by giving access to the direct Urca process. This process does not involve a neutron (or other) spectator and, if allowed, is an extremely efficient cooling mechanism ().

2. BUBBLE REARRANGEMENT OF THE FERMI SPHERE

A rearrangement of single-particle degrees of freedom takes place if the necessary condition for stability of the normal state of a Fermi liquid is violated. At $`T=0`$ this condition requires that the change of the ground-state energy $`E_0`$ remain positive for any admissible variation $`\delta n(p)`$ of the Landau quasiparticle distribution $`n(p)`$ away from the normal-state step-function distribution $`\theta (pp_F)`$. Formally,

$$\delta E_0=\xi (p,n(p))\delta n(p)\frac{\mathrm{d}^3p}{(2\pi )^3}>0,$$
(1)
where $`\xi (p,n(p))\epsilon (p,n(p))\mu `$ is the energy of a quasiparticle relative to the chemical potential $`\mu `$. The condition (1) fails if a depression with $`\xi <0`$ forms in the spectrum $`\xi (p)`$ at $`p>p_F`$; it likewise fails if there arises an elevation with $`\xi >0`$ at $`p<p_F`$. The rearrangement is precipitated when the density $`\rho `$ reaches a critical value $`\rho _{cF}`$ at which there emerges a bifurcation and a new root $`p=p_0`$ of the relation

$$\xi (p,n(p);\rho _{cF})=0,$$
(2)
which ordinarily serves merely to specify the Fermi momentum $`p_F`$.

The simplest kind of rearrangement of the momentum distribution $`n(p)`$ of quasiparticles of given spin projection retains the property that its values are restricted to 0 and 1, but the Fermi sea becomes doubly connected (). At densities exceeding the critical value $`\rho _{cF}`$, the normal-state distribution $`\theta (pp_F)`$ is modified by the presence of a “bubble,” or absence of particles, over the range $`p_i<p<p_f<p_F`$, with the inner surface $`p_i`$ located relatively close to the origin and the Fermi momentum $`p_F`$ readjusted to maintain the prescribed neutron density. The distance $`p_fp_i`$ between the two new Fermi surfaces lying interior to $`p_F`$ can be estimated using the formula $`\xi (pp_0,n(p);\rho )=\xi _0(\rho \rho _{cF})A(pp_0)^2`$, where $`A`$ and $`\xi _0`$ are positive constants. This formula embodies the essential properties that $`\xi (p)`$ is negative for any $`p<p_F`$ at $`\rho <\rho _{cF}`$ and that its maximum value first reaches zero at $`\rho =\rho _{cF}`$. Employing this parametrization in the relation (2), it is found that no bifurcation point exists for $`\rho <\rho _{cF}`$, whereas two solutions arise for $`\rho >\rho _{cF}`$, with the distance between the two new Fermi surfaces growing in proportion to $`\sqrt{\rho \rho _{cF}}`$.

Were such a rearrangement to occur in the neutron subsystem of neutron-star matter, the emergence of new Fermi surfaces situated at lower momenta would permit the direct nucleon Urca process to operate at a much lower density than hitherto considered possible (). In this process, the tandem reactions $`npe^{}\overline{\nu }_e`$ and $`pe^{}n\nu _e`$ are driven by thermal excitations. The condition of high degeneracy prevailing even in young neutron stars implies that these excitations remain close to the Fermi surfaces of the participants. For the direct Urca mechanism to contribute appreciably to neutron-star cooling, momentum conservation then demands satisfaction of the triangle inequalities $`|p_{Fp}p_{Fe}|p_{Fn}p_{Fp}+p_{Fe}`$ among the proton, electron, and neutron Fermi momenta. With the conventional singly-connected neutron Fermi sphere, these conditions are met only at very high baryon densities where the proton Fermi momentum $`k_{Fp}`$ reaches sufficiently large values. The best current estimates yield a threshold baryon density of at least five times the saturation density $`\rho _0`$ of symmetrical nuclear matter (). On the other hand, if the neutron subsystem undergoes a rearrangement that allows for thermal excitations of neutron quasiparticles at much lower momenta than $`p_{Fn}`$, specifically at $`p_i`$ and $`p_f`$ in the bubble rearrangement scenario, the triangle inequalities are much more easily satisfied and the direct Urca process is greatly facilitated.

We next provide a quantitative basis for this qualitative idea, by appealing to established methods of microscopic many-body theory. For a bifurcation point to arise in the solution of equation (2), both the spectrum $`\xi (p)`$ and the scalar component $`f(𝐩_1,𝐩_2;𝐤=0)`$ of the Landau amplitude of the quasiparticle interaction must depend strongly on momentum. The connection (, p. 37)

$$\frac{\xi _(p)}{𝐩}=\frac{\xi ^0(p)}{𝐩}+\frac{1}{2}\mathrm{Tr}_\sigma \mathrm{Tr}_{\sigma _1}f(𝐩,\sigma ,𝐩_1;\sigma _1,𝐤=0)\frac{n(p_1)}{𝐩_1}\frac{\mathrm{d}^3p_1}{(2\pi )^3}$$
(3)
between these two quantities, wherein $`\xi _0(p)`$ is the free spectrum and spin dependence is made explicit, is usually employed to find the effective mass $`M^{}`$ in terms of the first harmonic of the expansion of the amplitude $`f(\widehat{𝐩}_1\widehat{𝐩}_2;0)`$ in Legendre polynomials. We exploit this relation in the search for new roots of equation (2).

Within the integral in (3), the quasiparticle interaction amplitude $`f(𝐩_1,𝐩_2;𝐤=0)`$ may be replaced by $`(M^{}/M)\mathrm{\Gamma }^k(𝐩_1,𝐩_2)`$, the scattering amplitude $`\mathrm{\Gamma }^k(𝐩_1,𝐩_2)\mathrm{\Gamma }(𝐩_1,𝐩_2;𝐤=0)`$ modified by an effective-mass factor (). The requisite strong momentum dependences can arise if the system approaches a second-order phase transition that occurs at some critical density $`\rho _c`$ where a collective mode of frequency $`\omega _s(k)`$ collapses at wave number $`k=k_0`$ and the corresponding susceptibility diverges along with the scattering amplitude. Near such a soft-mode critical point, the singular part $`\mathrm{\Gamma }^s`$ of the scattering amplitude $`\mathrm{\Gamma }(𝐩_1,𝐩_2;𝐤)`$ at momentum transfer $`𝐤`$ may be expressed quite generally by ()

$$\frac{M^{}}{M}\mathrm{\Gamma }_{\alpha \kappa ;\beta \lambda }^s(𝐩_1,𝐩_2;𝐤,\rho \rho _c)=O_{\alpha \kappa }O_{\beta \lambda }D(k)+O_{\alpha \lambda }O_{\beta \kappa }D(|𝐩_1𝐩_2+𝐤|),$$
(4)
where $`O`$ denotes the vertex determining the structure of the collective-mode operator (e.g. $`O=(\sigma 𝐧)\tau `$ for the spin-isospin mode with pion quantum numbers, where $`𝐧`$ is a unit vector along the relevant momentum). In deriving (4), antisymmetry of the two-particle wave function under interchange of the coordinates and spins of the two particles has been invoked. The propagator $`D`$ may be parametrized as $`D(k)=\left[\beta ^2+\gamma ^2(k^2/k_0^21)^2\right]^1`$, where $`\beta (\rho )`$ measures the proximity to the phase transition point, with $`\beta (\rho \rho _c)\rho _c\rho `$ (cf. Dyugaev (1976)).

The essential messages of the preceding development are that the singular part $`\mathrm{\Gamma }^s(𝐩_1,𝐩_2;𝐤=0)D(|𝐩_1𝐩_2|)`$ of the scattering amplitude depends on the difference $`𝐩_1𝐩_2`$ and that as one approaches the soft-mode phase transition point this dependence becomes quite strong. We assume that the remaining contributions to $`\mathrm{\Gamma }(𝐩_1,𝐩_2;𝐤=0)`$ can be adequately incorporated by renormalization of the chemical potential $`\mu `$. Equation (3) is then easily integrated to produce an explicit expression suitable for calculation of the single-particle spectrum,

$$\xi (p)=\xi ^0(p)+\frac{1}{2}\mathrm{Tr}(O_{\alpha \lambda }O_{\lambda \alpha })D(|𝐩𝐩_1|)n(p_1)\frac{\mathrm{d}^3p_1}{(2\pi )^3}.$$
(5)
We now apply this equation to dense neutron matter in the vicinity of the second-order phase transition associated with neutral pion condensation, which is engendered by the softening of the spin-isospin mode having $`\pi ^0`$ quantum numbers (). It has been predicted that the collapse of this mode will take place at a neutron density $`\rho _c=\rho _{c\pi }`$ in the range (0.2 – 0.5) fm-3 (roughly 1–3 times $`\rho _0`$), depending on theoretical assumptions (). Unfortunately, there is as yet no definitive microscopic treatment of neutron-star matter from which one can extract or derive quantitatively reliable values for the input parameters $`\beta `$, $`\gamma `$, and $`k_0`$ of our model.

In this situation, a reasonable strategy is to perform calculations based on expression (5) for several choices of the parameters of the microscopic model. Substituting (5) into relation (2), one finds the critical density $`\rho _{cF}`$ for the onset of a bifurcation of the latter equation. For $`\rho >\rho _{cF}`$ this equation then determines two new momenta $`p_i`$ and $`p_f`$ where $`\xi (p)`$ vanishes, which delimit the bubble region of $`n(p)`$ and between which $`\xi (p)`$ is positive. Representative numerical results for the spectrum $`\xi (p)`$ are plotted in Fig. 1. Results for the phase diagram of dense neutron matter in the $`\rho /\rho _0`$ versus $`\beta ^2/m_\pi ^2`$ plane are displayed in Fig. 2. Different values of $`\gamma `$ are considered, while keeping the parameter $`k_0`$ fixed at the value $`0.9p_{Fn}`$ suggested by earlier numerical calculations ().

It is evident that variation of the parameters $`\beta `$, $`\gamma `$, and $`k_0`$ within sensible bounds can have strong effects on the phase diagram and therefore on the extent of the phase with rearranged quasiparticle occupation. Nevertheless, our numerical study has documented three salient features of the bubble rearrangement. First, the critical density $`\rho _{cF}`$ for the rearrangement is less than the critical density $`\rho _{c\pi }`$ for pion condensation. Since both phenomena stem from the strong momentum dependence of the Landau amplitude $`f(𝐩_1,𝐩_2;𝐤0)`$, rearrangement of the quasiparticle distribution may be regarded as a precursor of pion condensation. Second, the bifurcation point corresponding to formation of a hole bubble in the neutron momentum distribution is positioned at small momenta, $`p_0<0.2p_F`$, irrespective of the applicable value of $`\rho _{c\pi }`$. Third, the spectrum $`\xi (p)`$ shows a deep depression for $`p(0.50.6)p_F`$. And fourth, the ratios $`\rho _{cF}/\rho _{c\pi }`$ and $`p_0/p_F`$ are insensitive to the actual value taken by $`\rho _{c\pi }`$ within the usual range of theoretical predictions.

Analogous considerations apply to the proton subsystem of neutron-star matter, in which case one is dealing with the charged pion mode. Estimates () indicate that this mode is also softened in dense matter, with a critical density not far from that for neutral pion condensation. One may then argue that, under the influence of the strongly momentum-dependent external field provided by the neutron medium, protons will leave the old Fermi sphere and occupy states of relatively large momentum, $`p0.5p_{Fn}`$. The impact of this further rearrangement on the proton-neutron ratio and on the rate of neutrino cooling requires a separate analysis, which we defer.

Having laid the microscopic basis for a rearrangement of the neutron Fermi surface that creates a bubble at low momenta in the Fermi sea, we return to its most striking astrophysical implication. Beyond the bifurcation point, the triangle inequalities can now be satisfied without the conventional requirement ()) that the proton fraction exceed some 11–14%. Accordingly, the direct Urca process becomes active in the density regime just short of the threshold for pion condensation. At temperatures $`T`$ above the anticipated superfluid phase transition, the $`T`$-dependence of the resulting neutrino emissivity is determined through the usual expression

$`ϵ(npe^{}\overline{\nu })`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}2G_F^2(1+3g_A^2){\displaystyle \underset{i}{}}\epsilon _\nu n_n(1n_p)(1n_e)`$ (6)
$`\times \delta (𝐩_n𝐩_p𝐩_e𝐩_\nu )\delta (\epsilon _n\epsilon _p\epsilon _e\epsilon _\nu ),`$
derived in quasiparticle approximation in terms of the occupations $`n_i`$ and energies $`\epsilon _i`$ of the reacting particles. We do not consider a renormalization of the coupling constants $`G_F^2`$ and $`g_A^2`$ due to medium effects, and the emissivity takes the customary form ()

$$ϵ(npe^{}\overline{\nu })1.2\times 10^{27}\frac{M_n^{}}{M_n}\frac{M_p^{}}{M_p}\left(\frac{\mu _e}{100\text{MeV}}\right)T_9^6\text{erg}\text{cm}^3\text{sec}^1,$$
(7)
where $`\mu _e`$ is the electron chemical potential and the temperature $`T`$ is measured in multiples of $`10^9`$ K. Since the neutron and proton effective masses $`M_n^{}`$ and $`M_p^{}`$ remain $`T`$-independent, the bubble rearrangement serves only to turn on the direct Urca process at a lower density, without altering the familiar $`T^6`$ power-law behavior of the emissivity.

3. FERMION CONDENSATION

Interestingly enough, there exists a more radical scenario for rearrangement of the quasiparticle distribution, known as fermion condensation (). In this case, the occupancy $`n(p)`$ may be partial, i.e., it may lie between 0 and 1. At $`T=0`$, the new quasiparticle distribution $`n(p)`$ is to be found from the variational condition

$$\frac{\delta E_0[n(p)]}{\delta n(p)}=\mu ,p\mathrm{\Omega }.$$
(8)
The left-hand side of equation (8) is just the quasiparticle energy $`\epsilon (p)`$. Hence its coincidence with the chemical potential $`\mu `$ means that the spectrum of single-particle excitations will be dispersionless (i.e., the quasiparticle group velocity will vanish) throughout the entire momentum domain $`\mathrm{\Omega }`$ where $`\xi (p)=\epsilon (p)\mu =0`$, and not merely at isolated points as implied by equation (2). The family of quasiparticles having momenta $`p\mathrm{\Omega }`$ is called the fermion condensate because of a conspicuous analogy with the low-temperature Bose gas, in which the energy of condensate particles is also equal to the chemical potential $`\mu `$. A key signature of fermion condensation has been observed in strongly correlated electron systems (): flat portions of the single-particle spectrum have been seen experimentally in a number of high-temperature superconductors ().

The two types of rearrangement – bubble formation and fermion condensation – can compete with each other in the density regime just below the soft-mode phase transition. Numerical studies () demonstrate that fermion condensation wins the contest at nonzero $`T`$. Let us suppose this is the case in the neutron-star medium, while continuing to disregard nucleonic pairing phenomena. Consistency of the Fermi-Dirac form $`n(p,T)=\left\{\mathrm{exp}\left[\xi (p;n(p,T))/T\right]+1\right\}^1`$ for the quasiparticle distribution with the variational condition (8) requires that the spectrum $`\xi (p,T)`$ of the fermion-condensate phase grows linearly with $`T`$ at low temperature, implying an effective mass inversely proportional to $`T`$ ().

Although the rearranged momentum distribution derived from equation (8) differs from the bubble configuration, its structure will also admit thermal excitations at low neutron momenta. Hence we may again expect most neutron stars to contain a region of relatively moderate density, bounded below by $`\rho _{cF}`$ and above by $`\rho _{c\pi }`$, in which the direct Urca process operates vigorously. However, due to the new feature of a $`T`$-dependent neutron effective mass, $`M_n1/T`$, we may anticipate an enhancement of the neutrino emissivity relative to the standard result (), corresponding to a $`T^5`$ rather than a $`T^6`$ dependence on the temperature.

4. CONCLUSIONS

We have explored the possibility that an effective interaction with strong momentum dependence gives rise to a rearrangement of the neutron momentum distribution in neutron-star matter. Two plausible manifestations of this phase transformation – creation of a doubly-connected Fermi surface and fermion condensation – have been considered. Both open the prospect that direct nucleon Urca cooling is present in a density regime just below the threshold for pion condensation and consequently at a density much lower than previously estimated. If a fermion condensate is formed, the resulting neutrino emissivity is significantly larger than that generated by the direct Urca process in normal matter. Within the affected density range, it would therefore dominate all other proposed neutrino cooling mechanisms (). Future studies along this line will focus on temperatures below the superfluid transition and on the effect of the dramatically increased emissivity on neutrino opacity.

This research was supported by NSF Grants PHY-9602127 and PHY-9900713 (JWC and VAK) and by the McDonnell Center for the Space Sciences (VAK). We thank M. Baldo, M. Di Toro, and E. E. Kolomeitsev for fruitful discussions. DNV expresses his appreciation for hospitality and support provided by GSI Darmstadt. MVZ acknowledges the hospitality of INFN (Sezione di Catania).

FIGURE CAPTIONS

Fig. 1. The dimensionless neutron spectrum $`y_n(p)=\xi _n(p)/(p_F^2/2M)`$ at the critical densities $`\rho _{cF}`$ corresponding to three different sets of model parameters: (a) $`\gamma =1.25m_\pi `$, $`k_0=0.9p_{Fn}`$, $`\beta ^2=0.22m_\pi ^2`$ ($`\rho _{cF}1.19\rho _0)`$, (b) $`\gamma =1.25m_\pi `$, $`k_0=0.9p_{Fn}`$, $`\beta ^2=0.25m_\pi ^2`$ ($`\rho _{cF}1.76\rho _0)`$, (c) $`\gamma =1.25m_\pi `$, $`k_0=p_{Fn}`$, $`\beta ^2=0.13m_\pi ^2`$ ($`\rho _{cF}1.88\rho _0)`$. Two different positions of the bifurcation point, namely $`p_0=0`$ (for parameter sets (a) and (b)) and $`p_00.12p_{Fn}`$ (for set (c)), are indicated by arrows.

Fig. 2. Phase diagram of neutron matter in the variables $`\rho `$ (measured in $`\rho _0`$) and $`\beta ^2`$ (measured in $`m_\pi ^2`$), as calculated for $`k_0=0.9p_{Fn}`$ and four different values of $`\gamma `$, which (in $`m_\pi `$ units) label the corresponding the phase boundaries separating the bubble phase (upper left) from the normal phase (lower right).


text:
# 1 Introduction

## 1 Introduction

Accretion disks are inferred to exist in objects of very different scales: km to millions of km in low Mass X-ray Binaries (LMXB) and Cataclysmic Variables (CV), solar radius-to-AU scale in protostellar disks, and AU-to-parsec scales for the disks in Active Galactic Nuclei (AGN).

An interesting observational connection exists between accretion disks and jets (such as the spectacular jets from AGN and protostars), and outflows (the ‘CO-outflows’ from protostars and possibly the ‘broad-line-regions’ in AGN). Lacking direct (i.e. spatially resolved) observations of disks, theory has tried to provide models, with varying degrees of success. Uncertainty still exists with respect to some basic questions. In this situation, progress made by observations or modeling of a particular class of objects is likely to have direct impact for the understanding of other objects, including the enigmatic connection with jets.

In this lecture I concentrate on the more basic aspects of accretion disks, but an attempt is made to mention topics of current interest, such as magnetic viscosity, as well. Emphasis is on those aspects of accretion disk theory that connect to the observations of LMXB and CV’s. For other reviews on the basics of accretion disks, see Pringle (1981), Treves et al. (1988). For a more in-depth treatment, see the textbook by Frank et al. (1992).

## 2 Accretion: general

Gas falling into a point mass potential

$$\mathrm{\Phi }=\frac{GM}{r}$$

from a distance $`r_0`$ to a distance $`r`$ converts gravitational into kinetic energy, by an amount $`\mathrm{\Delta }\mathrm{\Phi }=GM(1/r1/r_0)`$. For simplicity, assuming that the starting distance is large, $`\mathrm{\Delta }\mathrm{\Phi }=GM/r`$. If the gas is then brought to rest, for example at the surface of a star, the amount of energy $`e`$ dissipated per unit mass is

$$e=\frac{GM}{r}(\mathrm{r}est)$$

or, if it goes into a circular Kelper orbit at distance $`r`$:

$$e=\frac{1}{2}\frac{GM}{r}(\mathrm{o}rbit).$$

The dissipated energy goes into internal energy of the gas, and into radiation which escapes to infinity (usually in the form of photons, but neutrino losses can also play a role).

### 2.1 Adiabatic accretion

Consider first the case when radiation losses are neglected. This is adiabatic accretion. For an ideal gas with constant ratio of specific heats $`\gamma `$, the internal energy per unit mass is

$$e=\frac{P}{(\gamma 1)\rho }.$$

With the equation of state

$$P=\rho T/\mu $$
(1)
where $``$ is the gas constant, $`\mu `$ the mean atomic weight per particle, we find the temperature of the gas after the dissipation has taken place (assuming that the gas goes into a circular orbit):

$$T=\frac{1}{2}(\gamma 1)T_{\mathrm{vir}},$$
(2)
where $`T_{\mathrm{vir}}`$, the virial temperature is given by

$$T_{\mathrm{vir}}=\frac{GM\mu }{r}.$$

In an atmosphere with temperature near $`T_{\mathrm{vir}}`$, the sound speed is close to the escape speed from the system, the hydrostatic pressure scale height is of the order of $`r`$, and such an atmosphere may evaporate on a relatively short time scale in the form of a stellar wind.

A simple example is spherical adiabatic accretion (Bondi, 1952). An important result is that such accretion is possible only if $`\gamma 5/3`$. The larger $`\gamma `$, the larger the temperature in the accreted gas (eq. 2), and beyond a critical value the temperature is too high for the gas to stay bound in the potential. A classical situation where adiabatic and roughly spherical accretion takes place is a supernova implosion: when the central temperature becomes high enough for the radiation field to start desintegrating nuclei, $`\gamma `$ drops and the envelope collapses onto the forming neutron star via a nearly static accretion shock. Another case are Thorne-Zytkow objects (e.g. Cannon et al. 1992), where $`\gamma `$ can drop to low values due to pair creation, initiating an adiabatic accretion onto the black hole.

Adiabatic spherical accretion is fast, taking place on the dynamical or free fall time scale

$$\tau _\mathrm{d}=r/v_\mathrm{K}=(r^3/GM)^{1/2},$$
(3)
where $`v_\mathrm{K}`$ is the Kepler orbital velocity.

When radiative loss becomes important, the accreting gas can stay cool irrespective of the value of $`\gamma `$, and Bondi’s critical value $`\gamma =5/3`$ plays no role. With such losses, the temperatures of accretion disks are usually much lower than the virial temperature. The optical depth of the accreting flow increases with the accretion rate $`\dot{M}`$. When the optical depth becomes large enough so that the photons are ‘trapped’ in the flow, the accretion just carries them in, together with the gas (Rees 1978, Begelman 1979). Above a certain critical rate $`\dot{M}_\mathrm{c}`$, accretion is therefore adiabatic.

### 2.2 The Eddington Limit

Objects of high luminosity have a tendency to blow their atmospheres away due to the radiative force exerted when the outward traveling photons are scattered or absorbed. Consider a volume of gas on which a flux of photons is incident from one side. Per gram of matter, the gas presents a scattering (or absorbing) surface area of $`\kappa `$ cm2. The force exerted by the radiative flux $`F`$ on one gram is $`F\kappa /c`$. The force of gravity pulling back on this one gram of mass is $`GM/r^2`$. The critical flux at which the two forces balance is

$$F_\mathrm{E}=\frac{c}{\kappa }\frac{GM}{r^2}$$
(4)
Assuming that the flux is spherically symmetric, this can be converted into a critical luminosity

$$L_\mathrm{E}=4\pi GMc/\kappa ,$$
(5)
the Eddington luminosity (e.g. Rybicki and Lightman, 1979). If the gas is fully ionized, its opacity is dominated by electron scattering, and for solar composition $`\kappa `$ is then of the order $`0.3`$ cm2/g (about a factor 2 lower for fully ionized helium, a factor up to $`10^3`$ higher for partially ionized gases). With these assumptions,

$$L_\mathrm{E}\mathrm{1.7\hspace{0.17em}10}^{38}\frac{M}{M_{}}\mathrm{e}rg/s\mathrm{4\hspace{0.17em}10}^4\frac{M}{M_{}}L_{}$$

If this luminosity results from accretion, it corresponds to the Eddington accretion rate $`\dot{M}_\mathrm{E}`$:

$$\frac{GM}{r}\dot{M}_\mathrm{E}=L_\mathrm{E}\dot{M}_\mathrm{E}=4\pi rc/\kappa .$$
(6)
Whereas $`L_\mathrm{E}`$ is a true limit that can not be exceeded by a static radiating object except by geometrical factors of order unity (see chapter 10 in Frank et al, 1992), no maximum exists on the accretion rate. For $`\dot{M}>\dot{M}_\mathrm{E}`$ the plasma is just swallowed whole, including the radiation energy in it (cf. discussion in the preceding section). With $`\kappa =0.3`$:

$$\dot{M}_\mathrm{E}\mathrm{1.3\hspace{0.17em}10}^{18}r_6\mathrm{g}/s\mathrm{2\hspace{0.17em}10}^8r_6M_{}\mathrm{y}r^1,$$

where $`r_6`$ is the radius of the accreting object in units of 10 km.

## 3 Accretion with Angular Momentum

When the accreting gas has a zonzero angular momentum with respect to the accreting object, it can not accrete directly. A new time scale appears, the time scale for outward transport of angular momentum. Since this is in general much longer than the dynamical time scale, much of what was said about spherical accretion needs modification for accretion with angular momentum.

Consider the accretion in a close binary consisting of a compact (white dwarf, neutron star or black hole) primary of mass $`M_1`$ and a main sequence companion of mass $`M_2`$. The mass ratio is defined as $`q=M_2/M_1`$ (note: $`q`$ is just as often defined the other way around).

If $`M_1`$ and $`M_2`$ orbit each other in a circular orbit and their separation is $`a`$, the orbital frequency $`\mathrm{\Omega }`$ is

$$\mathrm{\Omega }^2=G(M_1+M_2)/a^3.$$

The accretion process is most easily described in a coordinate frame that corotates with this orbit, and with its origin in the center of mass. Matter that is stationary in this frame experiences an effective potential, the Roche potential (Ch. 4 in Frank, King and Raine, 1992), given by

$$\varphi _\mathrm{R}(𝐫)=\frac{GM}{r_1}\frac{GM}{r_2}\frac{1}{2}\mathrm{\Omega }^2r^2$$
(7)
where $`r_{1,2}`$ are the distances of point $`𝐫`$ to stars 1,2. Matter that does not corotate experiences a very different force (due to the Coriolis force). The Roche potential is therefore useful only in a rather limited sense. For non-corotating gas intuition based on the Roche geometry is usually confusing. Keeping in mind this limitation, consider the equipotential surfaces of (7). The surfaces of stars $`M_{1,2}`$, assumed to corotate with the orbit, are equipotential surfaces of (7). Near the centers of mass (at low values of $`\varphi _\mathrm{R}`$) they are unaffected by the other star, at higher $`\mathrm{\Phi }`$ they are distorted and at a critical value $`\mathrm{\Phi }_1`$ the two parts of the surface touch. This is the critical Roche surface $`S_1`$ whose two parts are called the Roche lobes. Binaries lose angular momentum through gravitational radiation and a magnetic wind from the secondary (if it has a convective envelope). Through this loss the separation between the components decreases and both Roche lobes decrease in size. Mass transfer starts when $`M_2`$ fills its Roche lobe, and continues as long as the angular momentum loss from the system lasts. A stream of gas then flows through the point of contact of the two parts of $`S_1`$, the inner Lagrange point $`L_1`$. If the force acting on it were derivable entirely from (7) the gas would just fall in radially onto $`M_1`$. As soon as it moves however, it does not corotate any more and its orbit under the influence of the Coriolis force is different (Fig. 1).

Since the gas at $`L_1`$ is very cold compared with the virial temperature, its sound speed is small compared with the velocity it gets after only a small distance from $`L_1`$. The flow into the Roche lobe of $`M_1`$ is therefore highly supersonic. Such hypersonic flow is essentially ballistic, that is, the stream flows along the path taken by free particles.

Though the gas stream on the whole follows an orbit close to that of a free particle, a strong shock develops at the point where the orbit intersects itself. \[In practice shocks already develop shortly after passing the pericenter at $`M_1`$, when the gas is decelerated again. Supersonic flows that are decelerated by whatever means in general develop shocks (e.g. Courant and Friedrichs 1948, Massey, 1968). The effect can be seen in action in the movie published in Różyczka and Spruit, 1993\]. After this, the gas settles into a ring, into which the stream continues to feed mass. If the mass ratio $`q`$ is not too small this ring forms fairly close to $`M_1`$. An approximate value for its radius is found by noting that near $`M_1`$ the tidal force due to the secondary is small, so that the angular momentum of the gas with respect to $`M_1`$ is approximately conserved. If the gas continues to conserve angular momentum while dissipating energy, it settles into the minimum energy orbit with the specific angular momentum $`j`$ of the incoming stream. The radius of this orbit, the circularization radius $`r_\mathrm{c}`$ is determined from

$$(GM_1r_\mathrm{c})^{1/2}=j.$$

The value of $`j`$ is found by a simple integration of the orbit starting at $`L_1`$ and measuring $`j`$ at some point near pericenter. In units of the orbital separation $`a`$, $`r_\mathrm{c}`$ and the distance $`r_{\mathrm{L1}}`$ from $`M_1`$ to $`L_1`$ are functions of the mass ratio only. As an example for $`q=0.2`$, $`r_{\mathrm{L1}}0.66a`$ and the circularization radius $`r_\mathrm{c}0.16a`$. In practice the ring forms somewhat outside $`r_\mathrm{c}`$, because there is some angular momentum redistribution in the shocks that form at the impact of the stream on the ring.

The evolution of the ring depends critically on nature and strength of the angular momentum transport processes. If sufficient ‘viscosity’ is present it spreads inward and outward to form a disk.

At the point of impact of the stream on the disk the energy dissipated is a significant fraction of the orbital kinetic energy, hence the gas heats up to a significant fraction of the virial temperature. For a typical system with $`M_1=1M_{}`$, $`M_2=0.2M_{}`$ having an orbital period of 2 hrs, the observed size of the disk (e.g. Wood et al. 1989b, Rutten et al. 1992) $`r_\mathrm{d}/a0.3`$, the orbital velocity at $`r_\mathrm{d}`$ about 900 km/s, the virial temperature at $`r_\mathrm{d}`$ is $`10^8`$K. The actual temperatures at the impact point are much lower, due to rapid cooling of the shocked gas. Nevertheless the impact gives rise to a prominent ‘hot spot’ in many systems, and an overall heating of the outermost part of the disk.

## 4 Thin disks: properties

### 4.1 Flow in a cool disk is supersonic

Ignoring viscosity, the equation of motion in the potential of a point mass is

$$\frac{𝐯}{t}+𝐯𝐯=\frac{1}{\rho }P\frac{GM}{r^2}\widehat{𝐫},$$
(8)
where $`\widehat{𝐫}`$ is a unit vector in the spherical radial direction $`r`$. To compare the order of magnitude of the terms, choose a position $`r_0`$ in the disk, and choose as typical time and velocity scales the orbital time scale $`\mathrm{\Omega }_0^1=(r_0^3/GM)^{1/2}`$ and velocity $`\mathrm{\Omega }_0r_0`$. The pressure gradient term is

$$\frac{1}{\rho }P=\frac{}{\mu }T\mathrm{ln}P.$$

In terms of the dimensionless quantities

$$\stackrel{~}{r}=r/r_0,\stackrel{~}{v}=v/(\mathrm{\Omega }_0r_0),$$

$$\stackrel{~}{t}=\mathrm{\Omega }_0t,\stackrel{~}{}=r_0,$$

the equation of motion is then

$$\frac{\stackrel{~}{𝐯}}{\stackrel{~}{t}}+\stackrel{~}{𝐯}\stackrel{~}{}\stackrel{~}{𝐯}=\frac{T}{T_{\mathrm{vir}}}\stackrel{~}{}\mathrm{ln}P\frac{1}{\stackrel{~}{r}^2}\widehat{𝐫}.$$
(9)
All terms and quantities in this equation are of order unity by the assumptions made, except the pressure gradient term which has the coefficient $`T/T_{\mathrm{vir}}`$. If cooling is important, so that $`T/T_{\mathrm{vir}}1`$, the pressure term is negligible to first approximation, and vice versa. Equivalent statements are also that the gas moves hypersonically on nearly Keplerian orbits, and that the disk is thin, as is shown next.

### 4.2 Disk thickness

The thickness of the disk is found by considering its equilibrium in the direction perpendicular to the disk plane. In an axisymmetric disk, using cylindrical coordinates ($`\varpi ,\varphi ,z`$), measure the forces at a point $`𝐫_0`$ ($`\varpi ,\varphi ,0`$) in the midplane, in a frame rotating with the Kepler rate $`\mathrm{\Omega }_0`$ at that point. The gravitational acceleration $`GM/r^2\widehat{𝐫}`$ balances the centrifugal acceleration $`\mathrm{\Omega }_0^2`$$`\varpi `$$`\varpi `$$`\varpi `$at this point, but not at some distance $`z`$ above it because gravity and centrifugal acceleration work in different directions. Expanding both accelerations near $`𝐫_0`$, one finds a residual acceleration toward the midplane of magnitude

$$g_z=\mathrm{\Omega }_0^2z.$$

Assuming again an isothermal gas, the condition for equilibrium in the $`z`$ direction under this acceleration yields a hydrostatic density distribution

$$\rho =\rho _0(\varpi )\mathrm{exp}\left(\frac{z^2}{2H^2}\right).$$

$`H`$, the scale height of the disk, is given in terms of the isothermal sound speed $`c_\mathrm{s}=(T/\mu )^{1/2}`$ by

$$H=c_\mathrm{s}/\mathrm{\Omega }_0.$$

We define $`\delta H/r`$, the aspect ratio of the disk, and find that it can be expressed in several equivalent ways:

$$\delta =\frac{H}{r}=\frac{c_\mathrm{s}}{\mathrm{\Omega }r}=M^1=\left(\frac{T}{T_{\mathrm{vir}}}\right)^{1/2},$$

where $`M`$ is the Mach number of the orbital motion.

### 4.3 Viscous spreading

The shear flow between neighboring Kepler orbits in the disk causes friction due to viscosity. The frictional torque is equivalent to exchange of angular momentum between these orbits. But since the orbits are close to Keplerian, a change in angular momentum of a ring of gas also means it must change its disctance from the central mass. If the angular momentum is increased, the ring moves to a larger radius. In a thin disk angular momentum transport (more precisely a nonzero divergence of the angular momentum flux) therefore automatically implies redistribution of mass in the disk.

A simple example (Lüst 1952, see also Lynden-Bell and Pringle 1974) is a narrow ring of gas at some distance $`r_0`$. If at $`t=0`$ this ring is released to evolve under the viscous torques, one finds that it first spreads into an asymmetric hump with a long tail to large distances. As $`t\mathrm{}`$ the hump flattens in such a way that almost all the mass of the ring is accreted onto the center, while a vanishingly small fraction of the gas carries almost all the angular momentum to infinity. As a result of this asymmetric behavior essentially all the mass of a disk can accrete, even if there is no external torque to remove the angular momentum.

### 4.4 Observations of disk viscosity

Evidence for the strength of the angular momentum transport processes in disks comes from observations of variability time scales. This evidence is not good enough to determine whether the processes really have the same effect as a viscosity, but if this is assumed, estimates can be made of the magnitude of the viscosity.

Cataclysmic Variables give the most detailed information. These are binaries with white dwarf (WD) primaries and (usually) main sequence companions (for reviews see Meyer-Hofmeister and Ritter 1993, Cordova 1995, Warner 1995). A subclass of these systems, the Dwarf Novae, show semiregular outbursts. In the currently most developed theory, these outbursts are due to an instability in the disk (Smak 1971, Meyer and Meyer-Hofmeister 1981, for recent references see, King 1995, Hameury et al. 1998). The outbursts are episodes of enhanced mass transfer of the disk onto the primary, involving a significant part of the whole disk. The decay time of the burst is thus a measure of the viscous time scale of the disk (the quantitative details depend on the model, see Cannizzo et al. 1988, Hameury et al. 1998):

$$t_{\mathrm{visc}}=r_\mathrm{d}^2/\nu ,$$

where $`r_\mathrm{d}`$ is the size of the disk. With decay times on the order of days, this yields viscosities of the order $`10^{15}`$ cm2/s, about 14 orders of magnitude above the microscopic viscosity of the gas.

Other evidence comes from the inferred time scale on which disks around protostars disappear, which is of the order of $`10^7`$ years (Strom et al, 1993).

### 4.5 $`\alpha `$-parametrization

The process responsible for such a large viscosity has not been identified with certainty yet. Many processes have been proposed, some of which demonstrably work, though often not with an efficiency as high as the observations of CV outbursts seem to indicate. Other ideas, such as certain turbulence models, do not have much predictive power and are based on ad-hoc assumptions about hydrodynamic instabilities in disks. In order to compare the viscosities in disks under different conditions, one introduces a dimensionless vsicosity $`\alpha `$:

$$\nu =\alpha \frac{c_\mathrm{s}^2}{\mathrm{\Omega }},$$
(10)
where $`c_\mathrm{s}`$ is the isothermal sound speed as before. The quantity $`\alpha `$ was introduced by Shakura and Sunyaev (1973), as a way of parametrizing our ignorance of the angular momentum transport process (their definition is based on a different formula however, and differs by a constant of order unity).

How large can the value of $`\alpha `$ be, on theoretical grounds? As a simple model, let’s assume that the shear flow between Kepler orbits is unstable to the same kind of shear instabilities found for flows in tubes, channels, near walls and in jets. These instabilities occur so ubiquitously that the fluid mechanics community considers them a natural and automatic consequence (e.g. DiPrima and Swinney 1981, p144 2nd paragraph) of a high Reynolds number:

$$\mathrm{R}e=\frac{LV}{\nu }$$

where $`L`$ and $`V`$ are characteristic length and velocity scales of the flow. If this number exceeds about 1000 (for some forms of instability much less), instability and turbulence are generally observed. It has been argued (e.g. Zel’dovich 1981) that for this reason hydrodynamic turbulence is the cause of disk viscosity. Let’s look at the consequences of this assumption. If an eddy of radial length scale $`l`$ develops due to shear instability, it will rotate at a rate given by the rate of shear, $`\sigma `$, in the flow, here

$$\sigma =r\frac{\mathrm{\Omega }}{r}\frac{3}{2}\mathrm{\Omega }.$$

The velocity amplitude of the eddy is $`V=\sigma l`$, and a field of such eddies produces a turbulent viscosity of the order (leaving out numerical factors of order unity)

$$\nu _{\mathrm{turb}}=l^2\mathrm{\Omega }.$$
(11)
In compressible flows, there is a maximum to the size of the eddy set by causality considerations. The force that allows an instability to form an overturning eddy is the pressure, which transports information about the flow at the sound speed. The eddies formed by a shear instability can therefore not move faster than $`c_\mathrm{s}`$, hence their size does not exceed $`c_\mathrm{s}/\sigma H`$. At the same time, the largest eddies formed also have the largest contribution to the turbulent viscosity. Thus we should expect that the turbulent viscosity is given by eddies with size of the order $`H`$:

$$\nu H^2\mathrm{\Omega },$$

or

$$\alpha 1.$$

Does hydrodynamical turbulence along these lines exist in disks? Unfortunately, this question is still open, but current opinion is leaning toward the view that the angular momentum transport in sufficiently ionized disks is due a a small scale magnetic field (Shakura and Sunyaev 1973). This is discussed briefly in section 8.

## 5 Thin Disks: equations

Consider a thin (= cool, nearly Keplerian, cf. section 4.2) disk, axisymmetric but not stationary. Using cylindrical coordinates ($`r,\varphi ,z`$), (note that we have changed notation from $`\varpi `$ to $`r`$ compared with section 4.2) we define the surface density $`\mathrm{\Sigma }`$ of the disk as

$$\mathrm{\Sigma }=_{\mathrm{}}^{\mathrm{}}\rho dz2H_0\rho _0,$$
(12)
where $`\rho _0`$, $`H_0`$ are the density and scaleheight at the midplane. The approximate sign is used to indicate that the coefficient in front of $`H`$ in the last expression actually depends on details of the vertical structure of the disk. Conservation of mass, in terms of $`\mathrm{\Sigma }`$ is given by

$$\frac{}{t}(r\mathrm{\Sigma })+\frac{}{r}(r\mathrm{\Sigma }v_r)=0.$$
(13)
(derived by integrating the continuity equation over $`z`$). Since the disk is axisymmetric and nearly Keplerian, the radial equation of motion reduces to

$$v_\varphi ^2=GM/r.$$
(14)
The $`\varphi `$-equation of motion is

$$\frac{v_\varphi }{t}+v_r\frac{v_\varphi }{r}+\frac{v_rv_\varphi }{r}=F_\varphi ,$$
(15)
where $`F_\varphi `$ is the azimuthal component of the viscous force. By integrating this over height $`z`$ and using (13), one gets an equation for the angular momentum balance:

$$\frac{}{t}(r\mathrm{\Sigma }\mathrm{\Omega }r^2)+\frac{}{r}(r\mathrm{\Sigma }v_r\mathrm{\Omega }r^2)=\frac{}{r}(Sr^3\frac{\mathrm{\Omega }}{r}),$$
(16)
where $`\mathrm{\Omega }=v_\varphi /r`$, and

$$S=_{\mathrm{}}^{\mathrm{}}\rho \nu dz\mathrm{\Sigma }\nu .$$
(17)
The second approximate equality in (17) holds if $`\nu `$ can be considered independent of $`z`$. The right hand side of (16) is the divergence of the viscous angular momentum flux, and is derived most easily with a physical argument, as described in, e.g. Pringle (1981) or Frank et al. (1992)111If you prefer a more formal derivation, the fastest way is to consult Landau and Lifshitz (1959) chapter 15 (hereafter LL). Noting that the divergence of the flow vanishes for a thin axisymmetric disk, the viscous stress $`\sigma `$ becomes (LL eq. 15.3)
$$\sigma _{ik}=\eta \left(\frac{v_i}{x_k}+\frac{v_k}{x_i}\right),$$
where $`\eta =\rho \nu `$. This can be written in cylindrical or spherical coordinates using LL eqs. (15.15-15.18). The viscous force is
$$F_i=\frac{\sigma _{ik}}{x_k}=\frac{1}{\eta }\frac{\eta }{x_k}\sigma _{ik}+\eta ^2v_i,$$
Writing the Laplacian in cylindrical coordinates, the viscous torque is then computed from the $`\varphi `$-component of the viscous force by multiplying by $`r`$, and is then integrated over $`z`$..

Assume now that $`\nu `$ can be taken constant with height. For an isothermal disk ($`T`$ independent of $`z`$), this is equivalent to taking the viscosity parameter $`\alpha `$ independent of $`z`$. As long as we are not sure what causes the viscosity this is a reasonable simplification. Note, however, that recent numerical simulations of magnetic turbulence suggest that the effective $`\alpha `$, and the rate of viscous dissipation per unit mass, are higher near the disk surface than near the midplane. See the discussion in section 8. While eq (16) is still valid for rotation rates $`\mathrm{\Omega }`$ deviating from Keplerian (only the integration over disk thickness must be justifiable), we now use the fact that $`\mathrm{\Omega }r^{3/2}`$. Then Eqs. (13-16) can then be combined into a single equation for $`\mathrm{\Sigma }`$:

$$r\frac{\mathrm{\Sigma }}{t}=3\frac{}{r}[r^{1/2}\frac{}{r}(\nu \mathrm{\Sigma }r^{1/2})].$$
(18)
Under the same assumptions, eq. (15) yields the mass flux $`\dot{M}`$ at any point in the disk:

$$\dot{M}=2\pi r\mathrm{\Sigma }v_r=6\pi r^{1/2}\frac{}{r}(\nu \mathrm{\Sigma }r^{1/2}).$$
(19)
Eq. (18) is the standard form of the thin disk diffusion equation. An important conclusion from this equation is: in the thin disk limit, all the physics which determines the time dependent behavior of the disk enters through one quantitity only, the viscosity $`\nu `$. This is the main attraction of the thin disk approximation.

### 5.1 Steady Thin Disks

In a steady disk ($`/t=0`$) the mass flux $`\dot{M}`$ is constant through the disk and equal to the accretion rate onto the central object. From (19) we get the surface density distribution:

$$\nu \mathrm{\Sigma }=\frac{1}{3\pi }\dot{M}\left[1\beta \left(\frac{r_i}{r}\right)^{1/2}\right],$$
(20)
where $`r_i`$ is the inner radius of the disk and $`\beta `$ is a parameter appearing through the integration constant. It is related to the flux of angular momentum $`F_J`$ through the disk:

$$F_J=\dot{M}\beta \mathrm{\Omega }_ir_i^2,$$
(21)
where $`\mathrm{\Omega }_i`$ is the Kepler rotation rate at the inner edge of the disk. If the disk accretes onto an object with a rotation rate $`\mathrm{\Omega }_{}`$ less than $`\mathrm{\Omega }_i`$, one finds (Shakura and Sunyaev, 1973, Lynden-Bell and Pringle, 1974) that $`\beta =1`$, independent of $`\mathrm{\Omega }_{}`$. Thus the angular momentum flux (torque on the accreting star) is inward (spin-up) and equal to the accretion rate times the specific angular momentum at the inner edge of the disk. For stars rotating near their maximum rate ($`\mathrm{\Omega }_{}\mathrm{\Omega }_i`$) and for accretion onto magnetospheres, which can rotate faster than the disk, the situation is different (Sunyaev and Shakura 1977, Popham and Narayan 1991, Paczyński 1991, Bisnovatyi-Kogan 1993).

Accreting magnetospheres, for example, can spin down by interaction with the disk. This case has a surface density distribution (20) with $`\beta <1`$ (see also Spruit and Taam 1993). The angular momentum flux is then outward, and the accreting star spins down. This is possible even when the interaction between the disk and the magetosphere takes place only at the inner edge of the disk. Magnetic torques due interaction with the magetosphere may exist at larger distances in the disk as well, but are not necessary for creating an outward angular momentum flux. Recent numerical simulations of disk-magnetosphere interaction (Miller and Stone 1997) give an interesting new view of how such interaction may take place, and suggests it happens very differently from what is assumed in previous ‘standard’ models.

### 5.2 Disk Temperature

In this section I assume accretion onto not-too-rapidly rotating objects, so that $`\beta =1`$. The surface temperature of the disk, which determines how much energy it loses by radiation, is governed primarily by the energy dissipation rate in the disk, which in turn is given by the accretion rate. From the first law of thermodynamics we have

$$\rho T\frac{\mathrm{d}S}{\mathrm{d}t}=\mathrm{d}iv𝐅+Q_\mathrm{v},$$
(22)
where $`S`$ the entropy per unit mass, $`𝐅`$ the heat flux (including radiation and any form of ‘turbulent’ heat transport), and $`Q_\mathrm{v}`$ the viscous dissipation rate. For changes which happen on time scales longer than the dynamical time $`\mathrm{\Omega }^1`$, the left hand side is small compared with the terms on the right hand side. Integrating over $`z`$, the divergence turns into a surface term and we get

$$2\sigma _\mathrm{r}T_\mathrm{s}^4=_{\mathrm{}}^{\mathrm{}}Q_\mathrm{v}dz,$$
(23)
where $`T_\mathrm{s}`$ is the surface temperature of the disk, $`\sigma _\mathrm{r}`$ is Stefan-Boltzmann’s radiation constant $`\sigma _\mathrm{r}=a_\mathrm{r}c/4`$, and the factor 2 comes about because the disk has 2 radiating surfaces (assumed to radiate like black bodies). Thus the energy balance is local (for such slow changes): what is generated by viscous dissipation inside the disk at any radius $`r`$ is also radiated away from the surface at that position. The viscous dissipation rate is equal to $`Q_\mathrm{v}=\sigma _{ij}v_i/x_j`$, where $`\sigma _{ij}`$ is the viscous stress tensor (see footnote in section 5), and this works out222using, e.g. LL eq. 16.3 to be

$$Q_\mathrm{v}=9/4\mathrm{\Omega }^2\nu \rho .$$
(24)
Eq. (23), using (20) then gives the surface temperature in terms of the accretion rate:

$$\sigma _\mathrm{r}T_\mathrm{s}^4=\frac{9}{8}\mathrm{\Omega }^2\nu \mathrm{\Sigma }=\frac{GM}{r^3}\frac{3\dot{M}}{8\pi }\left[1\left(\frac{r_i}{r}\right)^{1/2}\right].$$
(25)
This shows that the surface temperature of the disk, at a given distance $`r`$ from a steady accreter, depends only on the product $`M\dot{M}`$, and not on the highly uncertain value of the viscosity. For $`rr_i`$ we have

$$T_\mathrm{s}r^{3/4}.$$
(26)
These considerations only tells us about the surface temperature. The internal temperature in the disk is quite different, and depends on the mechanism transporting energy to the surface. Because it is the internal temperature that determines the disk thickness $`H`$ (and probably also the viscosity), this transport needs to be considered in some detail for realistic disk models. This involves a calculation of the vertical structure of the disk. Because of the local (in $`r`$) nature of the balance between dissipation and energy loss, such calculations can be done as a grid of models in $`r`$, without having to deal with exchange of energy between neighboring models. Schemes borrowed from stellar structure computations are used (e.g. Meyer and Meyer-Hofmeister 1982, Pringle et al. 1986, Cannizzo et al. 1988).

An approximation to the temperature in the disk can be found when a number of additional assumptions is made. As in stellar interiors, the energy transport is radiative rather than convective at high temperatures. Assuming local thermodynamic equilibrium (LTE, e.g. Rybicki and Lightman 1979), the temperature structure of a radiative atmosphere is given, in the Eddington approximation by:

$$\frac{\mathrm{d}}{\mathrm{d}\tau }\sigma _\mathrm{r}T^4=\frac{3}{4}F.$$
(27)
The boundary condition that there is no incident flux from outside the atmosphere yields the approximate condition

$$\sigma _\mathrm{r}T^4(\tau =2/3)=F,$$
(28)
where $`\tau =_z^{\mathrm{}}\kappa \rho dz`$ is the optical depth at geometrical depth $`z`$, and $`F`$ the energy flux through the atmosphere. Assuming that most of heat is generated near the midplane (which is the case if $`\nu `$ is constant with height), $`F`$ is approximately constant with height and equal to $`\sigma _\mathrm{r}T_\mathrm{s}^4`$, given by (25). Eq (27) then yields

$$\sigma _\mathrm{r}T^4=\frac{3}{4}(\tau +\frac{2}{3})F.$$
(29)
Approximating the opacity $`\kappa `$ as constant with $`z`$, the optical depth at the midplane is $`\tau =\kappa \mathrm{\Sigma }/2`$. If $`\tau 1`$, the temperature at the midplane is then:

$$T^4=\frac{27}{64}\sigma _\mathrm{r}^1\mathrm{\Omega }^2\nu \mathrm{\Sigma }^2\kappa .$$
(30)
With the equation of state (1), valid when radiation pressure is small, we find for the disk thickness, using (20):

$`{\displaystyle \frac{H}{r}}=`$ $`(/\mu )^{2/5}\left({\displaystyle \frac{3}{64\pi ^2\sigma _\mathrm{r}}}\right)^{1/10}(\kappa /\alpha )^{1/10}(GM)^{7/20}r^{1/20}(f\dot{M})^{1/5}`$ (31)
$`=`$ $`\mathrm{5\hspace{0.17em}10}^3\alpha ^{1/10}r_6^{1/20}\left(M/M_{}\right)^{7/20}(f\dot{M}_{16})^{1/5},(P_rP)`$ (32)
where $`r_6=r/(10^6`$ cm), $`\dot{M}_{16}=\dot{M}/(10^{16}`$g/s), and

$$f=1\left(r_i/r\right)^{1/2}.$$

From this we conclude that: i) the disk is thin in X-ray binaries, $`H/r<0.01`$, ii) the disk thickness is relatively insensitive to the parameters, especially $`\alpha `$, $`\kappa `$ and $`r`$. It must be stressed, however, that this depends fairly strongly on the assumption that the energy is dissipated in the disk interior. If the dissipation takes place close to the surface \[such as in some magnetic reconnection models (Haardt et al. 1994, Di Matteo et al. 1999 and references therein)\], the internal disk temperature will be much closer to the surface temperature. The midplane temperature and $`H`$ are even smaller in such disks than calculated from (32).

The viscous dissipation rate per unit area of the disk, $`W_\mathrm{v}=(9/4)\mathrm{\Omega }^2\nu \mathrm{\Sigma }`$ \[cf. eq. 25)\] can be compared with the local rate $`W_\mathrm{G}`$ at which gravitational energy is liberated in the accretion flow. Since half the gravitational energy stays in the flow as orbital motion, we have

$$W_\mathrm{G}=\frac{1}{2\pi r}\frac{GM\dot{M}}{2r^2},$$
(33)
so that

$$W_\mathrm{v}/W_\mathrm{G}=3f=3[1(r_\mathrm{i}/r)^{1/2}].$$
(34)
At large distances from the inner edge, the dissipation rate is 3 times larger than the rate of gravitational energy release. This may seem odd, but becomes understandable when it is realized that there is a significant flux of energy through the disk associated with the viscous stress333See LL section 16. Integrating the viscous energy dissipation over the whole disk, one finds

$$_{r_i}^{\mathrm{}}2\pi rW_\mathrm{v}dr=\frac{GM\dot{M}}{2r_\mathrm{i}},$$
(35)
as expected. That is, globally, but not locally, half of the gravitational energy is radiated from the disk while the other half remains in the orbital kinetic energy of the accreted material.

What happens to this remaining orbital energy depends on the nature of the accreting object. If the object is a slowly rotating black hole, the orbital energy is just swallowed by the hole. If it has a solid surface, the orbiting gas slows down until it corotates with the surface, dissipating the orbital energy into heat in a boundary layer. Unless the surface rotates close to the orbital rate (‘breakup’), the energy released in this way is of the same order as the total energy released in the accretion disk. The properties of this boundary layer are therefore crucial for accretion onto neutron stars and white dwarfs. See also section 9.1 and Inogamov and Sunyaev (1999, and elsewhere in this volume).

### 5.3 Radiation pressure dominated disks

In the inner regions of disks in XRB, the radiation pressure can dominate over the gas pressure, which results in a different expression for the disk thickness. The total pressure $`P`$ is

$$P=P_r+P_g=\frac{1}{3}aT^4+P_g.$$
(36)
Defining a ‘total sound speed’ by $`c_\mathrm{t}^2=P/\rho `$ the relation $`c_\mathrm{t}=\mathrm{\Omega }H`$ still holds. For $`P_rP_g`$ we get from (30), with (25) and $`\tau 1`$:

$$cH=\frac{3}{8\pi }\kappa f\dot{M},$$

(where the rather approximate relation $`\mathrm{\Sigma }=2H\rho _0`$ has been used). Thus,

$$\frac{H}{R}\frac{3}{8\pi }\frac{\kappa }{cR}f\dot{M}=\frac{3}{2}f\frac{\dot{M}}{\dot{M}_\mathrm{E}},$$
(37)
where $`R`$ is the stellar radius and $`\dot{M}_\mathrm{E}`$ the Eddington rate for this radius. It follows that the disk becomes thick near the star, if the accretion rate is near Eddington (though this is mitigated somewhat by the decrease of the factor $`f`$). Accretion near the Eddington limit is evidently not geometrically thin any more. In addition, other processes such as angular momentum loss by ‘photon drag’ have to be taken into account.

### 5.4 Time scales in a disk

Three locally defined time scales play a role in thin disks. The dynamical time scale $`t_\mathrm{d}`$ is the orbital time scale:

$$t_\mathrm{d}=\mathrm{\Omega }^1=(GM/r^3)^{1/2}.$$
(38)
The time scale for radial drift through the disk over a distance of order $`r`$ is the viscous time scale:

$$t_\mathrm{v}=r/(v_r)=\frac{2}{3}\frac{rf}{\nu }=\frac{2f}{3\alpha \mathrm{\Omega }}(\frac{r}{H})^2,$$
(39)
(using (19 and (20), valid for steady accretion). Finally, there are thermal time scales. If $`E_\mathrm{t}`$ is the thermal energy content (enthalpy) of the disk per unit of surface area, and $`W_\mathrm{v}=(9/4)\mathrm{\Omega }^2\nu \mathrm{\Sigma }`$ the heating rate by viscous dissipation, we can define a heating time scale:

$$t_\mathrm{h}=E_\mathrm{t}/W_\mathrm{v}.$$
(40)
In the same way, a cooling time scale is defined by the energy content and the radiative loss rate:

$$t_\mathrm{c}=E_\mathrm{t}/(2\sigma _\mathrm{r}T_\mathrm{s}^4).$$
(41)
For a thin disk, the two are equal since the viscous energy dissipation is locally balanced by radiation from the two disk surfaces. \[In thick disks (ADAFs), this balance does not hold, since the advection of heat with the accretion flow is not negligible. In ADAFs, $`t_\mathrm{c}>t_\mathrm{h}`$ (see elsewhere in this volume)\]. Thus, we can replace both time scales by a single thermal time scale $`t_\mathrm{t}`$, and find, with (24):

$$t_\mathrm{t}=\frac{1}{W_\mathrm{v}}_{\mathrm{}}^{\mathrm{}}\frac{\gamma P}{\gamma 1}dz,$$
(42)
where the enthalpy of an ideal gas of constant ratio of specific heats $`\gamma `$ has been used. Leaving out numerical factors of order unity, this yields

$$t_\mathrm{t}\frac{1}{\alpha \mathrm{\Omega }}.$$
(43)
That is, the thermal time scale of the disk is independent of most of the disk properties and of the order $`1/\alpha `$ times longer than the dynamical time scale. This independence is a consequence of the $`\alpha `$-parametrization used. If $`\alpha `$ is not a constant, but dependent on disk temperature for example, the dependence of the thermal time scale on disk properties will become apparent again.

If, as seems likely from observations, $`\alpha `$ is generally $`<1`$, we have in thin disks the ordering of time scales:

$$t_\mathrm{v}t_\mathrm{t}>t_\mathrm{d}.$$
(44)
## 6 Comparison with CV observations

The number of meaningful quantitative tests between the theory of disks and observations is somewhat limited since in the absence of a theory for $`\nu `$, it is a bit meagre on predictive power. The most detailed information perhaps comes from modeling of CV outbursts.

Two simple tests are possible (nearly) independently of $`\nu `$. These are the prediction that the disk is geometrically quite thin (eq. 32) and the prediction that the surface temperature $`T_\mathrm{s}r^{3/4}`$ in a steady disk. The latter can be tested in a subclass of the CV’s that do not show outbursts, the nova-like systems, which are believed to be approximately steady accreters. If the system is also eclipsing, eclipse mapping techniques can be used to derive the brightness distribution with $`r`$ in the disk (Horne, 1985, 1993). If this is done in a number of colors so that bolometric corrections can be made, the results (e.g. Rutten et al. 1992) show in general a fair agreement with the $`r^{3/4}`$ prediction. Two deviations occur: i) a few systems show significantly flatter distributions than predicted, and ii) most systems show a ‘hump’ near the outer edge of the disk. The latter deviation is easily explained, since we have not taken into account that the impact of the stream heats the outer edge of the disk. Though not important for the total light from the disk, it is an important local contribution near the edge.

Eclipse mapping of Dwarf Novae in quiescence gives a quite different picture. Here, the inferred surface temperature profile is often nearly flat (e.g. Wood et al. 1989a, 1992). This is understandable however since in quiescence the mass flux depends strongly on $`r`$. In the inner parts of the disk it is small, near the outer edge it is close to its average value. With eq. (25), this yields a flatter $`T_\mathrm{s}(r)`$. The lack of light from the inner disk is compensated during the outburst, when the accretion rate in the inner disk is higher than average (see Mineshige and Wood 1989 for a more detailed comparison). The effect is also seen in the 2-dimensional hydrodynamic simulations of accretion in a binary by Różyczka and Spruit (1993). These simulations show an outburst during which the accretion in the inner disk is enhanced, between two episodes in which mass accumulates in the outer disk.

## 7 Comparison with LMXB observations: irradiated disks

In low mass X-ray binaries a complication arises because of the much higher luminosity of the accreting object. Since a neutron star is roughly 1000 times smaller than a white dwarf, it produces 1000 times more luminosity for a given accretion rate.

Irradiation of the disk by the central source leads to a different surface temperature than predicted by (25). The central source (star plus inner disk) radiates the total accretion luminosity $`GM\dot{M}/R`$ (assuming sub-Eddington accretion, see section 2). If the disk is concave, it will intercept some of this luminosity. If the central source is approximated as a point source the irradiating flux on the disk surface is

$$F_{\mathrm{irr}}=ϵ\frac{GM\dot{M}}{4\pi Rr^2},$$
(45)
where $`ϵ`$ is the angle between the disk surface and the direction from a point on the disk surface to the central source:

$$ϵ=\mathrm{d}H/\mathrm{d}rH/r.$$
(46)
The disk is concave if $`ϵ`$ is positive. We have

$$\frac{F_{\mathrm{irr}}}{F}=\frac{2}{3}\frac{ϵ}{f}\frac{r}{R},$$

where $`F`$ is the flux generated internally in the disk, given by (25). On average, the angle $`ϵ`$ is of the order of the aspect ratio $`\delta =H/r`$. With $`f1`$, and our fiducial value $`\delta \mathrm{5\hspace{0.17em}10}^3`$, we find that irradiation in LMXB dominates for $`r>10^9`$cm. This is compatible with observations (for reviews see van Paradijs and McClintock 1993), which show that the optical and UV are dominated by reprocessed radiation.

When irradiation by an external source is included in the thin disk model, the surface boundary condition of the radiative transfer problem, equation (28) becomes

$$\sigma _\mathrm{r}T_\mathrm{s}^4=F+(1a)F_{\mathrm{irr}},$$
(47)
where $`a`$ is the X-ray albedo of the surface, i.e. $`1a`$ is the fraction of the incident flux that is absorbed in the optically thick layers of the disk (photons absorbed higher up only serve to heat up the corona of the disk). The surface temperature $`T_\mathrm{s}`$ increases in order to compensate for the additional incident heat flux. The magnitude of the incident flux is sensitive to the assumed disk shape $`H(r)`$, as well as on the assumed shape (plane or spherical, for example) of the central X-ray emitting region. The disk thickness depends on temperature, and thereby also on the irradiation. It turns out, however, that this dependence on the irradiating flux is small, if the disk is optically thick, and the energy transport is by radiation (Lyutyi and Sunyaev 1976). To see this, integrate (27) with the modified boundary condition (47). This yields

$$\sigma _\mathrm{r}T^4=\frac{3}{4}F(\tau +\frac{2}{3})+\frac{(1a)F_{\mathrm{irr}}}{F}.$$
(48)
The irradiation adds an additive constant to $`T^4(z)`$. At the midplane, this constant has much less effect than at the surface. For the midplane temperature and the disk thickness to be affected significantly, it is necessary that

$$F_{\mathrm{irr}}/F>\tau .$$
(49)
The reason for this weak dependence of the midplane conditions on irradiation is the same as in radiative envelopes of stars, which are also insensitive to the surface boundary condition. The situation is very different for convective disks. As in fully convective stars, the adiabatic stratification then causes the conditions at the midplane to depend much more directly on the surface temparture. The outer parts of the disks in LMXB with wide orbits may be convective, and their thickness affected by irradiation.

In the reprocessing region of the disks of LMXB, the conditions are such that $`F<<F_{\mathrm{irr}}\tau F`$, hence we must use eq. (32) for $`H`$. This yields $`ϵ=(21/20)H/r\mathrm{5\hspace{0.17em}10}^3`$, and $`T_\mathrm{s}r^{0.5}`$, and we still expect the disk to remain thin.

From the paucity of sources in which the central source is eclipsed by the companion one deduces that the companion is barely or not at all visible from the inner disk, presumably because the outer parts of the disk are much thicker than expected from the above arguments. This is consistent with the observation that the characteristic modulation of the optical light curve due to irradiation of the secondary’s surface by the X-rays is not very strong in LMXB (with the exception of Her X-1, which has a large companion). The place of the eclipsing systems is taken by the so-called ‘Accretion Disk Corona’ (ADC) systems, where shallow eclipses of a rather extended X-ray source are seen instead of the expected sharp eclipses of the inner disk (for reviews of the observations, see Lewin et al. 1995). The conclusion is that there is an extended X-ray scattering ‘corona’ above the disk. It scatters a few per cent of the X-ray luminosity.

What causes this corona and the large inferred thickness of the disk ? The thickness expected from disk theory is a rather stable small number. To ‘suspend’ matter at the inferred height of the disk forces are needed that are much larger than the pressure forces available in an optically thick disk. A thermally driven wind, produced by X-ray heating of the disk surface, has been invoked (Begelman et al. 1983, Schandl and Meyer 1994). For other explanations, see van Paradijs and McClintock (1995). Perhaps a magnetically driven wind from the disk, such as seen in protostellar objects (e.g. Königl and Ruden 1993) can explain both the shielding of the companion and the scattering. Such a model would resemble magnetically driven wind models for the broad-line region in AGN (e.g. Emmering et al., 1992, Königl and Kartje 1994). A promising possibility is that the reprocessing region at the disk edge consists of matter ‘kicked up’ at the impact of the mass transfering stream (Meyer-Hofmeister et al.1997, Armitage and Livio 1998, Spruit et al. 1998). This produces qualitatively the right dependence of X-ray absorption on orbital phase in ADC sources, and the light curves of the so-called supersoft sources.

### 7.1 Transients

Soft X-ray transients (also called X-ray Novae) are believed to be binaries similar to the other LMXB, but somehow the accretion is episodic, with very large outbursts recurring on time scales of decades (sometimes years). There are many black hole candidates among these transients (see Lewin et al. 1995 for a review). As with the Dwarf Novae, the time dependence of the accretion in transients can in principle be exploited to derive information on the disk viscosity, assuming that the outburst is caused by an instability in the disk. The closest relatives of soft transients among the White Dwarf plus main sequence star systems are probably the WZ Sge stars (van Paradijs and Verbunt 1984, Kuulkers et al. 1996), which show (in the optical) similar outbursts with similar recurrence times (cf. Warner 1987, O’Donoghue et al. 1991). Like the soft transients, they have low mass ratios ($`q<0.1`$). For a given angular momentum loss, systems with low mass ratios have low mass transfer rates, so the speculation is that the peculiar behavior of these systems is somehow connected with a low mean accretion rate.

### 7.2 Disk Instability

The most developed model for outbursts is the disk instability model of Osaki (1974), Hōshi (1979), Smak (1971, 1984), Meyer and Meyer-Hofmeister (1981), see also King (1995), Osaki (1993). In this model the instability that gives rise to cyclic accretion is due to a temperature dependence of the viscous stress. In any local process that causes an effective viscosity, the resulting $`\alpha `$\- parameter will be a function of the main dimensionless parameter of the disk, the aspect ratio $`H/r`$. If this is a sufficiently rapidly increasing function, such that $`\alpha `$ is large in hot disks and low in cool disks, an instability results by the following mechanism. Suppose we start the disk in a stationary state at the mean accretion rate. If this state is perturbed by a small temperature increase, $`\alpha `$ goes up, and by the increased viscous stress the mass flux $`\dot{M}`$ increases. By (25) this increases the disk temperature further, resulting in a runaway to a hot state. Since $`\dot{M}`$ is larger than the average, the disk empties partly, reducing the surface density and the central temperature (eq. 30). A cooling front then transforms the disk to a cool state with an accretion rate below the mean. The disk in this model switches back and forth between hot and cool states. By adjusting $`\alpha `$ in the hot and cool states, or by adjusting the functional dependence of $`\alpha `$ on $`H/r`$, outbursts are obtained that agree reasonably with the observations of soft transients (Lin and Taam 1984, Mineshige and Wheeler, 1989). A rather strong increase of $`\alpha `$ with $`H/r`$ is needed to get the observed long recurrence times.

Another possible mechanism for instability has been found in 2-D numerical simulations of accretion disks (Blaes and Hawley 1988, Różyczka and Spruit 1993). The outer edge of a disk is found, in these simulations, to become dynamically unstable to a oscillation which grows into a strong excentric perturbation (a crescent shaped density enhancement which rotates at the local orbital period). Shock waves generated by this perturbation spread mass over most of the Roche lobe, at the same time the accretion rate onto the central object is strongly enhanced. This process is different from the Smak-Osaki-Hōshi mechanism, since it requires 2 dimensions, and does not depend on the viscosity (instead, the internal dynamics in this instability generates the effective viscosity that causes a burst of accretion).

### 7.3 Other Instabilities

Instability to heating/cooling of the disk can be the due to several effects. The cooling rate of the disk, if it depends on temperature in an appropriate way, can cause a thermal instability like that in the interstellar medium. Other instabilities may result from the dependence of viscosity on conditions in the disk. For a general treatment see Piran (1978), for a shorter discussion see Treves et al., 1988.

## 8 Sources of Viscosity

The high Reynolds number of the flow in accretion disks (of the order $`10^{11}`$ in the outer parts of a CV disk) would, to most fluid dynamicists, seem an amply sufficient condition for the occurrence of hydrodynamic turbulence. A theoretical argument against such turbulence often used in astrophysics (Kippenhahn and Thomas 1981, Pringle 1981) is that in cool disks the gas moves almost on Kepler orbits, which are quite stable (except for the orbits that get close to the companion or near a black hole). This stability is related to the known stabilizing effect that rotation has on hydrodynamical turbulence (Bradshaw 1969, for a discussion see Tritton 1992). Kippenhahn and Thomas also point out that the one laboratory experiment that comes close to the situation in accretion disks, namely the rotating Couette flow, does not become unstable for parameters like in disks (for the rather limited range in Reynolds numbers available). A (not very strong) observational argument is that hydrodynamical turbulence as described above would produce an $`\alpha `$ that does not depend on the nature of the disk, so that all objects should have the same value. This is unlikely to be the case. From the modeling of CV outbursts one knows, for example, that $`\alpha `$ probably increases with temperature (more accurately, with $`H/r`$, see previous section). Also, there are indications from the inferred life times and sizes of protostellar disks (Strom et al. 1993) that $`\alpha `$ may be rather small there, $`10^3`$, whereas in outbursts of CV’s one infers values of the order $`0.11`$.

The indeterminate status of the hydrodynamic turbulence issue is an annoying problem in disk theory. Direct 3-D numerical simulation of the hydrodynamics in accretion disks is possible, and so far has not shown the expected turbulence. In fact, Balbus and Hawley (1996), and Hawley et al (1999) argue, on the basis of such simulations and a physical argument, that disks are actually quite stable against hydrodynamic turbulence, as long as the specific angular momentum increases outward. \[Such heresy would not pass a referee in a fluid mechanics journal.\] If it is true that disks are stable to hydrodynamic turbulence it will be an uphill struggle to convince the fluid mechanics community, since it can always be argued that one should go to even higher Reynolds numbers to see the expected turbulence in the simulations or experiments.

The astrophysical approach has been to circumvent the problem by finding plausible alternative mechanisms that might work just as well. Among the processes that have been proposed repeatedly as sources of viscosity is convection due to a vertical entropy gradient (e.g. Kley et al. 1993), which may have some limited effect in convective parts of disks. Another class are waves of various kinds. Their effect can be global, that is, not reducible to a local viscous term because by traveling across the disk they can communicate torques over large distances. For example, waves set up at the outer edge of the disk by tidal forces can travel inward and by dissipating there can effectively transport angular momentum outward (e.g. Narayan et al. 1987, Spruit et al. 1987). A nonlinear version of this idea are selfsimilar spiral shocks, observed in numerical simulations (Sawada et al. 1987) and studied analytically (Spruit 1987). Such shocks can produce accretion at an effective $`\alpha `$ of $`0.01`$ in hot disks, but are probably not very effective in disks as cool as those in CV’s and XRB. A second non-local mechanism is provided by a magnetically accelerated wind originating from the disk surface (Blandford 1976, Bisnovatyi-Kogan and Ruzmaikin 1976, Lovelace 1976, Blandford and Payne 1982, for reviews see Blandford 1989, Blandford and Rees 1992, for an introduction see Spruit 1996). In principle, such winds can take care of all the angular momentum loss needed to make accretion possible in the absence of a viscosity (Blandford 1976, Königl 1989). The attraction of this idea is that magnetic winds are a strong contender for explaining the strong outflows and jets seen in protostellar objects and AGN. It is not yet clear however if, even in these objects, the wind is actually the main source of angular momentum loss.

In sufficiently cool or massive disks, selfgravitating instabilities of the disk matter can produce internal friction. Paczýnski (1978) has proposed that the resulting heating would limit the instability and keep the disk in a well defined moderately unstable state. The angular momentum transport in such a disk has been modeled by several authors (e.g. Ostriker et al. 1999). Disks in XRB are too hot for selfgravity to play a role.

### 8.1 magnetic viscosity

Magnetic forces can be very effective at transporting angular momentum. If it can be shown that the shear flow in the disk produces some kind of small scale fast dynamo process, that is, some form of magnetic turbulence, an effective $`\alpha O(1)`$ expected (Shakura and Sunyaev 1973, Eardley and Lightman 1975, Pudritz 1981, Meyer and Meyer-Hofmeister 1982). Numerical simulations of initially weak magnetic fields in accretion disks have now shown that this does indeed happen in sufficiently ionized disks (Hawley et al. 1995, Brandenburg et al. 1995, Armitage 1998). These show a small scale magnetic field with azimuthal component dominating (due to stretching by differential rotation). The effective $`\alpha `$’s are of the order 0.05. The angular momentum transport is due to magnetic stresses. The fluid motions induced by the magnetic forces contribute only little to the angular momentum transport. In a perfectly conducting plasma this turbulence can develop from an arbitrarily small initial field through magnetic shear instability (also called magetorotational instability, Velikhov 1959, Chandrasekhar 1961, Balbus and Hawley 1991, 1992). The significance of this instability is that it shows that at large conductivity accretion disks must be magnetic. The actual form of the highly time dependent small scale magnetic field which develops can only be found from numerical simulations.

### 8.2 viscosity in radiatively supported disks

A disk in which the radiation pressure $`P_\mathrm{r}`$ dominates must be optically thick (otherwise the radiation would escape). The radiation pressure then adds to the total pressure is larger than it would be, for a given temperature, if only the gas pressure were effective. If the viscosity is then parametrized by (10), it turns out (Lightman and Eardley, 1974) that the disk is locally unstable. An increase in temperature increases the radiation pressure, which increases the viscous dissipation and the temperature, leading to a runaway. This has raised the question whether the radiation pressure should be included in the sound speed that enters expression (10). If it is left out, a lower viscosity results, and there is no thermal-viscous runaway. Without knowledge of the process causing the effective viscous stress, this question can not be answered. Sakimoto and Coroniti (1989) have shown, however, that if the stress is due to some form of magnetic turbulence, it most likely scales with the gas pressure alone, rather than the total pressure. Now that it seems likely, from the numerical simulations, that the stress is indeed magnetic, there is reason to believe that in the radiation pressure-dominated case the effective viscosity will scale as $`\nu \alpha P_\mathrm{g}/(\rho \mathrm{\Omega })`$ (this case has not been studied with simulations yet). Nayakshin and Rappaport (1999) show that, depending on how the viscosity scales in the intermediate regime $`P_\mathrm{g}P_{\mathrm{rad}}`$, interesting cyclic behavior can occur akin to the ‘S-curve’ instability in CV disks (section 7.2).

## 9 Beyond thin disks

Ultimately, much of the progress in developing useful models of accretion disks will depend on detailed numerical simulations in 2 or 3 dimensions. In the disks one is interested in, there is usually a large range in length scales (in LMXB disks, from less than the $`10`$km neutron star radius to the more than $`10^5`$km orbital scale). Correspondingly, there is a large range in time scales that have to be followed. This not technically possible at present and in the foreseeable future. In numerical simulations one is therefore limited to studying in an approximate way aspects that are either local or of limited dynamic range in $`r,t`$ (for examples, see Hawley 1991, Różyczka and Spruit 1993, Armitage 1998). For this reason, there is still a need for approaches that relax the strict thin disk framework somewhat without resorting to full simulations. Due to the thin disk assumptions, the pressure gradient does not contribute to the support in the radial direction and the transport of heat in the radial direction is negligible. Some of the physics of thick disks can be included in a fairly consistent way in the ‘slim disk’ approximation (Abramowicz et al., 1988). The so-called Advection Dominated Accretion Flows (ADAFs) are related to this approach (for a review see Yi 1998, for an introduction Spruit, elsewhere in this volume).

### 9.1 Boundary layers

In order to accrete onto a star rotating at the rate $`\mathrm{\Omega }_{}`$, the disk matter must dissipate an amount of energy given by

$$\frac{GM\dot{M}}{2R}\left[1\mathrm{\Omega }_{}/\mathrm{\Omega }_k(R)\right]^2.$$
(50)
The factor in brackets measures the kinetic energy of the matter at the inner edge of the disk ($`r=R`$), in the frame of the stellar surface. Due to this dissipation the disk inflates into a ‘belt’ at the equator of the star, of thickness $`H`$ and radial extent of the same order. Equating the radiation emitted from the surface of this belt to (50) one gets for the surface temperature $`T_{\mathrm{sb}}`$ of the belt, assuming optically thick conditions and a slowly rotating star ($`\mathrm{\Omega }_{}/\mathrm{\Omega }_k1`$):

$$\frac{GM\dot{M}}{8\pi R^2H}=\sigma _\mathrm{r}T_{\mathrm{sb}}^4$$
(51)
To find the temperature inside the belt and its thickness, use eq. (29). The value of the surface temperature is higher, by a factor of the order $`(R/H)^{1/4}`$, than the simplest thin disk estimate (25, ignoring the $`(r/r_i)^{1/2}`$ factor). In practice, this works out to a factor of a few. The surface of the belt is therefore not very hot. The situation is quite different if the boundary layer is not optically thick (Pringle and Savonije 1979). It then heats up to much higher temperatures. Analytical methods to obtain the boundary layer structure have been used by Regev and Hougerat (1988), numerical solutions of the slim disk type by Narayan and Popham (1993), Popham (1997), 2-D numerical simulations by Kley (1991). These considerations are primarily relevant for CV disks; in accreting neutron stars, the dominant effects of radiation pressure have to be included. More analytic progress on the structure of the boundary layer between a disk and a neutron star and the way in which it spreads over the surface of the star is reported by Inogamov and Sunyaev (1999, see also elsewhere in this volume).


text:
# On the controversy over the stochastic density functional equations

## 1 Introduction

In supercooled liquids, due to the dense packing and strong correlation of the constituent particles, the nonvibrational diffusive motions take much more time than collisions. In other words, the momentum and the energy flow much more quickly via collisions through the system than the slowly decaying number density. Recently there has been a considerable effort to describe such slow dynamics in liquids . We shall in particular discuss one of the approaches, the following stochastic equations in terms of the number density field $`\rho (x,t)`$ \[2-5\]:

$`{\displaystyle \frac{\rho (x,t)}{t}}=L[\rho (x,t)]{\displaystyle \frac{\delta H(\{\rho \})}{\delta \rho }}|_{\rho (x,t)}+\xi (x,t),`$ (1)
or its equivalent, i.e. the Fokker-Planck equation for the probability distribution functional $`P(\{\rho \},t)`$:

$$\frac{P(\{\rho \},t)}{t}=𝑑x\frac{\delta }{\delta \rho (x)}L[\rho (x)]\left[T\frac{\delta }{\delta \rho (x)}+\frac{\delta H(\{\rho \})}{\delta \rho (x)}\right]P(\{\rho \},t).$$
(2)
In equations (1) and (2), the Hamiltonian $`H`$ is of the free energy functional form as

$$H(\{\rho \})=\frac{1}{2}𝑑x𝑑y\rho (x)V(xy)\rho (y)+T𝑑x\rho (x)\mathrm{log}\rho (x),$$
(3)
$`L[\rho (x)]`$ is the kinetic coefficient written as $`L[\rho (x)]=\rho (x,t)\mathrm{\Gamma }`$ with $`\mathrm{\Gamma }`$ being the mobility of particles, and $`\xi `$ is the divergence of a random force and its correlation function is given by

$`\xi (x,t)\xi (y,t^{})=2T_xL[\rho (x)]_y\delta (xy)\delta (tt^{}).`$ (4)
For the explicit representation of this averaging, see equation (9) below.

These equations have attractive features for studying the slow dynamics: One is the physically clear incorporation of the thermally activated hopping processes via the last random term $`\xi `$ on the right hand side (rhs) of equation (1). Also of interest is the density dependence of the kinetic coefficient, $`L\rho `$, that produces the nonlinear term of dynamic origin  and implies the relevance of these equations to describing dynamical heterogeneity .

Nevertheless the stochasticity for the number density variable is still controversial. To see this, let us first mention one of some attempts \[2-4\], presented by Dean , to justify the above stochastic equations (1) and (2): Consider here a canonical system of $`N`$ particles interacting via a pairwise potential $`V(x)`$ and surrounded by a thermal white noise heat bath. The $`i`$th particle then obeys the Langevin equation,

$$\frac{dX_i(t)}{dt}=\mathrm{\Gamma }\underset{j=1}{\overset{N}{}}_iV[X_i(t)X_j(t)]+\eta _i(t),$$
(5)
where the components of the noise $`\eta _i^\alpha `$ are taken to be uncorrelated as $`\eta _i^\alpha (t)\eta _i^\beta (t)=2T\mathrm{\Gamma }\delta _{ij}\delta ^{\alpha \beta }\delta (tt^{})`$. Dean shows, using the Ito prescription for the change of variables and summing over the $`i`$, that equation (5) is transformed to the equation of the density operator, $`\widehat{\rho }(x,t)=_i\widehat{\rho }_i(x,t)=_i\delta [xX_i(t)]`$:

$`\begin{array}{cc}\hfill {\displaystyle \frac{\widehat{\rho }(x,t)}{t}}=\widehat{\rho }(x,t){\displaystyle 𝑑y\widehat{\rho }(y,t)V(xy)}+T^2\widehat{\rho }(x,t)+\widehat{\xi }(x,t)& \\ \hfill \widehat{\xi }(x,t)={\displaystyle \underset{i}{}}[\widehat{\rho }_i(x,t)\eta _i(t)].& \end{array}`$ (8)
Since one finds

$`\widehat{\xi }(x,t)\widehat{\xi }(y,t)`$ $`{\displaystyle 𝑑\eta _i\widehat{\xi }(x,t)\widehat{\xi }(y,t)\mathrm{exp}\left(𝑑t\frac{\eta _i^2}{4T\mathrm{\Gamma }}\right)}`$ (9)
$`=2T_xL[\widehat{\rho }(x)]_x\delta (tt^{}),`$ (10)
equation (1) is verified so long as the operator $`\widehat{\rho }`$ reads $`\rho `$ of c-number.

As expected, though, Marconi and Tarazona (MT)  subsequently object to the last supposition: they claim that $`\rho `$ is to be defined by averaging $`\widehat{\rho }`$ over the noise as $`\rho _{av}=\widehat{\rho }`$, where the subscript $`av`$ is appended for emphasizing the procedure. As a consequence dynamical density functional equation becomes deterministic:

$$\frac{\rho _{av}(x,t)}{t}=𝑑y\widehat{\rho }(x,t)\widehat{\rho }(y,t)V(xy)+T^2\rho _{av}(x,t),$$
(11)
whereby the Boltzmann distribution of number density is assured as time-independent solution .

To settle such controversy over the stochastic density functional equations (1) and (2), this letter aims at complementing Dean’s argument from (5) to (10) so that the above criticism by MT may become invalid. Our strategy is twofold: First, in the next section, we demonstrate that standard manipulations  transform equation (8) of the density operator to the Fokker-Planck equation (2). Moreover we verify in section 3, with the help of the conditional grand canonical partition function, the static solution $`P_0(\{\rho \})`$ of (2):

$$P_0(\{\rho \})\mathrm{exp}(\beta H),$$
(12)
where $`\beta =T^1`$. In the final section, to clarify the connection between the stochastic and the deterministic equation, we confirm using the WKB-like approximation  to the Fokker-Planck equation (2) that the noise-averaged deterministic equation (11) corresponds to that for the saddle-point path of $`P(\{\rho \},t)`$; this reveals that MT’s argument produces only the mean-field equation and not the first member of the BBGKY hierarchy  including two-point equal-time correlation function.

## 2 From equation (8) to the Fokker-Planck equation (2)

Turning our attention to the functional space, we immediately find that the density operator $`\widehat{\rho }`$ may be directly mapped to the distribution functional $`P(\{\rho \},t)`$ as

$$P(\{\rho \},t)=\underset{x}{}\delta [\widehat{\rho }(x,t)\rho (x)],$$
(13)
not via the averaging $`\rho _{av}=\widehat{\rho }`$; essentially only this definition has dissolved the MT’s critique.

Let us then exhibit below that $`P(\{\rho \},t)`$ with equation (8) satisfies the Fokker-Planck equation (2). We first differentiate (13) with respect to time, and obtain

$`{\displaystyle \frac{P(\{\rho \},t)}{t}}={\displaystyle 𝑑x\frac{\widehat{\rho }(x,t)}{t}\frac{\delta }{\delta \widehat{\rho }(x,t)}\delta [\widehat{\rho }(x,t)\rho (x)]\underset{yx}{}\delta [\widehat{\rho }(y,t)\rho (y)]}`$
$`\begin{array}{cc}\hfill ={\displaystyle }dx[\widehat{\rho }(x,t){\displaystyle }dy\widehat{\rho }(y,t)V(xy)+T^2\widehat{\rho }(x,t)+\widehat{\xi }(x,t)]& \\ \hfill \times {\displaystyle \frac{\delta }{\delta \widehat{\rho }(x,t)}}\delta [\widehat{\rho }(x,t)\rho (x)]{\displaystyle \underset{yx}{}}\delta [\widehat{\rho }(y,t)\rho (y)].& \end{array}`$ (16)
We may replace $`\widehat{\rho }`$ by $`\rho `$ using the $`\delta `$–function, and hence equation (16) reads

$`\begin{array}{cc}\hfill {\displaystyle \frac{P(\{\rho \},t)}{t}}={\displaystyle 𝑑x\frac{\delta }{\delta \rho (x)}\left[\rho (x)𝑑y\rho (y)V(xy)+T^2\rho (x)\right]P(\{\rho \},t)}& \\ \hfill {\displaystyle 𝑑x\frac{\delta }{\delta \rho (x)}\widehat{\xi }(x,t)\underset{x}{}\delta [\widehat{\rho }(x,t)\rho (x)]}.& \end{array}`$ (19)
With use of the identity for an arbitrary function $`F(\{\eta _i^\alpha \})`$, $`F(\{\eta _i^\alpha \})\eta _i^\alpha (t)=2T\mathrm{\Gamma }\delta F(\{\eta _i^\alpha \})/\delta \eta _i^\alpha (t)`$, the last bracket term on the rhs of (19) is further transformed to

$`\widehat{\xi }(x,t){\displaystyle \underset{x}{}}\delta [\widehat{\rho }(x,t)\rho (x)]=2T\mathrm{\Gamma }{\displaystyle \underset{i,\alpha }{}}{\displaystyle \frac{\widehat{\rho }_i(x,t)}{x^\alpha }}{\displaystyle \frac{\delta \widehat{\rho }(y,t)}{\delta \eta _i^\alpha }}{\displaystyle \frac{\delta }{\delta \widehat{\rho }(y,t)}}{\displaystyle \underset{y}{}}\delta [\widehat{\rho }\rho ]`$
$`=T\mathrm{\Gamma }{\displaystyle \underset{i}{}}_x_y[\widehat{\rho }_i(x,t)\widehat{\rho }_i(y,t)]{\displaystyle \frac{\delta }{\delta \widehat{\rho }(y,t)}}{\displaystyle \underset{y}{}}\delta [\widehat{\rho }\rho ],`$ (20)
where the superscript $`\alpha `$ denotes the component of $`x`$ and $`\eta _i`$, and use has been made of

$$\frac{\delta \widehat{\rho }_i(y,t)}{\delta \eta _i^\alpha }\frac{1}{2}\frac{\widehat{\rho }_i(y,t)}{y_i^\alpha }$$
(21)
that is obtained from standard mathematical manipulation of the discretized Langevin equation . Also, noting that relation $`\widehat{\rho }_i(x,t)\widehat{\rho }_i(y,t)=\delta (xy)\rho _i(x,t)`$ gives

$$_x_y[\widehat{\rho }_i(x,t)\widehat{\rho }_i(y,t)]=_x\widehat{\rho }_i(x,t)_x\delta (xy)$$
(22)
and replacing $`\widehat{\rho }`$ by $`\rho `$ as before, the bracket term finally reads

$$\widehat{\xi }(x,t)\underset{x}{}\delta [\widehat{\rho }(x,t)\rho (x)]=T\mathrm{\Gamma }\rho (x,t)\frac{\delta P(\{\rho \},t)}{\delta \rho (x,t)}.$$
(23)
Equation (19) with this is none other than the Fokker-Planck equation (2).

Thus it has been demonstrated that the stochastic equation of density operator (8) leads to the Fokker-Planck equation (2) for the (c-number) density distribution functional \[or equation (1)\].

## 3 Verification of the equilibrium distribution functional (12)

We arrive at a time-independent solution of the Fokker-Planck equation in the large time limit: $`P_0(\{\rho \})=lim_t\mathrm{}P(\{\rho \},t)`$. Therefore it is plausible to suppose that the noise averaging in calculating $`P_0`$ becomes equivalent to the configurational one in equilibrium:

$$P_0(\{\rho \})\frac{1}{N!}\underset{i=1}{\overset{N}{}}𝑑X_i\underset{x}{}\delta [\widehat{\rho }(x,t)\rho (x)]\mathrm{exp}[\beta \underset{i,j}{}V(X_iX_j)].$$
(24)
The problem is then how to derive expression (12) from the above conditional partition function.

Let us once move to the grand canonical system where we are to consider

$$P_0^\mathrm{\Xi }=\underset{N=0}{\overset{\mathrm{}}{}}P_0\lambda ^N,$$
(25)
with $`\beta \mu =\mathrm{ln}\lambda `$ being chemical potential. Introducing the auxiliary field $`\psi `$ as $`\delta [\widehat{\rho }(x)\rho ]=𝑑\psi \mathrm{exp}[i\psi (\widehat{\rho }\rho )]`$, the configurational representation of $`P_0^\mathrm{\Xi }`$ given by (24) and (25) reads

$$P_0^\mathrm{\Xi }D\psi \mathrm{exp}\left[\beta 𝑑x𝑑y\frac{1}{2}\rho (x)V(xy)\rho (y)+i\rho (x)\psi (x)e^{i\psi (x)+\mu }\right],$$
(26)
where $`D\psi `$ is formally defined as $`_xd\psi (x)`$. Since there is no contribution to $`P_0^\mathrm{\Xi }`$ of principal quadratic fluctuation of the auxiliary field $`\psi `$ around the saddle point path $`\psi _{sp}`$ as shown elsewhere , Gaussian approximation for $`\psi `$ reduces the functional integral form (26) to

$$P_0^\mathrm{\Xi }\mathrm{exp}\left[\beta H+\beta 𝑑x\rho (x)+\mu \rho (x)\right],$$
(27)
as found from substituting $`\rho =e^{i\psi _{sp}+\mu }`$ into (26).

For returning to the canonical system, we have only to perform the following contour integral,

$$P_0=\frac{1}{2\pi i}𝑑\lambda \frac{P_0^\mathrm{\Xi }}{\lambda ^{N+1}},$$
(28)
where $`\lambda `$ is now a complex variable. This relation with use of the Cauchy’s integral theorem gives back the canonical form:

$`P_0\mathrm{exp}\left[\beta H+\beta {\displaystyle 𝑑x\rho (x)}\right]{\displaystyle \frac{1}{2\pi i}}{\displaystyle 𝑑\lambda \frac{1}{\lambda ^{1+[N{\scriptscriptstyle 𝑑x\rho (x)}]}}}`$
$`=\{\begin{array}{cc}e^{\beta (HN)}\hfill & \text{if }𝑑x\rho (x)=N\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}`$ (31)
The static solution (12) has been thus verified, and supplementing Dean’s discussion has been completed.

## 4 Discussion: connection with the deterministic equation (11)

Here we would like to first confirm the saddle-point path of the Fokker-Planck (2) from exploiting the WKB-like approach made in Ref. . Setting, similarly to the WKB approximation, that

$$P(\{\rho \},t)\mathrm{exp}\left[\beta \mathrm{\Phi }(\{\rho \},t)\right],$$
(32)
we obtain the Hamilton-Jacobi like equation:

$$\frac{\mathrm{\Phi }(\{\rho \},t)}{t}=𝑑x\frac{\delta \mathrm{\Phi }(\{\rho \},t)}{\delta \rho (x)}L[\rho (x)]\left[\frac{\delta \mathrm{\Phi }(\{\rho \},t)}{\delta \rho (x)}\frac{\delta H(\{\rho \})}{\delta \rho (x)}\right].$$
(33)
A short-cut way of deriving from this the most probable (or saddle-point) path $`\{\overline{\rho }\}`$ is to expand $`\mathrm{\Phi }`$ and $`H`$ around $`\{\overline{\rho }\}`$ as

$`\begin{array}{cc}\hfill \mathrm{\Phi }(\{\rho \},t)=\mathrm{\Phi }(\{\overline{\rho }\},t)+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑x𝑑y[\rho (x)\overline{\rho }(x)]\mathrm{\Phi }\mathrm{"}(xy)[\rho (y)\overline{\rho }(y)]}+\mathrm{}& \\ \hfill J(\{\rho \})=J(\{\overline{\rho }\})+{\displaystyle \frac{\delta J(\{\rho \})}{\delta \rho }}|_{\{\rho \}=\{\overline{\rho }\}}[\rho (x)\overline{\rho }(x)]+\mathrm{},& \end{array}`$ (36)
with $`J(\{\rho \})L[\rho (x)]\delta H(\{\rho \})/\delta \rho (x)`$. Substitution of these into equation (33) yields in $`𝒪[\rho (x)\overline{\rho }(x)]`$

$$\frac{\overline{\rho }(x,t)}{t}=L[\overline{\rho }(x,t)]\frac{\delta H(\{\rho \})}{\delta \rho }|_{\{\rho \}=\{\overline{\rho }\}}.$$
(37)
Since $`\{\overline{\rho }\}`$ is to be in accord with the noise averaged density $`\{\rho _{av}\}`$, equation (37) implies that $`\widehat{\rho }(x,t)\widehat{\rho }(y,t)=\rho _{av}(x,t)\rho _{av}(y,t)`$ for the first term on the rhs of (11), i.e., no spatial correlation of noise-averaging. In other words, the above derivation reveals that the noise-averaged equation (11) is not the first member of the dynamical BBGKY hierarchy unlike the proposal by MT, but is only the mean-field equation for the saddle-point path of $`P(\{\rho \},t)`$.

We have thus validated the stochastic density functional equations, which must be a powerful tool for the understanding of supercooled fluids and glasses, via proving the irrelevance of MT’s objection to Dean’s argument in three ways: (i) demonstrating that standard manipulations enable to replace with the c-number density field $`\rho `$ the corresponding operator variable $`\widehat{\rho }`$ in the stochastic equation (8) derived by Dean, (ii) verifying the static solution (12) of the Fokker-Planck equation for the density distribution functional with the help of the conditional grand canonical partition function, and (iii) pointing out that the noise averaged path satisfying the deterministic equation (11) merely corresponds to the saddle-point one.

The next problem is how to solve these dynamical density functional equations, stochastic or deterministic. In previous works , the static density functional theory has been exploited as input, and some justifications have been also described by Kawasaki and MT . However, the present discussion does not support these; from our point of view, what to suppose for incorporating the static theory remains an open problem.

We acknowledge the financial support from the Ministry of Education, Science, Culture, and Sports of Japan.

## References