Approaching the black hole by numerical simulations

TL;DR

This paper reviews how numerical simulations are crucial for understanding the complex physics near astrophysical black holes, highlighting recent progress and different simulation approaches for a non-expert audience.

Contribution

It provides a broad overview of numerical methods used in black hole astrophysics, emphasizing recent advancements and differences in simulation techniques.

Findings

Significant progress in numerical code development and computational power.

Magnetohydrodynamic simulations are emphasized.

Differences between radiation transfer and merger simulations are highlighted.

Abstract

Black holes represent extreme conditions of physical laws. Being predicted about a century ago, they are now accepted as astrophysical reality by most of the scientific community. Only recently more direct evidence of their existence has been found - the detection of gravitational waves from black hole mergers and of the shadow of a supermassive black hole in the center of a galaxy. Astrophysical black holes are typically embedded in an active environment which is affected by the strong gravity. When the environmental material emits radiation, this radiation may carry imprints of the black hole that is hosting the radiation source. In order to understand the physical processes that take place in the close neighbourhood of astrophysical black holes, numerical methods and simulations play an essential role. This is simply because the dynamical evolution and the radiative interaction are far too complex in order to allow for an analytic solution of the physical equations. A huge progress has been made over the last decade(s) in the numerical code development, as well as in the computer power that is needed to run these codes. This review tries to summarize the basic questions and methods that are involved in the undertaking of investigating the astrophysics of black holes by numerical means. It is intended for a non-expert audience interested in an overview over this broad field. The review comes along without equations and thus without a detailed expert discussion of the underlying physical processes or numerical specifics. Instead, it intends to illustrate the richness of the field and to motivate for further reading. The review puts some emphasis on magnetohydrodynamic simulations, but also touches radiation transfer and merger simulations, in particular pointing out differences in these approaches.

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date \TitleApproaching the black hole by numerical simulations \AuthorChristian Fendt* \AuthorNamesChristian Fendt

\corresCorrespondence: fendt@mpia.de

1 Introduction

For many people - the scientists and the public – black holes are considered as the most fascinating objects in the universe. In principle, black holes simply consist of highly compressed and invisible mass, solely described by its mass and angular momentum, and are, thus, physically and chemically less rich compared to e.g. the Solar atmosphere. However, the compactness of the accumulated mass distorts the surrounding space time in such a way that the physical processes occurring close to a black hole seem to behave in contradiction to the every-day experience.

A prime mystery is (and may remain) the interior of a black hole – nobody knows how the innermost parts of a black hole are structured as we do not yet understand the physics at work under such extreme conditions. Below the black hole horizon the time and space coordinates change their meaning, which may impact our common understanding of causality. In order to understand the very center - if it exists - the extreme physical conditions would ask for a new theory of physics, that is quantum gravity, the combination of general relativity and quantum physics.

Astrophysical black holes come in various manifestations - from stellar mass black holes that can sometimes be observed as micro-quasars (Mirabel & Rodríguez, 1999) to supermassive black holes that reside in the nuclei of active galaxies and weigh up to billions of Solar masses (Fabian, 2012; Netzer, 2013). Intermediate mass black holes weighing some 1000s of Solar masses are discussed but not yet confirmed. Recently, by their emission of gravitational waves, a new flavor of black holes has been discovered rather unexpectedly – black holes with masses of 30-50 Solar masses for which the origin is not yet known.

Astrophysical black holes do not exist alone in space, but are surrounded by other material that is moving with high speed and is also radiating. As a consequence, we expect that certain features of this radiation may carry along characteristic imprints of the black hole from where it originated.

This review wants to touch on three basic questions of black hole physics that can be addressed by using numerical simulations. We want to understand (i) how does the environmental material dynamically behave when it comes close to the black hole, (ii) how does that material look like for a distant observer, and (iii) what could be learned about the black hole parameters by comparing theoretical modeling and observational data.

In order to answer these questions, the use of numerical simulations is essential, if not inevitable. Fortunately, the availability of large computer clusters for high-performance codes has vastly increased over the recent years. In general, these codes do solve the time-dependent physical equations considering general relativistic space time. As a result, the time-evolution of the dynamics of the matter and the electromagnetic field and also the radiation field is provided. Depending on the astrophysical context, different numerical approaches are undertaken.

One option is to concentrate on the dynamics of the material around a black hole. This is typically the accretion disk that is surrounding the supermassive black holes in the center of an active galactic nuclei or the stellar mass black holes in micro quasars. It is this accretion disk that makes the black hole shine (Shakura & Sunyaev, 1973). Depending on the black hole mass and the physical radiation process, we may observe X-ray emission or UV light, emission and absorption lines, a Doppler broadening of these lines, and we may even use these lines to re-construct the disk structure from the line signal (reverberation mapping, (Blandford & McKee, 1982)). The same sources may eject jets - collimated and magnetized beams of high velocity that are either launched in the surrounding accretion disk (Blandford & Payne, 1982) or by a rotating black hole itself (Blandford & Znajek, 1977). For both processes a strong magnetic field is essential. To simulate of the (magneto)hydrodynamics of disks and outflows as well as their radiation if one of the prime targets of numerical investigations. What are the mass accretion and ejection rates, the speed and energetics of the outflows, or the accretion disk luminosity?

For a black hole located in a intensive radiation field, we may expect to observe its ”shadow" against the background light as the black hole absorbs all the radiation that enters the horizon. For a long time – that means before the recent detection of gravitational waves – for many colleagues only the observation of the black hole shadow would have provided strong evidence for the actual existence of black holes. Numerical simulations of ray-tracing in a strong gravitational field are essential for the understanding of a (yet-to-come) observational signal of such a shadow. The idea is to compare the numerically constructed images of the area close to the black hole with the observation of nearby supermassive black hole shadow. These experiments are currently undertaken and are focusing on the black holes with the largest angular diameter that is the one in the Galactic center and the one in the active galaxy M87 (see https://eventhorizontelescope.org/; https://blackholecam.org/).

Even isolated black holes may be observed if the path of light from distant sources is penetrating deep enough into the gravitational potential of a black hole that is located just between the source and the observer and thus becomes characteristically deflected (gravitational lensing; see (Einstein, 1936; Refsdal, 1964)).

Close binary compact objects (neutron stars or black holes) will finally merge and thereby strongly distort the space-time around them. This distortion propagates through space as a gravitational wave which could then be detected on earth. As this is extremely difficult to measure, only recently such waves could be detected by the LIGO gravitational wave interferometer (see https://www.ligo.caltech.edu). In order to interpret the observed merger signal, the numerical simulation of a merging binary system and the subsequent gravitational wave is a fundamental tool which can deliver the progenitor masses and even the distance to the gravitational wave source.

2 Numerical approaches

In this review we will discuss the different numerical methods that were invented to investigate the close environment of black holes, considering their increasing numerical complexity. We distinguish the different approaches in three branches – these are

(i) simulations mainly treating the dynamics of the black hole environment on a fixed metric,

(ii) simulations that ray-trace photons that are affected by strong gravity, and

(iii) simulations that solve the time-dependent Einstein equations.

Of course, a combination of theses approaches are possible, but will be even more demanding.

In the following we first want to briefly compare and summarize the different methods, before we go on and discuss a few technical details and exemplary results in the remaining sections.

2.1 Dynamics on fixed space-time

It is considerably the most straight-forward approach to numerically evolve the environment of a black hole on a fixed space-time, thus applying a fixed metric. We may understand this approach may as an analogy to the non-relativistic treatment of gravity considering a (hydrodynamic) test mass in a Newtonian gravitational potential. The equations to be solved are the general relativistic (magneto)hydrodynamic (MHD) equations, that is the equation of motion defined by the time evolution of the energy-momentum-tensor on a given metric, fixed in time. The general approach is to choose a metric and coordinate system that are convenient for the black hole system under consideration and then solve the physical equations of interest.

A typical astrophysical prospect could be to investigate the evolution of an accretion disk or torus that is orbiting a black hole. In particular, one may investigate the stability of disks and tori, or their mass loss by accretion or the ejection of disk winds or jets. A topical motivation of this application is to learn about the feedback of supermassive black holes for cosmological galaxy formation simulations, which apply parameterized models of the energy output from galactic black holes.

The approach just described is equivalent to a non-relativistic numerical treatment that solves the (non-relativistic) MHD equations in the gravitational potential of e.g. a star. In both cases, the surrounding masses are test masses affected by gravity. As such we cannot tell from pure (magneto-)hydrodynamic simulations their astrophysical masses or the densities of the objects affected by a fixed gravity. Typically, the numerical approach gives an accretion rate in code units with no direct relation to an astrophysical value. However, if radiation transport is included (see below) the density scale actually matters as absorption and emission coefficients depend on the exact physical values (see e.g. (Noble et al., 2009; Jiang et al., 2016). Disk luminosities derived from such kind of simulations could be directly compared to observations and can thus constrain the disk density and the mass fluxes.

Numerical simulations of the close environment of the black hole were pioneered by Wilson and collaborators (Wilson, 1972, 1977; Hawley et al., 1984a, b; De Villiers & Hawley, 2002). After some years of silence - probably due to the lack of necessary computer power - the field of general relativistic time-dependent simulations flourished just after the turn of the century with a number of groups working on different numerical codes (Koide et al., 1999; Komissarov, 2001; De Villiers & Hawley, 2002, 2003a; Gammie et al., 2003; Noble et al., 2006; McKinney, 2006; Del Zanna et al., 2007).

More physics has been added, e.g. radiation processes or a physical resistivity. Note, however, that for a number of microscopic processes a consistent general relativistic description is missing, which may potentially lead to causality issues in the simulation (see for example (Hiscock & Lindblom, 1983; Denicol et al., 2008) for the case of viscosity or conduction).

An interesting feature to note is that GR-MHD simulations (Gammie et al., 2003; Noble et al., 2006) usually apply (modified) Kerr-Schild coordinates. In Kerr-Schild coordinates the coordinate singularity at the horizon that is present in the Schwarzschild metric or the Kerr metric is avoided by a suitable coordinate transformation. These codes actually continue to calculate the physical evolution of the variables inside the horizon. Naturally, physical processes taking place inside the horizon cannot affect the (numerical) world outside the horizon. The physical conditions defined by the horizon still apply - a neat boundary condition, that is difficult to accomplish if the horizon itself is the numerical boundary and if that boundary implies physical singularities.

A general difficulty that MHD codes have, are areas of low matter density. The numerical time step decreases with increasing Alfvén speed and may ”kill" the simulation. While this is a problem for non-relativistic simulations applying strong magnetic fields and low densities, in relativistic simulations the Alfvén speed is limited by the speed of light and thus the time-stepping is limited by the light-crossing time over a grid cell.

Low densities together with a strong magnetic field impose severe problems with certain variable transformations that are essential in relativistic MHD codes. Therefore, relativistic simulations consider a numerical floor model for the critical areas of low density, keeping the gas density and the internal energy above a ”floor" value, a chosen (local) threshold. By that, the codes can avoid to evolve into an unphysical parameter space, for example negative densities or pressures.

Typically, the choice of the numerical floor density typically affects the most relativistic regions of the simulation, for example a Poynting flux-dominated Blandford-Znajek jet (see below). For such regions an option may be to apply general relativistic force-free electrodynamic simulations (McKinney, 2006; McKinney & Narayan, 2007). However, while such simulation are able to treat highly magnetized areas very well, they hardly can capture the gas dynamics.

An essential parameter for any simulation using grid codes is the numerical resolution. High spatial resolution is desirable in order to the resolve the proper physical evolution. On the other hand, high resolution in a explicit code slows down the simulation because of the small numerical time step needed to ensure a causally correct propagation across the grid. A typical astrophysical context here is the evolution of hydrodynamical instabilities and turbulence, for example the magneto-rotational instability (MRI, (Balbus & Hawley, 1991)). While the simulations have made great progress towards capturing the MRI and subsequent processes such as angular momentum transport or a direct dynamo, the simulations often fail to capture a significant dynamical range of the turbulent cascade and thus the energy dissipation on small scales. The is of course a problem that is not restricted to general relativity.

2.2 Dynamics of interacting compact objects

Compared to simulations performed on a static metric as discussed above, a more general approach is to solve Einstein equations numerically and time-dependent. This procedure is obviously needed when the mass distribution in the simulation domain changes substantially in time, leading to a similar change in the metric. A typical application is the compact object merger. No fixed metric can be chosen initially - the metric is time-dependent and results from the evolution of the mass distribution. Further, the change in the metric that is propagated away in space, is just known as the gravitational wave.

As for the approach with the fixed metric discussed above, also this approach has non-relativistic equivalent. The example may be a typical star formation simulation of a molecular cloud collapsing under self-gravity. In order to evolve the simulation, after each time step the gravitational potential of the mass distribution has to be integrated, before in the next time step the motion of masses is calculated using the updated gravitational potential.

Pioneering work on the theoretical framework for numerical general relativity considering e.g. the mass formula for black holes, a systematic approach to solve Einstein’s equations on dynamical space-times, or even a definition of ”numerical relativity" has been published during the 1970s (Smarr, 1973, 1976; Smarr & York, 1978; Eardley & Smarr, 1979).

It took some years until seminal development of the first 3D numerical codes designed to solve the Einstein equations for general vacuum space times succeeded providing the evolution for the first 3D black hole spacetime (Anninos et al., 1995). A major difficulty here that had to overcome is to avoid the singularity inside the horizon or to treat the boundary conditions properly for spacetimes. Another feature, the so-called apparent horizon becomes important - a kind of spatial boundary that at a certain time effectively acts as an event horizon. In a 3D simulations both horizons are time-dependent and may fluctuate as the simulation progresses. The difficulty is how to determine the mass inside them?

The problem of dealing with the initial data in case of rotating black holes or distorted non-rotating black holes was treated by (Brandt & Seidel, 1996). Constructing the initial data for a simulation is a complex task involving to solve Einsteins equations on chosen initial spatial coordinates that than has to be embedded in the 4-dimensional space time (see e.g. (Varma et al., 2018)). Another seminal step was to successfully treat the motion of a black hole through a numerical grid. The difficulty here are the curvature singularities contained within the black hole(s). One way to deal with it is the so-called black hole-excision by that all parts of teh black hole interior are excluded from the computational domain and only the exterior region is evolved (Cook et al., 1998).

2.3 Radiation - propagation of light close to a black hole

The third approach we are discussing is to numerically propagate photons along space-time, thus to do relativistic radiation transport. What has discussed above was mainly related to the dynamical evolution of matter (and the magnetic field). Simulations of the dynamics deliver the time-dependent distribution of e.g. mass density and velocity. What can be potentially observed, is the radiation that is emitted from this material. Thus, radiation transport will be the essential step that allows to compare numerical simulations of hydrodynamics with the observed radiation pattern.

The path of light that is passing by a black hole is governed by the strong gravity. The principles seem rather simple: the trajectory of light rays follows the null geodesics of the metric - the Schwarzschild metric for non-rotating black holes or the Kerr metric for rotating black holes. More complicated situations arise if the metric is time-dependent as for example in case of merging compact stars, or if the mass distribution that is emitting and absorbing the light is changing rapidly with time.

Figure 1 shows a schematic diagram of photon paths around a black hole. The photon trajectories can be bent around the black hole if a photon comes within the the co-called photon sphere. Figure 2 displays a 3-dimensional realization of a numerical ray-tracing procedure applied to a luminous accretion disk orbiting a Kerr black, as seen by a distant observer under different inclination angles. Obviously, a rich radiative structure emerges including parts of the disk that are located behind the black hole and become visible by gravity.

The same principle applies and is observationally confirmed also for weaker gravitational fields. The famous historical example is the deflection of light by stars predicted by Einstein (Einstein, 1911) and its confirmation observed during a Solar eclipse (Dyson et al., 1920). A modern application is the observation of gravitational lensing (Einstein, 1936; Refsdal, 1964) for a variety of sources. Examples are the discovery of planets by micro-lensing events in the Milky Way (Beaulieu et al., 2006), the observation of multiple images of distant quasars (Huchra et al., 1985), or the measure of shear and flexion in the large scale structure of the universe (see e.g. (Blandford & Narayan, 1992)).

Apart from the general relativistic effects discussed above, special relativistic Doppler beaming for moving sources (a jet or a disk) plays an essential role for the observed signal. The blue-shifted signal is Doppler-boosted, thus brighter, while the red-shifted is de-boosted for the observer.

In addition to the pure ray-tracing of light rays trough space-time, radiative transfer can be applied, dealing with absorption and emissivity along the path of light. This problem is substantially more complex as dealing with material opacities or scattering and can be solved in approximations. However, in general, the material properties such as densities, temperatures etc as well as the light sources can be provided by either a kinematic prescription or by post-processing of dynamical data resulting from an (magneto-)hydrodynamic simulation (Schnittman et al., 2006; Dexter et al., 2010; Vincent et al., 2011; Chan et al., 2013; Zhu et al., 2015; Dexter, 2016; Pu et al., 2016b) (see below).

3 Dynamics & feedback - (magneto)hydrodynamic simulations

Black holes are able to feedback to their environment – implying that mass, energy or angular momentum is taken away from the accompanying components or the black hole itself and is deposited further away from the black hole.

A typical context is the nuclear activity of galaxy such as quasars with extremely high luminosity that arises from the black hole accretion disk, or radio galaxies that may drive Mpc-long relativistic jets. These feedback processes are thought to be essential for cosmological structure formation and the formation of galaxies (Fabian, 2012). On a lower energy scales similar processes work in Galactic micro-quasars that are powered by stellar-mass black holes (Mirabel & Rodríguez, 1999). In a similar way, central black holes are believed to be the very origin of the explosive Gamma-ray bursts observed on cosmological distances (Kumar & Zhang, 2015).

A fundamental question that can be tackled by a numerical treatment is the mass accretion to the black hole. The physics of the accretion process determines the black hole luminosity and also the growth of the black hole mass (and angular momentum). Quite a number of general relativistic simulations of black hole accretion flows have been carried out over the last decades, both in the hydrodynamic limit (e.g. (Wilson, 1972; Hawley et al., 1984a, b; Hawley, 1991; Fragile & Anninos, 2005)), but also applying MHD, and some of them even considering some sort of radiative transport (e.g. (Gammie et al., 2003; Fragile et al., 2017; Noble et al., 2009; Sadowski et al., 2014; Morales Teixeira et al., 2014)).

The physical origin of magnetized jet outflows has been proposed by seminal papers by (Blandford & Znajek, 1977) and (Blandford & Payne, 1982). The first mechanism - known as Blandford-Znajek mechanism - works only for rotating black holes. Here, the frame dragging of space-time leads to the induction of a toroidal magnetic field from the poloidal magnetic field threading the black hole ergosphere. As a result, a strong Poynting flux can be generated, leading to a magnetically dominated jet tower that can be filled by lepton pairs that are produced by the strong radiation field of the accretion disk.

The second mechanism - known as Blandford-Payne mechanism - deals with the acceleration and collimation of disk winds by magneto-centrifugal forces. Many non-relativistic MHD simulations on this topic have been published (see e.g. (Ouyed & Pudritz, 1997)), but also a few relativistic simulations that show that this effect may also work in relativity (Porth & Fendt, 2010; Porth et al., 2011). However, it seems to be difficult to gain very high Lorentz factors for disk winds (Komissarov et al., 2009).

The launching question - the question how the disk wind material is lifted into the outflow - is more complex and requires to treat the accretion disk evolution simultaneously. Analytical (non-relativistic) theory has demonstrated that magnetic diffusivity (respectively, resistivity) is essential for launching (Ferreira, 1997) – that implies to re-distribute material between different magnetic flux surfaces and by that to feed the outflow with disk material.

Simulations of wind or jet launching from thick hot disks have been performed applying non-relativistic codes (Yuan et al., 2012, 2015), see also Yuan & Narayan (2014). Similarly, launching simulations of jets from thin disks have been published (see for example (Zanni et al., 2007; Sheikhnezami et al., 2012; Zhu & Stone, 2018) pointing again out the importance of the magnetic diffusivity for mass loading. From non-relativistic simulations we also know that the disk magnetization plays a leading role for governing the outflow properties (Stepanovs & Fendt, 2016).

General relativistic MHD simulations were pioneered by (Wilson, 1977), while first results were published by (Koide et al., 1999). Soon, further numerical codes became available boosting the number of publications in this field (De Villiers & Hawley, 2002; Gammie et al., 2003; De Villiers & Hawley, 2003a, b; McKinney & Gammie, 2004; De Villiers et al., 2005; Noble et al., 2006; McKinney, 2006; Noble et al., 2009; Narayan et al., 2012, 2012; Porth et al., 2017). Over the last decade constant and rapid progress in this field has been made.

Most often for the simulation a stationary hydrodynamic torus (Fishbone & Moncrief, 1976) is assumed as initial condition from which a (thin) disk may emerge. However, it looks like that some leading parameters of the flow evolution seem to be determined by this (probably not realistic) initial setup and are not "forgotten" during the simulation, so its use is maybe limited (see e.g. (Penna et al., 2013)). The torus-disk structure evolves under the magneto-rotational instability (Balbus & Hawley, 1991), leading to accretion of material towards the black hole and also to weak disk winds. Advection of magnetic flux along the disk towards the black hole leads to the accumulation of a strong axial poloidal flux that ”threads" the black hole.

The feasibility of the so-called Blandford-Znajek mechanism was numerically demonstrated by (Komissarov, 2001) by a pure electromagnetic simulations and later also by a number of MHD simulations in general relativity (see for example (Koide, 2003; Gammie et al., 2003; McKinney & Gammie, 2004; De Villiers et al., 2005; McKinney, 2005; Komissarov, 2005)).

Figure 3 shows the result from a 3D GR-MHD simulation investigating the stability of jets launched by a rotating black hole (McKinney & Blandford, 2009). The authors find jets with Lorentz factor of about 10 emerging along the black hole rotational axis with opening angles of about 10 degrees. The jet reaches about 1000 gravitational radii distance from the black hole without significant distortion by a 3D instability. These numbers are constrained by the numerical capability of the available codes. So far much higher relativistic jets cannot be modeled, meaning that the jet mass fluxes and the jet speed are limited by the floor model (see above). The jet stability is essentially affected by the numerical resolution and also by the conditions in the environment.

Of particular interest is the so-called plunging region between the inner disk radius (located approximately at the innermost stable orbit, ISCO) and the horizon. Of course, also the generation of jets driven by black hole rotation is governed by this area. This region is in particular sensitive to the frame dragging of a spinning black hole and also determines the energy conversion from the potential and orbital energy of the disk material. Examples for simulations emphasizing on the plunging region are (Penna et al., 2010) who investigate how the electromagnetic stresses within the ISCO depend on the black hole spin or the thickness of the accretion disk; or (Tchekhovskoy et al., 2011) who consider jet formation in from magnetically arrested disks. In (Tchekhovskoy et al., 2010) the interrelation between black hole spin and the jet power output applying their results to the long-standing problem of the radio loud/quiet dichotomy of AGN. In a similar approach (Tchekhovskoy & McKinney, 2012) find that black holes that rotate prograde compared to the disk may trap more magnetic flux then retrograde black holes leading to a higher jet efficiency.

The time-evolution of thin disks has been studied by (Noble et al., 2010) in particular considering the evolution of and the electromagnetic stress at the inner disk radius in the Schwarzschild case. This is important in order to understand the efficiency of the accretion process and the time scale for black hole growth. Tilted thin disks and the potential onset of precession due to the Bardeen-Petersen effect were investigated by (Morales Teixeira et al., 2014) in 3D simulations lasting up to 13,000 time units $t_{\rm g}=R_{\rm g}/c$ that are measured in light crossing time over the gravitational radius $R_{\rm g}$. Precession of the central black hole or the inner accretion disk may be seen as a time-dependent alignment of the jet outflow on much larger (pc-kpc) scale.

Also resistive GR-MHD have been applied, in particular to investigate the magnetic field generation by a dynamo process (Bucciantini & Del Zanna, 2013; Bugli et al., 2014), to compare the energy output from the disk wind and the Blandford-Znajek jet (Qian et al., 2017, 2018), or a stellar magnetized collapse to a black hole (Palenzuela et al., 2009; Dionysopoulou et al., 2013).

4 Radiation from simulated gas dynamics

From GR-MHD simulations of the black hole environment the time evolution of the gas dynamics and the magnetic field can be derived. However, what is observed in reality is the radiation that is emitted from this gas - observers obtain spectra or radiation images. In order to compare simulation results with observations, it is therefore essential to couple radiation transfer to the hydrodynamics of the simulation.

Most of the emission from the black hole source is expected to come from the innermost accretion disk. However, the light that is emitted in the plunging region during the final infall may actually carry most of potential (general) relativistic imprints, in particular imprints from the metric that is defined by the central black hole. By simulating these signals and comparing them to observations we may gain direct information about the black hole characteristics such as e.g. the black hole spin (see e.g. (Zhu et al., 2012)).

Combining MHD dynamics and radiation transfer is computationally highly expensive and yet numerically and physically still limited. A number of issues have to be considered. The method of choice would be to do radiation transfer modeling along with the hydrodynamics that is to follow the absorption and emission of photons by the gas along with the gas dynamics.

It is well known that the radiation transfer equations need to be closed. A closure relation is needed as radiation transfer is typically treated considering the moments such as radiation energy, radiation flux and radiation pressure. The local geometry of the radiation field needs to be specified by a radiation model (for example LTE) which closes the set of radiation transfer equations. A few astrophysical options for the closure problem have been established. Each method may have its advantages and disadvantages for a certain astrophysical problem. One option is the flux-limited diffusion, applied for a short mean free path of photons, thus for large optical depth close the the case of LTE. The application of the diffusion approximation is often coarse, however, radiative transfer would be very expensive. A well known closure is the Eddington approximation, typically applied in a plane-parallel geometry like a geometrically thin stellar photosphere. However, this is probably not a good choice for a turbulent disk-jet structure around a black hole.

The method of choice applied in most GR-MHD simulations is the so-called M1 closure that evolves the radiation flux independently from the radiation energy. The basic assumption is that the radiation field is locally a dipole aligned with the radiation flux. The M1 closure has been introduced in GR hydrodynamics by (Sadowski et al., 2013) and for special relativistic magnetized fluids by (Takahashi & Ohsuga, 2013). In contrast to the diffusion limit, M1 could handle "shadows". In (Sadowski et al., 2014) this approach was applied to investigate the impact of opacity on the accretion flow structure and evolution in case of super-Eddington luminosities. In a similar way, (McKinney et al., 2014) could model super-Eddington disks with accretion rates 100-200 times the Eddington accretion rate.

Another closure that has recently be proposed is the Monte-Carlo closure for general relativistic radiation hydrodynamics (Foucart, 2018). This has been applied in particular also to investigate the neutrino transport. Neutrino transport has been proven to be essential for the understanding the neutrino emission they radiate along with the gravitational waves when forming a black hole (Foucart et al., 2016).

A computationally less expensive approach that avoids the costly radiative transfer modeling is to consider a cooling function in the GR-MHD simulation. This method has been introduced in GR-MHD by (Noble et al., 2009) to calculate the radiative efficiency of the black hole disk. The disk radiation has been ray-traced to compute the flux that eventually is received by the observer. In (Noble et al., 2007) the mm and sub-mm flux for Sgr A* has been derived. Spectra and radiation maps were derived applying various combinations of the black hole spin parameter and the disk accretion rate in order to match the observed 1-mm flux of 4 Jy. In (Dexter et al., 2010) the authors compare the numerically derived disk luminosity to observations of Sgr A* deriving an accretion rate of $5^{+15}_{-2}\times 10^{-9}{M}_{\odot}\,{\rm yr}^{-1}$ as well as temperatures and inclination angle. Similarly, (Chan et al., 2015) derive spectra from GR-MHD modeling and compare them to the properties of Sgr A*. In particular, they show that by combining this information with the emission images certain models (those with strong funnel emission) can be ruled out. Heating and cooling of electrons have been included in recent GR-MHD simulations by (Ressler et al., 2015) and (Ryan et al., 2017).

A widely used approach to construct black hole radiation features is post-processing of the simulated MHD data. This is an option for example when considering the non-thermal radiation or optically thin cases. For non-thermal emission, however, the energy distribution of the emitting high energy particles cannot be derived self-consistently within the MHD approach, and must thus be prescribed. Depending on whether the particle energies are interrelated with the gas pressure, the gas density or the magnetic field energy, the resulting radiation maps may vary greatly (see e.g. (Porth et al., 2011)).

Post-processing of radiation features is always problematic when dealing with optically thick media. The optical depth of the material and also its opacity is not always known and post-processing remains an approximation. This is in particular an issue when dealing with thermal emission. Then, ideally, radiation transport should be done along with the MHD simulation (as discussed above), as the radiation field may affect the gas temperature that in turn governs the radiation field.

We mention another issue for this method that is linked to the coordinate system and time stepping generally used by GR-MHD codes. GR-MHD simulations are typically performed in (modified) Kerr-Schild coordinates, thus providing time slices - snapshots - in Kerr-Schild time. Observations - ray tracing - is typically done in the observers frame, for example in Boyer-Lindquist coordinates. The difficulty is that in every snapshot in Kerr-Schild coordinates - representing a 2D or 3D distribution of e.g. gas density or magnetic field - each pixel represents a different time in Boyer-Lindquist coordinates. Thus, performing a ray-tracing along the photon path that is smooth in time seems impossible.

In addition to that, ray-tracing across material that moving relativistically needs to consider the retarded times. The problem is that while the signal (emitted at a certain retarded time) is propagating with time through the gas, ”new" material is moving into that path of the photon (this motion is actually) calculated by the GR-MHD simulation. It is therefore be difficult to interpret post-processed radiation maps for highly time-variable sources such as jets or accretion disk coronae.

So far, ray tracing in GR-MHD-simulated gas dynamics circumvents this problem by assuming the so-called fast-light approximation that postulates that the light travels faster than the dynamics of the system evolves. Another option to get (post-processed) mock images of simulation data is to stack simulation snapshots and ray-trace along the time-averaged gas distribution. While this is a valid approach in order to obtain time-averaged radiation features such as a typical spectrum of the source, the method is obviously not applicable for highly time-variable phenomena.

In reality, post-processing of radiation transfer may not be fully consistent with the dynamics of the relativistic plasma, but at the moment, there seem to be no numerical tools at hand that can circumvent this problem in respect of the available CPU resources. Post-processing radiative transfer is widely applied (Schnittman et al., 2006; Dolence et al., 2009; Dexter et al., 2010; Chan et al., 2013; Zhu et al., 2015; Dexter, 2016; Pu et al., 2016a, b, 2017) and has recently been extended also for comptonization (Narayan et al., 2016).

As a specific example, we mention (Schnittman et al., 2006) who derive light curves of black hole accretion disks from radiative transfer based on GR-MHD simulations with the application on quasi-periodic oscillations of X-ray binaries. In (Gold et al., 2017) a comparison study of black hole shadow images for Sgr A* is presented demonstrating the impact of e.g. the inclination angle and the temperature prescription, scattering, or the stacking of images on the visibility of the shadow. This is shown in Figure 4 demonstrating the complexity of the parameter space and thus the difficulty to obtain ”true" radiation maps. The authors conclude that all these effects together may make it difficult to unambiguously detect the black hole shadow or to extract information about the space-time from its size and shape, before the plasma physics and dynamical properties of disk and jet are constrained (see also (Psaltis et al., 2015)).

Nevertheless, mock observations from simulated MHD simulations of the black hole environment are extremely important for understanding the observations of the "shadow" of nearby supermassive black holes (Sgr A∗ and M87) as proposed by the Event Horizon Telescope (EHT) initiative (see https://eventhorizontelescope.org/. Radiation transfer modeling is essential in order to disentangle the various processes at work at these spatial scales, such as disk and outflow dynamics, light bending, or absorption or scattering (Chan et al., 2015; Pu et al., 2017).

5 Compact object mergers

”For a brief moment, a binary black hole merger can be the most powerful astrophysical event in the visible Universe" (Keitel et al., 2017). The numerical treatment of compact object mergers is considerably the most difficult numerical task to date. While the dynamical evolution of the material around a single black hole can be treated on a fixed space-time, the evolution of (close) binary black holes is substantially more complex. As the metric changes with time, for this problem the Einstein equations needs to be solved numerically in a way that (i) allows for a well-posed initial-boundary-value problem, and that (ii) singularites inherent in the black hole spacetimes can be circumvented (see e.g. (Centrella et al., 2010; Sperhake, 2015)). The same challenges appear for mergers of binary neutron stars into a black hole.

In comparison to the simulations discussed in the last sections, most simulations of black hole mergers typically consider only gravity. Therefore, the results depend only on few parameters (black hole masses and spins), and do not rely on further assumptions on physical quantities as e.g. the accretion disk density or the magnetic field. Fundamental breakthrough in the simulation of binary black holes came around 2005 with seminal papers considering the three basic phases of a black hole merger - the orbit (only up to a single orbit in these early years), the merger, and the following ringdown (Brügmann et al., 2004; Pretorius, 2005; Baker et al., 2006; Campanelli et al., 2006).

Deriving the gravitational wave signal from a merger in its final stage requires to evolve the 3D Einstein equations. Further, the merger simulation has to run sufficiently long (Alcubierre et al., 2001; Diener et al., 2006). An early example of the gravitational wave signal and the apparent horizons of the merging black hole as derived by (Baker et al., 2006) is shown in Figure 5. These early modeling efforts succeeded with using excision to exclude singular regions within the horizons (see above). For the extraction of the wave signal from the numerical simulation of the merger – thus to extract the radiative part of the numerical solution of the Einstein equations – a number of methods have been developed. For a detailed discussion of the problems involved we refer to (Bishop & Rezzolla, 2016).

The field of compact merger simulations has evolved rapidly, benefiting also from the advance in CPU power. Some codes such as the ”Einstein Toolkit" have even been made publicly available (Löffler et al., 2012; Mösta et al., 2014). Hundreds of binary black hole simulations have been performed by now, allowing to do a proper statistics of gravitational wave forms and the gravitational wave luminosity or on the final black hole properties based on different initial mass ratios or black hole spins (see e.g. (Jani et al., 2016; Healy et al., 2017; Keitel et al., 2017). A recent, colorful 3D representation of the merging black holes emitting gravitational waves is shown in Figure 6.

As considerably one of the most fascinating events in recent science, the long lasting efforts in both theoretical predictions and detector experiments were culminating in 2015 when LIGO interferometer eventually detected a gravitational wave signal (Abbott et al., 2016a, b) that was perfectly matching the theoretically predicted waveform models for a black hole merger (Hannam et al., 2014; Taracchini et al., 2014) . The detection clearly indicates on a black hole merger of masses of $36^{+5}_{-4}{M}_{\odot}$ and $29^{+4}_{-4}{M}_{\odot}$ to a final black hole of mass $62^{+4}_{-4}{M}_{\odot}$ at a luminosity distance of $410^{+160}_{-180}$ Mpc, thereby radiating away gravitational waves111 Interested readers can actually redo the modeling of the GW 150914 binary black hole merger gravitational wave signal by using the ”Einstein Toolkit”, see https://einsteintoolkit.org/gallery/bbh/. of energy of $3.0^{+0.5}_{-0.5}{M}_{\odot}c^{2}$. Besides the fascinating confirmation of the existence of gravitational waves, the data also revealed a new and unexpected mass range for black holes. As of early 2019, gravitational wave detections has reached the number of 11 (for the actual status see https://www.ligo.caltech.edu/page/detection-companion-papers).

Another class of merger simulations consider mergers of neutron star binaries. This field has been pioneered by Shibata & Uryu (Shibata & Uryu, 2000; Shibata et al., 2003, 2005). Binary neutron star mergers are related to black hole numerics as these mergers do form a black hole. As a matter of fact, these kind of mergers may provide an observational signal in addition to the pure gravitational wave.

A particular difficulty in neutron star merger simulations compared to black hole merger simulations is that one has to considermaterial with a certain equation of state. However, despite a huge observational and theoretical progress, the equation of state for neutron stars is still one of the major unsolved problems in astrophysics (see e.g. (Steiner et al., 2013; Özel & Freire, 2016)).

Merger simulations investigating different nuclear equations of state have been performed (Hotokezaka et al., 2011), indicating for example that different kinds of final merger states can be realized that are strongly dependent on the equation of state applied. Simulations of this kind also allow to tackle the evolution of a massive torus surrounding the remnant black hole. In turn, the gravitational wave signal from binary neutron star mergers could be used to probe the extreme-density matter in these stars (Radice et al., 2017).

Examples for further specifics of neutron star mergers in comparison to black hole merges are tidal effects of the stellar mass distribution (Bernuzzi et al., 2015), a possible mass ejection and radiation (Sekiguchi et al., 2016), or neutrino emission (Foucart et al., 2016). Also, neutron star mergers have long been discussed as sources of short gamma-ray bursts, therefore simulations of such a scenario are essential to confirm this hypothesis (Bauswein & Stergioulas, 2015; Paschalidis, 2017). A particularly exciting example of this kind of simulations is the modeling (Shibata et al., 2017) of the gravitational wave signal GW170817 that was connected to the simultaneous detection of a gamma ray burst (Abbott et al., 2017).

A remaining step to do would be the combination of the two approaches discussed above, namely the combination of a merger simulation with the magnetohydrodynamics of the binary black hole environment, thus to combine the numerical approach discussed in this section (pure gravity) with the dynamical MHD simulations discussed before. This has been achieved by e.g. (Farris et al., 2012), who, for the first time, run a fully general relativistic MHD simulations of an equal-mass black-hole binary merger in a magnetized, circum-(black hole)-binary accretion disk (see also (Gold et al., 2014)). Such simulations can predict for example the change of accretion rate during the merger event. Essentially, they may also provide - in addition to the binary merger gravitational wave signal - the contemporary electromagnetic signal during the merger event.

6 Summary & outlook

This brief (and therefore incomplete) review intended to summarize the rapid development and the various applications of numerical methods in general relativity over the past decades.

The numerical approaches involving black holes can roughly be distinguished in simulations (i) treating the dynamics of the black hole environment on a fixed metric; or (ii) solving for time-dependent solutions of the Einstein equations that allow to to follow the merger of compact objects into a remnant black hole and the gravitational waves that are emitted; or (iii) following the photons that are affected by strong gravity, and thus allowing the see how black holes ”look alike". Even more demanding are combinations of these approaches such as obtaining mock observations of simulated dynamical data involving magnetohydrodynamics and radiation transfer.

The comfortable situation today is that quite a number of codes have been developed for these different approaches. In addition, also the computational resources to operate them have vastly improved. An essential point is to compare the existing codes and the robustness of their results concerning e.g. the accretion towards the black hole. This task is currently undertaken as a collaboration of nine groups and first results are published recently (Porth et al., 2019).

Besides fundamental questions of theoretical physics such as e.g. testing general relativity in the strong field limit, from the astrophysical perspective the most intriguing problems are maybe (i) the formation and early growth of (distant) supermassive black holes - involving most probably both, the merger scenario and also disk accretion; (ii) the mystery of gamma-ray bursts, or (iii) the inner structure (the equation of state) and the evolution of compact stars.

Essential observational input is expected from future gravitational wave signals measured by LIGO, as well as from highly resolved radio observations that are going to image the shadow of nearby supermassive black holes. Striking first images of a black hole shadow namely that of the nearby galaxy M87 and corresponding simulations in GR-MHD that are essential for the interpretation of the data have been published by the EHT collaboration very recently Event Horizon Telescope Collaboration (2019, 2019).

With all the advance in the observational data that will become available, the numerical modeling will be central for the understanding of these findings.

\conflictsofinterest

The author declares no conflict of interest.

\abbreviations

The following abbreviations are used in this manuscript:

[TABLE]

\reftitle

References