Exact solutions of a nonlinear theory, such as GR, provide invaluable insights into the physical properties of the theory. Finding such solutions analytically and through a direct attack, that is by inserting an educated ansatz into the field equations, can be a tour de force, and, in general, only leads to success if a large isometry group is assumed from the beginning for the spacetime geometry. For instance, assuming spherical symmetry, in vacuum, leads to a fairly simple problem, whose general solution is the Schwarzschild metric . This simplicity is intimately connected with the inexistence of a spherically symmetric mode for gravitational radiation in Einstein gravity, which means that, in vacuum, a spherically symmetric solution must be static, as recognized by Birkhoff . On the other hand, assuming only axial symmetry leads to a considerably more difficult problem, even under the additional assumption of staticity. This problem was first considered by Weyl  who unveiled a curious and helpful mapping from these solutions to axially symmetric solutions of Newtonian gravity in an auxiliary -dimensional flat space; under this mapping, a solution to the latter problem yields a solution to the vacuum Einstein equations: a Weyl solution. For instance, the Schwarzschild solution of mass can be recovered as a Weyl solution from the Newtonian gravitational field of an infinitely thin rod of linear density and length . As we shall discuss in Section 4.1.2, the generalization of Weyl solutions plays an important role in the construction of qualitatively new solutions to the higher-dimensional Einstein equations.
Within the axially symmetric family of solutions, the most interesting case from the astrophysical viewpoint is the solution for a rotating source, which could describe the gravitational field exterior to a rotating star or the one of a rotating BH. An exact solution of Einstein’s equations describing the exterior of a rotating star has not been found (rotating stars are described using perturbative and numerical approaches [728*]),5 but in the case of a rotating BH, such a solution does exist. To obtain this stationary, rather than static, geometry, the Weyl approach by itself is unhelpful and new methods had to be developed. These new methods started with Petrov’s work on the classification of the Weyl tensor types . The Weyl tensor determines four null complex ‘eigenvectors’ at each point, and the spacetime is called ‘algebraically special’ if at least two of these coincide. Imposing the algebraically special condition has the potential to reduce the complicated nonlinear PDEs in two variables, obtained for a vacuum axially symmetric stationary metric, to ordinary differential equations. Using the (then) recently shown Goldberg–Sachs theorem , Kerr eventually succeeded in obtaining the celebrated Kerr metric in 1963 . This family of solutions was generalized to include charge by Newman et al. – the Kerr–Newman solution  – and to include a cosmological constant by Carter . In Boyer–Lindquist coordinates, the Kerr–Newman-(A)dS metric reads:
At the time of its discovery, the Kerr metric was presented as an example of a stationary, axisymmetric (BH) solution. The outstanding importance of the Kerr metric was only realized some time later with the establishment of the uniqueness theorems [188, 654]: the only asymptotically flat, stationary and axisymmetric, electrovacuum solution to the Einstein equations, which is non-singular on and outside an event horizon is the Kerr–Newman geometry. Moreover, Hawking’s rigidity theorem [406*] made the axisymmetric assumption unnecessary: a stationary BH must indeed be axisymmetric. Although the stability of the Kerr metric is not a closed subject, the bottom line is that it is widely believed that the final equilibrium state of the gravitational collapse of an enormous variety of different stars is described by the Kerr geometry, since the electric charge should be astrophysically negligible. If true, this is indeed a truly remarkable fact (see, however, Section 4.2 for “hairy” BHs).
Even if we are blessed to know precisely the metric that describes the final state of the gravitational collapse of massive stars or of the merger of two BHs, the geometry of the time-dependent stages of these processes seems desperately out of reach as an exact, analytic solution. To understand these processes we must then resort to approximate or numerical techniques.
3 there are various motivations to consider generalizations of (or alternative theories to) four-dimensional electrovacuum GR with . A natural task is then to address the exact solutions of such theories. Here we shall briefly address the exact solutions in two different classes of modifications of Einstein electrovacuum gravity: i) changing the dimension, ; ii) changing the equations of motion, either by changing the right-hand side – i.e., theories with different matter fields, including non-minimally coupled ones –, or by changing the left-hand side – i.e., higher curvature gravity. We shall focus on relevant solutions for the topic of this review article, referring to the specialized literature where appropriate.
- Changing the number of dimensions: GR in . Exact solutions in higher-dimensional GR,
, have been explored intensively for decades and an excellent review on the subject is
Ref. [308*]. In the following we shall focus on the vacuum case.
The first classical result is the generalization of the Schwarzschild BH, i.e., the vacuum, spherically – that is – symmetric solution to the -dimensional Einstein equations (with or without cosmological constant), obtained by Tangherlini [740*] in the same year the Kerr solution was found. Based on his solution, Tangherlini suggested an argument to justify the (apparent) dimensionality of spacetime. But apart from this insight, the solution is qualitatively similar to its four-dimensional counterpart: an analog of Birkhoff’s theorem holds and it is perturbatively stable.
On the other hand, the existence of extra dimensions accommodates a variety of extended objects with reduced spherical symmetry – that is – surrounded by an event horizon, generically dubbed as -branes, where stands for the spatial dimensionality of the object [441*, 285*]. Thus, a point-like BH is a 0-brane, a string is a 1-brane and so on. The charged counterparts of these objects have played a central role in SMT, especially when charged under a type of gauge field called ‘Ramond–Ramond’ fields, in which case they are called -branes or simply -branes . Here we wish to emphasize that the Gregory–Laflamme instability discussed in Section 3.2.4 was unveiled in the context of -branes, in particular black strings [367*, 368*]. The understanding of the nonlinear development of such instability is a key question requiring numerical techniques.
The second classical result was the generalization of the Kerr solution to higher dimensions, i.e., a vacuum, stationary, axially – that is6 – symmetric solution to the -dimensional Einstein equations, obtained in 1986 by Myers and Perry [565*] (and later generalized to include a cosmological constant [351*, 350*]). The derivation of this solution was quite a technical achievement, made possible by using a Kerr-Schild type ansatz. The solution exhibits a number of new qualitative features, in particular in what concerns its stability. It has independent angular momentum parameters, due to the nature of the rotation group in dimensions. If only one of these rotation parameters is non-vanishing, i.e., for the singly spinning Myers–Perry solution, in dimensions there is no bound on the angular momentum in terms of the BH mass . Ultra-spinning Myers–Perry BHs are then possible and their horizon appears highly deformed, becoming locally analogous to that of a -brane. This similarity suggests that ultra-spinning BHs should suffer from the Gregory–Laflamme instability. Entropic arguments also support the instability of these BHs [305*] (see Section 7.4 for recent developments).
The third classical result was the recent discovery of the black ring in [307*], a black object with a non-simply connected horizon, having spatial sections that are topologically . Its discovery raised questions about how the results on uniqueness and stability of vacuum solutions generalized to higher-dimensional gravity. Moreover, using the generalization to higher dimensions of Weyl solutions  and of the inverse scattering technique , geometries with a non-connected event horizon – i.e., multi-object solutions – which are asymptotically flat, regular on and outside an event horizon have been found, most notably the black Saturn . Such solutions rely on the existence of black objects with non-spherical topology; regular multi-object solutions with only Myers–Perry BHs do not seem to exist , just as regular multi-object solutions with only Kerr BHs in are inexistent [574, 424].
Let us briefly mention that BH solutions in lower dimensional GR have also been explored, albeit new ingredients are necessary for such solutions to exist. vacuum GR has no BH solutions, a fact related to the lack of physical dimensionality of the would be Schwarzschild radius , where is the 3-dimensional Newton’s constant. The necessary extra ingredient is a negative cosmological constant; considering it leads to the celebrated Bañados–Teitelboim–Zanelli (BTZ) BH [68*]. In a BH spacetime was obtained by Callan, Giddings, Harvey and Strominger (the CGHS BH), by considering GR non-minimally coupled to a scalar field theory [156*]. This solution provides a simple, tractable toy model for numerical investigations of dynamical properties; for instance see [55, 54] for a numerical study of the evaporation of these BHs.
- Changing the equations: Different matter fields and higher curvature gravity.
The uniqueness theorems of four-dimensional electrovacuum GR make clear that BHs are selective objects. Their equilibrium state only accommodates a specific gravitational field, as is clear, for instance, from its constrained multipolar structure. In enlarged frameworks where other matter fields are present, this selectiveness may still hold, and various “no-hair theorems” have been demonstrated in the literature, i.e., proofs that under a set of assumptions no stationary regular BH solutions exist, supporting (nontrivial) specific types of fields. A prototypical case is the set of no-hair theorems for asymptotically flat, static, spherically symmetric BHs with scalar fields [546*]. Note, however, that hairy BHs, do exist in various contexts, cf. Section 4.2.
The inexistence of an exact stationary BH solution, i.e., of an equilibrium state, supporting (say) a specific type of scalar field does not mean, however, that a scalar field could not exist long enough around a BH so that its effect becomes relevant for the observed dynamics. To analyse such possibilities dynamical studies must be performed, typically involving numerical techniques, both in linear and nonlinear analysis. A similar discussion applies equally to the study of scalar-tensor theories of gravity, where the scalar field may be regarded as part of the gravitational field, rather than a matter field. Technically, these two perspectives may be interachanged by considering, respectively, the Jordan or the Einstein frame. The emission of GWs in a binary system, for instance, may depend on the ‘halo’ of other fields surrounding the BH and therefore provide smoking guns for testing this class of alternative theories of gravity.
Finally, the change of the left-hand side of the Einstein equations may be achieved by considering higher curvature gravity, either motivated by ultraviolet corrections to GR, i.e., changing the theory at small distance scales, such as Gauss–Bonnet  (in ), Einstein-Dilaton-Gauss–Bonnet gravity and Dynamical Chern–Simons gravity [602*, 27*]; or infrared corrections, changing the theory at large distance scales, such as certain models. This leads, generically, to modifications of the exact solutions. For instance, the spherically symmetric solution to Gauss–Bonnet theory has been discussed in Ref.  and differs from, but asymptotes to, the Tangherlini solution. In specific cases, the higher curvature model may share some GR solutions. For instance, Chern–Simons gravity shares the Schwarzschild solution but not the Kerr solution [27*]. Dynamical processes in these theories are of interest but their numerical formulation, for fully nonlinear processes, may prove challenging or even, apart from special cases (see, e.g., [265*] for a study of critical collapse in Gauss–Bonnet theory), ill-defined.
- : The essential results in higher-dimensional vacuum gravity are the Tangherlini [740*] and Myers–Perry [565*] BHs, the (vacuum) black -branes [441*, 285*] and the Emparan–Reall black ring . Solutions with multi-objects can be obtained explicitly in with the inverse scattering technique. Their line element is typically quite involved and given in Weyl coordinates (see  for a list and references). The Myers–Perry geometry with a cosmological constant was obtained in in Ref.  and for general and cosmological constant in [351, 350]. Black rings have been generalized, as numerical solutions, to higher in Ref. [472*]. Black -branes have been discussed, for instance, in Ref. [441, 285]. In the best known examples of BH solutions are, respectively, the BTZ  and the CGHS BHs .
- Changing the equations of motion: Hawking showed  that in Brans–Dicke gravity the only stationary BH solutions are the same as in GR. This result was recently extended by Sotiriou and Faraoni to more general scalar-tensor theories . Such type of no-hair statements have also been proved for spherically symmetric solutions in GR (non-)minimally coupled to scalar fields  and to the electromagnetic field ; but they are not universal: for instance, a harmonic time dependence for a (complex) scalar field or a generic potential (together with gauge fields) are ways to circumvent these results (see Section 4.2 and e.g. the BH solutions in ). BHs with scalar hair have also been recently argued to exist in generalized scalar-tensor gravity .
Given the complexity of the Einstein equations, it is not surprising that, in many circumstances, stationary exact solutions cannot be found in closed analytic form. In this subsection we shall very briefly mention numerical solutions to such elliptic problems for cases relevant to this review.
The study of the Einstein equations coupled to nonlinear matter sources must often be done numerically, even if stationarity and spatial symmetries – typically spherical or axisymmetry – are imposed.7 The study of numerical solutions of elliptic problems also connects to research on soliton-like solutions in nonlinear field theories without gravity. Some of these solitons can be promoted to gravitating solitons when gravity is included. Skyrmions are one such case . In other cases, the nonlinear field theory does not have solitons but, when coupled to gravity, gravitating solitons arise. This is the case of the Bartnik–McKinnon particle-like solutions in Einstein–Yang–Mills theory . Moreover, for some of these gravitating solitons it is possible to include a BH at their centre giving rise to “hairy BHs”. For instance, in the case of Einstein–Yang–Mills theory, these have been named “colored BHs” . We refer the reader interested in such gravitating solitons connected to hairy BHs to the review by Bizoń  and to the paper by Ashtekar et al. .
A particularly interesting type of gravitating solitons are boson stars (see [685*, 516*] for reviews), which have been suggested as BH mimickers and dark matter candidates. These are solutions to Einstein’s gravity coupled to a complex massive scalar field, which may, or may not, have self-interactions. Boson stars are horizonless gravitating solitons kept in equilibrium by a balance between their self-generated gravity and the dispersion effect of the scalar field’s wave-like character. All known boson star solutions were obtained numerically; and both static and rotating configurations are known. The former ones have been used in numerical high energy collisions to model particles and test the hoop conjecture [216*] (see Section 7.3 and also Ref. [599*] for earlier boson star collisions and  for a detailed description of numerical studies of boson star binaries). The latter ones have been shown to connect to rotating BHs, both for Myers–Perry BHs in AdS [275*] and for Kerr BHs [422*], originating families of rotating BHs with scalar hair. Crucial to these connections is the phenomenon of superradiance (see Section 7.5), which also afflicts rotating boson stars . The BHs with scalar hair branch off from the Kerr or Myers–Perry-AdS BHs precisely at the threshold of the superradiant instability for a given scalar-field mode , and display new physical properties, e.g., new shapes of ergo-regions .
The situation we have just described, i.e., the branching off of a solution to Einstein’s field equations into a new family at the onset of a classical instability, is actually a recurrent situation. An earlier and paradigmatic example – occurring for the vacuum Einstein equations in higher dimensions – is the branching off of black strings at the onset of the Gregory–Laflamme instability [367*] (see Section 3.2.4 and Section 7.2) into a family of non-uniform black strings. The latter were found numerically by Wiseman  following a perturbative computation by Gubser . We refer the reader to Ref.  for more non-uniform string solutions, to Refs. [11, 790] for a discussion of the techniques to construct these numerical (vacuum) solutions and to  for a review of (related) Kaluza–Klein solutions. Also in higher dimensions, a number of other numerical solutions have been reported in recent years, most notably generalizations of the Emparan–Reall black ring [474, 472, 473] and BH solutions with higher curvature corrections (see, e.g., [131, 475, 132]). Finally, numerical rotating BHs with higher curvature corrections but in , within dilatonic Einstein–Gauss–Bonnet theory, were reported in .
In the context of holography (see Section 3.3.1 and Section 7.8), numerical solutions have been of paramount importance. Of particular interest to this review are the hairy AdS BHs that play a role in the AdS-Condensed matter duality, by describing the superconducting phase of holographic superconductors. These were first constructed (numerically) in . See also the reviews [396, 436] for further developments.
Finally, let us mention, as one application to mathematical physics of numerical stationary solutions, the computation of Ricci-flat metrics on Calabi–Yau manifolds .