Foundations of Black Hole Accretion Disk Theory

Abstract

This review covers the main aspects of black hole accretion disk theory. We begin with the view that one of the main goals of the theory is to better understand the nature of black holes themselves. In this light we discuss how accretion disks might reveal some of the unique signatures of strong gravity: the event horizon, the innermost stable circular orbit, and the ergosphere. We then review, from a first-principles perspective, the physical processes at play in accretion disks. This leads us to the four primary accretion disk models that we review: Polish doughnuts (thick disks), Shakura-Sunyaev (thin) disks, slim disks, and advection-dominated accretion flows (ADAFs). After presenting the models we discuss issues of stability, oscillations, and jets. Following our review of the analytic work, we take a parallel approach in reviewing numerical studies of black hole accretion disks. We finish with a few select applications that highlight particular astrophysical applications: measurements of black hole mass and spin, black hole vs. neutron star accretion disks, black hole accretion disk spectral states, and quasi-periodic oscillations (QPOs).

1 Introduction

Because of its firm connection to black holes themselves, black hole accretion disk theory belongs to the realm of fundamental physics. However, the theory itself employs a complicated maze of fluid-dynamics results and several phenomenological estimates and guesses known only to specialists. This Living Review aims to give readers a useful guide, a “road map” if you will, through the unavoidable complexity of the subject.

Below we list fourteen “Destinations” on this road map and explain their logical connections. In our opinion, these Destinations are the key issues, or landmarks, of the theory being reviewed. The particular road map that we present is, of course, biased by our own ideas and research histories. However, we are confident that the grand landscape of this field will, nevertheless, shine through in the end.

We start by pointing out that black holes are one of the most remarkable inventions of the human mind. Their bizarre properties capture nearly everyone’s imagination, from Princeton string theorists to Hollywood movie makers. Initially, black holes only existed virtually, as weird theoretical and mathematical ideas. Models of them were studied with interest, but their real existence was questioned by many, including Albert Einstein himself. This situation changed in the latter part of the twentieth century, after unambiguous and robust detections were made of several astrophysical black hole “candidates” within our Galaxy and in many others. Thus far, two¹ classes of black holes have been observed by astronomers. In our theory-minded Living Review we do not give detailed descriptions of their observational properties. Instead, we stress their importance by starting our road map from the two classes of observed black holes:

• Destination 1: Quasars and other similar supermassive objects, which are collectively called “active galactic nuclei” (or AGN), having masses² in the range 10⁶ M_⊙ < M < 10⁹ M_⊙. They reside at centers of our and other galaxies. The “active” ones among them are the most powerful steady energy sources known in the universe. Many have radiant powers L in excess of their corresponding Eddington luminosities.³ The high efficiency of quasars, η = L/Mc² > 0.1, where Ṁ is the mass supply rate from accretion (“the accretion rate”), is puzzling. Black hole accretion disk theory predicts that L > L_Edd would imply small accretion efficiency η ≪ 0.1. However, the famous “Sołtan argument,” based on quasar counts, shows that on a long time average, t ∼ t_Hubble, quasars can have both L > L_Edd and η ∼ 0.1 [290, 254]. AGN are described in much greater detail in the authoritative monograph by Krolik [163].

• Destination 2: Microquasars and similar “stellar-mass” black holes, having M ∼ 10 M_⊙. The term “microquasars” was invented by Mirabel [206] to convey that these objects, in many regards, behave like scaled-down versions of quasars. A few tens of them have been found in our Galaxy as members of X-ray binaries [184]. The natural scaling (time) ∼ (mass) adds importance to observations of microquasar variability, because the same processes that takes hundreds of years (say) in quasars, takes only minutes in microquasars. Particularly interesting are spectral state changes, which occur on timescales ∼ 1 day (see [258] for a review), and quasi-periodic oscillations (QPOs), which have timescales as short as ∼ 1 ms [308, 184, 258]. Quasars and microquasars are among the most intriguing astrophysical objects ever discovered, and the goal of black hole accretion disk theory is, obviously, to explain their observed properties — but one hopes for much more! One hopes that observations of quasars and microquasars, together with their proper theoretical interpretation, would eventually test the very heart of black hole physics itself. When this happens, we may meaningfully constrain our knowledge of the fundamental properties of space and time.⁴

To accomplish this goal, black hole accretion disk theory must find ways to filter out from the observational data the parts that bear clear black hole signatures from the remaining dirty astrophysics parts. The logical structure of our Living Review attempts to reflect as closely as possible this main task of the theory. Thus, our review starts (in Section 2) with the pure gravity part, concentrating on the three particular signatures of black hole gravity that are crucial to black hole accretion disk theory:

• Destination 3: Event horizon. This is a sphere of radius ∼ GM/c² surrounding the black hole singularity, from within which nothing may emerge — a one-way membrane. Note that this means that black holes have no rigid surfaces. This is a unique signature of black holes; other relativistic features may be observable around non-black hole objects, specifically sufficiently compact neutron stars, but the event horizon is a defining property of black holes.

• Destination 4: Ergosphere. This is a region around a rotating black hole where spacetime itself is dragged along in the direction of rotation at a speed greater than the local speed of light in relation to the rest of the universe. In this region, negative energy states are possible, which means that the rotational energy of the black hole can be tapped through various manifestations of the “Penrose process” [242].

• Destination 5: Marginally stable orbit (also called the “innermost stable circular orbit” or ISCO). This is the smallest circle (r = r_ms) along which free particles may stably orbit around a black hole. No stable circular motion is possible for r < r_ms. This is a unique feature of relativity, as in Newtonian theory, orbits at all radii are possible.⁵

As fascinating as these destinations are, they are not easily probed in nature. The radiation we receive from quasars and microquasars comes not from the black holes themselves, but instead originates in the accretion disks which surround them⁶ (see Figure 1). In these accretion disks, angular momentum is gradually removed by some presumably (although not necessarily [48]) dissipative process, causing matter to spiral down into the black hole, converting its gravitational energy into heat, and then, by various processes, radiating this energy.⁷ The radiation subsequently leaks through the disk, escapes from its surface, and travels along trajectories curved (in space) by the strong gravity of the black hole, eventually reaching our telescopes. Bingo! That is what we are interested in this Living Review.

Figure 1

An artist’s rendition of a generic black hole accretion disk and jet. Inset figures include a time sequence of radio images from the jet in microquasar, GRS 1915+105 [204] and an optical image of the jet in quasar, M87 (Credit: J.A. Biretta et al., Hubble Heritage Team (STScI / AURA), NASA).

The preceding description of how black hole accretion works is well established in general, yet several rather crucial details remain either not sufficiently understood or are too complex to be studied even by the most powerful present-day computers. To deal with this seemingly hopeless situation, several purely phenomenological approaches have been adopted. Two of them should certainly be marked as Destinations, in and of themselves, on the road map, because they are among the very basic ingredients of the accretion process:

• Destination 6: Angular momentum transport and energy dissipation. The quasi-steady accretion of a particle of mass m through a Keplerian disk from a large outer radius, r_out, to an inner radius, r_in, requires that the particle give up an amount of energy⁸ ∼ 0.1 mc². To do this, the particle must also give up an amount of angular momentum ∼ (G M r_out)^1/2. We will see in Section 3.2 that viscous stresses within the fluid can facilitate this (mass transfer in, angular momentum transfer out, and energy dissipation). However, the stresses can not come from ordinary molecular viscosity, as this is much too weak in astrophysical accretion disks. Instead, the stresses likely come from turbulence that acts like an effective viscosity.

• Destination 7: Radiative processes and radiative transfer. These depend on the thermodynamic state of matter (electron density, ion density, temperature), its motion, and the magnetic field, but most importantly on whether the matter is opaque or transparent to radiation, i.e., whether its optical depth is large, τ ≫ 1, or small, τ ≪ 1. In the general case when the matter is optically thick τ > 1, the accretion disk can be quite luminous and also efficiently cooled by radiation. Accretion disks with τ < 1 are inefficiently cooled and thus less luminous.

With these complexities, it is useful to have both analytic theories and numerical simulations to facilitate progress in the field. Before getting into details, though, let us simply study the parameter space available before us. Let us first consider three generic types of physical processes that must occur in black hole accretion disks: “dynamical” processes with a characteristic time-scale, t_dyn ∼ 1/Ω, where Ω is the orbital angular velocity; “thermal” processes with a characteristic time-scale, ${t_{{\rm{th}}}} \sim c_{\mathcal S}^2/{\nu _\ast}{\Omega ^2}$, where c_s is the sound speed and ν_* is the kinematic viscosity; and “viscous” processes with a characteristic time-scale, t_vis ∼ r²/ν_*, where r is the radial distance from the black hole. One excepts that, typically,⁹

$${t_{{\rm{dyn}}}} \ll {t_{{\rm{th}}}} \ll {t_{{\rm{vis}}}}.$$

(1)

Because dynamical processes act much faster than thermal or viscous ones, in the first approximation one may consider only the dynamical structure of the disk. The disk dynamical structure is governed by three forces: gravity, which at a given place is fixed by the Kerr geometry; pressure; and rotation forces. The relative importance of each of these depend on the particular astrophysical situation — their values may either be “small” (dynamically unimportant) or “large” (dynamically important). Correspondingly, there are four “genuine” types of dynamical situations, as described in the upper half of Table 1. The four types of accretion disks we consider in this Living Review (thin, slim, thick, ADAF) are located in three of the four regions. The fourth region corresponds to a “free-fall” of dust.¹⁰ It is remarkable that the same types of accretion disks pop-up when one introduces other, very different, criteria for dividing the parameter space, such as accretion rate and opacity, as shown in lower half of Table 1. These four types of accretion disks are the next four Destinations on our road map. We characterize them in terms of: relative thickness, h = H/R; dimensionless accretion rate, ṁ = 0.1 Ṁc²/L_Edd; optical depth, τ; importance of advection, q = Q_adv/Q_rad, where Q represents an energy flux; importance of radiation pressure, β = P_gas/(P_gas + P_rad); location of inner edge, r_in; and accretion efficiency η.

• Destination 8: Thick Disk.

  $$h > 1,\;\dot m \gg 1,\;\tau \gg 1,\;q\sim 1,\;\beta \ll 1,\;{r_{{\rm{in}}}}\sim {r_{{\rm{mb}}}},\;\eta \ll 0.1$$

• Destination 9: Thin Disk.

  $$h \ll 1,\;\dot m < 1,\;\tau \gg 1,\;q = 0,\;\beta \sim 1,\;{r_{{\rm{in}}}} = {r_{{\rm{ms}}}},\;\eta \sim 0.1$$

• Destination 10: Slim Disk.

  $$h\sim 1,\dot m \gtrsim 1, \;\tau \gg 1,\;q\sim 1,\;\beta < 1,\;{r_{{\rm{mb}}}} < {r_{{\rm{in}}}} < {r_{{\rm{ms}}}},\;\eta < 0.1$$

• Destination 11: Advection-Dominated Accretion Flow (ADAF).

  $$h < 1,\;\dot m \ll 1,\;\tau \ll 1,\;q\sim 1,\;\beta = 1,\;{r_{{\rm{mb}}}} < {r_{{\rm{in}}}} < {r_{{\rm{ms}}}},\;\eta \ll 0.1$$

Table 1 Upper half: The four a priori possible types of dynamical states of accretion structures, corresponding to the division of the dynamical parameter space into fast-slow rotation and large-small pressure. Lower half: A different division of the parameter space, corresponding to high-low accretion rate and large-small opacity.

Interestingly, the parameter space of each of these types of accretion disks overlaps that of other solutions. For instance, ADAF solutions exist that have the same mass accretion rates as thin disk solutions [59]. In such cases, it is not clear how nature might choose one over the other.

In Sections 4–7, we will demonstrate that analytic (or at least semi-analytic) models exist for each of these types of accretion disks. However, this does not guarantee that they are stable. The issue of stability is our next destination.

• Destination 12: Stability. Stability analysis is important because the systematic differential rotation that is one of the defining characteristics of accretion disks is also a potential source of destabilizing energy. On the one hand, this may be essential, as the angular momentum transport and energy dissipation required for accretion may require disks to be mildly unstable. On the other hand, if a model is violently unstable, then the basic assumption of a “steady-state” would be violated.

Along with stability, we also look at the natural oscillation modes associated with accretion disks. The frequencies of some of these modes are tied to the properties of the black hole space-time, which is how these relate to the fundamental physics issues of interest in this Living Review.

• Destination 13: Oscillations. As with any finite distribution of fluid, accretion disks have natural oscillation modes associated with them. If these modes can be excited at appreciable amplitudes, they may be able to modify the observed light curve of the disk in measurable ways. This makes disk oscillations a leading candidate for explaining the quasi-periodic oscillations (QPOs) that we discuss in Section 12.4. Because of the close observational links between black hole accretion disks and jets [92, 94], we include jets as the final Destination of our review.

• Destination 14: Jets. Jets are narrow (opening angle < 5°), long (length > 10⁷ ly in the case of AGN), and fast (v > 0.9c) streams of matter emerging from very compact regions around the black hole, usually in opposite directions, presumably normal to the plane of the accretion disk. Jets can play a significant role in transporting energy and angular momentum away from the accretion disk [48]. They also play an important role in shaping the black hole’s environment far beyond the gravitational reach of the black hole itself, affecting galactic evolution, particle acceleration, and intragalactic ionization.

Going hand-in-glove with analytic models of accretion disks are direct numerical simulations. Although analytic theories have been extremely successful at explaining many general observational properties of quasars and microquasars, numerical simulations can be critically important in at least two respects: 1) as an extension of analytic work, by treating nonlinear perturbations and higher order coupling terms, and 2) in cases that are highly time variable or contain little symmetry, such that the prospects of finding an analytic solution are poor. There is also an important overlap region where various analytic and numerical methods are applicable and can be used to independently validate results. Because of these close connections between analytic and numerical work, we have dedicated Section 11 to the discussion and review of direct numerical simulation of black hole accretion disks.

We finish this Living Review with Section 12, which tries to make some connections between the concepts discussed in earlier sections and actual observational phenomena. We emphasize that we are not aiming to provide a comprehensive review of black hole observations, but rather to highlight a very small subset of these that are of particular relevance to our review.

Throughout this review, we adopt the − + + + metric signature and often use units where c = 1 = G. To make all physical quantities dimensionless, we sometimes also use the mass of the black hole as a unit, M = 1. We use the common Einstein summation convention, where repeated indices in a formula imply summation over the range of that index. We also follow the common convention where Greek (Latin) indices are used for four-(three-)dimensional tensor quantities.

2 Three Destinations in Kerr’s Strong Gravity

In this section, we briefly describe the three destinations within Kerr’s strong gravity that are most relevant to black hole accretion disk theory:

• Destination 3: Event Horizon: That radius inside of which escape from the black hole is not possible;

• Destination 4: Ergosphere: That radius inside of which negative energy states are possible (giving rise to the potentiality of tapping the energy of the black hole).

• Destination 5: Innermost Stable Circular Orbit (ISCO): That radius inside of which free circular orbital motion is not possible;

Our principal question is: Could accretion disk theory unambiguously prove the existence of the event horizon, ergosphere, and ISCO using currently available or future observations?

In realistic astrophysical situations involving astrophysical black holes (in particular quasars and microquasars), the black hole itself is uncharged, and the gravity of accretion disk is practically negligible. This means that the spacetime metric g_μν is given by the Kerr metric, determined by two parameters: total mass M_* and total angular momentum J_*. It is convenient to rescale them by

$$M = {{G{M_\ast}} \over {{c^2}}}$$

(2a)

$$a = {{{J_\ast}} \over {{M_\ast}c}},$$

(2b)

such that both M and a are measured in units of length.

In the standard spherical Boyer-Lindquist coordinates the Kerr metric takes the form [31],

$$\begin{array}{*{20}c} {{g_{tt}} = - \left({1 - {{2Mr} \over {{\varrho ^2}}}} \right),\quad \quad \quad \quad \quad \quad \quad} & {{g^{tt}} = - {{{{({r^2} + {a^2})}^2} - {a^2}\Delta {{\sin}^2}\theta} \over {{\varrho ^2}\Delta}},} \\ {{g_{t\phi}} = - {{2Mar{{\sin}^2}\theta} \over {{\varrho ^2}}},\quad \quad \quad \quad \quad \quad \;} & {{g^{t\phi}} = - {{2Mar} \over {{\varrho ^2}\Delta}},\quad \quad \quad \quad \quad \;\;\;} \\ {{g_{\phi \phi}} = \left({{r^2} + {a^2} + {{2M{a^2}r{{\sin}^2}\theta} \over {{\varrho ^2}}}} \right){{\sin}^2}\theta,} & {{g^{\phi \phi}} = {{\Delta - {a^2}{{\sin}^2}\theta} \over {\Delta {\varrho ^2}{{\sin}^2}\theta}},\quad \quad \quad \quad} \\ {{g_{rr}} = {{{\varrho ^2}} \over \Delta},\quad {g_{\theta \theta}} = {\varrho ^2},\quad \quad \quad \quad \quad \quad \;\;} & {{g^{rr}} = {\Delta \over {{\varrho ^2}}},\quad {g^{\theta \theta}} = {1 \over {{\varrho ^2}}},\quad \quad \quad} \\ \end{array}$$

(3)

where Δ = r² − 2 M r + a² and ϱ² = r² + a² cos² θ.

The Kerr metric depends neither on time t, nor on the azimuthal angle φ around the symmetry axis. These two symmetries can be expressed in a coordinate independent way by the two commuting Killing vectors η^μ = δ^μt and ξ^μ = δ^μφ,

$${\nabla _\mu}{\eta _\lambda} + {\nabla _\lambda}{\eta _\mu} = 0,\quad \quad {\nabla _\mu}{\xi _\lambda} + {\nabla _\lambda}{\xi _\mu} = 0,$$

(4a)

$${\eta ^\mu}{\nabla _\mu}{\xi ^\nu} = {\xi ^\mu}{\nabla _\mu}{\eta ^\nu}.$$

(4b)

Here ∇_μ denotes the covariant derivative,

$${\nabla _\mu}{\eta _\lambda} = {\partial _\mu}{\eta _\lambda} - \Gamma _{\mu \lambda}^\nu {\eta _\nu},\quad \quad \Gamma _{\mu \lambda}^\nu = {1 \over 2}{g^{\nu \kappa}}({\partial _\lambda}{g_{\mu \kappa}} + {\partial _\mu}{g_{\lambda \kappa}} - {\partial _\kappa}{g_{\lambda \mu}}),$$

(5)

and ∂_μ = ∂/∂x^μ denotes the standard partial derivative. Formulae for the Kerr metric (3) and all its non-zero Christoffell symbols $\Gamma _{\mu \lambda}^\nu$ (5), are available from [305].

In Boyer-Lindquist coordinates the t and φ components of the Kerr metric can be expressed as scalar products of the Killing vectors,

$${g_{tt}} = {\eta ^\mu}{\eta _\mu},\quad \quad {g_{t\phi}} = {\eta ^\mu}{\xi _\mu},\quad \quad {g_{\phi \phi}} = {\xi ^\mu}{\xi _\mu}.$$

(6)

From the Killing vectors one can also define the following constants of motion for a particle or photon with four-momentum p_μ

$${\rm{energy}}:\quad \mathcal{E} \equiv - {\eta ^\mu}{p_\mu} = - {p_t},$$

(7a)

$${\rm{angular}}\;{\rm{momentum}}:\quad \mathcal{L} \equiv {\xi ^\mu}{p_\mu} = {p_\phi},$$

(7b)

$${\rm{specific}}\,{\rm{angular}}\,{\rm{momentum}}:\quad \ell \equiv {{\mathcal L} \over {\mathcal E}} = - {{{p_\phi}} \over {{p_t}}} = - {{{u_\phi}} \over {{u_t}}},$$

(7c)

$${\rm{Carter}}\;{\rm{constant}}:\quad \mathcal{C} = {K^{\mu \nu}}{p_\mu}{p_\nu} = {({p_\theta})^2} + {{p_\phi ^2} \over {{{\sin}^2}\theta}}.$$

(7d)

The Carter constant ${\mathcal C}$ is connected to the Killing tensor K_μν, which exists in the Kerr metric. Killing tensors obey,

$${1 \over 6}[{\nabla _\kappa}{K_{\mu \nu}} + {\nabla _\mu}{K_{\nu \kappa}} + {\nabla _\nu}{K_{\kappa \mu}} + {\nabla _\kappa}{K_{\nu \mu}} + {\nabla _\nu}{K_{\mu \kappa}} + {\nabla _\mu}{K_{\kappa \nu}}] = 0.$$

(8)

The other coordinates often used in black hole accretion disk research are the Kerr-Schild coordinates, in which the metric takes the form,

$$\begin{array}{*{20}c} {{g_{\mu \nu}} = {\eta _{\mu \nu}} + f{k_\mu}{k_\nu},\quad \quad {\eta _{\mu \nu}} = {\rm{diag}}(- 1,1,1,1)} \\ {f = {{2M{r^3}} \over {{r^4} + {a^2}{z^2}}},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {{k_t} = 1,\;\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {{k_x} = {{rx + ay} \over {{r^2} + {a^2}}},\quad \quad \quad \quad \quad \;\;\quad \quad \quad \quad \quad \quad} \\ {{k_y} = {{ry - ax} \over {{r^2} + {a^2}}},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \;\;} \\ {{k_z} = {z \over r},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$

(9)

where k_ν = (k_t, k_x, k_y, k_z) is a unit vector, and r is given implicitly by the condition,

$${{{x^2} + {y^2}} \over {{r^2} + {a^2}}} + {{{z^2}} \over {{r^2}}} = 1.$$

(10)

Note that we have given the Kerr-Schild metric in its Cartesian form to prevent confusion with the spherical-polar Boyer-Lindquist coordinates. In keeping with this, unless specifically stated otherwise, the indices {t, φ, r, θ} will always refer to the Boyer-Lindquist coordinates in this review.

2.1 The event horizon

The mathematically precise, general, definition of the event horizon involves topological considerations [207]. Here, we give a definition which is less general, but in the specific case of the Kerr geometry is fully equivalent.

The Boyer-Lindquist coordinates split the Kerr spacetime into a “time” coordinate t and a three-dimensional “space,” defined as t = const hypersurfaces. This split may be done in a coordinate independent way, based on the Killing vectors which exist in the Kerr spacetime. Indeed, the family of non-geodesic observers N^μ with trajectories orthogonal to a family of 3-D spaces t = const is defined as,

$${\rm{ZAMO}}:\;\;{N^\mu} \equiv {e^{- \Phi}}{\tilde \eta ^\mu},\;\;{\tilde \eta ^\mu} = {\eta ^\mu} + \omega {\xi ^\mu},$$

(11a)

$${\rm{frame}}\;{\rm{dragging}}:\;\;\omega \equiv - {{{\eta ^\nu}{\xi _\nu}} \over {{\xi ^\mu}{\xi _\mu}}} = - {{{g_{t\phi}}} \over {{g_{\phi \phi}}}}.$$

(11b)

They are called zero-angular-momentum-observers (ZAMO), because for them, the angular momentum defined by (7b) is zero, ${\mathcal L} = {N_\phi} = 0$. The ZAMO observers provide the standard of rest in the 3-D space: objects motionless with respect to the ZAMO frame of reference occupy fixed positions in space.

We can also define a gravitational potential in the ZAMO frame:

$${\rm{potential}}:\;\;\Phi \equiv - {1 \over 2}\ln \left[ {{{{\xi ^\mu}{\xi _\mu}} \over {{{({\eta ^\nu}{\xi _\nu})}^2} - ({\eta ^\nu}{\eta _\nu})({\xi ^\mu}{\xi _\mu})}}} \right] = - {1 \over 2}\ln \vert{g^{tt}}\vert.$$

(12)

The primary reason to call Φ the gravitational potential is that, in Newton’s theory, the observer who stays still in space experiences an acceleration due to “gravity” g_μ, which equals the gradient of the gravitational potential. In the Kerr spacetime it is,

$${g_\mu} = {({a_\mu})_{{\rm{ZAMO}}}} \equiv {N^\nu}{\nabla _\nu}{N_\mu} = {\nabla _\mu}\Phi.$$

(13)

From (11a) one sees that at the surface 1/g^tt = 0, the vector ῆ^μ is null, ῆ^μῆ_μ = 0. Therefore, the ZAMO observers who provide the standard of rest, move on that surface with the speed of light. In order to stand still in this location, one must move radially out with the speed of light.¹¹ As it is clear from (3), 1/g^tt = 0 is equivalent to Δ = 0. The last equation has a double solution,

$$r = {r_H} = {r_G}\left({1 \pm \sqrt {1 - {a \over M}}} \right).$$

(14)

Note, that for r < r₊ the ZAMO “observers” are spacelike: standing still at a given radial location implies moving along a spacelike trajectory — i.e., faster than light. All trajectories that move radially out are also spacelike. Thus, the outer root r = r₊ of Eq. (14) defines the Kerr black hole event horizon: a null surface that surrounds a region from which nothing may escape. Outside the outer horizon (i.e., for r > r₊) the normalization of N^μ is non-singular, and therefore the gravitational potential (12) is a non-singular, well-defined quantity.

2.1.1 Detecting the event horizon

One may think of two general classes of astrophysical observations that could provide evidence for a black hole horizon. Arguments in the first class are indirect; they are based on estimating a dimensionless “compactness parameter”

$$\mathfrak{C} = {{({\rm{size}}\;{\rm{of}}\;{\rm{the}}\;{\rm{object}})} \over {(G/{c^2})({\rm{mass}}\;{\rm{of}}\;{\rm{the}}\;{\rm{object}})}}.$$

(15)

Arguments in the second class are more direct. They are based (in principle) on showing that some amount of radiation emitted by the source is lost inside the horizon.

Evidence based on estimating the compactness parameter: A source for which observations indicate ℭ ≈ 1 may be suspected of having an event horizon. Values ℭ ≈ 1 have indeed been found in several astronomical sources. In order to know ℭ, one must know mass and size of the source. The mass measurement is usually a direct one, because it may be based on an application of Kepler’s laws. In a few cases the mass measurement is remarkably accurate. For example, in the case of Sgr A*, the supermassive black hole in the center of our Galaxy, the mass is measured to be M = (4.3 ± 0.5) × 10⁶ M_⊙ [111].

Until recently, estimates of size were always indirect, and generally not accurate. They are usually based on time variability or spectral considerations. For the former, the measurement rests on the logic that if the shortest observed variability time-scale is Δt, then the size of the source cannot be larger than R = cΔt. For the latter, the argument goes like this: If the total radiative power L and the radiative flux F can be independently measured for a black-body source, then its size can be estimated from L = 4π R² F. Keep in mind that one must know the distance to the source in order to measure L. The flux can be estimated from F = aT⁴, where T is the temperature corresponding to the peak in the observed intensity versus frequency electromagnetic spectrum.

It is hoped that in the near future, the next generation of high-tech radio telescopes will be able to measure directly the size of “the light circle”, which is uniquely related to the horizon size (see Figure 2). For Sgr A*, at a distance of 8.28 ± 0.44 kpc [111], the event horizon corresponds to an angular size of ∼ 10 μas in the sky, making it an ideal target for near-future microarcsecond very long base interferometric techniques [77, 83]. Here the plan is to measure the black hole “shadow” or “silhouette.”

Figure 2

Silhouettes of Sgr A* calculated for four optically thin accretion structures, characterized by very different physical conditions. The display is intentionally reversed in black-and-white and saturated in order to better show the less luminous parts. Although “dirty astrophysics” makes the most prominent differences, effects of the “pure strong gravity” are also seen in the form of “the light circle”, a tiny almost circular feature at the center. Its shape and size depends only on the black hole mass and spin. Image reproduced by permission from [297], copyright by ESO.

Evidence based on the “no escape” argument: For accretion onto an object with a physical surface (such as a star), 100% of the gravitational binding energy released by accretion must be radiated away. This does not apply for a black hole since the event horizon allows for the energy to be advected into the hole without being radiated. This may allow for a black hole with an event horizon to be distinguished from another, similar-mass object with a surface, such as a neutron star. This argument was first developed by Narayan and collaborators [215, 216, 214]; we describe it in more detail in Section 12.2.

2.2 The ergosphere

In Newtonian gravity, angular momentum ℓ and angular velocity Ω are related by the formula ℓ = r²Ω, and therefore there is no ambiguity in defining a non-rotating frame as Ω = 0 = ℓ. However, in the Kerr geometry ℓ ∝ (Ω − ω), where ω = −g_tφ/g_φφ is the angular velocity of the frame dragging induced by the Lense-Thirring effect. Therefore, Ω = 0 does not imply ℓ = 0. This leads to two different standards of “rotational rest”: the Zero Angular Velocity Observer (ZAVO) and the Zero Angular Momentum Observer (ZAMO),

$${\rm{ZAVO}}\;{\rm{frame}}(\Omega = 0):\;{n^\mu} = {(- {\eta ^\nu}{\eta _\nu})^{- 1/2}}{\eta ^\mu},$$

(16)

$${\rm{ZAMO}}\;{\rm{frame}}(\ell = 0):\;\;{N^\mu} = {e^\Phi}({\eta ^\mu} + \omega {\xi ^\mu}){.}$$

(17)

These two frames rotate with respect to each other with the frame-dragging angular velocity ω = −g_tφ/g_φφ.

The ZAMO frame defines a local standard of rest with respect to the local compass of inertia: a gyroscope stationary in the ZAMO frame does not precess. Considering its kinematic invariants,¹² one sees that the ZAMO frame is non-inertial (a_μ ≠ 0), non-rigid (σ_μν ≠ 0, Θ = 0), and surface-forming (ω_μν = 0). The ZAMO vectors ῆ^μ and N^μ are time-like everywhere outside the horizon, i.e., outside the surface 1/g^tt = 0. This means that the energy of a particle or photon with a four-momentum p_μ measured by the ZAMO is positive, E_ZAMO = N^μp_μ > 0.

The ZAVO frame defines a global standard of rest with respect to distant stars: a telescope that points to a fixed star does not rotate in the ZAVO frame. Considering its kinematic invariants one sees that the ZAVO frame is non-inertial (a_μ ≠ 0), rigid (σ_μν = 0, Θ = 0), and not surface-forming (ω_μν ≠ 0). At infinity, i.e., for r → ∞, it is (η^νη_ν) = g_tt → −1, and therefore n^μ → η^μ. For this reason, η^μ is called the stationary observer at infinity. The ZAVO vectors η^μ and n^μ are timelike outside the region surrounded by the surface g_t = 0, called the ergosphere. Inside the ergosphere η^μ and n^μ are spacelike. This means that inside the ergosphere, the conserved energy of a particle (i.e., the energy measured “at infinity”), as defined by (7a), may be negative.

Penrose [242] considered a process in which, inside the ergosphere, a particle with energy E^∞ > 0 decays into two particles with energies $E_ + ^\infty > 0$ and $E_ - ^\infty = - \vert E_ - ^\infty \vert < 0$. The particle with positive energy escapes to infinity, and the particle with the negative energy gets absorbed by the black hole. Then, because $E_ + ^\infty = E - E_ - ^\infty = E + \vert E_ - ^\infty \vert > E$, one gets a net gain of positive energy at infinity. The source of energy in this Penrose process is the rotational energy of the black hole. Indeed, the angular momentum absorbed by the black hole, J^∞ = p_iξⁱ is necessarily negative (in the sense that J^∞ω_H < 0), which follows from

$${E_{{\rm{ZAMO}}}} \equiv - {p_i}{e^\Phi}({\eta ^i} + {\omega _H}{\xi ^i}) = {e^\Phi}(E_ - ^\infty - {\omega _H}{J^\infty}) > 0,$$

and thus

$${\omega _H}{J^\infty} < E_ - ^\infty < 0.$$

(18)

A more complete presentation of the Penrose process is made in [310]. At this time it appears the most likely realization of the Penrose process would be the Blandford-Znajek mechanism [49] for launching jets from quasars and microquasars. Observations suggest [255, 252], and simulations confirm [303, 304], that through this mechanism it is possible to extract more energy from the system than is being delivered by accretion. We discuss jets and the Blandford-Znajek mechanism more in Sections 10 and 11.7.

2.3 ISCO: the orbit of marginal stability

Particles (with velocity normalization u^μu_μ = −1) and photons (with velocity normalization u^μu_μ = 0) move freely on “geodesic” trajectories x^μ = x^μ(s), with velocities u^μ = dx^μ/ds, characterized by vanishing accelerations

$${a^\nu} \equiv {u^\mu}{\nabla _\mu}{u^\nu} = {{{d^2}{x^\nu}} \over {d{s^2}}} - \Gamma _{\;\;\kappa \mu}^\nu {{d{x^\kappa}} \over {ds}}{{d{x^\mu}} \over {ds}} = 0.$$

(19)

A constant of motion χ, such as those defined in Eq. (7), is conserved along a geodesic trajectory (19) in the sense that u^μΔ_μχ = 0.

Circular geodesic motion in the equatorial plane (θ =π/2) is of fundamental importance in black hole accretion disk theory. The four velocity corresponding to circular motion is defined by,

$${u^\mu} = A({\eta ^\mu} + \Omega {\xi ^\mu}),$$

(20)

where Ω = u^φ/u^t = dφ/dt is the angular velocity measured by the stationary observer (ZAVO, see Section 2.2), and the redshift factor, A = u^t, follows from u^μu^νg_μν = −1,

$$- {A^{- 2}} = {g_{tt}} + 2\Omega {g_{t\phi}} + {\Omega ^2}{g_{\phi \phi}}.$$

(21)

Other connections between these quantities that are particularly useful in our later calculations also follow from u^μu^νg_μν = −1:

$$\Omega = - {{\ell {g_{tt}} + {g_{t\phi}}} \over {\ell {g_{t\phi}} + {g_{\phi \phi}}}},\quad \ell = - {{\Omega {g_{\phi \phi}} + {g_{t\phi}}} \over {\Omega {g_{t\phi}} + {g_{tt}}}},\quad {u_t} = A({g_{tt}} + {\Omega _{gt\phi}}){.}$$

(22)

It is convenient to define the effective potential,

$${\mathcal{U}_{{\rm{eff}}}} = - {1 \over 2}\ln \vert{g^{tt}} - 2\ell {g^{t\phi}} + {\ell ^2}{g^{\phi \phi}}\vert,$$

(23)

because in terms of ${{\mathcal U}_{{\rm{eff}}}}$ and the rescaled energy ${{\mathcal E}^*} = \ln \;{\mathcal E}$, slightly non-circular motion, i.e., with V² = u^ru^rg_rr + u^θu^θg_θθ ≪ u^φu^φg_φφ, is characterized by the equation,

$${1 \over 2}{V^2} = {\mathcal{E}^\ast} - {\mathcal{U}_{{\rm{eff}}}},$$

(24)

which has the same form and the same physical meaning as the corresponding Newtonian equation. Therefore, exactly as in Newtonian theory, unperturbed circular Keplerian orbits are given by the condition of an extremum (minimum or maximum) of the effective potential (θ = π/2),

$$\left\{{{{\left({{{\partial {\mathcal{U}_{{\rm{eff}}}}} \over {\partial r}}} \right)}_\ell} = 0} \right\} \Rightarrow \left\{{\left({{{\partial {g^{tt}}} \over {\partial r}}} \right) - 2\ell \left({{{\partial {g^{t\phi}}} \over {\partial r}}} \right) + {\ell ^2}\left({{{\partial {g^{\phi \phi}}} \over {\partial r}}} \right) = 0} \right\}.$$

(25)

This quadratic equation for ℓ has two roots ℓ = ±ℓ_K (r, a), corresponding to “corotating” and “counterrotating” Keplerian orbits. Their explicit algebraic form is given in Eq. (35) in Section 2.5.

As in Newtonian theory, slightly non-circular orbits (with V ≠ 0 being either δṙ or $\delta \dot \theta$) are fully determined by the simple harmonic oscillator equations,

$$\delta \ddot r + \omega _r^2\delta r = 0,\quad \delta \ddot \theta + \omega _\theta ^2\delta \theta = 0,$$

(26)

where the radial ω_r and vertical ω_θ epicyclic frequencies are second derivatives of the effective potential,

$$\omega _r^2 = {\left({{{{\partial ^2}{\mathcal{U}_{{\rm{eff}}}}} \over {\partial {r_\ast}^2}}} \right)_\ell},\quad \omega _\theta ^2 = {\left({{{{\partial ^2}{\mathcal{U}_{{\rm{eff}}}}} \over {\partial {\theta _\ast}^2}}} \right)_\ell},$$

(27)

where ∂$\partial {x_\ast}^2 = - {g_{xx}}\partial {x^2}$. The epicyclic frequencies (27) are measured by the comoving observer. To get the frequencies ω*_x measured by the stationary “observer at infinity” (Section 2.2), one must rescale by the redshift factor ω*_x = Ax ω_x. Obviously, when (ω_r)² < 0, the epicyclic radial oscillations described by Eq. (26) are unstable — from Eq. (27) we see that they correspond to maxima of the effective potential. This happens for all circular orbits with radii less than r = r_ms (a), and this limiting radius is called ISCO, the innermost stable circular orbit.

Free circular orbits with r > r_ms are stable, while those with r < r_ms are not. Accordingly, accretion flows of almost free matter (i.e., with stresses insignificant in comparison with gravity or centrifugal effects), resemble almost circular motion for r > r_ms, and almost radial free-fall for r < r_ms. For thin disks, this transition in the character of the flow is expected to produce an effective inner truncation radius in the disk (see Section 5.3). The exceptional stability of the inner radius of the X-ray binary LMC X-3 [293], provides considerable evidence for such a connection and, hence, for the existence of the ISCO. The transition of the flow at the ISCO may also show up in the observed variability pattern, if variability is modulated by the orbital motion. In this case, one may expect that the there will be no variability observed with frequencies ν > ν_ISCO, i.e., higher than the Keplerian orbital frequency at ISCO, or that the quality factor for variability, Q ∼ ν/Δν will significantly drop at ν_ISCO. Several variants of this idea have been discussed [33, 34], and some observational evidence to support them has been presented (see Figure 3).

Figure 3

Evidence for the existence of the ISCO from data recorded by the Rossi X-ray Timing Explorer satellite from neutron star binary source 4U 1636–536 [33]. The source shows quasi-periodic oscillations (QPOs) with frequencies in the range 650 Hz < ν < 900 Hz. The sharp drop in the quality factor (bottom panel) seen at ∼ 870 Hz may be attributable to the ISCO [34].

2.4 The Paczyński-Wiita potential

For a non-rotating black hole (a = 0), the Kerr metric reduces to the Schwarzschild solution,

$${\rm{d}}{s^2} = - \left({1 - {{2M} \over r}} \right){\rm{d}}{t^2} + {\left({1 - {{2M} \over r}} \right)^{- 1}}{\rm{d}}{r^2} + {r^2}[{\rm{d}}{\theta ^2} + {\sin ^2}\theta \;{\rm{d}}{\phi ^2}].$$

(28)

Paczyński and Wiita [236] proposed a practical and accurate Newtonian model for a Schwarzschild black hole, based on the gravitational potential,

$${\Phi _{PW}} = - {{GM} \over {r - {r_S}}},\quad \quad {r_S} = {{2GM} \over {{c^2}}}.$$

(29)

The Paczyński-Wiita potential became a very handy tool for studying black hole astrophysics. It has been used in many papers on the subject and still has applicability today. The Schwarzschild and the Paczyński-Wiita expressions for the Keplerian angular momentum and locations of the marginally stable and marginally bound orbits (Section 2.3) are identical. Similar, though less commonly adopted, pseudo-Newtonian potentials have also been found for Kerr (rotating) black holes [23, 276, 213].

2.5 Summary: characteristic radii and frequencies

We end this section with a few formulae for the Kerr geometry that we will use elsewhere in this review.

Keplerian circular orbits exist in the region r > r_ph, with r_ph being the circular photon orbit. Bound orbits exist in the region r > r_mb, with r_mb being the marginally bound orbit, and stable orbits exist for r > r_ms, with r_ms being the marginally stable orbit (also called the ISCO — Section 2.3). The location of these radii, as well as the location of the horizon r_H and ergosphere r₀, are given by the following formulae [31]:

$${\rm{photon}}\;\;{r_{{\rm{ph}}}} = 2{r_G}\left\{{1 + \cos \left[ {{2 \over 3}{{\cos}^{- 1}}({a_\ast})} \right]} \right\},$$

(30)

$${\rm{bound}}\;\;{r_{{\rm{mb}}}} = 2{r_G}\left({1 - {{{a_\ast}} \over 2} + \sqrt {1 - {a_\ast}}} \right),$$

(31)

$${\rm{stable}}\;\;{r_{{\rm{ms}}}} = {r_G}\left\{{3 + {Z_2} - {{[(3 - {Z_1})(3 + {Z_1} + 2{Z_2})]}^{1/2}}} \right\},$$

(32)

$${\rm{horizon}}\;\;{r_{\rm{H}}} = {r_G}\left({1 + \sqrt {1 - a_\ast^2}} \right),$$

(33)

$${\rm{ergosphere}}\;\;{r_0} = {r_G}\left({1 + \sqrt {1 - a_\ast^2{{\cos}^2}\theta}} \right),$$

(34)

where Z₁ = 1 + (1 − a^2*)^1/3[(1 + a_*)^1/3 + (1 − a_*)^1/3], ${Z_2} = {(3a_\ast^2 + Z_1^2)^{1/2}}$, a_* = a/M, and r_G = GM/c² is the gravitational radius.

The Keplerian angular momentum ℓ_k and angular velocity Ω_K, and the angular velocity of frame dragging ω are given by,

$${\ell _K} = {\ell _K}(r,\;a) = {{{M^{1/2}}({r^2} - 2a{M^{1/2}}{r^{1/2}} + {a^2})} \over {{r^{3/2}} - 2M{r^{1/2}} + a{M^{1/2}}}},$$

(35)

$$\Omega _K^2 = {{GM} \over {{{({r^{3/2}} + a{M^{1/2}})}^2}}},$$

(36)

$$\omega = {{2aMr} \over {({a^2} + {r^2}){\varrho ^2} + 2{a^2}Mr{{\sin}^2}\theta}}.$$

(38)

The epicyclic frequencies measured “at infinity” are (here x = r/M),

$${(\omega _r^\ast)^2} = \Omega _K^2\left({1 - 6{x^{- 1}} + 8{a_\ast}{x^{- 3/2}} - 3a_\ast^2{x^{- 2}}} \right),$$

(39)

$${(\omega _\theta ^\ast)^2} = \Omega _K^2\left({1 - 4{a_\ast}{x^{- 3/2}} + 3a_\ast ^2{x^{- 2}}} \right).$$

(40)

Comparing the Keplerian and epicyclic frequencies and the characteristic radii between the Schwarzschild metric and the Paczyński-Wiita potential (Section 2.4), we find for the Schwarzschild metric,

$$\begin{array}{*{20}c} {\Omega _K^2 = {{GM} \over {{r^3}}} = {{(\omega _\theta ^\ast)}^2},\quad \quad {{(\omega _r^\ast)}^2} = \Omega _K^2(1 - 6{x^{- 1}}),} \\ {{\rm{and}}\quad \quad {x_{{\rm{ms}}}} = 6,\quad \quad {x_{{\rm{mb}}}} = 4,\quad \quad {x_{{\rm{ph}}}} = 3,\quad \quad} \\ \end{array}$$

(41)

and for the Paczyński-Wiita potential,

$$\begin{array}{*{20}c} {\Omega _K^2 = {{GM} \over {{r^3}{{(1 - 2{x^{- 1}})}^2}}} = {{(\omega _\theta ^\ast)}^2},\quad \quad {{(\omega _r^\ast)}^2} = \Omega _K^2{{(1 - 6{x^{- 1}})} \over {(1 - 2{x^{- 1}})}},} \\ {{\rm{and}}\quad {x_{{\rm{ms}}}} = 6,\quad {x_{{\rm{mb}}}} = 4.\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$

(42)

3 Matter Description: General Principles

Having provided a detailed description of the key signatures associated with a black hole spacetime, we now move into the mirkier realm of the accretion disk itself. We start from the fundamental conservation laws that govern the behavior of all matter, namely the conservation of rest mass and conservation of energy-momentum, stated mathematically as

$${\nabla _\mu}(\rho {u^\mu}) = 0,\quad \quad {\nabla _\mu}T_\nu ^\mu = 0.$$

(43)

Here ρ is the rest mass density, u^μ is the four velocity of matter, and $T_\nu ^\mu$ is the stress energy tensor describing properties of the matter. The conservation equations (43) are supplemented by numerous “material” equations, like the equation of state, prescriptions of viscosity, opacity, conductivity, etc. Several of them are phenomenological or simple approximations. Nevertheless, we can give a GEN-eral form of $T_\nu ^\mu$ that is relevant to accretion disk theory as a sum of FLU-id, VIS-cous, MAX-well, and RAD-iation parts, which may be written as,

$${(T_\nu ^\mu)_{{\rm{GEN}}}} = {(T_\nu ^\mu)_{{\rm{FLU}}}} + {(T_\nu ^\mu)_{{\rm{VIS}}}}{(T_\nu ^\mu)_{{\rm{MAX}}}} + {(T_\nu ^\mu)_{{\rm{RAD}}}},$$

(44)

$${(T_\nu ^\mu)_{{\rm{FLU}}}} = (\rho {u^\mu})(W{u_\nu}) + \delta _\nu ^\mu P,$$

(45)

$${(T_\nu ^\mu)_{{\rm{VIS}}}} = {\nu _\ast}\sigma _\nu ^\mu,$$

(46)

$${(T_\nu ^\mu)_{{\rm{MAX}}}} = {F^{\mu \alpha}}{F_{\alpha \nu}} - {1 \over 4}\delta _\nu ^\mu {F_{\alpha \beta}}{F^{\alpha \beta}},$$

(47)

$${(T_\nu ^\mu)_{{\rm{RAD}}}} = {4 \over 3}E{u^\mu}{u_\nu} + {u^\mu}{F_\nu} + {u_\nu}{F^\mu}.$$

(48)

Here W = enthalpy, $\delta _\nu ^\mu = {\rm{Kronecker}}$ delta tensor, P = pressure, ν_* = kinematic viscosity, $\sigma _\nu ^\mu = {\rm{shear}}$, F^μν = Faraday electromagnetic field tensor, E = radiation energy density, and F^μ = radiation flux. In the remainder of this section we describe these components one by one, including the most relevant details. Most models of accretion disks are given by steady-state solutions of the conservation equations (43), with particular choices of the form of the stress-energy tensor $T_\nu ^\mu$, and a corresponding choice of the supplementary material equations. For example, thick accretion disk models (Section 4) often assume ${(T_\nu ^\mu)_{{\rm{VIS}}}} = {(T_\nu ^\mu)_{{\rm{MAX}}}} = {(T_\nu ^\mu)_{{\rm{RAD}}}} = 0$, thin disk models (Section 5) assume ${(T_\nu ^\mu)_{{\rm{MAX}}}} = 0$, and most current numerical models (Section 11) assume ${(T_\nu ^\mu)_{{\rm{RAD}}}} = 0$.

3.1 The fluid part

The one absolutely essential piece of the stress-energy tensor for describing accretion disks is the fluid part, ${(T_\nu ^\mu)_{{\rm{FLU}}}} = (\rho {u^\mu})(W{u_\nu}) + \delta _\nu ^\mu P$. The fluid density, enthalpy, and pressure, as well as other fluid characteristics, are linked by the first law of thermodynamics, dU = T dS − P dV, which we write in the form,

$${\rm{d}}\epsilon = W\;{\rm{d}}\rho + nT\;{\rm{d}}S,$$

(49)

where U is the internal energy, T is the temperature, S is the entropy, and ϵ = ρc² +Π is the total energy density, with Π being the internal energy density, and

$$V = {1 \over n},\quad \quad U = {\Pi \over n},\quad \quad W = {{P + \epsilon} \over \rho}.$$

(50)

The equation of state is often assumed to be that of an ideal gas,

$$P = {\mathcal{R} \over \mu}\rho T,$$

(51)

with ${\mathcal R}$ being the gas constant and μ the mean molecular weight.

Sometimes we may wish to consider a two temperature fluid, where the temperature T_i and molecular weight μ_i of the ions are different from those of the electrons (T_e and μ_e). For such a case

$$P = {P_i} + {P_e} = {\mathcal{R} \over {{\mu _i}}}\rho {T_i} + {\mathcal{R} \over {{\mu _e}}}\rho {T_e}.$$

(52)

Two-temperature plasmas are critical in advection-dominated flows (discussed in Section 7). Two-temperature fluids are also important when one considers radiation [281], as it is the ions that are generally heated by dissipative processes in the disk, while it is generally the electrons that radiate. Ions and electrons normally exchange energy via Coulomb collisions. As this process is generally not very efficient, the electrons in the inner parts of accretion flows are usually much cooler than the ions (Coulomb collisions are not able to heat the electrons as fast as they radiate or advect into the black hole). However, there have been suggestions that more efficient processes may couple the ions and electrons [245], such as plasma waves [38] or kinetic instabilities [282]. At this point, this remains an open issue in plasma physics, so it is difficult to know how much heating electrons experience.

3.1.1 Perfect fluid

In the case of a perfect fluid, the whole stress-energy tensor (44) is given by its fluid part (45), and all other parts vanish, i.e., ${(T_\nu ^\mu)_{{\rm{GEN}}}} = {(T_\nu ^\mu)_{{\rm{FLU}}}}$. In this particular case, one can use ${\nabla _\mu}(T_\nu ^\mu {\eta ^\nu}) = {\eta ^\nu}{\nabla _\mu}(T_\nu ^\mu) + {T^{\mu \nu}}({\nabla _\mu}{\eta _\nu}) = 0 + 0 = 0$, and similarly derived ${\nabla _\mu}(T_\nu ^\mu {\xi ^\nu}) = 0$, to prove that

$$\mathcal{B} = - W(u\eta)\; = - W\;{u_t},\quad \quad \mathcal{J} = W(u\xi)\; = W\;{u_\phi},$$

(53)

are constants of motion. We can identify ${\mathcal B}$ as the Bernoulli function and ${\mathcal J}$ as the angular momentum. Their ratio is obviously also a constant of motion,

$$\ell = {\mathcal{J} \over {\mathcal B}} = - {{{u_\phi}} \over {{u_t}}},$$

(54)

identical in form with the specific angular momentum (7c), which is a constant of geodesic motion.

3.2 The stress part

In the stress part ${(T_\nu ^\mu)_{{\rm{VIS}}}} = {\nu _\ast}\sigma _\nu ^\mu$, the shear tensor $\sigma _\nu ^\mu$ is a kinematic invariant (cf. Footnote 12). It is defined as

$${\sigma _{\mu \nu}} \equiv {\left[ {{1 \over 2}({\nabla _\mu}{u_\nu} + {\nabla _\nu}{u_\mu}) - \Theta {g_{\mu \nu}}} \right]_ \bot}$$

(55)

where the symbol ⊥ denotes projection into the instantaneous 3-space perpendicular to u^μ in the sense that (X^μ)⊥ = X^α (δ^μα + u^μ u_α). The other kinematic invariants are vorticity,

$${\omega _{\mu \nu}} \equiv {1 \over 2}{({\nabla _\mu}{u_\nu} - {\nabla _\nu}{u_\mu})_ \bot}$$

(56)

and expansion,

$$\Theta \equiv {1 \over 3}({\nabla _\mu}{u^\mu})$$

(57)

In the standard hydrodynamical description (e.g. [168]), the viscous stress tensor, $S_\nu ^\mu$, is proportional to the shear tensor,

$$S_\nu ^\mu = {\nu _\ast}\rho \;\sigma _\nu ^\mu.$$

(58)

The rate of heat generation by viscous stress in a volume V is then given by

$${Q^ +} = \int {S_\nu ^\mu} {\sigma ^\nu}_\mu \;{\rm{d}}V.$$

(59)

In addition, the rates of viscous angular momentum and energy transport across a surface S, with a unit normal vector N^μ, are

$${\mathcal{J}_S} = \int {S_\nu ^\mu} {\xi ^k}{\nu _\ast}{N_\mu}{\rm{d}}S,\quad \quad {\mathcal{B}_S} = \int {S_\nu ^\mu} {\eta ^\nu}{N_\mu}{\rm{d}}S.$$

(60)

For the case of purely circular motion, where uⁱ = A(ηⁱ + Ωξⁱ), the kinematic invariants are

$$\Theta = 0,\quad \quad {\sigma _{\mu \phi}} = {1 \over 2}{A^3}{\Psi ^2}{\partial _\mu}\Omega, \quad \quad {\omega _{\mu \phi}} = {1 \over 2}{A^3}{\Psi ^{- 2}}{\partial _\mu}\ell,$$

(61)

where ${\Psi ^2} = g_{t\phi}^2 - {g_{tt}}{g_{\phi \phi}}$. It is a general property that (X^μ)⊥ u_μ = 0, and so for purely circular motion, one has,

$${\sigma _{\mu \nu}}{\eta ^\nu} = - \Omega {\sigma _{\mu \nu}}{\xi ^\nu},\quad \quad {\omega _{\mu \nu}}{\eta ^\nu} = - \Omega {\omega _{\mu \nu}}{\xi ^\nu}.$$

(62)

From Eqs. (60) and (62) one deduces that for purely circular motion, the rates of energy and angular momentum transport are related as

$${\mathcal{B}_S} = - {\Omega _S}{\mathcal{J}_S},$$

(63)

where Ω_S is the angular velocity averaged on the surface S. From this, one sees that as angular momentum is transported outward, additional energy is carried inward by the fluid.

3.2.1 The alpha viscosity prescription

As we mentioned in Section 1, the viscosity in astrophysical accretion disks can not come from ordinary molecular viscosity, as this is orders of magnitude too weak to explain observed phenomena. Instead, the source of stresses in the disk is likely turbulence driven by the magneto-rotational instability (MRI, described in Section 8.2). Even so, one can still parametrize the stresses within the disk as an effective viscosity and use the normal machinery of standard hydrodynamics without the complication of magnetohydrodynamics (MHD). This is sometimes desirable as analytic treatments of MHD can be very difficult to work with and full numerical treatments can be costly.

For these reasons, the Shakura-Sunyaev “alpha viscosity” prescription [279] still finds application today. It is an ad hoc assumption based on dimensional arguments. Shakura and Sunyaev realized that if the source of viscosity in accretion disks is turbulence, then the kinematic viscosity coefficient ν_* has the form,

$${\nu _\ast} \approx {l_0}{v_0},$$

(64)

where l₀ is the correlation length of turbulence and v₀ is the mean turbulent speed. Assuming that the velocity of turbulent elements cannot exceed the sound speed, v₀ < c_S, and that their typical size cannot be greater than the disk thickness, l₀ < H, one gets

$${\nu _\ast} = \alpha H{c_S},$$

(65)

where 0 < α < 1 is a dimensionless coefficient, assumed by Shakura and Sunyaev to be a constant.

For thin accretion disks (see Section 5) the viscous stress tensor reduces to an internal torque with the following approximate form [see Eqs. (55) and (58)]

$${\mathcal{T}_{r\phi}} \approx \rho {\nu _\ast}r{{\partial \Omega} \over {\partial r}}.$$

(66)

However, for thin disks, r(∂Ω/∂r) ≈ −Ω and c_S ≈ (P/ρ)^1/2 ≈ ΩH, so Shakura and Sunyaev argued that the torque must have the form Trφ = −αP. A critical question that was left unanswered was what pressure P one should consider: P_gas, P_rad, or P_Tot = P_gas + P_rad? This question has now been answered using numerical simulations [128], so that we now know the appropriate pressure to be P_Tot. Typical values of α estimated from magnetohydrodynamic simulations are close to 0.02 [122], while observations suggest a value closer to 0.1 (see [148] and references therein).

3.3 The Maxwell part

Magnetic fields may play many interesting roles in black hole accretion disks. Large scale magnetic fields threading a disk may exert a torque, thereby extracting angular momentum [48]. Similarly, large scale poloidal magnetic fields threading the inner disk, ergosphere, or black hole, have been shown to be able to carry energy and angular momentum away from the system, and power jets [49]. Weak magnetic fields can tap the differential rotation of the disk itself to amplify and trigger an instability that leads to turbulence, angular momentum transport, and energy dissipation (exactly the processes that are needed for accretion to happen) [26, 27].

In most black hole accretion disks, it is reasonable to assume ideal MHD, whereby the conductivity is infinite, and consequently the magnetic diffusivity is zero. Whenever this is true, magnetic field lines are effectively frozen into the fluid. A corollary to this is that parcels of fluid are restricted to moving along field lines, like “beads” on a wire. In ideal MHD, the Faraday tensor obeys the homogeneous Maxwell’s equation

$${\nabla _\mu}{(^\ast}F_\nu ^\mu) = 0,$$

(67)

where $^*F_\nu ^\mu$ is the dual. If we define a magnetic field 4-vector ${b^u} \equiv {u^v}{\;^*}{F_v^\mu }$, then using ${b^\mu}{u_\mu} = 0$ one can show that

$$^\ast F_\nu ^\mu = {b^\mu}{u_\nu} - {b_\nu}{u^\mu}.$$

(68)

Using this, it is easy to show that the spatial components of (67) give the induction equation

$${\partial _t}(\sqrt {- g} {B^i}) = - {\partial _j}[\sqrt {- g} ({B^i}{v^j} - {B^j}{v^i})],$$

(69)

while the time component gives the divergence-free constraint

$${\partial _i}(\sqrt {- g} {B^i}) = 0,$$

(70)

where Bⁱ = u^tbⁱ − uⁱb^t, and g is the 4-metric determinant.

3.3.1 The magneto-rotational instability (MRI)

We mentioned in Section 3.2 that a hydrodynamic treatment of accretion requires an internal viscous stress tensor of the form ${{\mathcal T}_{r\phi}} < 0$. However, we also pointed out that ordinary molecular viscosity is too weak to provide the necessary level of stress. Another possible source is turbulence. The mean stress from turbulence always has the property that ${{\mathcal T}_{r\phi}} < 0$, and so it can act as an effective viscosity. As we will explain in Section 8.2, weak magnetic fields inside a disk are able to tap the shear energy of its differential rotation to power turbulent fluctuations. This happens through a mechanism known as the magneto-rotational (or “Balbus-Hawley”) instability [26, 118, 27]. Although the non-linear behavior of the MRI and the turbulence it generates is quite complicated, its net effect on the accretion disk can, in principle, be characterized as an effective viscosity, possibly making the treatment much simpler. However, no such complete treatment has been developed at this time.

3.4 The radiation part

Radiation is important in accretion disks as a way to carry excess energy away from the system. In geometrically thin, optically thick (Shakura-Sunyaev) accretion disks (Section 5.3), radiation is highly efficient and nearly all of the heat generated within the disk is radiated locally. Thus, the disk remains relatively cold. In other cases, such as ADAFs (Section 7), radiation is inefficient; such disks often remain geometrically thick and optically thin.

In the optically thin limit, the radiation emissivity f has the following components: bremsstrahlung f_br, synchrotron f_synch, and their Comptonized parts f_br,C and f_synch,C. In the optically thick limit, one often uses the diffusion approximation with the total optical depth τ = τ_abs + τ_es coming from the absorption and electron scattering optical depths. In the two limits, the emissivity is then

$$f = \left\{{\begin{array}{*{20}c} {{f_{{\rm{br}}}} + {f_{{\rm{synch}}}} + {f_{{\rm{br}},{\rm{C}}}} + {f_{{\rm{synch}},{\rm{C}}}}} & {{\rm{optically}}\;{\rm{thin}}\;{\rm{case}}\;(\tau \ll 1)\;,} \\ {{{8\sigma T_{\rm{e}}^4} \over {3H\tau}}\quad \quad \quad \quad \quad \quad \quad \;} & {{\rm{optically}}\;{\rm{thick}}\;{\rm{case}}\;(\tau \gg 1)\;,} \\ \end{array}} \right.$$

(71)

where σ is the Stefan-Boltzmann constant. In the intermediate case one should solve the transfer equation to get reliable results, as has been done in [288, 66]. Often, though, the solution of the grey problem obtained by Hubeny [133] can serve reliably:

$$f = {{4\sigma T_e^4} \over H}{\left[ {{{3\tau} \over 2} + \sqrt 3 + {{4\sigma T_e^4} \over H}{{({f_{{\rm{br}}}} + {f_{{\rm{synch}}}} + {f_{{\rm{br}},{\rm{C}}}} + {f_{{\rm{synch}},{\rm{C}}}})}^{- 1}}} \right]^{- 1}}.$$

(72)

In sophisticated software packages like BHSPEC, color temperature corrections in the optically thick case (the “hardening factor”) are often applied [66].

In the remaining parts of this section we give explicit formulae for the bremsstrahlung and synchrotron emissivities and their Compton enhancements. These sections are taken almost directly from the work of Narayan and Yi [225]. Additional derivations and discussions of these equations in the black hole accretion disk context may be found in [299, 295, 225, 87].

3.4.1 Bremsstrahlung

Thermal bremsstrahlung (or free-free emission) is caused by the inelastic scattering of relativistic thermal electrons off (nonrelativistic) ions and other electrons. The emissivity (emission rate per unit volume) is f_br = f_ei + f_ee. The ion-electron part is given by [225]

$${f_{ei}} = {n_e}\bar n{\sigma _T}c{\alpha _f}{m_e}{c^2} \times \left\{{\begin{array}{*{20}c} {4\;{{\left({{{2{\theta _{\rm{e}}}} \over {{\pi ^3}}}} \right)}^{1/2}}(1 + 1{.}781\theta _e^{1.34})\quad} & {{\theta _e} < 1,} \\ {{{9{\theta _{\rm{e}}}} \over {2\pi}}[\ln (1{.}123{\theta _e} + 0{.}48) + 1{.}5]} & {{\theta _e} \geq 1.} \\ \end{array}} \right.$$

(73)

where n_e is the electron number density, ${\bar n}$ is the ion number density averaged over all species, σ_T = 6.62 × 10⁻²⁵ cm² is the Thomson cross section, α_f = 1/137 is the fine structure constant, θ_e = k_BT_e/m_ec² is the dimensionless electron temperature, and k_B is the Boltzmann constant. The electron-electron part is given by [225]

$${f_{ee}} = n_e^2cr_e^2{m_e}{c^2}{\alpha _f} \times \left\{{\begin{array}{*{20}c} {{{20} \over {9{\pi ^{1/2}}}}(44 - 3{\pi ^2})\theta _e^{3/2}(1 + 1{.}1{\theta _e} + \theta _e^2 - 1.25\theta _e^{5/2})} & {{\theta _e} < 1,} \\ {24{\theta _e}[\ln (1{.}1232{\theta _e}) + 1{.}28]\quad \quad \quad \quad \quad \quad \quad} & {{\theta _e} \geq 1.} \\ \end{array}} \right.$$

(74)

where ${r_e} = {e^2}/{m_e}{c^2}$ is the classical radius of the electron.

3.4.2 Synchrotron

Assuming the accretion environment is threaded by magnetic fields, the hot (relativistic) electrons can also radiate via synchrotron emission. For a relativistic Maxwellian distribution of electrons, the formula is [225]

$$f_{{\rm{synch}}}^ - = {{2\pi} \over {3{c^2}}}k{T_e}{d \over {dr}}\left[ {{{3eB\;\theta _e^2{x_M}} \over {4\pi {m_e}c}}} \right],$$

(75)

where e is the electric charge, B is the equipartition magnetic field strength, and x_M is the solution of the transcendental equation

$$\exp (1{.}8899x_M^{1/3}) = 2{.}49 \times {10^{- 10}}\left({{{4\pi {n_e}r} \over B}} \right){{x_M^{- 7/6} + 0{.}40x_M^{- 17/12} + 0{.}5316x_M^{- 5/3}} \over {\theta _e^3{K_2}(1/{\theta _e})}},$$

(76)

where the radius r must be in physical units and K₂ is the modified Bessel function of the second kind. This expression is valid only for θ_e > 1, but that is sufficient in most applications.

3.4.3 Comptonization

The hot, relativistic electrons can also Compton up-scatter the photons emitted via bremsstrahlung and synchrotron radiation. The formulae for this are [225]

$${f_{{\rm{br}},{\rm{C}}}} = {f_{{\rm{br}}}}\left\{{{\eta _1} - {{{\eta _1}{x_c}} \over {3{\theta _e}}} - {{3{\eta _1}[{3^{- ({\eta _3} + 1)}} - {{(3{\theta _e})}^{- ({\eta _3} + 1)}}]} \over {{\eta _3} + 1}}} \right\}$$

(77)

$${f_{{\rm{synch}},{\rm{C}}}} = {f_{{\rm{synch}}}}[{\eta _1} - {\eta _2}{({x_c}/{\theta _e})^{{\eta _3}}}].$$

(78)

Here η = 1 + η₁ + η₂(x/θ_e)^η3 is the Compton energy enhancement factor, and

$$\begin{array}{*{20}c} {x = {{h\nu} \over {{m_e}{c^2}}},\quad {x_c} = {{h{\nu _c}} \over {{m_e}{c^2}}},} & {{\eta _1} = {{{x_2}({x_1} - 1)} \over {1 - {x_1}{x_2}}},\quad \;\;} \\ {\quad \;\;{x_1} = 1 + 4{\theta _e} + 16\theta _e^2,} & {{\eta _2} = {{- {\eta _1}} \over {{3^{{\eta _3}}}}},\quad \quad \quad \quad} \\ {\quad \;\;{x_2} = 1 - \exp (- {\tau _{es}}),} & {{\eta _3} = - 1 - \ln {x_2}/\ln {x_1},} \\ \end{array}$$

(79)

where h is Planck’s constant and v_c is the critical frequency, below which it is assumed that the emission is completely self-absorbed and above which the emission is assumed to be optically thin.

4 Thick Disks, Polish Doughnuts, & Magnetized Tori

In this section we discuss the simplest analytic model of a black hole accretion disk — the “Polish doughnut.” It is simplest in the sense that it only considers gravity (Section 2), plus a perfect fluid (Section 3.1.1), i.e., the absolute minimal description of accretion. We include magnetized tori in Section 4.2, which allows for ${(T_\nu ^\mu)_{{\rm{MAX}}}} \neq 0$, but otherwise ${(T_\nu ^\mu)_{{\rm{VIS}}}} = {(T_\nu ^\mu)_{{\rm{MAX}}}}{(T_\nu ^\mu)_{{\rm{RAD}}}} = 0$ throughout this section.

4.1 Polish doughnuts

Paczyński and his collaborators developed, in the late 1970s and early 1980s, a very general method of constructing perfect fluid equilibria of matter orbiting around a Kerr black hole [139, 236, 235, 234]. They assumed for the stress energy tensor and four velocity,

$$\begin{array}{*{20}c} {{T^\mu}_\nu = {{({T^\mu}_\nu)}_{{\rm{FLU}}}} = \rho W{u^\mu}{u_\nu} + P{\delta ^\mu}_\nu,} \\ {\quad \quad \quad \quad \quad {u^\mu} = A({\eta ^\mu} + \Omega {\xi ^\mu}),\quad} \\ \end{array}$$

(80)

and derived from ∇_μT^μv = 0 that,

$${\nabla _\nu}\ln A - {{\ell {\nabla _\nu}\Omega} \over {1 - \ell \Omega}} = {1 \over \rho}{\nabla _\nu}P.$$

(81)

In the case of a barytropic fluid P = P(ϵ), the right-hand side of Eq. (81) is the gradient of a scalar function, and thus the left-hand side must also be the gradient of a scalar, which is possible if and only if

$$\ell = \ell (\Omega){.}$$

(82)

This statement is one of several useful integrability conditions, collectively called von Zeipel theorems, found by a number of authors [51, 29, 1, 156].

In real flows, the function ℓ = ℓ(Ω) is determined by dissipative processes that have timescales much longer than the dynamical timescale, and are not yet fully understood. Paczyński realized that instead of deriving ℓ = ℓ(Ω) from unsure assumptions about viscosity that involve a free function fixed ad hoc (e.g., by assuming α(r, θ) = const), one may instead assume the result, i.e., assume ℓ = ℓ(Ω). Assuming ℓ = ℓ(Ω) is not self-consistent, but neither is assuming α(r, θ) = const.

In Boyer-Lindquist coordinates, the equation for the equipressure surfaces, P = (r, θ) = const, may be written as r_P = r_P(θ), with the function r_P(θ) given by

$$- {{d{r_P}} \over {d\theta}} = {{{\partial _\theta}P} \over {{\partial _r}P}} = {{(1 - \ell \Omega){\partial _\theta}\ln A + \ell {\partial _\theta}\Omega} \over {(1 - \ell \Omega){\partial _r}\ln A + \ell {\partial _r}\Omega}}.$$

(83)

Using the expressions for A = A(r, θ), Ω = Ω(r, θ), and ℓ = ℓ(r, θ) from Section 2 (Eqs. 21 and 22), one can integrate Eq. (83) to get the equipressure surfaces. A description of how to do this for both Schwarzschild and Kerr black holes is given in [57]. Figure 4 illustrates the simplest (and important) case of ℓ = ℓ(Ω) = ℓ₀ = const.

Figure 4

In equilibrium, the equipressure surfaces should coincide with the surfaces shown by the solid lines in the right panel. Note the Roche lobe, self-crossing at the cusp. The cusp and the center, both located at the equatorial plane θ = π/2, are circles on which the pressure gradient vanishes. Thus, the (constant) angular momentum of matter equals the Keplerian angular momentum at these two circles, ℓ₀ = ℓ_K(r_cusp) = ℓ_K(r_center), as shown in the upper left panel. In this figure W refers to the effective potential. Image reproduced by permission from [98], copyright by RAS.

Another useful way to think about thick disks is from the relativistic analog of the Newtonian effective potential Φ,

$$\Phi - {\Phi _{{\rm{in}}}} = - \int\nolimits_0^P {{{dP} \over {\rho W}}},$$

(84)

where Φ_in is the potential at the boundary of the thick disk. For constant angular momentum ℓ, the form of the potential reduces to Φ = ln(−u_t). Provided ℓ > ℓ_ms, the potential Φ(r, θ) will have a saddle point Φ_cusp at r = r_cusp, θ = π/2. We can define the parameter ΔΦ = Φ_in − Φ_cusp as the potential barrier (energy gap) at the inner edge of the disk. If ΔΦ < 0, the disk lies entirely within its Roche lobe, whereas if ΔΦ > 0, matter will spill into the black hole even without any loss of angular momentum.

Before leaving the topic of Polish doughnuts, we should point out that, starting with the work of Hawley, Smarr, and Wilson [125], this simple, analytic solution has been the most commonly used starting condition for numerical studies of black hole accretion.

4.2 Magnetized Tori

Komissarov [156] was able to extend the Polish doughnut solution by adding a purely azimuthal magnetic field to create a magnetized torus. This is possible because a magnetic field of this form only enters the equilibrium solution as an additional pressure-like term. For example, the extended form of Eq. (81) is

$${\nabla _\nu}\ln A - {{\ell {\nabla _\nu}\Omega} \over {1 - \ell \Omega}} = {{{\nabla _\nu}P} \over \rho} + {{{\nabla _\nu}({\Psi ^2}{P_{{\rm{mag}}}})} \over {{\Psi ^2}\rho}},$$

(85)

where Ψ² = gt_φg_tφ − g_ttg_φφ and Eq. (84) becomes

$$\Phi - {\Phi _{{\rm{in}}}} = - \int\nolimits_0^P {{{dP} \over {\rho W}}} - \int\nolimits_0^{{{\tilde P}_{{\rm{mag}}}}} {{{d{{\tilde P}_{{\rm{mag}}}}} \over {\Psi \rho W}}},$$

(86)

where ${{\tilde P}_{{\rm{mag}}}} = {\Psi ^2}{P_{{\rm{mag}}}}$. Komissarov [156] gives a procedure for solving the case of a barotropic magnetized torus with constant angular momentum (ℓ = const.).

5 Thin Disks

Most analytic accretion disk models assume a stationary and axially symmetric state of the matter being accreted into the black hole. In such models, all physical quantities depend only on the two spatial coordinates: the “radial” distance from the center r, and the “vertical” distance from the equatorial symmetry plane z. In addition, the most often studied models assume that the disk is not vertically thick. In “thin” disks z/r ≪ 1 everywhere inside the matter distribution, and in “slim” disks (Section 6) z/r ≤1.

In thin and slim disk models, one often uses a vertically integrated form for many physical quantities. For example, instead of density ρ(r, z) one uses the surface density defined as,

$$\Sigma (r) = \int\nolimits_{- H(r)}^{+ H(r)} \rho (r,\;z){\rm{d}}z,$$

(87)

where z = ± H(r) gives the location of the surface of the accretion disk.

5.1 Equations in the Kerr geometry

The general relativistic equations describing the physics of thin disks have been derived independently by several authors [169, 7, 13, 106, 40]. Here we present them in the form used in [268]:

• (i) Mass conservation (continuity):

  $$\dot M = - 2\pi \Sigma {\Delta ^{1/2}}{V \over {\sqrt {1 - {V^2}}}},$$

  (88)

  where V is the gas radial velocity measured by an observer at fixed r who co-rotates with the fluid, and Δ has the same meaning as in Section 2.

• (ii) Radial momentum conservation:

  $${V \over {1 - {V^2}}}{{{\rm{d}}V} \over {{\rm{d}}r}} = {\mathcal{A} \over r} - {1 \over \Sigma}{{{\rm{d}}P} \over {{\rm{d}}r}},$$

  (89)

  where

  $$\mathcal{A} = - {{M\tilde A} \over {{r^3}\Delta \Omega _K^ + \Omega _K^ -}}{{(\Omega - \Omega _K^ +)(\Omega - \Omega _K^ -)} \over {1 - {{\tilde \Omega}^2}{{\tilde R}^2}}},$$

  (90)

  Ã= (r² + a²)² − a²Δsin² θ, Ω = u^φ/u^t is the angular velocity with respect to the stationary observer, $\tilde \Omega = \Omega - \omega$ is the angular velocity with respect to the inertial observer, $\Omega _K^ \pm = \pm {M^{1/2}}/({r^{3/2}} \pm a{M^{1/2}})$ are the angular frequencies of the co-rotating and counter-rotating Keplerian orbits, and $\tilde R = \tilde A/({r^2}{\Delta ^{1/2}})$ is the radius of gyration.

• (iii) Angular momentum conservation:

  $${{\dot M} \over {2\pi}}(\mathcal{L} - {\mathcal{L}_{in}}) = {{{{\tilde A}^{1/2}}{\Delta ^{1/2}}\gamma} \over r}\alpha \Pi,$$

  (91)

  where ${\mathcal L} = {u_\phi}$ is the specific angular momentum, γ is the Lorentz factor, Π = 2HP can be considered to be the vertically integrated pressure, is the standard alpha viscosity (Section 3.2.1), and ${{\mathcal L}_{in}}$ is the specific angular momentum at the horizon, which can not be known a priori. As we explain in the next section, it provides an eigenvalue linked to the unique eigensolution of the set of thin disk differential equations, once they are properly constrained by boundary and regularity conditions.

• (iv) Vertical equilibrium:

  $${\Pi \over {\Sigma {H^2}}} = {{{\mathcal{L}^2} - {a^2}({\mathcal{E}^2} - 1)} \over {2{r^4}}},$$

  (92)

  with ${\mathcal E} = - {u_t}$ being the conserved energy associated with the time symmetry.

• (v) Energy conservation:

  $$- {{\alpha \Pi \tilde A{\gamma ^2}} \over {{r^3}}}{{{\rm{d}}\Omega} \over {{\rm{d}}r}} - {{32} \over 3}{{\sigma {T^4}} \over {\kappa \Sigma}} = - {{\dot M} \over {2\pi r\rho}}{1 \over {{\Gamma _3} - 1}}\left({{{{\rm{d}}P} \over {{\rm{d}}r}} - {\Gamma _1}{P \over \rho}{{{\rm{d}}\rho} \over {{\rm{d}}r}}} \right),$$

  (93)

  where T is the temperature in the equatorial plane, k is the mean (frequency-independent) opacity,

  $$\begin{array}{*{20}c} {\quad \;\;{\Gamma _1} = {\beta ^\ast} + (4 - 3{\beta ^\ast})({\Gamma _3} - 1)\;,\quad} \\ {{\Gamma _3} = 1 + {{(4 - 3{\beta ^\ast})({\gamma _g} - 1)} \over {12(1 - \beta/{\beta _m})({\gamma _g} - 1) + \beta}},} \\ \end{array}$$

  β = P_gas/(P_gas + P_rad + P_mag), β_m = P_gas/(P_gas + P_mag), β* = β(4 − β_m)/3β_m, and γ_g is the ratio of specific heats of the gas.

5.2 The eigenvalue problem

Through a series of algebraic manipulations one can reduce the thin disk equations to a set of two ordinary differential equations for two dependent variables, e.g., the Mach number ${\mathcal M} = - V/{c_S} = - V\Sigma/P$ and the angular momentum ${\mathcal L} = {u_\phi}$. Their structure reveals an important point here,

$${{{\rm{d}}\ln \mathcal{M}} \over {{\rm{d}}\ln r}} = {{{\mathcal{N}_1}(r,\mathcal{M},\mathcal{L})} \over {\mathcal{D}(r,\mathcal{M},\mathcal{L})}}{{{\rm{d}}\ln \mathcal{L}} \over {{\rm{d}}\ln r}} = {{{\mathcal{N}_2}(r,\mathcal{M},\mathcal{L})} \over {\mathcal{D}(r,\mathcal{M},\mathcal{L})}}.$$

(94)

In order for this to yield a non-singular physical solution, the numerators ${{\mathcal N}_1}$ and ${{\mathcal N}_2}$ must vanish at the same radius as the denominator ${\mathcal D}$. The denominator vanishes at the “sonic” radius r_sonic where the Mach number is equal to unity, and the equation ${\mathcal D}(r,{\mathcal M},{\mathcal L}) = 0$ determines its location.

The extra regularity conditions at the sonic point ${{\mathcal N}_i}(r,{\mathcal M},{\mathcal L}) = 0$ are satisfied only for one particular value of the angular momentum at the horizon ${{\mathcal L}_{in}}$, which is the eigenvalue of the problem that should be found. For a given α the location of the sonic point depends on the mass accretion rate. For low mass accretion rates one expects the transonic transition to occur close to the ISCO. Figure 5 shows that this is indeed the case for accretion rates smaller than about 0.4 Ṁ_Edd, independent of α, where we use the authors’ definition of Ṁ_Edd = 16 L_Edd/c². At Ṁ = 0.4 Ṁ_Edd a qualitative change occurs, resembling a “phase transition” from the Shakura-Sunyaev behavior to a very different slim-disk behavior. For higher accretion rates the location of the sonic point significantly departs from the ISCO. For low values of α, the sonic point moves closer to the horizon, down to ∼ 4 M for α = 0.001. For α > 0.2 the sonic point moves outward with increasing accretion rate, reaching values as high as 8 M for α = 0.5 and Ṁ = 100 Ṁ_Edd. This effect was first noticed for small accretion rates by Muchotrzeb [212] and later investigated for a wide range of accretion rates by Abramowicz [10], who explained it in terms of the disk-Bondi dichotomy.

Figure 5

Location of the sonic point as a function of the accretion rate for different values of α, for a non-rotating black hole, a = 0, taking Ṁ_Edd = 16 L_Eddc². The solid curves are for saddle type solutions while the dotted curves present nodal type regimes. Image reproduced by permission from [9], copyright by ESO.

The topology of the sonic point is important, because physically acceptable solutions must be of the saddle or nodal type; the spiral type is forbidden. The topology may be classified by the eigenvalues λ₁, λ₂, λ₃ of the Jacobi matrix,

$$\mathcal{J} = \left[ {\begin{array}{*{20}c} {{{\partial \mathcal{D}} \over {\partial r}}} & {{{\partial \mathcal{D}} \over {\partial \mathcal{M}}}} & {{{\partial \mathcal{D}} \over {\partial \mathcal{L}}}} \\ {{{\partial {\mathcal{N}_1}} \over {\partial r}}} & {{{\partial {\mathcal{N}_1}} \over {\partial \mathcal{M}}}} & {{{\partial {\mathcal{N}_1}} \over {\partial \mathcal{L}}}} \\ {{{\partial {\mathcal{N}_2}} \over {\partial r}}} & {{{\partial {\mathcal{N}_2}} \over {\partial \mathcal{M}}}} & {{{\partial {\mathcal{N}_2}} \over {\partial \mathcal{L}}}} \\ \end{array}} \right].$$

(95)

Because $\det ({\mathcal J}) = 0$, only two eigenvalues λ₁, λ₂ are non-zero, and the quadratic characteristic equation that determines them takes the form,

$$2{\lambda ^2} - 2\lambda \;{\rm{tr}}(\mathcal{J}) - [{\rm{tr}}({\mathcal{J}^2}) - {\rm{t}}{{\rm{r}}^2}(\mathcal{J})] = 0.$$

(96)

The nodal-type solution is given by λ₁ λ₂ > 0 and the saddle type by λ₁ λ₂ < 0, as marked in Figure 5 with the dotted and the solid lines, respectively. For the lowest values of only the saddle-type solutions exist. For moderate values of (0. 1 ≤ α ≤ 0. 4) the topological type of the sonic point changes at least once with increasing accretion rate. For the highest solutions, only nodal-type critical points exist.

5.3 Solutions: Shakura-Sunyaev & Novikov-Thorne

Shakura and Sunyaev [279] noticed that a few physically reasonable extra assumptions reduce the system of thin disk equations (88)–(93) to a set of algebraic equations. Indeed, the continuity and vertical equilibrium equations, (88) and (92), are already algebraic. The radial momentum equation (90) becomes a trivial identity 0 = 0 with the extra assumptions that the radial pressure and velocity gradients vanish, and the rotation is Keplerian, $\Omega = \Omega _k^ +$. The algebraic angular momentum equation (91) only requires that we specify ${{\mathcal L}_{in}}$. The Shakura-Sunyaev model makes the assumption that ${{\mathcal L}_{in}} = {{\mathcal L}_k}({\rm{ISCO}})$. This is equivalent to assuming that the torque vanishes at the ISCO. This is a point of great interest that has been challenged repeatedly [164, 104, 25]. Direct testing of this hypothesis by numerical simulations is discussed in Section 11.4.

The right-hand side of the energy equation (93) represents advective cooling. This is assumed to vanish in the Shakura-Sunyaev model, though we will see that it plays a critical role in slim disks (Section 6) and ADAFs (Section 7). Because the Shakura-Sunyaev model assumes the rotation is Keplerian, $\Omega = \Omega _k^ +$, meaning Ω is a known function of r, the first term on the left-hand side of Eq. (93), which represents viscous heating, is algebraic. The second term, which represents the radiative cooling (in the diffusive approximation) is also algebraic in the Shakura-Sunyaev model.

In addition to being algebraic, these thin-disk equations are also linear in three distinct radial ranges: outer, middle, and inner. Therefore, as Shakura and Sunyaev realized, the model may be given in terms of explicit algebraic (polynomial) formulae. This was an achievement of remarkable consequences — still today the understanding of accretion disk theory is in its major part based on the Shakura-Sunyaev analytic model. The Shakura-Sunyaev paper [279] is one of the most cited in astrophysics today (see Figure 6), illustrating how fundamentally important accretion disk theory is in the field.

Figure 6

The number of citations to the Shakura & Sunyaev paper [279] is still growing exponentially, implying that the field of black hole accretion disk theory still has not reached saturation. Image reproduced from the SAO/NASA Astrophysics Data System, URL (accessed 9 Jan 2013): http://adsabs.harvard.edu/abs/1973A&A....24..337S.

The general relativistic version of the Shakura-Sunyaev disk model was worked out by Novikov and Thorne [229], with important extensions and corrections provided in subsequent papers [237, 265, 241]. Here we reproduce the solution, although with a more general scaling: m = M/M_⊙ and ṁ = Ṁc²/L_Edd.

Outer region: P = P_gas, k = K_ff (free-free opacity)

$$\begin{array}{*{20}c} {\quad \quad \quad \;F = [7 \times {{10}^{26}}\;{\rm{erg}}\;{\rm{c}}{{\rm{m}}^{- 2}}\;{{\rm{s}}^{- 1}}]({m^{- 1}}\dot m)r_\ast ^{- 3}{\mathcal{B}^{- 1}}{\mathcal{C}^{- 1/2}}\mathcal{Q},\quad \quad \quad \quad \quad \quad \quad \quad \quad \;} \\ {\quad \quad \quad \;\,\Sigma = [4 \times {{10}^5}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 2}}]({\alpha ^{- 4/5}}{m^{2/10}}\dot m_{0*}^{7/10})r_\ast ^{- 3/4}{\mathcal{A}^{1/10}}{\mathcal{B}^{- 4/5}}{\mathcal{C}^{1/2}}{\mathcal{D}^{- 17/20}}{\mathcal{E}^{- 1/20}}{\mathcal{Q}^{7/10}},} \\ {\quad \quad \quad \;H = [4 \times {{10}^2}\;{\rm{cm}}]\;({\alpha ^{- 1/10}}{m^{18/20}}{{\dot m}^{3/20}})r_\ast ^{9/8}{\mathcal{A}^{19/20}}{\mathcal{B}^{- 11/10}}{\mathcal{C}^{1/2}}{\mathcal{D}^{- 23/40}}{\mathcal{E}^{- 19/40}}{\mathcal{Q}^{3/20}},} \\ {\quad \quad \quad \,{\rho _0} = [4 \times {{10}^2}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 3}}]({\alpha ^{- 7/10}}{m^{- 7/10}}{{\dot m}^{11/20}})r_\ast ^{- 15/8}{\mathcal{A}^{- 17/20}}{\mathcal{B}^{3/10}}{\mathcal{D}^{- 11/40}}{\mathcal{E}^{17/40}}{\mathcal{Q}^{11/20}},} \\ {\quad \quad \quad \;\,T = [2 \times {{10}^8}\;{\rm{K}}]({\alpha ^{- 1/5}}{m^{- 1/5}}{{\dot m}^{3/10}})r_\ast ^{- 3/4}{\mathcal{A}^{- 1/10}}{\mathcal{B}^{- 1/5}}{\mathcal{D}^{- 3/20}}{\mathcal{E}^{1/20}}{\mathcal{Q}^{3/10}},\quad \quad \quad \quad \;} \\ {\beta/(1 - \beta) = [3]({\alpha ^{- 1/10}}{m^{- 1/10}}{{\dot m}^{- 7/20}})r_\ast ^{3/8}{\mathcal{A}^{- 11/20}}{\mathcal{B}^{9/10}}{\mathcal{D}^{7/40}}{\mathcal{E}^{11/40}}{\mathcal{Q}^{- 7/20}},\quad \quad \quad \quad \quad \quad} \\ {\quad \;{\tau _{f\,f}}/{\tau _{es}} = [2 \times {{10}^{- 3}}]({{\dot m}^{- 1/2}})r_\ast ^{3/4}{\mathcal{A}^{- 1/2}}{\mathcal{B}^{2/5}}{\mathcal{D}^{1/4}}{\mathcal{E}^{1/4}}{\mathcal{Q}^{- 1/2}},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$

(97)

where r_* = rc²/GM.

Middle region: P = P_gas, κ = κ_es (electron-scattering opacity)

$$\begin{array}{*{20}c} {\quad \quad \quad \;F = [7 \times {{10}^{26}}\;{\rm{erg}}\;{\rm{c}}{{\rm{m}}^{- 2}}\;{{\rm{s}}^{- 1}}]({m^{- 1}}\dot m)r_\ast ^{- 3}{\mathcal{B}^{- 1}}{\mathcal{C}^{- 1/2}}\mathcal{Q},\quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \;\,\Sigma = [9 \times {{10}^4}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 2}}]({\alpha ^{- 4/5}}{m^{1/5}}{{\dot m}^{3/5}})r_\ast ^{- 3/5}{\mathcal{B}^{- 4/5}}{\mathcal{C}^{1/2}}{\mathcal{D}^{- 4/5}}{\mathcal{Q}^{3/5}},\quad \quad \;} \\ {\quad \quad \quad \;H = [1 \times {{10}^3}\;{\rm{cm}}]\;({\alpha ^{- 1/10}}{m^{9/10}}{{\dot m}^{1/5}})r_\ast ^{21/20}\mathcal{A}{\mathcal{B}^{- 6/5}}{\mathcal{C}^{1/2}}{\mathcal{D}^{- 3/5}}{\mathcal{E}^{- 1/2}}{\mathcal{Q}^{1/5}},\;\;} \\ {\quad \quad \quad \,{\rho _0} = [4 \times {{10}^1}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 3}}]({\alpha ^{- 7/10}}{m^{- 7/10}}{{\dot m}^{2/5}})r_\ast ^{- 33/20}{\mathcal{A}^{- 1}}{{\mathcal B}^{3/5}}{\mathcal{D}^{- 1/5}}{\mathcal{E}^{1/2}}{\mathcal{Q}^{2/5}},} \\ {\quad \quad \quad \;\,T = [7 \times {{10}^8}\;{\rm{K}}]({\alpha ^{- 1/5}}{m^{- 1/5}}{{\dot m}^{2/5}})r_\ast ^{- 9/10}{\mathcal{B}^{- 2/5}}{\mathcal{D}^{- 1/5}}{\mathcal{Q}^{2/5}},\quad \quad \quad \quad \quad \;} \\ {\beta/(1 - \beta) = [7 \times {{10}^{- 3}}]({\alpha ^{- 1/10}}{m^{- 1/10}}{{\dot m}^{- 4/5}})r_\ast ^{21/20}{\mathcal{A}^{- 1}}{\mathcal{B}^{9/5}}{\mathcal{D}^{2/5}}{\mathcal{E}^{1/2}}{\mathcal{Q}^{- 4/5}},\quad \quad \;} \\ {\quad \;{\tau _{f\,f}}/{\tau _{es}} = [2 \times {{10}^{- 6}}]({{\dot m}^{- 1}})r_\ast ^{3/2}{\mathcal{A}^{- 1}}{\mathcal{B}^2}{\mathcal{D}^{1/2}}{\mathcal{E}^{1/2}}{\mathcal{Q}^{- 1}},\quad \quad \quad \quad \quad \quad \quad \quad \;} \\ \end{array}$$

(98)

Inner region: P = P_rad, κ = κ_es

$$\begin{array}{*{20}c} {\quad \quad \quad \;F = [7 \times {{10}^{26}}\;{\rm{erg}}\;{\rm{c}}{{\rm{m}}^{- 2}}\;{{\rm{s}}^{- 1}}]({m^{- 1}}\dot m)r_\ast ^{- 3}{\mathcal{B}^{- 1}}{\mathcal{C}^{- 1/2}}\mathcal{Q},\quad \quad \quad \quad \quad \;\;} \\ {\quad \quad \quad \;\,\Sigma = [5\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 2}}]({\alpha ^{- 1}}{{\dot m}^{- 1}})r_\ast ^{3/2}{\mathcal{A}^{- 2}}{\mathcal{B}^3}{\mathcal{C}^{1/2}}\mathcal{E}{\mathcal{Q}^{- 1}},\quad \quad \quad \quad \quad \quad \quad} \\ {\quad \quad \quad \,H = [1 \times {{10}^5}\;{\rm{cm}}]\;(\dot{m}){\mathcal{A}^2}{\mathcal{B}^{- 3}}{\mathcal{C}^{1/2}}{\mathcal{D}^{- 1}}{\mathcal{E}^{- 1}}\mathcal{Q},\quad \quad \quad \quad \quad \quad \quad \quad \;\;} \\ {\quad \quad \quad {\rho _0} = [2 \times {{10}^{- 5}}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 3}}]({\alpha ^{- 1}}{m^{- 1}}{{\dot m}^{- 2}})r_\ast ^{3/2}{\mathcal{A}^{- 4}}{\mathcal{B}^6}\mathcal{D}{\mathcal{E}^2}{\mathcal{Q}^{- 2}},\quad \quad \quad \;\;\;} \\ {\quad \quad \quad \;\,T = [5 \times {{10}^7}\;{\rm{K}}]({\alpha ^{- 1/4}}{m^{- 1/4}})r_\ast ^{- 3/8}{\mathcal{A}^{- 1/2}}{\mathcal{B}^{1/2}}{\mathcal{E}^{1/4}},\quad \quad \quad \quad \quad \quad \quad} \\ {\beta/(1 - \beta) = [4 \times {{10}^{- 6}}]({\alpha ^{- 1/4}}{m^{- 1/4}}{{\dot m}^{- 2}})r_\ast ^{21/8}{\mathcal{A}^{- 5/2}}{\mathcal{B}^{9/2}}\mathcal{D}{\mathcal{E}^{5/4}}{\mathcal{Q}^{- 2}},\quad \quad \quad \quad} \\ {\;\,{{({\tau _{f\,f}}{\tau _{es}})}^{1/2}} = [1 \times {{10}^{- 4}}]({\alpha ^{- 17/16}}{m^{- 1/16}}{{\dot m}^{- 2}})r_\ast ^{93/32}{\mathcal{A}^{- 25/8}}{\mathcal{B}^{41/8}}{\mathcal{C}^{1/2}}{\mathcal{D}^{1/2}}{\mathcal{E}^{25/16}}{\mathcal{Q}^{- 2}}}. \\ \end{array}$$

(99)

The radial functions ${\mathcal A}, \ldots, {\mathcal Q}$ that appear in Eqs. (97), (98), (99), are defined in terms of y = (r/M)^1/2 and a_* = a/M as [237]:

$$\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\quad \quad \mathcal{A} = 1 + a_ * ^2{y^{ - 4}} + 2a_ * ^2{y^{ - 6}},\quad \quad \quad }&{\mathcal{B} = 1 + {a_ * }{y^{ - 3}},\quad \quad \quad \quad } \\ {\quad \quad \mathcal{C} = 1 - 3{y^{ - 2}} + 2{a_ * }{y^{ - 3}},\quad \quad \quad \quad }&{\mathcal{D} = 1 - 2{y^{ - 2}} + a_ * ^2{y^{ - 4}},\quad } \\ {\quad \quad \mathcal{E} = 1 + 4a_ * ^2{y^{ - 4}} - 4a_ * ^2{y^{ - 6}} + 3a_ * ^4{y^{ - 8}}}&{{\mathcal{Q}_0} = \frac{{1 + {a_ * }{y^{ - 3}}}}{{y{{(1 - 3{y^{ - 2}} + 2{a_ * }{y^{ - 3}})}^{1/2}}}},} \end{array}} \\ {\mathcal{Q} = {\mathcal{Q}_0}\left[ {y - {y_0} - \frac{3}{2}{a_ * }\ln \left( {\frac{y}{{{y_0}}}} \right) - \frac{{3{{({y_1} - {a_ * })}^2}}}{{{y_1}({y_1} - {y_2})({y_1} - {y_3})}}\ln \left( {\frac{{y - {y_1}}}{{{y_0} - {y_1}}}} \right)} \right]\quad \quad \quad \quad \;} \\ {\quad - {\mathcal{Q}_0}\left[ {\frac{{3{{({y_2} - {a_ * })}^2}}}{{{y_2}({y_2} - {y_1})({y_2} - {y_3})}}\ln \left( {\frac{{y - {y_2}}}{{{y_0} - {y_2}}}} \right) - \frac{{3{{({y_3} - {a_ * })}^2}}}{{{y_3}({y_3} - {y_1})({y_3} - {y_2})}}\ln \left( {\frac{{y - {y_3}}}{{{y_0} - {y_3}}}} \right)} \right]} \end{array}$$

Here y₀ = (r_ms/M)^1/2, and y₁, y₂, and y₃ are the three roots of y³ − 3y + 2a_* = 0; that is

$$\begin{array}{*{20}c} {{y_1} = 2\cos [({{\cos}^{- 1}}{a_\ast} - \pi)/3],} \\ {{y_2} = 2\cos [({{\cos}^{- 1}}{a_\ast} + \pi)/3],} \\ {{y_3} = - 2\cos [({{\cos}^{- 1}}{a_\ast})/3].\quad} \\ \end{array}$$

for the numerical solutions reproduced in Figure 7, the opacities were assumed to be κ_es = 0.34 cm² g⁻¹ and ${\kappa _{ff}} = 6.4 \times {10^{22}}{\rho _{{\rm{cgs}}}}T_{\rm{K}}^{- 7/2}{\rm{c}}{{\rm{m}}^2}{g^{- 1}}$, where ρ_cgs is the density in g cm⁻³ and T_K is the temperature in Kelvin.

Figure 7

The innermost part of the disk. In the Shakura-Sunyaev and Novikov-Thorne models, the locations of the maximum pressure (a.k.a. the center) r_center and the cusp r_cusp, as well as the sonic radius r_sound, are assumed to coincide with the ISCO. Furthermore, the angular momentum is assumed to be strictly Keplerian outside the ISCO and constant inside it. In real flows, r_center ≠ r_cusp ≠ r_sound ≠ ISCO, and angular momentum is super-Keplerian between r_cusp and r_center. Image reproduced by permission from [9], copyright by ESO.

The Shakura-Sunyaev and Novikov-Thorne solutions are only local solutions; this is because they do not take into account the full eigenvalue problem described in Section 5.2. Instead, they make an assumption that the viscous torque goes to zero at the ISCO, which makes the model singular there. For very low accretion rates, this singularity of the model does not influence the electromagnetic spectrum [298], nor several other important astrophysical predictions of the model. However, in those astrophysical applications in which the inner boundary condition is important (e.g., global modes of disk oscillations), the Novikov-Thorne model is inadequate. Figure 7 illustrates a few ways in which the model fails to capture the true physics near the ISCO.

6 Slim Disks

The Shakura-Sunyaev and Novikov-Thorne models of thin disks assume that accretion is radiatively efficient. This assumption means that all the heat generated by viscosity at a given radius is immediately radiated away. In other words, the viscous heating is balanced by the radiative cooling locally and no other cooling mechanism is needed. This assumption can be satisfied as long as the accretion rate is small. At some luminosity (L ≈ 0.3 L_Edd), however, the radial velocity is large, and the disk is thick enough, to trigger another mechanism of cooling: advection. It results from the fact that the viscosity-generated heat does not have sufficient time to transform into photons and leave the disk before being carried inwards by the motion of the gas. The higher the luminosity, the more significant advective cooling is. At the highest luminosities, it becomes comparable to the radiative cooling (see Figure 8), and the standard, thin disk approach can no longer be applied.

Figure 8

The advection factor (ratio of advective to radiative cooling) profiles for ṁ = 0.01, 1.0 and 10.0 (here, ṁ = Ṁc²/16 L_Edd). Profiles for α = 0.01 and 0.1 are presented with solid black and dashed red lines, respectively. The fraction f_adv/(1 + f_adv) of heat generated by viscosity is carried along with the flow. In regions with f_adv < 0 the advected heat is released. Image reproduced by permission from [269].

The problem of accretion with an additional cooling mechanism has to be treated in a different way than radiatively efficient flows. Without the assumptions of radiative efficiency and Keplerian angular momentum, it is no longer possible to find an analytic solution to the system of equations presented in Section 5.1. Instead, one has to solve a two-dimensional system of ordinary differential equations (94) with a critical point — the radius at which the gas velocity exceeds the local speed of sound (the sonic radius). This was first done in the pseudo-Newtonian limit by Abramowicz [8], who forged the term “slim disks”. It has since been done using a fully relativistic treatment by Beloborodov [40]. Recently, Sadowski [268] constructed slim disk solutions for a wide range of parameters applicable to X-ray binaries.

These slim disks are in some sense more physical than thin disks and offer a more general set of solutions, while in the limit of low accretion rates they converge to the standard thin disk solutions. Slim disks are more physical in that they extend down to the black hole horizon, as opposed to thin disks that formally terminate at the ISCO. Slim disks are more general in that they may rotate with an angular momentum profile significantly different than the Keplerian one — the higher the accretion rate, the more significant the departure (see Figure 9). The disk thickness also increases with the accretion rate. For rates close the Eddington limit, the maximal H/R ratio reaches 0.3. Finally, the flux emerging from the slim disk surface is modified by the advection. At high luminosities, a large fraction of the viscosity-generated heat is advected inward and released closer to the black hole or not released at all. As a result, the slope of the radial flux profile changes, and radiation is even emitted from within the ISCO (see Figure 10). Due to the increasing rate of advection, the efficiency of transforming gravitational energy into radiative flux decreases with increasing accretion rate. Despite highly super-Eddington accretion rates, the disk luminosity remains only moderately super-Eddington (see Figure 11). The Eddington luminosity may be exceeded because the geometry of the flow is not spherical and the classical definition of this quantity does not apply — most of the accretion takes place in the equatorial plane while the radiation escapes vertically. Thus, the radiation is not capable of stopping the inflow, though it may cause outflows from the surface.

Figure 9

Profiles of the disk angular momentum (u_φ) for α = 0.01 (left) and α = 0.1 (right panel) for different accretion rates (as a reminder, ṁ = Ṁc²/16 L_Edd), showing the departures from the Keplerian profile. These plots are for a non-rotating black hole. Image reproduced by permission from [270], copyright by ESO.

Figure 10

Flux profiles for different mass accretion rates in the case of a non-rotating black hole and two values of α: 0.01 (black solid), 0.1 (red dashed lines). For each value of α there are five lines corresponding to the following mass accretion rates: ṁ = 0.01, 0.1, 1.0, 2.0 and 10.0 (as a reminder, Ṁ = Ṁc²/16 L_Edd). The black hole mass is 10M_⊙. Image reproduced by permission from [270], copyright by ESO.

Figure 11

Top panel: Luminosity vs accretion rate for three values of black hole spin (a_* = a/M = 0.0, 0.9, 0.999) and two values of α = 0.01 (black) and 0.1 (red line). Bottom panel: efficiency of accretion η = (L/L_Edd)/(⊙/⊙_Edd) (here Ṁ_Edd = 16L_Edd/c²). Image reproduced by permission from [269].

7 Advection-Dominated Accretion Flows (ADAFs)

The ADAF, or advection-dominated accretion flow, solution also involves advective cooling. In fact, it carries it to an extreme — nearly all of the viscously dissipated energy is advected into the black hole rather than radiated. Unlike the slim disk solution, which is usually invoked at high luminosities, the ADAF applies when the luminosity (and generally the mass accretion rate) are low.

Because of their low efficiency, ADAFs are much less luminous than the Shakura-Sunyaev thin disks. The solutions tend to be hot (close to the virial temperature), optically thin, and quasispherical (see Figure 12). Their spectra are non-thermal, appearing as a power-law, often with a strong Compton component. This makes them a good candidate for the Hard state observed in X-ray binaries (discussed in Section 12.3).

Figure 12

Profiles of temperature, optical depth, ratio of scale height to radius, and advection factor (the ratio of advective cooling to turbulent heating) of a hot, one-T ADAF (solid lines). The parameters are M = 10M_⊙, Ṁ = 10⁻⁵ L_Edd/c², α = 0.3, and β = P_gas/(P_gas + P_mag) = 0.9. The outer boundary conditions are R_out = 10³R_S, T = 10⁹ K, and v/c_s = 0.5. Two-T solutions with the same parameters and δ = 0.5 (dashed lines) and 0.01 (dot-dashed lines) are also shown for comparison, where δ is the fraction of the turbulent viscous energy that directly heats the electrons. Image reproduced by permission from [321], copyright by AAS.

ADAFs were formally introduced in the Newtonian limit through a series of papers by Narayan and Yi [223, 224, 225], followed closely by Abramowicz [6, 7] and others [106], although the existence of this solution had been hinted at much earlier [134, 257]. In the same spirit as we gave the equations for the Novikov-Thorne solution in Section 5.3 for thin disks, we report the self-similar ADAF solution found by Narayan and Yi [224]. Again we present the solution with the following scaling: r_* = rc²/GM, m = M/M_⊙ and ṁ = Ṁc²/L_Edd.

$$\begin{array}{*{20}c} {v = [ - 3.00 \times {{10}^{10}}\;{\rm{cm}}\;{{\rm{s}}^{- 1}}]\alpha {c_1}r_\ast ^{- 1/2},\quad \quad \quad \quad \quad \quad \quad} \\ {\Omega = [2.03 \times {{10}^5}\;{{\rm{s}}^{- 1}}]{c_2}{m^{- 1}}r_\ast ^{- 3/2},\quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\,c_S^2 = [9.00 \times {{10}^{20}}\;{\rm{cm}}\;{{\rm{s}}^{- 2}}]{c_3}r_\ast ^{- 1},\quad \quad \quad \quad \quad \quad \quad \quad \quad \;\;} \\ {\;\rho = [1.07 \times {{10}^{- 5}}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 3}}]{\alpha ^{- 1}}c_1^{- 1}c_3^{- 1/2}{m^{- 1}}\dot mr_\ast ^{- 3/2},\quad \quad \quad \;\;} \\ {\;P = [9.67 \times {{10}^{15}}\;{\rm{g}}\;{\rm{c}}{{\rm{m}}^{- 1}}\;{{\rm{s}}^{- 2}}]{\alpha ^{- 1}}c_1^{- 1}c_3^{1/2}{m^{- 1}}\dot mr_\ast ^{- 5/2},\quad \quad \quad} \\ {\;B = [4.93 \times {{10}^8}\;{\rm{G}}]{\alpha ^{- 1/2}}{{(1 - {\beta _m})}^{1/2}}c_1^{- 1/2}c_3^{1/4}{m^{- 1/2}}{{\dot m}^{1/2}}r_\ast ^{- 5/4},} \\ {{q^ +} = [2.94 \times {{10}^{21}}\;{\rm{erg}}\;{\rm{c}}{{\rm{m}}^{- 3}}\;{{\rm{s}}^{- 1}}]\epsilon ^{\prime}c_3^{1/2}{m^{- 2}}\dot mr_\ast ^{- 4},\quad \quad \quad \quad \;\;} \\ {{\tau _{{\rm{es}}}} = [1.75]{\alpha ^{- 1}}c_1^{- 1}\dot mr_\ast ^{- 1/2},\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ \end{array}$$

(100)

where v is the radial infall velocity and q⁺ is the viscous dissipation of energy per unit volume. The constants c₁, c₂, and c₃ are given by

$$\begin{array}{*{20}c} {{c_1} = {{(5 + 2\epsilon ^{\prime})} \over {3{\alpha ^2}}}g(\alpha, \; \epsilon ^{\prime})\quad \quad \quad} \\ {{c_2} = {{\left[ {{{2\epsilon ^{\prime}(5 + 2\epsilon ^{\prime})} \over {9{\alpha ^2}}}g(\alpha, \epsilon^{\prime})} \right]}^{1/2}}} \\ {{c_3} = {{2(5 + 2\epsilon ^{\prime})} \over {9{\alpha ^2}}}g(\alpha, \; \epsilon^{\prime})\;,\quad \quad} \\ \end{array}$$

where

$$\begin{array}{*{20}c} {\quad \quad \;\epsilon ^{\prime} = {1 \over {{f_{{\rm{adv}}}}}}\left({{{5/3 - {\gamma _g}} \over {{\gamma _g} - 1}}} \right)\quad \quad} \\ {g(\alpha, \; \epsilon ^{\prime}) \equiv {{\left[ {1 + {{18{\alpha ^2}} \over {{{(5 + 2\epsilon ^{\prime})}^2}}}} \right]}^{1/2}} - 1,} \\ \end{array}$$

and the parameter f_adv represents the fraction of viscously dissipated energy which is advected. The remaining amount, 1 − f_adv, is radiated locally.

The rapid advection in ADAFs generally has two effects: 1) dissipated orbital energy can not be radiated locally before it is carried inward and 2) the rotation profile is generally no longer Keplerian, although Abramowicz [6] found solutions where the dominant cooling mechanism was advection, even when the angular momentum profile was Keplerian. Fully relativistic solutions of ADAFs have also been found numerically [13, 41]. Further discussion of ADAFs is given in the review article by Narayan and McClintock [219].

8 Stability

Having reviewed some of the main analytic models of accretion disks, it is important now to discuss the issue of stability. Since all analytic models presume steady-state solutions, such models are only useful if the resulting solutions are stable. One reason to suspect accretion disks may not be stable is that the systematic differential rotation that is a signature feature of accretion is a potential source of energy, and therefore, of instability. Another is that some level of instability may be essential in accretion disks as it can provide a pathway to the kind of sustained turbulence anticipated by Shakura and Sunyaev (see Section 3.2.1).

8.1 Hydrodynamic stability

Within ideal hydrodynamics, local linear stability of an axisymmetric rotating flow is guaranteed if the Høiland criterion is satisfied [302]:

$${1 \over {{R^3}}}{{\partial {\ell ^2}} \over {\partial R}} - {1 \over {{C_p}\rho}}\nabla p\cdot\nabla S > 0,$$

(101a)

$${{\partial p} \over {\partial z}}\left({{{\partial {\ell ^2}} \over {\partial R}}{{\partial S} \over {\partial z}} - {{\partial {\ell ^2}} \over {\partial z}}{{\partial S} \over {\partial R}}} \right) < 0,$$

(101b)

where C_p is the specific heat at constant pressure and R is the cylindrical radius (see [275] for the criterion for relativistic stars). This criterion can be easily understood in two limits: For non-rotating equilibria (e.g., a non-rotating star), the criterion reduces to the Schwarzschild criterion (∂S/∂r > 0) that the entropy must not increase toward the interior (for stability against convection). Provided this is true, local fluid elements will simply oscillate under stable buoyancy forces. To see the effects of rotation, we can consider an equilibrium that has constant entropy everywhere. Then the Høiland criterion reduces to the Rayleigh criterion (dℓ/dR > 0): the specific angular momentum must not decrease outward. Physically, if one perturbs a fluid element radially outward, it conserves its own specific angular momentum. If the ambient specific angular momentum decreases outward, then the fluid element will be rotating too fast to stay in its new position, and centrifugal forces will push it further outward. Stability would be a fluid element that oscillates at the local epicyclic frequency.

As it turns out, the Høiland criterion is a huge disappointment for understanding why turbulence might exist in accretion disks. This is because it indicates that accretion disks with rotation profiles that do not differ too much from Keplerian should be strongly stable!

8.1.1 Papaloizou-Pringle Instability (PPI)

The Høiland criterion is only a local stability criterion. Flows can be locally stable, yet have global instabilities. An example of this occurs in the Polish doughnut solution (Section 4). Papaloizou and Pringle [238] showed that this solution is marginally stable with respect to local axisymmetric perturbations yet unstable to low-order nonaxisymmetric modes. As with all global instabilities, the existence of the Papaloizou-Pringle instability (PPI) is sensitive to the assumed boundary conditions [44]. In cases where the disk overflows its potential barrier (Roche lobe) and accretes through pressure-gradient forces across the cusp, the PPI is generally suppressed [117].

8.1.2 Runaway instability

Another instability associated with the Polish doughnut is the runaway instability [5]. If matter is overflowing its Roche lobe and accreting onto the black hole, then one of two evolutionary tracks are possible: (i) As the disk loses material it contracts inside its Roche lobe, slowing the mass transfer and resulting in a stable situation, or (ii) as the black hole mass grows, the cusp moves deeper inside the disk, causing the mass transfer to speed up, leading to the runaway instability. Recent numerical simulations show that, while this instability grows very fast, on timescales of a few orbital periods, over a wide range of disk-to-black hole mass ratios when ℓ = const., i.e., a constant specific angular momentum profile [98], it is strongly suppressed whenever the specific angular momentum of the disk increases with the radial distance as a power law, ℓ ∝ r^p [63]. Even values of p much smaller than the Keplerian limit (p = 1/2) suffice to suppress this particular instability. [This is equivalent to angular velocity profiles, Ω ∝ r^−q, with q > 3/2.]

8.2 Magneto-rotational instability (MRI)

Although it had long been suspected that some sort of MHD instability might provide the necessary turbulent stresses to make accretion work, the nature of this instability remained a mystery until the rediscovery of the magneto-rotational instability by Balbus and Hawley [26, 118, 27]. Originally discovered by Velikhov [309], and generalized by Chandrasekhar [58], in the context of vertically magnetized Couette flow between differentially rotating cylinders, the application of this instability to accretion disks was originally missed.

The instability itself can be understood through a simple mechanical model. Consider two particles of gas connected by a magnetic field line. Arrange the particles such that they are initially located at the same cylindrical distance from the black hole but with some vertical separation. Give one of the particles (say the upper one) a small amount of extra angular momentum, while simultaneously taking away a small amount of angular momentum from the lower one. The upper particle now has too much angular momentum to stay where it is and moves outward to a new radius. The lower particle experiences the opposite behavior and moves to a smaller radius. In the usual case where the angular velocity of the flow drops off with radius, the upper particle will now be orbiting slower than the lower one. Since these two particles are connected by a magnetic field line, the differing orbital speeds mean the field line will get stretched. The additional tension coming from the stretching of the field line provides a torque, which transfers angular momentum from the lower particle to the upper one. This just reinforces the initial perturbation, so the separation grows and angular momentum transfer is enhanced. This is the fundamental nature of the instability.

In more concrete terms, consider a disk threaded with a vertical magnetic field B_z and having an Alfvén speed $\upsilon _A^2 = B_z^2/(4\pi \rho)$. The dispersion relation for perturbations of a fluid quantity δX ∼ exp[i(kz − ωt)] is [27]

$${\omega ^4} - (2k{v_A} + \omega _r^2){\omega ^2} + k{v_A}(k{v_A} + rd{\Omega ^2}/dr) = 0{.}$$

(102)

This equation has an unstable solution (ω² < 0), if and only if, kv_A + rdΩ²/dr < 0. Since

$${{\partial \Omega} \over {\partial r}} < 0,$$

(103)

in accretion disks, the instability criteria can generally be met for weakly magnetized disks. More specifically, the MRI exists for intermediate magnetic field strengths. In terms of the natural length scale of the instability (∼ v_A/Ω), the field strength is constrained at the upper limit by the requirement that the unstable MRI wavelength fit inside the vertical thickness of the disk (v_A/Ω ≲ H). This corresponds to field energy densities that are less than the thermal pressure, i.e., b² < P_gas. At the lower end, dissipative processes set a floor on the relevant length scales, and hence, field strengths.

If the conditions for the instability are met, the fastest-growing mode, which dominates the early evolution, has the form of a “channel flow” involving alternating layers of inward- and outward-moving fluid. The amplitude of this solution grows exponentially until it becomes unstable to three-dimensional “parasitic modes” that feed off the gradients of velocity and magnetic field provided by the channel flow. The flow rapidly reaches a state of magnetohydrodynamic turbulence [118, 119]. This instability can be self-sustaining through a nonlinear dynamo process [52] — nonlinear because the motion that sustains or amplifies the magnetic field is driven by the field itself through the MRI. A more complete description of the linear and non-linear evolution of the MRI is provided in the review article by Balbus and Hawley [27]. A general relativistic linear analysis is presented in [20].

8.3 Thermal and viscous instability

It was realized by Shakura and Sunyaev themselves [280], as well as other authors [176, 287], that the Shakura-Sunyaev solution (Section 5.3) should be thermally and viscously unstable for disks in which radiation pressure dominates (when the opacity is governed by electron scattering). The most general and elegant arguments are presented in the classic paper by Piran [246]. This discovery started a long debate, which continues unresolved to this day. A recent update is provided in [60].

To understand the thermal instability better, we consider a disk cooling through radiative diffusion. The local emergent flux at radius r is given by

$${F^ -} = {{ac{T^4}} \over \tau},$$

(104)

where a is the radiation density constant, T is a measure of the disk interior temperature, and τ is half the total vertical optical depth. If the opacity is dominated by electron scattering, then τ ∼ κ_TΣ/2, where Σ is the surface density of the disk, which is constant on the time scales of the instabilities (very much less than the radial flow time scale). The Thomson opacity κ_T is also constant, being independent of temperature, provided there is already sufficient ionization. Therefore, the optical depth is independent of temperature and the cooling rate per unit area is F⁻ ∝ T⁴.

The dissipation rate per unit area is

$${F^ +}\sim rH{\mathcal{T}_{r\phi}}{{d\Omega} \over {dr}}.$$

(105)

Vertical hydrostatic equilibrium implies that the disk half thickness $H \sim 2P/(\Omega _K^2\Sigma)$, so that Eq. (105) becomes

$${F^ +}\sim {{2rP{\mathcal{T}_{r\phi}}} \over {\Omega _K^2\Sigma}}{{d\Omega} \over {dr}}.$$

(106)

In the radiation pressure dominated inner region, P ≃ aT⁴/3, so that Eq. (106) plus the Shakura-Sunyaev assumption ${{\mathcal T}_{r\phi}} = - \alpha P$ imply that F⁺ ∝ T⁸! Hence a perturbative increase in temperature increases both the local cooling and heating rates, but the heating rate increases much faster, leading to a thermal runaway.

Note, though, that this argument only applies when the viscous stress is proportional to the total pressure P_Tot (α being the proportionality constant). For some time it seemed that a plausible way to avoid this instability was to argue that the stress is proportional instead to the gas pressure P_gas. Recent numerical simulations, though, of the magneto-rotational instability in radiation-pressure dominated disks have shown that the stress is, in fact, proportional to the total pressure [128]. Interestingly, these simulations exhibit no sign of the predicted thermal instability. Most observations also argue against the existence of this instability. In the case of accretion onto black holes, the instability is supposed to set in for luminosities in excess of L > 0.01 L_Edd. However, during outbursts, many stellar-mass black hole sources cross this limit both during their rise to peak luminosity and on their decline to quiescence, showing no dramatic symptoms (although they do undergo state changes, as described in Section 12.3). On the contrary, observations suggest that disks in black hole X-ray binaries are stable up to at least L ≈ 0.5L_Edd [80]. Certainly there is no evidence for the sensational behavior anticipated by some models [172, 300].

9 Oscillations

Even when analytic disk solutions are stable against finite perturbations, it is often the case that these perturbations will, nevertheless, excite oscillatory behavior. Oscillations are a common dynamical response in many fluid (and solid) bodies. Here we briefly explore the nature of oscillations in accretion disks. This topic is particularly relevant to understanding the physical mechanisms that may be behind quasi-periodic oscillations (or QPOs, which are discussed in Section 12.4).

There are a number of local restoring forces available in accretion disks to drive oscillations. Local pressure gradients can drive oscillations via sound waves. Buoyancy forces can act through gravity waves. The Coriolis force can operate through inertial waves. Surface waves can also exist, with the restoring force given by the local effective gravity.

Of particular interest are families of low order modes that may exist in various accretion geometries. Such modes will tend to have the largest amplitudes and produce more easily observed changes than their higher-order counterparts. Here we briefly review a couple relatively simple examples for the purpose of illustration. More details can be found in the references given.

9.1 Dynamical oscillations of thick disks

A complete analysis of the spectrum of modes in thick disks has not yet been done. Some progress has been made by considering the limiting case of a slender torus, where slender here means that the thickness of the torus is small compared to its radial separation from the central mass (i.e., the torus has a small cross-sectional area). In this limit, the complete set of modes have been determined for the case of constant specific angular momentum in a Newtonian gravitational potential [43]. A more general analysis of slender torus modes is given in [45].

Any finite, hydrodynamic flow orbiting a black hole, such as the Polish doughnuts described in Section 4, is susceptible to axisymmetric, incompressible modes corresponding to global oscillations at the radial (σ = ω_r) and vertical (σ = ω_θ) epicyclic frequencies. Other accessible modes are found by solving the relativistic Papaloizou-Pringle equation [4]

$$\begin{array}{*{20}c} {{1 \over {{{(- g)}^{1/2}}}}\left\{{{\partial _i}\left[ {{{(- g)}^{1/2}}{g^{ij}}{f^n}{\partial _j}W} \right]} \right\} - ({m^2}{g^{\phi \phi}} - 2m\omega {g^{t\phi}} + {\omega ^2}{g^{tt}}){f^n}W} \\ {\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad = - {{{{({u^t})}^2}{{(\omega - m\Omega)}^2}} \over {c_{s0}^2}}{f^{n - 1}}W,} \\ \end{array}$$

(107)

where the function f is defined by

$${P \over \rho} = {{{P_0}} \over {{\rho _0}}}f(r,\;\theta)\;,\quad \quad f(r,\;\theta) = 1 - {1 \over {nc_{s0}^2}}\left[ {{{\mathcal{E}_0^2} \over 2}({\mathcal{U}_{{\rm{eff}}}} - {\mathcal{U}_{{\rm{eff}},0}}) + \Phi} \right],$$

(108)

m is the azimuthal wave number, n is the polytropic index (assuming an equation of state of the form P = Kρ^1+1/n), c_s0 is the sound speed at the pressure maximum r₀, and g is the determinant of the metric. The necessary boundary condition comes from requiring that the Lagrangian pressure perturbation vanish at the unperturbed surface (f = 0):

$$\Delta P = (\delta P + {\xi ^\alpha}{\nabla _\alpha}P) = 0,$$

(109)

where ξ^α is the Lagrangian displacement vector. Eqs. (107) and (109) describe a global eigenvalue problem for the modes of Polish doughnuts, with three characteristic frequencies: the radial epicyclic frequency ω_r, the vertical epicyclic frequency ω_θ, and the characteristic frequency of inertial modes κ.

In the slender-torus limit, one can write down analytic expressions for a few of the lowest order modes [45] besides just the radial and vertical epicyclic ones. In terms of local coordinates measured from the equilibrium point,

$$x \equiv g_{rr0}^{1/2}\left({{{r - {r_0}} \over {{r_0}}}} \right)\quad {\rm{and}}\quad y \equiv g_{\theta \theta 0}^{1/2}\left({{{\pi/2 - \theta} \over {{r_0}}}} \right),$$

(110)

an eigenfunction of the form W = Axy, for some constant A, yields two modes with eigenfrequencies [45]

$$\bar \sigma _0^2 = {1 \over 2}\left\{{\omega _r^2 + \omega _\theta ^2 \pm {{[{{(\omega _r^2 + \omega _\theta ^2)}^2} + 4\kappa _0^2\omega _\theta ^2]}^{1/2}}} \right\},$$

(111)

where ${{\bar \sigma}_0} = {\sigma _0}/{\Omega _0}$. The positive square root gives a surface gravity mode that has the appearance of a cross(×)-mode. The negative square root gives a purely incompressible inertial (c-) mode, whose poloidal velocity field represents a circulation around the pressure maximum.

An eigenfunction of the form W = A + Bx² + Cy², with arbitrary constants A, B, C, also has two modes with eigenfrequencies [45]

$$\begin{array}{*{20}c} {\bar \sigma _0^2 = {1 \over {2n}}\left\{{(2n + 1)(\omega _r^2 + \omega _\theta ^2) - (n + 1)\kappa _0^2} \right.\quad \quad \quad \quad \quad \quad \quad \quad \quad} \\ {\left. {\quad \quad \pm {{[{{((2n + 1){{(\omega _\theta ^2 - \omega _r^2)}^2} + (n + 1)\kappa _0^2)}^2} + 4(\omega _r^2 - \kappa _0^2)\omega _\theta ^2]}^{1/2}}} \right\}.} \\ \end{array}$$

(112)

The upper sign results in an acoustic mode with the velocity field of a breathing mode. This mode is comparable to the acoustic mode in the incompressible Newtonian limit for ℓ = const., while in the Keplerian limit, the mode frequency becomes that of a vertical acoustic wave. The lower sign corresponds to a gravity wave, with a velocity field reminiscent of a plus(+)-mode. The poloidal velocity fields of all these lowest order modes are illustrated in Figure 13.

Figure 13

Poloidal velocity fields (δu_x, δu_y) of the lowest order, non-trivial thick disk modes. Image reproduced by permission from [45].

9.2 Diskoseismology: oscillations of thin disks

To analyze oscillations of thin disks, one can express the Eulerian perturbations of all physical quantities through a single function δW ∝ δP/ρ. Since the accretion disk is considered to be stationary and axisymmetric, the angular and time dependence are factored out as δW = W(r, z)e^i(mφ−σt), where the eigenfrequency σ(r, z) = ω − mΩ. Here it is useful to assume that the variation of the modes in the radial direction is much stronger than in the vertical direction, z = r cos θ. The resulting two (separated) differential equations for the functional amplitude W(r, z) = W_r(r)W_y(r, y) are given by [244, 311, 312]

$${{{{\rm{d}}^2}{W_r}} \over {{\rm{d}}{r^2}}} - {1 \over {({\omega ^2} - \omega _r^2)}}\left[ {{{\rm{d}} \over {{\rm{d}}r}}({\omega ^2} - \omega _r^2)} \right]{{{\rm{d}}{W_r}} \over {{\rm{d}}r}} + {{{{({u^t})}^2}{g_{rr}}} \over {{c_s}}}({\omega ^2} - \omega _r^2)\left({1 - {\Psi \over {{{\tilde \omega}^2}}}} \right){W_r} = 0,$$

(113)

$$(1 - {y^2}){{{{\rm{d}}^2}{W_y}} \over {{\rm{d}}{y^2}}} - 2gy{{{\rm{d}}{W_y}} \over {{\rm{d}}y}} + 2g[{\tilde \omega ^2}{y^2} + \Psi (1 - {y^2})]{W_y} = 0.$$

(114)

The radial eigenfunction, W_r, varies rapidly with r, while the vertical eigenfunction, W_y, varies slowly. Here y = (z/H)[γ_g/(γ_g − 1)]¹/² is the re-scaled vertical coordinate, γ_g is the adiabatic index, $\tilde \omega (r) = {\omega _r}/\omega \theta$ is the ratio of the epicyclic frequencies from Section 2.3, Ψ is the eigenvalue of the (WKB) separation function, and g = (P/P_c)/(ρ/ρ_c), where P_c and ρ_c are the midplane values. The radial boundary conditions depend on the type of mode and its capture zone (see below). Oscillations in thin accretion disks, then, are described in terms of Ψ(r, σ), along with the angular, vertical, and radial mode numbers (number of nodes in the corresponding eigenfunction) m, j, and n, respectively. Modes oscillate in the radial range where

$$({\omega ^2} - \omega _r^2)\left({1 - {\Psi \over {{{\tilde \omega}^2}}}} \right) > 0$$

(115)

outside the inner radius, r > r_i.

p-modes are inertial acoustic modes defined by Ψ < ῶ² and are trapped where ω² > ωr², which occurs in two zones. The inner p-modes are trapped between the inner disk edge and the inner “Lindblad” radius, i.e., r_i < r < r₋, where gas is accreted rapidly. The outer p-modes occur between the outer Lindblad radius and the outer edge of the disk, i.e., r₊ < r < r_o. The Lindblad radii, r₋ and r₊, occur where ω = ω_r. The outer p-modes are thought to be more consequential as they produce stronger luminosity modulations [233]. In the corotating frame these modes appear at frequencies slightly higher than the radial epicyclic frequency. Pressure is the main restoring force of p-modes.

g-modes are inertial gravity modes defined by Ψ > ῶ². They are trapped where ${\omega ^2} < \omega _r^2$ in the zone r₋ < r < r₊ given by the radial dependence of ω_r, i.e., g-modes are gravitationally captured in the cavity of the radial epicyclic frequency and are thus the most robust among the thin-disk modes. Since this is the region where the temperature of the disk peaks, g-modes are also expected to be most important observationally [244]. In the corotating frame these modes appear at low frequencies. Gravity is their main restoring force.

c-modes are corrugation modes defined by Ψ = ῶ². They are non-radial (m = 1) and vertically incompressible modes that appear near the inner disk edge and precess slowly around the rotational axis. These modes are controlled by the radial dependence of the vertical epicyclic frequency. In the corotating frame they appear at the highest frequencies.

All modes have frequencies ∝ 1/M. Upon the introduction of a small viscosity (ν ∝ α, α ≪ 1), most of the modes grow on a dynamical timescale t_dyn, such that the disk should become unstable. However, evidence for these modes has so far mostly been lacking in MRI turbulent simulations (see Section 11.6). This leaves their relevance in some doubt.

10 Relativistic Jets

Although the main focus of this review is on black hole accretion disk theory, we note that there has long been a strong observational connection between accreting black holes and relativistic jets across all scales of black hole mass. For supermassive black holes this includes quasars and active galactic nuclei; for stellar-mass black holes this includes microquasars. However, the theoretical understanding of disks and jets has largely proceeded separately and the physical link between the two still remains uncertain. Therefore, we present only a few brief comments on the subject in this review. More complete discussions of the theory of relativistic jets may be found in [193]. A review of their observational connection to black holes is given in [205].

In Section 2.2, we described the “Penrose process,” whereby rotational energy may be extracted from a black hole and carried to an observer at infinity. To briefly recap, Penrose [242, 243] imagined a freely falling particle with energy E^∞ disintegrating into two particles with energies $E_ - ^\infty < 0$ and $E_ + ^\infty > 0$. Then, the particle with negative energy $E_ - ^\infty$ falls into the black hole, and the other one escapes to infinity. Clearly, $E_ + ^\infty > {E^\infty}$, so that there is a net gain of energy.

It was first suggested by Wheeler at a 1970 Vatican conference and soon after by others [183, 96] that such a Penrose process may explain the energetics of superluminal jets commonly seen emerging from quasars and other black hole sources. However, a number of authors [31, 314, 161] showed that for $E_ + ^\infty$ to be greater than E^∞, the disintegration process must convert most of the rest mass energy of the infalling particle to kinetic energy, in the sense that, in the center-of-mass frame, the $E_ - ^\infty$ particle must have velocity v > c/2. The argument of Wald [314] is powerful, short and elegant, so we give it here in extenso.

Let p^ν = mu^ν be the four momentum of a particle with the mass m. We assume that in the ZAMO frame (Section 2.2) the particle has a four velocity of the form, u^ν = Γ(N^ν + VΛ^ν), where Λ^ν is a timelike-unit vector (for simplicity we assume Λ^vη_ν = 0), V is the particle 3-velocity in the ZAMO frame, and Γ⁻² = 1 − V². If the disintegration fragments move in the directions ±Λ^ν (which one may prove is energetically most favorable), then the four velocities of these fragments in the center-of-mass frame are,

$$u_ \pm ^\nu = \gamma ({u^\nu} \pm v{\tau ^\nu}),$$

(116)

$${\tau ^\nu} = \Gamma (V{N^\nu} + {\Lambda ^\nu}){.}$$

(117)

The form of (117) follows from the Lorentz transformation {N^ν,Λ^ν}→ {u^ν,τ^ν}. Multiplying (116) by η_ν gives

$${{E_ \pm ^\infty} \over {{m_ \pm}}} = \gamma \left({{{{E^\infty}} \over m}} \right) \pm \gamma v{\left[ {{{\left({{{{E^\infty}} \over m}} \right)}^2} + {1 \over {{g^{tt}}}}} \right]^{1/2}}.$$

(118)

Since 1/g^tt < 1 and realistic particles have $E/m > 1/\sqrt 3$, the condition $E_ - ^\infty < 0$ necessarily requires ν > 1/2. Such highly relativistic disintegration events are not generally seen in nature. To make matters worse, from the upper limit of (118), Wald deduced that the presence of the black hole limits the energy increase to a maximal factor of $1 + \sqrt 2$. Thus, he concluded [314]: “The Penrose mechanism cannot serve as a useful energy source for astrophysical processes. In no case can one obtain energies which are greater by a significant factor than those which already could be obtained by a similar breakup process without the presence of the black hole.”

Replacing particle disintegration with particle collision does not help, even though the center-of-mass energy of such a collision happening arbitrarily close to the horizon of the maximally rotating Kerr black hole may be arbitrarily large [247, 28]. This is because the Wald limit of $1 + \sqrt 2$ still holds [39]. It would seem that even under idealized conditions, the maximal energy of a particle escaping via the Penrose process is only a modest factor above the total initial energy [39].

Therefore, we consider a general matter distribution, described by an unspecified stress-energy tensor $T_\nu ^\mu$. In this case, the energy flux in the ZAMO frame is Eⁱ = −T^ikN^k, and the energy absorbed by the black hole is

$$E = - \int {{T^i}_k} {n^k}d{N_i} > 0,$$

(119)

where f dN_i is the surface integral over the horizon. The inequality sign follows from the fact that the locally measured energy must be positive. The above integral may by transformed into

$$0 < E = - \int {{e^\Phi}} {T^i}_k({\eta ^k} + \omega {\xi ^k})d{N_i} = e_H^\Phi (E_ - ^\infty - {\omega _H}{J^\infty}){.}$$

(120)

As in the classic Penrose process, the necessary condition for the energy gain is:

$${\omega _H}{J^\infty} < E_ - ^\infty < 0.$$

(121)

Thus, in a way fully analogous to the Penrose process for particles (18), one may say that if the energy at infinity increases because the black hole absorbed negative-at-infinity energy, then the black hole rotation must also slow down by absorbing matter with negative angular momentum.

Blandford and Znajek [49] made the brilliant discovery that an electromagnetic form of the Penrose process may work. In their model, the energy for the jet is extracted from the spin energy of the black hole via a torque provided by magnetic field lines that thread the event horizon or ergosphere. The estimated luminosity of the jet is given by [178] (although see [303] for higher order expressions that apply when a/M ∼ 1)

$${L_{{\rm{BZ}}}} = {1 \over {32}}\omega _F^2B_ \bot ^2r_H^2{(a/M)^2},$$

(122)

where $\omega _F^2 \equiv {\Omega _F}({\Omega _H} - {\Omega _F})/\Omega _H^2$ is a measure of the effect of the angular velocity of the field Ω_F relative to that of the hole ${\Omega _H} \equiv a/(r_H^2 + {a^2})$, B_⊥ is the magnetic field normal to the horizon, and r_H is the radius of the event horizon (14). In this model, the only purpose of the disk is to act as the current sheet which continually provides magnetic field to the black hole. This last point led to one of the main objections to the Blandford-Znajek model: Ghosh and Abramowicz [109] argued on astrophysical grounds that accretion disks simply cannot feed the required fields into the black hole. However, recent work by Rothstein and Lovelace [267] has countered this claim and suggested that indeed the disk can serve this role. There are also more fundamental reservations with the Blandford-Znajek model, some of which are presented in [252, 155]. Such claims and counter-claims were for many years characteristic of the uncertainty in the theory of relativistic jets (see [157] for a discussion). However, direct numerical simulations may be helping to clarify the picture, as we discuss in Section 11.7. Plus, there is now observational evidence suggesting a possible connection between black hole spin and jet power, exactly as predicted by the Blandford-Znajek model [220], although again there are countering claims [93]

11 Numerical Simulations

In simulating accretion disks around black holes, there are a number of challenging issues. First, there is quite a lot of physics involved: relativistic gravity, hydrodynamics, magnetic fields, and radiation being the most fundamental. Then there is the issue that accretion disks are inherently multi-dimensional objects. The computational expense of including extra dimensions in a numerical simulation is not a trivial matter. Simply going from one to two dimensions (still assuming axisymmetry for a disk) increases the computational expense by a very large factor (10² or more). Going to three-dimensions and relaxing all symmetry requirements increases the computational expense yet again by a similar factor. Simulations of this size have only become feasible within the last decade and still only with a subset of the physics one is interested in and usually with a very limited time duration.

Another hindrance in simulating accretion disks is the very large range of scales that can be present. In terms of a grid based code, a disk with a scale height of H/R requiring N_z zones to resolve in the vertical direction at some radius R_in, would require something of the order N_z/(H/R) zones to cover each factor of R_in that is treated in the radial direction. The azimuthal direction in a full three-dimensional simulation would require a comparable number of zones to what is used in the radial direction. Given that a very large calculation by today’s standards is 10⁹–10¹⁰ zones, we can see that treating very thin disks (H/R ≲ 0.01) in a full three-dimensional simulation, even with a very modest number of zones in the vertical direction (N_z ≥ 50) would take all the resources one could muster. However, two-dimensional (axisymmetric) simulations of thin disks and three-dimensional simulations of thick, slim, or moderately thin disks have now become quite common.

The approach of numericists in many ways parallels that of theorists, so we structure this section much the same as the first part of this Living Review. That is, after a brief introduction to numerical techniques, we discuss how the various components of the matter description (à la Section 3) are implemented in numerical simulations. We then review a few special cases illustrating how analytic models (Sections 4–7) and numerical simulations can complement one another. We finish this numerical section with a few topics of special interest.

11.1 Numerical techniques

11.1.1 Computational fluid dynamics codes

There are many numerical codes available today that include relativistic hydrodynamics or MHD that are, or can be, used to simulate accretion disks. A partial list includes: Cosmos++ [18], ECHO [74], HARM [105], and RAISHIN [208] plus the codes developed by Koide et al. [152], De Villiers and Hawley [71], Komissarov [154], Antón et al. [19], and Anderson et al. [16]. This represents tremendous growth in the field as prior to 1998 there were only a couple such codes, largely derived from the early work of Wilson [318] and later developments by Hawley [126, 125]. Those early codes were only available to a handful of researchers. This was simply a reflection of the fact that the computational resources available to most researchers at that time were insufficient to make use of such codes, so there was little incentive to develop or acquire them. This has clearly changed.

Today there are primarily two types of numerical schemes being used to simulate relativistic accretion disks, differentiated by how they treat discontinuities (shocks) that might arise in the flow: the artificial viscosity scheme, still largely based upon formulations developed by Wilson [318]; and Godunov-type approaches, using exact or approximate Riemann solvers. Both are based on finite difference representations of the equations of general relativistic hydrodynamics, although the latter tend to also incorporate finite-volume representations. The artificial viscosity schemes have the advantages that they are more straightforward to implement and easier to extend to include additional physics. More importantly, they are computationally less expensive to run than Godunov schemes. Furthermore, recent work by Anninos and collaborators [17, 18] has shown that variants of this scheme can be made as accurate as Godunov schemes even for ultrarelativistic flows, which historically was one of the principal weaknesses of the artificial viscosity approach. Godunov schemes, on the other hand, appreciate the advantage that they are fully conservative, and therefore, potentially more accurate. They also require less tuning since there are no artificial viscosity parameters that need to be set for each problem. More thorough reviews of these two methods, with clear emphasis on the High-Resolution Shock-Capturing variant of the Godunov approach, are provided in the Living Reviews by Martí & Müller [182] and Font [97]. Other numerical schemes, such as smooth-particle hydrodynamics (SPH) and (pseudo-)spectral methods are less well developed for work on relativistic accretion disks.

11.1.2 Global vs. shearing-box simulations

Along with settling on a numerical scheme, a decision must also be made whether or not to try to treat the disk as a whole or to try to understand it in parts. The latter choice includes “shearing-box” simulations, the name coming from the type of boundary conditions one imposes on the domain to mimic the shear that would be present in a real disk [121]. The obvious advantage of treating the disk in parts is that you circumvent the previously noted problem of the large range in scales in the disk by simply ignoring the large scales. Instead one treats a rectangular volume generally no larger than a few vertical scale heights on a side and in some cases much smaller. In this way, for a moderate number of computational zones, one can get as good, or often much better, resolution over the region being simulated than can be achieved in global simulations. The obvious disadvantage is that all information about what is happening on larger scales is lost. By construction, no structures larger than the box can be captured and the periodic boundaries impose some artificial conditions on the simulations, such as no background gradients in pressure or density and no flux of material into or out of the box, meaning the surface density remains fixed. Perhaps most important in the context of a review on relativity is that the box is treated as a local patch within the disk, with a scale smaller than that of the curvature of the metric. Thus most general relativistic effects are not, and by construction need not be, included in shearing box simulations. Since this is a Living Review in Relativity, we will not dwell much on this class of simulations. Still, there are some significant highlights that should be mentioned.

By far the most extensive work using shearing box simulations has been aimed at better understanding the magneto-rotational instability (MRI). Starting from the earliest “proof-of-concept” simulations [121, 296], shearing box simulations have been used to demonstrate various properties of the saturated state of MRI turbulence [271, 272, 173, 283] and much about the vertical structure of MRI turbulent disks [46]. Shearing box simulations have also proven valuable in studying radiation-dominated disks. For example, Turner [307] studied the vertical structure of radiation dominated accretion disks and Hirose et al. [129] studied the gas-pressure dominated case, both including radiation using flux-limited diffusion. Hirose and collaborators subsequently used shearing box simulations to demonstrate that radiation-dominated disks are thermally stable [128], though they had long been held not to be [176, 280, 246]. As we said, though, there are limits imposed by the finite size of the shearing box. For instance, the viscous, Lightman-Eardley instability [176] can only be studied through global simulations.

11.2 Matter description in simulations

The minimum physics required for a global simulation of an accretion disk are gravity and hydrodynamics (assuming the disk is dense enough for the continuum approximation to hold). Since many disks have masses that are small compared to the mass of the central compact object, the self-gravity of the disk can often be ignored. Therefore, in the next three sections, gravity will simply mean that of the central black hole.

11.2.1 Hydrodynamics + gravity

The first researcher to develop and use numerical algorithms for simulating relativistic accretion flows was Wilson [316], who considered the spherical infall of material with a non-zero specific angular momentum toward a Kerr black hole using the full metric, although restricted to two spatial dimensions. Wilson was able to confirm the additional centrifugal support that the infalling material experienced due to the rotating black hole. For sufficiently high values of angular momentum, the material could not immediately accrete onto the black hole, instead forming a fat disk of the type described in Section 4.

In many ways, Wilson’s pioneering work was at least a decade ahead of its time and there was consequently a lull in activity until Hawley and Smarr collaborated with Wilson to revive his work [126, 125]. As an example of how analytic and numerical work can complement each other, it is worth noting that one of the test cases they used for their new two-dimensional relativistic hydrodynamics code was based on the analytic theory for relativistic thick disks, which had been worked out in the time since Wilson’s original simulations. Using the analytic theory, they constructed a series of disks with different (constant) specific angular momenta, from ℓ < ℓ_ms to ℓ > ℓ_mb. The ℓ > ℓ_mb case yielded a static solution as expected and confirmed the ability of their code to accurately evolve such a solution in multi-dimensions over a dynamical time. The ℓ < ℓ_mb cases showed greater time variability and illustrated the power of direct numerical simulations to extend our understanding of black hole accretion.

11.2.2 Magnetohydrodynamics + gravity

Magnetic fields can play many important roles in relativistic accretion disks, from providing local viscous stresses through turbulence that results from the magnetorotational instability (Section 8.2), to providing a mechanism for launching and confining jets (Section 10). Thus, the inclusion of magnetic fields in numerical simulations of relativistic accretion disks is important. Although the techniques for doing this were first described in [317], relatively little work was done in this area until fairly recently. Perhaps the two most important contributions so far from relativistic MHD simulations of accretion disks have been: 1) elucidating the behavior of MRI turbulent disks in the vicinity of black holes, and 2) exploring the many interrelations between magnetic fields and relativistic jets.

In at least one case, the inclusion of MHD + gravity is sufficient to adequately capture the dynamics of a real black hole accretion disk. Sgr A*, the black hole at the center of the Milky Way galaxy has such an anemically low luminosity [194] that numerical simulations that ignore radiation can safely be applied [76], as has been done by many authors [112, 227, 211, 75]. In most other cases, radiative processes must somehow be accounted for.

11.2.3 Radiation-Magnetohydrodynamics + Gravity

Probably the most glaring shortcoming of almost all numerical simulations of accretion disks and many other phenomena in astrophysics to date is the unrealistic treatment of radiation, which is most often simply ignored. This is not due to a lack of appreciation of its importance on the part of numericists, but simply a reflection of the fact that there are very few efficient ways to treat radiation computationally in multiple dimensions.

The optically-thin limit is an exception. In this case, the radiative cooling only enters the energy and momentum conservation equations as a source term

$${\nabla _\mu}{T^{\mu \nu}} = f{u^\nu},$$

(123)

where f is the cooling function of the gas (Section 3.4). The subsequent propagation of the radiation through the matter can be ignored or post-processed. Even so, including radiative cooling can be critically important to capturing the true dynamics of the disk [101, 76], as shown in Figure 14.

Figure 14

Pseudo-color plots of log(T) with contours of logρ from four different general relativistic MHD simulations. The simulations all begin with the same initial conditions, but have different energy conservation and cooling treatments: The upper-left panel conserves internal energy and ignores cooling; the upper-right panel conserves internal energy and includes cooling; the lower-left panel conserves total energy and ignores cooling; the lower-right panel conserves total energy and includes cooling. The very different end states illustrate the importance of properly capturing thermodynamic processes. Image reproduced by permission from [101], copyright by AAS.

In the optically-thick limit, some progress can be made by employing the flux-limited diffusion approximation [174]. Global simulations of disks including a flux-limited diffusion treatment of radiation (but using pseudo-Newtonian gravity) are now being done by Ohsuga and collaborators [232, 231]. A few steps toward the goal of relativistic radiation MHD simulations of black hole accretion disks have also been taken in recent years. A method for treating optically thick accretion using a conservative, Godunov scheme was developed by Farris and collaborators [91]. The same basic method has now been used to examine both Bondi [100] and Bondi-Hoyle [326, 266] accretion. Simulations of accretion disks, though, must await the generalization of this method to treat radiation both in the optically thick and thin limits.

11.2.4 Evolving GRMHD

In most simulations of accretion disks around black holes, the self-gravity of the disk is ignored. In many cases this is justified as the mass of the disk is often much smaller than the mass of the black hole. This is also much simpler as it allows one to treat gravity as a background condition, either through a Newtonian potential or a relativistic metric (the so-called Cowling approximation in relativity). However, there are plausible astrophysical scenarios in which this approximation is not valid. Two of the more interesting are: 1) a tidally disrupted neutron star accreting onto a stellar-mass black hole; and 2) an overlying stellar envelope accreting onto a nascent black hole during the final dying moments of a massive star. Interestingly these scenarios are currently the most popular models of gamma-ray bursts [82, 319], which are the most powerful explosions in the Universe since the Big Bang and are observed at a rate of about one per day.

Treating a self-gravitating disk in general relativity requires a code that can simultaneously evolve the spacetime metric (by solving the Einstein field equations) and the matter, magnetic, and radiation fields. Codes capable of doing this have very recently become available, divided into the following “flavors”: self-gravity + hydro [24]; self-gravity + MHD [81, 286, 110, 88]; and self-gravity + radiation MHD [91]. More detailed discussion of these methods and their application to problems in astrophysics is provided in the Living Review by Font [97].

The importance of evolving the gravitational field in the case of black hole accretion is well illustrated in the case of the runaway instability, which we described in Section 8.1.2. That instability has now been simulated directly, including evolving the spacetime metric [209, 210, 160]. These studies have confirmed earlier results that the runaway instability does not develop.

Another application of these types of codes that is relevant to our review is the reexamination of the Papaloizou-Pringle instability (discussed in Section 8.1.1). Through this instability, a massive, self-gravitating Polish doughnut orbiting around a black hole could become a strong emitter of large amplitude, quasiperiodic gravitational waves [149].

11.3 Polish doughnuts (thick) disks in simulations

As we already mentioned, it is computationally less expensive to simulate thick disks than thin, so it comes as no surprise that most of the numerical work for nearly four decades, starting with the work of Wilson [316], has been on thick disks. In this section we touch on a few highlights from this area of research

Some of this work has focused on confirming predictions from analytic theory. Igumenshchev and Beloborodov [137] used two-dimensional relativistic hydrodynamic simulations to confirm that: (i) The structure of the innermost disk region strongly depends on the black hole spin, and (ii) the mass accretion rate is proportional to the size of the energy gap between the inner edge of the disk and the cusp (see Figure 15), as predicted by Kozlowsk and collaborators [162].

Figure 15

Time-average mass accretion rate ṁ = Ṁc² /L_Edd as a function of the energy gap ΔΦ for models with a = 0 (circles), a/M = 0.9 (squares), and a/M = −0.9 (triangles). The bars show the variability of ṁ. The lines represent the predicted dependencies $\dot m \propto {(\Delta \Phi)^{{\gamma _g}/({\gamma _g} - 1)}}$, where γ_g = 4/3 is the adiabatic index. Image reproduced by permission from [137], copyright by RAS.

Hawley [116, 117] studied the nonlinear evolution of the non-axisymmetric Papaloizou-Pringle instability (Section 8.1.1). He discovered that for radially slender tori, the principal mode of the PPI saturates with the formation of ellipsoidal overdensity regions he termed “planets” (see Figure 16 for an example). In cases where several wavenumbers were unstable, multiple planets would begin to form but then nonlinear mode-mode coupling would cause them to merge. For wide tori, particularly when the tori extended beyond their Roche limits, overflowed their cusps, and accreted into the hole, the PPI saturated at low amplitude or was prevented from growing altogether. The initial work only considered Schwarzschild black holes, but all of the results were later qualitatively confirmed for Kerr black holes [69].

Figure 16

Equatorial slice through hydrodynamic tori at saturation of the Papaloizou-Pringle instability showing formation of significant non-axisymmetric (m = 1) overdensity clumps. The density contours are linearly spaced between ρ_max and 0.0. This figure represents models A3p (left) and B3r (right) of [69]. Image reproduced by permission, copyright by AAS.

Another interesting note is that most of the global numerical simulations of MRI turbulent disks have used a Polish doughnut as the starting condition. This was first done by De Villiers and Hawley [70] (animations can be viewed at [68, 115]) and later by a number of other groups [105, 18]. As expected from linear analysis and earlier Newtonian simulations, the Maxwell stresses produced by the MRI-driven turbulence naturally force the final distribution of the angular momentum to be nearly Keplerian (see Figure 17), independent of the initial matter distribution.

Figure 17

Specific angular momentum ℓ as a function of radius at t = 0 (thin line) and at t = 10.0 orbits (thick line). The individual plots are labeled by model. In each case the Keplerian distribution for a test particle, ℓ_Kep, is shown as a dashed line. Image reproduced by permission from [72], copyright by AAS.

Although the magnetic field plays a crucial role through the action of the MRI, its amplitude saturates at a relatively weak level, always remaining subthermal in the disk (see Figure 18). On the other hand, simulations show that magnetic fields do tend to dominate regions outside the disk, including in the hot, magnetized corona that sandwiches the disk and in the evacuated, highly magnetized funnel region. The funnel region and the boundary between the funnel and the corona are where jets and outflows are generically observed in these simulations (Section 11.7).

Figure 18

Color contours of the ratio of azimuthally averaged magnetic to gas pressure, P_mag/P_gas. The scale is logarithmic and is the same for all panels; the color maps saturate in the axial funnel. The body of the accretion disk is identified with overlaid density contours at 10⁻², 10^−1.5, 10⁻¹, and 10^−0.5 of ρ_max(t = 0). The individual plots are labeled by model. In all cases, the magnetic pressure is low (P_mag/P_gas ≪ 1) in the disk, comparable to gas pressure (P_mag/P_gas ∼ 1) in the corona above and below the disk, and high (P_mag/P_gas ≫ 1) in the funnel region. Image reproduced by permission from [72], copyright by AAS.

Despite the increased complexity of the disk when MRI-driven turbulence is at play, many features and properties remain strikingly similar to what is revealed in hydrodynamic-only studies. For instance, a commonly seen feature in global, non-radiative MHD simulations of black hole accretion disks is an “inner torus” [120]. The inner torus is predominantly a hydrodynamic structure, as it is largely gas pressure supported. Generically, this inner torus consists of a sequence of equidensity and equipressure surfaces with a pressure maximum at ∼ 10r_G, a cusp at ∼ 3r_G, and a roughly parabolic, evacuated, though magnetized, funnel along the rotation axis. This structure is remarkably similar to the Polish doughnut model, as shown in Figure 19. Thus, many hydrodynamic torus results may retain relevance even in MRI turbulent disks [253].

Figure 19

On the left, equidensity contours calculated from an analytic Polish doughnut. On the right, equidensity contours from a numerical MHD simulation (model 90h from [99]). Note, though, that the contours on the left are linearly spaced, while those on the right are logarithmically spaced. Thus, the gradients represented on the left are shallower than those on the right. Image reproduced by permission from [253], copyright by ESO.

11.4 Novikov-Thorne (thin) disks in simulations

Early numerical simulations of black hole accretion disks nearly all focused on thick disks. This was mostly for computational convenience as thicker disks require fewer resources than thin ones. In recent years, though, the resources have become available to start testing thinner disks. This has enabled researchers to begin rigorously testing the Novikov-Thorne model (Section 5.3), especially the assumption that the internal torque of the disk vanishes at the ISCO. This is important as some researchers have argued that magnetic fields might nullify this hypothesis by maintaining stresses inside the ISCO [164, 104, 25].

Initial global MHD simulations did show some variations from the Novikov-Thorne model at the ISCO at the level of ∼ 10% [123, 226], while others [240] showed much smaller deviations ≲ 3%. Penna and collaborators [240] discussed possible reasons for the discrepancies, and suggested that the results are, in fact, in better agreement than they may have initially appeared. Their conclusion is that, in fact, the Novikov-Thorne model matches numerical results quite well (see Figure 20, plus animations can be viewed at [239]). Any remaining differences are small enough to have negligible effect on common applications of the Novikov-Thorne model, such as measuring black hole spin via spectral fitting [166, 327] (Section 12.1).

Figure 20

Left: Time-averaged rest mass density in the r−z plane for four GRMHD simulations with a = 0 and various disk thicknesses. The dashed vertical line marks the ISCO. The disk opening angle, h = H/r, and effective Shakura-Sunyaev viscosity, α, are reported in each panel. The top three panels have h ≪ α and the inner edge of the disk is located outside the ISCO. The bottom panel has h ≪ α and the density increases monotonically down to the event horizon. Figure from [241]. Right: Various fluxes as functions of radius for a numerical Novikov-Thorne disk simulation. Top: Mass accretion rate. Second panel: Accreted specific angular momentum. Solid line is simulation data; dashed line gives Novikov-Thorne solution; dotted line is ISCO value. Note that the specific angular momentum does not drop significantly inside the ISCO, unlike for thick disks, such as in Figure 17. Third panel: The “nominal” efficiency, which is the total loss of specific energy from the fluid. Bottom panel: Specific magnetic flux. The near constancy of this quantity inside the ISCO is an indication that magnetic stresses are not significant in this region. Image reproduced by permission from [240].

We mention that since present day relativistic numerical simulations do not treat the radiation in geometrically thin, optically thick disks, the simulations of Novikov-Thorne disks have so far employed an ad hoc cooling prescription [278, 226]. This prescription conveniently makes the same assumption that the Novikov-Thorne model does: that all energy dissipated as heat in the disk is radiated away locally on roughly the orbital timescale. This is probably reasonable for an appropriate range of mass accretion rates, though it will be good to test this assumption with future global radiation-MHD simulations.

11.5 ADAFs in simulations

A lot of recent simulation work has been focused on exploring more ADAF-like flows under the action of realistic MRI turbulence [222, 322, 320]. To some extent, this is an extension of the Polish-doughnut simulations of the last decade, yet goes beyond it in at least two important respects: 1) the simulations cover a significantly larger spatial range (a few hundreds of GM/c² versus a few tens); and 2) the simulations explore much longer temporal evolution (hundreds of thousands of GM/c³ versus tens of thousands). The result of this is that the simulations are able to explore steady-state accretion out to much larger radii (∼ 100 GM/c² as opposed to ∼ 10 GM/c²). More can be expected from this work in the coming years.

11.6 Oscillations in simulations

Hydrodynamic simulations have also been used extensively to study the natural oscillation modes of relativistic disks, particularly as they might relate to QPOs (Section 12.4). Rezzolla and collaborators [263, 264] identified a mode they referred to as a p-mode that occurred in a near 3:2 ratio with the radial epicyclic mode, which was subsequently confirmed through numerical simulations [325, 324] (animations can be viewed at [262]). The 3:2 ratio of these modes is important as the highest-frequency QPOs in black hole low-mass X-ray binaries are observed to occur in this ratio. Later, Blaes and collaborators [45] identified a different pair of modes, the vertical epicyclic and axisymmetric breathing mode, that have a near 3:2 ratio for a broader range of parameters. Similar to [325, 324], they also demonstrated that these modes could be identified in numerical simulations.

A number of MHD simulations have also been performed which claim to observe QPOs. Kato [145] performed a global MRI simulation in a pseudo-Newtonian potential and claimed to see QPOs in the power spectrum in a measure of the luminosity associated with radial infall. Machida and Matsumoto [179] have reported the formation of a one-armed spiral within the inner torus of a global MRI simulation, and suggested that this could be responsible for the low frequency QPO in the hard state. Schnittman, Krolik, and Hawley [274] have found tentative evidence for high frequency QPO features in light curves generated by coupling a general relativistic ray tracing code to a global GRMHD simulation with an assumed emission model. Henisey and collaborators [127] found similar tentative evidence in a sample of tilted disk simulations, possibly confirming earlier suggestions that disk tilt may be an important mechanism for driving high frequency QPOs in black hole accretion disks [143, 95]. Most recently Dolence and collaborators [78] report evidence for near-infrared and X-ray QPOs in numerical simulations of Sgr A*.

Aside from these few interesting examples, though, no robust QPO signal has generically emerged in MHD simulations [22, 261]. It is unclear what the implications of this are. It may be that some missing physics, such as radiation transport, plays a fundamental role in exciting QPOs. This is an important open problem in numerical simulations of black hole accretion disks.

11.7 Jets in simulations

In recent years several numerical simulations have demonstrated the generation of jets self-consistently from simulations of disks. A few researchers have claimed to produce jetted outflows from purely hydrodynamic interactions. Nobuta and Hanawa [228], for instance, were able to achieve jetted outflows driven by shock waves when infalling gas with high specific angular momentum collided with the centrifugal barrier. However, such simulations usually require very special starting conditions and the jets tend to be transient features of the flow; therefore, it is unlikely these are related to the well-collimated expansive jets observed in many black hole systems.¹³

Among MHD simulations, some distinction should be drawn between those that impose large scale magnetic fields that extend beyond the domain of the simulation and those that, at least initially, impose a self-contained magnetic field. In the first case, the disk is often treated merely as a boundary condition for the evolution of large-scale magnetic field. Shibata and Uchida [284, 285] used this technique to show that jets could be driven both by a pinching of the fields due to radial inflow through the disk and by the j × B force caused when twisted fields unwind.

Since 1999, many MHD disk simulations, usually starting with locally confined, weak magnetic fields, have shown some sort of magnetically driven bipolar outflow [153, 146, 73, 124, 191, 189]. Since these jets are propagating at the Alfvén speed in a magnetically dominated region (β_m ≪ 1), they proceed supersonically relative to the gas. However, some authors claim the outflow is only a transient state [146] and most do not achieve the very high Lorentz factors observed in astrophysical jets [73]. It is also unclear whether or not the jets in these simulations are connected with the Blandford-Znajek process [73] (Section 10), although see [188]. McKinney [189] and Komissarov [155] suggest that part of the problem is that the jets arising from MHD simulations of disks are often contaminated by the application of numerical floors on density and energy and through numerical diffusion, both of which can lead to anomalous mass loading of the jet. In their own carefully constructed simulations, they independently found that the magnetic structures exhibited by their simulated jets are, in fact, consistent with the expectations of the Blandford-Znajek model, as illustrated in Figure 21 (animations can be viewed at [187]).

Figure 21

On the left is a schematic diagram of the Blandford-Znajek mechanism [49] for an assumed parabolic field distribution. On the right is the result of a numerical simulation from [158] showing a very similar structure. Images reproduced by permission; copyright RAS.

Tchekhovskoy and collaborators [303] have used GRMHD simulations of black hole accretion disks plus jets to investigate possible explanations for the observed radio loud/quiet dichotomy in AGN.¹⁴ For a black hole surrounded by a thin disk, the Blandford-Znajek mechanism predicts the luminosity of the jet should scale as ${L_{BZ}} \propto \Omega _H^2$. This limits the range of power expected for realistic AGN spins to a factor of a few tens at most — too small to explain the observed differences. Thicker disks, on the other hand, produce jets whose luminosity can scale as ${L_{BZ}} \propto \Omega _H^4$ or even $\Omega _H^6$ [303], providing sufficient range to perhaps explain the dichotomy.

There has also been some work recently trying to understand a different jet dichotomy — one that is observed in black hole low-mass X-ray binaries (LMXBs). These exhibit radio emission (associated with jets) whose properties change with the observed X-ray spectral and, to a less well determined extent, temporal properties of the accretion disk [92, 94]. Briefly, compact, steady jets are observed in the Low/Hard state, whereas jets appear absent in the High/Soft state. Fragile and collaborators [102] used numerical simulations to rule out disk scale height as the controlling factor, suggesting instead that perhaps the jets are intimately connected with the corona or failed MHD wind. Alternatively, it could be that magnetic field topology is the key factor [136, 190].

11.8 Highly magnetized accretion in simulations

Recently, work has begun to focus on highly magnetized disk configurations, for which some modes of the MRI may be suppressed. One motivation is that these might provide an alternate explanation for the Low/Hard state [35] (Section 12.3).

One way a highly magnetized state could come about is as a result of a thermal instability in an initially hot, thick, weakly-magnetized disk [179, 101]. Machida and collaborators [179] simulated this process for an optically thin disk, assuming bremsstrahlung cooling and pseudo-Newtonian gravity. They found that, indeed, the densest inner regions of the disk collapse down to a cool, thin, magnetically-supported structure. Fragile and Meier [101] extended these results by including more cooling processes and using a general relativistic MHD code.

Another way to achieve a highly magnetized state is to have the disk “drag” the magnetic field in from some distant region [42, 291]. If the field has a consistent net flux, then doing so will necessarily increase the strength of the field near the hole due simply to the smaller area through which the flux must thread. Recent numerical simulations have demonstrated this effect, first in the pseudo-Newtonian limit [135] (Figure 22), and more recently in relativistic simulations [304, 192]. These examples are interesting because they have led directly to a phenomenon known as a “magnetically arrested” accretion state [218], where the accumulation of magnetic field near the black hole is sufficient to temporarily halt the inflow of matter. This state has been shown to produce jet efficiencies η = (Ṁ − Ė)/〈M〉 > 1 [304], where

$$\dot M(t) = - \int {\sqrt {- g}} \rho {u^r}\;{\rm{d}}\theta \;{\rm{d}}\phi$$

(124)

is the mass accretion rate and

$$\dot E(t) = - \int {\sqrt {- g}} T_t^r\;{\rm{d}}\theta \;{\rm{d}}\phi$$

(125)

is the energy flux, both taken at the black hole event horizon r_H, and the angle brackets indicate a time-averaged quantity. This can only happen if at least some of the energy powering the jet is being extracted from the rotational energy of the black hole itself!

Figure 22

Top: Distributions of density in the meridional plane at different simulation times, showing a magnetically arrested state (left) and a non-arrested state (right). Bottom: Snapshot of magnetic field lines at the same simulation times. Image reproduced by permission from [135], copyright by AAS.

12 Selected Astrophysical Applications

12.1 Measurements of black-hole mass and spin

Astrophysical black holes are not charged, and thus are characterized only by their mass and spin. Measurements of black hole mass are generally straightforward, requiring only the observation of an orbital companion and an application of Kepler’s laws. Current mass estimates for stellar-mass black holes (in particular microquasars) are reviewed by McClintock & Remillard [184] and for supermassive black holes by Kormendy & Richstone [159] (see also [196] for somewhat more recent data).

Nevertheless, there remain several fundamental questions connected with the mass of black holes, and some of them are directly connected to accretion disk theory. One of them is the question of ultra-luminous X-ray sources (ULXs) [61, 181]. ULXs are powerful X-ray sources located outside of galactic nuclei, which have luminosities in excess of the Eddington limit for M = 10² M_⊙, assuming isotropic emission. The huge luminosities of ULXs lead some to conclude that they are accreting intermediate-mass black holes, with masses Mulx > 10² M_⊙ (e.g., [165, 198]). Others think that they are stellar mass black holes with M_ULX ∼ 10 M_⊙, either exhibiting beamed emission [108] or surrounded by disks that are somehow able to produce highly super-Eddington luminosities (e.g., [147, 36]). At this time, at least one ULX (ESO 243-49/HLX-1) has been convincingly demonstrated to have an intermediate mass (∼ 500 M_⊙) [90, 67].

As difficult as it is to nail down the mass on some of these systems, it is even more difficult to measure black hole spin, even though it plays a direct, and important, role in accretion disk physics. One obvious example of the role of spin is the dependence of r_ms, the coordinate radius of the ISCO, on spin. In the symmetry plane of the black hole r_ms = 6r_G for a_* = a/M = 0 (non-rotating), 1r_G for a_* = 1 (maximal prograde rotation), and 9r_G for a_* = −1 (maximal retrograde rotation) (Section 2.5). It is believed that the inner edge of the accretion disk will be similarly affected. In addition, for rapidly rotating black holes, the Blandford-Znajek mechanism, and similar processes which depend on spin, may account for a fair share of the global energetics, comparable to that of accretion itself (see, e.g., [155, 189] and references therein). Therefore, measuring the spin of accreting black holes is integrally tied up in understanding black hole accretion generally.

Four methods for determining black hole spin have been proposed in the literature. With references to some of the earliest results, they are: 1) fitting the continuum spectra of microquasars observed in the thermally dominant state using disk emission models [65, 185, 197, 277]; 2) fitting observed relativistically broadened iron line profiles with theoretical models [142, 315, 201, 202, 260]; 3) matching observed QPO frequencies to those predicted by theoretical models [62, 11, 259, 306]; and 4) analyzing the “shadow” a black hole makes on the surface of an accretion disk [301]. The first three methods have been the most commonly applied to date. Although there have been some glaring discrepancies in the spin estimates published to date (e.g., one group claiming that Cyg X-1 has a near-zero spin, a_* = 0.05±0.01 [200], and then later claiming it has a near-maximal spin, a_* > 0.95 [89]), there appears to be a settling of values in recent years and a growing confidence in the methods, particularly the continuum fitting. There are still some concerns, however. For the continuum fitting method, the main issue is that the inclination of the X-ray emitting region must be measured by some independent means. This is because the effect of the inclination on the spectrum is degenerate with the effect of spin [175], so both can not be accounted for within the continuum fitting method. In cases where such an independent measure is available (e.g. [292]), the continuum fitting method appears robust. For the relativistically-broadened iron line method, there are difficulties in properly estimating the extent of the “red” wing, which is most directly related to the spin of the black hole, and in modeling the hard X-ray source photons and the disk ionization, both of which strongly affect the reflection spectrum. The reviews by Remillard, McClintock, and collaborators [258, 186] give more complete introductions to the topic of measuring black hole spin, with emphasis mainly on the continuum fitting method.

12.2 Black hole vs. neutron star accretion disks

Little of what we have said so far has depended on whether the central compact object is a black hole or neutron star, provided only that the neutron star is compact enough to lie inside the inner radius of the disk r_in. In this case, its presence will not be noticed by the disk except through its gravity, which will be practically the same as for a black hole (an exception would be if the neutron star is strongly magnetized [113]). However, this does not mean that accreting black hole and neutron star sources will be indistinguishable, as we have not yet fully addressed the question of what happens to energy advected past r_in. For optically thick, geometrically thin Shakura-Sunyaev disks (Section 5.3), a significant fraction of the gravitational energy liberated by advection is radiated by the gas prior to it passing through r_in. Thus, the total luminosity of thin disks will not depend sensitively on the nature of the central object. However, this is not the case for the ADAF solution (Section 7), for which much of the thermal energy gained by the gas from accretion is carried all the way in to the central object. Narayan and his collaborators [225, 215, 195, 107] have convincingly argued that this may allow observers to distinguish between black hole and neutron star sources.

The key is that, for black hole sources, advection through the event horizon allows the excess thermal energy to be effectively absorbed without ever radiating. For neutron star sources, on the other hand, the presence of a hard surface ensures that the excess energy of accretion is released upon impact and must be radiated to infinity. This implies that for systems in the ADAF state, a black hole source should be significantly less luminous that a neutron star one with the same mass accretion rate [225]. Perhaps a more important point is that the range of luminosities should be wider for a black hole source than for a neutron star one [215]. This is because, while the luminosity goes as L ∝ ṁ for all neutron star states and for black holes in a high accretion state, it goes as L ∝ ṁ² for black holes in the ADAF state for which ṁ < 10⁻² − 10⁻¹. Furthermore, since the luminosity is also proportional to the mass of the central object L ∝ M, at the highest accretion rates a black hole source should be more luminous than a neutron star one due to its higher mass (Figure 23). Another key point to this argument is that neutron stars can independently and reliably be confirmed if they display type I bursts, which are thermonuclear flashes occurring in material accumulated on the surface of the neutron star [141]. Thus one can compare known neutron star sources against suspected black hole sources. This has now been done in a number of otherwise similar sources and Narayan’s expectations have indeed been confirmed [215, 195, 107] (a recent example is shown in Figure 23). This provides compelling observational evidence for the existence of black hole event horizons, although this falls short of being a proof [12]. This topic is discussed further in the review article by Narayan and McClintock [219].

Figure 23

Left: Luminosity as a function of accretion rate for neutron star and black hole sources, illustrating that a wider range of luminosities are expected for black holes. Image reproduced by permission from [215], copyright by AAS. Right: Recent data showing that neutron star sources (open symbols) are systematically more luminous than black hole sources (filled symbols) in analogous spectral states. Image reproduced by permission from [171], copyright by Elsevier.

12.3 Black-hole accretion disk spectral states

Black hole accretion disks, particularly in X-ray binaries, exhibit complex spectra composed of both thermal and nonthermal components. During outbursts, the relative strengths of these components change frequently in concert with changes in luminosity and the characteristics of the radio features (i.e., jets). Astronomers have developed a set of empirical spectral classification states to broadly characterize these observations. These spectral states likely reveal important information about the underlying physical state of the system; therefore it is worth summarizing the states and their observed properties here.

Probably the easiest state to connect with a theoretical model is the “High/Soft” (HS) or “Thermally Dominant” state. As the name implies, the spectrum in this state is dominated by the thermal component. This state is best explained as ∼ 1 keV thermal emission from a multi-temperature accretion disk, as predicted by the Shakura-Sunyaev (thin) disk model (Section 5.3). However, most sources spend the majority of their lifetimes in the “Low/Hard” (LH) or even “quiescent” state. The quiescent state is characterized by exceptionally low luminosity and a hard, non-thermal spectrum (photon index Γ = 1.5–2.1). As the luminosity increases the sources usually enters the Hard state. Here the 2–10 keV intensity is still comparatively low and the spectrum is still nonthermal. This spectrum is best fit with a power law of photon index Γ ∼ 1.7 (220 keV). In this state, the thermal component is either not detected or appears much cooler, indicating the thin disk may truncate further out than in the Thermally Dominant state, although see [199] for claims that the thin disk extends all the way to the ISCO even in the Hard state. Observations suggest the region interior to the thin disk may be filled with a hot (presumably thick), optically-thin plasma, which accounts for the nonthermal part of the spectrum. This is the picture suggested by the “truncated disk model” [85, 86, 84], which pictures the Hard state as a truncated thin (Shakura-Sunyaev) disk adjoined with an inner thick (ADAF-like) flow. This model has shown tremendous phenomenological success [79], although other models for the Hard state still abound [323, 101, 136]. The Hard state is also linked with observations of a persistent radio jet that is not seen in other states (see Section 10). The final spectral state, which is referred to as the “Very High” (VH) or “Steep Power Law” state, is characterized by the appearance of high-frequency QPOs (Section 12.4) and the presence of both disk and powerlaw components, each of which contributes substantial luminosity. In this state, the powerlaw component is observed to be steep (Γ ∼ 2.5), giving the state its name. This state is sometimes associated with intermittent jets.

These states, their distinguishing observational properties, and a sampling of observations from various black hole X-ray binaries is presented in the review by McClintock and Remillard [184].

12.4 Quasi-Periodic Oscillations (QPOs)

Einstein’s general theory of gravity has never been tested in its strong field limit, characteristic of the region very near black holes (or neutron stars), i.e., a few gravitational radii away from these sources. Soon VLBI measurements may be able to resolve these scales for the supermassive black hole at the center of our galaxy. However, for most sources, resolution in time seems to be a more practical approach. To zeroth order, the light curves from accreting black holes vary in a chaotic manner, resembling a loud noise. However, Fourier analysis of the light curve reveals stinkingly regular patterns buried in the noise. For Galactic black hole and neutron star sources, these quasi-periodic oscillations (QPOs) have frequencies of a few hundred Hz.

Frequencies in the range 100–1000 Hz formally correspond to orbital frequencies a few gravitational radii away from a stellar-mass object. The focus on orbital frequencies is further motivated by the stability of observed QPO frequencies over very long periods of time. For example, QPOs of 300 Hz and 450 Hz were observed from the microquasar GRO J1655–40 during its 1996 outburst and again almost nine years later during its 2005 outburst. This strongly suggests that the frequencies cannot depend on quantities such as magnetic field, density, temperature, or accretion rate, as these all vary greatly in time. The only parameters of a black hole accretion system that do not vary over a nine year period are the mass and the spin of the central black hole. Thus, the oscillation frequencies must only depend on these two parameters, and only frequencies connected to orbital motion have the property that they depend only on mass and spin. Thus, the possible frequencies are: the Keplerian frequency, the two epicyclic frequencies (as originally suggested by Kluźniak and Abramowicz [11]), the Lense-Thirring frequency (as originally suggested by Stella and Vietri [294]), and their combinations (e.g., [145, 143]).

In several microquasars the detected high-frequency QPOs come not as single oscillations, but as part of a pair. Furthermore, Abramowicz and Kluźniak [11] noticed that they are commensurable, being most often in a 2/3 ratio, as shown in Table 2. This suggests that a resonance may be at work. Twin peak QPOs in the kilohertz range have been also detected from binaries containing accreting neutron stars. These neutron star QPOs show a similar, though less obvious, 2/3 ratio.

Table 2 Frequency ratio of the “twin peak” QPOs in all four microquasars where they have been detected.

It was realized [150] that the behavior of the observed QPO frequencies and amplitudes in both neutron star binaries and microquasars is typical for a certain type of non-linear resonance. Indeed, it may be observed when a properly tuned spring is attached to a pendulum with a properly chosen length. Such a system oscillates in two “modes”: the pendulum mode and the spring mode.¹⁵ Mostly due to efforts of Rebusco and Horák [256, 130, 131, 132], a mathematical resonance model was developed to describe an arbitrary system oscillating in two modes near a 2/3 non-linear resonance. The model’s predictions for the frequency and amplitude behaviors are strikingly similar to the ones observed in X-ray binary QPOs, suggesting that they may indeed be explained as a non-linear resonance of two modes of oscillation. The model does not, however, explain what these modes are, how they are excited, nor what energy reservoir they tap. Only when these questions are answered satisfactorily could one say that the QPO puzzle is solved.¹⁶

As a final note, a crucial discovery by Barret and collaborators [32, 33, 34] concerning the behavior of the quality factor of twin peak QPOs proves that they are disk oscillations and cannot be explained by kinematic (Doppler) effects due to the presence “hot spots” on the accretion disk surface. These effects are, however, important in modulating the QPO’s signal [55]. It is not clear, though, whether the oscillations are explained by the discoseismic modes discussed in Section 9.2 (see, e.g., [313] for references).

12.5 The case of Sgr A*

As already mentioned in Section 2.1.1, there is a very good chance that the first direct evidence for a black hole event horizon will come from Sgr A*, the compact, supermassive object at the center of the Milky Way. There have already been a number of strong, indirect arguments in favor of the black hole nature of Sgr A* [54, 53], but no direct evidence yet. Sgr A* is also of interest because it represents a unique case of black hole accretion, having by far the lowest (scaled) mass accretion rate and radiative efficiency of any known source. We anticipate greatly expanding this section in the near future as new results become available.

13 Concluding Remarks

Since the early 1970s, the study of black hole accretion disks has yielded to remarkable successes. And yet, as in most fields of research, each step forward has been met by new questions. In this article, we have tried to give a tour of some of the successes, such as the many disk models (thick, thin, slim, ADAF, …) that have given a firm foundation on which to work, the many studies of disk instabilities and oscillations that help us to understand the ways in which real disks can deviate from the simplistic models, and the numerical simulations that come as close as possible to an experimental test-bed for black hole accretion. We have also tried to indicate what we believe are some of the most pressing challenges of the day, including matching our theoretical knowledge to actual observed phenomena such as black hole spectral states, quasi-periodic oscillations, and relativistic jets. There are also observational challenges to find direct evidence of black hole event horizons and definitively constrain black hole spins. We hope, as we continue to update this Living Review, to be able to report on future discoveries in these areas, just as we expect to report new puzzles we have yet to encounter.

Clearly our tour has been incomplete. For instance, despite their prominent role in nature, e.g., as an AGN feedback mechanism, outflows are not accounted for in any of the four main accretion models we presented, as the models all assume that the accretion rate is constant with radius. In reality, outflows may be triggered by any of three mechanisms: thermal, radiative, or centrifugal. Thermal winds are expected to result from heating of the outer regions of an accretion disk by its hot inner region. Radiative winds are driven by radiative flux acting on line opacities. Centrifugal acceleration of particles can take place along magnetic field lines which are sufficiently inclined to the disk plane.

The ADIOS (adiabatic inflow-outflow solution) [47] is a generalization of the ADAF solution that actually has much of the mass, energy, and angular momentum of the accretion “disk” carried away in the form of winds, rather than being advected into the black hole as in a normal ADAF. However, the argument behind the ADIOS model is flawed [14]. Blandford and Begelman [47] claim that black hole accretion flows with small radiative efficiency must necessarily experience strong outflows because the matter all has a positive Bernoulli constant. Yet a positive Bernoulli constant is only a necessary, but not sufficient, condition for outflows. For example, the classical Bondi accretion solution has ${\mathcal B} > 0$ everywhere and yet experiences no outflows. Furthermore, low efficiency accretion flows do not really have a positive Bernoulli constant everywhere, as boundary conditions (ignored in [47]) will impose some regions of negative Bernoulli constant. A more recent look at the ADIOS solution is presented in [37].

There are several other analytic and semi-analytic models of accretion disks. Some are closely related to the models we have already discussed. For example, the CDAF (convection-dominated accretion flow) [217, 221] is another variant on the ADAF, in which long-wavelength convective instabilities transport angular momentum inward and energy outward. Other models relax some of the standard assumptions about accretion disks. For example, Bardeen and Petterson [30] and others [167, 273, 177] relaxed the assumption that disks are axisymmetric by considering tilted accretion disks, acted on by the Lense-Thirring precession of the central (rotating) black hole. Then there is the exact solution for stationary, axisymmetric non-circular, accretion flows found by Kluźniak and Kita [151].

On the numerical side, too, there are many interesting accretion configurations that have been identified, but are not included in this review. Some examples include quasi-spherical (low angular momentum) accretion flows [250, 249] and convection-dominated disks [138].

For those who wish for more details or a different perspective, we can recommend several excellent text books and review articles devoted, partially or fully, to black hole accretion disks. We recall here some of the most often quoted. The oldest, but still very useful and informative, is the classic review by Pringle [248]. The most authoritative text book on accretion is Accretion Power in Astrophysics by Frank, King and Raine [103]. Two monographs devoted to black hole accretion disks are: Black-Hole Accretion Disks by Kato, Fukue and Mineshige [144], and Theory of Black Hole Accretion Disks by Abramowicz, Björnsson, and Pringle [3]. Lasota [170] also wrote an excellent non-technical Scientific American article on black hole accretion in microquasars. Finally, there is a nice series of lecture notes by Ogilvie available on the web [230].