Numerical Relativity and High Energy Physics: Recent Developments

We review recent progress in the application of numerical relativity techniques to astrophysics and high-energy physics. We focus on some developments that took place within the"Numerical Relativity and High Energy Physics"network, a Marie Curie IRSES action that we coordinated, namely: spin evolution in black hole binaries, high-energy black hole collisions, compact object solutions in scalar-tensor gravity, superradiant instabilities and hairy black hole solutions in Einstein's gravity coupled to fundamental fields, and the possibility to gain insight into these phenomena using analog gravity models.


Introduction
The last global meeting of the Numerical Relativity and High Energy Physics network -a Marie Curie International Research Staff Exchange Scheme (IRSES) partnership (2012)(2013)(2014)(2015) funded by the European Union and coordinated by the authors of this paper -started in Belém (Brazil) on 28 September 2015. Everyone in attendance, as well as the large majority of the scientific community, was unaware that a major breakthrough in science had just taken place: precisely two weeks earlier, the LIGO/Virgo collaboration observed the first gravitational-wave (GW) signal from the merger of two black holes (BHs). 1 The detection relied on decades of technological efforts to perform an apparently impossible measurement, corresponding to displacements that are ∼10 3 times smaller than the atomic nucleus. The unambiguous interpretation of the signal observed by Advanced LIGO as a BH binary coalescence was also the result of a decades-long effort: it took over 40 years to numerically solve Einstein' equations of general relativity (GR) and to understand the behavior of BH binaries through their inspiral, merger and ringdown.
A toolbox of powerful techniques became available after the numerical relativity breakthrough that took place in 2005. [2][3][4] This naturally led to a community effort looking for applications of these tools both in astrophysics 5 and beyond. 6,7 The area at the interface between numerical relativity and high-energy physics has made impressive strides in the recent past. A comprehensive review is beyond the scope of this paper, where we will focus on our own efforts in this emerging research field.
The plan of the paper is as follows. We will start in Sec. 2 by considering BH collisions, first in astrophysics (paying particular attention to some recent developments concerning spin dynamics), and then in the context of fundamental and high-energy physics (with applications to large extra dimension scenarios and to the gauge gravity duality). In Sec. 3, we will look at compact objectsi.e. BHs and neutron stars (NSs) -in alternative theories of gravity, focusing on models with scalar degrees of freedom in the gravitational sector: tensor-(multi)scalar theories, Horndeski gravity and Einstein-dilaton-Gauss-Bonnet gravity. In Sec. 4, we will consider GR minimally coupled to fundamental scalar and tensor fields and present some remarkable results obtained in the last few years in these simple models, including new types of numerical BH solutions that defied common lore. The existence of these BHs with scalar or Proca hair is intimately related with the complex phenomenon of superradiance, that can occur for rotating and charged BHs. Numerical relativity techniques have been (and will be) instrumental in probing the dynamics of these objects. Section 5 looks at many of these phenomena (in particular those discussed in Sec. 4) from a different perspective: that of analog gravity models. We close with some brief remarks.
several analytical and numerical studies. [22][23][24][25][26][27][28][29][30][31] Quasi-circular inspirals remain the most likely and best understood scenario, so here, we shall concentrate on this case. We will also focus on binaries in the framework of GR, but we note that there are preliminary explorations of scalar radiation from BH binaries in scalar-tensor theory, triggered either by nonasymptotically flat boundary conditions 32 or by a nonvanishing potential. 33 The estimation of source parameters in GW observations employs a method called matched filtering where the data stream is compared with a catalog of theoretically predicted GW templates 34 ; for an application of this technique to hybrid waveforms constructed out of numerical relativity and PN calculations see for example Ref. 35. A main challenge for the theoretical community is the generation of such template catalogs covering with high accuracy the whole range of BH binary masses and spins. Given the high computational cost of numerical relativity simulations, this construction typically stitches together post-Newtonian and numerical relativity waveforms, [36][37][38] or employs numerical simulations to calibrate free parameters in analytic prescriptions such as the effective-one-body model. [39][40][41][42] BH binaries with generic spins will undergo spin precession during which the orbital plane changes orientation. The modeling of these systems is significantly more involved than that of their nonprecessing counterparts, but benefits enormously from the presence of three distinctly different timescales. If we denote by r the separation of the two constituent BHs, these orbit around each other on the orbital timescale t orb ∝ r 3/2 , while the spin directions change on the precession timescale t pre ∝ r 5/2 , and the emission of GWs reduces the separation r on the radiation reaction timescale t RR ∝ r 4 . At sufficiently large separation r, this implies the hierarchy t orb t pre t RR . The first inequality has been used to derive orbit-averaged evolution equations for the individual spin vectors S i (i = 1, 2) fromṠ i = Ω i × S i , where the precession frequency depends on the orbital angular momentum L and the S i , but not on the separation vector r. [43][44][45] This quasi-adiabatic approach has been combined with some additional simplifications for the precession dynamics in order to construct template banks for precessing binaries. These techniques include a single effective spin model, modifications to the stationary-phase approximation, or the use of nonprecessing templates modulated through an effective precession parameter. [46][47][48][49][50][51] Orbit-averaged PN calculations have also been employed in the discovery of spin-orbit resonances 52 and for predictions of the final spins and recoil in BH binary mergers. [53][54][55] The success of the orbit-averaging procedure relies heavily on the analytic solutions for Keplerian orbits that are employed in the averaging over the orbital timescale. Until recently, no analogous analytic solution was known for the precession equations, so that the second inequality of the above hierarchy, t orb t RR , has not been brought to the same level of fruition. This picture changed with the identification of analytic solutions on the precessional timescale. 56,57 Consider for this purpose a BH binary with orbital angular momentum L, individual masses m i and spin vectors S i , i = 1, 2 and mass ratio q = m 2 /m 1 ≤ 1. For fixed mass ratio, the system is described by nine parameters, three each for S 1 , S 2 and L. Conservation of the spin magnitudes S i reduces this number to seven. On the precession timescale, the total angular momentum J ≡ S 1 + S 2 + L as well as the magnitude L are also conserved at the PN orders considered here, leaving three numbers to determine the state of the binary. A convenient choice for these variables is given by the angles θ i between the individual spins and the orbital angular momentum vector and the angle ∆Φ between the projections of the individual spins onto the orbital plane: cf. e.g. Fig. 1 in Ref. 57. One further variable can be eliminated through a convenient choice of a noninertial frame. Finally, the projected effective spin defined by 58,59 is conserved by the orbit-averaged spin-precession equations at 2PN order, and even under radiation reaction at 2.5PN order. Spin precession at this order is therefore described in terms of a single evolution variable, conveniently chosen to be the magnitude of the total spin S ≡ |S 1 + S 2 |. For a BH binary with specified parameters and separation, i.e. fixed values m i , S i , L, J, ξ, the precession is described completely in terms of the variable S. The set of physically allowed systems can then be represented as the area inside a closed loop constructed from two "effective potentials" ξ ± (S) in the (S, ξ) plane. The functions ξ ± (S) are determined by the physical constraints on the spin and angular momenta and are given in closed analytic form by Eq. (14) in Ref. 57. For a given value of ξ inside the range compatible with these constraints, this implies that S oscillates between two extrema S ± , where S min ≤ S − ≤ S + ≤ S max . a All remaining variables of the binary can be obtained from S through Eqs. (20) in Ref. 57 for θ 1 , θ 2 and ∆Φ which, in turn, determine all other physical variables.
A particularly intriguing consequence of this formulation of the spin precession dynamics is that all binaries fall into one of three morphologies, which are best characterized by the behavior of the angle ∆Φ on the precession time scale. As the variable S oscillates inside its allowed range, ∆Φ either (i) librates around 0, (ii) librates around π, or (iii) circulates through the entire range ∆Φ ∈ [−π, π].
a The resonance configurations of Ref. 52 correspond to the maximal and minimal allowed values of ξ in this area, ξmax and ξ min , at which S − = S + and, hence, S remains constant on the precession timescale; cf. Fig. 2 in Ref. 57. As the binary inspirals on the much larger radiation reaction timescale, the orbital and total angular momentum evolve, and the binary may undergo phase transitions between these morphologies. The inspiral of the binary under GW emission can be modeled in a remarkably efficient manner at 1PN order, if we express the binary separation r in terms of the orbital angular momentum L given by the Newtonian expression L = η(rM 3 ) 1/2 , where M ≡ m 1 + m 2 and η ≡ m 1 m 2 /(m 1 + m 2 ) is the symmetric mass ratio parameter. The evolution of the total angular momentum J averaged over a precession cycle is then given by 57 where The inspiral of precessing BH binaries is thus modeled in terms of a single ordinary differential equation (4) which, thanks to the precession-averaging procedure, can furthermore be solved numerically using much larger timesteps than possible in a formulation using only orbit-averaged variables. By suitably compactifying the variables involved, accurate numerical evolutions from infinite separations become possible at drastically reduced computational cost. The formalism truly bridges the gap between astrophysical BH separations and the regime close to merger, where numerical relativity predictions for BH kicks are valid. Further applications of the formalism identified a precessional instability of binaries where the spin of the more (less) massive BH is (anti) aligned with the orbital angular momentum, 60 and highlighted how the precessional morphology may carry a memory of the astrophysical processes that formed the binary. 57,61 Preliminary studies of the potential of present and future GW detectors to determine the morphologies in BH observations are encouraging, except for highly symmetric binaries, where precessional effects are suppressed. 62,63

Black hole collisions, fundamental and high-energy physics
Even 100 years after its publication, GR still confronts us with some of the most important questions in contemporary physics. As described in Sec. 4, astrophysical and cosmological observations suggest the presence of an enigmatic dark sector which appears to dominate the gravitational dynamics of much of our universe. At an even more fundamental level, GR predicts the limits of its own range of validity. Seminal work by Hawking and Penrose [64][65][66] demonstrated that gravitational collapse in the framework of GR leads to singularities under generic initial conditions. The appearance of infinities in the theory is expected to be cured by a future

1641022-6
Numerical relativity and high energy physics: Recent developments quantum theory of gravitation. Quite remarkably, however, relativity appears to have a built-in protection mechanism against the potentially fatal consequences of spacetime singularities; according to Penrose's cosmic censorship conjecture, [67][68][69] the singularities do not appear in naked form for physically realistic, generic initial data, but instead are cloaked inside horizons which causally disconnect the exterior spacetime from being influenced by the singularity. b A natural question then arising in the study of BH collisions is whether these can violate the conjecture as for example in the form of a super extremal Kerr BH generated in ultrarelativistic collisions. According to Thorne's hoop conjecture, 70 BH horizons should furthermore form (in D = 4 spacetime dimensions) whenever a physical system of mass M gets compacted inside a region with circumference 2πR s , where R s = 2M is the Schwarzschild radius associated with M . The conjecture has also been generalized to higher dimensions. 71 A particularly intriguing consequence arising from the hoop conjecture is the possibility of BH formation in proton-proton collisions at colliders such as the LHC, or in cosmic-ray showers hitting the Earth's atmosphere. In the trans-Planckian regime, where the colliding partons can be approximated as classical particles, the hoop conjecture predicts formation of a BH if the boost parameter γ satisfies γ c 4 R/(4Gm 0 ), i.e. if two particles of rest mass m 0 with center-of-mass energy M = 2γm 0 get compacted inside a volume of radius ∼ R. Taking the radius to be given by the de Broglie wavelength hc/M associated with the center-of-mass energy, the condition for BH formation becomes 72 (up to factors of order unity) M E P = c 5 /G, i.e. the center-of-mass energy of the collision must exceed the Planck energy. At the four-dimensional (4D) standard-model value E p ∼ 10 19 GeV, experimental tests of BH formation are clearly out of the range of present and foreseeable colliders. The so-called TeV gravity scenarios involving large or warped extra dimensions, 73-76 however, provide an appealing explanation of the hierarchy problem of physics and may lower the effective Planck energy to values as low as the TeV regime, which would allow for the possibility to form BHs in particle collisions at the LHC 77,78 ; for reviews see Refs. 79-81. Simulations of BH collisions can provide important information about the cross-section and energy loss through GWs which form key input for the Monte Carlo generators employed in the analysis of experimental data. 77,82 Yet another rich area of applying BH studies has emerged in the context of the gauge-gravity duality, often also referred to as the AdS/CFT correspondence, [83][84][85] which states the equivalence between string theory in asymptotically AdS spacetimes (times a compact space) and conformal field theories living on the AdS boundary. The duality provides a new approach to the (notoriously difficult) modeling of physical systems in strongly coupled gauge theories in terms of classical spacetimes, b This hypothesis is sometimes referred to as weak cosmic censorship and independent of the strong version that conjectures the inextendibility of the maximal cauchy development of generic compact or asymptotically flat initial data. often involving BHs, that are one dimension higher and asymptote to AdS at infinite radius.
We will next review some of the recent developments in BH modeling in the context of these topics, but we note that there are several more comprehensive reviews on these subjects. 6,7,11,86 The hoop conjecture has been tested in the context of high-speed collisions of scalar-field 72 and fluid-ball configurations. 87,88 In all these simulations, BH formation is observed at high velocities, consistent with the prediction of the hoop conjecture. In high-energy collisions, most of the center-of-mass energy is present in the form of the colliding object's kinetic energy. At sufficiently large velocities, this energy will be large enough to meet the conditions of the hoop conjecture and the colliding objects will appear as BHs. The hoop conjecture thus supports the idea that parton-parton collisions can be well modeled by colliding BHs, the GR analog of point particles. High-energy collisions have been studied most comprehensively in D = 4 spacetime dimensions and revealed a number of intriguing features. Numerical simulations of equal-mass, nonspinning BHs colliding head-on predict 89 that in the ultrarelativistic limit a fraction of 14 ± 3% (recently confirmed and refined to 13 ± 1% by the RIT group 90 ) of the total energy is radiated in GWs. This value is about half of Penrose's upper limit, 91,92 and in good agreement with the value of 16.4% obtained in second-order perturbative calculations on a background composed of two superposed Aichelburg-Sexl shock waves. [93][94][95][96][97][98][99][100][101] In grazing collisions, the BHs are allowed to approach each other with a nonzero impact parameter b, and the outcome of the collision depends on whether this parameter exceeds a threshold value or not. This scattering threshold b scat separates configurations that result in the formation of a single BH (b < b scat ) or in the constituents scattering off each other to infinity (b > b scat ), and was shown to be approximately given as a function of the collision speed v and the center-of-mass energy M c 2 by the remarkably simple formula 102 Numerical simulations 103 furthermore identified the presence of zoom-whirl orbits 104 in a regime where b is close to a critical value b * b scat . The energy released in GWs in these grazing collisions can be enormous, exceeding 35% of the total energy. Extrapolation to the speed of light, however, demonstrates that the maximum energy saturates at about half of the total (i.e. kinetic) energy. 105 The remaining kinetic energy, instead, ends up as rest mass, either in the single BH resulting from merger or in the two constituents in scattering configurations. Simulations of rotating BHs 105 also demonstrate that the impact of the spin on the collision dynamics is washed out in the limit v → c. We find here another confirmation of the "matter does not matter" conjecture already encountered in scalar-field and fluid-ball collisions 72,87,88 : Ultrarelativistic collisions are dominated by the kinetic energy, so that the internal structure of the colliding objects becomes irrelevant for the collision process. Further evidence for this conjecture has recently been obtained in numerical studies of unequal-mass BH collisions. Head-on collisions of this type emit ∼13% of the total mass in the form of GWs in the ultrarelativistic limit, in excellent agreement with the equal-mass result mentioned above. 106 It remains to be seen whether this picture remains intact as electric charge is added to the colliding particles. Collisions of electrically charged BHs have so far considered only the low-velocity regime, and revealed qualitatively similar dynamics as for the case of neutral BHs. 107 As intuitively expected, however, the collision is slowed down in the case of equal electric charges, reducing the energy radiated in electromagnetic and GWs as the charge-to-mass ratio approaches the critical value Q/M = 1. Conversely, the collision is accelerated by the additional attractive force between the BHs if they carry opposite charges, which increases the GW radiation by a factor of ∼2.7 as the charge-to-mass ratio |Q|/M increases from 0 to 0.99. 108 The electromagnetic radiation becomes dominant in these collisions at |Q|/M ∼ 0.37 and exceeds its GW counterpart by a factor ∼5.8 when |Q|/M = 0.99.
The high-energy collision of particles has also attracted considerable interest in a more astrophysical context through the collisional Penrose process. 109 Particle collisions near rapidly rotating BHs could in principle lead to arbitrarily large center-of-mass energies, 110 but there are several caveats on the astrophysical viability of this process. [111][112][113] The significance of such collisions is limited, in particular, by the redshift experienced by particles escaping from the near-horizon area to observers far away from the source. 114,115 An interesting possibility is that one of the colliding particles could have outgoing radial momentum: this can happen either because the particle reaches a turning point in the orbit, 116 or by allowing at least one outgoing particle to be generated close to the BH via previous collisions. 117 In both cases, the efficiency of the process (i.e. the ratio of the escaping particle's energy to the sum of the pre-collision particles' energies) can reach values as large as ∼13.9. 118,119 Collisions of BHs in D > 4 spacetime dimensions are not as well understood as the D = 4 case, mostly because of difficulties arising in the numerical stability of the simulations. The most extensive exploration of BH collisions in D = 5 was able to determine the scattering threshold in grazing collisions at velocities up to v ≈ 0.6 c. 120 This study used the so-called modified Cartoon method [121][122][123][124][125] and identified regions of exceptionally high curvature above the Planck regime that are not hidden inside a BH horizon. These regions with curvature radius below the Planck length are realized during the close encounter of two BHs in scattering configurations. The emission of GWs in D = 5 has been analyzed in head-on collisions of nonspinning BHs starting from rest. 126,127 It predicts E rad /M = (0.089 ± 0.006)% for equal masses and a monotonic decrease, in good agreement with point particle approximations, as the mass ratio is lowered from q = 1 to smaller values. The numerical method of this work is based on a dimensional reduction by isometry, 128 analogous to the Geroch decomposition. 129 Results from the two different codes have been compared in D = 5, demonstrating excellent agreement. 130 This work also provided the first estimate in D = 6, where equal-mass head-on collisions yield E rad /M = (0.081 ± 0.004)%. All these studies assume rotational symmetry in planes involving the extra dimensions such that the effective computational domain remains three-dimensional (3D). These symmetry assumptions still accommodate most of the scenarios relevant in the context of testing TeV gravity or fundamental properties of BH spacetimes.
Perturbative calculations based on superposed shock waves have also been extended to D ≥ 5 dimensions, including nonzero impact parameters. 92,131,132 These calculations predict a significant increase of the threshold impact parameter for formation of a common apparent horizon relative to the 4D case; see, in particular, Table 2 in Ref. 132. Extension of the work by d'Eath and Payne for the head-on case to D ≥ 5 resulted in a remarkably simple expression at first perturbative order for the energy fraction radiated in GWs 96,97 : E rad /M = 1/2 − 1/D. This result, originally obtained in a numerical study, has more recently been confirmed analytically to be exact at first-order. 100 The modeling of shock-wave and BH collisions is significantly more complicated in asymptotically AdS spacetimes. In contrast to the asymptotically flat case, the outer boundary now causally affects the interior. The use of large but finite computational domains, as often done in the modeling of asymptotically flat spacetimes, does not adequately capture this effect, so that incorporation of spacelike and/or null infinity, e.g. through compactification, becomes mandatory. The singular behavior of the metric components at the boundary furthermore requires careful numerical handling to avoid the generation of nonassigned numbers in the simulations. See, for example, Refs. 133 and 134 for a discussion of the numerical methods developed to handle these issues. Over the past decade, substantial progress has been made in overcoming these difficulties, achieving the first collision of BHs in an asymptotically AdS spacetime 135 (see also Ref. 136 for an earlier toy model). Numerical simulations of shock waves and BHs have clearly demonstrated their capacity for obtaining new insight into the strongly coupled regime of gauge theories. These studies have addressed in particular the thermalization of quarkgluon plasma in the heavy-ion collisions performed for example at the Brookhaven RHIC collider. Estimates for the thermalization time, i.e. the time after which the plasma reaches thermal equilibrium, obtained through the AdS/CFT correspondence are ∼0.35 fm/c, in good agreement with experimental data. [137][138][139][140] As in the asymptotically flat case, point-particle and perturbative calculations provide a convenient tool complementing full numerical relativity studies. 141,142 For a more extended discussion of BH and shock wave collisions in AdS, we refer the reader to Ref. 7.
Finally, we remark that BH collisions in asymptotically de Sitter spacetimes have been considered as tests of the cosmic censorship conjecture. The numerical simulations support the conjecture. 143

Compact Objects in Modified Theories of Gravity
In this section, we will discuss isolated compact objects (BHs and NSs) in modified theories of gravity. We will focus on one of the most natural and best studied extensions of GR: scalar-tensor gravity, in which one or more scalar degrees of freedom are included in the gravitational sector through a nonminimal coupling. [144][145][146][147][148][149] We will discuss compact objects in these theories at increasing levels of complexity, starting from the "standard" Bergmann-Wagoner formulation (Sec. 3.1) and then considering extensions to multiple scalar fields (Sec. 3.2), Horndeski gravity (Sec. 3.3) and Einstein-dilaton-Gauss-Bonnet gravity (Sec. 3.4).

"Bergmann-Wagoner " scalar-tensor theories
The most general action of scalar-tensor gravity, at most quadratic in derivatives of the fields and with one scalar field, was studied by Bergmann and Wagoner. 150,151 The action of this theory can be written (with an appropriate field redefinition) as: where U (φ) and ω(φ) are arbitrary functions of the scalar field φ, and S M is the action of the matter fields Ψ. When U (φ) = 0 and ω(φ) = ω BD is constant, the theory reduces to (Jordan-Fierz-)Brans-Dicke gravity. [152][153][154] The Bergmann-Wagoner theory (8) can be expressed in a different form through a scalar field redefinition ϕ = ϕ(φ) and a conformal transformation of the metric g µν → g µν = A −2 (ϕ)g µν . In particular, fixing A(ϕ) = φ −1/2 , the Jordan-frame action (8) transforms into the Einstein-frame action where g and R are the determinant and Ricci scalar of g µν , respectively, and the potential V (ϕ) ≡ A 4 (ϕ)U (φ(ϕ)). The price paid for the minimal coupling of the scalar field in the gravitational sector is the nontrivial coupling in the matter sector of the action: particle masses and fundamental constants depend on the scalar field. The actions (8) and (9) are just different representations of the same theory, 155,156 so it is legitimate (and customary) to choose the conformal frame in which calculations are simpler. For instance, in vacuum the Einstein-frame action (9) formally reduces to the GR action minimally coupled with a scalar field. It may then be necessary to change the conformal frame when extracting physically meaningful statements (since the scalar field is minimally coupled to matter in the Jordan frame, test particles follow geodesics of the Jordan-frame metric, not of the Einsteinframe metric). The relation between Jordan-frame and Einstein-frame quantities is 157 Many phenomenological studies neglect the scalar potential, setting U (φ) = 0 or V (ϕ) = 0. c If the potential vanishes, the theory is characterized by a single function α(ϕ). The expansion of this function around the asymptotic value ϕ 0 can be written in the form The choice α(ϕ) = α 0 = const. (i.e. ω(φ) = const.) corresponds to Brans-Dicke theory. A more general formulation, proposed by Damour and Esposito-Farèse, is parametrized by α 0 and β 0 . 159,160 Another simple variant is massive Brans-Dicke theory, in which α(ϕ) is constant, but the potential is nonvanishing and has the form so that the scalar field has a mass m 2 s ∼ U (φ 0 ). Since, the scalar field ϕ in the action (9) is dimensionless, the function α(ϕ) and the constants α 0 , β 0 are dimensionless as well.
In the Einstein-frame, the field equations are where is the Einstein-frame stress-energy tensor of matter fields, and T = g µν T µν is its trace. Equation (11b) shows that α(ϕ) determines the strength of the coupling of the scalar fields to matter. 144,161 Astrophysical observations set bounds on the parameter space of scalar-tensor theories. In the case of Brans-Dicke theory, the best observational bound (α 0 < 3.5 × 10 −3 ) comes from the Cassini measurement of the Shapiro time delay. 162 An interesting feature of scalar-tensor gravity is the prediction of characteristic physical phenomena which do not occur in GR. Even though, we know from observations that α 0 1 and that GR deviations are generally small, these phenomena may lead to observable consequences. The best known example is the fact that compact binary systems in scalar-tensor gravity emit dipolar gravitational radiation. 163,164 Dipolar gravitational radiation is "pre-Newtonian", i.e. it occurs at lower PN order than quadrupole radiation, and it does not exist in GR. In the more general case with β 0 = 0, the phenomenon of spontaneous scalarization (described below) can lead in principle to macroscopic modifications in the structure of NSs, significantly c This approximation corresponds to neglecting the cosmological term, the mass of the scalar field and possible self-interactions. In an asymptotically flat spacetime, the scalar field tends to a constant φ 0 at spatial infinity, corresponding to a minimum of U (φ). Taylor expanding U (φ) around φ 0 yields, at the lowest orders, a cosmological constant and a mass term for the scalar field. 151,158 affecting the amount of dipolar radiation emitted by a binary system. Therefore the best constraints in the (α 0 , β 0 ) plane come from observations of NS-NS and NS-WD binary systems. 165 Observations of compact binary systems also constrain massive Brans-Dicke theory, leading to exclusion regions in the (α 0 , m s ) plane. 158

Spontaneous scalarization in compact stars
An interesting feature of scalar-tensor theories is the existence of nonperturbative NS solutions in which the scalar field amplitude is finite even for α 0 1: this phenomenon, known as spontaneous scalarization, 159,160 may significantly affect the mass and radius of a NS, and therefore the orbital motion of a compact binary system, even far from coalescence. A simple way to illustrate the principle behind spontaneous scalarization is by taking the limit in which the scalar field ϕ is a small perturbation around a GR solution. 166 Expanding around the constant value ϕ 0 to first-order inφ ≡ ϕ − ϕ 0 1, the field equations in the Einstein-frame (11) read Here, we have assumed analyticity around ϕ ∼ ϕ 0 and we have used Eq. (10). It is clear from Eq. (13b) that α 0 controls the effective coupling between the scalar and matter. Various observations, such as weak-gravity constraints and tests of violations of the strong equivalence principle, require α 0 to be negligibly small when the scalar tends to its asymptotic value. 160,165,167 This implies that a configuration in which the scalar ϕ ≈ ϕ 0 and α 0 ≈ 0 should be at least an approximate solution in most viable scalar-tensor theories. With α 0 = 0, any background GR solution solves the field equations above at first-order in the scalar field. At this order, the Klein-Gordon equation reads Thus, the coupling of the scalar field to matter is equivalent to an effective x νdependent mass. Depending on the sign of β 0 T , the effective mass squared can be negative. Because typically d −T ≈ ρ > 0 , this happens when β 0 < 0. When µ 2 s < 0 in a sufficiently large region inside the NS, scalar perturbations of a GR equilibrium solution develop a tachyonic instability associated with an exponentially growing mode, which causes the growth of scalar hair in a process similar to ferromagnetism. 159,160 Spherically symmetric NSs develop spontaneous scalarization for β 0 −4.35. 170 Detailed investigations of stellar structure, 160,171 numerical simulations of collapse [172][173][174][175] and stability studies 170,176 confirmed that spontaneously scalarized configurations would indeed be the end-state of stellar collapse in these theories. d Some nuclear equations of state (EOSs) allow for positive T in the NS interior, with potentially interesting phenomenological implications. 168,169 This is subject, however, to the collapse process reaching a sufficient level of compactness. Recent simulations of supernova core collapse identified clear signatures of spontaneous scalarization when the collapse ultimately formed a BH, but not in case of neutron-star end products, as these were not compact enough. Further studies using more elaborate treatment of the microphysics and/or relaxing symmetry assumptions are needed to determine how generic a feature this is in core collapse scenarios. Finally, spontaneously scalarized configurations may also be the result of semiclassical vacuum instabilities. [177][178][179][180] The nonradial oscillation modes of spontaneously scalarized, nonrotating stars were studied by various authors. [181][182][183][184] The bottom line is that the oscillation frequencies can differ significantly from their GR counterparts if spontaneous scalarization modifies the equilibrium properties of the star (e.g. the mass-radius relation) by appreciable amounts. However, current binary pulsar observations yield very tight constraints on spontaneous scalarization -implying in particular that β −4.5 -and the oscillation modes of scalarized stars for viable theory parameters are unlikely to differ from their GR counterparts by any measurable amount. Note, however, that the binary pulsar constraints on β apply to the case of massless ST theories. For massive scalars, much larger (negative) β and correspondingly stronger effects on the structure and dynamics of compact objects may be possible. 185 Spinning NSs at first-order in the Hartle-Thorne slow-rotation approximation were studied by Damour and Esposito-Farèse 160 and later by Sotani. 187 At firstorder in rotation, the scalar field only affects the moment of inertia, mass and radius of the NS. Second-order calculations 186 are necessary to compute corrections to the spin-induced quadrupole moment, tidal and rotational Love numbers, as well as higher order corrections to the NS mass and to the scalar charge. Figure 1 shows representative examples of the properties of NSs in a scalar-tensor theory with spontaneous scalarization at second-order in the rotation parameter.
Rapidly rotating NSs in scalar-tensor theories were recently constructed 188 by extending the RNS code. 189 Scalarization effects are stronger -and deviations from GR are larger -for rapidly spinning NSs. 190,191 Therefore, despite the tight binary pulsar bounds, it is still possible that spontaneous scalarization may occur in rapidly rotating stars.
Old, isolated NSs, as well as the NSs whose inspiral and merger, we expect to observe with GW detectors, are expected to be rotating well below their massshedding limit. However, these considerations may not apply just before merger, where the rotational frequencies of each NS may approach the mass-shedding limit. In these conditions, numerical simulations have also recently revealed the possibility of "dynamical scalarization" -a growth of the scalar field that may significantly affect the waveform near merger, and potentially be detectable. [192][193][194][195][196] A more exotic mechanism to amplify the effects of scalarization is anisotropy in the matter composing the star. 197 Nuclear matter may be anisotropic at very high densities, where the nuclear interactions must be treated relativistically and phase transitions (e.g. to pion condensates or to a superfluid state) may occur. For example, Nelmes and Piette 198 recently considered NS structure within the Skyrme model -a low energy, effective field theory for quantum chromodynamics (QCD) -finding significant anisotropic strains for stars with mass M 1.5M (see also work by Adam et al. 199,200 ). The effect of anisotropy is shown in Fig. 2. For illustration, in the figure we adopt a very simple model developed in the 70's by Bowers and Liang,201 where the degree of anisotropy is parametrized by a parameter λ BL . The left panel shows the critical threshold for scalarization as a function of stellar compactness for several "ordinary" (isotropic) EOSs: the EOS has almost no effect on the critical threshold for scalarization, which is always around β = −4. 35  that the critical β for scalarization (and, as it turns out, also the effects of scalarization on macroscopic NS properties) increases (decreases) when the tangential pressure is bigger (smaller) than the radial pressure. An interesting feature of the Bowers-Liang models is that it allows for stellar configurations with compactness approaching the Schwarzschild limit r = 2M . Yagi and Yunes used this observation to study the recently found "I-Love-Q" universal relations -which relate bulk NS properties such as the moment of inertia, spininduced quadrupole moment and tidal deformability in an EOS-independent wayas NSs approach the BH limit. 202-204

Black hole hair ?
The phenomenology of scalar-tensor theory in vacuum spacetimes, such as BH spacetimes, is less interesting. When the matter action S M can be neglected, the Einstein-frame formulation of the theory is equivalent to GR minimally coupled to a scalar field. BHs in Bergmann-Wagoner theories satisfy the same no-hair theorem as in GR, and thus the stationary BH solutions in the two theories coincide. [205][206][207] Moreover, dynamical (vacuum) BH spacetimes satisfy a similar generalized no-hair theorem: the dynamics of a BH binary system in Bergmann-Wagoner theory with vanishing potential are the same as in GR, 144 up to at least 2.5 PN order for generic mass ratios 208 and at any PN order in the extreme mass-ratio limit. 166 If there is more than one massive real scalar field, however, or at least one massive complex scalar field, the situation concerning stationary BH solutions can be very different: axisymmetric, hairy BHs do exist , 148,209,210 as will be reviewed in Sec. 4. Tensor-multi-scalar theories have indeed received more attention in the recent literature, as we now discuss.

Tensor-multi-scalar theories
A natural generalization of the Bergmann-Wagoner formulation (8) consists in including more than one scalar field coupled with gravity. The action of tensormulti-scalar (TMS) gravity 144,211 is: where F and V are functions of the N scalar fields φ a (a = 1 . . . N). The presence of multiple fields allows for a wider range of kinetic terms, i.e. terms quadratic in the fields' first derivatives. These are conveniently formulated in terms of a manifold (separate from the spacetime manifold) spanned by the scalar fields themselves. The scalar fields live on this manifold (the target space) with metric γ ab (φ). The action (15) is invariant not only under spacetime diffeomorphisms, but also under target-space diffeomorphisms, i.e. scalar field redefinitions. TMS theories are more complex than theories with a single scalar field, since the geometry of the target space can affect the dynamics.
The simplest extension of a ST theory with a single real scalar field is a theory with two real scalar fields. If we work, equivalently, with a single complex scalar ϕ instead, the action reduces to Hereafter, we assume that the potential vanishes, i.e. V (ϕ,φ) = 0, and that the target space is maximally symmetric. Upon stereographic projection and field redefinition 211 the target-space metric can be written as where Ö is the radius of curvature of the target-space geometry. For a spherical geometry we have Ö 2 > 0, for a hyperbolic geometry Ö 2 < 0, and in the limit Ö → ∞ the geometry is flat. The function A(ϕ,φ) determines the scalar-matter coupling. What enters the field equations is actually the function κ, defined as Without loss of generality, we assume that far away from the source the field vanishes: ϕ ∞ = 0. We then expand the function log A in a series about ϕ = 0: where β 0 is real, while α * and β * 1 are in general complex numbers. Redefine β * 1 ≡ β 1 e iθ , where θ is chosen such that β 1 is real. Then, after defining α * ≡ αe iθ/2 and 1641022-17 E. Berti et al. a new field ψ ≡ ϕe iθ/2 , the field equations become where log A(ψ,ψ) = αψ +ᾱψ + 1 2 and in the second line, we have split the field ψ into real and imaginary parts: . The structure of this theory when α = 0 is determined by three real parameters: β 0 + β 1 , β 0 − β 1 and the target-space curvature Ö 2 . When α = 0, two further parameters (|α| and arg α) are necessary to define the theory. This two-scalar model is the simplest generalization of the spontaneous scalarization model by Damour and Esposito-Farèse. 159 Note that the quantity |α| 2 ≡ αᾱ ≡ Re[α] 2 + Im[α] 2 is strongly constrained by observations, similarly to the single-scalar case. However, in TMS theories α is a complex quantity and its argument, arg α, is unconstrained in the weak-field regime. When α = 0, the conformal coupling at second-order in ψ reduces to log A(ψ,ψ) = 1 2 β 0 ψψ + 1 4 Compact stars in theories with α = 0 and α = 0 are rather different. When α = 0 and β 1 = 0, the theory is invariant under the symmetries ψ →ψ and ψ → −ψ. In this case, we only found solutions where either the real or the imaginary part of the scalar field has a nontrivial profile; these theories are effectively equivalent to single-scalar theories.
When α = 0 the situation is more interesting, as shown in Fig. 3. Introduction of α ∈ R partially breaks this symmetry down to conjugation only, whereas introduction of α ∈ C fully breaks the symmetry. Now GR configurations are not solutions of the field equations. In particular, a constant (or zero) scalar field does not satisfy Eq. (21) when T = 0. Therefore, it is not surprising that when α = 0 we can find solutions with two nontrivial scalar profiles even when β 0 = β 1 = 0. A more interesting question is whether there are stellar configurations in which both scalar fields have a large amplitude. These "biscalarized" solutions are absent in the α = 0 case, but as it turns out they exist when α = 0. For concreteness, in the figure we set |α| = 10 −3 , so that the theory is compatible with experimental bounds from binary pulsar observations. We set 1/Ö = 0 (i.e. for simplicity we consider a flat target space), we fix β 0 = −5, and we vary arg α in the range [0, 2π] in steps of π/6. As shown in Fig. 3   , however, remains satisfied, so that scalarized models should cluster close to the Re[ψ 0 ]-axis. This is indeed observed in the bottom panels of Fig. 3. A more detailed investigation of the phenomenology of these models is underway.

Horndeski theories
Besides the obvious addition of one or more scalar field(s), a second possibility to generalize scalar-tensor theories of the Bergmann-Wagoner type has recently attracted a great deal of attention. The theory in question was first formulated by Horndeski,212 and it is the most general single-scalar theory with second-order field equations. In "modern" notation, the action of Horndeski gravity can be written in terms of Galileon interactions 213 as where The functions G i = G i (φ, X) depend only on the scalar field φ and its kinetic energy, X = −∂ µ φ∂ µ φ/2. For brevity, we have also defined the shorthand notation . An attractive feature of Horndeski gravity is its generality. The theory includes a broad spectrum of phenomenological dark energy models, as well as modified gravity theories with a single scalar degree of freedom. Some important special limits of the theory are listed below: (1) GR corresponds to choosing G 4 = 1/2 and G 2 = G 3 = G 5 = 0.
(2) When G 4 = F (φ) and all other G i 's are zero, we recover a scalar-tensor theory with nonminimal coupling of the form F (φ)R. Therefore Brans-Dicke theory and f (R) gravity are special cases of Horndeski gravity. (3) A theory that we will consider in some detail below, namely Einstein-dilaton-Gauss-Bonnet (EdGB) gravity, has the action where V (φ) is the scalar potential, ξ(φ) is a coupling function and is the Gauss-Bonnet invariant. This theory can be recovered with the choices 214,215

1641022-20
Numerical relativity and high energy physics: Recent developments where we have defined ξ (n) ≡ ∂ n ξ/∂φ n . 214 (4) A theory with nonminimal derivative coupling between the scalar field φ and the Einstein tensor G µν (the "John" Lagrangian in the language of the so-called "Fab Four" model 216,217 ), with action can be constructed by setting where Λ 0 , η, ζ and β are constants. Incidentally, a coupling of the form G µν φ µ φ ν can also be obtained by setting G 5 = −φ and integrating by parts. 218 (5) The Lagrangian L 2 corresponds to the k-essence field. [219][220][221] For this reason, some of the literature denotes the function G 2 by the letter K. Because of the generality of Horndeski gravity, a comprehensive review of compact objects would inevitably have to discuss important subclasses that have been studied for a long time. 5 A more specific review of compact objects in the subclasses (4)-(6) can be found in this same volume, 223 and some examples are also discussed in another review. 148 In the next section, we complement these reviews focusing on recent work in EdGB gravity.

Einstein-dilaton-Gauss-Bonnet gravity
In EdGB gravity 224 (see Sec. 3 (27) is coupled with a scalar field. e The resulting action of EdGB gravity, Eq. (26), is a special case of Horndeski gravity, as discussed in Sec. 3.3. With the choice ξ(φ) = αφ this theory is shift-symmetric, and it has been shown 225 that it is the only shift-symmetric Horndeski gravity theory in which the no-hair theorems do not hold.

.3), the Gauss-Bonnet invariant
EdGB gravity can also be seen as belonging to a different class of modified gravity: that of quadratic gravity theories, 226,227 in which quadratic curvature terms are included in the action. EdGB gravity holds a special place as the only quadratic gravity theory with equations of motion of second differential order. Other theories of quadratic curvature gravity (e.g. dynamical Chern-Simons gravity 228,229 ) have equations of motion of higher differential order, and are then subject to Ostrogradski's instability. 230 In order to avoid this instability, they should be considered as e The normalization of the scalar is different from those in Secs. 3.1 and 3.2, by a factor 2. effective theories, obtained as truncations of more general theories. In other words, EdGB gravity is consistent for any value of the coupling constant, while other quadratic gravity theories should only be considered in the weak-coupling limit. Note also that the EdGB term without the coupling to a scalar field would be trivial, since R 2 GB is a total derivative. Including a quadratic curvature term in the action is an interesting modification of GR, for a variety of reasons. First of all, this is the simplest way to modify the strong-curvature regime of gravity, and, second, it is also a way to circumvent no-hair theorems (see the discussion in Secs. 3.1 and 4 for different ways to grow BH hair). Moreover, quadratic curvature terms can make the theory renormalizable. 231 In particular, the EdGB term naturally arises in low-energy effective string theories 232,233 when ξ(φ) = αe φ /4.
Hereafter, we consider EdGB gravity with ξ(φ) = αe φ /4. The first BH solution of this theory, derived about 20 years ago by Kanti et al., 224 is a numerical solution describing a spherically symmetric BH. The solution has scalar hair, i.e. a nontrivial configuration of the scalar field, but only secondary hair (the scalar charge D is determined by the mass M , and hence is not a free parameter. It can be shown that Kanti's solution only exists for 224,234 where M is the BH mass. The best observational bound on the coupling parameter is α 47 M 2 . 235 This bound is weaker than the theoretical bound (31) for BHs with M 8.2 M . 236 In recent years, numerical solutions describing slowly rotating 234 and rapidly rotating 237 BHs have been derived. These solutions describe stationary BHs for all values of the mass and the spin, and for all values of the coupling parameter in the allowed range (31). However, these solutions require a numerical integration for each set of parameters. In order to devise and implement observational tests based on astrophysical or GW observations (for instance, for Monte Carlo data analysis), an approximate, analytical solution can be more useful than numerical solutions.
Analytical BH solutions in EdGB theory have been determined as perturbative expansions in the dimensionless coupling parameter ζ = α/M 2 and the dimen- 239 and finally at order O(ζ 7 ,ā 5 ). 240 For a slowly rotating BH, the solution derived in Ref. 240 reproduces the most relevant geodesic quantities (the ISCO location and the epicyclic frequencies) within 1%, for the entire allowed range (31) of the coupling parameter.
Astrophysical observations from the near-horizon region of BHs can allow tests of GR against modified theories (such as EdGB gravity) which predict deviations in the strong-field, high-curvature regime. Indeed, near the horizon of stellar-mass BHs the spacetime curvature is very large, and (for sufficiently large values of α) BH solutions in EdGB theory may be significantly different from the Kerr solution. These deviations can affect observable quantities, such as the quasi-periodic

1641022-22
Numerical relativity and high energy physics: Recent developments oscillations (QPOs) observed in the X-ray flux of accreting BHs. 236 Indeed, in many astrophysical models the frequencies of these QPOs are appropriate combinations of the epicyclic frequencies of the (near-horizon) BH geodesics, in which the strongfield regime of gravity is manifest. Therefore, future large-area X-ray telescopes such as LOFT 241 could set constraints on the coupling parameter α. For instance, the detection of two QPO triplets from a BH with M = 5.3 M andā = 0.2 by a detector having the LOFT design sensitivity could exclude the range ζ 0.4 with 3σ confidence. 236

Implications of Superradiant Instabilities for Fundamental Physics and Astrophysics
In the previous section, we focused on specific modifications of Einstein's gravity and on the different physical consequences, as well as compact object solutions, that arise in these models. Perhaps somewhat more surprisingly, even within Einstein's gravity, considering simple fundamental fields that satisfy the energy conditions can also lead to new types of compact objects, with interesting physical consequences.
In this section, we will review these recent developments, that are related to the complex phenomenon known as superradiance. 242 Other relevant recent developments related to superradiance, in different physical and mathematical directions to the ones reviewed here, can be found, for instance, in Refs. 243-248.

Setup
Einstein's GR minimally coupled to fundamental fields, such as massive scalars or vectors, is described by the Lagrangian We have defined κ = 16π and F µν ≡ ∇ µ A ν − ∇ ν A µ is the Maxwell tensor. Both the scalar and vector fields are assumed to be complex, for reasons that will become clear soon. The mass m B of the bosons under consideration is related to the mass parameters above through µ S,V = m B / . By "fundamental" we mean fields which are not effective descriptions of other microscopic degrees of freedom. For most of the analysis below, however, the true nature of these fields (i.e. whether they are truly fundamental or rather a coarse-grained representation of more fundamental degrees of freedom), is not relevant. Each of them is completely equivalent to two real scalar or vector fields, but some of our considerations below apply equally well to one or many real scalar and vector fields. The theories represented by this action are relevant for several reasons. Because they are simple, they can be thought of as proxies for more complex interactions, of which they would be faithful models in certain regimes (presumably when higher order interactions are negligible). Fundamental bosons also play a key role in particle physics. For instance they could describe the axion or axionlike particles, originally 1641022-23 E. Berti et al. intended to solve the strong-CP problem in QCD, which recently gained prominence as dark-matter candidates. [249][250][251] In this context, self-gravitating solutions of fundamental fields allow us to understand and study quantitatively the growth of dark matter structures and their clustering inside compact stars. 252,253 Whether or not they form a significant component of dark matter, minimally coupled fundamental fields should obey the equivalence principle and freely fall in the same way as standard model fields. Thus, the most promising channel to look for their imprints consists of gravitational interactions.

Superradiance and superradiant instabilities
Fundamental fields in the presence of gravity display of course a panoply of interesting effects, such as the critical phenomena identified in Choptuik's seminal study. 254 In strong gravitational fields, one of the most peculiar is superradiance, i.e. the amplification of low-frequency waves scattering off rotating BHs. 242,255,256 Superradiance is required by the second law of thermodynamics, and is akin to tidal acceleration in planetary dynamics. 257 Superradiance is active for low-frequency, bosonic fields satisfying the superradiance condition ω < mΩ (33) with m an integer azimuthal number and Ω the angular frequency of the BH. The amplitude of the superradiant amplification of any incident wave depends on the rotation Ω, on the wave frequency ω and on the field being scattered. 242,258 Superradiant mechanisms can trigger instabilities in spacetimes that are able to confine the fluctuations. In such cases, the wave is forced to bounce back and forth, being repeatedly amplified upon interaction with the BH, and leading to exponential growth of linearized fluctuations. This mechanism is called a black hole bomb, [259][260][261][262] and leads to instabilities in truly confined spacetimes like anti-de Sitter. [263][264][265][266][267] It is interesting that the same mechanism also makes Kerr BHs unstable under massive, scalar-field fluctuations, 268-271 vector-field fluctuations [272][273][274] or even tensor-field perturbations. 275 Physically, massive states prevent full leakage to infinity and act as an effective barrier for low-frequency waves.

Hairy black holes bifurcating from the Kerr solution
Since Kerr BHs are unstable against sufficiently low frequency modes of a massive bosonic field, a relevant question is: what is the endpoint of the instability? While this is still an open question (but see Sec. 4.3.4), a relevant observation is the existence of stationary, asymptotically flat BH solutions of the model (32), which are regular on and outside the event horizon and for which the horizon is in equilibrium with a nontrivial scalar or vector field condensate. Moreover these BHs are continuously connected with the Kerr solution, and as such they have been dubbed Kerr BHs with scalar 209,210,276 or Proca hair. 277 They are manifestly related to the phenomenon of superradiance, as they exist at the threshold of the inequality (33), and

1641022-24
Numerical relativity and high energy physics: Recent developments they are likely to play a role either as endpoints or as long-lived intermediate states in the development of the superradiant instability of Kerr BHs in the presence of massive scalar or vector fields.
The existence of these hairy BH solutions raises three immediate questions: (1) "How is it possible that stationary, asymptotically flat BH solutions different from Kerr exist in the very simple model (32), in view of the well-known no-hair theorems that apply to this model (in particular the pioneering theorems due to Bekenstein for the scalar 278 and Proca 279,280 cases)?" (see also Ref. 148 for a review of no-hair theorems applying to the scalar case). (2) "If these hairy BH solutions are continuously connected to the Kerr solution, then there must be an imprint of their existence when we consider the corresponding matter fields on the Kerr background as test fields. Is it so?" (3) "Do these BHs trivialize in the limit of vanishing horizon or is there some residual gravitating configuration in this limit?" We shall tackle each of these questions in the following three subsections.

Circumventing no-hair theorems
The answer to (1) is simple and enlightening: theorems have assumptions and assumptions can be dropped. In the present case, an assumption in many of the no-hair theorems, including those of Bekenstein, is that the metric and the matter field share the same symmetries. This is not necessary: the spacetime and the energy-momentum tensor should share the same symmetries, but not the matter field itself. This apparently innocuous observation allows us to circumvent the simplest Bekenstein-type no-hair theorems, but observe that it is a necessary but not sufficient ingredient. The reason will become clear in the following. The metric ansatz that has been successfully used for finding (nonextremal) Kerr BHs with scalar 209 and Proca 277 hair reads: where N (r) ≡ 1 − r H /r and r H is a constant. The metric is completely determined by four functions of the spheroidal coordinates (r, θ). These coordinates reduce to prolate spheroidal coordinates (rather than the more familiar oblate spheroidal ones, obtained in the flat-spacetime limit of the Boyer-Lindquist form of the Kerr metric) in an appropriate Minkowski limit. 277  is the event horizon. Observe that the line element (34) admits two independent Killing vector fields: k ≡ ∂ t and m = ∂ ϕ .
On the other hand, the "matter" ansatz that has been used to find the hairy BHs is for the scalar case, 209 and for the Proca case. 277 The two constant parameters ω, m are the frequency and azimuthal quantum number, with ω ∈ R + , m ∈ Z/{0}. An immediate observation is that the matter fields are not invariant under the two aforementioned Killing vector fields: but the corresponding energy-momentum-tensors are Thus Bekenstein-type theorems are inapplicable and the absence of hair is no longer guaranteed, but the left-and right-hand sides of the Einstein equations still have the same symmetries. The ansatz (34), in combination with (35) or (36), yields axially symmetric solutions. One may wonder whether BH solutions could also exist in the much simpler spherically symmetric case, obtained by taking W = 0 and F 1 = F 2 in (34) and m = 0 and φ = φ(r) in (35); H 2 = H 3 = 0 and V = V (r), H 1 = H 1 (r) in (36), respectively. In that case, however, it was shown for both the scalar case 281 and the Proca case 277 that no BH solutions exist. Thus, as mentioned above, symmetry (of the metric) noninheritance by the matter fields is a necessary but not sufficient ingredient. A further ingredient is necessary; this can be seen by answering question (2) above.

Stationary clouds and the threshold of superradiance
The answer to question (2) above is "yes". A test field analysis shows the existence of stationary, everywhere regular (on and outside the horizon) solutions of the scalar 209,[282][283][284] or Proca field 277 on the Kerr BH spacetime: stationary scalar or Proca clouds around Kerr BHs. The existence of these stationary clouds is intimately related to superradiance, as we now illustrate for the scalar case.

1641022-26
Numerical relativity and high energy physics: Recent developments Here, M and a are the ADM mass and ADM angular momentum per unit mass of the Kerr solution, and ∆ ≡ r 2 − 2M r + a 2 . Λ is the separation constant, that reduces to the familiar ( + 1) in the Schwarzschild limit. The angular equation defines the spheroidal harmonics. To tackle the radial wave equation, looking for bound state solutions, one requires exponentially decaying solutions towards spatial infinity and a purely ingoing boundary condition on the horizon (in a frame co-rotating with the horizon). Then, one finds in general that the frequency ω is complex: ω = R(ω) + iI(ω). For R(ω) = mΩ H , however, I(ω) = 0 and thus one finds truly stationary states with a real frequency. This condition is interpreted as a zero mode of the superradiant instability, which sets in for R(ω) < mΩ H yielding I(ω) > 0.
This bound state problem becomes particularly simple and elegant for extremal Kerr BHs. 282 In this case the radial equation above, generically of confluent Heun type, reduces to the confluent hypergeometric type, precisely the same equation one finds for the Hydrogen atom (without spin). In this problem, the quantization condition can be interpreted as a condition on the background parameters. Thus, the corresponding stationary clouds -labelled by three quantum numbers (n, , m), where the first is the number of nodes of the radial function and the last two are the spheroidal harmonic indices -can only exist in a subspace of Kerr solutions, actually an one-dimensional existence line, for fixed quantum numbers. This conclusion changes, however, when the test scalar field is allowed to have self-interactions. 285 The Proca case is similar in spirit, but more involved technically, since the Proca potentials do not decouple and no separation of variables has been observed. 272,273 To summarize: the answer to question (1) showed that there is a breach in the wall; the answer to (2) shows that there is indeed something beyond the wall.

Solitonic limits and phenomenology
The construction of Kerr BHs with scalar and Proca hair adapted the technology already in use for (rotating) boson stars. 286,287 Scalar boson stars can be constructed with the ansatz (34)- (35) taking r H = 0, and thus will be a limiting case of the corresponding Kerr BHs with scalar hair. f The Einstein-Klein-Gordon system of equations yields five coupled nonlinear PDEs for the five unknown functions plus two "constraint" equations (which are differentially related to the remaining ones).
f To construct the scalar boson and Proca stars, it is useful to rescale the function W as W/r in (34) and the function H 1 as H 1 /r in (36).
These equations can be solved by a Newton-Raphson relaxation method. 210 Likewise, the Einstein-Proca system, taking the ansatz (34)- (36), yields eight coupled nonlinear PDEs for the eight unknown functions plus two "constraint" equations. Solutions regular on and outside r = r H can be found, and they correspond to Kerr BHs with Proca hair. 277 The r H = 0 limit yields rotating Proca stars, 288 spin-1 cousins of the aformentioned (scalar) boson stars. These observations answer question (3) above.
The exploration of the physical and phenomenological properties of these new families of hairy BHs connected to the Kerr solutions is ongoing research. For the scalar case, it has been noted that the hairy BHs can have quadrupoles and orbital frequency at the ISCO quite distinct from the Kerr case. 209,210 Particularly striking are the BH shadows that have been obtained for some examples of Kerr BHs with scalar hair, with remarkably different shapes and sizes from the Kerr case. 289 These shadows can be partly understood regarding the hairy BHs as composites of a boson star with a horizon, a perspective that can also explain, for instance, the ergoregion structure of these spacetimes. 290 Generalizations of the hairy BHs to include selfinteractions of the matter field have been considered. 291,292 It is likely that similar generalizations are possible in the Proca case.
Finally, let us remark that Myers-Perry BHs with scalar hair have been found in D = 5. These are also anchored to a similar condition between the frequency of the scalar field and the angular velocity of the horizon. 293,294 In D = 5 asymptotically flat spacetimes, vacuum Myers-Perry BHs are not, however, afflicted by superradiant instabilities of massive scalar fields. As such, when the scalar field is set equal to zero, the hairy solutions do not reduce to vacuum Myers-Perry solutions: even though the local geometry can become arbitrarily close to that of the vacuum solutions, there is always a mass gap. A generalization of these solutions including higher curvature terms has also been constructed. 295

Can hairy BHs form?
The existence of Kerr BHs with scalar and Proca hair is theoretically interesting, and it presents us with a rich landscape of previously unknown BH solutions in GR. These solutions require complex bosonic fields, but extremely long-lived solutions exist even for real fields. These solutions describe BHs surrounded by a "cloud" of massive bosons. 296 Are these solutions relevant for astrophysics? The answer to this question depends on two main issues: (A) the existence of massive (and very light) bosonic fields in Nature, and (B) the formation mechanism of these solutions and their stability properties.
Question (A) is an open issue. Question (B) has been studied in a specific scenario. The development of the superradiant instability of massive scalars and vectors has recently been addressed taking into account gravitational radiation, superradiant growth and the effects of a putative accretion disk around the BH, 297 but in an adiabatic approximation (rather than a fully nonlinear numerical evolution).
Assuming that the bosonic cloud is formed through the development of the superradiant instability, it was shown that, within the previous approximations, (i) the observation of supermassive BHs would show gaps in the Reggae-plane, corresponding to BHs which quickly become unstable due to superradiant effects; (ii) the bosonic cloud never backreacts significantly on the geometry; and even though a hairy BH can effectively form, it does not depart significantly from the Kerr geometry. 297 Progress on question (B) has also been achieved using a different toy model: a Reissner-Nordström BH enclosed in a cavity. This system is afflicted by the superradiant instability of bosonic fields (not necessarily massive, since the trapping mechanism is now provided by the cavity) and it was observed that superradiant instabilities in this system -at the test-field level -grow much faster than for Kerr BHs, occurring even for spherically symmetric modes. 262,298-300 These two features make the system tractable with current numerical relativity technology, allowing us to perform fully nonlinear evolutions of the superradiant instability. 301 The simulations showed that the final states of these unstable BHs are indeed hairy BHs at the threshold of superradiance, which can be regarded as the charged counterparts (in this context) of Kerr BHs with scalar hair. 302 Similar results have also been obtained for superradiantly unstable charged BHs in anti-de-Sitter spacetime. 267 Finally, an orthogonal process for the formation of these hairy BHs could be starting from the solitonic limit, rather than the BH limit. It would be nonhairy and very interesting to understand if unstable rotating scalar boson (or Proca) stars could develop into hairy BHs, and how hairy these BHs would be. This is an open issue.

Analog Gravity
Analog models of gravity are a useful tool to investigate kinematical aspects of curved spacetimes in condensed matter systems. 303,304 Analogues have been presented in many contexts, like Bose-Einstein condensates, optical media, and fluids. [305][306][307] Here, we will give emphasis to the progress made in the latter context. Indeed, many interesting physical properties of sonic analogues of BHs have been recently studied, like, for instance, absorption and scattering phenomena, 308-313 as well as quasi-normal modes (QNMs). 314-317

Acoustic analogues
Propagation of sound waves in an ideal fluid, under certain considerations, may be described using the Klein-Gordon equation for a massless scalar field Φ in an effective curved spacetime, namely where g µν are the covariant components of the effective metric (g µν being its contravariant components), with determinant g. We should emphasize that g µν is a 1641022-29 E. Berti et al. function of the local properties of the fluid and, in general, it is not a solution of Einstein's equations.

Canonical acoustic hole
The simplest acoustic analogue to a BH is the so-called "canonical acoustic hole". It consists of a spherically symmetric steady flow of an irrotational barotropic fluid (considered also as incompressible and inviscid), presenting a sink at the origin. It may be described by the following line element: Here where r H is the radius of the sonic event horizon, inside which the radial velocity exceeds the speed of sound c s in the fluid. The canonical acoustic hole is an analogue of the Schwarzschild BH.

Draining bathtub
An acoustic analogue of a rotating BH is the so-called "draining bathtub", whose line element may be written as Here and the constants C and D > 0 stand for the circulation and the draining, respectively. This effective geometry corresponds to the one experienced by sound waves propagating in a fluid with flow velocity The draining bathtub has an ergoregion (defined by the supersonic flow condition |v| ≥ c s ) within the radius r dbt e = √ C 2 + D 2 /c s , and a sonic horizon (defined by v ·r = c s ) at radius r dbt H = D/c s . 310

Hydrodynamic vortex
By setting D = 0 in Eq. (46), we are left with a purely circulating flow, which, in a (3+1)-dimensional setup, may be associated to the following line element

1641022-30
Numerical relativity and high energy physics: Recent developments where r hv e ≡ |C|/c s is the outer boundary of the ergoregion. This effective spacetime is the so-called "hydrodynamic vortex".
In the remainder of this section, we will set c s ≡ 1.

Ergoregion instabilities in acoustic systems
Ergoregion instabilities 318 in acoustic systems have been recently studied for the hydrodynamic vortex, both for incompressible 319 and for compressible fluids. 320 Here, we will review the investigation of instabilities of the hydrodynamic vortex composed by an incompressible fluid. 319 (The numerical results exhibited here are obtained for higher values of the azimuthal number m, complementing the ones exhibited in Ref. 319).
Using the line element (47) in the Klein-Gordon equation (41), and assuming a decomposition of the field Φ as we find the ordinary differential equation where the effective potential V hv m (r) is given by where m is an integer number related with the angular momentum, and ω is the frequency of the perturbation. Solutions describing QNMs can be obtained from Eq. (49), considering the asymptotic behavior at large radial distances u ωm (r → ∞) ∼ exp(iωr) (51) and a boundary condition of Neumann type at r = r min (close to the center of the vortex): This condition is related to a cutoff on the radial velocity increment, i.e. [ ∂Φ ∂r ] r=rmin = 0 [cf. Eq. (48)]. 319 Ergoregion instabilities are present in acoustic systems possessing an ergoregion but not an event horizon. These instabilities are developed inside the ergoregion, i.e. for r < |C|. 319 In order to obtain the QNM frequencies of the hydrodynamic vortex, we can use different numerical techniques to integrate Eq. (49) in the frequency domain. Some results for QNM frequencies are exhibited in Table 1, for different  Table 1. QNM frequencies ω for different values of the azimuthal number m and circulation C = 0.5, obtained numerically from estimates via the DI and CF methods. We impose the asymptotic behavior given by Eq. (51) and a boundary condition of Neumann type, represented by Eq. (52), at r min = 0.51 (outside the ergoregion) and r min = 0.25 (inside the ergoregion).  values of the azimuthal number m and r min , obtained using two different frequencydomain methods, namely the direct integration (DI) and continued fraction (CF) methods. Real and imaginary parts of the QNM frequencies are plotted in Fig. 4, as functions of r min , for different values of the azimuthal number m, obtained using the CF method, considering a circulation C = 0.5. From the results exhibited in Table 1, it can be seen that as the azimuthal number m increases, the magnitude of the real (imaginary) part of the QNM frequencies increases (decreases). Moreover, from the plots exhibited in Fig. 4, we find that as r min decreases, the magnitude of the real and imaginary parts of the QNM frequencies increase (decrease) for unstable (stable) modes. This behavior of the imaginary part can be clearly seen in the inset of the right panel of Fig. 4.

Acoustic clouds
As analogues to the clouds around rotating BHs, described in Sec. 4.3.2, we may have acoustic clouds around the draining bathtub. Taking advantage of the symmetries of the draining bathtub spacetime, characterized by Eq. (44), we can search for

1641022-32
Numerical relativity and high energy physics: Recent developments solutions of the Klein-Gordon equation (41), assuming the separation of variables Φ m (r, φ, t) = e i(mφ−ωt) ζ m (r). (53) The radial function ζ m (r) obeys the ordinary differential equation Using the tortoise coordinate, defined by we can rewrite Eq. (54) as the Schrödinger-like equation where u m ≡ √ rζ m , and we have defined the effective potential Considering the asymptotic limit of Eq. (56), we find the solutions ζ ωm (r) ∼ e −i(ω−ωc)r * , for r → r H , e iωr * , for r → ∞.
In order to have clouds we must choose ω = ω c ≡ mΩ H , and enclose the system inside a "barrier" located at r = r 0 . At the "barrier", we impose suitable boundary conditions, usually chosen to be of Dirichlet or Neumann type. 319 In Fig. 5, we analyze the behavior of the acoustic clouds by plotting the values of the frontier location r 0 as a function of the angular velocity at the horizon Ω H , for different choices of the azimuthal number m. We see that, for a fixed position r 0 of the barrier, the acoustic clouds occur for smaller values of Ω H , as we increase the value of m. Three-dimensional plots of the radial and azimuthal profiles of acoustic clouds are shown in Fig. 6, for Dirichlet (left panel) and Neumann (right panel) boundary conditions.

Concluding Remarks
The fantastic conceptual and formal elegance of Einstein's gravity hides a tremendous complexity when it is applied to realistic, dynamical systems. Quite often, all hope of finding elegant analytic solutions is lost. Then, to tackle this complexity, one needs to resort to numerical solutions. This necessity is now well understood by the scientific community and with the current available techniques, together with the ones under development, there is a strong belief that a lot can be learned about the most elegant physical theory -or generalizations thereof -in the strong field, dynamical regime. We foresee interesting times ahead.