List of Figures

	Figure 1: Range of various approximation tools (“UR” stands for ultra-relativistic). NR is mostly limited by resolution issues and therefore by possible different scales in the problem.
	Figure 2: Illustration of two hypersurfaces of a foliation $Σt$ . Lapse $α$ and shift $βμ$ are defined by the relation of the timelike unit normal field $μ n$ and the basis vector $∂t$ associated with the coordinate $t$ . Note that $⟨dt,αn ⟩ = 1$ and, hence, the shift vector $β$ is tangent to $Σt$ .
	Figure 3: $D$ -dimensional representation of head-on collisons for spinless BHs, with isometry group $SO (D − 2)$ (left), and non-head-on collisons for BHs spinning in the orbital plane, with isometry group $SO (D − 3 )$ (right). Image reproduced with permission from [841], copyright by APS.
	Figure 4: Illustration of mesh refinement for a BH binary with one spatial dimension suppressed. Around each BH (marked by the spherical AH), two nested boxes are visible. These are immersed within one large, common grid or refinement level.
	Figure 5: Illustration of singularity excision. The small circles represent vertices of a numerical grid on a two-dimensional cross section of the computational domain containing the spacetime singularity, in this case at the origin. A finite region around the singularity but within the event horizon (large circle) is excluded from the numerical evolution (white circles). Gray circles represent the excision boundary where function values can be obtained through regular evolution in time using sideways derivative operators as appropriate (e.g., [630]) or regular update with spectral methods (e.g., [677, 678]), or through extrapolation (e.g., [703, 723]). The regular evolution of exterior grid points (black circles) is performed with standard techniques using information also on the excision boundary.
	Figure 6: Illustration of the conjectured mass-scaling relation (172 ). The data refer to three separate one-parameter variations of the pulse shape (171 ). The constants $αi$ and $βi$ are chosen to normalize the ranges of the abscissa and place the data point corresponding to the smallest BH in each family at the origin. Image reproduced with permission from [212], copyright by APS.
	Figure 7: Left panel: Embedding diagram of the AH of the perturbed black string at different stages of the evolution. The light (dark) lines denote the first (last) time from the evolution segment shown in the corresponding panel. Right panel: Dimensionless Kretschmann scalar $𝒦2$ at the centre of mass of a binary BH system as a function of the (areal) coordinate separation between the two BHs in a $D = 5$ scattering, in units of $Rg = RS$ . Images reproduced with permission from (left) [511] and from (right) [587], copyright by APS.
	Figure 8: Snapshots of the rest-mass density in the collision of fluid balls with boost factor $γ = 8$ (upper panels) and $γ = 10$ (lower panels) at the initial time, shortly after collision, at the time corresponding to the formation of separate horizons in the $γ = 10$ case, and formation of a common horizon (for $γ = 10$ ) and at late time in the dispersion ( $γ = 8$ ) or ringdown ( $γ = 10$ ) phase. Image reproduced with permission from [288], copyright by APS.
	Figure 9: Instability against BH formation in AdS (left panel) and Minkowski enclosed in a cavity (right panel). In both panels, the horizontal axis represents the amplitude of the initial (spherically symmetric) scalar field perturbation. The vertical axis represents the size of the BH formed. Perturbations with the largest plotted amplitude collapse to form a BH. As the amplitude of the perturbation is decreased so does the size of the BH, which tends to zero at a first threshold amplitude. Below this energy, no BH is formed in the first generation collapse and the scalar perturbation scatters towards the boundary. But since the spacetime behaves like a cavity, the scalar perturbation is reflected off the boundary and re-collapses, forming now a BH during the second generation collapse. At smaller amplitudes a second, third, etc, threshold amplitudes are found. The left (right) panel shows ten (five) generations of collapse. Near the threshold amplitudes, critical behavior is observed. Images reproduced with permission from (left) [108] and from (right) [537], copyright by APS.
	Figure 10: (a) and (b): $+$ and $×$ modes of gravitational waveform (solid curve) from an unstable six-dimensional BH with $q = 0.801$ as a function of a retarded time defined by $t − r$ , where $r$ is the coordinate distance from the center. Image reproduced with permission from [700], copyright by APS.
	Figure 11: Evolution of a highly spinning BH ( $a∕M = 0.99$ ) during interaction with different frequency GW packets, each with initial mass $≈ 0.1M$ . Shown (in units where $M = 1$ ) are the mass, irreducible mass, and angular momentum of the BH as inferred from AH properties. Image reproduced with permission from [289], copyright by APS.
	Figure 12: Massive scalar field (nonlinear) evolution of the spacetime of an initially non-rotating BH, with $M μ = 0.29$ . Left panel: Evolution of a spherically symmetric $l = m = 0$ scalar waveform, measured at $rex = 40M$ , with $M$ the initial BH mass. In addition to the numerical data (black solid curve) we show a fit to the late-time tail (red dashed curve) with $t−0.83$ , in excellent agreement with linearized analysis. Right panel: The dipole signal resulting from the evolution of an $l = m = 1$ massive scalar field around a non-rotating BH. The waveforms, extracted at different radii $rex$ exhibit pronounced beating patterns caused by interference of different overtones. The critical feature is however, that these are extremely long-lived configurations. Image reproduced with permission from [588], copyright by APS.
	Figure 13: BH trajectories in grazing collisions for $γ = 1.520$ and three values of the impact parameter corresponding to the regime of prompt merger (solid, black curve), of delayed merger (dashed, red curve), and scattering (dotted, blue curve). Note that for each case, the trajectory of one BH is shown only; the other BH’s location is given by symmetry across the origin.
	Figure 14: Total energy radiated in GWs (left panel) and final dimensionless spin of the merged BH (right panel) as a function of impact parameter $b$ for the same grazing collisions with $γ = 1.520$ . The vertical dashed (green) and dash-dotted (red) lines mark $b∗$ and $bscat$ , respectively. Image reproduced with permission from [720], copyright by APS.
	Figure 15: Left panels: Scattering threshold (upper panel) and maximum radiated energy (lower panel) as a function of $v$ . Colored “triangle” symbols pointing up and down refer to the aligned and antialigned cases, respectively. Black “circle” symbols represent the thresholds for the nonspinning configurations. Right panel: Trajectory of one BH for a delayed merger configuration with anti-aligned spins $j = 0.65$ . The circles represent the BH location at equidistant intervals $Δt = 10M$ corresponding to the vertical lines in the inset that shows the equatorial circumference of the BH’s AH as a function of time.
	Figure 16: The (red) plus and (blue) circle symbols mark scattering and merging BH configurations, respectively, in the $b − v$ plane of impact parameter and collision speed, for $D = 5$ spacetime dimensions.
	Figure 17: Energy fluxes for head-on collisions of two BHs in $D = 5$ spacetime dimensions, obtained with two different codes, HD-Lean [841, 797] (solid black line) and SacraND [820, 587] (red dashed line). The BHs start off at an initial coordinate separation $d∕RS = 6.47$ . Image adapted from [796].
	Figure 18: Trajectories of BHs immersed in a scalar field bubble of different amplitudes. The BH binary consists of initially non-spinning, equal-mass BHs in quasi-circular orbit, initially separated by $11M$ , where $M$ is the mass of the binary system. The scalar field bubble surrounding the binary has a radius $r0 = 120M$ and thickness $σ = 8M$ . Panels $A, B, C$ correspond to $φ0 = 0(GR ),1∕80,1 ∕40$ and a zero potential amplitude $λ$ . Panel $D$ corresponds to $3 2 φ0 = 1∕80, 4πλ = 10 M$ . Image reproduced with permission from [410], copyright by IOP. All rights reserved.
	Figure 19: Numerical results for a BH binary inspiralling in a scalar-field gradient $−7 M σ = 10$ . Left panel: dependence of the various components of the scalar radiation $Re(φ11)∕(M σ)$ on the extraction radius (top to bottom: $112M$ to $56M$ in equidistant steps). The dashed line corresponds instead to $10− 3Im (φ )∕(M σ) 11$ at the largest extraction radius. This is the dominant mode and corresponds to the fixed-gradient boundary condition, along the $z$ -direction, at large distances. Right panel: time-derivative of the scalar field at the largest and smallest extraction radii, rescaled by radius and shifted in time. Notice how the waveforms show a clean and typical merger pattern, and that they overlap showing that the field scales to good approximation as $1∕&tidle;r$ . Image reproduced with permission from [92], copyright by APS.
	Figure 20: The dominant quadrupolar component of the gravitational $ψ4$ scalar for an equal-mass, non-spinning NS binary with individual baryon masses of $1.625M ⊙$ . The solid (black) curve refers to GR, and the dashed (red) curve to a scalar-tensor theory with $− 5 β∕(4π) = − 4.5, φ0 = 10$ . Image reproduced with permission from [73], copyright by APS.
	Figure 21: Left panel: Collision of two shock waves in AdS₅. The energy density $4 ℰ ∕μ$ is represented as a function of an (advanced) time coordinate $v$ and a longitudinal coordinate $z$ . $μ$ defines the amplitude of the waves. Right panel: Evolution of the scalar field in an unstable RN-AdS BH. $z$ is a radial coordinate and the AdS boundary is at $z = 0$ . Due to the instability of the BH, the scalar density grows exponentially for $0 < tTc ≲ 6$ . Then, the scalar density approaches some static function. Images reproduced with permission from (left) [207], copyright by APS and (right) [562], copyright by SISSA.
	Figure 22: Left: Elementary cells for the 8-BH configuration, projected to $3 R$ . The marginal surface corresponding to the BH at infinity encompasses the whole configuration. Note that the 8 cubical lattice cells are isometric after the conformal rescaling. Right: Several measures of scaling in the eight-BH universe, as functions of proper time $τ$ , plotted against a possible identification of the corresponding FLRW model (see Ref. [86] for details). All the quantities have been renormalized to their respective values at $τ = 0$ . Images reproduced with permission from [86], copyright by IOP. All rights reserved.

	Abstract
	Acronyms
	Notation and conventions
1	Prologue
2	Milestones
3	Strong Need for Strong Gravity
	3.1	Astrophysics
	3.2	Fundamental and mathematical issues
	3.3	High-energy physics
4	Exact Analytic and Numerical Stationary Solutions
	4.1	Exact solutions
	4.2	Numerical stationary solutions
5	Approximation Schemes
	5.1	Post-Newtonian schemes
	5.2	Spacetime perturbation approach
	5.3	The zero-frequency limit
	5.4	Shock wave collisions
6	Numerical Relativity
	6.1	Formulations of the Einstein equations
	6.2	Higher-dimensional NR in effective “3 + 1” form
	6.3	Initial data
	6.4	Gauge conditions
	6.5	Discretization of the equations
	6.6	Boundary conditions
	6.7	Diagnostics
7	Applications of Numerical Relativity
	7.1	Critical collapse
	7.2	Cosmic censorship
	7.3	Hoop conjecture
	7.4	Spacetime stability
	7.5	Superradiance and fundamental massive fields
	7.6	High-energy collisions
	7.7	Alternative theories
	7.8	Holography
	7.9	Applications in cosmological settings
8	Conclusions
	Acknowledgements
	References
	Footnotes
	Figures