## Abstract

Einstein’s theory of general relativity affords an enormously successful description of gravity. The theory encodes the gravitational interaction in the metric, a tensor field on spacetime that satisfies partial differential equations known as the Einstein equations. This review introduces some of the fundamental concepts of numerical relativity—solving the Einstein equations on the computer—in simple terms. As a primary example, we consider the solution of the general relativistic two-body problem, which features prominently in the new field of gravitational wave astronomy.

The basic equations of general relativity are the Einstein equations, first published in 1915 (*1*). However, even today there are large gaps in our understanding of the physics implied by the Einstein equations. Stated in general terms, a major goal of research in general relativity is to solve the Einstein equations for the physical situations of interest. Fundamental analytic solutions of the Einstein equations include the flat Minkowski spacetime known from special relativity, the Schwarzschild and Kerr spacetimes describing single black holes, and the simple Big Bang cosmologies. Also predicted by general relativity are gravitational waves, which for weak fields can be obtained as analytic solutions of the linearized Einstein equations. However, the few known analytic solutions describe only very special situations, and approximation methods fail in the regime where the nonlinear, strong-field effects of relativity play a crucial role. If we are interested in the truly relativistic regime, we must turn to computer simulations to obtain numerical solutions to the full Einstein equations.

Solving the full Einstein equations on the computer is the subject of numerical relativity, which could also be called computational general relativity. Computers also play a role in algebraic computations and in approximation schemes, and such calculations are important topics in numerical relativity. But the distinguishing feature of numerical relativity is that, in principle, the Einstein equations in full generality can and must be solved numerically.

Numerical relativity spans a large range of different topics including mathematical general relativity, astrophysics, numerical methods for partial differential equations, computer programming, and simulation science. Current research in numerical relativity is in a transition from a self-contained topic in theoretical physics to a physical theory with numerous connections to observational astronomy (*2*, *3*). Gravitational wave astronomy holds much promise for the future, as recognized by the 2017 Nobel Prize in Physics, and numerical relativity is providing key theoretical predictions and analysis tools for the ongoing gravitational wave observations (*4*).

## The general relativistic two-body problem

As a primary application of numerical relativity, we consider the gravitational two-body problem. The two-body problem in Newtonian gravitational physics can be formulated for two point masses moving in their mutual gravitational field. A particular solution of the Newtonian two-body problem is a Keplerian elliptical orbit. However, in Einsteinian gravity, such orbital motion generates gravitational waves that carry away energy and momentum. Binary orbits therefore decay, and the motion of the two bodies follows an inward spiral that eventually terminates with the collision and merger of the two objects. In most astrophysical situations, the energy loss due to the emission of gravitational waves is so small that a binary orbit decays only on time scales of millions or billions of years. However, for compact objects such as neutron stars or black holes in very tight binaries, general relativistic effects such as gravitational wave emission play a major role (*5*).

Research in this field seeks to provide a theoretical framework for the physics of binary black holes, neutron stars, and gravitational waves. Such an endeavor must rely on numerical simulations in general relativity and general relativistic hydrodynamics. But a reasonably complete framework still requires substantial progress in numerical relativity and related fields. Currently there are serious limitations in our ability to model the entire range of relevant physics, from the nuclear physics of neutron star matter to the large-scale, strong-gravity effects encountered in binary neutron star mergers (*6*). The different dynamical phases of the binary evolution—known as the inspiral, the merger, and the evolution of the remnant—are accompanied by characteristic gravitational wave signatures (Fig. 1). For binaries involving at least one neutron star, depending on the specifics of the system, key features include the disruption of the star(s) before merger, the formation of a hypermassive neutron star, the prompt or delayed collapse to a black hole, the dynamics of the accretion torus plus the central merged object, and the creation of unbound material, the ejecta. Before discussing simulations of these systems, we introduce the mathematical foundation of numerical relativity.

## Mathematical foundation

Combining space and time into spacetime can be considered a triumph of human thought, allowing us to perceive the true nature of relativistic and gravitational physics (*7*). However, this does not mean that we cannot or need not consider space and time separately. Somewhat ironically, after working hard to unify space and time, the mathematical setup of numerical relativity starts by splitting spacetime again into space and time and by making gauge (coordinate) choices (*8*) in order to reformulate the Einstein equations as a well-posed mathematical problem.

General relativity is the theory of a metric tensor on a four-dimensional manifold, plus matter described by additional tensor fields. A manifold ℳ endowed with a metric *g _{ab}* is called a spacetime (ℳ,

*g*). The metric measures lengths, here in four dimensions. The infinitesimal line element(1)provides a generalization of the Pythagorean theorem. Repeated indices are summed over, following the Einstein summation convention. The metric is symmetric (

_{ab}*g*=

_{ab}*g*), has a Lorentz signature of – + + +, and there exists an inverse metric

_{ba}*g*defined by , where is the identity matrix. A special example is the Minkowski line element

^{ab}*ds*

^{2}= –

*c*

^{2}

*dt*

^{2}+

*dx*

^{2}+

*dy*

^{2}+

*dz*

^{2}, where

*c*is the speed of light,

*t*is the time coordinate, and

*x*,

*y*, and

*z*are spatial coordinates. The components of the Minkowski metric are constants, but in general

*g*is a field with nonconstant components.

_{ab}The field equations of general relativity are the Einstein equations,(2)where *G _{ab}* is the Einstein tensor, which depends on the metric and its first and second derivatives, and

*T*is the stress-energy tensor constructed from the matter fields Φ and in general also from the metric. For numerical implementations, the first step is to write the Einstein equations as a well-posed system of partial differential equations (PDEs) for the metric. Equation 2 represents 10 coupled, nonlinear PDEs for the 10 independent components of the metric, but without further adjustments these equations are in no known sense hyperbolic (i.e., well-posed as an initial value problem).

_{ab}The differential operator acting on the metric in the Einstein equations is given by the Ricci tensor,(3)The first term by itself, *g ^{cd}*∂

*∂*

_{c}*, which is often denoted*

_{d}g_{ab}*□*

^{g}*g*(where □ is the d’Alembert operator), has the form of the principal part of a simple hyperbolic wave equation, but note that the metric appears in two places: as the wave field

_{ab}*g*and as the inverse metric

_{ab}*g*in the wave operator. The other second-derivative terms are not standard wave operators. The best we can say about the complete principal part in Eq. 3 is that it is quasi-linear in the metric; that is, it is linear in the highest-order derivatives but with coefficients that depend (nonlinearly) on the variable itself. The lower-order terms are quite involved as well, with typical terms of the form

^{cd}*g*

^{–1}

*g*

^{–1}∂

*g*∂

*g*. Approaching the problem in this way makes it difficult to recognize that these equations are describing the simple geometric concept of curvature and that there is a time evolution being defined. Further, their well-posedness properties are quite unclear.

The so-called 3+1 decomposition—for example, in the form of Arnowitt, Deser, and Misner (ADM) (*8*)—assumes that the manifold (at least locally) allows a split into time and space, ℳ = *R* × Σ. Physics is then describable by time-dependent tensors on three-dimensional hypersurfaces Σ, which correspond to *t* = constant slices of ℳ, resulting in a “foliation” of spacetime in terms of three-dimensional spaces. Geometrically, we obtain a normal vector *n ^{a}* to Σ that allows the decomposition of tensors in directions normal and tangential to the hypersurfaces. These are the time-like and space-like directions, respectively. For concreteness, we can assume coordinates

*x*= (

^{a}*t*,

*x*) = (

^{i}*t*,

*x*,

*y*,

*z*) with a time coordinate

*x*

^{0}=

*t*and spatial coordinates

*x*, where

^{i}*i*= 1, 2, 3 and

*a*= 0, 1, 2, 3.

Decomposing the Einstein equations is accomplished by projecting *G _{ab}* and

*R*in the directions parallel and orthogonal to

_{ab}*n*. We discover that the differential operator Eq. 3 leads to two types of equations: (i) evolution equations containing time derivatives, and (ii) four constraint equations that are essentially elliptic equations, highlighting the indeterminate type of Eq. 3. The constraints are the Hamiltonian constraint and the momentum constraints. The latter are reminiscent of the Gauss law constraint of electrodynamics, where the divergence of the electric field gives the charge density.

^{a}Given the evolution and constraint equations as PDEs, we still must choose a spatial and temporal domain with boundary conditions. For problems in astrophysics such as the two-body problem of black holes and neutron stars, we consider isolated systems where at large distances gravitational fields become weak and spacetime becomes asymptotically flat (in contrast to typical cosmological models). Because gravity is universally attractive and long-range, it is not natural to restrict a system to a finite box, especially given that the goal is to compute waves traveling to infinity. Nonetheless, a typical configuration for numerical simulations is a finite-size spatial domain (e.g., a sphere) with boundary conditions at some finite radius that implement the proper fall-off of the fields and an outgoing-wave boundary (*9*).

Features unique to numerical relativity are various aspects of black hole spacetimes, in particular the causal structure associated with black hole event horizons and the possibility of spacetime singularities. This latter aspect can be viewed as the problem of specifying additional boundaries that represent black holes within the simulation domain.

## Building blocks of numerical relativity

To define a particular strategy to solve the Einstein equations, we consider the following building blocks that define the anatomy of a numerical relativity simulation, with a focus on the compact binary problem. The following items are certainly relevant to any evolution problem in computational physics: initial data, evolution, analysis, and numerics. We must specify the initial conditions, integrate the equations of motion to obtain the evolved data, and perform an analysis of the evolved data to extract physical information. The numerical treatment of each item may require the implementation of specific numerical techniques.

### Evolution

*Formulation: *Choose one of many inequivalent formulations (i.e., choose variables and rewrite the Einstein equations to obtain a well-posed initial value problem). Choose the order of time and space derivatives; make structural choices about the gauge and the constraints.

Reformulating the Einstein equations as a well-posed initial value problem has been the subject of much research (*10*, (*11*). To give an example, the result of the generalized harmonic gauge (GHG) formulation (*12*) can be cast in a standard first-order PDE form as(4)Here, the state vector *u*^{μ} collects all 10 components *g _{ab}*, the 40 first derivatives ∂

*, and a few additional fields depending on the formulation. In addition, there may be variables for the matter fields. Equation 4 for GHG is strongly and even symmetric hyperbolic (*

_{c}g_{ab}*10*, (

*11*). To give an example, the result of the generalized harmonic gauge (GHG) formulation (

*12*) and is well suited for numerical implementation. Another standard way to proceed is the classic ADM formulation that makes the geometry of the 3+1 decomposition in time and space more explicit. Basic variables are the 3-metric

*g*and the extrinsic curvature

_{ij}*K*, which is essentially the first time derivative of the metric (

_{ij}*8*). The ADM equations are only weakly hyperbolic and are not suitable for numerics. However, closely related systems, the so-called BSSN and Z4c formulations (

*10*,

*11*), are strongly hyperbolic. Most current simulations in numerical relativity rely on either the GHG or BSSN/Z4c families of formulations.

*Constraint propagation:* Maintain the constraints during evolution. Perform free evolutions and monitor the convergence of the constraints; use, for example, constraint damping to maintain the constraints explicitly.

Analytically, the constraints propagate; that is, if they are satisfied initially, they remain satisfied during a well-posed evolution. Numerically, even small rounding errors can trigger divergence from the constraint-satisfying solution, which can lead to a catastrophic failure of the simulation. How the constraints are controlled is a distinguishing feature of each formulation. A key ingredient in stable binary black hole evolutions (*13*) is the constraint-damping scheme (*14*). The Z4c formulation improves BSSN in the way the Hamiltonian constraint is treated, which leads to improved conservation of mass for neutron star simulations (*15*). Apart from instabilities, constraint violations in 3+1 relativity signify a problem with four-dimensional covariance. The 3+1 decomposition breaks covariance of the full theory by choosing a foliation, but the constraints ensure that four-dimensional covariance is maintained.

*Gauge:* Choose a coordinate condition—for example, in terms of lapse and shift or in terms of gauge source functions. Construct coordinates that avoid physical and coordinate singularities and are suitable for the black hole problem.

The main point about the gauge choice is that not only do we have the freedom to choose coordinates, but it is necessary to choose nontrivial coordinates. For example, even for the simplest black hole spacetimes, a foliation can fail by running into the physical singularity, and the hypersurface (or slice) may become badly distorted by slice stretching when points start falling into the black hole. The topic of how to dynamically construct good coordinates that lead to stable evolutions, cover spacetime with a regular foliation, avoid coordinate singularities, and avoid physical singularities inside black holes has become its own area of interest. In that context, the 3+1 decomposition is about “spacetime engineering” because we not only evolve the metric variables, but also build up the spacetime slice by slice in coordinates that are constructed dynamically during the evolution. The GHG formulation relies on the harmonic gauge to obtain hyperbolicity (*12*). For BSSN, the moving puncture gauge is essential to obtain long-term black hole evolutions, preventing slice stretching (*16*) and allowing the black hole punctures to move freely (*17*, *18*).

*Boundary conditions:* Specify outer boundary conditions appropriate for outgoing waves and asymptotic flatness. Specify inner boundaries for black holes; choose between black hole excision and black hole punctures. Handle coordinate patch boundaries.

Although approximate boundary conditions are possible, for a clean treatment, strong or symmetric hyperbolicity is required for well-posedness (*11*). We can then specify boundary conditions in terms of the ingoing and outgoing characteristic fields. The outer boundary conditions in numerical relativity tend to be substantially more complicated than the Einstein equations themselves, because outgoing-wave boundaries are typically constructed by taking additional derivatives of the right sides of the equations (*10*). The Einstein equations share with other nonlinear wave equations the feature that there is backscattering by waves off themselves (and, for binary systems, also due to the gradient in the gravitational well). This is a fundamental problem for boundaries at finite radius, because in principle we must account for all future backscattering from outside the domain. Consistent boundaries at finite radius have only been addressed quite recently, considering the long history of the Einstein equations (*19*, *20*).

### Initial data

*Formulation:* Rewrite the constraints as elliptic equations, identifying suitable free and dependent variables.

To give an indication of what the formulation of the constraints entails (*8*), consider a conformal rescaling of the metric, , which conveniently transforms the Hamiltonian constraint into a scalar elliptic equation for ψ. We have the freedom to specify a conformal metric, , that is not physical because it does not solve the constraints, but by solving the elliptic equation for ψ we find a physical solution *g _{ab}* that solves the constraints. The conformal transverse-traceless decomposition (

*8*) is widely used for the full set of constraints; for neutron star initial data in particular, the conformal thin-sandwich construction (

*21*,

*22*) is used, where typically an additional elliptic equation is added to initialize the gauge condition.

*Physical content: *Solve the constraints for data that contain multiple black holes or neutron stars with arbitrary mass, spin, and momentum.

Because the constraints are nonlinear, we cannot simply “add up” the metric tensors of, for example, two Schwarzschild black holes to obtain binary data, although that can be a useful approximate initial guess. As a result, some aspects of the initial data construction are indirect. For example, we can start with two single black hole solutions for particular masses, which will be combined to form a binary. But solving the constraints for the binary leads to a change in the individual masses of the black holes because of the conformal rescaling. In some cases we have to perform evolutions to determine whether the initial data were constructed appropriately for a particular dynamical situation.

There is a growing variety of initial data constructions for binaries that correspond to the variety of physical configurations. For black holes, there are excision-type data, where the interior of black holes is removed (*23*, *24*). Alternatively, black hole puncture data handle the black hole interior with a coordinate singularity at a point (*25*), which sometimes is called automatic excision. The thin-sandwich formulation is well suited for quasi-equilibrium data of black holes and/or neutron stars, which, for example, can approximate the state of a binary system during a quasi-circular inspiral (*26*). Only quite recently have methods been developed for neutron stars that generalize the quasi-equilibrium, quasi-circular construction to eccentric orbits (*27*) and to neutron stars with spin (*28*) (Fig. 2). Solving the constraints for electromagnetic field configurations is another recent topic of investigation (*29*).

### Analysis

*Black holes and neutron stars:* Determine all physical parameters during evolution. Find horizons of black holes. Analyze the rich phenomenology of neutron star mergers with the remnant, torus, jets, and ejecta. Connect to multi-messenger astronomy.

In any binary simulation, a wide range of detailed information is of interest, especially when matter is involved. The “relativity” in general relativity means, however, that many quantities have no direct physical meaning. In general, any tensor component (such as *g _{tt}* or

*g*) is not meaningful by itself; we have to construct proper gauge-invariant quantities. For example, mass and spin must be carefully defined because their local meaning at a point is problematic. For black holes, special methods are required to find the event horizon, which is a global concept in spacetime and therefore expensive to compute. See Fig. 3 for examples. Instead, black hole excision relies on the apparent horizon [e.g., (

_{xy}*13*)].

*Gravitational waves: *Compute gravitational wave emission; control numerical and systematic errors. Produce gravitational wave templates in a form that is ready to use for gravitational wave detectors. Treat both waveform prediction and waveform analysis.

Gravitational waves are propagating variations in the metric tensor, and the challenge is to separate the physical waves from various coordinate effects. In the weak-field limit, we can define gravitational waves as small perturbations around a background metric, and a first-order gauge-invariant formalism can be used to eliminate leading-order gauge effects (*30*). Such methods are applicable because we assume that the detectors are located far from the source where an asymptotically flat background is available. In simulations, the numerical grids often include extra patches for the far zone [e.g., (*9*)], possibly at lower resolution (see below).

A major effort in numerical relativity is directed toward obtaining accurate waveforms with controlled error bars for long time intervals. For the signal-to-noise ratio of current observations, a sufficiently accurate waveform model may begin with a post-Newtonian approximation (assuming nonrelativistic speeds) for the initial inspiral, matched to 10 to 20 orbits up to and including the merger from numerical simulations of the full Einstein equations. Initially, the goal was to filter the signal out of the noise by matching against theoretical waveforms. However, as the quality of the signals is improving, the main goal of gravitational wave astronomy is to estimate unknown source parameters. For example, we need detailed waveform models to distinguish black hole mergers from neutron star mergers, determine masses and spins, etc. The first detection of gravitational waves by Advanced LIGO (*2*) was accompanied by a theory paper describing how the properties of GW150914 were deduced from the observational data (*4*). Only by combining data with theory was it possible to arrive at the interpretation of GW150914 as the signature of a binary black hole merger, with specific parameters and credibility intervals. Two families of models were used, the EOBNR and Phenom families of waveforms (*2*) (Fig. 4). To analyze the data stream from the detectors, various parametrized waveform models are being developed for high-speed template matching (e.g., reduced-order surrogate models) (*31*).

### Numerics

*Discretization:* Choose a discretization in space and time. Introduce adaptive mesh refinement (AMR) in space and time to efficiently represent different physical length scales. Choose coordinate patches and transformations to adapt coordinates to the underlying physics.

Once a suitably hyperbolic form of the PDEs of general relativity has been derived, we have access to several standard discretizations from applied mathematics. The recent trend has been toward high-order discretizations, with different choices for the geometry and the matter fields. In vacuum or where the matter is smooth, the geometry is smooth as well. For smooth metrics, fourth- to eighth-order finite differencing in space is applied routinely, as well as pseudospectral methods for exponential convergence. Neutron star matter is represented by general relativistic fluids, and handling relativistic shocks becomes important. Several high-resolution shock-capturing (HRSC) fifth-order methods are available (*6*), as is work on smoothed particle hydrodynamics (*32*, *33*).

The physics of a binary involves several physical scales. The wavelength of gravitational waves near merger is about 100 times the size of the black holes, and the simulation domain is typically chosen to be at least 1000 times the size of the black holes to accommodate several wave cycles. Simulations in three spatial dimensions therefore become several orders of magnitude more efficient with AMR, often of the Berger-Oliger type with refinements not just in space, but also in time. Many codes use several coordinate patches to transition from two (or more) central objects to spherical shells near the outer boundary.

*Scientific computing:* Implement parallel algorithms for high-performance computing. Invest in professional software engineering for a collaborative computational infrastructure.

Numerical relativity has been very successful with the hybrid MPI (message passing interface) plus OpenMP (open multiprocessing) or a similar parallelization strategy. Still, a typical numerical relativity simulation for a binary coalescence, representing just a single data point in a template catalog, may take roughly 1 month on 1000 to 10,000 cores of a supercomputer. The numerical relativity community is working on improving the efficiency of these methods, including spectral methods and improved AMR schemes, which tend to be a bottleneck for massive parallelism. Most efforts in numerical relativity are group efforts with a long-term investment in an evolving code base. These efforts include SpEC (*34*), SACRA (*35*), Whisky/THC (*36*), Pretorius (*37*), HAD (*38*), BAM (*39*), and the community code Einstein Toolkit (*40*). Some codes approximate general relativity but provide more advanced neutron star physics (*32*, *33*). Although similar in some regards—after all, the same or similar physics is studied—the different projects vary greatly in the range and the specifics of the physics modules, the flexibility and extensibility of the codes, the level of software optimization, and the collaboration and code-sharing models.

The main challenge common to all these projects is that they are implementing a “moving target,” as formulations and basic equations are still changing and more physics is added to the simulations. Simultaneously, they must handle the trend in technology toward massively parallel computers and heterogeneous hardware, which is challenging given the complex algorithms required for numerical relativity.

## Short history of binary simulations

The first simulations of black holes in vacuum were attempted in 1964 (*41*). By the 1970s, many concepts of the 3+1 ADM formulation had been brought into numerical relativity (*42*), which led to the seminal numerical work on head-on (axisymmetric, 2+1-dimensional) black hole collisions and gravitational waves (*43*, *44*). It took until the early 1990s (*45*, *46*) to revisit the head-on collision with improved numerics, which confirmed the early results on gravitational waves (*46*). Numerical relativity in 3+1 dimensions began in 1995 with the evolution of a Schwarzschild black hole on a Cartesian grid (*47*) and the evolution of gravitational waves (*48*), followed by the first fully 3+1-dimensional simulation of a black hole binary (*49*, *50*). All the early black hole simulations mentioned so far were numerically unstable, with barely enough evolution time to start with two separate black holes that promptly merged. The first full orbit was achieved in 2004 (*51*). In 2005–2006, the last missing ingredients for long-term stable black hole evolutions were found in two different approaches, one based on a harmonic gauge formulation and excision (*13*) and the other based on the BSSN formulation and black hole punctures (*17*, *18*, *51*). By 2010, the robustness and flexibility of these methods had been established. Improvements in the formulations, the boundary conditions, etc., are still ongoing today (*11*, *12*, *51*).

Neutron star simulations were pursued in parallel with the black hole simulations. The Valencia formalism of general relativistic hydrodynamics (GRHD), now the primary approach, was developed in the 1990s (*52*). The first fully general relativistic binary neutron star simulations were published in 2000 (*53*), with enormous progress in many groups since then. As far as the geometry of general relativity matters in these simulations, it turns out that the methods established for stable black hole simulations carry over to neutron star simulations (gauge, boundaries, initial data formulation, etc.). However, GRHD introduces its own challenge of relativistic shocks, and the range of different physics phenomena makes this a much more complex problem than black holes in vacuum.

## Outlook

Numerical relativity is developing rapidly in several directions, and we highlight a few representative examples.

### High-order methods

High-order methods to address the ever-increasing demand for even more accurate and detailed simulations are a major topic of current research. Among the different high-order methods to solve partial differential equations, the discontinuous Galerkin (DG) method has emerged in recent years as a particularly successful general-purpose paradigm (*54*). It can be argued that the DG spectral-element method subsumes several of the key advantages of traditional finite-element and finite-volume methods. In particular, the DG method works with element-local stencils, which is a great advantage for parallelization and the construction of complicated grids. Furthermore, DG methods offer easy access to *hp*-adaptivity, where both the size of the computational elements (or cells) and the order of the polynomial approximation within each element can be adapted to the problem.

There are three major efforts to use DG methods for general relativity and/or GRHD (*55*–*57*). The first simulations of a single neutron star were achieved recently (*55*, *58*), and simple binaries are a work in progress. With regard to high-order approximations, there is no doubt that if exponentially convergent spectral methods such as DG (or pseudospectral methods) are applicable, they will constitute a big improvement over finite-difference approximations, which give only polynomial convergence. High-order methods can provide breakthroughs by reaching accuracies that make new physics possible (e.g., for magnetic field amplification due to small-scale turbulence) or by reducing numerical errors to make gravitational wave analysis possible. Viewed differently, we can reach a given error criterion with much lower computational resources, making simulations feasible that are otherwise too computationally expensive.

### Multi-physics

The spectacular first observation of both gravitational waves (*3*) and electromagnetic radiation (*59*, *60*) from a neutron star merger represents the beginning of multi-messenger astronomy including gravitational waves. To model such systems, we need to perform “multi-physics” simulations.

Modeling electromagnetic fields in GRHD can be accomplished by coupling the Maxwell equations to the GRHD equations, for which the prevalent approach has been ideal magnetohydrodynamics (IMHD). The assumption of IMHD is that the fluid has zero resistivity, but for the merger—and in particular for the fields surrounding the remnant with torus and ejecta—the quality of that approximation is unclear. Resistive magnetohydrodynamics (RMHD) is expected to be important for realistic models of plasma instabilities and magnetic reconnection. Apart from unknown physics, the mathematical character of the RMHD equations may be problematic (*61*, *62*). There are only a few general relativistic simulations with RMHD [e.g., (*61*, *63*, *64*)]. Developing a proper theory of resistive relativistic plasmas is a large project in itself (*65*).

The microphysical equation of state of neutron stars remains unknown and is also a target for numerical models and for observations. Investigations may include 20 or more different equations of state in an attempt to cover all sensible proposals. Even determining just one parameter—the existence of neutron stars with 2.0 solar masses (*66*, *67*)—provided a strong constraint. In principle, gravitational wave observations can do much better, gleaning information from the inspiral and the merger. Although inspiral signals will show rather systematic long-time effects (*68*–*70*), one of the grand challenges will be to disentangle the much more messy merger signal (*71*).

Standard merger models predict strong heating of the neutron star matter, which is expected to lead to an enormous amount of neutrino emission with luminosity on the order of 10^{54} erg s^{–1}. This burst of energy plays a role in models of short gamma ray bursts (*72*) and also for the ejecta, which in turn affects heavy-element production and macro- or kilonovae (*73*). However, currently the high dimensionality of such radiative transport problems (3+1 spacetime plus 3 for the radiative transport) is prohibitive, leading to a wide array of approximations with variable applicability (*74*, *75*). A coherent picture for neutrino physics in binary mergers is still lacking but should be a part of multi-messenger astrophysics.

### Beyond current astrophysics

Numerical relativity has a large number of applications outside the area of compact binaries and gravitational waves (*76*, *77*). Topics include gravitational collapse with surprising critical phenomena, boson stars and other exotic matter, and cosmological simulations. Going beyond classical general relativity, the field of numerical relativity for alternative gravity theories and gravity in higher dimensions is wide open.

## Conclusion

The next decade is sure to see numerical relativity grow in terms of computational power and applicability to different physical scenarios. The detailed theoretical models for black hole and neutron star binaries that are the target of research in numerical relativity are closely linked to the observation of gravitational waves. Numerical relativity, in combination with the highly anticipated future observations of gravitational waves, is expected to provide entirely new insights into extreme gravity and extreme matter.

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**I gratefully acknowledge the joint work evident from the list of references. Without my collaborators, this review would not have been possible.

**Funding:**Supported in part by DFG/NSF grant BR 2176/5-1.

**Author contributions:**B.B. is responsible for the entire manuscript.

**Competing interests:**None.

**Data and materials availability:**There are no new data in this review.