of freedom (gravitational waves). An entire subfield of gravitational physics, *numerical relativity*, spent nearly five decades looking for ways to solve Einstein’s equations on computers (see [@baumgarte:10book] for a comprehensive introduction). In general relativity, changes in the gravitational field propagate at the speed of light. Hence, time evolution equations must be solved. This is done by splitting 4D spacetime into 3D spatial slices that are evolved in the time direction. In the simplest way of writing the equations (the so-called Arnowitt-Deser-Misner \[ADM\] formulation), they form a system of 12 partial differential evolution equations, 4 gauge variables that must be specified (and evolved in time or recalculated on each slice), and 4 elliptic *constraint* equations without time derivatives. The ADM formulation has poor numerical stability properties. These lead to violations of the constraint equations and numerical instabilities that make long-term evolution impossible.

    It took until the 2000s for numerical relativity to find formulations of Einstein’s equations and gauge choices that together lead to stable long-term evolutions. In some cases, well-posedness and strong or symmetric hyperbolicity can be proven. The equations are typically evolved time-explicitly with straightforward high-order (fourth and higher) finite difference schemes or with multi-domain pseudospectral methods.

    Since numerical relativity 