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  (see  for a review of instabilities in astrophysical accretion disks). Since the pioneering work by Shakura and Sunyaev , 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.  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 ). 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  and references therein). The importance of magnetic fields and their consequences for the stability properties of this solution are critically discussed in .
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., ). 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 , 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  and relativistic gravity . 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 . In  (see also ), 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  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 .
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  (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 . The code contains a parametrized viscosity based on the -model . 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  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  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. . 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.,  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 , 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,  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  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.  (animations can be visualized at the website listed in ). Furthermore, recent simulations with non-constant angular momentum disks and rotating black holes , 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 , 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 . The Kerr metric is described in Boyer-Lindquist coordinates, with a radial tortoise coordinate to enhance the resolution near the horizon.
In , 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.  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  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 . 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. . In this work, Wilson's formulation of the hydrodynamic equations was adopted. The integration algorithm was borrowed from  with the transport terms finite-differenced following the prescription given in . 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  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 , 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 , motivated by the non-stationary patterns found in Newtonian simulations (see, e.g., ). 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 , 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 , 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., ). 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 . Papadopoulos and Font  used the limiting case for Schwarzschild black holes, i.e., the inhomogeneous Bardeen-Press equation . 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  for a detailed explanation).
One-dimensional numerical simulations of a self-gravitating perfect fluid accreting onto a black hole were presented in , 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  performed the simulations adopting an ingoing null foliation of a spherically symmetric black hole spacetime . 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.  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  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 ). Zanotti et al.  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
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to firstname.lastname@example.org