Disk accretion theory is primarily based on the study of (viscous) stationary flows and their stability properties through linearized perturbations thereof. A well-known example is the solution consisting of isentropic constant angular momentum tori, unstable to a variety of non-axisymmetric global modes, discovered by Papaloizou and Pringle [223] (see [18] for a review of instabilities in astrophysical accretion disks). Since the pioneering work by Shakura and Sunyaev [253], thin disk models, parametrized by the so-called -viscosity, in which the gas rotates with Keplerian angular momentum which is transported radially by viscous stress, have been applied successfully to many astronomical objects. The thin disk model, however, is not valid for high luminosity systems, as it is unstable at high mass accretion rates. In this regime, Abramowicz et al. [3] found the so-called slim disk solution, which is stable against viscous and thermal instabilities. More recently, much theoretical work has been devoted to the problem of slow accretion, motivated by the discovery that many galactic nuclei are under-luminous (e.g., NGC 4258). Proposed accretion models involve the existence of advection-dominated accretion flows (ADAF solution; see, e.g., [202, 200]) and advection-dominated inflow-outflow solutions (ADIOS solution [33]). The importance of convection for low values of the viscosity parameter is currently being discussed in the so-called convection-dominated accretion flow (CDAF solution; see [130] and references therein). The importance of magnetic fields and their consequences for the stability properties of this solution are critically discussed in [17].

For a wide range of accretion problems, the Newtonian theory
of gravity is adequate for the description of the background
gravitational forces (see, e.g., [98]). Extensive experience with Newtonian astrophysics has shown
that explorations of the relativistic regime benefit from the use
of model potentials. In particular, the Paczynski-Wiita
pseudo-Newtonian potential for a Schwarzschild black hole [217], gives approximations of general relativistic effects with an
accuracy of
outside the
*marginally stable*
radius (at
*r*
= 6
*M*, or three times the Schwarzschild radius). Nevertheless, for
comprehensive numerical work a three-dimensional formalism is
required, able to cover also the maximally rotating hole. In
rotating spacetimes the gravitational forces cannot be captured
fully with scalar potential formalisms. Additionally, geometric
regions such as the ergosphere would be very hard to model
without a metric description. Whereas the bulk of emission occurs
in regions with almost Newtonian fields, only the observable
features attributed to the inner region may crucially depend on
the nature of the spacetime.

In the following, we present a summary of representative
*time-dependent*
accretion simulations in relativistic hydrodynamics. We
concentrate on multi-dimensional simulations. In the limit of
spherical accretion, exact stationary solutions exist for both
Newtonian gravity [43] and relativistic gravity [179]. Such solutions are commonly used to calibrate time-dependent
hydrodynamical codes, by analyzing whether stationarity is
maintained during a numerical evolution [123,
163,
78,
244,
21].

Wilson's formulation has been extensively used in time-dependent numerical simulations of thick disk accretion. In a system formed by a black hole surrounded by a thick disk, the gas flows in an effective (gravitational plus centrifugal) potential, whose structure is comparable to that of a close binary. The Roche torus surrounding the black hole has a cusp-like inner edge located at the Lagrange point , where mass transfer driven by the radial pressure gradient is possible [1]. In [123] (see also [120]), Hawley and collaborators studied, in the test fluid approximation and axisymmetry, the evolution and development of nonlinear instabilities in pressure-supported accretion disks formed as a consequence of the spiraling infall of fluid with some amount of angular momentum. The code used explicit second-order finite difference schemes with a variety of choices to integrate the transport terms of the equations (i.e., those involving changes in the state of the fluid at a specific volume in space). The code also used a staggered grid (with scalars located at the cell centers and vectors at the cell boundaries) for its suitability to finite difference the transport equations in a stable numerical way. Discontinuous solutions were stabilized with artificial viscosity terms.

With a three-dimensional extension of the axisymmetric code of Hawley et al. [122, 123], Hawley [121] studied the global hydrodynamic non-axisymmetric instabilities in thick, constant angular momentum accretion gas tori orbiting around a Schwarzschild black hole. Such simulations showed that the Papaloizu-Pringle instability saturates in a strong spiral pressure wave, not in turbulence. In addition, the simulations confirmed that accretion flows through the torus could reduce and even halt the growth of the global instability. Extensions to Kerr spacetimes have been recently reported in [62].

Yokosawa [313,
314], also using Wilson's formulation, studied the structure and
dynamics of relativistic accretion disks and the transport of
energy and angular momentum in magneto-hydrodynamical accretion
onto a rotating black hole. In his code, the hydrodynamic
equations are solved using the Flux-Corrected-Transport (FCT)
scheme [46] (a second-order flux-limiter method that avoids oscillations
near discontinuities by reducing the magnitude of the numerical
flux), and the magnetic induction equation is solved using the
constrained transport method [80]. The code contains a parametrized viscosity based on the
-model [253]. The simulations revealed different flow patterns inside the
marginally stable orbit, depending on the thickness
*h*
of the accretion disk. For thick disks with
,
being the radius of the event horizon, the flow pattern becomes
turbulent.

Igumenshchev and Beloborodov [131] have performed two-dimensional relativistic hydrodynamical simulations of inviscid transonic disk accretion onto a Kerr black hole. The hydrodynamical equations follow Wilson's formulation, but the code avoids the use of artificial viscosity. The advection terms are evaluated with an upwind algorithm that incorporates the PPM scheme [58] to compute the fluxes. Their numerical work confirms analytic expectations: (i) The structure of the innermost disk region strongly depends on the black hole spin, and (ii) the mass accretion rate goes as , being the potential barrier between the inner edge of the disk and the cusp, and the adiabatic index.

Furthermore, thick accretion disks orbiting black holes may be
subjected to the so-called
*runaway instability*, as first proposed by Abramowicz et al. [2]. Starting with an initial disk filling its Roche lobe, so that
mass transfer is possible through the cusp located at the
Lagrange point, two evolutions are feasible when the mass of the
black hole increases: (i) Either the cusp moves inwards towards
the black hole, which slows down the mass transfer, resulting in
a stable situation, or (ii) the cusp moves deeper inside the disk
material. In the latter case, the mass transfer speeds up,
leading to the runaway instability. This instability, whose
existence is still a matter of debate (see, e.g., [87] and references therein), is an important issue for most models
of cosmic GRBs, where the central engine responsible for the
initial energy release is such a system consisting of a thick
disk surrounding a black hole.

In [87], Font and Daigne presented time-dependent simulations of the runaway instability of constant angular momentum thick disks around black holes. The study was performed using a fully relativistic hydrodynamics code based on HRSC schemes and the conservative formulation discussed in Section 2.1.3 . The self-gravity of the disk was neglected and the evolution of the central black hole was assumed to be that of a sequence of Schwarzschild black holes of varying mass. The black hole mass increase is determined by the mass accretion rate across the event horizon. In agreement with previous studies based on stationary models, [87] found that by allowing the mass of the black hole to grow the disk becomes unstable. For all disk-to-hole mass ratios considered (between 1 and 0.05), the runaway instability appears very fast on a dynamical timescale of a few orbital periods (), typically a few and never exceeding , for a black hole and a large range of mass fluxes ().

An example of the simulations performed by [87] appears in Figure
7
. This figure shows eight snapshots of the time-evolution of the
rest-mass density, from
*t*
=0 to
. The contour levels are linearly spaced with
, where
is the maximum value of the density at the center of the initial
disk. In Figure
7
one can clearly follow the transition from a quasi-stationary
accretion regime (panels (1) to (5)) to the rapid development of
the runaway instability in about two orbital periods (panels (6)
to (8)). At
, the disk has almost entirely disappeared inside the black hole
whose size has, in turn, noticeably grown.

Extensions of this work to marginally stable (or even stable) constant angular momentum disks are reported in Zanotti et al. [239] (animations can be visualized at the website listed in [236]). Furthermore, recent simulations with non-constant angular momentum disks and rotating black holes [86], show that the instability is strongly suppressed when the specific angular momentum of the disk increases with the radial distance as a power law, . Values of much smaller than the one corresponding to the Keplerian angular momentum distribution suffice to supress the instability.

The most promising processes for producing relativistic jets like those observed in AGNs, microquasars, and GRBs involve the hydromagnetic centrifugal acceleration of material from the accretion disk [34], or the extraction of rotational energy from the ergosphere of a Kerr black hole [225, 36]. Koide and coworkers have performed the first MHD simulations of jet formation in general relativity [141, 140, 142]. Their code uses the 3+1 formalism of general relativistic conservation laws of particle number, momentum, and energy, and Maxwell equations with infinite electric conductivity. The MHD equations are numerically solved in the test fluid approximation (in the background geometry of Kerr spacetime) using a finite difference symmetric scheme [61]. The Kerr metric is described in Boyer-Lindquist coordinates, with a radial tortoise coordinate to enhance the resolution near the horizon.

In [142], the general relativistic magneto-hydrodynamic behaviour of a
plasma flowing into a rapidly rotating black hole (*a*
=0.99995) in a large-scale magnetic field is investigated
numerically. The initial magnetic field is uniform and strong,
,
being the mass density. The initial Alfvén speed is
. The simulation shows how the rotating black hole drags the
inertial frames around in the ergosphere. The azimuthal component
of the magnetic field increases because of the azimuthal twisting
of the magnetic field lines, as is depicted in Figure
8
. This frame-dragging dynamo amplifies the magnetic field to a
value which, by the end of the simulation, is three times larger
than the initial one. The twist of the magnetic field lines
propagates outwards as a torsional Alfvén wave. The magnetic
tension torques the plasma inside the ergosphere in a direction
opposite to that of the black hole rotation. Therefore, the
angular momentum of the plasma outside receives a net increase.
Even though the plasma falls into the black hole, electromagnetic
energy is ejected along the magnetic field lines from the
ergosphere, due to the propagation of the Alfvén wave. By total
energy conservation arguments, Koide et al. [142] conclude that the ultimate result of the generation of an
outward Alfvén wave is the magnetic extraction of rotational
energy from the Kerr black hole, a process the authors call the
MHD Penrose process. Koide and coworkers argue that such a
process can be applicable to jet formation, both in AGNs and
GRBs.

We note that, recently, van Putten and Levinson [292] have considered, theoretically, the evolution of an accretion torus in suspended accretion, in connection with GRBs. These authors claim that the formation of baryon-poor outflows may be associated with a dissipative gap in a differentially rotating magnetic flux tube supported by an equilibrium magnetic moment of the black hole. Numerical simulations of non-ideal MHD, incorporating radiative processes, are necessary to validate this picture.

Numerical simulations have extended the simplified analytic
models and have helped to develop a thorough understanding of the
hydrodynamic accretion scenario, in its fully three-dimensional
character (see, e.g., [245,
27] and references therein). The numerical investigations have
revealed the formation of accretion disks and the appearance of
non-trivial phenomena such as shock waves and tangential (*flip-flop*) instabilities.

Most of the existing numerical work has used Newtonian
hydrodynamics to study the accretion onto non-relativistic
stars [245]. For compact accretors such as neutron stars or black holes,
the correct numerical modeling requires a general relativistic
hydrodynamical description. Within the relativistic, frozen star
framework, wind accretion onto ``moving'' black holes was first
studied in axisymmetry by Petrich et al. [226]. In this work, Wilson's formulation of the hydrodynamic
equations was adopted. The integration algorithm was borrowed
from [277] with the transport terms finite-differenced following the
prescription given in [123]. An artificial viscosity term of the form
, with
*a*
a being constant, was added to the pressure terms of the
equations in order to stabilize the numerical scheme in regions
of sharp pressure gradients.

An extensive survey of the morphology and dynamics of relativistic wind accretion past a Schwarzschild black hole was later performed by [91, 90]. This investigation differs from [226] in both the use of a conservative formulation for the hydrodynamic equations (see Section 2.1.3) and the use of advanced HRSC schemes. Axisymmetric computations were compared to [226], finding major differences in the shock location, opening angle, and accretion rates of mass and momentum. The reasons for the discrepancies are related to the use of different formulations, numerical schemes, and grid resolution, much higher in [91, 90]. Non-axisymmetric two-dimensional studies, restricted to the equatorial plane of the black hole, were discussed in [90], motivated by the non-stationary patterns found in Newtonian simulations (see, e.g., [27]). The relativistic computations revealed that initially asymptotic uniform flows always accrete onto the hole in a stationary way that closely resembles the previous axisymmetric patterns.

In [219], Papadopoulos and Font presented a procedure that simplifies
the numerical integration of the general relativistic
hydrodynamic equations near black holes. This procedure is based
on identifying classes of coordinate systems in which the black
hole metric is free of coordinate singularities at the horizon
(unlike the commonly adopted Boyer-Lindquist coordinates),
independent of time, and admits a spacelike decomposition. With
those coordinates the innermost radial boundary can be placed
*inside*
the horizon, allowing for an unambiguous treatment of the entire
(exterior) physical domain. In [94,
95] this approach was applied to the study of relativistic wind
accretion onto rapidly rotating black holes. The effects of the
black hole spin on the flow morphology were found to be confined
to the inner regions of the black hole potential well. Within
this region, the black hole angular momentum drags the flow,
wrapping the shock structure around. An illustrative example is
depicted in Figure
9
. The left panel of this figure corresponds to a simulation
employing the Kerr-Schild form of the Kerr metric, regular at the
horizon. The right panel shows how the accretion pattern would
look if it were the computation performed using the more common
Boyer-Lindquist coordinates. The transformation induces a
noticeable wrapping of the shock around the central hole. The
shock would wrap infinitely many times before reaching the
horizon. As a result, the computation in these coordinates would
be much more challenging than in Kerr-Schild coordinates.

In [220], such conclusions were corroborated with numerical simulations of the gravitational radiation emitted during the accretion process of an extended object onto a black hole. The first-order deviations from the exact black hole geometry were approximated using curvature perturbations induced by matter sources whose nonlinear evolution was integrated using a (nonlinear) hydrodynamics code (adopting the conservative formulation of Section 2.1.3 and HRSC schemes). All possible types of curvature perturbations are captured in the real and imaginary parts of the Weyl tensor scalar (see, e.g., [55]). In the framework of the Newman-Penrose formalism, the equations for the perturbed Weyl tensor components decouple, and when written in the frequency domain, they even separate [285]. Papadopoulos and Font [220] used the limiting case for Schwarzschild black holes, i.e., the inhomogeneous Bardeen-Press equation [23]. The simulations showed the gradual excitation of the black hole quasi-normal mode frequency by sufficiently compact shells.

An example of these simulations appears in the movie of
Figure
10
. This movie shows the time evolution of the shell density (left
panel) and the associated gravitational waveform during a
complete accretion/collapse event. The (quadrupolar) shell is
parametrized according to
with
and
. Additionally,
denotes a logarithmic radial (Schwarzschild) coordinate. The
animation shows the gradual collapse of the shell onto the black
hole. This process triggers the emission of gravitational
radiation. In the movie, one can clearly see how the burst of the
emission coincides with the most dynamical accretion phase, when
the shell crosses the peak of the potential and is subsequently
captured by the hole. The gravitational wave signal coincides
with the quasinormal ringing frequency of the Schwarzschild black
hole, 17
*M*
. The existence of an initial burst, separated in time from the
*physical*
burst, is also noticeable in the movie. It just reflects the
gravitational radiation content of the initial data (see [220] for a detailed explanation).

One-dimensional numerical simulations of a self-gravitating
perfect fluid accreting onto a black hole were presented
in [222], where the effects of mass accretion during the gravitational
wave emission from a black hole of growing mass were explored.
Using the conservative formulation outlined in Section
2.2.2
and HRSC schemes, Papadopoulos and Font [222] performed the simulations adopting an ingoing null foliation of
a spherically symmetric black hole spacetime [221]. Such a foliation penetrates the black hole horizon, allowing
for an unambiguous numerical treatment of the inner boundary. The
essence of non-spherical gravitational perturbations was captured
by adding the evolution equation for a minimally coupled massless
scalar field to the (characteristic) Einstein-perfect fluid
system. The simulations showed the familiar damped-oscillatory
radiative decay, with both decay rate and frequencies being
modulated by the mass accretion rate. Any appreciable increase in
the horizon mass during the emission reflects on the
instantaneous signal frequency
*f*, which shows a prominent negative branch in the
evolution diagram. The features of the frequency evolution
pattern reveal key properties of the accretion event, such as the
total accreted mass and the accretion rate.

Recently, Zanotti et al. [317] have performed hydrodynamical simulations of constant angular momentum thick disks (of typical neutron star densities) orbiting a Schwarzschild black hole. Upon the introduction of perturbations, these systems either become unstable to the runaway instability [87] or exhibit a regular oscillatory behaviour resulting in a quasi-periodic variation of the accretion rate as well as of the mass quadrupole (animations can be visualized at the website listed in [236]). Zanotti et al. [317] have found that the latter is responsible for the emission of intense gravitational radiation whose amplitude is comparable or larger than the one expected in stellar core collapse. The strength of the gravitational waves emitted and their periodicity are such that signal-to-noise ratios can be reached for sources at or , respectively, making these new sources of gravitational waves potentially detectable.

Numerical Hydrodynamics in General Relativity
José A. Font
http://www.livingreviews.org/lrr-2003-4
© Max-Planck-Gesellschaft. ISSN 1433-8351 Problems/Comments to livrev@aei-potsdam.mpg.de |