. 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 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,
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 . For large amplitudes , 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 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
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  or different gravitational theories . 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. , whereas recently the first 3 + 1 simulations of the collapse of minimally coupled scalar fields were reported . 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 .
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 processes19 – 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 . 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  have shown that the 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  but these processes seem to require fine-tuned initial conditions, such as in spherically symmetric collapse or in the critical collapse  discussed in Section 7.1.
A claim of cosmic censorship violation in 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 ; therefore, fine-tuning of initial conditions may be required for the formation of a naked singularity.
Another suggestion that Planckian scale curvature becomes visible in a classical evolution in 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 were reported. For sufficiently small initial velocities () it is possible to find the threshold impact parameter such that the BHs merge into a (spinning) BH for or scatter off to infinity for . For high velocities, however, only a lower bound on the impact parameter for scattering and an upper bound on the impact parameter for merger could be found, since simulations with crashed before the final outcome could be determined (cf. Figure 16* below). Moreover, an analysis of a scattering configuration with and , 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 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 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., states that when the mass of a system (in dimensions) gets compacted into a region whose circumference in every direction has radius , a horizon – and thus a BH – forms (for a generalization in , see ). 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 and radius (corresponding to the de Broglie wavelength of the process) collide with a boost parameter , the mass-energy in the centre-of-mass frame is . The threshold radius of Thorne’s hoop is then and the condition translates into a bound on the boost factor, .
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 (where the “hoop” critical boost factor 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 . Here, the fluid balls are modeled as two superposed Tolman–Oppenheimer–Volkoff “stars” with a polytropic equation of state.
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 that does not result in horizon formation (upper panels) and one at 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  has shown similar results. Therein it has been found that a BH forms when the compactness of the star is , i.e., for . Type-I critical behaviour has also been identified, with BH formation for initial masses above a critical value scaling as .
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 . For BH spacetimes, however, mode stability does not guarantee linear stability, cf. the discussion in . 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 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 : , where and 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  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 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 . 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 . 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 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.
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, . 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 , 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 -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 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 [700*], but also in [701*, 700*] – for which a slightly different argument for instability was given in  – 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 , where are the standard mass and angular momentum parameters of the Myers–Perry solution, the onset values for the instability were found to be: . We remark that the corresponding values found in  for the Gregory–Laflamme instability in 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 by a linear analysis in Ref. [273*]; such linear analysis suggests that for , however, the single spinning Myers–Perry BH is linearly stable.20
Results for the “gravitational waveforms” (see Eqs. (43) – (44) in Ref. [700*] for a precise definition of these quantities) for an initial dimensionless 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 (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  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 , 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 . 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., ).
There are several reasons to consider extensions of GR with minimally, or non-minimally coupled massive scalar fields with mass parameter . 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 introduce both an extra scale in the problem and a potential barrier at distances , 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*][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 , and an incoming quadrupolar GW packet. Their results are summarized in Figure 11*, for three different wavepacket frequencies, (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 and spin of the BH decrease for the superradiant wavepacket with ), and that nonlinear results agree quantitatively with linear predictions for small wavepacket amplitudes . 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.
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 .
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 , 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.
- 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 for carefully tuned scalar field mass. However, current
long-term simulations are able to extract GWs induced by the scalar cloud .
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 , 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 for the dominant quadrupolar emission. These are called floating orbits, and were verified at linearized level [165, 164].21 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 is the charge of the field and 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 . 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 ). Together with the fact that even -waves, i.e., 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.  (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  – 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 and the field charge ; but whereas the former may change in a nonlinear evolution, the latter is conserved .
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 , 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” , 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 that separates collisions resulting in the formation of a single BH () from scattering encounters (), as a function of the number of spacetime dimensions and the collision velocity in the center-of-mass frame or boost parameter .
- 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 , scattering parameter , 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 ? (See discussion in Section 3.2.3 about this limit.)
These issues are presently rather well understood through NR simulations in spacetime dimensions but remain largely unanswered for the important cases .
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 . 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 . 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 of spacetime dimensions, (ii) the Lorentz factor or, equivalently, the collision velocity , and (iii) the impact parameter , where and 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 in dimensions and was analysed in Ref.  varying the boost parameter in the range . 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  zero-frequency limit (see Section 5.3)175*) yields 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.
Grazing collisions in four dimensions represent two-parameter studies, where the boost factor and the impact parameter 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 : (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 , the scattering threshold that we have already mentioned above and the threshold of immediate merger . This threshold has been identified in numerical BH simulations by Pretorius & Khurana  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 ; see also [413, 355]. This zoom-whirl-like behaviour has also been identified in high-energy grazing collisions in . The three different regimes are illustrated in Figure 13* which shows the BH trajectories for and , and . For this boost factor, the thresholds are given by and . For impact parameters close to the threshold values and , grazing collisions can generate enormous amounts of GWs. This is shown in the left panel of Figure 14*. Starting from in the head-on limit , the radiated energy increases by more than an order of magnitude to for . 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 . By fitting their numerical results, Shibata et al.  have found an empirical relation for the scattering threshold given by
Grazing collisions of spinning and non-spinning BHs have been compared in Ref. . The initial configurations for these simulations have been chosen with factors up to 2.49 and equal spin for both BHs of dimensionless magnitude and 0.65 aligned or anti-aligned with the orbital angular momentum . 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 at low collision speeds. At velocities above 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 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. ). 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 . 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 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 , the emitted GW energy decreases by about 3 orders of magnitude while the electromagnetic wave energy reaches a maximum at , and drops towards 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 .
An extended study of BH collisions using various analytic approximation techniques including geodesic calculations and the ZFL has been presented in Berti et al. ; see also [94*] for a first exploration in higher dimensions. Weak scattering of BHs in , which means large scattering parameters , and for velocities , has been studied by Damour et al.  using NR as well as PN and EOB calculations. Whereas PN calculations start deviating significantly from the NR results for , the NR calibrated EOB model yields good agreement in the scattering angle throughout the weak scattering regime.
BH collisions in 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.  for has been successfully generalized to higher in Ref.  and provides solutions with comparable accuracy as in .
A systematic exploration of the scattering threshold in dimensions has been performed in Ref. [587*]. By superposing non-rotating, boosted single BH initial data, they have evolved grazing collisions up to . 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 and the impact parameter . The simulations show a decrease of the scattering threshold at increasing velocity up to , similar to the case in the upper left panel of Figure 15*. At larger , 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.
The GW emission in BH mergers in 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 respectively (note that in 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 . 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 ) agree within numerical uncertainties . 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 .
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 . 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.
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 is the Ricci scalar associated to the metric , and and are arbitrary functions (see e.g., [92*, 824*] and references therein). The matter fields minimally coupled to are collectively denoted by . This form of the action corresponds to the choice of the so-called “Jordan frame”, where the matter fields obey the equivalence principle. For , 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 77*) enlarged to allow for a generic self-interaction potential (which could include the mass term). The label denotes quantities constructed from the Einstein-frame metric . In the Einstein frame the scalar field is minimally coupled to gravity, but any matter field is coupled to the metric .
In vacuum, this action leads to the following equations of motion:
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 , the equations of motion in the Einstein frame (180*) – (181*) admit GR (with ) 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 , at 2.5 PN order by Mirshekari and Will  and to all orders in the point particle limit in Ref. .
Thus, at least one of the following three ingredients are necessary to generate interesting dynamics in scalar-tensor theories:
- Nontrivial potential and initial conditions. Healy et al. studied an equal-mass BH binary in an inflation-inspired potential with nontrivial initial conditions on the scalar given by [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.
- 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*.
- Matter. When matter is present, new effects (due to the coupling of matter to the effective metric
) can dominate the dynamics and wave emission. For example, it has been shown that, for
, 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
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 ) 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. .
- 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 , 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.
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.
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 AdS5 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 AdS5 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 . 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 AdS5 have been presented in Ref. . 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 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. , 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  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 . 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].
Some initial applications of NR methods addressing specific issues in cosmology have been reviewed in the Living Reviews article by Anninos , 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. ). For further details, we refer the reader to Ref. .
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.  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.
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. ) 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 . The puncture locations in that work projected down to 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.  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. , 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 . For an earlier NR code developed to address inhomogeneous cosmologies see Ref. .
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 . These take into account interactions of dark matter, scalar perturbations, GWs, magnetic fields and turbulent plasma. Finally, Ref.  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.