Numerical relativity begin in the mid-1960s and had a major breakthrough in 2005. The LIGO gravitational wave observatory had started collecting data in 2002, so there was a strong impetus to be able to match theoretical simulations of merging black holes to observations. This paid off in 2016 when LIGO made its first detection.
Einstein's equations are ten coupled, non-linear, second-order partial differential equations in four dimensions... a formidable computational challenge!
The first problem was hardware limitations. As late as 1995, physicists could not even numerically solve the equations for the simple, analytically-known, spherically symmetric Schwarzschild metric, due to the complications of dealing with the singularity. The supercomputers of that time did not have sufficient memory and computational power to perform accurate calculations of 3D spacetimes.
In a situation without any spatial symmetry, the number of 3D grid points in the discretized equations are enormous if you want to have decent resolution. But within a few years, progress was made on head-on collisions of black holes, exploiting the cylindrical symmetry. Eventually hardware reached the point where it was no longer the bottleneck, even in situations where symmetry could not be exploited. But a long series of other computational challenges had to be overcome.
The first was formulating the equations in a way that made them a well-posed initial boundary-value problem with satisfactory numerical stability. The Arnowitt-Deser-Misner (ADM) "3+1" formalism had been around since 1959. This is a Hamiltonian approach in which spacetime is foliated into a spacelike 3D slices, each with its own internal 3D metric and extrinsic curvature, which evolve in time. It reduces the Einstein equations to twelve coupled first-order-in-time evolution equations (six for the 3-metric, six for the extrinsic curvature), plus four constraint equations. This formalism was suitable for numerically evolving an initial spacetime slice forward in time, but keeping numerical errors from building up was problematic because the equations were only "weakly hyperbolic".
The Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formalism, developed between 1987 and 1999 overcame this problem by reformulating the ADM equations to make them "strongly hyperbolic", a condition which makes for much better numerical stability.
Next, since general relativity is a gauge theory, another challenge was the question of which of various possible gauges was best for doing calculations. It turned out to be highly non-trivial to find gauge conditions that ensure numerically stable evolutions, but eventually a family of gauges called generalized harmonic gauge (GHG) proved suitable.
The question of formulating appropriate data for initial conditions was a difficult one. Not only did the initial data need to be physically correct -- for example, to describe two orbiting black holes, each with spin -- but it also needed to satisfy the four constraint equations.
Dealing with spatial boundary conditions at infinity were another roadblock. Far from the two black holes, spacetime must take the form of outgoing gravitational radiation. The numerical solution must ensure that there is no gravitational radiation coming in from infinity.
Mesh refinement turned out to be necessary to handle the various distance scales that black holes have, all the way from their horizon to the wave zone. And this mesh refinement had to be implemented in a way that could be parallelized on multiple processors.
Extracting physical results, such as gravitational wave waveforms, in a gauge-invariant manner from the numerical simulation was non-trivial.
Dealing with the singularity of each hole, and of the merged hole, was a major problem. Two different techniques were developed. In the "excision" technique, proposed in the late 1990's, a region around the singularity, but inside the horizon, is simply not evolved, since nothing happening inside the hole can affect the outside.
The second technique, called the "puncture method", divided the solution into an analytical part containing the singularity, and a numerically-constructed part that was singularity-free. But at first, the puncture containing the singularity remained at fixed coordinates even as the holes moved, resulting in the coordinate system getting stretched and warped to the point that numerical instabilities arose.
The breakthrough of 2005 was allowing the punctures to move through the coordinate system to control the numerical instabilities. After that point, spacetimes for merging black holes could be accurately simulated.
Forty years of hard work had brought the field of numerical relativity to maturity!
This post was based on two sources: "Numerical relativity" (https://en.wikipedia.org/wiki/Numerical_relativity) and "The numerical relativity breakthrough for binary black holes" (https://arxiv.org/abs/1411.3997).
