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 nonaxisymmetric 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, parameterized by the 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 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., [285, 284]) and advection-dominated inflow-outflow solutions (ADIOS solution ). The importance of convection for low values of the viscosity parameter has also been discussed in the 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 . Magnetic effects have also been long advocated in connection with outflows from black-hole–accretion-disk systems, particularly when the mass accretion rate is much larger than the Eddington rate. Recent Chandra observations of the stellar-mass black-hole binary GRO J1655-40  provide strong evidence supporting the magnetic nature of disk accretion onto black holes and its associated magnetic-powered winds. Theoretically, MHD outflows (and relativistic jets) are believed to be driven by energy flow through the accretion disk associated with the magnetic torque. The most convincing argument for the existence of the torque is the presence of the MRI, which not only drives the turbulent transport of angular momentum within the disk, self-regulating the accretion process, but also amplifies the strength of a seed vertical field on a dynamic timescale given by the inverse of the angular frequency of the disk. Extensive numerical work has been carried out in the last few years for a better understanding of such complex dynamics (see below).
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 Paczyński–Wiita pseudo-Newtonian potential for a Schwarzschild black hole  gives approximations of general-relativistic effects with an accuracy of 10 – 20% outside the marginally stable radius (at ). Nevertheless, for comprehensive numerical work, a three-dimensional formalism is required, that is also able to cover the maximally-rotating black 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 (as the relativistic X-ray emission lines) 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 and MHD. We concentrate on multidimensional 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 hydrodynamic codes, by analyzing whether stationarity is maintained during a numerical evolution [171, 237, 117, 339, 36]. Correspondingly, similar tests have been proposed for the GRMHD equations, such as the spherically-symmetric accretion solution of a perfect fluid onto a Schwarzschild black hole in the presence of an (unphysical) radial magnetic field () or the equatorial magnetized inflow solution in the region between a Kerr black-hole horizon and the marginally stable orbit, derived by  (see also ). These analytic solutions have been used in the literature to assess numerical codes [149, 86, 108, 24, 91].
Pioneering numerical efforts in the study of black-hole accretion [411, 171, 167, 168] made use of the frozen star paradigm of a black hole. In this framework, the time “slicing” of the spacetime is synchronized with that of asymptotic observers far from the black hole. Within this approach, Wilson  first investigated numerically the time-dependent accretion of inviscid matter onto a rotating black hole. This was the first problem to which his formulation of the hydrodynamic equations, as presented in Section 2.1.2, was applied. Wilson used an axisymmetric hydrodynamic code in cylindrical coordinates to study the formation of shock waves and X-ray emission in the strong-field regions close to the black-hole horizon, and was able to follow the formation of thick accretion disks during the simulations.
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 cell centers and vectors at the cell boundaries) for its suitability to difference the transport equations. Discontinuous solutions were stabilized with artificial viscosity terms.
With a three-dimensional extension of the axisymmetric code of Hawley et al. [170, 171], Hawley  studied the global hydrodynamic nonaxisymmetric 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 are reported in .
Igumenshchev and Beloborodov  performed two-dimensional relativistic hydrodynamic simulations of inviscid transonic disk accretion onto a Kerr black hole. The hydrodynamic 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 regarding the strong dependence of the structure of the innermost disk region on the black-hole spin, and the functional form of the mass accretion rate with the potential barrier between the inner edge of the disk and the cusp, , with being the adiabatic index.
Thick accretion disks orbiting black holes, on the other hand, may be subject to the 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: either (i) 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 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–black-hole mass ratios considered (between 1 and 0.05), the runaway instability appears very fast on a dynamic timescale of a few orbital periods (1 100), typically a few 10 ms and never exceeding 1 s, for a black hole and a large range of mass fluxes (). Extensions of this work to marginally stable (or even stable) constant angular momentum disks were reported in Zanotti et al.  (animations can be visualized at the website listed in ).
An example of the simulations performed by  appears in Figure 10. This figure shows eight snapshots of the time-evolution of the rest-mass density, from 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 10 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.
On the other hand, recent simulations with nonconstant angular momentum disks and rotating black holes [128, 82], show that the runaway instability is strongly suppressed in such cases. In the models built and evolved by [128, 82] the angular momentum of the disks is assumed to increase outwards with the radial distance according to a power law, , where and are constant, such that corresponds to prograde disks (with respect to the black-hole rotation) and to retrograde disks. The results of  show that nonconstant angular momentum disks are dramatically stabilized for very small values of the angular momentum slope , much smaller than the Keplerian value . In light of these results further ingredients need to be incorporated into the codes to finally discern the likelihood of the instability, namely magnetic fields and self-gravity.
The impact of magnetic fields on the dynamics of accretion disks around black holes has long been recognized, but only recently have the first systematic GRMHD simulations become possible [86, 87, 85, 149, 200, 173, 24, 202, 251, 368, 256, 269]. We note, however, that as early as 1977 Wilson  developed the first numerical GRMHD scheme to study the axisymmetric evolution of plasma near a Kerr black hole. A key idea under intense numerical scrutiny currently is that the combined presence of a magnetic field and of a differentially rotating Keplerian disk, can lead to the magneto-rotational instability , generating effective viscosity and transporting angular momentum outwards through MHD turbulence. The dynamo generated by the MRI amplifies the magnetic field on dynamic timescales and, in turn, the associated magnetic torque is believed to be responsible for the generation of strong MHD outflows (see below). The sustained level of activity in the numerical modelling of magnetized accretion disks has also been accompanied by advances in the analytic front. In this respect, Komissarov  has derived an analytic solution for an axisymmetric, stationary, barotropic torus with constant distribution of specific angular momentum and a toroidal magnetic field configuration, which can be used to test multidimensional GRMHD codes in the strong gravity regime.
Yokosawa [429, 430], using Wilson’s formulation, studied the structure and dynamics of relativistic accretion disks and the transport of energy and angular momentum in MHD accretion on to a rotating black hole. In his code, the hydrodynamic equations are solved using the 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 parameterized viscosity based on the -model . The simulations revealed different flow patterns inside the marginally-stable orbit, depending on the thickness, , of the accretion disk. For thick disks with , being the radius of the event horizon, the flow pattern becomes turbulent.
The evolution of the MRI in both the Schwarzschild and Kerr metrics has been numerically investigated, either in axisymmetry or in three dimensions, in a number of recent papers [86, 149, 87, 173]. The simulations have been performed with two of the codes listed in Table 2, namely HARM  and the code of De Villiers and Hawley . We note that both codes not only implement different formulations of the GRMHD equations but they also use entirely different numerical schemes (conservative and incomplete Riemann solver the former, nonconservative and artificial viscosity the latter). Furthermore, the code of  uses Kerr–Schild coordinates (modified to increase the resolution where the disk lies), regular at the event horizon (see below), while the code of  employs the more traditional Boyer–Lindquist coordinate system (singular at the horizon). We note that the initial models in all these investigations are purely hydrodynamic, constant–angular-momentum, thick disks in equilibrium, with an ad hoc poloidal magnetic field distribution superposed. The simulations, carried out for a large sample of parameter values (such as black-hole spin or initial magnetization of the plasma) reveal the following basic evolution of the accretion process: first, there is an initial rapid growth of the toroidal magnetic field driven by shear. This is followed by the development of the MRI until saturation. The accretion process then continues through the quasi-stationary evolution of the disk by sustained MHD turbulence, which redistributes the density and angular momentum of the disk. An example of such an evolution is depicted in Figure 11, taken from . This figure shows the first (left) and last snapshot (after about 8 orbital periods) of the logarithm of the density field for a magnetized, constant–angular-momentum torus around a Kerr black hole with . The exponential growth of the MRI happens within the first orbital period, after which the magnetic field has been amplified to high enough levels to distort the torus and initiate the accretion process.
De Villiers et al.  are able to identify five generic features in the accretion flow: a turbulent and dense disk body with a roughly constant ( and being the pressure scale height and radius, respectively), an inner disk consisting of an inner torus and a plunging inflow, a magnetized coronal envelope, an evacuated axial funnel, and a jet-like outflow along the funnel wall. Overall, they found a decrease in the accretion rate into the black hole with increasing spin parameter. The accretion rate is highly variable in all models, which makes the characteristics of the inner torus highly variable as well, which plays a key role in the accretion flow into the black hole. The outflows are found to be launched from a region in the coronal envelope above the inner torus, which is closer to the horizon as the black-hole spin increases, resulting in more powerful jets.
Three-dimensional numerical work has also addressed recently the evolution of accretion disks that are misaligned (tilted) with respect to the rotation axis of a Kerr black hole. Studies have been carried out both for hydrodynamic flows  and for magnetized flows . These studies are motivated by the increasing evidence that misaligned black holes may actually exist in nature, evidence collected both from observational data as from theoretical arguments (see  and references therein). The simulations have been done using the cosmos and cosmos++ codes listed in Tables 1 and 2, respectively. For hydrodynamic flows it was found that Lense–Thirring precession causes the disks to warp differentially, yet the precession timescale is not shorter than other dynamic timescales. The late-time evolution of the disks showed solid-body precession at rates consistent with some low frequency quasi-periodic oscillations (QPOs). Similarly, the MRI turbulent tilted-disk simulations have not shown indication of a Bardeen–Petterson effect at large.
Recent work on disk accretion has also focused on improving the microphysics of the tori, as in the axisymmetric magnetohydrodynamic simulations of  for neutrino-cooled accretion tori around rotating black holes. (See also  for comparisons between a simple ideal fluid EOS and the relativistic Synge EOS in GRMHD accretion flows.) The EOS used takes into account neutronization, the nuclear statistical equilibrium of a gas of free nucleons and -particles, black body radiation, and a relativistic Fermi gas (neutrinos, electrons, and positrons), along with several cooling processes. It is found that the neutrino luminosity is for a torus of mass and a black hole of , irrespective of its spin, a figure compatible with current estimates in short-hard GRBs. However, the luminosity decreases with decreasing torus mass, and, in light of recent state-of-the-art simulations of neutron-star binary mergers (see below), the formed black-hole–torus systems have generally small torus masses (). On the positive side, however, the simulations predict a strong variability in the neutrino luminosity (on timescales of a few ms), due to the MRI-turbulent accretion flow, which could help explain the observed variability of GRB light curves.
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 by magnetic processes [316, 50]. The Blandford–Payne model  is based upon the existence of a large-scale poloidal magnetic field passing through the disk, whose matter dynamics is dominated by magnetic tension. Following angular-momentum–conservation considerations, accelerated outflows are produced whenever the magnetic field lines have a sufficient angle with respect to the disk. Correspondingly, in the Blandford–Znajek mechanism  the black hole is seen as a magnetized rotating conductor whose rotational energy can be efficiently extracted by means of a magnetic torque, giving rise to a Poynting flux jet. On the other hand, in the magnetic Penrose process  the magnetic field lines across the ergosphere are twisted by frame dragging. The twist of the lines propagates outwards as a torsional Alfvén wave train carrying electromagnetic energy and making the total energy of the plasma near the black hole decrease to negative values. In addition, the accretion of this plasma by the black hole reduces the black-hole rotational energy. As in the Blandford–Znajek process, the type of outflows generated by this process is also a Poynting flux jet. The viability of all these various mechanisms has been investigated numerically in the last few years.
Koide et al. performed the first MHD simulations of jet formation in general relativity [196, 195, 197, 194, 288]. 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 [197, 194], in particular, the GRMHD behavior of a plasma flowing into a rapidly-rotating black hole () in a large-scale magnetic field was 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 azimuthal twisting of the magnetic field lines, as depicted in Figure 12. 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 (but see ).
The first time-dependent GRMHD simulations of the magnetically-dominated monopole magnetospheres of black holes were performed by , who found that the numerical solution evolves towards a stable steady-state solution, which is very close to the corresponding force-free solution found by Blandford & Znajek . A similar result was reported by  for weakly magnetized accretion disks around Kerr black holes, who found numerical solutions consistent with the Blandford–Znajek solution in the low-density funnel region around the black hole. The numerical simulations of  showed for the first time the development of an ultrarelativistic particle wind, Poynting dominated, from a rotating black hole.
Further numerical work on the Blandford–Znajek mechanism was pursued by Komissarov in  using initial models similar to those employed by Koide et al. [197, 194]. The long-term solution shows properties that are significantly different from those of the initial transient phase studied by [197, 194]. In particular, no regions of negative hydrodynamic ‘energy at infinity’ are seen inside the ergosphere and the MHD Penrose process does not operate. Similar conclusions were drawn by , which highlights the fact that the models of [197, 194] were evolved for too short a time to observe unbound mass outflows along the funnel region (but see  for longer evolutions). The results of Komissarov  indicate that the rotational energy of the black hole continues to be extracted via the purely electromagnetic Blandford–Znajek mechanism. This effect also operates when the black-hole spin is the maximum allowed, (extreme Kerr), as has been recently shown by . Perhaps the most important result of these axisymmetric models has been not finding strong relativistic outflows from the black hole, contrary to observational data indicating the existence of strongly relativistic jets in numerous scenarios (but see below).
Ultrarelativistic outflows have not been found either in the comprehensive study of Hawley et al. [87, 173, 88, 210, 169]. Their third paper on the series , as well as , are devoted to studying the formation of unbound outflows in the simulations. Such unbound outflows, which are produced self-consistently from accretion disks that do not have large-scale magnetic fields as an initial condition, do not require rotation in the black hole, yet their strength is greatly enhanced by the black-hole spin. If the spin of the black hole is high, as in their almost extreme Kerr (Kepler disk) model KDJ with , the energy ejected in the Poynting-flux–dominated outflows can be tens of percent of the accreted rest mass. At low spin, kinetic energy and enthalpy of the matter dominate the outflow energetics. These basic results are in qualitative agreement with axisymmetric simulations performed by  using the conservative, shock-capturing scheme HARM. The jets formed in the three-dimensional simulations have two major components: a matter-dominated outflow, which moves at moderate speed () along the centrifugal barrier surrounding an evacuated axial funnel, and a highly-relativistic, Poynting-flux–dominated jet within the funnel. This jet is accelerated and collimated by magnetic and gas pressure forces in the inner torus and the surrounding corona. A critical discussion of the physical origin of the relativistic Poynting jet formed in these numerical simulations has been conducted by Punsly , who points out the importance of incorporating resistive MHD reconnection in the modelling to gain further insight.
Axisymmetric models of jet formation allow one to extend the computational time of the simulations at an affordable cost. This is the case for the simulations reported by , which extended the earlier work of  by studying jets from a collapsar GRB model (see below) until and out to radial distances of . For this study a new version of the GRMHD code HARM was used. It was found that, at such large spatial and temporal scales, the Poynting-flux–dominated jet is accelerated by continuous mass-loading from the disk to reach bulk Lorentz factors of and maximum terminal Lorentz factors of and it collimates to narrow half-opening angles of . This remarkable result is in concordance with the values observed in jets associated with GRBs, AGN, or X-ray binary systems.
At stellar scales, relativistic thick disks (or tori) orbiting around stellar-mass black holes may form after the gravitational collapse of the core of a rotating massive star () or after a neutron-star–binary merger. Such an accreting torus–plus–black-hole system is the most likely candidate to be the mechanism powering gamma-ray bursts. Part of the large amounts of energy released by accretion is deposited in the low-density region along the rotation axis, where the heated gas expands in a jet-like fireball. In connection with GRBs, van Putten and Levinson  have considered the theoretical evolution of an accretion torus in suspended accretion. 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 nonideal MHD, incorporating radiative processes, are necessary to validate this picture.
Numerical simulations of relativistic jets propagating through progenitor stellar models of collapsars have been presented in [9, 435] (see also  for GRMHD simulations of jet formation in such context). The collapsar scenario, proposed by , is currently the most favored model for explaining long duration GRBs. The hydrodynamic simulations performed by  employ the three-dimensional code genesis  with a 2D spherical grid and equatorial-plane symmetry. The gravitational field of the black hole is described by the Schwarzschild metric, and the relativistic hydrodynamic equations are solved in the test-fluid approximation using a Godunov-type scheme. Aloy et al.  showed that the jet, initially formed by an ad hoc energy deposition of a few within a cone around the rotation axis, reaches the surface of the collapsar progenitor intact, with a maximum Lorentz factor of .
Similarly, the genesis code has also been used to perform hydrodynamic simulations of relativistic jets in the neutron-star–binary merger scenario. Figure 13 shows a snapshot at time s of an axisymmetric simulation of the launch and evolution of a relativistic outflow driven by energy deposition through neutrino-antineutrino annihilation in the vicinity of a black-hole–plus–accretion-torus system (model B1 of ). The left panel displays the logarithm of the Lorentz factor in the entire computational domain, while the right panel focuses on the innermost regions closer to the central black hole, depicting the logarithm of the rest-mass density. Since in this simulation the system is the remnant of a neutron-star–binary merger, it naturally provides a baryon-free channel along the rotation axis through which the outflow can reach ultrarelativistic velocities. Indeed, terminal Lorentz factors of about 1000 are attained for initial energy deposition rates of some erg/s per stereo-radian. As shown in Figure 13 the outflow is strongly collimated due to the presence of the accretion torus, which, in turn, makes the gamma-ray emission very weak outside of the polar cones. The calculations reported by  predict that about one out of a hundred neutron-star mergers produces a detectable GRB, as its ultrarelativistic jet happens to be sent along the line of sight 2. It is worth noting that the simulations performed by  have also been used by  to illustrate a new effect in relativistic hydrodynamics that can act as an efficient booster in jets and help explain the observed ultrarelativistic speeds. This effect, concealed in the subtleties of the relativistic Riemann problem, is purely hydrodynamic and occurs when large velocities tangential to a discontinuity are present in the flow, yielding Lorentz factors or larger in flows with moderate initial Lorentz factors.
The term “wind” or hydrodynamic accretion refers to the capture of matter by a moving object under the effects of the underlying gravitational field. The canonical astrophysical scenario in which matter is accreted in such a nonspherical way was suggested originally by Bondi, Hoyle, and Lyttleton [175, 58], who studied, using Newtonian gravity, the accretion on to a gravitating point mass moving with constant velocity through a nonrelativistic gas of uniform density. The matter flow inside the accretion radius, after being decelerated by a conical shock, is ultimately captured by the central object. Such a process applies to describe mass transfer and accretion in compact X-ray binaries, in particular in the case in which the donor (giant) star lies inside its Roche lobe and loses mass via a stellar wind. This wind impacts on the orbiting compact star forming a bow-shaped shock front around it. This process is also believed to occur during the common envelope phase in the evolution of a binary system.
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., [341, 42, 125] and references therein). The numerical investigations have revealed the formation of accretion disks and the appearance of nontrivial 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 nonrelativistic stars . More recently, Newtonian MHD studies have also investigated the spin-down of moving magnetars . For compact accretors such as neutron stars or black holes, the correct numerical modeling requires a general-relativistic (magneto-)hydrodynamic 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 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 performed by [133, 132]. This investigation differs from  in both the use of a conservative formulation for the hydrodynamic equations and the use of Godunov-type 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 [133, 132]. Nonaxisymmetric two-dimensional studies, restricted to the equatorial plane of the black hole, were discussed in , motivated by the nonstationary patterns found in Newtonian simulations (see, e.g., [42, 125]). The relativistic computations revealed that initially-asymptotic uniform flows always accrete onto the black 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 hydrodynamics 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 [135, 136] 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 14. 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 were the computation performed using the more common Boyer–Lindquist coordinates. The transformation induces a noticeable wrapping of the shock around the central black 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. We note that such “horizon-penetrating” coordinates (or variations thereof) are currently employed in most of the state-of-the-art codes developed to study accretion onto black holes [149, 201, 140, 396].
Semi-analytical studies of finite-sized collections of dust, shaped in the form of stars or shells and falling isotropically onto a black hole are available in the literature [282, 165, 350, 302, 318]. These early investigations approximate gravitational collapse by a dust shell of mass falling into a Schwarzschild black hole of mass . These studies have shown that for a fixed amount of infalling mass, the gravitational radiation efficiency is reduced compared to the point particle limit (), due to cancellations of the emission from distinct parts of the extended object.
Get Flash to see this player.
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. Similar studies based on a metric formalism have recently been discussed in [275, 276].
An example of the simulations of  appears in the movie of Figure 15. 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 parameterized 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 dynamic accretion phase, when the shell crosses the peak of the potential and is subsequently captured by the black hole. The gravitational wave signal coincides with the quasinormal ringing frequency of the Schwarzschild black hole, . 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; this “junk” radiation is also a common inevitable feature in the recent black-hole–binary simulations ).
One-dimensional numerical simulations of a self-gravitating perfect fluid accreting onto a black hole are presented in , where the effects of mass accretion during the gravitational wave emission from a black hole of growing mass are 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 nonspherical gravitational perturbations was captured by adding to the (characteristic) Einstein–perfect-fluid system, the evolution equation for a minimally-coupled massless scalar field. 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, , 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. To highlight the recent advances accomplished in numerical relativity it is worth it linking this result with the 3D simulations of black-hole formation performed by , which have provided information on the ring-down gravitational signal through which the dynamic black hole relaxes to its stationary state (see Section 5.1.2).
In a series of recent papers [434, 433, 269] it has also been shown that upon the introduction of perturbations, stable, relativistic, accretion tori manifest a long-term oscillatory behavior lasting for tens of orbital periods. The associated changes in the mass-quadrupole moment of such disks render these objects promising sources of high-frequency, detectable gravitational radiation for ground-based interferometers and advanced resonant bar detectors, particularly for Galactic systems and when the average disk density is close to nuclear matter density. If the disks are instead composed of low-density material stripped from the secondary star in Low-Mass X-Ray Binaries (LMXBs), their oscillations could help explain the high-frequency QPOs observed in the spectra of X-ray binaries. Indeed, such HFQPOs can be explained in terms of -mode oscillations of a small-size torus orbiting around a stellar-mass black hole .
The studies reported in the existing papers have considered both Schwarzschild and Kerr black holes, as well as constant and nonconstant (power-law) distributions of the specific angular momentum of the disks. The most recent investigation  has accounted for magnetic fields in the tori, which, in addition to satisfying a polytropic EOS and having a constant distribution of the specific angular momentum, were built with a nonzero toroidal magnetic field component, for which equilibrium configurations exist . A representative sample of initial models was built and the dependence of their dynamics on the strength of the magnetic field was investigated in long-term GRMHD evolutions in timescales of 100 orbital periods. Overall, the dynamics of the magnetized tori considered was found to be strikingly similar to that found in the purely hydrodynamic case by [434, 433]. It should be noted that, since the axisymmetric initial models do not have a poloidal magnetic field component, the MRI, which could potentially alter the conclusions, can not develop.
This work is licensed under a Creative Commons Attribution-Noncommercial-No Derivative Works 2.0 Germany License.