Applications of Numerical Relativity

7 Applications of Numerical Relativity

Numerical relativity was born out of efforts to solve the two-body problem in GR, and aimed mainly at understanding stellar collapse and GW emission from BH and NS binaries. There is therefore a vast amount of important results and literature on NR in astrophysical contexts. Because these results fall outside the scope of this review, we refer the interested reader to Refs. [631 , 592 , 191 , 715 , 18 , 617 , 429 ] and to the relevant sections of Living Reviews¹⁸ for (much) more on this subject. Instead, we now focus on applications of NR outside its traditional realm, most of which are relatively recent new directions in the field.

7.1 Critical collapse

The nonlinear stability of Minkowski spacetime was established by Christodoulou and Klainerman, who showed that arbitrarily “small” initial fluctations eventually disperse to infinity [219 ]. On the other hand, large enough concentrations of matter are expected to collapse to BHs, therefore raising the question of how the threshold for BH formation is approached.

Choptuik performed a thorough investigation of this issue, by evolving initial data for a minimally coupled massless scalar field [212 *]. Let the initial data be described by a parameter $p$ which characterizes the initial scalar field wavepacket. For example, in Choptuik’s analysis, the following family of initial data for the scalar field $Φ$ was considered,

where $′ Φ = ϕ$ ; therefore any of the quantities $ϕ0,r0,δ,q$ is a suitable parameter $p$ .

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.

The evolution of such initial data close to the threshold of BH formation is summarized in Figure 6 *. Fix all but one parameter, say the scalar field amplitude $p = ϕ0$ . For large amplitudes $ϕ0$ , a large BH is formed. As the amplitude of the initial data decreases, the mass of the formed BH decreases, until a critical threshold amplitude $ϕ0∗$ is reached below which no BH forms and the initial data disperses away (consistently with the nonlinear stability of Minkowski). Near the threshold, BHs with arbitrarily small masses can be created, and the BH mass scales as

It was found that $γ ≈ 0.37$ is a universal (critical) exponent which does not depend on the initial data, or in other words it does not depend on which of the parameters $ϕ0,r0,δ,q$ is varied (but it may depend on the type of collapsing material).

The BH threshold in the space of initial data for GR shows both surprising structure and surprising simplicity. In particular, critical behavior was found at the threshold of BH formation associated with universality, power-law scaling of the BH mass, and discrete self-similarity, which bear resemblance to more familiar statistical physics systems. Critical phenomena also provide a route to develop arbitrarily large curvatures visible from infinity (starting from smooth initial data) and are therefore likely to be relevant for cosmic censorship (see Section 7.2), quantum gravity, astrophysics, and our general understanding of the dynamics of GR.

Choptuik’s original result was extended in many different directions, to encompass massive scalar fields [125 , 586 *], collapse in higher dimensions [344 ] or different gravitational theories [265 ]. Given the difficulty of the problem, most of these studies have focused on 1 + 1 simulations; the first non spherically (but axially) symmetric simulations were performed in Ref. [9 ], whereas recently the first 3 + 1 simulations of the collapse of minimally coupled scalar fields were reported [411 ]. The attempt to extend these results to asymptotically AdS spacetimes would uncover a new surprising result, which we discuss below in Section 7.4. A full account of critical collapse along with the relevant references can be found in a Living Reviews article on the subject [381 ].

7.2 Cosmic censorship

As discussed in Section 3.2.1, an idea behind cosmic censorship is that classical GR is self-consistent for physical processes. That is, despite the fact that GR predicts the formation of singularities, at which geodesic incompleteness occurs and therefore failure of predictablity, such singularities should be – for physical processes¹⁹ – causally disconnected from distant observers by virtue of horizon cloaking. In a nutshell: a GR evolution does not lead, generically, to a system GR cannot tackle. To test this idea, one must analyze strong gravity dynamics, which has been done both using numerical evolutions and analytical arguments. Here we shall focus on recent results based on NR methods. The interested reader is referred to some historically relevant numerical [692 , 361 ] and analytical [218 , 653 ] results, as well as to reviews on the subject [768 , 88 *, 649 , 461 *] for further information.

The simplest (and most physically viable) way to violate cosmic censorship would be through the gravitational collapse of very rapidly rotating matter, possibly leading to a Kerr naked singularity with $a > M$ . However, NR simulations of the collapse of a rotating NS to a BH [568 , 63 , 348 ] have shown that when the angular momentum of the collapsing matter is too large, part of the matter bounces back, forming an unstable disk that dissipates the excess angular momentum, and eventually collapses to a Kerr BH. Simulations of the coalescence of rapidly rotating BHs [769 , 412 ] and NSs [465 ] have shown that the $a > M$ bound is preserved by these processes as well. These simulations provide strong evidence supporting the cosmic censorship conjecture. Let us remark that analytical computations and NR simulations show that naked singularities can arise in the collapse of ideal fluids [461 ] but these processes seem to require fine-tuned initial conditions, such as in spherically symmetric collapse or in the critical collapse [375 ] discussed in Section 7.1.

A claim of cosmic censorship violation in $D > 4$ spacetime dimensions was made in the context of the nonlinear evolution of the Gregory–Laflamme (see Section 3.2.4) instability for black strings. In Ref. [511 *] long-term numerical simulations were reported showing that the development of the instability leads to a cascade of ever smaller spherical BHs connected by ever thinner black string segments – see Figure 7 *, left (top, middle and bottom) panel for a visualization of the (first, second and third) generations of spherical BHs and string segments. Observe, from the time scales presented in Figure 7 *, that as viewed by an asymptotic observer, each new generation develops more rapidly than the previous one. The simulations therefore suggest that arbitrarily thin strings, and thus arbitrarily large curvature at the horizon, will be reached in finite asymptotic time. If true, this system is an example where a classical GR evolution is driving the system to a configuration that GR cannot describe, a state of affairs that will presumably occur when Planck scale curvatures are attained at the horizon. The relevance of this example for cosmic censorship, may, however, be questioned, based on its higher dimensionality and the lack of asymptotic flatness: cosmic strings with horizons require the spacetime dimension to be greater or equal to five and the string is infinitely extended in one dimension. In addition, the simulations of [511 *] assume cylindrical symmetry, and cylindrically symmetric matter configurations are unstable [174 ]; therefore, fine-tuning of initial conditions may be required for the formation of a naked singularity.

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.

Another suggestion that Planckian scale curvature becomes visible in a classical evolution in $D = 5$ GR arises in the high-energy scattering of BHs. In Ref. [587 *], NR simulations of the scattering of two non-spinning boosted BHs with an impact parameter $b$ were reported. For sufficiently small initial velocities ( $v ≲ 0.6c$ ) it is possible to find the threshold impact parameter $bscat = bscat(v)$ such that the BHs merge into a (spinning) BH for $b < bscat$ or scatter off to infinity for $b > bscat$ . For high velocities, however, only a lower bound on the impact parameter for scattering $bC = bC(v)$ and an upper bound on the impact parameter for merger $bB = bB(v)$ could be found, since simulations with $bB < b < bC$ crashed before the final outcome could be determined (cf. Figure 16 * below). Moreover, an analysis of a scattering configuration with $v = 0.7$ and $b = bC$ , shows that very high curvature develops outside the individual BHs’ AHs, shortly after they have reached their minimum separation – see Figure 7 * (right panel). The timing for the creation of the high curvature region, i.e., that it occurs after the scattering, is in agreement with other simulations of high energy collisions. For instance, in Refs. [216 *, 288 *] BH formation is seen to occur in the wake of the collision of non-BH objects, which was interpreted as due to focusing effects [288 *]. In the case of Ref. [587 *], however, there seems to be no (additional) BH formation. Both the existence and significance of such a high curvature region, seemingly uncovered by any horizon, remains mysterious and deserves further investigation.

In contrast with the two higher-dimensional examples above, NR simulations that have tested the cosmic censorship conjecture in different $D = 4$ setups, found support for the conjecture. We have already mentioned simulations of the gravitational collapse of rotating matter, and of the coalescence of rotating BH and NS binaries. As we discuss below in Section 7.6, the high-energy head-on collisions of BHs [719 *], boson stars [216 *] or fluid particles [288 *, 647 *] in $D = 4$ result in BH formation but no naked singularities. A different check of the conjecture involves asymptotically dS spacetimes [837 *]. Here, the cosmological horizon imposes an upper limit on the size of BHs. Thus one may ask what is the outcome of the collision of two BHs with almost the maximum allowed size. In Ref. [837 *] the authors were able to perform the evolution of two BHs, initially at rest with the cosmological expansion. They observe that for all the (small) initial separations attempted, a cosmological AH, as viewed by an observer at the center of mass of the binary BH system, eventually forms in the evolution, and both BH AHs are outside the cosmological one. In other words, the observer in the center of mass loses causal contact with the two BHs which fly apart rather than merge. This suggests that the background cosmological acceleration dominates over the gravitational attraction between ‘large’ BHs. It would be interesting to check if a violation of the conjecture can be produced by introducing opposite charges to the BHs (to increase their mutual attraction) or give them mutually directed initial boosts.

7.3 Hoop conjecture

The hoop conjecture, first proposed by K. Thorne in 1972 [750 ], states that when the mass $M$ of a system (in $D = 4$ dimensions) gets compacted into a region whose circumference in every direction has radius $R ≲ Rs = 2M$ , a horizon – and thus a BH – forms (for a generalization in $D > 4$ , see [449 ]). This conjecture is important in many contexts. In high-energy particle collisions, it implies that a classical BH forms if the center-of-mass energy significantly exceeds the Planck energy. This is the key assumption behind the hypothesis – in the TeV gravity scenario – of BH production in particle accelerators (see Section 3.3.2). In the trans-Planckian regime, the particles can be treated as classical objects. If two such “classical” particles with equal rest mass $m0$ and radius (corresponding to the de Broglie wavelength of the process) $R$ collide with a boost parameter $γ$ , the mass-energy in the centre-of-mass frame is $M = 2γm 0$ . The threshold radius of Thorne’s hoop is then $Rs = 4γm0$ and the condition $R ≲ Rs = 2M = 4γm0$ translates into a bound on the boost factor, $γ ≳ γh ≡ R ∕(4m0 )$ .

Even though the hoop conjecture seems plausible, finding a rigorous proof is not an easy task. In the last decades, the conjecture has mainly been supported by studies of the collision of two infinitely boosted point particles [256 , 286 *, 817 , 818 ], but it is questionable that they give an accurate description of an actual particle collision (see, e.g., the discussion in Ref. [216 *]). In recent years, however, advances in NR have made it possible to model trans-Planckian collisions of massive bodies and provided more solid evidence in favor of the validity of the hoop conjecture.

The hoop conjecture has been first addressed in NR by Choptuik & Pretorius [216 *], who studied head-on collisions of boson stars in four dimensions (see Section 4.2). The simulations show that the threshold boost factor for BH formation is $∼ 1∕3γh$ (where the “hoop” critical boost factor $γh$ is defined above), well in agreement with the hoop conjecture. These results have been confirmed by NR simulations of fluid star collisions [288 *], showing that a BH forms when the boost factor is larger than $∼ 0.42γh$ . Here, the fluid balls are modeled as two superposed Tolman–Oppenheimer–Volkoff “stars” with a $Γ = 2$ polytropic equation of state.

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.

The simulations furthermore show that for boosts slightly above the threshold of BH formation, there exists a brief period where two individual AHs are present, possibly due to a strong focusing of the fluid elements of each individual star caused by the other’s gravitational field. These results are illustrated in Figure 8 * which displays snapshots of a collision for $γ = 8$ that does not result in horizon formation (upper panels) and one at $γ = 10$ that results in a BH (lower panels). We remark that the similarity of the behaviour of boson stars and fluid stars provides evidence supporting the “matter does not matter” hypothesis discussed below in Section 7.6. A similar study of colliding fluid balls [647 ] has shown similar results. Therein it has been found that a BH forms when the compactness of the star is $m ∕R ≳ 0.08γ− 1.13 0$ , i.e., for $γ1.13 ≳ R ∕(12m ) = 1∕3γ 0 h$ . Type-I critical behaviour has also been identified, with BH formation for initial masses $m0$ above a critical value scaling as $−1.0 ∼ γ$ .

7.4 Spacetime stability

Understanding the stability of stationary solutions to the Einstein field equations, or generalisations thereof, is central to gauge their physical relevance. If the corresponding spacetime configuration is to play a role in a given dynamical process, it should be stable or, at the very least, its instabilities should have longer time scales than those of that dynamical process. Following the evolution of unstable solutions, on the other hand, may unveil smoking guns for establishing their transient existence. NR provides a unique tool both for testing nonlinear stability and for following the nonlinear development of unstable solutions. We shall now review the latest developments in both these directions, but before doing so let us make a remark. At the linear level, typical studies of space-time stability are in fact studies of mode stability. A standard example is Whiting’s study of the mode stability for Kerr BHs [780 ]. For BH spacetimes, however, mode stability does not guarantee linear stability, cf. the discussion in [236 ]. We refer the reader to this reference for further information on methods to analyse linear stability.

Even if a spacetime does not exhibit unstable modes in a linear analysis it may be unstable when fully nonlinear dynamics are taken into account. A remarkable illustration of this possibility is the turbulent instability of the AdS spacetime reported in Ref. [108 *]. These authors consider Einstein gravity with a negative cosmological constant $Λ$ and minimally coupled to a massless real scalar field $ϕ$ in $D = 4$ spacetime dimensions. The AdS metric is obviously a solution of the system together with a constant scalar field. Linear scalar-field perturbations around this solution generate a spectrum of normal modes with real frequencies [150 ]: $ωN L = 2N + 3 + ℓ$ , where $N ∈ ℕ0$ and $L, ℓ$ are the AdS length scale and total angular momentum harmonic index, respectively. The existence of this discrete spectrum is quite intuitive from the global structure of AdS: a time-like conformal boundary implies that AdS behaves like a cavity. Moreover, the fact that the frequencies are real shows that the system is stable against scalar-field perturbations at linear level.

The latter conclusion dramatically changes when going beyond linear analysis. Setting up spherically symmetric Gaussian-type initial data of amplitude $𝜖$ , Bizoń and Rostworowski [108 *] made the following observations. For large $𝜖$ the wave packet collapses to form a BH, signalled by an AH at some radial coordinate. As $𝜖$ is made smaller, the AH radius also decreases, reaching zero size at some (first) threshold amplitude. This behaviour is completely analogous to that observed in asymptotically flat spacetime by Choptuik [212 ] and discussed in Section 7.1; in fact, the solutions obtained with this threshold amplitude asymptote – far from the AdS boundary – to the self-similar solution obtained in the $Λ = 0$ case. For amplitudes slightly below the first threshold value, the wave packet travels to the AdS boundary, where it is reflected, and collapses to form a BH upon the second approach to the centre. By further decreasing the amplitude, one finds a second threshold amplitude at which the size of the AH formed in this second generation interaction decreases to zero. This pattern seems to repeat itself indefinitely. In [108 *] ten generations of collapse were reported, as shown in Figure 9 * (left panel). These results were confirmed and extended in subsequent work [142 ]. If indeed the pattern described in the previous paragraph repeats itself indefinitely, a remarkable conclusion is that, no matter how small the initial amplitude is, a BH will form in AdS after a time scale $𝒪 (𝜖−2)$ . A corollary is then that linear analysis misses the essential physics of this problem; in other words, the evolution always drives the system away from the linear regime.

The central property of AdS to obtain this instability is its global structure, rather than its local geometry. This can be established by noting that a qualitatively similar behaviour is obtained by considering precisely the same dynamical system in Minkowski space enclosed in a cavity [537 *], see Figure 9 * (right panel). Moreover, the mechanism behind the instability seems to rely on nonlinear interactions of the field that tend to shift its energy to higher frequencies and hence smaller wavelenghts. This process stops in GR since the theory has a natural cutoff: BH formation.

It has since been pointed out that collapse to BHs may not be the generic outcome of evolutions in AdS [145 *, 66 , 274 , 538 , 539 ]. For example, in Ref. [145 *] “islands of stability” were discovered for which the initial data, chosen as a small perturbation of a boson star, remain in a nonlinearly stable configuration. In Ref. [586 *], the authors raised the possibility that some of the features of the AdS instability could also show up in asymptotically flat spacetimes in the presence of some confinement mechanism. They observed that the evolution of minimally coupled, massive scalar wavepackets in asymptotically flat spacetimes can also lead to collapse after a very large number of “bounces” off the massive effective potential barrier in a manner akin to that discovered in AdS. Similarly, in some region of the parameter space the evolution drives the system towards nonlinearly stable, asymptotically flat “oscillatons” [688 *, 326 ]. Nevertheless, for sufficiently small initial amplitudes they observe a $t−3∕2$ decay of the initial data, characteristic of massive fields, and showing that Minkowski is nonlinearly stable. The “weakly turbulent” instability discovered in AdS is stimulating research on the topic of turbulence in GR. A full understanding of the mechanism(s) will require further studies, including collapse in non-spherically symmetric backgrounds, other forms of matter and boundary conditions, etc.

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.

We now turn to solutions that display an instability at linear level, seen by a mode analysis, and to the use of NR techniques to follow the development of such instabilities into the nonlinear regime. One outstanding example is the Gregory–Laflamme instability of black strings already described in Section 7.2. Such black strings exist in higher dimensions, $D ≥ 5$ . It is expected that the same instability mechanism afflicts other higher-dimensional BHs, even with a topologically spherical horizon. A notable example are Myers–Perry BHs. In $D ≥ 6$ , the subset of these solutions with a single angular momentum parameter have no analogue to the Kerr bound, i.e., a maximum angular momentum for a given mass. As the angular momentum increases, they become ultra-spinning BHs and their horizon becomes increasingly flattened and hence resembling the horizon of a black $p$ -brane, which is subject to the Gregory–Laflamme instability [305 *]. It was indeed shown in Ref. [272 *, 271 , 270 ], by using linear perturbation theory, that rapidly rotating Myers–Perry BHs for $7 ≤ D ≤ 9$ are unstable against axisymmetric perturbations. The nonlinear growth of this instability is unknown but an educated guess is that it may lead to a deformation of the pancake like horizon towards multiple concentric rings.

A different argument – of entropic nature – for the instability of ultra-spinning BHs against non-axisymmetric perturbations was given by Emparan and Myers [305 *]. Such type of instability has been tested in $6 ≤ D ≤ 8$ [700 *], but also in $D = 5$ [701 *, 700 *] – for which a slightly different argument for instability was given in [305 ] – by evolving a Myers–Perry BH with a non-axisymmetric bar-mode deformation, using a NR code adapted to higher dimensions. In each case sufficiently rapidly rotating BHs are found to be unstable against the bar-mode deformation. In terms of a dimensionless spin parameter $q ≡ a∕μ1∕(D−3)$ , where $μ,a$ are the standard mass and angular momentum parameters of the Myers–Perry solution, the onset values for the instability were found to be: $D = 5, q = 0.87;D = 6, q = 0.74;D = 7, q = 0.73;D = 8, q = 0.77$ . We remark that the corresponding values found in [272 ] for the Gregory–Laflamme instability in $7 ≤ D ≤ 9$ are always larger than unity. Thus, the instability triggered by non-axisymmetric perturbations sets in for lower angular momenta than the axisymmetric Gregory–Laflamme instability. Moreover, in [700 *], long-term numerical evolutions have been performed to follow the nonlinear development of the instability. The central conclusion is that the unstable BHs relax to stable configurations by radiating away the excess angular momentum. These results have been confirmed for $D = 6,7$ by a linear analysis in Ref. [273 *]; such linear analysis suggests that for $D = 5$ , however, the single spinning Myers–Perry BH is linearly stable.²⁰

Results for the $D = 6$ “gravitational waveforms” $h+,×$ (see Eqs. (43) – (44) in Ref. [700 *] for a precise definition of these quantities) for an initial dimensionless $q = 0.801$ are shown in Figure 10 *. The early stage shows an exponential increase of the amplitude, after which a saturation phase is reached, and where angular momentum is being shed through GW emission. After this stage, exponential decay of the oscillations ensues. This bar-mode instability discovered – before linear perturbation analysis – in Refs. [701 , 700 *] has recently been seen at linear level for BHs having all angular momentum parameters equal in Refs. [395 , 310 ].

Another spacetime instability seen at linear level is the superradiant instability of rotating or charged BHs in the presence of massive fields or certain boundary conditions. This instability will be discussed in detail in the next section.

Concerning the nonlinear stability of BH solutions, the only generic statement one can produce at the moment is that hundreds of NR evolutions of binary Kerr or Schwarzschild BHs in vacuum, over the last decade, lend empirical support to the nonlinear stability of these solutions. One must remark, however, on the limitations of testing instabilities with NR simulations. For instance, fully nonlinear dynamical simulations cannot probe – at least at present – extremal Kerr BHs; they are also unable to find instabilities associated with very high harmonic indices $ℓ,m$ (associated with very small scales), as well as instabilities that may grow very slowly. Concerning the first caveat, it was actually recently found by Aretakis that extremal RN and Kerr BHs are linearly unstable against scalar perturbations [42 , 41 , 43 ], an observation subsequently generalised to more general linear fields [530 ] and to a nonlinear analysis [563 *, 44 ]. This growth of generic initial data on extremal horizons seems to be a very specific property of extremal BHs, in particular related to the absence of a redshift effect [563 ], and there is no evidence a similar instability occurs for non-extremal solutions.

To conclude this section let us briefly address the stability of BH interiors already discussed in Section 3.2.2. The picture suggested by Israel and Poisson [621 , 622 ] of mass inflation has been generically confirmed in a variety of toy models – i.e., not Kerr – by numerical evolutions [151 , 152 , 153 , 393 , 56 , 448 , 57 ] and also analytical arguments [233 ]. Other numerical/analytical studies also suggest the same holds for the realistic Kerr case [388 , 386 , 387 , 531 ]. As such, the current picture is that mass inflation will drive the curvature to Planckian values, near or at the Cauchy horizon. The precise nature of the consequent singularity, that is, if it is space-like or light-like, is however still under debate (see, e.g., [235 ]).

7.5 Superradiance and fundamental massive fields

There are several reasons to consider extensions of GR with minimally, or non-minimally coupled massive scalar fields with mass parameter $μS$ . As mentioned in Section 3.1.4, ultra-light degrees of freedom appear in the axiverse scenario [49 , 50 ] and they play an important role in cosmological models and also in dark matter models. Equally important is the fact that massive scalar fields are a very simple proxy for more complex, realistic matter fields, the understanding of which in full NR might take many years to achieve.

At linearized level, the behavior of fundamental fields in the vicinities of non-rotating BHs has been studied for decades, and the main features can be summarized as follows:

(i) A prompt response at early times, whose features depend on the initial conditions. This is the counterpart to light-cone propagation in flat space.

(ii) An exponentially decaying “ringdown” phase at intermediate times, where the BH is ringing in its characteristic QNMs. Bosonic fields of mass $μSℏ$ introduce both an extra scale in the problem and a potential barrier at distances $∼ 1∕μS$ , thus effectively trapping fluctuations. In this case, extra modes appear which are quasi-bound states, i.e., extremely long-lived states effectively turning the BH into a quasi-hairy BH [76 , 794 *, 280 *, 656 , 604 , 135 ].

(iii) At late times, the signal is dominated by a power-law fall-off, known as “late-time tail” [633 , 506 , 209 , 491 ]. Tails are caused by backscattering off spacetime curvature (and a potential barrier induced by massive terms) and more generically by a failure of Huygens’ principle. In other words, radiation in curved spacetimes travels not only on, but inside the entire light cone.

When the BH is rotating, a novel effect can be triggered: superradiance [827 *, 828 *, 83 *, 164 *]. Superradiance consists of energy extraction from rotating BHs, and a transfer of this energy to the interacting field [164 *]. For a monochromatic wave of frequency $ω$ , the condition for superradiance is [827 , 828 , 83 *, 164 *]

where $m$ is the azimuthal harmonic index and $ΩH$ is the angular velocity of the BH horizon. If, in addition, the field is massive, a “BH bomb-type” mechanism can ensue [241 , 268 , 843 , 171 *] leading to an instability of the spacetime and the growth of a scalar condensate outside the BH horizon [605 , 815 , 794 *, 280 , 588 *, 164 *]. The rich phenomenology of scenarios where fundamental fields couple to gravity motivated recent work on the subject, where full nonlinear evolutions are performed [588 *, 289 *, 410 *, 92 *]. East et al. [289 *] have performed nonlinear scattering experiments, solving the field equations in the generalized harmonic formulation, and constructing initial data representing a BH with dimensionless spin $a∕M = 0.99$ , and an incoming quadrupolar GW packet. Their results are summarized in Figure 11 *, for three different wavepacket frequencies, $M ω = 0.75,0.87,1$ (note that only the first is superradiant according to condition (173 *)). The wavepackets carry roughly 10% of the spacetime’s total mass. These results confirm that low frequency radiation does extract mass and spin from the BH (both the mass $MBH$ and spin $JBH$ of the BH decrease for the superradiant wavepacket with $M ω = 0.75$ ), and that nonlinear results agree quantitatively with linear predictions for small wavepacket amplitudes [745 ]. To summarize, superradiance is confirmed at full nonlinear level, providing a rigorous framework for the complex dynamics that are thought to arise for massive fields around rotating BHs.

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.

Self-interacting scalars can give rise to stable or very long-lived configurations. For example, self-interacting complex scalar fields can form boson stars for which the scalar field has an oscillatory nature, but the metric is stationary [458 , 685 , 516 *, 533 ]. Real-valued scalars can form oscillating solitons or “oscillatons”, long-lived configurations where both the scalar field and the metric are time-dependent [688 , 689 , 594 , 586 *]. Dynamical boson star configurations were studied by several authors [516 *], with focus on boson star collisions with different velocities and impact parameters [597 , 599 , 216 *]. These are important for tests of the hoop and cosmic censorship conjectures, and were reviewed briefly in Sections 7.2 and 7.3. For a thorough discussion and overview of results on dynamical boson stars we refer the reader to the Living Reviews article by Liebling and Palenzuela [516 ].

The first steps towards understanding the nonlinear interaction between massive fields and BHs were taken by Okawa et al. [588 *, 585 *], who found new ways to prescribe, and evolve, constraint-satisfying initial data, analytically or semi-analytically, for minimally coupled self-interacting scalar fields [588 *, 586 ]. This construction was reviewed in Section 6.3. In Ref. [588 *], the authors used this procedure to generate initial data and to evolve wavepackets of arbitrary angular shape in the vicinity of rotating BHs. Their results are summarized in Figure 12 *. Spherically symmetric initial data for massless fields reproduce previous results in the literature [382 ], and lead to power-law tails of integer index. The mass term adds an extra scale and a barrier at large distances, resulting in characteristic late-time tails of massive fields.

Full nonlinear results from Ref. [588 *] are reproduced in Figure 12 *, and agree with linearized predictions. Higher multipoles “feel” the centrifugal barrier close to the light ring which, together with the mass barrier at large distances, provides a confining mechanism and gives rise to almost stationary configurations, shown in the right panel of Figure 12 *. The beating patterns are a consequence of the excitation of different overtones with similar ringing frequency [794 , 588 *]. These “scalar condensates” are extremely long-lived and can, under some circumstances, be considered as adding hair to the BH. They are not however really stationary: the changing quadrupole moment of the “scalar cloud” triggers the simultaneous release of gravitational radiation [588 *, 816 , 585 *]. In fact, gravitational radiation is one of the most important effects not captured by linearized calculations. These nontrivial results extend to higher multipoles, which display an even more complex behavior [588 *, 585 ].

Although only a first step towards understanding the physics of fundamental fields in strong-field gravity, these results are encouraging. We expect that with more robust codes and longer simulations one will be able to fully explore the field, in particular, the following features.

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.

⋅

Superradiant instability and its saturation. The timescales probed in current nonlinear simulations are still not sufficient to unequivocally observe superradiance with test scalar fields. The main reason for this is the feebleness of such instabilities: for scalar fields they have, at best, an instability timescale of order $7 10 M$ for carefully tuned scalar field mass. However, current long-term simulations are able to extract GWs induced by the scalar cloud [588 ].

The biggest challenge ahead is to perform simulations which are accurate enough and last long enough to observe the scalar-instability growth and its subsequent saturation by GW emission. This will allow GW templates for this mechanism to finally be released.

Due to their simplicity, scalar fields are a natural candidate to carry on this program, but they are not the only one. Massive vector fields, which are known to have amplification factors one order of magnitude larger, give rise to stronger superradiant instabilities, and might also be a good candidate to finally observe superradiant instabilities at the nonlinear level. We note that the development of the superradiant instability may, in some special cases, lead to a truly asymptotically flat, hairy BH solution of the type recently discussed in [422 *].

⋅

Turbulence of massive fields in strong gravity. Linearized results indicate that the development of superradiant instabilities leaves behind a scalar cloud with scalar particles of frequency $ω ∼ m Ω$ , in a nearly stationary state. This system may therefore be prone to turbulent effects, where nonlinear terms may play an important role. One intriguing aspect of these setups is the possibility of having gravitational turbulence or collapse on sufficiently large timescales. Such effects were recently observed in “closed” systems where scalar fields are forced to interact gravitationally for long times [108 , 537 , 145 ]. It is plausible that quasi-bound states are also prone to such effects, but in asymptotically flat spacetimes.

⋅

Floating orbits. Our discussion until now has focused on minimally coupled fundamental fields. If couplings to matter exist, new effects are possible: a small object (for example, a star), orbiting a rotating, supermassive BH might be able to extract energy and angular momentum from the BH and convert it to gravitational radiation. For this to happen, the object would effectively stall at a superradiant orbit, with a Newtonian frequency $Ω = 2μS$ for the dominant quadrupolar emission. These are called floating orbits, and were verified at linearized level [165 , 164 ].²¹ Nonlinear evolutions of systems on floating orbits are extremely challenging on account of all the different extreme scales involved.

⋅

Superradiant instabilities in AdS. The mechanism behind superradiant instabilities relies on amplification close to the horizon and reflection by a barrier at large distances. Asymptotically AdS spacetimes provide an infinite-height barrier, ideal for the instability to develop [171 , 166 , 172 *, 755 , 169 *].

Because in these backgrounds there is no dissipation at infinity, it is both possible and likely that new, non-symmetric final states arise as a consequence of the superradiant instability [172 , 168 , 514 , 275 *]. Following the instability growth and its final state remains a challenge for NR in asymptotically AdS spacetimes.

⋅

Superradiant instabilities of charged BHs. Superradiant amplification of charged bosonic fields can occur in the background of charged BHs, in quite a similar fashion to the rotating case above, as long as the frequency of the impinging wave $ω$ obeys

where $q$ is the charge of the field and $ΦH$ is the electric potential on the BH horizon. In this case both charge and Coulomb energy are extracted from the BH in a way compatible with the first and second law of BH thermodynamics [83 ]. In order to have a recurrent scattering, and hence, an instability, it is not enough, however, to add a mass term to the field [339 , 432 , 430 , 261 , 671 ]; but an instability occurs either by imposing a mirror like boundary condition at some distance from the BH (i.e., a boxed BH) or by considering an asymptotically AdS spacetime. In Refs. [421 , 262 ] it has been established, through both a frequency and a time domain analysis, that the time scales for the development of the instability for boxed BHs can be made much smaller than for rotating BHs (in fact, arbitrarily small [431 ]). Together with the fact that even $s$ -waves, i.e., $ℓ = 0$ modes, can trigger the instability in charged BHs, makes the numerical study of the nonlinear development of this type of superradiant instabilities particularly promising. One should be aware, however, that there may be qualitative differences in both the development and end-point of superradiant instabilities in different setups. For instance, for AdS and boxed BHs, the end-point is likely a hairy BH, such as those constructed in Ref. [275 ] (for rotating BHs), since the scalar field cannot be dissipated anywhere. This applies to both charged and rotating backgrounds. By contrast, this does not apply to asymptotically flat spacetimes, wherein rotating (but not charged) superradiance instability may occur. It is then an open question if the system approaches a hairy BH – of the type constructed in [422 ] – or if the field is completely radiated/absorbed by the BH. Concerning the development of the instability, an important difference between the charged and the rotating cases may arise from the fact that a similar role, in Eqs. (173 *) and (174), is played by the field’s azimuthal quantum number $m$ and the field charge $q$ ; but whereas the former may change in a nonlinear evolution, the latter is conserved [169 ].

7.6 High-energy collisions

Applications of NR to collisions of BHs or compact matter sources near the speed of light are largely motivated by probing GR in its most violent regime and by the modelling of BH formation in TeV gravity scenarios. The most important questions that arise in these contexts can be summarized as follows.

Does cosmic censorship still apply under the extreme conditions of collisions near the speed of light? As has already been discussed in Section 7.2, numerical simulations of these collisions in four dimensions have so far identified horizon formation in agreement with the censorship conjecture. The results of higher-dimensional simulations are still not fully understood, cf. Section 7.2.
Do NR simulations of high-energy particle collisions provide evidence supporting the validity of the hoop conjecture? As discussed in Section 7.3, NR results have so far confirmed the hoop conjecture.
In collisions near the speed of light, the energy mostly consists of the kinetic energy of the colliding particles such that their internal structure should be negligible for the collision dynamics. Furthermore, the gravitational field of a particle moving at the speed of light is non-vanishing only near the particle’s worldline [16 ], suggesting that the gravitational interaction in high-energy collisions should be dominant at the instant of collision and engulfed inside the horizon that forms. This conjecture has sometimes been summarized by the statement that “matter does not matter” [216 ], and is related to the hoop conjecture discussed above. Do NR simulations of generic high-energy collisions of compact objects support this argument in the classical regime, i.e., does the modelling of the colliding objects as point particles (and, in particular, as BHs) provide an accurate description of the dynamics?
Assuming that the previous question is answered in the affirmative, what is the scattering threshold for BH formation? This corresponds to determining the threshold impact parameter $bscat$ that separates collisions resulting in the formation of a single BH ( $b < bscat$ ) from scattering encounters ( $b > bscat$ ), as a function of the number of spacetime dimensions $D$ and the collision velocity $v$ in the center-of-mass frame or boost parameter $√ -----2 γ = 1∕ 1 − v$ .
How much energy and momentum is lost in the form of GWs during the collision? By conversion of energy and momentum, the GW emission determines the mass and spin of the BH (if formed) as a function of the spacetime dimension $D$ , scattering parameter $b$ , and boost factor $γ$ of the collision. Collisions near the speed of light are also intriguing events to probe the extremes of GR; in particular what is the maximum radiation that can be extracted from any collision and does the luminosity approach Dyson’s limit $dE ∕dt ≲ 1$ [157 ]? (See discussion in Section 3.2.3 about this limit.)

These issues are presently rather well understood through NR simulations in $D = 4$ spacetime dimensions but remain largely unanswered for the important cases $D ≥ 5$ .

The relevance of the internal structure of the colliding bodies has been studied in Ref. [716 *], comparing the GW emission and scattering threshold in high-energy collisions of rotating and non-rotating BHs in $D = 4$ . The BH spins of the rotating configurations are either aligned or anti-aligned with the orbital angular momentum corresponding to the so-called hang-up and anti-hang-up cases which were found to have particularly strong effects on the dynamics in quasi-circular BH binary inspirals [160 ]. In high-energy collisions, however, this (anti-)hang-up effect disappears; the GW emission as well as the scattering threshold are essentially independent of the BH spin at large collision velocities (cf. Figure 15 * which will be discussed in more detail further below). These findings suggest that ultra-relativistic collisions are indeed well modelled by colliding point-particles or BHs in GR. In the center-of-mass frame, and assuming that the two particles have equal mass, the collisions are characterized by three parameters. (i) The number $D$ of spacetime dimensions, (ii) the Lorentz factor $γ$ or, equivalently, the collision velocity $v$ , and (iii) the impact parameter $b = L∕P$ , where $L$ and $P$ are the initial orbital angular momentum and the linear momentum of either BH in the center-of-mass frame.

The simplest set of configurations consists of head-on collisions with $b = 0$ in $D = 4$ dimensions and was analysed in Ref. [719 ] varying the boost parameter in the range $1.07 ≤ γ ≤ 3$ . In agreement with the cosmic censorship conjecture, these collisions always result in the formation of a single BH that settles into a stationary configuration through quasi-normal ringdown. The total energy radiated in the form of GWs is well modelled by the following functional form predicted by Smarr’s [707 ] zero-frequency limit (see Section 5.3)

where $E∞$ is a free parameter that corresponds to the fraction of energy radiated in the limit $γ → ∞$ . Fitting the numerical results with Eq. (175 *) yields $E ∞ = 14 ± 3 %$ which is about half of Penrose’s upper limit [611 , 286 ]. Observe the good agreement with the second order result in Eq. (38 *) as discussed in Section 5.4.

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.

Grazing collisions in four dimensions represent two-parameter studies, where the boost factor $γ$ and the impact parameter $b$ are varied, and have been investigated in Refs. [697 *, 720 *]. At fixed Lorentz boost, such grazing collisions exhibit three distinct regimes as the impact parameter is increased from the head-on limit $b = 0$ : (i) prompt mergers, (ii) delayed mergers, and (iii) the scattering regime where no common horizon forms. These regimes are marked by two special values of the impact parameter $b$ , the scattering threshold $bscat$ that we have already mentioned above and the threshold of immediate merger $b∗$ . This threshold has been identified in numerical BH simulations by Pretorius & Khurana [632 ] as marking the onset of a regime where the two BHs whirl around each other prior to merging or scattering off for a number of orbits proportional to $log|b − b∗|$ ; see also [413 , 355 ]. This zoom-whirl-like behaviour has also been identified in high-energy grazing collisions in [720 ]. The three different regimes are illustrated in Figure 13 * which shows the BH trajectories for $γ = 1.520$ and $b∕M = 3.24$ , $3.29$ and $3.45$ . For this boost factor, the thresholds are given by $b∗∕M ≈ 3.25$ and $bscat∕M ≈ 3.35$ . For impact parameters close to the threshold values $b∗$ and $bscat$ , grazing collisions can generate enormous amounts of GWs. This is shown in the left panel of Figure 14 *. Starting from $E ∕M ≈ 2.2%$ in the head-on limit $b = 0$ , the radiated energy increases by more than an order of magnitude to $≳ 25%$ for $∗ b < b < bscat$ . These simulations can also result in BHs spinning close to the extremal Kerr limit as is shown in the right panel of the figure which plots the dimensionless final spin as a function of $b$ . By fitting their numerical results, Shibata et al. [697 ] have found an empirical relation for the scattering threshold given by

For such a value of the impact parameter, they observe that the dimensionless final spin of the merged BH is given by $J fin = (0.6 ± 0.1)J$ , where $J ∕M 2 = 1.25 ± 0.03$ is the initial angular momentum of the system.

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.

Grazing collisions of spinning and non-spinning BHs have been compared in Ref. [716 ]. The initial configurations for these simulations have been chosen with $γ$ factors up to 2.49 and equal spin for both BHs of dimensionless magnitude $χ = 0.85$ and 0.65 aligned or anti-aligned with the orbital angular momentum $L$ . This set has been complemented with collisions of non-spinning BHs covering the same range in $γ$ . The scattering threshold and the energy radiated in GWs in these simulations are shown in the left panel of Figure 15 *. As expected from the hang-up effect, aligned (anti-aligned) spins result in a smaller (larger) value of the scattering threshold $b scat$ at low collision speeds. At velocities above $∼ 80%$ of the speed of light, however, this effect is washed out and, in agreement with the matter-does-not-matter conjecture mentioned above, the collision dynamics are barely affected by the BH spins. Furthermore, the scattering threshold determined for non-spinning BHs agrees very well with the formula (176 *). As demonstrated in the bottom panel of the figure, the radiated energy is barely affected by the BH spin even in the mildly relativistic regime. The simulations also suggest an upper limit of the fraction of kinetic energy that can be converted into GWs. Extrapolation of the data points in Figure 15 * to $v = 1$ predicts that at most about half of the total energy can be dissipated in GWs in any four dimensional collision. The other half, instead, ends up as rest mass inside the common horizon formed in merging configurations or is absorbed by the individual BHs during the close encounter in scattering processes (the result of this extrapolation is consistent with the calculation in Ref. [376 ]). This is illustrated in the right panel of Figure 15 *, where the trajectory of one BH in a delayed merger configuration with anti-aligned spins is shown. The circles, with radius proportional to the horizon mass, represent the BH location at intervals $Δt = 10M$ . During the close encounter, (i) the BH grows in size due to absorption of gravitational energy and (ii) slows down considerably.

Collisions of BHs with electric charge have been simulated by Zilhão et al. [838 , 839 *]. For the special case of BHs with equal charge-to-mass ratio $Q ∕M$ and initially at rest, constraint-satisfying initial data are available in closed analytic form. The electromagnetic wave signal generated in these head-on collisions reveals three regimes similar to the pattern known for the GW signal, (i) an infall phase prior to formation of a common horizon, (ii) the nonlinear merger phase where the wave emission reaches its maximum and (iii) the quasi-normal ringdown. As the charge-to-mass ratio is increased towards $Q ∕M ≲ 1$ , the emitted GW energy decreases by about 3 orders of magnitude while the electromagnetic wave energy reaches a maximum at $Q ∕M ≈ 0.6$ , and drops towards $0$ in both the uncharged and the extreme limit. This behaviour of the radiated energies is expected because of the decelerating effect of the repulsive electric force between equally charged BHs. For opposite electric charges, on the other hand, the larger collision velocity results in an increased amount of GWs and electromagnetic radiation [839 ].

An extended study of BH collisions using various analytic approximation techniques including geodesic calculations and the ZFL has been presented in Berti et al. [93 ]; see also [94 *] for a first exploration in higher dimensions. Weak scattering of BHs in $D = 4$ , which means large scattering parameters $b∕M ∼ 10$ , and for velocities $v ≈ 0.2$ , has been studied by Damour et al. [245 ] using NR as well as PN and EOB calculations. Whereas PN calculations start deviating significantly from the NR results for $b∕M ≲ 10$ , the NR calibrated EOB model yields good agreement in the scattering angle throughout the weak scattering regime.

BH collisions in $D ≥ 5$ spacetime dimensions are not as well understood as their four-dimensional counterparts. This is largely a consequence of the fact that NR in higher dimensions is not yet that robust and suffers more strongly from numerical instabilities. Such complications in the higher-dimensional numerics do not appear to cause similar problems in the construction of constraint satisfying initial data. The spectral elliptic solver originally developed by Ansorg et al. [39 ] for $D = 4$ has been successfully generalized to higher $D$ in Ref. [836 ] and provides solutions with comparable accuracy as in $D = 4$ .

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.

A systematic exploration of the scattering threshold in $D = 5$ dimensions has been performed in Ref. [587 *]. By superposing non-rotating, boosted single BH initial data, they have evolved grazing collisions up to $v ≲ 0.8$ . Their results are summarized in Figure 16 *, where scattering (merging) BH collisions are marked by “plus” and “circle” symbols in the plane spanned by the collision velocity $v$ and the impact parameter $b$ . The simulations show a decrease of the scattering threshold at increasing velocity up to $v ≈ 0.6$ , similar to the $D = 4$ case in the upper left panel of Figure 15 *. At larger $v$ , the threshold cannot yet be determined because simulations with near critical impact parameter become numerically unstable. By monitoring the Kretschmann scalar at the point of symmetry between the two BHs, a large curvature regime was furthermore identified in [587 *], as discussed in Section 7.2.

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 *].

The GW emission in BH mergers in $D > 4$ has so far only been studied for collisions starting from rest. In Refs. [797 , 793 , 796 *] the wave signal was extracted using the Kodama-Ishibashi formalism discussed in Section 6.7.2; the GW emission contains about 0.09%, 0.08% of the center-of-mass energy for equal-mass binaries in $D = 5,6$ respectively (note that in $D = 4$ only 0.055% of the center-of-mass-eanergy goes into GWs), and decreases with the mass ratio. The dependency of the radiated energy and momentum on the mass ratio is well modelled by point particle calculations [94 ]. A comparison of the predicted GW emission in higher-dimensional collisions using two different numerical codes with different formulations of the Einstein equations, namely those discussed in Sections 6.2.1 and 6.2.2, has been presented in Ref. [796 *]. The predictions from the two codes using respectively the Kodama–Ishibashi formalism (cf. Section 6.7.2) and a direct extraction through the metric components (cf. Section IV B 1 in [700 ]) agree within numerical uncertainties [796 ]. This result, illustrated in Figure 17 *, represents an important validation of both the numerical evolution techniques and the diagnostics of the simulations, along with the first estimate of the emitted energy for head-on collisions in $D = 6$ .

The main challenges for future numerical work in the field of high-energy collisions are rather evident. For applications in the analysis of experimental data in the context of TeV gravity scenarios (cf. Section 3.3.2), it will be vital to generalize the results obtained in four dimensions to $D ≥ 5$ . Furthermore, it is currently not known whether the impact of electric charge on the collision dynamics becomes negligible at high velocities, as suggested by the matter-does-not-matter conjecture, and as is the case for the BH spin.

7.7 Alternative theories

As discussed in Section 6.1.7, one of the most straightforward extensions of Einstein’s theory is obtained by the addition of minimally coupled scalar fields. When the scalar couples to the Ricci scalar however, one gets a modification of Einstein gravity, called scalar-tensor theory. In vacuum, scalar-tensor theories are described by the generic action in Eq. (5 *), where $R$ is the Ricci scalar associated to the metric $gμν$ , and $F (ϕ),Z (ϕ)$ and $U (ϕ)$ are arbitrary functions (see e.g., [92 *, 824 *] and references therein). The matter fields minimally coupled to $gμν$ are collectively denoted by $Ψm$ . This form of the action corresponds to the choice of the so-called “Jordan frame”, where the matter fields $Ψm$ obey the equivalence principle. For $F = ϕ,Z = ωBD ∕ϕ,U = 0$ , the action (5 *) reduces to the standard Brans–Dicke theory.

The equations of motion derived from the action (5 *) are second-order and the theory admits a well-posed initial-value problem [670 *]. These facts turn scalar-theories into an attractive alternative to Einstein’s equations, embodying at least some of the physics one expects from an ultimate theory of gravity, and have been a major driving force behind the efforts to understand scalar-tensor theories from a NR point of view [696 , 679 , 680 , 582 , 410 *, 92 *, 73 *, 670 ]. In fact, scalar-tensor theories remain the only alternative theory to date where full nonlinear dynamical evolutions of BH spacetimes have been performed.

Scalar-tensor theories can be recast in such a way as to be formally equivalent, in vacuum, to GR with a minimally coupled scalar field, i.e., to the theory described previously in Section 6.1.7. This greatly reduces the amount of work necessary to extend NR to these setups. The explicit transformations that recast the previous action in the “Einstein frame” are [244 ]

U gEμν = F(ϕ )gμν, V = ----------, (177 ) 16πF (ϕ )2 1 ∫ [3 F′(ϕ)2 4πGZ (ϕ)]1∕2 φ (ϕ) = √------ dϕ -------2-+ --------- . (178 ) 4πG 4 F (ϕ) F (ϕ)

The Einstein-frame action is then

which is the action for a minimally coupled field (77 *) enlarged to allow for a generic self-interaction potential (which could include the mass term). The label $E$ denotes quantities constructed from the Einstein-frame metric $gEμν$ . In the Einstein frame the scalar field is minimally coupled to gravity, but any matter field $Ψm$ is coupled to the metric $F− 1gEμν$ .

In vacuum, this action leads to the following equations of motion:

1 REμν − -gEμνRE − 8πT μν = 0, (180 ) 2 ′ ∇ μ∇ μφ − V (φ) = 0, (181 )

where the energy-momentum tensor of the scalar field is determined by

In summary, the study of scalar-tensor theories of gravity can be directly translated, in vacuum, to the study of minimally coupled scalar fields. For a trivial potential $V = const$ , the equations of motion in the Einstein frame (180 *) – (181 *) admit GR (with $φ = const$ ) as a solution. Because stationary BH spacetimes in GR are stable, i.e., any scalar fluctuations die away rather quickly, the dynamical evolution of vacuum BHs is expected to be the same as in GR. This conclusion relies on hand-waving stability arguments, but was verified to be true to first PN order by Will and Zaglauer [787 ], at 2.5 PN order by Mirshekari and Will [549 ] and to all orders in the point particle limit in Ref. [824 ].

Thus, at least one of the following three ingredients are necessary to generate interesting dynamics in scalar-tensor theories:

⋅: Nontrivial potential $V$ and initial conditions. Healy et al. studied an equal-mass BH binary in an inflation-inspired potential $2 V = λ (φ2 − φ20) ∕8$ with nontrivial initial conditions on the scalar given by $φ = φ0 tanh (r − r0)∕ σ$ [410 *]. This setup is expected to cause deviations in the dynamics of the inspiralling binary, because the binary is now accreting scalar field energy. The larger the initial amplitude of the field, the larger those deviations are expected to be. This is summarized in Figure 18 *, where the BH positions are shown as a function of time for varying initial scalar amplitude.

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.
⋅: Nontrivial boundary conditions. As discussed, GR is recovered for constant scalar fields. For nontrivial time-dependent boundary conditions or background scalar fields, however, nontrivial results show up. These boundary conditions could mimic cosmological scenarios or dark matter profiles in galaxies [434 , 92 *]. Reference [92 *] modelled a BH binary evolving nonlinearly in a constant-gradient scalar field. The scalar-field gradient induces scalar charge on the BHs, and the accelerated motion of each BH in the binary generates scalar radiation at large distances, as summarized in Figure 19 *.

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.

The scalar-signal at large distances, shown in the right panel of Figure 19 *, mimics the inspiral, merger and ringdown stages in the GW signal of an inspiralling BH binary.
⋅: Matter. When matter is present, new effects (due to the coupling of matter to the effective metric $F −1gEμν$ ) can dominate the dynamics and wave emission. For example, it has been shown that, for $β ≡ ∂2φ (ln F(φ )) ≲ − 4$ , NSs can “spontaneously scalarize,” i.e., for sufficiently large compactnesses the GR solution is unstable. The stable branch has a nonzero expectation value for the scalar field [242 ].

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.

Scalarized matter offers a rich new phenomenology. For example, the dynamics and GW emission of scalarized NSs can be appreciably different (for given coupling function $− 1 F (ϕ )$ ) from the corresponding GR quantities, as shown by Palenzuela et al. [73 *, 596 *] and summarized in Figure 20 *. Strong-field gravity can even induce dynamical scalarization of otherwise GR stars during inspiral, offering new ways to constrain such theories [73 , 596 ].

The application of NR methods to the understanding of alternative theories of gravity and tests of GR is still in its infancy. Among various possible directions, we point out the following.

⋅: Understanding the well posedness of some theory, in particular those having some motivation from fundamental physics, as for example Einstein-Dilaton-Gauss–Bonnet and Dynamical Chern–Simons gravity [602 , 27 ]. A study on the well posedness of the latter has recently been presented in Ref. [263 ].
⋅: Building initial data describing interesting setups for such theories. Unless the theory admits particularly simple analytic solutions, it is likely that initial data construction will also have to be done numerically. Apart from noteworthy exceptions, such as Gauss–Bonnet gravity in higher dimensions [814 ], initial data have hardly been considered in the literature.

Once well-posedness is established and initial data are constructed, NR evolutions will help us understanding how these theories behave in the nonlinear regime.

7.8 Holography

Holography provides a fascinating new source of problems for NR. As such, in recent years, a number of numerical frameworks have been explored in asymptotically AdS spacetimes, as to face the various pressing questions raised within the holographic correspondence, cf. Section 3.3.1. At the moment of writing, no general purpose code has been reported, comparable to existing codes in asymptotically flat spacetime, which can evolve, say, BH binaries with essentially arbitrary masses, spins and momenta. Progress has occured in specific directions to address specific issues. We shall now review some of these developments emphasizing the gravity side of the problems.

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.

As mentioned in Section 3.3.1 an important problem in the physics of heavy ion collisions is to understand the “early thermalization problem”. In Ref. [207 *], the gauge/gravity duality was used to address this issue. On the gravity side, the problem at hand was to study a head-on collision of two shock waves in asymptotically AdS₅ spacetime. The numerical scheme was to perfom a null (characteristic) evolution. By choosing a specific metric ansatz, it was possible to unveil in the nonlinear Einstein equations a nested linear structure: the equations can be integrated as linear ordinary differential equations if an appropriate sequence is chosen. The AdS boundary condition was implemented by an adequate radial expansion near the boundary and the initial data consisted of two well-separated planar shocks, with finite thickness and energy density, moving toward each other. In this setup an AH is always found (even before the collision) and excision was performed by restricting the computational domain to start at this horizon. The evolution of the two shock waves is displayed in Figure 21 * (left panel). By following the evolution and using the gauge/gravity dictionary, the authors reported that the total time required for apparent thermalization was 0.35 fm/c. This is within the same order of magnitude as the thermalization scale obtained from accelerator data, already discussed in Section 3.3.1. A discussion of numerical approaches using null evolutions applied to asymptotically AdS problems can be found in [208 *]. Other recent applications of shock wave collisions in AdS₅ to describe phenomenological properties of heavy ion collisions can be found in Refs. [800 , 655 , 757 ].

Time-plus-space decompositions have also been initiated, both based on a generalized harmonic evolution scheme [70 *] and in an ADM formulation [416 ]. In particular the latter formulation seems very suited for extracting relevant physical quantities for holography, such as the boundary time for the thermalization process discussed in Section 3.3.1.

Evolutions of BHs deformed by a scalar field in AdS₅ have been presented in Ref. [70 ]. The evolution leads the system to oscillate in a (expected) superposition of quasi-normal modes, some of which are nonlinearly driven. On the boundary, the dual CFT stress tensor behaves like that of a thermalized $𝒩 = 4$ super-Yang–Mills fluid, with an equation of state consistent with conformal invariance and transport coefficients that match holographic calculations at all times. Similar conclusions were reached in Ref. [417 ], where the numerical scheme of Ref. [207 *, 208 *] was used to study the isotropization of a homogeneous, strongly coupled, non-Abelian plasma by means of its gravity dual, comparing the time evolution of a large number of initially anisotropic states. They find that the linear approximation seems to work well even for initial states with large anisotropies. This unreasonable effectiveness of linearized predictions hints at something more fundamental at work, perhaps a washing out of nonlinearities close to the horizon. Such effects were observed before in asymptotically flat spacetimes, for example the already mentioned agreement between ZFL (see Section 5.3) or close limit approximation predictions (see Section 5.2) and full nonlinear results.

Also of interest for accelerator physics, and the subject of intense work in recent years, are holographic descriptions of jet-quenching, i.e., the loss of energy of partons as they cross strongly coupled plasmas produced in heavy ion collisions [2 , 199 , 200 ]. Numerical work using schemes similar to that of Refs. [207 , 208 ] have been used to evolve dual geometries describing the quenches [143 , 144 , 204 ]; see also [319 ] for numerical stationary solutions in this context.

Another development within the gauge/gravity duality that gained much attention, also discussed in Section 3.3.1, is related to condensed matter physics. In asymptotically AdS spacetimes, a simple theory, say, with a scalar field minimally coupled to the Maxwell field and to gravity admits RN-AdS as a solution. Below a critical temperature, however, this solution is unstable against perturbations of the scalar field, which develops a tachyonic mode. Since the theory admits another set of charged BH solutions, which have scalar hair, it was suggested that the development of the instability of the RN-AdS BHs leads the system to a hairy solution. From the dual field theory viewpoint, this corresponds to a phase transition between a normal and a superconducting phase. A numerical simulation showing that indeed the spacetime evolution of the unstable RN-AdS BHs leads to a hairy BH was reported in [562 ]. Therein, the authors performed a numerical evolution of a planar RN-AdS BH perturbed by the scalar field and using Eddington–Finkelstein coordinates. A particular numerical scheme was developed, adapted to this problem. The development of the scalar field density is shown in Figure 21 *. The initial exponential growth of the scalar field is eventually replaced by an approach to a fixed value, corresponding to the value of the scalar condensate on the hairy BH.

Finally, the gauge/gravity duality itself may provide insight into turbulence. Turbulent flows of CFTs are dual to dynamical BH solutions in asymptotically AdS spacetimes. Thus, urgent questions begging for answers include how and when do turbulent BHs arise, and what is the (gravitational) origin of Kolmogorov scaling observed in turbulent fluid flows. These problems are just now starting to be addressed [185 , 12 , 13 ].

7.9 Applications in cosmological settings

Some initial applications of NR methods addressing specific issues in cosmology have been reviewed in the Living Reviews article by Anninos [36 ], ranging from the Big Bang singularity dynamics to the interactions of GWs and the large-scale structure of the universe. The first of these problems – the understanding of cosmological singularities – actually motivated the earlier applications of NR to cosmological settings, cf. the Living Reviews article [88 *]. The set of homogeneous but anisotropic universes was classified by Bianchi in 1898 into nine different types (corresponding to different independent groups of isometries for the 3-dimensional space). Belinskii, Khalatnikov and Lifshitz (BKL) proposed that the singularity of a generic inhomogeneous cosmology is a “chaotic” spacelike curvature singularity, and that it would behave asymptotically like a Bianchi IX or VIII homogeneous cosmological model. This is called BKL dynamics or mixmaster universe. The accuracy of the BKL dynamics has been investigated using numerical evolutions in Refs. [87 , 558 , 662 ], and the BKL sensitivity to initial conditions in various references (see for instance Ref. [89 ]). For further details, we refer the reader to Ref. [88 ].

More recently, NR methods have been applied to the study of bouncing cosmologies, by studying the evolution of adiabatic perturbations in a nonsingular bounce [801 *]. The results of Ref. [801 ] show that the bounce is disrupted in regions of the universe with significant inhomogeneity and anisotropy over the background energy density, but is achieved in regions that are relatively homogeneous and isotropic. Sufficiently small perturbations, consistent with observational constraints, can pass through the nonsingular bounce with negligible alteration from nonlinearity.

In parallel, studies of “bubble universes”, in which our universe is one of many nucleating and growing inside an ever-expanding false vacuum, have also been made with NR tools. In particular, Refs. [765 , 764 ] investigated the collisions between bubbles, by computing the cosmological observables arising from bubble collisions directly from the Lagrangian of a single scalar field.

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.

Applications of NR in more standard cosmological settings are still in their infancy, but remarkable progress has been achieved. One of these concerns the impact of cosmic inhomogeneities on the value of the cosmological constant and the acceleration of the universe. In other words, how good are models of homogeneous and isotropic universes – the paradigmatic Friedmann–Lemaître–Robertson–Walker (FLRW) geometry – when we know that our universe has structure and is inhomogeneous?

Studies of this (long-standing, see for instance Ref. [520 ]) question within NR have considered the evolution of BH lattices (the BHs mimicking strong, self-gravitating inhomogeneities) [86 *, 805 *]. In Ref. [86 *] the authors explicitly constructed and evolved a three-dimensional, fully relativistic, eight-BH lattice with the topology of $S3$ . The puncture locations in that work projected down to $R3$ are shown in the left panel of Figure 22 * (one of the punctures is projected out to infinity, see Ref. [86 *] for further details). The evolution of this 8-BH configuration is summarized in the right panel of Figure 22 *, showing the (minimal) proper distance between neighbouring surfaces and the proper length of each cell’s edges. These quantities are then compared against a reference FLRW closed model with spatial slices of spherical topology. The comparison procedure is not straightforward, but adopting the procedure of Ref. [86 ] it yields good agreement.

The effects of local inhomogeneities have been investigated in Ref. [804 *] using different initial data, describing an expanding inhomogeneous universe model composed of regularly aligned BHs of identical mass. The evolution of these initial data also indicates that local inhomogeneities do not significantly affect the global expansion law of the universe, despite the fact that the inhomogeneities themselves are extremely nonlinear [804 , 805 ]. Similar conclusions were reached in Ref. [834 ], where the ADM formalism is used to develop a practical scheme to calculate a proposed domain averaging effect in an inhomogeneous cosmology within the context of numerical large-scale structure simulations. This study finds that in the weak-field, slow-motion limit, the proposed effect implies a small correction to the global expansion rate of the universe. In this limit, their simulations are always dominated by the expanding underdense regions, hence the correction to the energy density is negative and the effective pressure is positive. The effects of strong gravity in more general scenarios are yet to be understood [834 ]. For an earlier NR code developed to address inhomogeneous cosmologies see Ref. [426 ].

More complex NR codes aimed at understanding cosmological evolutions are currently being developed. NR simulations of large scale dynamical processes in the early universe have recently been reported [345 ]. These take into account interactions of dark matter, scalar perturbations, GWs, magnetic fields and turbulent plasma. Finally, Ref. [837 ] considers the effect of (extreme) cosmological expansion on the head-on collision and merger of two BHs, by modelling the collision of BHs in asymptotically dS spacetimes.

	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