# LECTURE NOTES ON ACCRETION DISK PHYSICS

Philip J. Armitage  Stony Brook University & Center for Computational
Astrophysics, Flatiron Institute

###### Abstract

These notes introduce and review some of the physical principles underlying
the theory of astrophysical accretion, emphasizing the central roles of
angular momentum transport, angular momentum loss, and radiative cooling in
determining the structure and evolution of accretion flows. Additional topics
covered include the effective viscous theory of thin disks, classical
instabilities of disk structure, the evolution of warped or eccentric disks,
and the basic properties of waves within disks.

###### Contents

  1. I Introduction
  2. II Angular momentum transport in disk accretion
    1. II.1 Thin disk structure
    2. II.2 Microphysical viscosity in accretion disks
    3. II.3 The shearing sheet
    4. II.4 Magnetorotational instability
      1. II.4.1 The MRI dispersion relation
    5. II.5 Self-gravity
      1. II.5.1 Dispersion relation
      2. II.5.2 Toomre $Q$
      3. II.5.3 Angular momentum transport or fragmentation
    6. II.6 MHD disk winds
      1. II.6.1 Blandford-Payne disk winds
  3. III Effective viscous theory of accretion disks
    1. III.1 One dimensional time-dependent disk evolution
      1. III.1.1 1D evolution equation
      2. III.1.2 Steady solutions
      3. III.1.3 Irradiated disks
      4. III.1.4 Green’s function solution
      5. III.1.5 Self-similar solution
    2. III.2 The $\alpha$-prescription
    3. III.3 Time scales
    4. III.4 $\alpha$-model disks
    5. III.5 Self-gravitating $\alpha$-disks
    6. III.6 The values of $\alpha$
    7. III.7 The Shakura-Sunyaev solution
  4. IV Energetics of disk accretion
    1. IV.1 Eddington limit
    2. IV.2 Radiative efficiency of thin disks
      1. IV.2.1 Salpeter time
    3. IV.3 Electron-ion coupling in low density plasmas
    4. IV.4 Radiatively inefficient accretion models
      1. IV.4.1 The limiting ADAF solution
      2. IV.4.2 Variations on a theme of radiatively inefficient accretion
    5. IV.5 Hyperaccretion
  5. V Waves in disks
    1. V.1 Modified sound waves
    2. V.2 The linear wave equation
    3. V.3 Density and inertial waves
  6. VI Warped and eccentric disks
    1. VI.1 Warped disks
      1. VI.1.1 Viscous limit: heuristic derivation
      2. VI.1.2 Viscous evolution equation
      3. VI.1.3 Bardeen-Petterson effect
      4. VI.1.4 Wave-like evolution equation
      5. VI.1.5 Additional dynamical considerations
    2. VI.2 Eccentric disks
  7. VII Classical disk instabilities
    1. VII.1 Viscous instability
    2. VII.2 Thermal instability
      1. VII.2.1 The S-curve and limit cycles
    3. VII.3 Radiation warping instability
    4. VII.4 Other conjectured instabilities
  8. VIII Geometry of accretion
    1. VIII.1 Spherical (Bondi) accretion
    2. VIII.2 Supersonic accretion from an external medium

##  I Introduction

Accretion is central to astronomical systems as diverse as protoplanetary
disks, X-ray binaries, Active Galactic Nuclei, Gamma-Ray Bursts, and Tidal
Disruption Events. There are very many important differences between these
systems, but what they have in common is that they all involve accretion,
usually from a disk, on to a central gravitating object. Because the specific
angular momentum in a Keplerian potential increases outward as
$h\propto\sqrt{r}$, the physics problem posed by accretion is to understand
how gas in a rotating flow loses angular momentum and moves inward. This
problem is generally considered to be solved in principle, though important
aspects remain unclear. The astronomical problem is to figure out how the
observed manifestations of accretion, in the form of resolved images, spectra
and time variability, derive from the underlying physics. Many facets of this
problem are certainly not solved.

Indispensable references for the student of accretion are Jim Pringle’s review
of viscous disk theory [197], Steve Balbus and John Hawley’s review of
turbulence and angular momentum transport processes [14], and the textbook by
Juhan Frank, Andrew King and Derek Raine [72]. Gordon Ogilvie’s lecture
notes111 http://www.damtp.cam.ac.uk/user/gio10/accretion.html are extremely
useful, and I follow his treatment of several classical problems here. My goal
for these notes is to give an accessible introduction that follows a modern
point of view, in which the key physical processes of angular momentum
transport and radiative efficiency motivate older work on effective viscous
disk theory. Also included are a smorgasbord of results that are hard to find
outside the primary literature, extensive (but still incomplete) references,
and some editorializing on my part.

##  II Angular momentum transport in disk accretion

The angular velocity $\Omega$ of a particle in a circular orbit at distance
$r$ from a point mass $M$ in Newtonian gravity is,

$\Omega=\Omega_{\rm K}=\sqrt{\frac{GM}{r^{3}}},$ (1)

where $G$ is the gravitational constant and the subscript “K” identifies this
as being characteristic of Keplerian orbits. The specific angular momentum
(i.e. the angular momentum per unit mass) is an increasing function of orbital
radius,

$h_{\rm K}=r^{2}\Omega_{\rm K}=\sqrt{GMr}.$ (2)

This elementary property of orbits in Newtonian gravity leads immediately to
the central problem of accretion physics. Very frequently gas that is bound to
stars and compact objects has specific angular momentum that exceeds that of a
circular orbit grazing the object’s surface, and notwithstanding the
complications wrought by general relativity, pressure forces, and so forth,
often it orbits in a flattened disk-like structure with $h\approx h_{\rm K}$.
An individual parcel of gas must then lose angular momentum before it can move
to an orbit at smaller $r$ and be accreted. There are two logical
possibilities for how this can happen. One possibility is that the gas parcel
exchanges angular momentum with another parcel, so that one moves in while the
other moves out to conserve angular momentum. Invariably the viscosity of the
gas is too small to lead to appreciable angular momentum exchange on the time
scales inferred for astrophysical systems, so to realize this possibility we
require that the fluid flow exhibits an instability that leads to turbulence
and enhanced transport. The other possibility is that the disk is not an
isolated system and can lose angular momentum as a whole. This can occur if
the disk supports a magnetized outflow that exerts a torque back on the disk
to remove angular momentum.

In this Section we will first show that for an important class of disks that
are highly flattened, known as “thin disks”, the assumption that $h\approx
h_{\rm K}$ is easily justified. We justify neglecting the regular fluid
viscosity that results from particle-particle interactions, which is
negligible in most cases of interest. We then consider the linear stability of
magnetohydrodynamic (MHD) Keplerian shear flows, along with the case where the
disk is non-magnetized but self-gravitating. Finally, we will derive a simple
condition for when a thin disk, threaded by a large-scale magnetic field, can
launch an MHD wind.

###  II.1 Thin disk structure

Consider a disk of gas that orbits in the $z=0$ plane in cylindrical polar co-
ordinates $(r,\phi,z)$. We assume that the disk is accreting slowly enough
that it is in approximate hydrostatic balance, and that the mass of the disk
is negligible compared to the central mass $M$.

Figure 1: The vertical structure of a thin disk is determined by the
hydrostatic balance between the vertical component of gravity from the central
source and a pressure gradient.

To work out the vertical structure of the disk, we look at the balance of
forces in the vertical direction (shown in Figure 1). The vertical component
of gravity,

$g_{z}=g\sin(\theta)=\frac{GM}{\left(r^{2}+z^{2}\right)^{3/2}}z,$ (3)

is balanced by a pressure gradient. The vertical structure is determined by
the equation of vertical hydrostatic equilibrium,

$\frac{{\rm d}P}{{\rm d}z}=-\rho g_{z}.$ (4)

In different but reasonable physical circumstances the pressure might have
important contributions from gas, radiation, and magnetic fields. The simplest
case is when the pressure is dominated by gas pressure, and the vertical
temperature structure is isothermal. This is roughly appropriate when the disk
is optically thick and heated from outside. The equation of state is then,

$P=\rho c_{s}^{2},$ (5)

where $c_{s}$, the sound speed, is not a function of $z$. If we further assume
that $z\ll r$, the equation of hydrostatic equilibrium becomes,

$c_{s}^{2}\frac{{\rm d}\rho}{{\rm d}z}=-\Omega_{\rm K}^{2}\rho z.$ (6)

Solving this equation we find that,

$\rho(z)=\rho_{0}\exp\left[-\frac{\Omega_{\rm
K}^{2}z^{2}}{2c_{s}^{2}}\right]=\rho_{0}\left[-\frac{z^{2}}{2h^{2}}\right].$
(7)

The density at the disk mid-plane (at $z=0$) is $\rho_{0}$, and the second
equality serves to define the disk scale height,

$h\equiv\frac{c_{s}}{\Omega_{\rm K}}.$ (8)

It is often convenient to work in terms of the surface density $\Sigma$,
defined as,

$\Sigma\equiv\int_{-\infty}^{\infty}\rho(z){\rm d}z.$ (9)

For the density profile given by equation (7) the relation between the mid-
plane and surface densities is,

$\rho_{0}=\frac{1}{\sqrt{2\pi}}\frac{\Sigma}{h}.$ (10)

The above analysis is justified if $h/r=c_{s}/v_{\rm K}\ll 1$, i.e. if the
disk sound speed is much smaller than the Keplerian orbital velocity. Such
disks are described as being geometrically thin.

A geometrically thin disk, by definition, has a mid-plane sound speed that is
small compared to the Keplerian orbital velocity. In this limit, radial
pressure gradients do not affect the angular velocity profile of the disk at
leading order. The radial component of the momentum equation determines the
azimuthal velocity $v_{\phi}$ of the gas,

$\frac{v_{\phi}^{2}}{r}=\frac{GM}{r^{2}}+\frac{1}{\rho}\frac{{\rm d}P}{{\rm
d}r}.$ (11)

Approximating ${\rm d}P/{\rm d}r\sim P/r\sim\rho c_{s}^{2}/r$ we get,

$v_{\phi}^{2}=v_{\rm K}^{2}\left[1-{\cal
O}\left(\frac{h}{r}\right)^{2}\right].$ (12)

Deviations from Keplerian velocity are thus second order in $h/r$, and we can
imagine a toy disk model in which the angular velocity is Keplerian, the mid-
plane density and temperature are specified functions of radius, and the
vertical density profile is gaussian. In the right circumstances—when the disk
is geometrically thin, gas pressure dominated, and isothermal in the vertical
direction—such a model will be a decent approximation.

On occasion, one needs a true two-dimensional solution for an axisymmetric
disk in hydrostatic equilibrium. Several such solutions are known. Lin  _et
al._ [128] and Bate  _et al._ [20], for example, quote a solution for a disk
with a polytropic equation of state near the mid-plane, and an isothermal
atmosphere. This model can be useful as a background state when studying wave
propagation in disks. Fishbone and Moncrief [69] provide a solution for a
fluid torus around a black hole. Many numerical simulations of black hole
accretion are initialized with this, or similar, analytic models.

###  II.2 Microphysical viscosity in accretion disks

To accrete, gas in a Keplerian disk configuration that satisfies hydrostatic
and rotational equilibrium has to lose angular momentum. This can happen in
two ways. Gas in the disk can lose angular momentum, for example when a
magnetic field threading the disk exerts a torque on the disk surface.
Alternatively or additionally, angular momentum can be redistributed within
the disk, such that the inner part of the disk loses angular momentum and
accretes while the outer part expands to conserve angular momentum.

A rotating fluid redistributes angular momentum due to (microscopic)
viscosity, but this process is too slow to be astrophysically relevant in
essentially all disks. For an ionized gas of cosmic composition the kinematic
viscosity at temperature $T$ and density $\rho$ is [236],

$\nu=1.6\times 10^{-15}\frac{T^{5/2}}{\rho\ln\Lambda}\ {\rm cm^{2}\ s^{-1}}.$
(13)

Here, $\Lambda$ is the Coulomb logarithm for proton-proton scattering. (See
Balbus and Henri [15] for a discussion and partial derivation of this result.)
It is fiendishly hard to measure the density and temperature of most
astrophysical disks precisely, but there are many observational and
theoretical ways to get an estimate which is good enough for our purposes. For
dwarf novae (accreting white dwarfs in mass transfer binary systems), for
example, we can appeal to state-of-the-art numerical simulations by Hirose
_et al._ [91]. At a radius of $r\approx 10^{10}\ {\rm cm}$, they obtain
characteristic temperatures $T\approx 3\times 10^{4}\ {\rm K}$ and densities
$\rho\approx 10^{-7}\ {\rm g\ cm^{-3}}$. The corresponding viscosity is
$\nu\approx 2500/\ln\Lambda\ {\rm cm^{2}\ s^{-1}}$. Anticipating somewhat
results from §III.3 (or on dimensional grounds) we construct a time scale
$\tau\sim r^{2}/\nu$ assuming that the viscosity acts diffusively. With these
numbers, $\tau\sim 10^{9}\ {\rm yr}$. Observationally, however, dwarf nova
disks are seen to evolve dramatically on time scales of just days [for an
example with nice Kepler data, see 41]. Pretty clearly viscosity due to small-
scale kinetic processes in the gas is not responsible. Similar arguments apply
to protoplanetary disks and disks in AGN.

###  II.3 The shearing sheet

The low level of microphysical viscosity in accretion disks motivates study of
macroscopic instabilities that generate turbulence and angular momentum
transport. The most important of these is the magnetorotational instability
[MRI; 12], which we will discuss in §II.4. Our analysis will be based on the
shearing sheet approximation [81], which is frequently used in both analytic
and numerical studies of disks. It is a fluid cousin of Hill’s equations [90],
developed for the study of lunar motion.

To develop the shearing sheet we start with the inviscid momentum equation. In
an inertial frame this is just,

$\frac{\partial{\bf v}}{\partial t}+{\bf v}\cdot\nabla{\bf v}=-\frac{\nabla
P}{\rho}-\nabla\Phi.$ (14)

The velocity is ${\bf v}$, $P$ is the pressure, $\rho$ is the density, and
$\Phi$ is the gravitational potential. Equivalently, in terms of the
Lagrangian or material derivative,

$\frac{D{\bf v}}{Dt}=-\frac{\nabla P}{\rho}-\nabla\Phi.$ (15)

The right hand side expresses the forces acting on a fluid element, so
transforming to a frame rotating at constant angular velocity $\Omega$
requires adding the Coriolis and centrifugal forces in the same way as for
point mass dynamics,

$\frac{D{\bf v}}{Dt}=-\frac{\nabla P}{\rho}-\nabla\Phi-2{\bf\Omega}\times{\bf
v}+r\Omega^{2}{\hat{\bf r}}.$ (16)

Here $\hat{\bf r}$ is a unit vector in the radial direction. No approximation
has yet been made.

Figure 2: Geometry of the local or shearing sheet model for accretion disks.
The shearing sheet is a Cartesian $(x,y)$ representation of a small patch of
the disk, centered at $r=r_{0}$, that has an angular velocity
$\Omega=\Omega_{0}$. The disk shear is approximated by a linear function
across the sheet.

Figure 2 shows the geometry of the shearing sheet. The idea is to consider a
local patch of the disk, centered at $r_{0}$, that co-rotates with the fluid
at angular velocity $\Omega_{0}$. The dynamics is modeled in a Cartesian co-
ordinate system, within which disk quantities are replaced by constants or via
first-order Taylor expansion about $r_{0}$. Mathematically, assume that the
angular velocity can be approximated as a power-law,

$\Omega(r)=\Omega_{0}\left(\frac{r}{r_{0}}\right)^{-q}.$ (17)

$q=3/2$ corresponds to the Keplerian case (equation 1). We define a co-
rotating co-ordinate system,

$\displaystyle x$ $\displaystyle=$ $\displaystyle r-r_{0},$ (18)
$\displaystyle y$ $\displaystyle=$ $\displaystyle r\phi-\Omega_{0}t,$ (19)

and sum up the radial component of the gravitational and centrifugal forces,
keeping only the first order term in $x/r_{0}$,

$\displaystyle-\nabla\Phi+r\Omega_{0}^{2}{\hat{\bf r}}$ $\displaystyle=$
$\displaystyle-r\Omega_{0}^{2}\left(\frac{r}{r_{0}}\right)^{-q}+r\Omega_{0}^{2},$
(20) $\displaystyle\approx$ $\displaystyle 2q\Omega_{0}^{2}x.$ (21)

Adding in the vertical gravitational acceleration, with a value appropriate to
the center of the shearing sheet (using equation 6), the momentum equation in
the shearing sheet approximation is,

$\frac{D{\bf v}}{Dt}=-\frac{\nabla P}{\rho}-2{\bf\Omega}_{0}\times{\bf
v}+2q\Omega_{0}^{2}x{\hat{\bf x}}-\Omega_{0}^{2}z{\hat{\bf z}}.$ (22)

$\hat{\bf x}$ and $\hat{\bf z}$ are unit vectors in the $x$ (radial) and $z$
(vertical) directions respectively. This three-dimensional version is called a
shearing box.

Figure 3: Illustration of how radial boundary conditions are imposed in
computational studies that use the shearing box approximation [89].
Conceptually, we imagine that the simulation domain is bordered radially by
identical copies of itself, which shear azimuthally over time. Practically,
these boundary conditions are often implemented using time-dependent copying
of the fluid quantities in the active zones into one or more layers of “ghost”
zones, interpolating as necessary when the required azimuthal shift is not an
integer number of grid cells.

In addition to analytic applications, the shearing box is a useful tool for
local numerical simulations of accretion disks. The key point is that because
the angular velocity is constant across the box, the dynamical time scale
$\Omega^{-1}$ is not a function of the radial co-ordinate $x$. This makes it
possible to construct self-consistent shearing periodic radial boundary
conditions. The construction is illustrated in Figure 3 [89]. The active
computational domain is bordered radially by copies of itself, which shift
azimuthally over time at a rate that reflects the shear across the domain. For
a domain of radial extent $L_{x}$, the radial boundary conditions can be
written as [89],

$\displaystyle f(x,y,z)$ $\displaystyle=$ $\displaystyle
f(x+L_{x},y-q\Omega_{0}L_{x}t,z),$ (23) $\displaystyle v_{y}(x,y,z)$
$\displaystyle=$ $\displaystyle
v_{y}(x+L_{x},y-q\Omega_{0}L_{x}t,z)+q\Omega_{0}L_{x},$ (24)

where the first expression applies to all flow variables apart from the
azimuthal velocity $v_{y}$, which requires a boost to account for the shear
across the box.

The shearing box is a powerful tool. For analytic calculations, it captures
much of the essential small-scale dynamics of disks, while being easier to
deal with than global equations. For numerical work, on the other hand,
equation 22 is not much easier to solve than its global counterpart. The
advantage for computation studies derives mainly from the ability to use
shearing periodic radial boundary conditions, which are (relatively) non-
coercive. If, instead, one attempts to simulate a small cylindrical patch of
disk, it proves hard to devise boundary conditions which do not change the
flow dynamics. The shearing box can be readily adapted to include magnetic
fields (we will do so shortly), an energy equation, and radiation fields [92,
100]. It can also be extended to model disks whose background states are
eccentric [179] or warped [190].

The shearing box approximation can also set traps for the unwary. At the most
basic level, although we defined $x=r-r_{0}$, going to the shearing sheet
introduces a symmetry between $+x$ and $-x$, such that there is no physical
sense in which one direction is “inward” and the other “outward”. This means,
for example, that while an MRI simulation may develop sustained turbulence and
measurable stresses, the gas does not actually move radially. More
perniciously, this symmetry means that a shearing box threaded by a net
vertical magnetic field can develop a wind solution in which the gas flows to
$+x$ (say) above the mid-plane, and to $-x$ below the mid-plane. This sort of
solution is unphysical in a global context, but there is nothing to prevent it
developing within a shearing box [e.g. 9]. Another difficulty comes from
attempts to model radial gradients in disk properties, or a vertical gradient
in $\Omega$. These gradients are physically important for various disk
instabilities, but they are hard to incorporate consistently within the
shearing box formalism [153, 121].

###  II.4 Magnetorotational instability

The magnetorotational instability (MRI) is the most important process that can
generate turbulence and angular momentum transport in accretion disks [12,
14]. The MRI is a local, linear instability, that is present in accretion
disks provided that,

  * (i)

There is differential rotation, with ${\rm d}\Omega/{\rm d}r<0$. Except in
boundary layers, all disks satisfy this condition.

  * (ii)

“Weak” magnetic fields are present. The MRI can be damped by microphysical
viscosity for very weak fields [too weak to be relevant to any disks apart
perhaps from those around primordial stars; 115], and stabilized (though not
necessarily entitrely shutdown) by magnetic tension effects in disks where the
magnetic pressure is larger than the thermal pressure [192, 54].

  * (iii)

The ionization fraction is high enough to couple the magnetic field to the
fluid. This is the most important caveat. Although very low ionization
fractions suffice, protoplanetary disks can be so neutral as to call into
question the importance of the MRI [75].

These prerequisites are weak, and most astrophysical disks satisfy them.
Accordingly, we will skip over the tricky and rich question of whether non-
magnetic non-linear222Keplerian flows are linearly stable by the Rayleigh
criterion, as the specific angular momentum is an increasing function of
radius. fluid instabilities would exist in notional non-magnetic disks, and
proceed directly to the MHD case. Lyra and Umurhan [141] review potential
fluid instabilities in protoplanetary disk systems where the viability of the
MRI as a transport mechanism remains doubtful.

Figure 4: The operation of the magnetorotational instability (MRI) for a
perturbed vertical magnetic field in ideal MHD. The MRI is a fluid
instability, but the dynamics is closely related (and in some cases
mathematically equivalent) to that of two orbiting point masses connected by a
spring [13].

The physical origin of the MRI is shown in Figure 4, for the conceptually
simplest case of a disk threaded by a weak, initially vertical magnetic field.
Consider a radial perturbation to the magnetic field—exaggerated in the
cartoon—that results in linkage between an inner (blue) fluid element and an
outer (green) one. As the disk rotates, differential rotation causes the fluid
elements, initially at the same azimuth, to become azimuthally separated. This
separation is opposed by magnetic tension, which leads to a force that acts to
decrease the angular momentum of the inner fluid element and increase that of
the outer element. Because angular momentum is an increasing function of
radius in Keplerian disks, this transfer of angular momentum causes the fluid
elements to separate further radially, signalling an instability.

The existence of an instability in magnetized rotating flows was demonstrated
by Velikhov [246] and by Chandrasekhar [46], but the importance of this
instability for accretion disks was not recognized. Safronov [216] came close,
but no cigar. The MRI was rediscovered and applied to accretion disks in
breakthrough work by Balbus and Hawley [12].

####  II.4.1 The MRI dispersion relation

The dispersion relation for the MRI can be derived in several distinct ways.
Here, following Fromang [73], we work through a rather explicit calculation of
a differentially rotating disk containing a purely vertical magnetic field,
whose stability we assess in the shearing sheet approximation. The disk has a
power-law angular velocity profile, $\Omega\propto r^{-q}$, and a uniform
vertical magnetic field $B_{0}$. Radial and vertical variations in density,
and the vertical component of gravity, are ignored. The equation of state is
isothermal, $P=\rho c_{s}^{2}$, with $c_{s}$ a constant. As with any linear
stability analysis the goal is to find out whether infinitesimal perturbations
to the equilibrium are stable, or whether instead they exhibit exponentially
growth, implying instability.

The relevant equations are the continuity equation, the induction equation in
the ideal magnetohydrodynamic (MHD) limit, and the momentum equation
(including MHD forces) in the shearing box approximation,

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho{\bf v})$
$\displaystyle=$ $\displaystyle 0,$ (25) $\displaystyle\frac{\partial{\bf
v}}{\partial t}+({\bf v}\cdot\nabla){\bf v}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho}\nabla P+\frac{1}{4\pi\rho}(\nabla\times{\bf
B})\times{\bf B}-$ (26) $\displaystyle 2{\Omega}_{0}\times{\bf
v}+2q\Omega_{0}^{2}x{\hat{\bf x}},$ $\displaystyle\frac{\partial{\bf
B}}{\partial t}$ $\displaystyle=$ $\displaystyle\nabla\times({\bf v}\times{\bf
B}).$ (27)

The equilibrium state has uniform density, $\rho=\rho_{0}$, and uniform
magnetic field ${\bf B}=(0,0,B_{0})$. In this setup there are no pressure or
magnetic forces, and the initial velocity field is set by the balance of the
Coriolis and centrifugal terms,

$2\Omega_{0}\times{\bf v}=2q\Omega_{0}^{2}x{\hat{\bf x}}.$ (28)

The initial velocity field is then,

${\bf v}=\left(0,-q\Omega_{0}x,0\right).$ (29)

As expected given that we are working in the shearing sheet, there is a linear
shear around radius $r_{0}$.

To determine the stability of this system, we write the fluid variables as the
sum of their equilibrium values plus perturbations. For the MRI, it suffices
to consider a perturbation which depends only on $z$ and $t$. For the
velocity, for example, we take,

$\displaystyle v_{x}$ $\displaystyle=$ $\displaystyle v_{x}^{\prime}(z,t),$
(30) $\displaystyle v_{y}$ $\displaystyle=$
$\displaystyle-q\Omega_{0}x+v_{y}^{\prime}(z,t),$ (31) $\displaystyle v_{z}$
$\displaystyle=$ $\displaystyle v_{z}^{\prime}(z,t),$ (32)

and do similarly for the magnetic field and density. These expressions are
substituted into the continuity, momentum and induction equations, discarding
any terms that are quadratic in the primed quantities. This procedure leads to
seven equations, but to derive the MRI we need only a subset of them, namely
the $x$ and $y$ components of the momentum and induction equations. The four
linearized equations are,

$\displaystyle\frac{\partial v_{x}^{\prime}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{B_{0}}{4\pi\rho_{0}}\frac{\partial
B_{x}^{\prime}}{\partial z}+2\Omega_{0}v_{y}^{\prime},$ (33)
$\displaystyle\frac{\partial v_{y}^{\prime}}{\partial
t}-q\Omega_{0}v_{x}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{B_{0}}{4\pi\rho_{0}}\frac{\partial
B_{y}^{\prime}}{\partial z}-2\Omega_{0}v_{x}^{\prime},$ (34)
$\displaystyle\frac{\partial B_{x}^{\prime}}{\partial t}$ $\displaystyle=$
$\displaystyle B_{0}\frac{\partial v_{x}^{\prime}}{\partial z},$ (35)
$\displaystyle\frac{\partial B_{y}^{\prime}}{\partial t}$ $\displaystyle=$
$\displaystyle B_{0}\frac{\partial v_{y}^{\prime}}{\partial
z}-q\Omega_{0}B_{x}^{\prime}.$ (36)

These linearized differential equations are then converted into algebraic
equations by taking the perturbations to have the form,

$B_{x}^{\prime}=\bar{B}_{x}^{\prime}e^{i(\omega t-kz)},$ (37)

where $\omega$ is the frequency of a perturbation whose vertical wave-number
is $k$. The time derivatives give us factors of $i\omega$, while the spatial
derivatives give us $ik$. The four algebraic equations we end up with are,

$\displaystyle i\omega v_{x}^{\prime}$ $\displaystyle=$ $\displaystyle-
ik\frac{B_{0}B_{x}^{\prime}}{4\pi\rho_{0}}+2\Omega_{0}v_{y}^{\prime},$ (38)
$\displaystyle i\omega v_{y}^{\prime}$ $\displaystyle=$ $\displaystyle-
ik\frac{B_{0}B_{y}^{\prime}}{4\pi\rho_{0}}+(q-2)\Omega_{0}v_{x}^{\prime},$
(39) $\displaystyle i\omega B_{x}^{\prime}$ $\displaystyle=$ $\displaystyle-
ikB_{0}v_{x}^{\prime},$ (40) $\displaystyle i\omega B_{y}^{\prime}$
$\displaystyle=$ $\displaystyle-
ikB_{0}v_{y}^{\prime}-q\Omega_{0}B_{x}^{\prime}.$ (41)

(Here we’ve dropped bars on the variables.) The final step is to eliminate the
perturbation variables between the equations, leaving us with the dispersion
relation for the MRI. It takes the form,

$\displaystyle\omega^{4}-\omega^{2}\left[2k^{2}v_{A}^{2}+2(2-q)\Omega_{0}^{2}\right]+$
$\displaystyle k^{2}v_{A}^{2}\left[k^{2}v_{A}^{2}-2q\Omega_{0}^{2}\right]=0,$
(42)

where $v_{A}^{2}=B_{0}^{2}/(4\pi\rho_{0})$ is the Alfvén speed in the initial
state.

Figure 5: The MRI dispersion relation in ideal MHD is plotted for different
rotation laws. The lower red curve shows the $q=3/2$ case appropriate for a
Keplerian disk. The blue curves are plotted, from top downwards, for $q=-1/2$,
$q=0$, $q=1/2$, and $q=1$. Instability ($\omega^{2}/\Omega_{0}^{2}<0$) is
present at sufficiently small $k$ for $q>0$.

The MRI dispersion relation is a quadratic in $\omega^{2}$. It is plotted as a
function of $kv_{A}$, for different values of the shear parameter $q$, in
Figure 5. For the $q=3/2$ case that characterizes Keplerian disks
$\omega^{2}<0$, implying instability, for all scales larger than some minimum
(i.e. for $kv_{A}$ smaller than some maximum). Graphically, it can be seen
that the instability disappears for $q\leq 0$. Requiring that $\omega^{2}<0$
the instability criterion is,

$k^{2}v_{A}^{2}+\frac{{\rm d}\Omega^{2}}{{\rm d}\ln r}<0.$ (43)

In the limit where the vertical field $B_{0}\rightarrow 0$, the condition for
instability is just that the angular velocity decreases outward. This is
distinct from the Rayleigh stability criterion for a strictly unmagnetized
fluid, which instead requires that the specific angular momentum $r^{2}\Omega$
of the fluid decrease outward for an instability [see e.g. 201].

A number of other important quantities can be straightforwardly derived from
equation (42). Solving for ${\rm d}\omega^{2}/{\rm d}(kv_{A})=0$ on the
unstable branch of the dispersion relation yields the scale with the fastest
MRI linear growth rate. For a Keplerian disk the answer is,

$(kv_{A})_{\rm max}=\frac{\sqrt{15}}{4}\Omega_{0}.$ (44)

The growth rate at this spatial scale is,

$|\omega_{\rm max}|=\frac{3}{4}\Omega_{0}.$ (45)

This is a fast growth rate indeed! The most unstable MRI modes grow by a
factor $\exp(3\pi/2)\sim 10^{2}$ per orbit.

Returning to equation (42) we can set $\omega^{2}=0$ and find $k_{\rm crit}$,
the largest wavenumber (smallest spatial scale) which is unstable.
Specializing again to the Keplerian case we find,

$k_{\rm crit}v_{A}=\sqrt{3}\Omega_{0}.$ (46)

This result can be used to approximately quantify our earlier assertion that
the MRI can be stabilized if the vertical field is too strong. For a strong
field, the smallest unstable scale $\lambda_{\rm crit}$ may exceed the
vertical thickness of the disk, leading to a situation where the MRI is
stabilized on account of the disk being too thin to support any unstable modes
at all. For an estimate, we can set,

$\lambda_{\rm crit}=\frac{2\pi}{k_{\rm crit}}=2h,$ (47)

which leads to an estimate of the stabilizing field as,

$B_{0,{\rm max}}^{2}=\frac{12}{\pi}\rho c_{s}^{2}.$ (48)

Defining the plasma $\beta$ parameter to be the ratio of the gas pressure to
the magnetic pressure,

$\beta=\frac{8\pi P}{B_{0}^{2}},$ (49)

and making use of $h=c_{s}/\Omega_{0}$, we can re-express equation (48) as,

$\beta=\frac{2\pi^{2}}{3}.$ (50)

Having ignored vertical density gradients—which by definition are significant
on scales of the disk scale height $h$—this can only be an estimate. With that
caveat, the conclusion is that disks are unstable to the MRI provided that the
vertical magnetic field is significantly sub-thermal, in the sense of the
magnetic pressure being smaller than the thermal pressure.

Somewhat surprisingly, relaxing the severe simplifications that we have made
in the above analysis—a local shearing sheet, with a simple field geometry,
simple perturbations, and no consideration of the flow energetics—does not
materially change the answer. The review by Balbus and Hawley [14], and more
recent work by Latter  _et al._ [120], are good places to start for more
complete analyses of the MRI.

Figure 6: Rendering of the current density from a local shearing box
simulation of MRI disk turbulence. Based on simulations from Simon  _et al._
[231].

The linear analysis says nothing about whether the MRI leads to physically
significant levels of turbulence and angular momentum transport in accretion
disks. For this, we must turn to simulations, in either shearing box [89, 38]
or global [6] geometries. A local simulation of the MRI for a vertically
stratified isothermal disk [231] is shown in Figure 6. The most basic
simulation results, known now for more than twenty years, include the fact
that the MRI can act as a dynamo, sustaining turbulence and magnetic fields
against dissipation in a domain with no external currents. Angular momentum is
transported primarily due to the action of MHD (Maxwell) stresses. Stronger
levels of turbulence and transport occur if the disk is threaded by a net
vertical magnetic field. There remain open questions, of uncertain physical
significance, regarding the convergence of the MRI in local domains with
simple physics [at the least, very high resolution is needed; 212], but the
bulk of recent work has turned to models with more complete treatments of
thermal, radiative, relativistic, and plasma physics. Jiang  _et al._ [99],
for example, use radiation magnetohydrodynamic simulations to study the global
evolution of disks around supermassive black holes.

###  II.5 Self-gravity

In addition to the MRI, there are other linear instabilities of accretion disk
flow. Most of these either occur under less general circumstances than the
MRI, or have much lower growth rates, making them less important. Self-gravity
is a partial exception. The dynamics of a fluid disk that is massive and / or
cold is strongly modified by self-gravity if the Toomre $Q$ parameter,

$Q=\frac{c_{s}\kappa}{\pi G\Sigma}<1.$ (51)

Here $\kappa$ is the epicyclic frequency, $\kappa^{2}\equiv
4\Omega^{2}+2r\Omega{{\rm d}\Omega}/{{\rm d}r}$. A similar result [and the
eponymous work; 243] applies to particle disks, for which the sound speed in
the equation should be replaced by the velocity dispersion. Conditions where
$Q\sim 1$ are likely to occur at sufficiently large radii in Active Galactic
Nuclei (AGN) disks [229], and during star formation. Gravitational instability
can in turn lead to angular momentum transport, or fragmentation [for a
review, see 112].

We can gain intuition into where $Q$ comes from using a simple time scale
argument. Pressure will prevent the collapse of a patch of the disk, with
scale $\Delta r$, when the sound-crossing time $\Delta r/c_{s}$ is small
compared to the free-fall time $\sqrt{\Delta r^{3}/G\Delta r^{2}\Sigma}$.
(Ignoring factors of the order of unity.) Equating these time scales the
minimum scale of collapse can be estimated as $\Delta r\sim
c_{s}^{2}/G\Sigma$. On larger scales, collapse can be stopped by shear if the
free-fall time is long compared to the time scale for radial shear to separate
neighboring fluid elements. In a Keplerian disk the shear time scale is
$\sim\Omega^{-1}$. For a disk that is marginally unstable the minimum scale
set by pressure equals the maximum scale set by shear. This condition implies
that marginal stability occurs when $c_{s}\Omega/G\Sigma\sim 1$, as quoted
above.

####  II.5.1 Dispersion relation

Proceeding more formally [201] we consider a razor-thin circular gas disk with
uniform surface density $\Sigma_{0}$ and sound speed $c_{s}$ in the $z=0$
plane. In cylindrical polar coordinates $(r,\phi,z)$ the density of the disk
is given by,

$\rho_{0}(r,\phi,z)=\Sigma_{0}\delta(z),$ (52)

where $\delta(z)$ is a Dirac delta-function. The velocity field is,

${\bf v}_{0}(r,\phi,z)=(0,r\Omega,0).$ (53)

The angular velocity $\Omega(r)$ is not required to be Keplerian, but for
circular orbits the centrifugal force must balance gravity. If the
gravitational potential is $\Phi_{0}$,

$\Omega^{2}r=-\frac{{\rm d}\Phi_{0}}{{\rm d}r}.$ (54)

The assumptions of initially constant density and sound speed mean that
pressure gradient forces do not enter into the calculation of the equilibrium
state.

The stability analysis proceeds in a similar way as for the MRI, except that
we now ignore magnetic fields while including both the gravity of the central
object and the self-gravity of the disk itself. The equations are the
continuity and momentum equations, together with Poisson’s equation for the
gravitational field,

$\displaystyle\frac{\partial\Sigma}{\partial t}+\nabla\cdot(\Sigma{\bf v})$
$\displaystyle=$ $\displaystyle 0,$ (55) $\displaystyle\frac{\partial{\bf
v}}{\partial t}+({\bf v}\cdot\nabla){\bf v}$ $\displaystyle=$
$\displaystyle-\frac{\nabla p}{\Sigma}-\nabla\Phi,$ (56)
$\displaystyle\nabla^{2}\Phi$ $\displaystyle=$ $\displaystyle 4\pi G\rho.$
(57)

We will also make use of a two dimensional sound speed,

$c_{\rm s}^{2}\equiv\frac{{\rm d}p}{{\rm d}\Sigma},$ (58)

defined in terms of the pressure $p$ and surface density $\Sigma$ in the usual
way.

We consider infinitesimal axisymmetric perturbations to the equilibrium state,

$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\Sigma_{0}+\Sigma_{1}(r,t),$ (59) $\displaystyle p$
$\displaystyle=$ $\displaystyle p_{0}+p_{1}(r,t),$ (60) $\displaystyle\Phi$
$\displaystyle=$ $\displaystyle\Phi_{0}+\Phi_{1}(r,t),$ (61)
$\displaystyle{\bf v}$ $\displaystyle=$ $\displaystyle{\bf
v}_{0}+\left[v_{r}(r,t),\delta v_{\phi}(r,t),0\right]\\!,$ (62)

that have a spatial and temporal dependence given by (using the surface
density as an example),

$\Sigma_{1}(r,t)\propto\exp[{\rm i}(kr-\omega t)].$ (63)

Here $k$ is the spatial wavenumber of the perturbation (related to the
wavelength via $\lambda=2\pi/k$) and $\omega$ is the temporal frequency.
Making one further approximation, we assume that for the perturbations of
interest,

$kr\gg 1.$ (64)

This amounts to considering disturbances that are small compared to the radial
extent of the disk.

We now substitute the expressions for the surface density, pressure,
gravitational potential and velocity into the fluid equations, discarding any
terms we encounter that are quadratic in the perturbed quantities. For the
continuity equation this yields,

$-\\!{\rm i}\omega\Sigma_{1}+v_{r}\Sigma_{0}\left(\frac{1}{r}+{\rm
i}k\right)=0,$ (65)

which simplifies further in the local limit ($kr\gg 1$) to

$-\\!\omega\Sigma_{1}+kv_{r}\Sigma_{0}=0.$ (66)

Deriving the analogous algebraic equations from the momentum equation requires
us to express the convective operator $({\bf v}\cdot\nabla){\bf v}$ in
cylindrical coordinates. This takes the form,

$\displaystyle({\bf v}\cdot\nabla){\bf v}=\left[v_{r}\frac{\partial
v_{r}}{\partial r}+\frac{v_{\phi}}{r}\frac{\partial
v_{r}}{\partial\phi}+v_{z}\frac{\partial v_{r}}{\partial
z}-\frac{v_{\phi}^{2}}{r},\right.$ $\displaystyle\quad v_{r}\frac{\partial
v_{\phi}}{\partial r}+\frac{v_{\phi}}{r}\frac{\partial
v_{\phi}}{\partial\phi}+v_{z}\frac{\partial v_{\phi}}{\partial
z}+\frac{v_{r}v_{\phi}}{r},$ $\displaystyle\quad v_{r}\frac{\partial
v_{z}}{\partial r}+\frac{v_{\phi}}{r}\frac{\partial
v_{z}}{\partial\phi}+v_{z}\frac{\partial v_{z}}{\partial z}\Bigg{]}.$ (67)

With this in hand, the momentum equation reduces to,

$\displaystyle\displaystyle-{\rm i}\omega v_{r}-2\Omega\delta
v_{\phi}=-\frac{1}{\Sigma_{0}}\frac{{\rm d}p_{1}}{{\rm d}r}-\frac{{\rm
d}\Phi_{1}}{{\rm d}r},$ (68) $\displaystyle\displaystyle-\,{\rm i}\omega\delta
v_{\phi}+v_{r}\left[\Omega+\frac{\rm d}{{\rm
d}r}\left(r\Omega\right)\right]=0,$ (69)

where the two equations come from the radial and azimuthal components
respectively.

The next step is to relate the perturbations in pressure and gravitational
potential expressed on the right-hand-side of equation (68) to perturbations
in the surface density. For the pressure term this is straightforward.
Equation (58) implies that,

$\frac{1}{\Sigma_{0}}\frac{{\rm d}p_{1}}{{\rm d}r}=\frac{1}{\Sigma_{0}}c_{\rm
s}^{2}{\rm i}k\Sigma_{1}.$ (70)

Dealing with the potential perturbations requires more work. Starting from the
linearized Poisson equation,

$\nabla^{2}\Phi_{1}=4\pi G\Sigma_{1}\delta(z),$ (71)

we write out the Laplacian explicitly and simplify making use of the fact that
for short wavelength perturbations $kr\gg 1$. This yields a relation between
the density and potential fluctuations,

$\frac{{\rm d}^{2}\Phi_{1}}{{\rm d}z^{2}}=k^{2}\Phi_{1}+4\pi
G\Sigma_{1}\delta(z).$ (72)

For $z\neq 0$ the only solution to this equation that remains finite for large
$|z|$ has the form,

$\Phi_{1}=C\exp[-|kz|],$ (73)

where $C$ remains to be determined. To do so we integrate the Poisson equation
vertically between $z=-\epsilon$ and $z=+\epsilon$,

$\int_{-\epsilon}^{+\epsilon}\nabla^{2}\Phi_{1}{\rm
d}z=\int_{-\epsilon}^{+\epsilon}4\pi G\Sigma_{1}\delta(z){\rm d}z.$ (74)

Noting that both $\partial^{2}\Phi_{1}/\partial x^{2}$ and
$\partial^{2}\Phi_{1}/\partial y^{2}$ are continuous at $z=0$, whereas
$\partial^{2}\Phi_{1}/\partial z^{2}$ is not, we obtain,

$\left.\frac{{\rm d}\Phi_{1}}{{\rm d}z}\right|_{-\epsilon}^{+\epsilon}=4\pi
G\Sigma_{1}.$ (75)

Taking the limit $\epsilon\rightarrow 0$ we find that $C=-2\pi
G\Sigma_{1}/|k|$, and hence that the general relation between potential and
surface density fluctuations on the $z=0$ plane is,

$\Phi_{1}=-\frac{2\pi G\Sigma_{1}}{|k|}.$ (76)

Taking the radial derivative,

$\frac{{\rm d}\Phi_{1}}{{\rm d}r}=-\frac{2\pi{\rm i}kG\Sigma_{1}}{|k|},$ (77)

which allows us to eliminate the potential from the right-hand-side of
equation (68) in favor of the surface density. The result is,

$-\\!{\rm i}\omega v_{r}-2\Omega\delta v_{\phi}=-\frac{1}{\Sigma_{0}}c_{\rm
s}^{2}{\rm i}k\Sigma_{1}+\frac{2\pi{\rm i}kG\Sigma_{1}}{|k|}.$ (78)

Finally, we are ready to derive the functional relationship between $\omega$
and $k$ (the dispersion relation). Eliminating $v_{r}$ and $\delta v_{\phi}$
between equations (66), (69), and (78) we find that,

$\omega^{2}=\kappa^{2}+c_{\rm s}^{2}k^{2}-2\pi G\Sigma_{0}|k|,$ (79)

where the epicyclic frequency $\kappa$ is defined as,

$\kappa^{2}\equiv 4\Omega^{2}+2r\Omega\frac{{\rm d}\Omega}{{\rm d}r}.$ (80)

In a Keplerian potential $\kappa^{2}=\Omega^{2}$.

####  II.5.2 Toomre $Q$

Figure 7: The dispersion relation for a self-gravitating accretion disk in the
Keplerian limit is plotted for different values of the Toomre $Q$ parameter,
$Q=0.75$, $Q=1$ and $Q=1.5$. Instability is present at some scales for $Q<1$.
The blue dot-dashed lines show the contributions from pressure (dominant at
large $|k|$), rotation (dominant as $|k|\rightarrow 0$, and self-gravity
(destabilizing).

Let us consider, for simplicity, a Keplerian disk. Noting that
$h=c_{s}/\Omega$, we can write the dispersion relation in the form,

$\frac{\omega^{2}}{\Omega^{2}}=1-\frac{2h|k|}{Q}+h^{2}k^{2},$ (81)

where $Q=c_{s}\Omega/(\pi G\Sigma)$ is the dimensionless Toomre $Q$ parameter
that we previously deduced using a time scale argument. It is plotted in
Figure 7 for different values of $Q$. Each of the terms on the right-hand-side
has a simple physical interpretation. The constant term describes the effect
of rotation, which stabilizes all scales but which is dominant at small $h|k|$
(i.e. at large spatial scales). The quadratic term describes the effect of
pressure, which has the opposite tendency and preferentially stabilizes small
spatial scales. The linear term, which describes self-gravity, is negative and
thus destabilizing. The strength of self-gravity is fully parameterized by the
value $Q$, with small enough $Q$ leading to $\omega^{2}<0$ and unstable modes.
It is easy to verify that the criterion for instability corresponds to $Q<1$.

####  II.5.3 Angular momentum transport or fragmentation

The onset of self-gravity in an accretion disk can lead to angular momentum
transport or to fragmentation. Loosely speaking, fragmentation is the outcome
for disks that either cool too quickly [78, 209], or that are fed mass from
infall at too high a rate [113]. Kratter and Lodato [112] discuss the
quantitative thresholds, determined from simulations, for these outcomes.

###  II.6 MHD disk winds

Several physical processes can lead to mass loss in a disk wind from
geometrically thin accretion disks, including,

  * •

A vertical gradient of thermal pressure. A sufficiently strong thermal
pressure can develop, even in a thin disk, if high energy radiation heats the
disk’s surfaces to a temperature where $c_{s}\sim v_{\rm K}$. The hot gas then
has positive total energy and can escape to infinity. Thermal or
photoevaporative winds may be important in AGN [23] and in protoplanetary
disks [16], where they contribute to the dispersal of the disk at the end of
the protoplanetary disk phase [5, 66].

  * •

Radiation pressure acting on spectral lines. The physics of line-driven disk
winds is an extension of the accepted theory for mass loss from massive stars
[43]. For accretion disks, line-driving is efficient for AGN disks that are
strong emitters of ultraviolet radiation [202]. (Radiation absorbed by the
continuum can also be important, but typically only under conditions where the
disk is geometrically thick.)

  * •

MHD acceleration, which in general involves contributions from centrifugal
acceleration [30] along poloidal magnetic field lines (i.e. in the
$(r,\theta)$ plane of spherical polar co-ordinates), and from a gradient of
toroidal magnetic pressure [139].

In addition to these processes, which in principle can drive a wind from a
broad range of disk radii, there are others whose applicability is limited to
the innermost disk region. Disk magnetic fields interacting with a spinning
black hole can extract spin energy via the Blandford-Znajek effect [31], while
for accreting objects with a material surface the interaction between the
inner disk and a magnetosphere can generate outflow [230]. Inner disk
processes are implicated in the formation of well-collimated jets in accreting
systems.

Of these mass loss processes, magneto-centrifugal disk winds are particularly
interesting because they can extract angular momentum as well as mass from the
underlying disk flow. The basic idea is that gas, accelerating outwards along
a poloidal field line, rotates with the angular velocity of the gas at the
field line’s foot point in the disk. The departing gas thus gains angular
momentum, which is removed from the disk by a magnetic torque that can be
considered to act on the disk’s surface. The angular momentum loss leads to
inflow of gas in the disk, independent of internal angular momentum transport
processes such as the MRI or self-gravity.

The clearest description of the physics of MHD disk winds that I’m aware of
can be found in Spruit [237]. Here, we content ourselves with a derivation of
one of the basic properties of a Blandford and Payne [30] disk wind—the
existence of a critical inclination angle for a poloidal disk-threading field
to launch a cold wind.

####  II.6.1 Blandford-Payne disk winds

Figure 8: Illustration [adapted from the pedagogical review by 237] of a
magnetized disk wind driven by centrifugal force. Left panel: the region above
the disk surface is force-free. If magnetic field lines threading the disk are
inclined by a large enough angle with respect to the vertical, centrifugal
force can accelerate gas along the field lines even in the absence of a
pressure gradient. The dynamics is equivalent to a mechanical system of a bead
(the gas) on a rigid wire (the magnetic field line) that rotates with the
angular velocity of the disk at the foot point. Right panel: the geometry for
calculating the critical inclination angle for cold magneto-centrifugal wind
launching.

The geometry of a Blandford-Payne wind [30] wind is illustrated in Figure 8.
We envisage a Keplerian disk threaded by a large-scale poloidal magnetic
field, in the limit of ideal MHD. Within the disk the energy density in the
magnetic field, $B^{2}/8\pi$, is usually smaller than $\rho c_{s}^{2}$, the
thermal energy. Due to flux conservation, however, the energy in the vertical
field component, $B_{z}^{2}/8\pi$, is roughly constant with height for $z<r$,
while the gas pressure decreases rapidly (for a thin isothermal disk, as a
gaussian with a scale height $h\ll r$). This leads to a region above the disk
surface where magnetic forces dominate. The magnetic force per unit volume can
be written as the sum of a magnetic pressure gradient and a force due to
magnetic tension,

$\frac{{\bf J}\times{\bf
B}}{c}=-\nabla\left(\frac{B^{2}}{8\pi}\right)+\frac{{\bf B}\cdot\nabla{\bf
B}}{4\pi},$ (82)

where the current,

${\bf J}=\frac{c}{4\pi}\nabla\times{\bf B}.$ (83)

In the disk wind region where magnetic forces dominate, the requirement that
they exert a finite acceleration on the low density gas can only be satisfied
if the force approximately vanishes, i.e. that,

${\bf J}\times{\bf B}\approx 0.$ (84)

The structure of the magnetic field in the magnetically dominated region is
then described as being “force-free”, and in the disk wind case (where $B$
changes slowly with $z$) the field lines must be approximately straight to
ensure that the magnetic tension term is small. If the field lines support a
wind, the force-free structure persists up to where the kinetic energy density
in the wind, $\rho v^{2}$, first exceeds the magnetic energy density. This
criterion defines the Alfvén surface. Beyond the Alfvén surface, the inertia
of the gas in the wind is sufficient to bend the field lines, which wrap up
into a spiral structure as the disk below them rotates.

Magneto-centrifugal driving can launch a wind from the surface of a cold gas
disk if the magnetic field lines are sufficiently inclined to the disk normal.
The critical inclination angle in ideal MHD can be derived via an exact
mechanical analogy. To proceed, we note that in the force-free region the
magnetic field lines are (i) basically straight lines, and (ii) enforce rigid
rotation out to the Alfvén surface at an angular velocity equal to that of the
disk at the field line’s footpoint. The geometry is shown in Figure 8. We
consider a field line that intersects the disk at radius $r_{0}$, where the
angular velocity is $\Omega_{0}=\sqrt{GM/r_{0}^{3}}$, and that makes an angle
$\theta$ to the disk normal. We define the spherical polar radius $r$, the
cylindrical polar radius $\varpi$, and measure the distance along the field
line from its intersection with the disk at $z=0$ as $s$. In the frame co-
rotating with $\Omega_{0}$ there are no magnetic forces along the field line
to affect the acceleration of a wind; the sole role of the magnetic field is
to constrain the gas to move along a straight line at constant angular
velocity. Following this line of argument, the condition for the acceleration
of the wind can be described in terms of an effective potential,

$\Phi_{\rm eff}(s)=-\frac{GM}{r(s)}-\frac{1}{2}\Omega_{0}^{2}\varpi^{2}(s),$
(85)

that is the sum of the gravitational potential and the centrifugal potential
in the rotating frame.

Figure 9: The effective potential for a cold magnetocentrifugally driven disk
wind. Acceleration occurs for field lines that are inclined by $30^{\circ}$ or
more away from the vertical.

Written out explicitly the effective potential is,

$\displaystyle\Phi_{\rm
eff}(s)=-\frac{GM}{(s^{2}+2sr_{0}\sin\theta+r_{0}^{2})^{1/2}}$
$\displaystyle-\frac{1}{2}\Omega_{0}^{2}\left(r_{0}+s\sin\theta\right)^{2}.$
(86)

This function is plotted in Figure 9 for various values of the angle $\theta$.
If we consider first a vertical field line ($\theta=0$) the effective
potential is a monotonically increasing function of distance $s$. For modest
values of $\theta$ there is a potential barrier defined by a maximum at some
$s=s_{\rm max}$, while for large enough $\theta$ the potential decreases
monotonically from $s=0$. In this last case purely magneto-centrifugal forces
suffice to accelerate a wind off the disk surface, even in the absence of any
thermal effects. We compute the critical inclination angle $\theta_{\rm
crit}$, defined as the minimum angle that allows magneto-centrifugal wind
driving, via the condition,

$\left.\frac{\partial^{2}\Phi_{\rm eff}}{\partial s^{2}}\right|_{s=0}=0.$ (87)

Evaluating this condition, we find,

$\displaystyle 1-4\sin^{2}\theta_{\rm crit}$ $\displaystyle=$ $\displaystyle
0$ $\displaystyle\Rightarrow\theta_{\rm crit}$ $\displaystyle=$ $\displaystyle
30^{\circ},$ (88)

as the minimum inclination angle from the vertical needed for unimpeded wind
launching in ideal MHD.

The rigid rotation of the field lines interior to the Alfvén surface means
that gas being accelerated along them increases its specific angular momentum.
The magnetic field, in turn, applies a torque to the disk that removes a
corresponding amount of angular momentum. If a field line, anchored to the
disk at radius $r_{0}$, crosses the Alfvén surface at (cylindrical) radius
$r_{A}$, it follows that the angular momentum flux is,

$\dot{L}_{w}=\dot{M}_{w}\Omega_{0}r_{A}^{2},$ (89)

where $\dot{M}_{w}$ is the mass loss rate in the wind. Removing angular
momentum at this rate from the disk results in a local accretion rate
$\dot{M}=2\dot{L}_{w}/\Omega_{0}r_{0}^{2}$. The ratio of the disk accretion
rate to the wind loss rate is,

$\frac{\dot{M}}{\dot{M}_{w}}=2\left(\frac{r_{A}}{r_{0}}\right)^{2}.$ (90)

If $r_{A}$ substantially exceeds $r_{0}$ (by a factor of a few, which is
reasonable for detailed disk wind solutions) a relatively weak wind can carry
away enough angular momentum to support a much larger accretion rate.

##  III Effective viscous theory of accretion disks

The MRI, self-gravity, and disk winds are physical processes that lead to
angular momentum transport and evolution of accretion disks. In many cases, it
is not possible to simulate these processes with sufficient fidelity, or over
enough radial range and for long enough, to compare against observations. It
can therefore be necessary to turn to viscous disk theory, historically
developed much earlier, which can be used to model long-term disk evolution
[225, 140].

The basic assumption of an effective viscous disk theory is that angular
momentum transport within the disk can be represented approximately as a fluid
viscosity, using some parameterization that does not explicitly involve B. The
continuity and momentum equations are then,

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho{\bf
v}\right)$ $\displaystyle=$ $\displaystyle 0,$ (91)
$\displaystyle\frac{\partial{\bf v}}{\partial t}+{\bf v}\cdot\nabla{\bf v}$
$\displaystyle=$ $\displaystyle-\frac{\nabla
P}{\rho}-\nabla\Phi+\frac{1}{\rho}\nabla\cdot{\bf T}.$ (92)

Here $\rho$ is density, ${\bf v}$ is velocity, $P$ is pressure, $\Phi$ is
gravitational potential, and ${\bf T}$ is stress (represented by a tensor). We
further assume that we are dealing with a geometrically thin disk for which
the pressure gradient term is small. We can then make progress using
simplified versions of the equations for mass and momentum conservation.

###  III.1 One dimensional time-dependent disk evolution

The fluid equations (92) apply generally. We first specialize to the case of a
geometrically thin, circular, and planar disk, and derive a one-dimensional
(in radius $r$) evolution equation. We then look for time-independent and
time-dependent solutions.

####  III.1.1 1D evolution equation

In cylindrical polar co-ordinates $(r,\phi,z)$ the continuity equation is,

$\frac{\partial\rho}{\partial t}+\frac{1}{r}\frac{\partial}{\partial
r}\left(r\rho v_{r}\right)+\frac{1}{r}\frac{\partial}{\partial\phi}\left(\rho
v_{\phi}\right)+\frac{\partial}{\partial z}\left(\rho v_{z}\right)=0.$ (93)

Integrating over $\phi$ $[0,2\pi]$ and over $z$ $[-\infty,\infty]$ the first
term becomes,

$\frac{\partial}{\partial t}\int_{0}^{2\pi}\int_{-\infty}^{\infty}\rho
dzd\phi=\frac{\partial}{\partial t}\left(2\pi\Sigma\right),$ (94)

where $\Sigma(r,t)$ is the azimuthally averaged surface density. We are
working toward an evolution equation for this quantity. Integrating the second
term,

$\frac{1}{r}\frac{\partial}{\partial
r}\int_{0}^{2\pi}\int_{-\infty}^{\infty}r\rho
v_{r}dzd\phi=\frac{1}{r}\frac{\partial{\cal{F}}}{\partial r},$ (95)

where $\cal{F}$, the radial mass flux,

${\cal{F}}=\int_{0}^{2\pi}\int_{-\infty}^{\infty}r\rho v_{r}dzd\phi=2\pi
r\Sigma\bar{v}_{r},$ (96)

can be written in terms of the surface density and the density-weighted radial
velocity $\bar{v}_{r}$. On integration, the third and fourth terms vanish (the
latter assuming that there is no mass loss from the surfaces of the disk)
leaving,

$\frac{\partial\Sigma}{\partial t}+\frac{1}{r}\frac{\partial}{\partial
r}\left(r\Sigma\bar{v}_{r}\right)=0,$ (97)

as the one-dimensional version of the continuity equation.

Dealing with the momentum equation in the same way, we can simplify the
algebra by assuming at the outset that the disk is axisymmetric such that the
orbital velocity,

$v_{\phi}=r\Omega,$ (98)

depends only on the angular velocity $\Omega=\Omega(r)$ of circular orbits in
a fixed gravitational potential $\Phi$. The only non-zero terms in the $\phi$
component of the momentum equation then come from ${\bf v}\cdot\nabla{\bf v}$
and $\nabla{\bf\cdot T}$. Looking up the forms for these in cylindrical polar
co-ordinates the surviving terms are,

$\displaystyle\left.{\bf v}\cdot{\nabla{\bf v}}\right|_{\phi}$
$\displaystyle=$ $\displaystyle v_{r}\frac{\rm dv_{\phi}}{\rm
dr}+\frac{v_{\phi}v_{r}}{r},$ (99)
$\displaystyle\left.\frac{1}{\rho}\nabla\cdot{\bf T}\right|_{\phi}$
$\displaystyle=$ $\displaystyle\frac{1}{r^{2}\rho}\frac{\partial}{\partial
r}\left(r^{2}T_{r\phi}\right)+\frac{1}{\rho}\frac{\partial T_{\phi
z}}{\partial z}.$ (100)

Multiplying these expressions through by $r\rho$ the azimuthal component of
the momentum equation takes the form,

$\rho v_{r}\frac{{\rm d}h}{{\rm d}r}=\frac{1}{r}\frac{\partial}{\partial
r}\left(r^{2}T_{r\phi}\right)+r\frac{\partial T_{\phi z}}{\partial z},$ (101)

where the specific angular momentum $h$ is defined as,

$h\equiv r^{2}\Omega.$ (102)

Finally we multiply equation (101) through again by $r$, and integrate over
$\phi$ and $z$. If $T_{\phi z}$ vanishes as $z\rightarrow\pm\infty$ the result
is,

${\cal{F}}\frac{{\rm d}h}{{\rm d}r}=-\frac{\partial{\cal G}}{\partial r},$
(103)

where ${\cal F}$ is given by equation (96) and,

${\cal{G}}=-\int_{0}^{2\pi}\int_{-\infty}^{\infty}r^{2}T_{r\phi}dzd\phi,$
(104)

is the viscous torque.

Getting to this point from the conservation laws expressed in equation (92)
requires only some rather transparent assumptions: axisymmetry, a time-
independent potential, and no mass or angular momentum loss from the disk
surfaces. One could stop there, and consider $T_{r\phi}$ to be the key
quantity whose dependence on disk conditions needs to be determined.
Conventionally, however, we instead write the torque in terms of an effective
fluid viscosity that follows a Navier-Stokes form. For a fluid with viscosity
$\mu$ and bulk viscosity $\mu_{\rm b}$ we have,

${\bf T}=\mu\left[\nabla{\bf v}+\left(\nabla{\bf
v}\right)^{T}\right]+\left(\mu_{\rm
b}-\frac{2}{3}\mu\right)\left(\nabla\cdot{\bf v}\right){\bf I}.$ (105)

To leading order the divergence of a thin disk velocity field vanishes, so we
don’t have to worry about bulk viscosity at all. The $r\phi$ component of the
stress is,

$T_{r\phi}=\mu r\frac{{\rm d}\Omega}{{\rm d}r}.$ (106)

Defining the kinematic viscosity (later just “the viscosity”) $\nu$ as,

$\nu=\frac{1}{2\pi\Sigma}\int_{0}^{2\pi}\int_{-\infty}^{\infty}\mu dzd\phi,$
(107)

the viscous torque has a fairly intuitive form that is the product of the
circumference, the viscous force per unit length, and the lever arm,

${\cal G}=-2\pi r\cdot\nu\Sigma r\frac{{\rm d}\Omega}{{\rm d}r}\cdot r.$ (108)

Equation (103) is then,

$\Sigma\bar{v}_{r}\frac{{\rm d}h}{{\rm d}r}=\frac{1}{r}\frac{\rm d}{{\rm
d}r}\left(\nu\Sigma r^{3}\frac{{\rm d}\Omega}{{\rm d}r}\right).$ (109)

Given the aforementioned assumptions, this equation expresses angular momentum
conservation for a viscous fluid in a disk geometry.

Eliminating the radial velocity $\bar{v}_{r}$ between equation (97) and
equation (109) we obtain,

$\frac{\partial\Sigma}{\partial t}=-\frac{1}{r}\frac{\partial}{\partial
r}\left[\left(\frac{{\rm d}h}{{\rm d}r}\right)^{-1}\frac{\partial}{\partial
r}\left(\nu\Sigma r^{3}\frac{{\rm d}\Omega}{{\rm d}r}\right)\right].$ (110)

This form is valid for an arbitrary (fixed) profile of angular velocity and
angular momentum in the disk. Very often we are interested in the case of a
disk that orbits a Newtonian point mass $M$. In that limit,

$\displaystyle\Omega$ $\displaystyle=$ $\displaystyle\Omega_{\rm
K}=\sqrt{GM/r^{3}},$ (111) $\displaystyle h$ $\displaystyle=$
$\displaystyle\sqrt{GMr}.$ (112)

The radial velocity is given by equation (109) as,

$\bar{v}_{r}=-\frac{3}{\Sigma r^{1/2}}\frac{\rm d}{{\rm d}r}\left(\nu\Sigma
r^{1/2}\right),$ (113)

and the evolution equation has the form,

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

The surface density thus evolves according to a diffusive partial differential
equation, whose precise character depends upon the nature of the viscosity.
The equation is linear if $\nu\neq f(\Sigma)$, though there is no general
reason for this to be the case.

####  III.1.2 Steady solutions

Steady solutions to equation (114) are easily derived. Setting
$\partial\Sigma/\partial t=0$ and integrating we find that,

$\nu\Sigma=c_{1}+c_{2}r^{-1/2}.$ (115)

Determining the constants of integration takes a little more work, and the
right answer depends on the physics of the disk being modeled. It’s easiest to
start from equation (109). Noting that the accretion rate $\dot{M}$ is,

$\dot{M}=-2\pi r\Sigma\bar{v}_{r},$ (116)

and that $\dot{M}$ must be constant for a steady solution, we integrate
equation (109). The result is,

$-\frac{\dot{M}}{2\pi}h=\nu\Sigma r^{3}\frac{{\rm d}\Omega}{{\rm d}r}+{\rm
const}.$ (117)

For Keplerian (point mass) forms for $h$ and $\Omega$, the constant term is
negligible at large $r$. We have, as $r\rightarrow\infty$,

$\nu\Sigma\simeq\frac{\dot{M}}{3\pi}.$ (118)

The surface density of the disk is inversely proportional to the viscosity.

To get at the second constant of integration, we note that the constant
appearing in equation (117) is proportional to an angular momentum flux
$\dot{M}h$. To obtain the standard form of the steady disk solution we assume
that at some radius $r=\tilde{r}$ the first term on the right-hand-side of
equation (117), which is proportional to the viscous torque ${\cal G}$,
vanishes. If $h$ and $\Omega$ are given by Keplerian expressions, we then
find,

$\nu\Sigma=\frac{\dot{M}}{3\pi}\left[1-\sqrt{\frac{\tilde{r}}{r}}\right].$
(119)

This is the solution for a steady-state disk subject to a zero-torque boundary
condition at $r=\tilde{r}$. Classically, this boundary condition can be
physically justified for a disk around a slowly rotating, non-magnetized star,
with $\tilde{r}\simeq r_{*}$, the stellar radius [196], and for a disk around
a black hole, in which case $\tilde{r}$ can be identified with the innermost
stable circular orbit [17, 225, 184]. However, in neither situation is the
justification watertight [for the black hole case see, e.g.; 76, 4].

The heating rate per unit volume in the disk is given by,

$q_{+}=T_{r\phi}r\frac{{\rm d}\Omega}{{\rm d}r}=\mu\left(r\frac{{\rm
d}\Omega}{{\rm d}r}\right)^{2}=\frac{9}{4}\mu\Omega_{\rm K}^{2},$ (120)

where the last equality applies only for a point-mass Keplerian potential.
Integrating over $z$, the heating rate per unit surface area in the disk plane
is,

$Q_{+}=\int_{-\infty}^{\infty}\frac{9}{4}\mu\Omega_{\rm
K}^{2}dz=\frac{9}{4}\nu\Sigma\Omega_{\rm K}^{2}.$ (121)

This heat may result in an increase in the temperature of the disk, and it may
be transported radially by the disk flow. If these effects are negligible and
the energy is radiated locally, the disk effective temperature is,

$2\sigma T_{\rm eff}^{4}=\frac{9}{4}\nu\Sigma\Omega_{\rm K}^{2},$ (122)

where $\sigma$ is the Stefan-Boltzmann constant and the factor of two comes
from the fact that the disk radiates from both its upper and lower surfaces.
This equation does not require that the disk be in a steady state. If the disk
is in a steady state, however, with the profile given by equation (119), the
temperature distribution is,

$T_{\rm eff}^{4}=\frac{3GM\dot{M}}{8\pi\sigma
r^{3}}\left[1-\sqrt{\frac{\tilde{r}}{r}}\right].$ (123)

The steady state temperature profile does not depend on the viscosity.

The temperature profile for a steady viscous disk is not what you get from a
fully local toy model in which available gravitational potential energy is
lost as radiation at every radius. To see this, suppose that mass $\Delta m$
at radius $r$ moves to $r-\Delta r$ while remaining on a near circular orbit.
Half of the liberated potential energy goes into increased kinetic energy, so
the energy available to heat up the gas is $GM\Delta m\Delta r/2r^{2}$. If the
time scale for mass to move inward is $\Delta t$, and the heat is radiated
uniformly from the annulus with total surface area $4\pi r\Delta r$, the
expected temperature profile would be $T_{\rm eff}^{4}=GM\dot{M}/(8\pi\sigma
r^{3})$. At large $r$, this differs by a factor of three from the disk profile
given by equation (123). Viscous torques cause a significant radial
redistribution of energy in accretion disks.

####  III.1.3 Irradiated disks

Figure 10: Geometry for computing the radial temperature profile of a disk
primarily heated by irradiation from a central source, rather than by internal
dissipation.

The temperature profile given by equation (123) applies if the dominant source
of disk heating is internal dissipation. It is also possible for the dominant
source to be external irradiation, for example if the accreting object is a
star. The temperature profile in this limit depends upon the shape of the
disk, which determines the fraction of stellar photons that are absorbed at
each radius. The simplest case is a flat disk in the equatorial plane, that
absorbs all the incident stellar photons and re-emits the energy locally as a
single temperature blackbody.

To compute the resulting $T_{\rm eff}(r)$, consider a surface in the plane of
the disk at distance $r$ from a star of radius $R_{*}$. The star is assumed to
be a sphere of constant brightness $I_{*}$. Setting up spherical polar
coordinates such that the axis of the coordinate system points to the center
of the star, as shown in Figure 10, the stellar flux passing through this
surface is

$F=\int I_{*}\sin\theta\cos\phi{\rm d}\Omega,$ (124)

where ${\rm d}\Omega$ represents the element of solid angle. We count the flux
coming from the top half of the star only (and equate that to radiation from
only the top surface of the disk), so the limits on the integral are,

$\displaystyle-\pi/2<$ $\displaystyle\phi$ $\displaystyle\leq\pi/2$
$\displaystyle 0<$ $\displaystyle\theta$
$\displaystyle<\sin^{-1}\left(\frac{R_{*}}{r}\right).$ (125)

Substituting ${\rm d}\Omega=\sin\theta{\rm d}\theta{\rm d}\phi$, the integral
for the flux is,

$F=I_{*}\int_{-\pi/2}^{\pi/2}\cos\phi{\rm
d}\phi\int_{0}^{\sin^{-1}(R_{*}/r)}\sin^{2}\theta{\rm d}\theta,$ (126)

which evaluates to,

$\displaystyle
F=I_{*}\left[\sin^{-1}\left(\frac{R_{*}}{r}\right)-\left(\frac{R_{*}}{r}\right)\sqrt{1-\left(\frac{R_{*}}{r}\right)^{2}}\right].$
(127)

For a star with effective temperature $T_{*},$ the brightness
$I_{*}=(1/\pi)\sigma T_{*}^{4}$, with $\sigma$ the Stefan–Boltzmann constant
[e.g. 213]. Equating $F$ to the one-sided disk emission $\sigma T_{\rm
eff}^{4}$ we obtain a radial temperature profile,

$\displaystyle\frac{T_{\rm
eff}^{4}}{T_{*}^{4}}=\frac{1}{\pi}\left[\sin^{-1}\left(\frac{R_{*}}{r}\right)-\left(\frac{R_{*}}{r}\right)\sqrt{1-\left(\frac{R_{*}}{r}\right)^{2}}\right].$
(128)

Integrating over radii, we obtain the total disk luminosity,

$\displaystyle L_{\rm disk}$ $\displaystyle=$ $\displaystyle
2\times\int_{R_{*}}^{\infty}2\pi r\sigma T_{\rm eff}^{4}{\rm d}r$ (129)
$\displaystyle=$ $\displaystyle\frac{1}{4}L_{*}.$

A flat disk that extends all the way to the stellar equator intercepts a
quarter of the stellar flux.

The temperature profile given by equation (128) is approximately a power-law
at large radii. Expanding the right-hand-side in a Taylor series in the limit
that $(R_{*}/r)\ll 1$ (i.e. far from the stellar surface), we obtain,

$T_{\rm eff}\propto r^{-3/4},$ (130)

as the limiting temperature profile of a thin, flat, passive disk. For fixed
molecular weight $\mu$ this in turn implies a sound speed profile

$c_{\rm s}\propto r^{-3/8},$ (131)

and a scaling of the geometric thickness with radius,

$\frac{h}{r}\propto r^{1/8}.$ (132)

An irradiated disk therefore flares (i.e. has a concave shape) toward larger
radii. If the disk does flare then the outer regions intercept a larger
fraction of stellar photons, leading to a higher temperature. As a
consequence, a temperature profile $T_{\rm eff}\propto r^{-3/4}$ is the
steepest profile we would expect to obtain for a passive disk.

Irradiation is frequently important for protoplanetary disks, with standard
models (that include self-consistent treatments of disk flaring) having
effective temperature profiles close to $T_{\rm eff}\propto r^{-1/2}$ [105,
48]. It can also be important in high energy accretion environments, for
example in X-ray binaries where irradiation of the outer disk by X-rays from
the inner disk can dominate the local thermal balance [62].

####  III.1.4 Green’s function solution

Assume for simplicity that the viscosity $\nu(\Sigma,r,\ldots)$ is a constant.
The surface density of the disk $\Sigma(r,t)$ then obeys the equation,

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

To solve this equation we first manipulate it into the standard form of a
Bessel’s equation333We’re working toward the famous solution found by Lynden-
Bell and Pringle [140], but here following Gordon Ogilvie’s notes on
“Accretion Disks” from Part III of the Cambridge Mathematical Tripos. The
solution strategy is still not all that obvious, though you might note that we
have a diffusion equation in cylindrical co-ordinates, which is analogous to a
classical example of Bessel’s equation—heat diffusion in a cylinder.. We look
for a solution in which the variables are separated, and modes have a decaying
time dependence,

$\Sigma(r,t)=r^{\beta}\sigma(r)\exp[-\lambda t].$ (134)

Here $\lambda>0$ and by writing the spatial dependence as $r^{\beta}\sigma(r)$
we have given ourselves a free parameter in $\beta$. Substituting,

$-\lambda r^{\beta}\sigma=\frac{3\nu}{r}\frac{\rm d}{{\rm
d}r}\left[r^{1/2}\frac{\rm d}{{\rm d}r}\left(\sigma
r^{\beta+1/2}\right)\right].$ (135)

After evaluating the derivatives and dividing through by $r^{\beta-2}$ we
have,

$r^{2}\frac{{\rm d}^{2}\sigma}{{\rm
d}r^{2}}+\left(2\beta+\frac{3}{2}\right)r\frac{{\rm d}\sigma}{{\rm
d}r}+\beta\left(\beta+\frac{1}{2}\right)\sigma+\frac{\lambda}{3\nu}r^{2}\sigma=0.$
(136)

Defining $k^{2}\equiv\lambda/(3\nu)$ and using the freedom to choose
$\beta=-1/4$,

$r^{2}\frac{{\rm d}^{2}\sigma}{{\rm d}r^{2}}+r\frac{{\rm d}\sigma}{{\rm
d}r}+\left(k^{2}r^{2}-\frac{1}{16}\right)\sigma=0.$ (137)

This is in the form of Bessel’s equation, which has a general solution,

$\sigma=c_{1}J_{1/4}(kr)+c_{2}Y_{1/4}(kr),$ (138)

where $c_{1}$ and $c_{2}$ are constants and $J_{1/4}$ and $Y_{1/4}$ are Bessel
functions of the first and second kinds respectively. The term involving
$Y_{1/4}(kr)$ implies a non-zero torque as $r\rightarrow 0$, so in the case of
a point mass that does not spin up the disk material $c_{2}=0$. The solution
is therefore,

$\Sigma\propto r^{-1/4}J_{1/4}(kr)\exp[-3\nu k^{2}t].$ (139)

The properties of Bessel functions allow us to write a general initial
condition for the surface density in the form,

$\Sigma(r,0)=\int_{0}^{\infty}g(k)r^{-1/4}J_{1/4}(kr)dk.$ (140)

The time-dependent solution will then be,

$\Sigma(r,t)=\int_{0}^{\infty}g(k)r^{-1/4}J_{1/4}(kr)\exp[-3\nu k^{2}t]dk.$
(141)

The problem is thus solved provided that we can determine the decomposition of
the initial surface density into Bessel functions, given by $g(k)$.

To determine $g(k)$ we make use of the Fourier-Bessel (or Hankel) transform
pair. The textbook definition of this pair is,

$\displaystyle g(k)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}f(r)J_{m}(kr)rdr,$ (142) $\displaystyle f(r)$
$\displaystyle=$ $\displaystyle\int_{0}^{\infty}g(k)J_{m}(kr)kdk.$ (143)

Writing equation (140) in this form,

$r^{1/4}\Sigma(r,0)=\int_{0}^{\infty}k^{-1}g(k)J_{1/4}(kr)kdk,$ (144)

the inverse transform is,

$k^{-1}g(k)=\int_{0}^{\infty}s^{1/4}\Sigma(s,0)J_{1/4}(ks)sds.$ (145)

Substituting in equation (141) the general solution is,

$\displaystyle\Sigma(r,t)=r^{-1/4}\int_{0}^{\infty}\int_{0}^{\infty}\Sigma(s,0)J_{1/4}(ks)J_{1/4}(kr)$
$\displaystyle\times\exp[-3\nu k^{2}t]s^{5/4}kdsdk.$ (146)

We express this in the form,

$\Sigma(r,t)=\int_{0}^{\infty}G(r,s,t)\Sigma(s,0)ds,$ (147)

where,

$\displaystyle
G(r,s,t)=s^{5/4}r^{-1/4}\int_{0}^{\infty}J_{1/4}(ks)J_{1/4}(kr)$
$\displaystyle k\exp[-3\nu k^{2}t]dk,$ (148)

is the Green’s function. This integral evaluates to,

$G(r,s,t)=\frac{r^{-1/4}s^{5/4}}{6\nu t}I_{1/4}\left(\frac{rs}{6\nu
t}\right)\exp\left[-\frac{(r^{2}+s^{2})}{12\nu t}\right],$ (149)

with $I_{1/4}$ being a modified Bessel function. It is illustrative to
consider the solution for an initial ring of gas orbiting at $s=r_{0}$. Taking
the initial condition as,

$\Sigma(s,0)=\frac{m}{2\pi r_{0}}\delta(s-r_{0}),$ (150)

with $\delta$ being a Dirac delta function, the solution follows immediately
from equations (147) and (148). It can be written compactly in terms of
dimensionless variables $x$ and $\tau$,

$\displaystyle x$ $\displaystyle=$ $\displaystyle\frac{r}{r_{0}},$ (151)
$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\frac{12\nu}{r_{0}^{2}}t.$
(152)

In terms of these variables,

$\Sigma(x,\tau)=\frac{m}{\pi
r_{0}^{2}}\frac{x^{-1/4}}{\tau}I_{1/4}\left(\frac{2x}{\tau}\right)\exp\left[-\frac{(1+x^{2})}{\tau}\right].$
(153)

The solution described by this equation is shown in Figure 11. It displays
asymmetric diffusion, with mass flowing toward $r=0$ while the conserved
angular momentum is carried by a vanishing fraction of the mass toward
$r=\infty$. Although derived under the restrictive and normally unrealistic
assumption that $\nu$ is constant, these properties are qualitative features
of viscous disk evolution in the case where there is zero-torque at the inner
edge of the disk.

Figure 11: The evolution of a thin ring of gas, initially at $r=r_{0}$, under
the action of a constant viscosity (equation 153). The curves are plotted at
dimensionless times that are multiples of two, starting at $\tau=0.008$ and
going up to $\tau=0.512$.

Although we will not discuss the details here, time-dependent solutions can
also be derived that dispense with the zero-torque assumption [206, 168].
Particularly simple decretion disk solutions [aspects of which were already
discussed in §2.5 of 140] exist if one assumes that an external torque
maintains a $v_{r}=0$ boundary condition at a finite radius $r_{\rm in}$
[198]. It is also possible to derive a relativistic version of the disk
evolution equation [11].

####  III.1.5 Self-similar solution

Lynden-Bell and Pringle [140] also derived a self-similar solution to the disk
evolution equation (equation 114), for the case where the viscosity is a
power-law function of radius,

$\nu\propto r^{\gamma}.$ (154)

If a disk with characteristic size $r_{1}$ at $t=0$ has a surface density
profile of the form,

$\Sigma(t=0)=\frac{C}{3\pi\nu_{1}\tilde{r}^{\gamma}}\exp\left[{-\tilde{r}^{(2-\gamma)}}\right],$
(155)

where $C$ is a constant, $\tilde{r}\equiv r/r_{1}$, and
$\nu_{1}\equiv\nu(r_{1})$, then the time-dependent solution is,

$\displaystyle\Sigma(\tilde{r},T)=\frac{C}{3\pi\nu_{1}\tilde{r}^{\gamma}}T^{-(5/2-\gamma)/(2-\gamma)}\exp\left[-\frac{\tilde{r}^{(2-\gamma)}}{T}\right],$
(156)

where,

$\displaystyle T$ $\displaystyle\equiv$ $\displaystyle\frac{t}{t_{\rm s}}+1,$
(157) $\displaystyle t_{\rm s}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{3(2-\gamma)^{2}}\frac{r_{1}^{2}}{\nu_{1}}.$ (158)

This solution has proved to be quite useful for comparing theoretical models
of viscous disk evolution to data [e.g. 88]. It can be generalized to the case
where disk evolution is driven by a combination of viscous transport and MHD
winds [239].

###  III.2 The $\alpha$-prescription

A predictive model for disk evolution follows from equation (114) if we can
write down how the stress, or equivalently the viscosity, depends on
properties of the disk. Shakura and Sunyaev [225] advanced physical arguments
in favor of the form,

$T_{r\phi}=-\alpha P,$ (159)

where $P$ is the pressure and $\alpha$ is a dimensionless parameter444Shakura
and Sunyaev [225] first introduce $\alpha$ with an expression involving the
magnetic field ($H$ in their notation), similar but not identical to our
equation (162). The famous version is their equation (1.2),
$T_{r\phi}=-\alpha\rho c_{s}^{2}$, where the sound speed includes
contributions from both gas and radiation pressure.. This ansatz is known as
the “$\alpha$ prescription”. Using equation (106) and writing $\nu=\mu/\rho$,
for a Keplerian disk an equivalent form is,

$\nu=\frac{2}{3}\alpha c_{s}h\simeq\alpha c_{s}h.$ (160)

In keeping with the approximate nature of this exercise, the factor of two-
thirds is typically ignored and the $\alpha$-prescription written as just
$\nu=\alpha c_{s}h$.

To order of magnitude, the microphysical viscosity of a fluid can be written
in terms of the thermal velocity of the molecules $v_{\rm th}$ and the mean-
free-path $l$ as,

$\nu\sim v_{\rm th}l.$ (161)

By analogy, one can view equation (160) as describing an effective viscosity
due to turbulent eddies whose speed scales with the sound speed and whose size
scales with the disk thickness. Since turbulent velocities exceeding the sound
speed would cause shocks and rapid dissipation, and isotropic eddies could not
significantly exceed $h$, this argument bounds $\alpha<1$. This argument is
not terribly useful, as physical mechanisms for disk turbulence do not yield
turbulent structures that look much like eddies on scales of the order of $h$.
It is more instructive to follow Shakura and Sunyaev’s original train of
thought, and express $\alpha$ in terms of fluid (Reynolds) and magnetic
(Maxwell) stresses in a turbulent fluid [14],

$\alpha=\left\langle\frac{\delta v_{r}\delta
v_{\phi}}{c_{s}^{2}}-\frac{B_{r}B_{\phi}}{4\pi\rho
c_{s}^{2}}\right\rangle_{\rho},$ (162)

where the angle brackets indicate a density weighted average over space and
time. The first term is the Reynolds stress from correlated fluctuations in
the radial velocity and perturbed azimuthal velocity, the second is the
Maxwell stress from MHD turbulence.

###  III.3 Time scales

For a thin disk we can express several relevant time scales as simple
functions of the disk properties. The dynamical time scale,

$t_{\rm dyn}=\frac{1}{\Omega},$ (163)

is evidently $1/2\pi$ of the orbital period. The time scale for establishing
vertical hydrostatic equilibrium is the sound crossing time across a scale
height, $t_{\rm hydro}\sim h/c_{s}$ Using $h=c_{s}/\Omega$ we have,

$t_{\rm hydro}\sim\frac{h}{c_{s}}\sim\frac{1}{\Omega}\sim t_{\rm dyn}.$ (164)

Vertical hydrostatic equilibrium is thus established in a circular disk on the
shortest possible time scale. In an eccentric disk, however, where the
gravitational potential experienced by a fluid element varies around the
orbit, the approximate equality between $t_{\rm hydro}$ and $t_{\rm dyn}$
means that hydrostatic equilibrium is not established and interesting coupled
dynamics between the radial and vertical structure is possible [179].

The thermal time scale is the time scale on which the disk would cool if
heating processes were instantaneously cut off. The thermal energy per unit
surface area of the disk is $U\sim\Sigma c_{s}^{2}$. Using equation (121) and
equation (160),

$t_{\rm th}=\frac{U}{Q_{+}}\sim\frac{1}{\alpha\Omega}.$ (165)

The thermal time scale is the shortest time scale on which we expect the
emission from an annulus of the disk, heated by viscous-like dissipation, to
change.

The viscous time scale is the time scale on which redistribution of angular
momentum leads to gas inflow. If the surface density is not in a steady-state,
it is also the time scale over which the surface density evolves. Starting
from equation (114), we can estimate the viscous time scale by writing the
evolution equation in a form that more closely resembles a prototypical one-
dimensional diffusion equation. Defining,

$\displaystyle X$ $\displaystyle\equiv$ $\displaystyle 2r^{1/2},$ (166)
$\displaystyle f$ $\displaystyle\equiv$ $\displaystyle\frac{3}{2}\Sigma X,$
(167)

and assuming that $\nu$ is a constant, the evolution equation is,

$\frac{\partial f}{\partial t}=D\frac{\partial^{2}f}{\partial X^{2}},$ (168)

with diffusion coefficient $D$ given by,

$D=\frac{12\nu}{X^{2}}.$ (169)

The diffusion time scale across scale $\Delta X$ implied by equation (168) is
$(\Delta X)^{2}/D$. Returning to the original variables, the time scale over
which viscosity smooths out surface density gradients on radial scale $\Delta
r$ is,

$\tau_{\rm visc}\sim\frac{(\Delta r)^{2}}{\rm\nu}.$ (170)

If the disk has size $r$, the surface density can evolve on a time scale,

$\tau_{\rm visc}\approx\frac{r^{2}}{\nu}.$ (171)

Using the $\alpha$-prescription, we obtain,

$t_{\rm visc}\sim\frac{1}{\alpha\Omega}\left(\frac{h}{r}\right)^{-2}.$ (172)

For a thin disk the viscous time scale is substantially longer than the
thermal time scale, and absent special circumstances thermal equilibrium is
maintained in the vertical direction while the surface density evolves due to
accretion.

Figure 12: Estimates of the dynamical, thermal and viscous time scales at
$r=10GM/c^{2}$ for disks around supermassive black holes of different masses.
The disk aspect ratio is taken to be $h/r=0.1$. The colored bands for the
thermal and viscous time scales show the range obtained assuming
$10^{-2}\leq\alpha\leq 0.1$.

The time scale hierarchy,

$t_{\rm dyn}\sim t_{\rm hydro}\ll t_{\rm th}\ll t_{\rm visc},$ (173)

is a generic property of geometrically thin disks, and is one of the main
reasons why thin disk theory is internally consistent and useful [e.g. 197].
Figure 12 gives a concrete example of these time scales at ten gravitational
radii around supermassive black holes of various masses, for $h/r=0.1$ and
different assumed values for $\alpha$. For a $10^{9}\ M_{\odot}$ black hole,
for example, the dynamical time scale in the inner disk is of the order of a
day, the thermal time scale is around a month, and the viscous time scale is
around ten years.

###  III.4 $\alpha$-model disks

Adopting the $\alpha$-prescription (§III.2) the dependence of $\nu$ on the
local disk conditions and on $\alpha$, $\nu(\alpha,\Sigma,\Omega_{\rm K})$,
can be determined. With this function in hand, equation (114) can be solved
(usually numerically) for the time-dependent evolution of an arbitrary initial
surface density profile. Steady-state solutions (usually analytic) for the
surface density profile $\Sigma(r,\dot{M})$ can also be found. These are
“$\alpha$-model” or “Shakura-Sunyaev” disk solutions. Recall that none of this
effort is necessary if our only concern is the profile of the disk effective
temperature in steady state, as that is given by equation (123) independent of
the form of the viscosity.

A toy example shows how $\alpha$-model disks are constructed. Assume, for no
reason other than simplicity, that the vertical structure of the disk is
isothermal. The effective temperature $T_{\rm eff}$ is then the only
temperature characterizing the disk at some radius, and the viscosity can be
derived from a triplet of already-introduced equations,

$\displaystyle\nu$ $\displaystyle=$ $\displaystyle\alpha c_{s}h,$ (174)
$\displaystyle h$ $\displaystyle=$ $\displaystyle\frac{c_{s}}{\Omega_{\rm
K}},$ (175) $\displaystyle 2\sigma T_{\rm eff}^{4}$ $\displaystyle=$
$\displaystyle\frac{9}{4}\nu\Sigma\Omega_{\rm K}^{2}.$ (176)

The sound speed is related to the temperature through,

$c_{s}^{2}=\frac{k_{\rm B}T_{\rm eff}}{\mu m_{p}},$ (177)

where $k_{\rm B}$ is Boltzmann’s constant, $m_{p}$ is the mass of the proton,
and $\mu$ is the mean molecular weight in units of $m_{p}$. Using this, we
eliminate $T_{\rm eff}$ from equations (174)-(176) and obtain an expression
for the viscosity,

$\nu=\frac{9^{1/3}}{2\sigma^{1/3}}\left(\frac{k_{\rm B}}{\mu
m_{p}}\right)^{4/3}\alpha^{4/3}\Omega_{\rm K}^{-2/3}\Sigma^{1/3}.$ (178)

For fixed central mass $M$, the predicted viscosity scales as
$\nu\propto\alpha^{4/3}r\Sigma^{1/3}$. The equation for the evolution of the
disk surface density, equation (114), would be non-linear with this viscosity.
In steady-state, $\nu\Sigma\propto\dot{M}$ (equation 119), so away from the
inner boundary the disk surface density would scale as
$\Sigma\propto\dot{M}^{3/4}r^{-3/4}$. We have not specified $\alpha$, but as
long as this can be taken to be fixed (determined, perhaps, from simulations
of physical angular momentum transport mechanisms or inferred from
observations of time-dependent disk systems) we have a full solution for the
evolution of geometrically thin disks.

The toy model given above captures the spirit of $\alpha$-model disks, but
it’s not quite the full story. Real disks will not be vertically isothermal.
This extra complexity can be captured at different levels of approximation,

  * (i)

At each radius, characterize the disk’s vertical structure in terms of a
central temperature $T_{c}$ as well as an effective temperature $T_{\rm eff}$.
We can derive a relation between $T_{c}$ and $T_{\rm eff}$ by considering how
energy is transported vertically within the disk. Then, assuming that the
sound speed that enters into the expression for the viscosity ($\nu=\alpha
c_{s}^{2}/\Omega_{\rm K}$) corresponds to the central temperature, we can
proceed as before and derive the functional form of the viscosity. This is
described as a “one-zone” model for the disk vertical structure.

  * (ii)

Alternatively, we can write down and solve (numerically) differential
equations for the vertical disk structure, $\rho(z)$, $T(z)$, in a manner
directly analogous to the equations of stellar structure. This approach
requires a point-by-point specification of the stress, for which we could
adopt the original Shakura and Sunyaev [225] form, $T_{r\phi}=-\alpha P$, or
something else perhaps derived from simulations. Because of the separation
between the thermal and viscous time scales in a geometrically thin disk, it
is normally consistent to solve for the vertical structure separately from the
radial structure. This is sometimes called a “1+1D” disk model.

Both of these approaches are well-defined. Whether the additional complexity
of a vertical integration leads to a physically more realistic model, however,
is an open question. Unlike in the case of stellar structure—where the energy
transport processes and rate of nuclear energy generation are quite well-
known—the physical processes entering into calculations of disk vertical
structure are uncertain. I wouldn’t ascribe much physical reality to modest
differences between disk models with differing formal degrees of
approximation.

To write down a version of the one zone equations, suitable for deducing the
properties of steady-state $\alpha$-disks, we need only a relation between the
central temperature $T_{c}$ and the effective temperature $T_{\rm eff}$. As in
stellar structure, energy can be transported from the hot interior to the
cooler photosphere by radiative diffusion or by convection555Turbulent
transport or transport by waves are also in principle possible. In fact, they
may well be important, and we ignore them here only because they are harder to
capture in simple analytic formulae.. Consider the limit of an optically thick
disk with radiative transport. The vertical energy flux is [e.g. 213],

$F(z)=-\frac{16\sigma T^{3}}{3\kappa_{\rm R}\rho}\frac{{\rm d}T}{{\rm
d}T}=\sigma T_{\rm eff}^{4},$ (179)

where $\kappa_{\rm R}$ is the Rosseland mean opacity (with units of ${\rm
cm}^{2}\ {\rm g}^{-1}$). In equating the flux to a constant we have assumed
that energy dissipation is strongly concentrated at $z=0$. Noting that the
increment of optical depth ${\rm d}\tau=\kappa_{\rm R}{\rm d}z$, we integrate
from the mid-plane to the photosphere,

$-\frac{16}{3}\int_{T_{c}}^{T_{\rm eff}}T^{3}{\rm d}T=T_{\rm
eff}^{4}\int_{0}^{z_{\rm ph}}{\rm d}\tau.$ (180)

If the disk is sufficiently optically thick that $T_{c}\gg T_{\rm eff}$ then
we obtain,

$\left(\frac{T_{c}}{T_{\rm eff}}\right)^{4}\simeq\frac{3}{4}\tau,$ (181)

as the relation between the central and photospheric disk conditions.

With this expression in hand, we can write down a set of equations that
determine the steady-state radial structure of $\alpha$-model disks in the one
zone approximation. The disk is specified by the central mass $M$, accretion
rate $\dot{M}$, and innermost radius $\tilde{r}$, where a zero-torque inner
boundary condition is imposed. The variables to be determined are the mid-
plane density $\rho_{0}$, pressure $P$, temperature $T_{c}$, sound speed
$c_{s}$, surface density $\Sigma$, scale height $h$, optical depth to the mid-
plane $\tau$, Rosseland mean opacity $\kappa_{\rm R}$, and viscosity $\nu$. We
assume an $\alpha$-prescription, with $\alpha$ a constant, and approximate the
opacity as a power-law function of the central density and temperature.
Collecting together previous results (equations 10, 8, 5, 123, 160, and 119),
and adding in an equation of state together with some straightforward
definitions, the set is,

$\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi}}\frac{\Sigma}{h},$ (182) $\displaystyle h$
$\displaystyle=$ $\displaystyle\frac{c_{s}}{\Omega_{\rm K}},$ (183)
$\displaystyle c_{s}^{2}$ $\displaystyle=$ $\displaystyle\frac{P}{\rho_{0}},$
(184) $\displaystyle P$ $\displaystyle=$ $\displaystyle\frac{k_{\rm B}}{\mu
m_{p}}\rho_{0}T_{c}+\frac{4\sigma}{3c}T_{c}^{4},$ (185) $\displaystyle
T_{c}^{4}$ $\displaystyle=$ $\displaystyle\frac{9\dot{M}\Omega_{\rm
K}^{2}}{32\pi\sigma}\left[1-\sqrt{\frac{\tilde{r}}{r}}\right]\tau,$ (186)
$\displaystyle\nu\Sigma$ $\displaystyle=$
$\displaystyle\frac{\dot{M}}{3\pi}\left[1-\sqrt{\frac{\tilde{r}}{r}}\right],$
(187) $\displaystyle\nu$ $\displaystyle=$ $\displaystyle\alpha c_{s}h,$ (188)
$\displaystyle\tau$ $\displaystyle=$ $\displaystyle\frac{1}{2}\kappa_{\rm
R}\Sigma,$ (189) $\displaystyle\kappa_{\rm R}$ $\displaystyle=$
$\displaystyle\kappa_{0}\rho_{0}^{a}T_{c}^{b}.$ (190)

Up to a few not-so-important numerical factors, these are the standard
equations used to determine thin disk structure in the Newtonian limit.
Additional discussion of them can be found in Frank  _et al._ [72].

As written, the mid-plane pressure is the sum of a gas pressure component and
one due to radiation pressure. Usually one or the other of these pressure
sources is much larger than the other, with as a rule of thumb radiation
pressure dominating in black hole disks close to the innermost stable circular
orbit, and gas pressure dominating otherwise666Under some circumstances,
specifically when the disk is threaded by a net magnetic field, magnetic
pressure, $P_{B}=B^{2}/8\pi$, may also contribute to vertical support against
gravity [9, 218, 261, 158]. This possibility is not normally considered in
classical disk models.. Dropping either gas or radiation pressure, we can
solve the set of equations for the steady-state disk structure, verifying
after the fact that we dropped the right one. Away from the inner boundary,
the solutions take the form of power-laws, e.g. $\Sigma\propto
r^{w}M^{x}\dot{M}^{y}\alpha^{z}$, with power-law indices that depend upon the
source of pressure and upon the functional form of the opacity. Because radius
and central mass only enter the equations combined in the form of the
Keplerian angular velocity, the solutions only depend on $\Omega_{\rm K}$.

###  III.5 Self-gravitating $\alpha$-disks

In most cases, and in particular when the source of angular momentum transport
is the MRI, the $\alpha$-disk equations in no way determine $\alpha$ (though
the formalism would be inconsistent for $\alpha$ values large enough to induce
supersonic inflow). Self-gravitating disks can be an exception. Their
stability is a function of the Toomre $Q=c_{s}\Omega_{\rm K}/(\pi G\Sigma)$
(§II.5, specializing to a Keplerian disk). Suppose, somewhat reasonably, that
the angular momentum transport rate is a function of the local disk conditions
and increases rapidly as $Q$ drops below some critical value $Q_{0}\sim 1$.
Under these assumptions, the combination of self-gravitating transport and
local thermal equilibrium can establish a stable feedback loop that maintains
$Q\approx Q_{0}$. If $Q>Q_{0}$, reduced transport leads to reduced heating,
lowering $c_{s}$ to re-establish $Q=Q_{0}$. The reverse happens for $Q<Q_{0}$.
Self-gravitating disks are then expected to be everywhere marginally stable,
with $Q\approx Q_{0}$. This type of model was introduced by Paczynski [183].

Imposing $Q=Q_{0}$ in addition to the usual set of $\alpha$-disk equations has
an important consequence: $\alpha$ is no longer a free parameter but rather a
specified function of the local disk conditions [78, 124]. The evolution of
the disk in this limit is then fully determined. Rafikov [205], and references
therein, detail such “gravito-turbulent” disk models. They are useful provided
that non-local angular momentum transport and mass infall (a complication that
often accompanies self-gravity in astrophysically relevant settings) can be
consistently ignored.

###  III.6 The values of $\alpha$

Although a great deal of effort has been expended over the years in trying to
determine “the” value of $\alpha$, it should be clear from the discussion so
far that this is an illusory quest. Even to the extent that $\alpha$ provides
a good parameterization of the strength of accretion disk turbulence, its
value ought to depend upon the MHD properties of the disk, on the strength of
self-gravity, and so on. The consensus theoretical expectation is that for
disks that are well-described by ideal MHD, turbulence in the absence of net
magnetic flux yields [56, 231],

$\alpha_{\rm ZNF}\simeq 0.01-0.02.$ (191)

There remains some uncertainty about how well converged this computational
result is [212]. In the presence of a net vertical magnetic field $B_{z}$
transport is stronger, with local simulation results indicating that [89, 9,
218],

$\alpha_{\rm VF}\approx\beta_{z}^{-0.5},$ (192)

where $\beta_{z}$ is defined as the ratio of the gas pressure to the magnetic
pressure in the net field at the disk mid-plane,

$\beta_{z}\equiv\frac{P_{\rm gas}}{P_{\rm B}}.$ (193)

This scaling implies that disks with moderately strong vertical fields,
$\beta_{z}\lesssim 10^{2}$, are strongly turbulent and generate mid-plane
toroidal fields with magnetic pressure comparable to the gas pressure. Such
magnetically dominated or “magnetically elevated” disks have stability
properties that differ in interesting ways from standard Shakura-Sunyaev disks
[25].

For self-gravity, local gravito-turbulent models predict that $\alpha$ scales
with the local cooling time as [78],

$\alpha_{\rm SG}\propto\frac{1}{\Omega_{\rm K}t_{\rm cool}},$ (194)

with an upper limit set by fragmentation at $\alpha_{\rm SG}\sim 0.1$ [78,
209]. As with the MRI, pinning down these numbers precisely from numerical
simulations is none too easy a task.

Observationally, $\alpha$ can be estimated in systems where variability
exposes the viscous time scale (equation 170), with the most important example
being dwarf novae, whose disks show limit cycle behavior due to thermal
instability (see §VII.2). Dwarf nova outbursts can be well-described as time-
dependent $\alpha$-disks [154, 21, 155, 233]. Under well-ionized conditions,
modeling of these systems suggests $\alpha\approx 0.1$ [109, 86]. The inferred
larger value of $\alpha$, as compared to predictions from simplified MRI
simulations, may be due to convection in dwarf nova disks [91]. The strength
of turbulence in protoplanetary disks can be constrained more directly from
observations of molecular line broadening [95]. Such analyses suggest that
much lower levels of turbulence (in some cases only upper limits are obtained)
occur in very weakly ionized disks [71, 70].

###  III.7 The Shakura-Sunyaev solution

Thin disk solutions depend upon whether the main source of pressure is gas or
radiation, and upon the opacity under conditions encountered near the disk
mid-plane. For disks around compact objects (black holes, neutron stars, and
white dwarfs) two opacity regimes cover most conditions of interest. At high
temperatures, electron scattering dominates. For plasma with a typical
astrophysical distribution of elements, the opacity is,

$\kappa_{\rm es}=0.34\ {\rm cm^{2}\ g^{-1}}.$ (195)

At lower temperatures, free-free opacity, which can be approximated using
Kramers’ law, applies,

$\kappa_{\rm ff}=6.4\times 10^{22}\rho T^{-7/2}\ {\rm cm^{2}\ g^{-1}}.$ (196)

(Here, $\rho$ and $T$ are understood to be expressed in c.g.s. units.)
Kramers’ law remains valid down to the temperature where hydrogen recombines,
at $T\approx 1-2\times 10^{4}\ {\rm K}$. At lower temperatures, which can be
encountered at large radii in AGN disks and which are typical of
protoplanetary disks, molecules, dust and ice grains provide the opacity.

The properties of Shakura-Sunyaev disk solutions are not terribly intuitive,
but one important result—the disk thickness in the innermost radiation
pressure dominated region—is readily derived. The vertical flux of momentum
carried by radiation, $F_{z}/c$, is equal to (using equation 123),

$\frac{F_{z}}{c}=\frac{\sigma T_{\rm eff}^{4}}{c}=\frac{3}{8\pi}\Omega_{\rm
K}^{2}\dot{M}\left[1-\sqrt{\frac{\tilde{r}}{r}}\right].$ (197)

At high enough temperatures the force that the radiation exerts per unit mass
on the gas is $F_{z}\sigma_{\rm T}/c$, where $\sigma_{\rm T}=6.7\times
10^{-25}\ {\rm cm^{2}}$ is the Thomson cross-section appropriate for electron
scattering. Setting this equal to the vertical acceleration due to the gravity
of the central object, which in the Newtonian limit is just $\Omega_{\rm
K}^{2}z$, the scale height in the radiation pressure / electron scattering
dominated regime is,

$h=\frac{3\sigma_{\rm T}}{8\pi
c}\left[1-\sqrt{\frac{\tilde{r}}{r}}\right]\dot{M}.$ (198)

Away from the inner boundary, the scale height ($h$ itself, not $h/r$) is
predicted to be constant with radius, with a value that is proportional to the
accretion rate.

Without further ado, we quote without derivation the Shakura and Sunyaev [225]
disk solution, which takes the form of piece-wise power-laws corresponding to
three regimes,

  * •

An inner disk, dominated by radiation pressure and electron scattering
opacity.

  * •

A middle disk, dominated by gas pressure and electron scattering opacity.

  * •

An outer disk, dominated by gas pressure and free-free opacity.

The extent of the outer disk is limited by the validity of the free-free
opacity formula—which as noted already fails at low temperature—and / or by
the onset of other physical processes such as gravitational instability. The
model is conceptually just the solution of equations (182-176), though there
is some added subtlety that arises from the fact that electron scattering is a
true scattering process that does not alter the energies of either photons or
electrons. This means that even at high temperatures, sub-dominant absorption
opacity plays a critical role in thermalizing the emergent radiation. There
are also order unity numerical factors that differ between our equation set
and that of the original Shakura and Sunyaev [225] paper.

The Shakura and Sunyaev [225] solutions can usefully be expressed in terms of
dimensionless variables for mass, accretion rate, and radius,

$\displaystyle m$ $\displaystyle\equiv$ $\displaystyle\frac{M}{M_{\odot}},$
(199) $\displaystyle\dot{m}$ $\displaystyle\equiv$
$\displaystyle\frac{\dot{M}}{3\times 10^{-8}(M/M_{\odot})\ M_{\odot}\ {\rm
yr^{-1}}},$ (200) $\displaystyle r^{\prime}$ $\displaystyle\equiv$
$\displaystyle\frac{r}{3r_{\rm S}},$ (201)

where $r_{\rm S}$, the Schwarzschild radius, is given by,

$r_{\rm S}=\frac{2GM}{c^{2}}.$ (202)

These scaling are evidently intended for black hole accretion problems. The
Schwarzschild radius is the radius of the event horizon for a non-rotating
black hole, which has an innermost stable orbit at $r=3r_{\rm S}$. The
accretion rate scaling corresponds roughly to the Eddington limit, which is
also most directly relevant to black hole and other energetic accretion
environments (see IV.1). Nonetheless, these are fundamentally Newtonian
solutions, which can be rewritten with scalings more appropriate to, e.g.,
white dwarfs, without difficulty. Denoting the inner disk with a subscript
$i$, the scale height, surface density, and central temperature are,

$\displaystyle h_{i}$ $\displaystyle=$ $\displaystyle 3.2\times
10^{6}\dot{m}mf\ {\rm cm},$ $\displaystyle\Sigma_{i}$ $\displaystyle=$
$\displaystyle 4.6\alpha^{-1}\dot{m}^{-1}{r^{\prime}}^{3/2}f^{-1}\ {\rm g\
cm^{-2}},$ $\displaystyle T_{ci}$ $\displaystyle=$ $\displaystyle 2.3\times
10^{7}\alpha^{-1/4}m^{-1/4}{r^{\prime}}^{-3/4}\ {\rm K}.$ (203)

In these expressions, $f\equiv\left(1-{r^{\prime}}^{-1/2}\right)$. For the
middle disk, denoted with a subscript $m$,

$\displaystyle h_{m}$ $\displaystyle=$ $\displaystyle 1.2\times
10^{4}\alpha^{-1/10}\dot{m}^{1/5}m^{9/10}{r^{\prime}}^{21/20}f^{1/5}\ {\rm
cm},$ $\displaystyle\Sigma_{m}$ $\displaystyle=$ $\displaystyle 1.7\times
10^{5}\alpha^{-4/5}\dot{m}^{3/5}m^{1/5}{r^{\prime}}^{-3/5}f^{3/5}\ {\rm g\
cm^{-2}},$ $\displaystyle T_{ci}$ $\displaystyle=$ $\displaystyle 3.1\times
10^{8}\alpha^{-1/5}\dot{m}^{2/5}m^{-1/5}{r^{\prime}}^{-9/10}f^{2/5}\ {\rm K}.$
(204)

For the outer disk, denoted with a subscript $o$,

$\displaystyle h_{o}$ $\displaystyle=$ $\displaystyle 6.1\times
10^{3}\alpha^{-1/10}\dot{m}^{3/20}m^{9/10}{r^{\prime}}^{9/8}f^{3/20}\ {\rm
cm},$ $\displaystyle\Sigma_{o}$ $\displaystyle=$ $\displaystyle 6.1\times
10^{5}\alpha^{-4/5}\dot{m}^{7/10}m^{1/5}{r^{\prime}}^{-3/4}f^{7/10}\ {\rm g\
cm^{-2}},$ $\displaystyle T_{co}$ $\displaystyle=$ $\displaystyle 8.6\times
10^{7}\alpha^{-1/5}\dot{m}^{3/10}m^{-1/5}{r^{\prime}}^{-3/4}f^{3/10}\ {\rm
K}.$ (205)

The transition radii between the regimes ($r_{mi}^{\prime}$ for the inner to
middle disk, $r^{\prime}_{om}$ for the middle to outer disk) are given
implicitly by solving,

$\displaystyle r^{\prime}_{mi}$ $\displaystyle=$ $\displaystyle
150\left(\alpha
m\right)^{2/21}\dot{m}^{16/21}\left(1-{r^{\prime}}_{mi}^{-1/2}\right)^{16/21},$
(206) $\displaystyle r^{\prime}_{om}$ $\displaystyle=$ $\displaystyle
6.3\times 10^{3}\dot{m}^{2/3}\left(1-{r^{\prime}}_{om}^{-1/2}\right)^{2/3}.$
(207)

The dependence on $\alpha$ is weak, so these radii are mostly dependent on the
accretion rate and black hole mass (bearing in mind that the latter can vary
across many orders of magnitude).

Figure 13: Transition radii $r^{\prime}_{mi}$ (red: between the inner
radiation pressure dominated disk and the middle gas pressure / electron
scattering dominated disk) and $r^{\prime}_{om}$ (blue: between the middle gas
pressure / electron scattering dominated disk and the outer gas pressure /
free-free opacity dominated disk), plotted for different black hole masses in
the Shakura-Sunyaev solution. The solid lines assumes $\dot{m}=1$, the dashed
lines $\dot{m}=0.1$, in both cases for $\alpha=0.1$.

Figure 13 shows the dependence of $r^{\prime}_{mi}$ and $r^{\prime}_{om}$ as a
function of black hole mass. For a stellar mass black hole with $M=10\
M_{\odot}$, and $\dot{m}=1$, $r^{\prime}_{mi}\approx 140$ (i.e. $840\
GM/c^{2}$) and $r^{\prime}_{om}\approx 6250$. For a supermassive mass black
hole with $M=10^{7}\ M_{\odot}$, again assuming $\dot{m}=1$,
$r^{\prime}_{mi}\approx 540$ and $r^{\prime}_{om}\approx 6250$. In dimensional
units, in the supermassive case the transition from radiation to gas pressure
occurs at about $5\times 10^{15}\ {\rm cm}$, while the transition from
electron scattering to free-free opacity is at about $6\times 10^{16}\ {\rm
cm}$ (0.02 pc). A lower accretion rate of $\dot{m}=0.1$ moves both radii in by
a factor of 5-6. Qualitatively we observe that (1) all three Shakura-Sunyaev
regions are predicted to be present for reasonable parameter choices, and that
(2) radiation pressure is relatively more important for disks around
supermassive black holes as compared to stellar mass examples.

Figure 14: The geometric thickness of Shakura-Sunyaev disks as a function of
dimensionless disk radius $r^{\prime}=r/3r_{\rm S}$ is shown for supermassive
($m=10^{7}$, red) and stellar mass ($m=10$, green) black holes. The solid
lines show high accretion rate solutions ($\dot{m}=1$), the dashed lines a
lower accretion rate ($\dot{m}=0.1$). In all cases $\alpha=0.1$. The geometric
thickness of the inner radiation pressure dominated region depends strongly on
the accretion rate, and can be quite large for high accretion rates. The
middle and outer disk regions are much thinner, especially for disks around
supermassive black holes, where values of $h/r\sim{\rm few}\times 10^{-3}$ are
characteristic. (The small discontinuities in the plotted curves are due to
approximations and do not have any physical significance.)

Figure 14 shows the radial dependence of the predicted geometric thickness of
some selected Shakura-Sunyaev disk solutions. Although these are “thin” disk
solutions, the region where radiation pressure dominates (and $h$ is constant)
is not actually thin at all for $\dot{m}\sim 1$. Values of $h/r$ significantly
in excess of 0.1 are predicted at $r\sim 10-20\ GM/c^{2}$. The resulting
violation of the assumptions underlying thin disk theory is remedied in slim
accretion disk models [1], which should really be used to give a consistent
treatment of this region. Conversely, the gas pressure dominated middle and
outer regions of the Shakura-Sunyaev solution are quite thin, especially in
the case of supermassive black holes. For $M\gtrsim 10^{7}\ M_{\odot}$ we
expect $10^{-3}<h/r<10^{-2}$ in these regions. As a consequence, disk self-
gravity (§II.5, III.5) is predicted to become important for disk masses that
are much smaller than the mass of the black hole. The onset of self-gravity
and the likelihood of ensuing fragmentation, in turn, has far-reaching
consequences for the radial extent of AGN disks, for the formation of stars
and compact objects within the accretion disk, and for how supermassive black
holes are fuelled and grow [229, 84, 108, 124].

Novikov and Thorne [173] generalized the Newtonian Shakura-Sunyaev thin disk
solution to include relativity. The Novikov-Thorne solution does not introduce
any novelties in its treatment of the disk physics, but the proper inclusion
of all the relativistic effects is just as tricky as you might expect. There
are plenty of opportunities for making mistakes (even the original authors, no
slouches when it comes to relativity, didn’t get it quite right). I recommend
Abramowicz and Fragile [2] as a source for the explicit solution, and as a
starting point for reading the literature on Novikov-Thorne disks.

##  IV Energetics of disk accretion

The thin disk solutions are predicated on two assumptions: that the energy
released by accretion can be radiated on a time scale that is short compared
to the local inflow time scale, and that the outgoing radiation has a
negligible impact on the flow dynamics. We now turn to what happens when these
assumptions fail. Radiation becomes dynamically important when the luminosity
of the disk (or that of the central object) reaches the Eddington limit.
Radiative cooling ceases to be efficient both when the accretion rate is very
low, due to plasma physics effects (the regime of radiatively inefficient
accretion), and when it is very high, due to photon trapping (the regime of
hyperaccretion).

###  IV.1 Eddington limit

The Eddington limit is the luminosity at which radiation pressure from a
central point source balances gravity, curtailing spherically symmetric
accretion. Noting that photons of energy $E$ carry momentum $p=E/c$, the
momentum flux at distance $r$ from an isotropic source with luminosity $L$ is
$L/(4\pi cr^{2}$). If the opacity of the gas is $\kappa$, the outward
radiative force per unit mass of gas is,

$f_{\rm rad}=\frac{\kappa L}{4\pi cr^{2}}.$ (208)

Equating to the inward force per unit mass of gravity,

$f_{\rm grav}=\frac{GM}{r^{2}},$ (209)

the Eddington limiting luminosity is,

$L_{\rm Edd}=\frac{4\pi cGM}{\kappa}.$ (210)

At sufficiently high temperatures the opacity is due to Thomson (electron)
scattering, and $\kappa=\sigma_{\rm T}/m_{\rm H}$, where $\sigma_{\rm
T}=6.7\times 10^{-25}\ {\rm cm^{2}}$ is the Thomson cross-section and $m_{\rm
H}=1.66\times 10^{-24}\ {\rm g}$ is the mass of a hydrogen atom. (The
radiative force acts on the electrons, but these are tightly coupled to the
protons electrostatically.) Under these conditions, numerically,

$\displaystyle L_{\rm Edd}$ $\displaystyle=$ $\displaystyle\frac{4\pi cm_{\rm
H}GM}{\sigma_{\rm T}},$ (211) $\displaystyle\simeq$ $\displaystyle 1.2\times
10^{38}\left(\frac{M}{M_{\odot}}\right)\ {\rm erg\ s^{-1}},$ (212)
$\displaystyle\simeq$ $\displaystyle 3.2\times
10^{4}\left(\frac{M}{M_{\odot}}\right)\ L_{\odot}.$ (213)

A modest correction is needed for Thomson scattering on gas that is not pure
hydrogen, with a larger one being in order if the surrounding gas is cool and
the opacity is due to dust. Dust opacity can be relevant to accretion during
massive star formation, and to the dynamics of gas at relatively large
distances from supermassive black holes in AGN.

It is hopefully obvious that the Eddington limit, derived assuming spherical
symmetry, is best-regarded as a characteristic luminosity above which
radiative forces are guaranteed to matter for surrounding accretion flows. As
with the Pirate’s Code in Pirates of the Caribbean, the Eddington limit is
more what you call a guideline than an actual rule.

###  IV.2 Radiative efficiency of thin disks

The Newtonian estimate for the luminosity of accretion at rate $\dot{M}$ onto
a star with mass $M$ and radius $R_{*}$ is just,

$L=\frac{GM\dot{M}}{R_{*}}.$ (214)

For a black hole the presence of an event horizon means that a fraction of the
potential energy liberated by accretion can, in principle, be “lost” across
the horizon, augmenting the mass of the hole as measured at large distance.
The luminosity associated with black hole accretion then depends on the mode
of accretion, and on some uncertain MHD physics close to the hole. The
standard estimate is based upon the assumption that the black hole accretes
from a disk that is geometrically thin, with a zero-torque boundary condition
at an innermost radius that is close to the hole.

Figure 15: The horizon radius, and that of the innermost stable circular orbit
for prograde (solid) and retrograde (dashed) orbits, is plotted as a function
of the dimensionless black hole spin parameter $a_{*}$.

The fiducial estimate for the radiative efficiency of thin disks relies on
various properties of Kerr black holes [106]. For a rotating, uncharged, black
hole, we define the dimensionless spin parameter $a_{*}$ in terms of the mass
$M$ and angular momentum $J$ via,

$a_{*}\equiv\frac{cJ}{GM^{2}}.$ (215)

The spin is limited to $0\leq a_{*}<1$ (or $-1<a_{*}<1$ using negative values
to denote orbits that are counter-rotating with respect to the spin). Amongst
the many important properties of the Kerr metric are the radius of the event
horizon,

$r_{\rm H}=\left(1+\sqrt{1-a_{*}^{2}}\right)\frac{GM}{c^{2}},$ (216)

and the radii of the innermost stable circular orbits, which differ depending
on whether the orbits are co-rotating or counter-rotating in the equatorial
plane of the hole [19, 228],

$\displaystyle r_{\rm ISCO}$ $\displaystyle=$
$\displaystyle\left[3+Z_{2}\mp\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})}\right]\frac{GM}{c^{2}},$
$\displaystyle Z_{1}$ $\displaystyle\equiv$ $\displaystyle
1+\left(1-a_{*}^{2}\right)^{1/3}\left[\left(1+a_{*}\right)^{1/3}+\left(1-a_{*}\right)^{1/3}\right],$
$\displaystyle Z_{2}$ $\displaystyle\equiv$
$\displaystyle\sqrt{3a_{*}^{2}+Z_{1}^{2}}.$ (217)

These quantities are plotted in Figure 15. The innermost stable circular orbit
lies at $6GM/c^{2}$ for a non-rotating (Schwarzschild) black hole, and tends
towards $GM/c^{2}$ and $9GM/c^{2}$ for co-rotating and counter-rotating orbits
respectively as $a_{*}\rightarrow 1$. The variation of $r_{\rm ISCO}$ with
$a_{*}$ provides the physical underpinning for various black hole spin
estimates derived from electromagnetic observables [208].

Figure 16: The nominal radiative efficiency of thin disk accretion is plotted
as a function of the dimensionless black hole spin parameter $a_{*}$. Negative
values of $a_{*}$ denote accretion from equatorial counter-rotating orbits,
positive values accretion from equatorial co-rotating orbits. For $a_{*}=0$
the efficiency $\eta=1-\sqrt{8/9}=0.0572$, for $a_{*}=-1$ we have
$\eta=1-\sqrt{25/27}=0.0377$, while for $a_{*}=1$ (this value of the spin is
not quite physically realizable) $\eta=1-\sqrt{1/3}=0.423$.

Armed with these results, we can define the fiducial estimate of the radiative
efficiency properly. Assume that a thin disk extends from large radius, where
the mass accretion rate is $\dot{M}$, down to the innermost stable circular
orbit, where a zero-torque boundary condition is applied. Gas interior to the
disk—between $r_{\rm ISCO}$ and $r_{\rm H}$—is assumed to neither radiate nor
to exert any feedback effects on the disk flow. The radiative efficiency
$\eta$, defined via,

$L=\eta\dot{M}c^{2},$ (218)

is then fully specified by the value of the binding energy of particle orbits
at $r_{\rm ISCO}$. Defining the ancillary variable $E=1-\eta$ the efficiency
is found by solving the equation [228],

$a_{*}=\mp\frac{4\sqrt{2}\left(1-E^{2}\right)^{1/2}-2E}{3\sqrt{3}\left(1-E^{2}\right)}.$
(219)

The resulting function $\eta(a_{*})$ is shown in Figure 16. For $a_{*}=0$ the
efficiency is 5.7% (i.e. that fraction of the rest mass of the accreting gas
is radiated from the disk), while as $a_{*}\rightarrow 1$ the limiting
efficiency is 42%.777Although mechanically it is possible to spin up a Kerr
black hole arbitrarily close to $a_{*}=1$, a black hole spun up by disk
accretion necessarily consumes disk photons on orbits that counteract the spin
up. Thorne [241] found that this effect limits the maximum spin to
$a_{*}\simeq 0.998$, with a radiative efficiency $\eta\simeq 0.3$.

####  IV.2.1 Salpeter time

The results for the Eddington limit and the radiative efficiency of disk
accretion can be combined to give an estimate of how rapidly black holes can
grow from thin disk accretion. To do so, assume that the accretion rate is
limited to the value that would yield an Eddington-limited luminosity. For a
black hole of mass $M_{\rm BH}$,

$\eta\dot{M}c^{2}=\frac{4\pi cm_{\rm H}G}{\sigma_{\rm T}}M_{\rm BH}.$ (220)

Subtracting the rest-mass equivalent of the energy lost via radiation,
accretion at rate $\dot{M}$ increases the mass of the black hole according to,

$\dot{M}_{\rm BH}=\left(1-\eta\right)\dot{M}.$ (221)

A black hole that always grows as fast as it can—at the Eddington limit—then
obeys,

$\frac{{\rm d}M_{\rm BH}}{{\rm
d}t}=\frac{4\pi\left(1-\eta\right)}{\eta}\frac{m_{\rm H}G}{c\sigma_{\rm
T}}M_{\rm BH}.$ (222)

At fixed radiative efficiency (i.e. at fixed spin) the result is exponential
growth,

$M_{\rm BH}=M_{0}\exp\left[\frac{t}{t_{\rm S}}\right],$ (223)

where $M_{0}$ is some initial mass and $t_{\rm S}$, the Salpeter time [217],
is a characteristic time scale,

$t_{\rm S}=\frac{\eta}{4\pi\left(1-\eta\right)}\frac{c\sigma_{\rm T}}{m_{\rm
H}G}.$ (224)

For $\eta=0.1$, the Salpeter time is about 50 Myr.

The Salpeter time is a useful characteristic time scale for how fast black
holes accreting from geometrically thin disks grow. You will often see it
referenced as part of an argument about supermassive black hole formation,
which goes as follows. We observe luminous quasars, likely hosting black holes
with masses $M_{\rm BH}\sim 10^{9}\ {M_{\odot}}$, at redshifts exceeding $z=7$
[160, 247]. If we assume (say) that the seed for such a quasar formed at
$z=20$, the time available between $z=20$ and $z=7$ for it to grow is only
approximately $5.9\times 10^{8}\ {\rm yr}$. This is 11.6 Salpeter times (for
$\eta=0.1$), while 18 e-foldings are needed to grow from a $10\ M_{\odot}$
stellar mass black hole to a $10^{9}\ M_{\odot}$ supermassive black hole. The
observed early growth of massive black holes in the Universe appears to pose
problems, or at least to provide strong constraints on the masses of black
hole seeds [which may not be stellar mass black holes at all; 28].

This argument, in my opinion, is rather tired from overuse, and in its strong
form requires undue faith in the Eddington limit being a strict limit.
However, it does justify the weaker statement that the existence of high
redshift quasars requires a high duty cycle of prior accretion, at a rate high
enough to imply that radiative forces are important.

A separate argument, credited to Soltan [234], takes off from the observation
that the mass accumulated in supermassive black holes, and the total amount of
energy radiated from accretion during their growth, are two sides of the same
coin. If supermassive black holes grow primarily via thin disk accretion, the
total mass in black holes per comoving Mpc3 at $z=0$ is related to the
integral of the AGN luminosity function over redshift, and an appropriate
comparison between the two constrains the radiative efficiency [255].

###  IV.3 Electron-ion coupling in low density plasmas

The validity of thin disk solutions for black hole accretion is bounded above
at $\dot{m}\sim 1$, both by the onset of radiatively driven outflows as the
Eddington limit is exceeded and by the importance of radial advection of
energy as $h/r$ becomes larger. In the other direction, thin disk solutions
remain internally consistent for $\dot{m}\ll 1$, but they are not unique. A
hot, geometrically thick disk solution also exists below a critical accretion
rate. The physical origin of the hot solution is tied to an asymmetry in the
microphysics of electron-ion plasmas. Heating, due to small-scale dissipation
of turbulent energy, gives energy predominantly to the ions, while cooling,
due to processes such as free-free emission and synchrotron radiation, is much
more efficient for electrons. At low density, the time scale for Coulomb
collisions to transfer energy from the ions to the electrons becomes long—in
some cases longer than the time scale for the gas to accrete on to the black
hole. Unable to cool efficiently, the accretion flow becomes (or remains)
geometrically thick [227, 97, 207]. It resembles a torus or a doughnut more
than a disk.

Figure 17: Setup for computing the transverse momentum imparted by an
electron-proton encounter in the impulse approximation.

The time scale for electron-proton thermal equilibration can be computed to
better accuracy than we need via elementary methods. Step one is to calculate
the transverse momentum that is imparted when an electron, with mass $m_{e}$,
flies by a proton, with mass $m_{p}$, at velocity $v_{e}$ and impact parameter
$b$. The setup is shown in Figure 17. We work in the impulse approximation,
and compute the momentum change along the unperturbed straight-line trajectory
of the electron. The force perpendicular to the direction of motion is,

$F_{\perp}=\frac{e^{2}}{d^{2}}\cos\theta=\frac{e^{2}b}{\left(b^{2}+x^{2}\right)^{3/2}}.$
(225)

The total momentum imparted in the perpendicular direction as a result of the
encounter is just the integral of the force over time,

$\Delta p_{\perp}=\int_{-\infty}^{\infty}F_{\perp}(t){\rm
d}t=\frac{e^{2}}{b^{2}v_{e}}\int_{-\infty}^{\infty}\left(1+x^{2}/b^{2}\right)^{-3/2}{\rm
d}x.$ (226)

A quick trignometric substitution, or an even quicker trip to Wolfram Alpha,
and we have the answer,

$\Delta p_{\perp}=\frac{2e^{2}}{bv_{e}}.$ (227)

With this result in hand we proceed to step two, which is to sum up the effect
of many encounters occurring with a range of impact parameters $b_{\rm
min}<b<b_{\rm max}$. The mean-square momentum change is the lowest order
quantity that does not vanish. If the number density of particles is $n$,

$\displaystyle\left\langle\frac{\rm d}{{\rm d}t}\left(\Delta
p_{\perp}\right)^{2}\right\rangle$ $\displaystyle=$ $\displaystyle\int_{b_{\rm
min}}^{b_{\rm max}}2\pi bnv_{e}\left(\frac{2e^{2}}{bv_{e}}\right)^{2}{\rm d}b$
(228) $\displaystyle=$ $\displaystyle\frac{8\pi ne^{4}}{v_{e}}\int_{b_{\rm
min}}^{b_{\rm max}}\frac{{\rm d}b}{b}.$

The integral gives us $\ln(b_{\rm max}/b_{\rm min})$, which goes by the name
of the Coulomb logarithm and is written as $\ln\Lambda$. The same quantity
crops up, for the same reasons, in discussions of relaxation in gravitating
point-mass systems such as globular clusters. In a plasma, $b_{\rm max}$
should be set to be the Debye length [e.g. 242],

$\lambda_{D}=\left(\frac{kT}{4\pi ne^{2}}\right)^{1/2},$ (229)

because this is the spatial scale on which charges are effectively screened.
The lower limit, $b_{\rm min}$, is either the impact parameter where a single
deflection $\Delta\theta\sim 1$ radian (as in the gravitational case), or the
particle’s de Broglie wavelength, whichever is the larger [236]. All this can
be calculated properly as a function of density and temperature, but for our
purposes the details are neither very interesting nor very important. The
answer depends on $b_{\rm max}/b_{\rm min}$ only logarithmically, and taking
$\ln\Lambda\approx 20$ is a good enough approximation.

For the final step, we switch from thinking about momentum to thinking about
energy. The electron, with kinetic energy $(1/2)m_{e}v_{e}^{2}$, changes the
energy of the proton by an amount $(\Delta p_{\perp})^{2}/2m_{p}$ in a single
encounter. Considering all encounters, the equilibration time is,

$t_{ep}=\frac{m_{e}v_{e}^{2}}{\langle{\rm d}/{\rm d}t(\Delta
p_{\perp})^{2}/m_{p}\rangle}=\frac{m_{e}m_{p}v_{e}^{3}}{8\pi
ne^{4}\ln\Lambda}.$ (230)

Let us assume that the typical electron velocity is,

$v_{e}=\left(\frac{3k_{\rm B}T_{e}}{m_{e}}\right)^{1/2}.$ (231)

Then, noting that the Thomson cross-section [213],

$\sigma_{\rm T}=\frac{8\pi}{3}\frac{e^{4}}{m_{e}^{2}c^{4}},$ (232)

and defining a parameter expressing how relativistic the electrons are,

$\theta_{e}\equiv\frac{k_{\rm B}T_{e}}{m_{e}c^{2}},$ (233)

equation (230) becomes,

$t_{ep}=\frac{\sqrt{3}}{\sigma_{\rm
T}cn\ln\Lambda}\left(\frac{m_{p}}{m_{e}}\right)\theta_{e}^{3/2}.$ (234)

The rather sloppy methods we have employed here get us surprisingly close to
the answer obtained by more diligent investigators. A two-temperature
electron-proton plasma, where the two species have Maxwellian velocity
distributions, comes to a common temperature on a time scale888In the plasma
context this result is usually cited as Spitzer [236], but Spitzer’s textbook
passes the buck to his paper on stellar dynamics [235] for the details of the
calculation. That’s fine because, up to the question of the correct value of
$\ln\Lambda$, the gravitational and plasma calculations are the same.,

$t_{ep}=\frac{\sqrt{2\pi}}{2\sigma_{\rm
T}cn\ln\Lambda}\left(\frac{m_{p}}{m_{e}}\right)\left(\theta_{e}+\theta_{p}\right)^{3/2},$
(235)

where $\theta_{p}\equiv k_{\rm B}T_{p}/m_{p}c^{2}$. When
$\theta_{e}\gg\theta_{p}$ this differs from our rough and ready result by only
a modest numerical pre-factor.

###  IV.4 Radiatively inefficient accretion models

A geometrically thick disk solution is possible if the density is low enough
that the time scale for electron-ion thermal equilibration exceeds the time
scale for gas to flow into the black hole. This can be written as an upper
limit on the accretion rate, in units of the accretion rate that would
generate an Eddington limiting luminosity, that depends only on $\alpha$. To
obtain this result, we consider a quasi-spherical but rotationally supported
accretion flow, such that $n=n(r)$ only but $v_{r}$ is set by the local rate
of angular momentum transport. The continuity equation is,

$\dot{M}=-4\pi r^{2}v_{r}\rho.$ (236)

Substituting the order of magnitude estimates, $v_{r}=-r/t_{\rm visc}$ and
$t_{\rm visc}\sim 1/(\alpha\Omega_{\rm K})$ (remember that $h/r\sim 1$, by
assumption), and using $\rho=nm_{\rm H}$, gives,

$n(r)\sim\frac{\dot{M}}{4\pi m_{\rm H}\alpha\Omega_{\rm K}r^{3}}.$ (237)

We now substitute this expression for $n$ into equation (235) and impose the
physical condition that,

$t_{ep}\gtrsim t_{\rm visc}.$ (238)

Although the electrons cool to lower temperatures than the protons, they’re
also much less massive, so we can assume that $\theta_{e}>\theta_{p}$ and drop
the latter. The inequality gives,

$\dot{M}\lesssim\frac{\sqrt{2\pi}}{2\ln\Lambda}\frac{4\pi m_{\rm
H}GM\alpha^{2}}{\sigma_{\rm
T}c}\left(\frac{m_{p}}{m_{e}}\right)\theta_{e}^{3/2},$ (239)

which is approximately as clear as mud. The expression simplifies when we note
that many of the quantities on the right-hand-side also appear in the formula
for the Eddington limit (equation 213). Let us define an Eddington mass
accretion rate via,

$L_{\rm Edd}=\eta\dot{M}_{\rm Edd}c^{2},$ (240)

such that,

$\dot{M}_{\rm Edd}=\frac{4\pi m_{\rm H}GM}{\eta c\sigma_{\rm T}}.$ (241)

(Note that some authors define $\dot{M}_{\rm Edd}$ without the factor of
$\eta$, or fixing $\eta=0.1$.) The inequality then simplifies to,

$\frac{\dot{M}}{\dot{M}_{\rm
Edd}}\lesssim\sqrt{\frac{\pi}{2}}\frac{\eta}{\ln\Lambda}\left(\frac{m_{p}}{m_{e}}\right)\theta_{e}^{3/2}\alpha^{2}.$
(242)

It’s not easy to give an elementary estimate for $\theta_{e}$, but adopting
reasonable values for the parameters, $\theta_{e}=0.1$, $\eta=0.1$, and
$\ln\Lambda=20$, the result is,

$\frac{\dot{M}}{\dot{M}_{\rm Edd}}\lesssim 0.4\alpha^{2}.$ (243)

Our argument has been a simplified version of that given in Mahadevan and
Quataert [145], but there are multiple routes to this result. If $\alpha\sim
0.1$, the conclusion is that a geometrically thick, two-temperature accretion
flow, is a consistent possibility provided that the accretion rate, measured
in Eddington units, is below $10^{-3}$–$10^{-2}$.

The way in which plasma microphysics and flow macrophysics combine to give a
simple result for the critical accretion rate is quite appealing. It should be
obvious that our derivation is only good to an order of magnitude, but there
are also other caveats. The fact that a hot two-temperature flow can exist
below a critical accretion rate does not mean that it must—in principle a
denser thin disk could be a valid solution below the critical rate. There has