{'fraction_non_alphanumeric': 0.04206971727381524, 'fraction_numerical': 0.023275934353251464, 'mean_word_length': 4.350878981930905, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 37, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Black hole dynamics in generic spacetimes

The dynamics of black hole spacetimes play a crucial role in astrophysics, high energy physics and fundamental physics. In this thesis I have investigated the dynamics of black holes in generic spacetimes by extending established numerical relativity methods to higher dimensional or non-asymptotically flat spacetimes. Additionally, I have explored BH spacetimes perturbed my massive scalar fields. By developing and further improving NR techniques I have been able to push our knowledge to new grounds.

First of all, I would like to warmly thank my PhD supervisor Prof. Vítor Cardoso for many illuminating discussions and all the physics that I've learned in the last four years. I am grateful for your constant support and encouragement and the many opportunities that you have made possible. It has been a great pleasure to work with you. You are an exceptional supervisor and I could not have wished for a better one. Thank you, Vítor.
I am grateful to my closest collaborators Leonardo Gualtieri, Carlos Herdeiro, Andrea Nerozzi, Ulrich Sperhake, and Miguel Zilhão, for our many fruitful discussions on this exciting project. It has been a great pleasure to collaborate with you.
I would like to thank my long-time collaborator David Hilditch for our many enlightening discussions and interesting projects.
I am indebted to Andrea Nerozzi and Paolo Pani for carefully proof-reading this thesis. It is my pleasure to thank all members of CENTRA and especially our gravity group who provided an exceptional and inspiring working environement. Special thanks to Dulce Conceição and Rita Sousa who have been a constant support and never got tired to help me through the burocracy jungle. I would like to thank Sergio Almeida for building up and maintaining our cluster "Baltasar" and for being a great aid in all computer related problems.
I am thankful to my office-mates Andrea Nerozzi, Jorge Rocha and Antonino Flachi for our many inspiring discussions and for creating such a great working environment.
I would like to thank Emanuele Berti, Terence Delsate, Akihiro Ishibashi, Hirotada Okawa, Paolo Pani and Jan Steinhoff for many fruitful and interesting dicussions.
I acknowlegde the kind hospitality of the TAPIR group at Caltech, the gravitational physics group at OleMiss, the University of Mississippi, the Numerical Relativity group at the Friedrich-Schiller University Jena, and the Yukawa Institute for Theoretical Physics at Kyoto University, where parts of this work have been accomplished.
I warmly thank my family, Kerstin and Johannes and Manfred and Helga Witek, for your constant support, inspiration and unconditional love. I would not stand where I am today without you.
I am happy to thank my long-standing friends Tine, Anke and Stephan for your unconditional friendship and support.

Introduction
Black holes (BHs) are among the most intriguing predictions of General Relativity (GR) and alternative theories of gravity. They play a key role in many areas of physics reaching from astrophysics, cosmology and gravitational wave astronomy over high energy physics and the gauge/gravity duality to fundamental mathematical questions concerning their existence and stability properties. Of utmost interest are the dynamics of BHs in the fully non-linear, strong-field regime of GR. Because of the complex structure of Einstein's equations, a set of D(D + 1)/2 coupled, non-linear partial differential equations (PDEs), this regime can only be explored by employing Numerical Relativity (NR) techniques. Since the outstanding breakthrough in 2005 and 2006 [10][11][12] NR has grown into a mature research field capable to model a plethora of exciting phenomena (see e.g. [13][14][15][16][17][18][19][20] for recent reviews).
In this thesis I have explored BH dynamics in generic spacetimes, including BH collisions in four and higher dimensional asymptotically flat spacetimes, BH inspirals in a box, thus mimicking asymptotically Anti-de Sitter (AdS) spacetimes, and the stability of BHs against massive scalar field perturbations. In the following I give a brief overview and motivation for each of these topics.
Astrophysics and gravitational wave astronomy Since the first discovery of BHs as solutions to Einstein's Equations by Schwarzschild in 1916 they have outgrown their status of merely exotic theoretical predictions of GR and are nowadays known to make up an important class of astronomical objects. Supermassive BHs (SMBHs) with 10 6 − 10 9 M are expected to be hosted at the center of most galaxies [21][22][23]; in fact there is strong observational evidence for our very own SMBH in Sagittarius A * with M ∼ 4.2 · 10 6 M in the centre of the Milky Way. Their solar-mass counterparts with 3 − 30M are part of the galaxies population [24][25][26]. Their motion is dominated by gravitational interactions and the collision of two BHs is among the most violent scenarios in Nature, releasing a burst of energy in form of gravitational radiation. Moreover, gravitational waves emitted throughout such a collision provide the unique opportunity to explore and test the strong-field regime of GR or extensions thereof -an as yet still open issue in gravitational physics. Once detected, gravitational waves will open up a completely new window for astronomical observations, thus offering invaluable astrophysical insight into strong-field phenomena of gravity that are accessable in no other way. First detections of gravitational waves are expected within the next decade. Gravitational wave observatories, such as pulsar timing arrays [27][28][29] or ground-based interferometric GW detectors such as LIGO [30], Virgo [31] and GEO600 [32] are in place. The second generation, advanced LIGO-VIRGO detector network is expected to be up and running in 2015 [33,34] while KAGRA [35,36] is under construction. Third generation ground-based detectors such as the Einstein Telescope [37] are under development and space-based detectors in the spirit of the eLISA/NGO [38] operation might become available in the future. Although BH binaries are suspected to be among the strongest sources for GWs, their detection poses an extremely challenging task. The expected GW signal arriving at Earth is extremely weak; the dimensionless strain, or amplitude, is expected to be of the order h = 10 −20 . For the LIGO detectors with an arm length of L ∼ 3km this strain translates into a variation in length of ∆L ∼ 10 −17 m -two orders of magnitude smaller than the size of a proton. In order to facilitate the measurement of these very weak signals buried under all kinds of noise, such as seismic noise from the Earth, thermic noise of the test masses and radiation pressure from the laser light, matched-filtering techniques are mandatory. This technique aims at filtering real signals out of the noise background via the comparison with theoretical templates. Therefore a detailed theoretical understanding of the GW signal from the inspiral, merger and ringdown phase is of utmost importance. While the early stages of an inspiraling BH binary, when v/c 1, can be modelled well by analytic post-Newtonian (PN) and effective one-body (EOB) approaches and the final, after-merger ringdown phase is approximated well by perturbative approaches, both (semi-analytic) methods break-down in the highly dynamical merger phase. In order to accurately model the actual collision between the two BHs, the fully non-linear system of Einstein's equations has to be solved by employing numerical relativity methods.
Since the remarkable breakthrough in 2005 [10] followed by [11,12] shortly afterwards, NR has grown into a mature research field now able to evolve numerous configurations of BH binaries.
In Chapter 2, I will give a brief overview on the most commonly employed techniques in NR in 4-dimensional, asymptotically flat spacetimes -the setup of choice to simulate astrophysical BH configurations. In particular, I will briefly summarize the 3+1-splitting of spacetime into 3-dimensional spatial hypersurfaces that will be evolved in time. Subsequently I review the 3 + 1-decomposition of Einstein's equations and their particular formulation yielding (numerically) long-term stable simulations, the construction of initial data and the extraction of physically meaningful quantities.
In Chapter 3 I present numerical simulations of a particular BH binary setup. Specifically, I have evolved head-on collisions of unequal mass BHs with a mass ratio varying from q = 1, . . . , 1/100. Thus, I have been able to fill the gap between point particle calculations valid in the small mass ratio regime q 1 and a fully numerical modelling of the intermediate to comparable mass ratio regime.
High energy physics In recent years a further exciting application of BH physics in the highly dynamical non-linear regime has come into focus of the NR community: the intriguing possibility of BH production in high-energy particle colliders or interactions of ultra-high energetic cosmic rays with the atmosphere. If particles collide at very high energies such that the center of mass energy is well beyond the Planck scale, gravity becomes the dominant interaction. Then, due to Thorne's hoop conjecture [39] and "no-hair" theorem like arguments [40], the particular nature of these particles should not be important to describe their collision or, figuratively speaking, matter does not matter . A fully non-linear, numerical study by Choptuik & Pretorius [41] and East & Pretorius [42] has provided compelling evidence for this statement. The authors have simulated ultra-relativistic particle collisions by the collision of, respectively, bosonic solitons and fluid balls with generic initial boost and observed BH formation above a critical boost parameter of γ c = 2.9 in the former and γ c = 8.5 in the latter case. Therefore, the trans-Planckian scattering of particles is well described by scattering processes of highly boosted black holes. This fact gains a lot of importance in the context of so-called TeV gravity scenarios. In order to solve the hierarchy problem, higher dimensional theories of gravity have been proposed, in which the fundamental Planck scale is of the order of TeV [43][44][45][46][47]. This offers the exciting possibility that BHs could be produced, e.g., at the Large Hadron Collider (LHC) or in ultra high energy cosmic ray interactions with the atmosphere [48][49][50][51][52][53][54][55][56][57].
Recent searches for micro-BHs at Atlas [58] and CMS [59] put new constraints on BH production and excluded semi-classical or quantum BHs with M ∼ 4.3 . . . 6.3T eV [60].
Of particular importance is the computation of the energy released in form of gravitational radiation during the collision, which in turn determines the mass and angular momentum of the newly born BH. A second, crucial parameter is the cross-section of BH collision to provide estimates on the BH production rate. Lower bounds on the BH mass have been obtained from area theorem arguments [61,62]. Estimates on the radiated energy have been provided by computations of (highly) relativistic point particles falling into a higher dimensional BH [63,64], as well as shock wave collisions describing the ultra-relativistic regime [65][66][67][68]. Instead, I have aimed at complementary calculations of the energy emitted throughout the collision of BHs in higher dimensional spacetimes. In particular, I have attempted the challenging goal to evolve numerically collisions of comparable mass BHs in the fully dynamical, non-linear regime, because proton-proton collisions are the most likely type of collisions yielding BH formation. In Chapter 4 I present the framework based on the dimensional reduction and the numerical implementation enabling to explore higher dimensional BH spacetimes numerically. In Chapter 5 I summarize our results of these simulations, in particular benchmark tests of the code as well as head-on collisions of equal-mass BHs in D = 5 and D = 6 spacetime dimensions and of unequal-mass BHs with comparable mass ratios in D = 5.
Gauge/gravity duality Another fascinating application of BH physics is provided by the gauge/gravity duality [69,70]. Loosely speaking, this correspondence relates a theory of gravity in an Anti-de Sitter (AdS) spacetime to a gauge theory "living" on the boundary of this AdS spacetime. For example, it has been found that the Hawking-Page phase transition of black holes in AdS, i.e., the transistion between thermally stable and unstable BHs [71], is dual to the deconfinement/ confinement transition in a strongly coupled, thermal gauge theory [70]. Therefore, the analysis of BHs in AdS provides an appealing tool to study particles in a gauge theory. Eventually, these calculations may find applications in the description and understanding of high energy collisions of heavy ions in the Relativistic Heavy Ion Collider (RHIC). In the past couple of years first (fully) numerical works on the gravity side, i.e., of asymptotically AdS spacetimes have been presented. These studies include: (i) simulations of shock wave collisions in AdS [72,73], (ii) evolutions of scalar fields in AdS which has lead to turbulent-like instabilities and BH formation [74] or time-periodic solutions [75] and (iii) numerical studies of prompt BH formation through the ringdown phase which results in a static BH [76].
Because of the active role played by the AdS boundary, perturbations of (small) Kerr-AdS BHs are expected to give rise to superradiant instabilities [77][78][79]. It has been shown in [78] that this superradiant instability can be understood as "BH bomb" like instability first suggested by Press & Teukolsky [80] and, in a different setup, by Zel'dovich [81,82].
The key idea is that a wave packet with frequency ω impinging on a Kerr BH gets amplified due to superradiant scattering if its frequency satisfies the superradiant condition ω < mΩ H , where m is the azimuthal wave quantum number and Ω H the angular momentum of the BH horizon. Now, in the original gedanken experiment [80][81][82], the system has been surrounded by a perfectly reflecting cavity, at which the wave is reflected. Upon subsequent amplification due to superradiant scattering and reflection at the mirror, the field grows without bound and triggers the system to become unstable.
Whereas this setup at first glance appears somewhat artificial, such a "mirror" is naturally provided by the AdS boundary, thus yielding the superradiant instability of small Kerr-AdS BHs [77][78][79].
Furthermore it has been stated [78] that the global BH-AdS spacetime can be mimicked by an asymptotically flat BH spacetime in a finite-size box , which plays the role of the active boundary.
Because asymptotically AdS spacetime pose a challenging problem I have started with a toy model in order to study the effects of the boundary. Therefore, I have exploited the aforementioned relation and performed numerical simulations in 3 + 1-dimensions of BH binaries surrounded by a mirror-like box. The framework and results are presented in Chapter 6. The goal of this project has been twofold: (i) it has represented first fundamental steps towards actual simulations of fully dynamical BH-AdS spacetime; and (ii) it has offered the unique potential to study, fully non-linearly, a "BH bomb" like setup [80] and the interaction of a rotating BH with gravitational perturbations.

Stability of BHs in astrophysical systems
The "BH bomb" scenario in real physical setups cannot only be achieved by asymptotically AdS spacetimes but also by massive fields surrounding Kerr BHs. In this case, the mass term µ S effectively acts as a mirror if ω ≤ µ S , where ω is the frequency of the field [83][84][85][86][87][88][89][90].
Whereas this effect is very weak for known standard model particles in astrophysical environments, it might become astrophysically significant if there exist fields with small, but non-vanishing mass. Fascinating candidates are so-called axions. Axions are ultralight bosonic states emerging in string theory compactifications, which have not been ruled out by experiments. Recently the "axiverse" scenario has been suggested, implying the existence of a plethora of axionic particles covering a mass range from 10 −33 eV ≤ µ S ≤ 10 −8 eV [91][92][93]. The existence of such ultra-light bosonic states might have a vast number of astrophysical implications and observational signatures. These implications include modifications of the cosmic microwave background (CMB), or modification of a gravitational wave signal due to the creation of "gravitational atoms" -bosonic clouds in the vicinity of BHs. If the axion mass is in the range 10 −22 eV ≤ µ S ≤ 10 −10 eV , the time scales for the superradiant instability becomes astrophysically significant and its interaction with astrophysical BHs is expected to result in observable signatures [91][92][93][94][95][96][97][98]. Thus, axions have the fascinating potential to facilitate high precision BH physics [91][92][93].
With this plethora of applications in mind, I have investigated massive scalar fields in 4-dimensional, asymptotically flat BH backgrounds, which I present in Chapter 7. The main motivation to revisit the superradiant instabilities of massive scalar fields, that has been studied extensively in the frequency domain [80,83,87,88,94,97,99,100], in the time domain has been twofold: (i) By developing the code presented in Chapter 7 to evolve massive scalar fields in BH backgrounds, I have accomplished first important steps towards more generic setups. These scenarios include non-linear self-interaction terms of the scalar field as well as a fully non-linear evolution modelling backreaction effects and are subject to future investigations. (ii) Before the most recent studies by Yoshino & Kodama [97] and Dolan [101] there has been only one attempt to investigate the massive scalar field instability in the time domain by Strafuss & Khanna [102]. Some of their conclusions, however, are in stark contradiction to previous computations in the frequency domain [87,88]. I have found a mechanism that has the potential to explain these contracdictory results. I will discuss the underlying effect in detail in Chapter 7 based on Ref. [9].

Review of Numerical Relativity in D = 4
The basis of most NR schemes is the reformulation of Einstein's equations as time evolution problem. At the core of this formulation lies the splitting of spacetime into spacelike hypersurfaces labelled by a time parameter. In this chapter I give a brief review of the key concepts of the 3 + 1-approach NR. For the sake of illustration I focus on 4dimensional spacetimes and will discuss modifications of the approach to For a detailed overview and summary of the subject I refer the reader to [13][14][15][16][103][104][105][106][107] and references therein.
Whereas spacetime is usually considered as an entity, where space and time are on equal footing, it is mandatory to split them in order to investigate the time evolution of a physical system with some given initial configuration. There are various alternative ways to perform this splitting, such as the characteristic approach, in which the spacetime is foliated by null hypersurfaces (see e.g. [108] for a review) or the generalized harmonic formulation, in which the Einstein's equations are formulated as a set of wave equations together with harmonic coordinates [109,110]. The most commonly employed approach in NR is the 3 + 1-approach first introduced by Arnowitt, Deser and Misner [111] and later reformulated by York [103]. In this approach the equations of motion (EoM) are formulated as a Cauchy problem, i.e., as an initial (boundary) value problem.
In this chapter, I will summarize the 3 + 1-decomposition of the spacetime and Einstein's equations in Sec. 2.1. In Sec. 2.2 I give an overview on techniques to construct initial data for BHs spacetimes within this framework. In Sec. 2.3 I review the most common gauge choices describing the employed coordinate system. In Sec. 2.4 I will summarize the generalized Baumgarte-Shapiro-Shibata-Nakamura scheme of the 3 + 1form of Einstein's equations, which represents a well-posed initial value problem of the system. Finally, in Sec. 2.5, I summarize techniques to extract physically meaningful quantities.

3+1-decomposition 2.1.1. Foliation
Any globally hyperbolic (here 4-dimensional) spacetime M can be foliated into purely spatial, 3-dimensional hypersurfaces Σ t labeled by the timelike parameter t, as depicted in Fig. 2.1. Then, the geometry of this region of spacetime is determined by the induced 3-metric γ ij , which measures the proper distance dl 2 = γ ij dx i dx j within the spatial hypersurface Σ t , the lapse function α, which measures the elapsed proper time as measured by an observer moving along the normal vector n µ between two hypersurfaces (Eulerian observer), and the shift vector β µ = (0, β i ) which determines the relative velocity be- : Illustration of the foliation of a spacetime M into spatial hypersurfaces Σ t and Σ t+δt . Depicted is also the vector t µ = αn µ + β µ pointing from a point p ∈ Σ t to a point p ∈ Σ t+δt on the next hypersurface.
tween the Eulerian observer and lines of constant spatial coordinates. The lapse function α and the shift vector β i encode the information about the coordinates of the spacetime and are a manifestion of the coordinate degree of freedom in GR, thus freely specifiable. The vector t µ pointing from a point p in Σ t to a point p on the next slice Σ t+δt is constructed from where n µ is the vector normal to the spatial hypersurface Σ t with n µ =(−α, 0, 0, 0) , n µ = 1 α (1, −β i ) , n µ n µ = −1 . (2.2) Note, that we have choses the last relation such that n µ is timelike. In other words, the lapse function α relates the coordinate time t that labels the slices of the foliation to the time measured by an Eulerian observer. Then, the expression of the 4-acceleration a µ of this observer in terms of the lapse writes a µ =n ν ∇ ν n µ = 1 α D µ α = D µ ln α . (2.3) The line element of the 4-dimensional spacetime is obtained from From the line element (2.4) I can read of the components of the spacetime metric g µν in terms of the 3-metric γ ij , the lapse α and shift β i (2.5) The 3-metric γ ij is defined as the metric induced on the spatial hypersurface Σ t as γ µν =g µν + n µ n ν , γ µν = g µν + n µ n ν .

Intrinsic and Extrinsic curvature
The intrinsic curvature of the spatial hypersurface Σ t is determined by the 3-dimensional Riemann tensor (3) R k lij . It measures the non-commutativity of two successive (spatial) covariant derivatives D i associated with the 3-metric γ ij and is expressed by the Ricci identity for any spatial vector v k ∈ Σ t . In contrast, the extrinsic curvature K ij describes how a spatial slice Σ t is embedded into the spacetime manifold M, i.e., how the direction of the normal vector n µ changes as it is transported along the hypersurface Σ t , as illustrated in Fig. 2.2. Then, the extrinsic curvature is defined as K µν = − γ λ µ ∇ λ n ν = −∇ µ n ν − n µ n λ ∇ λ n ν = −∇ µ n ν − n µ a ν . (2.14) One can show that the extrinsic curvature is symmetric K µν = K νµ and purely spatial, i.e., K µν n µ = K µν n ν = 0. Hence, from now on I will only consider its spatial components K ij . Furthermore, the extrinsic curvature K ij is related to the Lie derivative of the 3metric γ ij along the normal vector n µ . If I apply the definition of the Lie derivative of a tensor field T µ 1 ...µp ν 1 ...νq along a vector u λ given by L u T µ 1 ...µp ν 1 ...νq =u λ ∇ λ T µ 1 ...µp ν 1 ...νq − T λ...µp ν 1 ...νq ∇ λ u µ 1 − ... + T µ 1 ...µp λ...νq ∇ ν 1 u λ + ... , (2.15) to the spatial metric I obtain L αn γ µν =αn λ ∇ λ γ µν + γ λν ∇ µ (αn λ ) + γ µλ ∇ ν (αn λ ) = −2αK µν . (2.16) This expression provides yet another interpretation of the extrinsic curvature as the "velocity" or "time derivative" of the spatial metric γ ij as seen by an Eulerian observer. Because the relation (2.16) has been derived only from geometrical concepts it is purely kinematic. The dynamics of the system will come into the game by considering the 3+1decomposition of the EoM provided by Einstein's equations, which yields an evolution equation for the extrinsic curvature K ij . : Illustration of a spatial hypersurface Σ t . The normal vector n µ changes its direction when it is displaced on the the hypersurface thus yielding a non-vanishing extrinsic curvature.

Projections of the Riemann tensor and the Gauss-Codazzi equations
Before applying the 3+1-split to Einstein's equations I first consider the different projections of the 4-dimensional Riemann tensor (4) R µ νρσ . The possible, non-zero projections of the 4-dimensional Riemann tensor (4) R µ νρσ include the full (spatial) projection onto the hypersurface, the contraction of one index with the vector normal to the hypersurface and the spatial projection of the three remaining indices as well as the projection twice along the normal and onto the hypersurface. All remaining combinations of projections vanish due to the symmetries of the Riemann tensor. First, I consider the second (spatial) derivative of a purely spatial vector v µ , which can be rewritten as which determines the Riemann tensor fully projected onto the spatial hypersurface. Next, I contract the Riemann tensor once with the vector normal to the spatial hypersurface Σ t and project the remaining three free indices onto the hypersurface. By employing the Ricci identity, Eq. (2.13), and the definition of the extrinsic curvature, Eq. (2.14), I obtain which yields Finally, I contract the Riemann tensor twice with the normal vector n µ and project it twice onto the hypersurface Σ t . This operation leads to where (4) R µν and (4) R are the 4-dimensional Ricci tensor and scalar, respectively, and T µν is the stress-energy tensor. In order to perform the 3 + 1-decomposition of Einstein's equations I write the EoM (2.22) in the form By contracting Eq. (2.23a) twice with the vector n µ normal to the spatial hypersurface, i.e., E 1,µν n µ n ν = 0, and inserting the Gauss relation (2.18) I get the Hamiltonian constraint where ρ = T µν n µ n ν is the energy density. If I project Eq. (2.23a) once onto the hypersurface Σ t and once perpendicular to it, i.e., γ λ µ E 1,λν n ν = 0, and substitute the 4-dimensional, projected Ricci tensor by the Gauss-Codazzi relation (2.20) I obtain the momentum constraint where j µ = −γ λ µ T λν n ν is the energy momentum flux. One can show straight-forwardly that the momentum density is purely spatial j µ n µ = −γ λ µ T λν n ν n µ = 0. Finally, I consider the fully spatial projection of Eq. (2.23b), γ κ µ γ λ ν E 2,κλ = 0, and insert Eqs. (2.18) and (2.21) to eliminate the Ricci tensor. Performing this operation yields where S ij = γ κ i γ λ j T κλ is the purely spatial matter stress tensor. By re-organizing the terms and applying Eq. (2.1) I obtain the time evolution equation for the extrinsic curvature In summary, the full set of Einstein's equations (2.22) in 3 + 1-form are provided by the time evolution equations for the 3-metric, Eq. (2.28), and extrinsic curvature, Eq. (2.27), as well as the Hamiltonian and momentum constraints, Eqs. (2.24) and (2.25), respectively. In the NR community, this particular set of equations is usually referred to as ADM equations, although they are not in the Hamiltonian form as originally introduced by Arnowitt, Deser & Misner [111] but rather in the form derived by York [103]. In the following sections I will only consider 3-dimensional quantities and will therefore drop the superscript (3) [...].

Initial Data
The 3+1-decomposition of Einstein's equations results in a set of evolution equations for the spatial metric γ ij and the extrinsic curvature K ij , Eqs. (2.28) and (2.27), respectively, and a set of constraints to the system which have to be satisfied for all times. In a free evolution scheme, the Hamiltonian and momentum constraints, Eqs. (2.24) and (2.25), are solved only to provide initial data for the 3-metric and curvature, thus representing the initial configuration for the physical system of interest. Although conceptually straightforward, the initial data construction is a non-trivial task. The constraint equations are a set of four coupled, elliptic PDEs which, in general, are difficult to solve. In the following I briefly summarize one of the most common procedures to solve the constraint equations, namely the conformal decomposition first introduced by Lichnerowicz and York [112][113][114]. For further methods and more detailed discussions I refer the reader to [105,106,115] and references therein.

York-Lichnerowicz decomposition and CTT approach
In the conformal transverse traceless approach for the construction of initial data I first conformally decompose the (physical) 3-metric γ ij into the conformal metricγ ij which I assume to be given and the conformal factor ψ. Note, that with this rescalingγ = detγ ij = 1. The extrinsic curvature is split into its trace K and tracefree part A ij Additionally, the tracefree part A ij of the extrinsic curvature is conformally rescaled according to whereD i andR denote the covariant derivative and Ricci scalar with respect to the conformal metricγ ij . The momentum constraint (2.25) becomes In order to solve the momentum constraint (2.33), I further split the traceless part of the extrinsic curvatureÂ ij into a transverse and longitudinal part [103,116] whereÂ ij T T is transverse and traceless, i.e., The longitudinal part (LX) ij is a conformal Killing form associated with the conformal metric acting on a vector field X i according to Then, the constraints reduce tõ Assuming that the conformal metricγ ij , the trace of the extrinsic curvature K and the transverse traceless part of the extrinsic curvatureÂ ij T T are given, the constraint equations (2.37) and (2.38) can be solved for the conformal factor ψ and the vector X i and the physical quantities (γ ij , K ij ) can be reconstructed on the initial time slice.

Initial data for multiple BH systems -the puncture approach
In this section I apply the previously discussed approach to solve the constraint equations to the case of BH spacetimes in which I are most interested in the present work. Therefore I consider the 3 + 1-Einstein equations in vacuum, implying that the matter variables (ρ, j i , S ij ) vanish. The currently most popular method to implement BH configurations is the moving puncture approach [11,117,118]. The initial data shows a Brill-Lindquist wormhole topology [119], depicted in Fig. 2.3. The spatial slice consists of two or more copies of R 3 with points r i removed. Thus, one obtains several asymptotic flat regions connected by Einstein-Rosen bridges and each asymptotically flat end is compactified and identified with a single point r i on R 3 . These coordinate singularities at the points r i are referred to as "puncture". In more practical terms, one uses an isotropic radial coordinate r iso = m 2 4r , such that the sphere r iso = m 2 remains invariant. Puncture data enables us to represent BH spacetimes without excision. The singularities are not cut off the numerical grid during the evolution but are avoided by the gauge choice. In the case of spacetimes containing N BHs each singular point represents infinity in a different asymptotically flat region, so that our universe is connected with N different universes through Einstein-Rosen brigdes (wormholes). Notice, that Eqs. (2.37) and (2.38) still constitute a system of four coupled PDEs. Therefore, I apply further simplifications to the initial data problem. In particular, I consider Brill-Lindquist type initial data, which describes a system of N BHs momentarily at rest on the initial time slice, and Bowen-York initial data, which sets up a system of N BHs with linear momentum and spin. Brill-Lindquist type wormhole topology of a spacetime containing two BHs. In this case our universe is joined to two distinct univeres via wormholes.
Brill-Lindquist type initial data I first consider time symmetric initial data, implying K ij = 0, in a conformally flat spacetime, i.e.,γ ij = η ij . Then, the momentum constraint is trivially satisfied and the Hamiltonian constraint reduces to a Laplace equation in flat space for the conformal factor ψ 0 at t = 0 whereD i is the flat-space covariant derivative. The boundary condition for an asymptotically flat spacetime requires ψ = 1 at infinity. Then, the ansatz for a single Schwarzschild BH is where m is the bare mass parameter of the BH. Because the Laplace equation (2.39) is linear, I obtain the solution representing a spacetime of N BHs momentarily at rest by superposition where m i is the bare mass parameter of the i th BH located at r i . This solution is known as Brill-Lindquist initial data [119,120].
Bowen-York initial data Whereas Brill-Lindquist type data sets up time symmetric data representing a spacetime of N BHs momentarily at rest, this assumption is relaxed in the case of Bowen-York initial data. Thus, it will set up more realistic configurations of BHs with linear momentum P i and spin S i . In this approach I assume K ij = 0, employ the flat metric as the conformal metric and impose the maximal slicing condition, i.e., where ψ 0 denotes the conformal factor on the initial time slice. Then, the Hamiltonian and momentum constraints decouple. Specifically, the latter becomes A ij P S(a) = 3 2r 2 P i a n j + P j a n i − (g ij − n i n j )P k a n k + 3 r 3 ikl S ka n l n j + jkl S ka n l n i .
(2.45) P i a and S ia are identified as the momentum and the spin of the a th BH, respectively. This can be seen by computing the ADM linear and angular momentum at infinity (see, e.g., [106] and references therin) These integrals are performed over spheres of constant radius r, where u i denotes the normal vector to the sphere. If I now consider the conformal factor at infinity ψ ∞ = 1 and insert Eq. (2.45), I find that P i a and S ia are indeed the linear and angular momenta of the spacetime.
The Hamiltonian constraint is given by Eq. (2.32) with K = 0 andÂ ij computed according to Eq. (2.44). The solution is given by the ansatz [117] where u is a regular function which is C 2 at the puncture and C ∞ elsewhere. The singular behaviour is absorbed in the term ψ BL . Inserting ansatz (2.47) into the Hamiltonian constraint (2.32) yields an elliptic equation for the function u. After finding a solution for u, the conformal factor and, thus, the physical metric γ ij can be reconstructed whereas the extrinsic curvature is given by Eqs. (2.42) and (2.44), describing the initial configuration of a system of N BHs each with linear momentum P i and spin S i .

Gauge conditions
As we have seen in Sec. 2.1, the coordinates, which are encoded by the lapse function α and the shift vector β i can be choosen freely. Nevertheless, for the sake of applications such as numerical simulations they need to be specified in order to close the system of equations. In the following we give a brief review of the most common choices of the slicing condition fixing the lapse function α and the shift condition prescribing β i .

Choice of foliation
Geodesic Slicing The simplest option for the choice of foliation is given by the geodesic slicing condition α = 1. This choice of the lapse function α yields a vanishing 4acceleration a µ of the Eulerian observer as can been seen from Eq. (2.3) implying that the worldlines of this observer are geodesics. However, from a practical viewpoint this choice of the foliation is a poor one because in a Schwarzschild spacetime every Eulerian observer starting from rest will freely fall into the singularity at a finite time, thus resulting in a break-down of any numerical simulation.
Maximal Slicing In order to avoid the pathologies of the geodesic slicing condition we impose a condition on the trace of the extrinsic curvature rather than on the lapse function α. Specifically, we choose K = 0, which is known as maximal slicing condition.
In strong field regions this condition will tend to hold back the evolution of the slice such that the proper time increases "more slowly". Whereas the entire spacetime outside the event horizon is covered by this foliation the slices pile up near the BH region such that they never reach the singularity. This behaviour is referred to as singularity avoidance. Consequently, the proper time of an Eulerian observer between two hypersurfaces tends to zero with increasing (coordinate) time t. This implies that α goes to zero as t goes to infinity and is called collapse of the lapse.
Harmonic Slicing The harmonic condition for spacetime coordinates (x µ ) is given by g x µ = 0, where g = ∇ µ ∇ µ is the d'Alembertian associated with the spacetime metric g µν . Imposing the harmonic condition only on the time coordinate x 0 while the remaining spatial coordinates are allowed to be choosen freely yields the harmonic slicing condition Like the maximal slicing condition the harmonic slicing also proves to be singularity avoiding.
A further succesful approach is achieved by employing generalized harmonic coordinates, in which the wave equation for each coordinate is allowed to contain a source term, i.e. g x µ = H µ . These gauge-driving source terms H µ can either be specified algebraically or evolved in such a way that hyperbolicity is preserved [110,[122][123][124]. The first successful numerical simulations of orbiting BH binaries considered a source term for the time coordinate that effectively kept the lapse close unity, while the spatial coordinates remained harmonic [10,110]. This was accomplished by evolving the source term itself, according to where ξ 1,2 are constants.
1 + log-Slicing The currently most common slicing condition used in many 3 + 1 numerical relativity codes is the 1 + log-slicing condition defined by In normal coordinates, i.e., for vanishing shift vector β i = 0, this equation can be rewritten as and yields hence the naming 1 + log-slicing. This slicing condition has strong singularity avoiding properties and has been found to mimic the maximal slicing condition while simultaneously being computationally cheaper [125].

Generalization of Slicing Conditions -Bona-Massó Family of Slicing Conditions
The slicing conditions introduced so far have been generalized by Bona et al. [126] and later further modified by [127] to where f (α) is an arbitrary, positive funtion of the lapse α and K 0 is the value of the trace of the extrinisic curvature on the initial time slice. In the case of f (α) = 0 and α = 1 the Bona-Massó slicing condition (2.53) reduces to the geodesic slicing. The harmonic slicing corresponds to the choice f (α) = 1, whereas the prescription f (α) = 2 α leads to the 1 + log-slicing.

Choice of spatial coordinates
Minimal Distortion The intention of the minimal distortion condition is the construction of a shift condition that minimizes the change of the conformal metric in time [128,129]. Therefore one introduces the distortion tensor Q ij that is defined as Because Q ij is tracefree it does not take into account the change of the volume of this domain but only its shear. Then, minimizing the integral of Q ij Q ij over the spatial hypersurface with respect to the shift gives the minimal distortion condition which implies D j ∂ tγ ij = 0. In the weak field limit this condition includes the standard transverse traceless gauge of the linearized theory [128]. Smarr and York [128] introduced the so-called radiation gauge which combines the minimal distortion with the maximal slicing condition K = 0, thus describing radiative spacetimes.
If we introduce the conformal connection functionΓ i =γ jkΓi jk = −∂ jγ ij that will play a key role in the BSSN formulation introduced in Sec. 2.4 below, we can relate the minimal distortion condition to the so-called Γ-freezing condition modulo lower order terms ∼ ∂ kγij .
Gamma-Driver Conditions Because the minimal distortion and Γ-freezing conditions are of elliptic type, they are computationally expensive to solve at each timestep. Therefore, modifications thereof have been considered, that are of hyperbolic type and can straightforwardly be implemented as time evolution equations. Of particular relevance for numerical implementations using the BSSN evolution system is the Γ-driver condition. The Γ-driver condition is a modification of the Γ-freezing condition first introduced in [127] and is given by whereΓ i =γ jkΓi jk = −∂ jγ ij is the conformal connection function, ξ Γ is a positive scalar function, η β is a damping parameter and B i is an auxiliary function. It has been found crucial to add the dissipation term η β B i in order to avoid strong oscillations in the shift [127]. This system of equations can be integrated to [130] where we have added the advection term β k ∂ k β i . The Γ-driver condition in one of the presented forms together with the 1 + log-slicing, Eq. (2.50), are known as moving puncture approach or also puncture gauge. The puncture gauge has been widely used in NR codes employing the so-called BSSN evolution system, described in the next section.

Generalized BSSN formulation
For a long time, numerical methods based on the ADM equations (2.24), (2.25), (2.27), and (2.28), suffered from a number of problems, most prominently from being numerically unstable (in 3 + 1) allowing only for short evolutions of BH spacetimes. In fact, a detailed analysis of the properties of the underlying PDE system reveals that the ADM equations are in general only weakly hyperbolic and therefore exhibit an ill-posed initial value problem (see e.g. [106,107,131] and references therein). Therefore, a large variety of alternative or modified versions of the evolution equations, suitable choices of the coordinate gauge and different treatments of the BH singularities have been investigated. A plethora of different well-posed formulations of the 3 + 1-Einstein's equations have been proposed, including the BSSN system [132,133] and generalizations thereof [134,135], as well as reductions to a first order BSSN version [136], the Nagy-Ortiz-Reula (NOR) [137] and Z4 scheme [138,139], a conformally decomposed version of Z4 [140][141][142][143], the Bona-Massó formulation [126,144,145], the Kidder-Scheel-Teukolsky family [146,147] and a Hamiltonian formulation [148], to name but a few. The currently most prominent and succesful techniques include the Generalized Harmonic (GH) formulation employed in the original breakthrough in 2005 [10,110], in which Einstein's equations are cast into a set of wave equations together with the generalized harmonic gauge, and the moving puncture technique which yielded first succesful simulations of BH binary spacetimes shortly afterwards [11,12]. The latter method makes use of the BSSN formulation of the 3 + 1-Einstein's equations [132,133] and the puncture gauge [11,12,127].
In this section we wish to review the key aspects of a generalized BSSN formalism introduced and studied in great detail in [134] based on the BSSN system [132,133] and modifications first considered in [135]. Whereas the BSSN equations traditionally are introduced as a conformal decomposition and readjustment of the ADM equations (see e.g. [127,133,149]), we here follow a slightly different approach [134]. Specifically, we consider the addition of physical and definition-differential constraints which alter the characteristic structure of the ADM equations, thus guaranteeing well-posedness of the system, accompanied by a conformal decomposition changing the evolved variables into a form convenient for numerical simulations.

Constraint addition
First we define the differential constraint which will be related to the relevant BSSN variables below. The ADM equations, providing evolution equations for the physical 3-metric γ ij , the extrinsic curvature K ij and f i , are readjusted to [131,134] ∂ t γ ij =[ADM] , (2.60a)

Conformal decomposition and densitization
The generalized BSSN system [134] employs conformal variables as dynamical variables where χ andγ ij are the conformal factor and metric, γ is the determinant of the physical 3-metric,K andÃ ij are the densitized trace and tracefree part of the extrinsic curvature andΓ i is the conformal connection function. The densitization parameter n K parametrizes the change of conformal variables. Note, that by construction the determinant of the conformal metricγ = 1. The key difference between the original BSSN evolution scheme and the generalized version is the densitization of the curvature variables. Notice, that we recover the BSSN system if the parameter n K vanishes. The definition of the conformal variables implies a set of additional algebraic constraints representing the requirements thatÃ ij remains tracefree and the determinant of the conformal metricγ = 1. The definition constraint, Eq. (2.59), becomes The evolution equations for the dynamical variables (2.61) are obtained by taking the time derivative of their definitions, substituting the evolution equations (2.60) and rewriting everything in terms of the conformal quantities. This procedure yields the evolution where [. . .] tf denotes the tracefree part with respect to the physical metric. The Ricci tensor is given by whereD i is the covariant derivative with respect to the conformal metricγ ij . The Christoffel symbols Γ i jk andΓ i jk associated with the physical and conformal metric, respectively, are related via The second derivative of the lapse in Eqs. (2.64) is given by The physical constraints, Eqs. (2.24) and (2.25), become (2.68b)

Gauge conditions
The system of evolution equations (2.64) is closed by the choice of gauge conditions for the lapse function α and shift vector β i . The BSSN system has been proven to be a wellposed, numerically robust formulation of the initial value problem when combined with the moving puncture approach. The moving puncture approach employs the 1 + logslicing condition, Eq. (2.50), and the Γ-driver shift condition, in the version of Eqs. (2.57) or Eq. (2.58). In the generalized BSSN system we consider the densitized lapse function Q, related to the lapse function α via as dynamical variable [134]. n Q is the densitization parameter and as such parametrizes the particular choice of the dynamical variable. The 1 + log-slicing condition (2.50) becomes One of the key goals in performing binary BH evolutions in astrophysical and other applications is the extraction of gravitational radiation. In order to gain information about gravitational waves two distinct methods have been developed: (i) The approach first introduced by Regge, Wheeler and Zerilli [150,151] considers perturbations of the Schwarzschild metric. From these metric perturbations one can compute a gauge invariant master function which encodes the information about gravitational radiation.
(ii) The second approach is based on the tetrad formalism introduced by Newman and Penrose [152]. In particular one computes the (coordinate) gauge independent Weyl scalars 1 , which are contractions of the Weyl tensor with the principal null directions of the tetrad. In particular, one can show that in a specific choice of this tetrad the Weyl scalars Ψ 4 and Ψ 0 carry information about out-and ingoing gravitational radiation [14,106,[153][154][155] . Note also, that Teukolsky's analysis of the stability of the Kerr BH [156] employs the Newman-Penrose formalism.
In this section we restrict ourselves to the summary of the Newman-Penrose approach [152]. We will study an extension of the Regge-Wheeler-Zerilli formalism to higher dimensional spacetimes in detail in Sec. 4.6 and provide a connection to the four dimensional case in this context. The relation between both methods has been investigated, e.g., in Refs. [3,[157][158][159].
Nowadays, the most common approach to extract information about gravitational radiation in numerical simulations of 4-dimensional BH binaries is the computation of the Weyl scalar Ψ 4 . It has the asymptotic property of being equal to the outgoing radiation if the complex null-tetrad is chosen properly. In Chapter 6 we intend to study BH binaries surrounded by a mirror-like box, thus mimicking BH-AdS spacetimes. This toy model has the same local geometry as vacuum models and for a sufficiently large box size we expect the standard tetrad to yield a good approximation to the gravitational wave information. Because of the imposed reflecting boundary conditions we also deal with ingoing waves in our simulations, once they are reflected from the box boundary. Therefore, we will also consider the Weyl scalar Ψ 0 in order to account for the ingoing contribution.
For self-containedness, Appendix A reviews the electromagnetic decomposition of the Weyl tensor and in particular the construction of the relevant quantities for our study, Ψ 0 and Ψ 4 .
To be explicit, we define a spherical coordinate system centered at the center-of-mass of the binary with orthonormal basis (r,θ,φ). The coordinates are chosen such that the azimuthal axis is aligned with the orbital angular momentum and the binary orbits are in the direction of increasing azimuthal coordinate. Our definitions and notation are the same as in [160,161]. To define our complex null-tetrad, we use the timelike unit vector normal to a given hypersurfacen and the radial unit vectorr to define an ingoing (k) and outgoing null vector ( ) by We define the complex null vector m and its complex conjugate by In terms of this tetrad, we define Ψ 0 and Ψ 4 as where C αβγδ is the Weyl tensor. To relate Ψ 0 and Ψ 4 to the amplitudes of the gravitational waves, we note that in the transverse-traceless (TT) gauge, assuming the functional form f (t ± r) for the ingoing or outgoing waves, we have A 'dot' denotes derivative with respect to the argument (t ± r). Following standard conventions [106,153,154], we take the h + and h × polarisations of the gravitational waves to be given byḧ Then, we find that in vacuum regions of the spacetime, for outgoing wavesḧ while for ingoing wavesḧ +,× =ḧ +,× (t + r) The fact that Ψ 4 (Ψ 0 ) are zero for ingoing (outgoing) waves is consistent with the leading order in perturbation theory. The correct expression and in particular its dependence on the radial coordinate is given by the solution of the Teukolsky master equation at large distances from the source, which states that for outgoing waves [156] Ψ 0 ≈ e iω(t−r) while for ingoing waves We decompose the resulting NP scalar NP scalarss Ψ 4 (Ψ 0 ) into modes by projection onto spherical harmonics of spin-weight s = −2 (s = 2) according to spin-weighted spherical harmonics [162,163] defined as where d l ms is the Wigner d-function and where C 1 = max(0, m − s) and C 2 = min(l + m, l − s). Here M is the ADM mass of the system, computed from the initial data and assuming this is an asymptotically flat spacetime.
In the numerical code, the null-tetrad is constructed from a Cartesian orthonormal triad (u, v, w) and the timelike vectorn is orthonormal to t = constant hypersurfaces. The spacetime is evolved in time t using Cartesian coordinates x, y, z. In practice, we compute the Newman-Penrose scalars Ψ 0 and Ψ 4 using the electromagnetic decomposition of the Weyl tensor according to Eqs. (A.7) on the entire Cartesian grid. Then, they are interpolated onto coordinate spheres of different extraction radii r ex with a uniform distribution of points in (θ, φ). All the waveform related data from the simulations presented in the course of this work are taken from such samplings of Ψ 0 (t, r = r ex , θ, φ) and Ψ 4 (t, r = r ex , θ, φ). A more detailed description is given in Appendix A.
Given the Newman-Penrose scalar Ψ 4 , we can compute the radiated energy, linear and angular momentum from the radiation content [106,164], and references therein: where i = (− sin θ cos φ, − sin θ sin φ, − cos θ) . The definitions above are based on time integrals which start in the infinite past (at retarded time t = −∞), and thus capture the complete gravitational wave signal. Starting the time integrations at t = −∞ corresponds to the limit of infinite extraction radius on the initial time slice -the slice would then extend all the way to spatial infinity, no part of the waveform would be lost, and it would take an infinite time for the waves to reach the extraction sphere. This situation cannot be handled with the current numerical codes; we therefore work with finite extraction radii.
The mass and angular momentum of the final BH can be estimated from balance arguments. Given the parameters (P y i , d) in the Bowen-York initial data, we straightforwardly calculate the total initial angular momentum as if the initial spin of each BH is zero. Ansorg's TwoPunctures [165] initial data solver directly provides the total ADM mass M of the system and we obtain radiated energy and angular momentum E rad and J rad from the gravitational wave signal. In case of a merger, this gives us the final angular momentum and mass of the BH The dimensionless spin parameter of the final BH follows directly from (2.85)

Apparent horizon properties
A further method to characterize the evolution of a BH spacetime is provided by the properties of the apparent horizon (AH) of the final BH itself. Since this relies only on local quantities, it does not depend upon the spacetime being asymptotically flat. In order to monitor the mass and spin of the final BH, we use Thornburg's Apparent Horizon Finder AHFinderDirect [166,167] in different ways, which also allow us to obtain uncertainty estimates: For comparison we also compute the spin of the final BH from the two following estimates: 2. We measure the ratio C r (j) = C p /C e of polar to equatorial circumference of the final BH [169]. If we assume the final object to be a Kerr BH, this ratio is C r = 1 − β 2 E(β 2 ), where β 2 ≡ j 2 M/(2r + ), E(β 2 ) is a complete elliptic integral and r + /M = 1 + 1 − j 2 . This expression can be inverted to find the dimensionless spin parameter, j Cr , of the final BH.
3. The equatorial circumference of a Kerr BH is C e = 4πM . Therefore 2πA AH /C 2 e = 1 + 1 − j 2 AH , where A AH is the area of the apparent horizon. Thus, the AH area and the equatorial circumference can be used to estimate the spin of the final BH from [170] (2.88) 3. Collisions of unequal mass black holes and the point particle limit

Introduction
Black holes play a key role in a variety of processes in astrophysics, gravitational wave physics and high-energy physics. Following the 2005 breakthroughs [10][11][12], NR has been an essential tool in the modeling of BH binaries in the strong-field regime. At the same time it has become clear that detailed studies of BH systems often involve a close interplay between fully non-linear numerical simulations and (semi)-analytic approximation techniques of various types. For example, the generation of gravitational wave (GW) template banks for use in the analysis of observational data from operational laser interferometric GW detectors such as LIGO, VIRGO, GEO600 and KAGRA or futureplanned observatories such as the Einstein Telescope requires the combination of NR with post-Newtonian or other techniques; see Refs. [14,16,33,[171][172][173][174][175][176] and references therein. Post-Newtonian studies have also played an important role in the guidance of numerical investigations of the BH recoil, most notably in the discovery of the socalled superkicks and their possible suppression due to spin alignment [177][178][179][180][181][182][183]. In the context of high-energy collisions of BHs, linearization tools such as the zero-frequency limit or point particle calculations provide valuable insight into the scattering threshold and GW emission of BH collisions in four and higher-dimensional spacetimes [184]. Of particular relevance for future spaceborne missions in the spirit of LISA, are extrememass-ratio inspirals (EMRIs). They represent a particularly difficult challenge to NR and their modeling relies heavily on perturbative methods and self-force calculations; see Refs. [185][186][187][188][189][190] and references therein. With the above as motivation, it is vital to obtain a detailed understanding of the range of validity of the various types of approximation methods. At the same time, these methods provide valuable tools to calibrate the accuracy of numerically generated solutions to the Einstein equations. In this chapter, based on [1], we provide such a study for the case of a classical calculation in GR, the head-on infall of a point-particle (PP) into a BH [191].
In recent years, NR has started probing the intermediate mass-ratio regime (i) by evolving the final orbits of (approximately) quasi-circular inspirals of BH binaries with mass-ratio q ≡ m 2 /m 1 = 1/10 [192,193]; (ii) by comparing numerical results with perturbative calculations employing the fully numerical BH trajectories for mass ratios up to q = 1/20 [194]; (iii) by evolving the last orbits of a quasi-circular inspiral with mass ratios up to q = 1/100 [195][196][197]. Recently, calculations of EMRIs with q = 1/10 6 in full GR have been presented [198] Nevertheless, in this chapter based on Ref. [1] I focus on head-on collision of BHs with q ≥ 1/100, for two reasons: (i) the lower computational cost due to the higher degree of spacetime symmetry and the absence of the lengthy inspiral phase and (ii) the availability of high-precision results in the PP limit.
In our study we will make extensive use of the calculation by Davis et al. [191] who model the collision of a small object of mass m with a BH of mass M m in the PP limit. In the original calculation the particle was falling from rest at infinity, and the total radiated energy was found to be E rad PP =0.0104 This setting has been generalized to arbitrary initial distance and boost, in which case initial data and consequent spurious radiation play a role [184,[199][200][201][202]. Fully numerical results for BH head-on collisions obtained in the equal and comparable mass regime have been compared with PP predictions and results obtained in the close-limit approximation [203] by Anninos et al [204,205]. These studies demonstrated agreement for the radiated energy and linear momentum bearing in mind the accuracies achievable at the time. The waveforms presented therein, however, exhibit a significant signal starting at t − R ex = 0; see for example Fig. 1 in [205]. This contribution most likely arises from spurious radiation inherent in the initial data due to the small initial separation and its impact on the quantitative conclusions is not entirely clear.
At the time when the present research was carried out no comparisons between PP calculations and fully numerical results for mass ratios in a truly perturbative regime have been available (but note Ref. [198]). By simulating BH binaries up to a mass ratio of q = 1/100 we fill this gap and identify those aspects of the PP predictions which describe BH dynamics well in general and which only hold in the extreme mass-ratio limit. From a different point of view, the agreement with the PP calculations represents an important validation of the fully numerical calculations in the regime of high-mass ratios. In this context we emphasize that we are able to accurately extract radiated GW energies of the order of 10 −6 M and linear momenta corresponding to recoil velocities of a few dozens of m/s from binary BH simulations. We note, however, that even smaller amounts of energy have been extracted from general relativistic simulations of stellar core collapse; see e. g. [159].

Numerical Setup and Analysis Tools
The numerical simulations of unequal-mass BH collisions starting from rest have been performed with the Lean code [161,206]. The Lean code is based on the Cactus computational toolkit [207,208] and uses the Carpet mesh refinement package [209,210], the apparent horizon finder AHFinderDirect [166,167] and the TwoPuncture initial data solver [165]. The 3 + 1 Einstein's equations are evolved using the BSSN [132,133] formulation, together with the moving puncture approach [11,12], discussed in Sec. 2.4. The gauge conditions are determined by the puncture gauge, i.e., the 1+log slicing and Γ driver shift condition [127] (cf. Sec. 2.3). The systems are set up using Brill-Lindquist initial data, summarized in Sec. 2.2. We have evolved BH binaries with mass ratios q ≡ m 2 /m 1 = 1, 1/2, 1/3, 1/4, 1/10 and 1/100, where m i is the bare mass parameter of the i-th BH.
We use the Newman-Penrose scalar Ψ 4 to measure gravitational radiation at extraction radii R ex , chosen in a range of 40 M to 90 M from the center of the collision. We decompose Ψ 4 into multipoles ψ lm using s = −2 spin-weighted spherical harmonics −2, −2 Y lm , according to Eq. (2.79). Due to the symmetry properties of the systems under consideration, the only non-vanishing multipoles all have m = 0 in a suitably chosen frame, and are purely real, corresponding to a single polarization state h + . In the equal-mass limit, the additional symmetry causes all multipoles with odd l to vanish identically. The energy spectrum and luminosity of the radiation are given by respectively, where a hat denotes the Fourier transform and ψ l0 is evaluated on a sphere at infinity.

Simulations and uncertainties
We have performed a series of simulations of head-on collisions with mass ratio ranging from q = 1 to q = 1/100 with initial coordinate separation d and proper horizonto-horizon separation L as given in Table 3.1. We describe the grid setup used for these simulations in terms of the number n rl of refinement levels, the radius R of the computational domain, the resolution H used in the wave extraction zone, typically the third refinement level counted from the outside, the radius r in units of the smaller hole's mass m 2 of the innermost refinement level centered on the individual punctures and the resolution h/m 2 of the innermost refinement level. The values for these parameters are summarized for all mass ratios in Table 3.2. Because of the vastly different length scales in case of the small mass ratios q = 1/10 and q = (1/100) we include, respectively, 2 or 5 additional refinement levels around the smaller BH. Thus we are able to evolve the system at feasible computational costs. Our results are affected by three main sources of uncertainties: finite extraction radius, discretization and, for small initial separations of the binary, spurious initial radiation. We reduce the error arising from finite extraction radius by measuring the waveform components at several radii, and fitting them to an expression of the form ψ lm (r, t) = ψ   lm (t) is the quantity reported throughout this work and used to calculate related quantities, such as the radiated energy. The uncertainty in this extrapolated value is estimated by performing a second fit including also a quadratic term ψ (2) lm /r 2 , and taking the difference between the first-   : Grid setup used for the different mass ratios q. The number of refinement levels is given by n rl , R is the radius of the computational domain, H the resolution in the wave extraction zone, r the radius of the innermost refinement box around the individual punctures and h the resolution used on that level. The additional low and high resolution for q = 1/4 and q = 1/100 have been used for the convergence studies. and second-order fits. The resulting uncertainty increases as we decrease the mass ratio q and is 1 − 4 % for the total radiated energy and the l = 2 waveform and energy, and 3 − 5 % for the subdominant multipoles and the radiated linear momentum. In order to estimate the discretization error of our simulations, we have performed a convergence analysis for models (q = 1/4, L = 16.53 M ) and (q = 1/100, L = 9.58 M ) using the three resolutions listed in Table 3.2. The resulting convergence plots for the l = 2 multipole of the wave signal is shown in Fig. 3.1 and demonstrates convergence between second and fourth order. With regard to the analysis below, we note in particular that the q = 1/100 case exhibits second order convergence in the plunge-merger signal around t − R ex ≈ 40 M but is close to fourth-order convergence for the remainder of the waveform. Bearing in mind that the plunge-merger transition represents the most dynamic part of the evolution and that the second-order ingredients in the code are associated with the prolongation of grid functions at the refinement boundaries in time, this observation is compatible with the numerical discretization. We observe similar convergence properties for the l = 3 multipole, but overall convergence close to fourth-order for the radiated energy and linear momentum, presumably because the accumulated errors are dominated by the fourth-order contributions observed for most of the signal. The resulting numerical uncertainties for q = 1/100 are about 10 % in the waveform for the plunge-merger transition and 5 % for the remainder of the signal as well as 6 % for the radiated energy and 8 % for the linear momentum emitted in gravitational waves. We note that in both cases the discretization error leads to an overestimate of the radiated quantities. For q = 1/4 we observe significantly smaller uncertainties in the range of 2 % for all quantities.
Finally, we comment on the unphysical gravitational radiation inherent in the conformally flat puncture initial data. In order to extract physically meaningful information, one has to separate the spurious radiation from the radiation generated by the collision itself. This is done by "waiting" for the spurious radiation to radiate off the computational domain, and then discarding the early, contaminated part of the wave signal. For small values of the initial separation, however, the binary will merge before the spurious radiation has had enough time to leave the system, and physical and unphysical contributions to the wave signal partially overlap and cannot be cleanly distinguished. For our set of simulations, this problem arises only in the case q = 1/100, L = 9.58 M , where it introduces an additional error of about 2 % to the radiated energy and momentum.

Results
All collisions summarized in Table 3.1 result in the formation of a single BH plus gravitational radiation, i. e. there is no indication of violation of the cosmic censorship conjecture. The final BH is born distorted, and eventually rings down to a Schwarzschild solution via emission of quasinormal modes [211].
We illustrate the l = 2 and l = 3 wave signal in Fig. 3.2 for the l = 2 and l = 3 multipoles obtained for the mass ratios q = 1/4 (top), q = 1/10 (center) and q = 1/100 (bottom). In each panel the solid (black) curves represent the PP prediction for infall from infinity whereas the dotted (red) and dash-dotted (blue) curves show the numerical results for different values of the finite initial separation. To leading order,  [199].
the gravitational radiation output of BH collisions scales with the square of the reduced mass µ ≡ M η of the system, where η = q/(q + 1) 2 is the dimensionless, symmetric mass ratio [191]. For comparison of the numerical results with PP predictions, we therefore rescale the former by the corresponding powers of η, quadratic for energy and linear for the waveforms in Fig. 3.2.
The waveforms show interesting features. For small initial separations, the early part of the waveform is contaminated by "spurious" radiation; cf. the dotted (red) curve in the top and bottom panels of Fig. 3.2. As the initial separation increases, however, this problem disappears, because the longer infall duration of the binary provides sufficient time for the unphysical radiation to propagate off the grid; cf. the dash-dotted (blue) curves. A closer inspection of the q = 1/100 case yields excellent agreement between the numerical and PP predictions except for the plunge-merger transition around t ≈ 0 in the figure. From the discussion in Sec. 3.3, however, we recall that the discretization error is particularly large in this regime. In fact, for the q = 1/100 model studied in Sec. 3.3, a second-order Richardson extrapolation predicts about a 10 % reduction in the amplitude around the first strong maximum in the l = 2 waveform which is very close in magnitude and sign to the deviation of the numerical from the PP result. As demonstrated by the upper central panel in Fig. 3.2, we find equally good agreement of the numerical l = 2 multipole with PP predictions for the less extreme mass ratio q = 1/10 and only a small deviation for the larger mass ratio q = 1/4 (upper top panel in Fig. 3.2). Our findings thus confirm over a wide range of mass ratios the observation by Ref. [204], that there is a weak dependence of the re-scaled waveforms on the mass ratio. The l = 3 mode, on the other hand, is a good discriminator between high-and low-mass ratios. This behavior was qualitatively expected, as higher multipoles are suppressed in the equal-mass case; by symmetry the l = 3 mode is absent when the masses are equal. It is interesting, however, that even for what one might call a small mass ratio, q = 1/10, higher multipoles are still visibly suppressed.
The total amount of energy radiated in gravitational waves during the collision depends on the initial separation of the holes. As discussed in Anninos et al. [204], two effects contribute to increasing the GW energy at larger initial separations; (i) there is more time to radiate GWs during the infall and (ii) the infalling velocity is larger. In practice, the second effect is found to be dominant. Anninos et al. have accounted for both contributions by defining One can write the corrections to the radiation emission With the above as motivation, we have fitted our results to a 1/L dependence, of the form with E rad ∞ the radiated energy for infinite initial separation. The results are summarized in Table 3.3. We remind the reader that L stands for proper initial separation between the holes. We also note that the results in Table 3.3 are normalized by η 2 . For comparison, we also show in the last entry of the table the results obtained in the PP limit, within a linearized calculation. This study was done by Lousto and Price [199] using the same type of initial data; we have used their Table I to obtain the behavior shown in Table 3.3 above. We note that already for q = 1/10 and q = 1/100 our results are in good agreement with PP calculations. We remind the reader, however, that in the q = 1/10 case there is a larger deviation in the l = 3 modes.
With the extrapolation above one gets an estimate for the total radiation of two BHs merging from infinite initial separation. A best fit of this number as function of mass ratio yields E rad ∞ M η 2 =0.0110 − 0.0088η (3.4) In the PP limit, when η → 0, this agrees with the classical PP calculation, Eq. (3.1) to within 6%, i.e., within the numerical uncertainties. Overall, the results in Table 3.1 demonstrate that we are able to accurately measure amounts of order E rad ∼ 10 −6 M in these fully nonlinear evolutions. The amount of spurious radiation in the initial data is also consistent with predictions from linearized gravity. Lousto and Price performed a detailed analysis of the amount of spurious radiation in the infall of PPs into massive BHs, using the same type of initial data [199]. Using their Table I for L > 11, we find that the amount of spurious radiation varies with L according to E rad /(M η 2 ) ∼ 0.15(L/M ) −2.5 . For q = 1/100, for instance, we obtain E rad /(M η 2 ) = 0.26(L/M ) −2.55 . Thus, we find good agreement in the decay power (roughly −2.5) and also in the proportionality coefficient.
If two BHs with different masses collide head-on, the remnant BH will recoil with respect to the center-of-mass frame, due to the emission of energy and momentum carried by gravitational waves. Based on PN tools, we have fit our results to [212] v where v f ∞ is a normalized recoil velocity for infinite initial separation. The normalized recoil velocity v f ∞ is shown in Table 3.3. The point particle limit was considered in Ref. [213], who obtained v f ∞ = 263km/s 1 . We note this is not a trivial agreement: unlike energy calculations, momentum involves interference with higher (typically highly suppressed) multipoles. Overall, our results agree well in the limit of small mass-ratios with the point particle limit. It is interesting to note in this context that for both, radiated energy and linear momentum, the numerical results exceed those obtained from the point particle limit by about 6 %. This value agrees in sign and magnitude with the discretization error obtained for the q = 1/100 simulation in Sec. 3.3. We therefore consider the discretization error the dominant source of the remaining discrepancies.

Conclusions
The simulation of dynamical, interacting BHs has a tremendous potential to provide answers to some of the most fundamental questions in physics. Recent developments in experimental and theoretical physics make this a pressing issue. We refer, in particular, to the prominent role of BHs in the gauge-gravity duality, in TeV-scale gravity or even on their own as solutions of the field equations [2]. Recent work along these lines includes the successful simulation and understanding of the collision of two BHs at close to the speed of light in four-dimensional spacetime [214][215][216][217], the low energy collisions in higher spacetime dimensions [2][3][4], BH scattering in five dimensions [218], stability studies in higher dimensions [219][220][221] and BH evolutions in non asymptotically flat spacetimes [8,76,222].
We have shown here that NR is capable of simulating dynamical BHs close to the regime of validity of linear calculations, and to make contact with (semi-)analytic approximation techniques. For this purpose we have evolved head-on collisions of nonspinning BH binaries over a range of mass ratios from q = 1 to q = 1/100. We obtain radiated energies decreasing from about 5.5 × 10 −4 for q = 1 to 10 −6 for q = 1/100. The recoil reaches a maximum of about 4 km/s near q = 3 and decreases towards 26 m/s for q = 1/100. In the limit of small mass ratios and extrapolating our results to infinite initial separation, we find the numerical values for radiated energy and linear momentum to be ≈ 6 % larger than the point-particle predictions. This discrepancy agrees rather well in sign and magnitude with the discretization error obtained from a convergence study of our q = 1/100 simulations. It thus appears likely that a significant part of the remaining differences can be attributed to the discretization error which mirrors the computational demands of numerical BH binary simulations with such small mass ratios.
With regard to the waveforms, the most remarkable result is the suppression of odd l multipoles. While we observe good agreement between numerical and point-particle results for the l = 2 mode, already for q = 1/10, the numerically calculated l = 3 multipole is visibly suppressed for this case and only agrees well with the PP limit for q = 1/100.
Overall, the good agreement for waveforms and radiated energy and momenta for the case q = 1/100 demonstrates that numerical techniques are capable of bridging the gap between linear analysis and the fully non-linear regime of general relativity.

Black hole collisions in higher
dimensional spacetimes -Framework
Of particular interest are the dynamics of BHs in higher dimensional spacetimes -a quite recent research branch in NR. Perturbative methods are not sufficient to tackle most outstanding issues in the strong-field regime, and the employment of numerical relativity methods is inevitable. Over the past few years, our group -among few others worldwide -has started to explore highly non-linear, dynamical scenarios, such as (headon) collisions of BHs in D ≥ 5 dimensional spacetimes and calculated the associated emission of gravitational radiation. These calculations provide improved estimates for the amount of energy and angular momentum emitted in form of gravitational waves and allows to compute (by energy balance arguments) the mass and spin of the newly created BH. In the context of TeV gravity scenarios this information is of utmost importance for BH event generators, such as Blackmax, Catfish, Charybdis2 or TrueNoir [50,[237][238][239][240]. Conceptually, such a NR code is based either on an extension of the GHG formulation or on the explicit splitting of space and time, thus evolving the dynamical system in time (see e.g. [16][17][18][19][20] for reviews on recent progress in higher dimensional NR). However, a full-blown (D − 1) + 1-approach, analogous to the 3 + 1-approach presented in Chapter 2, is computationally too expensive for currently available computational resources. Therefore, the problem of dynamically evolving higher dimensional spacetimes has been simplified to axissymmetric or effectively 3 + 1-dimensional setups [2, 3, 218-221, 235, 241-244].
In this chapter I will review the formalism and framework (developed by our group) adopted to perform numerical simulations of BH collisions in higher dimensional spacetimes. For a summary of the numerical results I refer the reader to Chapter 5. The presented summary is based on the publications [2][3][4][5] as well as on as yet unpublished work in progress [6].
The starting point for our approach are D-dimensional spacetimes with an isometry group. If this isometry group is sufficiently large, it permits the dimensional reduction of the problem to a 4-dimensional model. Specifically, I obtain 4-dimensional GR coupled, in the most general case, to scalar and gauge fields. The additional (D−4) spatial dimensions manifest themselves only in these additional fields emerging from the dimensional reduction. We focus on two types of models in particular, which are generalisations of axial symmetry to higher dimensional spacetimes: (i) D ≥ 5 dimensional vacuum spacetimes with an SO(D −2) isometry group, and (ii) D ≥ 6 dimensional vacuum spacetimes with an SO(D − 3) isometry group. In order to end up with an effectively 4-dimensional model I employ only part of this symmetry: Specifically, I perform a dimensional reduction by isometry on a (D −4)-sphere that has an SO(D −3) ⊂ SO(D −2) isometry group which allows us to investigate head-on collisions of BHs in D ≥ 5 as well as BH collisions with impact parameter and with spinning BHs in D ≥ 6, as long as all the dynamics are restricted to a single plane. These classes include the most interesting physical configurations relevant to accelerator and cosmic ray physics in the context of TeV-scale gravity, and to fundamental questions concerning the properties of higher-dimensional black objects such as their stability and phase diagrams.
Subsequently, I cast the effective 4-dimensional model into a time evolution problem by employing the 3 + 1-split, discussed in Sec. 2.1. The resulting ADM-like equations are then reformulated in the generalized BSSN evolution scheme [132][133][134], as described in Sec. 2.4. The system is completed by specifying the gauge using the moving puncture approach [11,12] adapted to our setup.
In the following sections I will present the framework for NR simulations of higher dimensional spacetimes, including the dimensional reduction to an effectively 4-dimensional model in Sec. 4.2, the specific formulation of the time evolution problem in Sec. 4.3, the setup of initial data in Sec. 4.5 and the technique used to extract gravitational radiation in Sec. 4.6. Additionally, I show that the generalized BSSN formulation of the evolution equations of the effectively 3+1-dimensional GR model coupled to a scalar field together with the modified moving puncture gauge exhibits a strongly hyperbolic PDE system in Sec. 4.4. The presented formalism has been implemented in the HD-Lean code [2,3] which is an extension of Sperhake's original Lean code [161].

4 + (D − 4)-form and dimensional reduction
The key ingredient of our NR framework for higher dimensional spacetimes is the reduction from D-dimensional vacuum GR to an effectively 4-dimensional model. The isometry group of D-dimensional Minkowski spacetime is ISO(1, D − 1) while BH solutions of GR generically break this symmetry into a subgroup. For instance, the isometry group of a D-dimensional Schwarzschild-Tangherlini BH [245] is SO(D − 1) × R, whereas the isometry group for a head-on collision of two non-rotating BHs is SO(D −2). In fact, neither the time direction nor the direction of the collision correspond to symmetries, but a rotation of the remaining (D − 2) spatial directions leaves the spacetime invariant. The spacetime can then be considered as the semi-direct product of a 3-dimensional spacetime (3) M with the sphere S D−3 = SO(D − 2)/SO(D − 3) and we take advantage of this symmetry to reduce the spacetime dimensionality. This can be accomplished by writing the D-dimensional vacuum Einstein's equations where M, N = 0, ..., (D − 1) are spacetime indices, in a coordinate system which makes the symmetry manifest and decomposing Einstein's equations with respect to this symmetry, thus providing a lower dimensional interpretation (in the spirit of Kaluza-Klein reduction). Note, however, that we are not performing a compactification but a dimensional reduction by isometry, as first proposed by Geroch [246]. Then, the extra dimensions appear as source terms of the lower dimensional Einstein's equations, defined on the lower dimensional manifold. Specifically, we use the SO(D − 3) symmetry, thus reducing the spacetime on a (D − 4)-sphere which yields a 4-dimensional manifold. The original proposal by Geroch [246] considered the symmetry space SO(2). A 5-dimensional extension with the same symmetry space has been derived in [247]. A generalisation to coset manifolds (like the sphere S n ) was given by Cho in [248,249]. We now describe in detail the reduction from a D-to a 4-dimensional spacetime. In order to highlight the particular classes of BH binaries we are able to study with this framework, it is convenient to begin this discussion with the isometry group of the S D−3 sphere, i.e., with the 3 + (D − 3)-split. The general D-dimensional spacetime metric can be written in the form where ābc are the structure constants of SO(D − 2). Because the fibre has the minimal dimension necessary to accommodate (D −3)(D −2)/2 independent Killing vector fields, we assume, without loss of generality, that the Killing vector fields have components exclusively along the fibre: ξā = ξīā∂ī. Furthermore, we normalise the Killing vectors so that they only depend on the coordinates of the fibre, i.e. ∂μξīā = 0. Then, Eq. (4.3) yields the conditions These expressions can be interpreted either as Lie derivatives of rank-2 tensors defined on the D-dimensional spacetime, or as Lie derivatives of a rank-2 tensor, a vector and a scalar, which are defined on S D−3 . Condition (4.5a) implies because Ωīj admits the maximal number of Killing vector fields and thus must be the metric on a maximally symmetric space at each xμ. Due to (4.4) this space must be the denotes the metric on an S D−3 sphere with unit radius. Eq. (4.5c) implies because the Killing vector fields ξā act transitively on the fibre and therefore the base space metric must be independent of the fibre coordinates. Because Eq. (4.5b) is equivalent to and there exist no non-trivial vector fields on S D−3 for D ≥ 5 that commute with all Killing vector fields on the sphere. We remark that (4.8) corresponds to the statement, expressed in [248] in group theoretical language, that the gauge group for a theory reduced on a coset space G/H is the normaliser of H in G; in the case of a sphere, where G = SO(D −2) and H = SO(D −3), the normaliser vanishes and then there are no "gauge vectors", i.e., no non-vanishing metric components gμī. If the normaliser of H in G is non-vanishing, such metric components appear, and with the dimensional reduction they yield gauge fields in the reduced theory. For example, in the case of head-on collision in D = 4 the isometry space is SO(2) and scalar as well as vector fields emerge in the reduced theory [246]. In D > 4, the isometry space is SO(D − 2)/SO(D − 3) and the dimensional reduction only yields a single scalar field. In the remainder of this work we focus on this subclass of spacetimes. Because we intend to employ and modify an existing 3 + 1-NR code, we now perform the 4 + (D − 4)-split of the D-dimensional spacetime. The metric on a unit sphere S D−3 can always be written in terms of the line element dΩ D−4 on a unit sphere S D−4 where θ ∈ [0, π] is a polar-like coordinate. Now we introduce 4-dimensional coordinates, x µ = (xμ, θ), µ = 0, . . . , 3, and define a 4-dimensional metric  as well as a new (conformal) factor The last equality holds for transformations into coordinates adapted to axial symmetry. Then, the most general D ≥ 5-dimensional metric compatible with SO(D − 2) isometry is The geometry (4.13) has only a manifest SO(D −3) symmetry. Note, that the additional scalar field λ ∼ y 2 (in coordinates adapted to the axial symmetry) arouses problems for the numerical implementation when y = 0. Therefore we will introduce regular variables in Sec. 4.3.2. The treatment of potentially pathogolical terms is presented in Appendix C.1.
Performing the dimensional reduction of the D-dimensional vacuum Einstein's equations (4.1) on a (D − 4)-sphere and considering the discussed symmetries leads to the 3 + 1-dimensional Einstein's equations coupled to the scalar field λ. If we consider SO(D − 2) as the full isometry group, the emergent scalar field is completely determined by the 3 + 1-dimensional geometry via Eq. (4.12). Instead, considering the full isometry group SO(D − 3) allows us to explore BH collisions with impact parameter and with spin as long as it is restricted to the same 2-plane. We depict both considered scenarios in Fig. 4.1. Taking metric (4.13), which describes general spacetimes with SO(D − 3) isometry in D ≥ 6 as follows from the discussion of (4.9), the D-dimensional vacuum Einstein's equations (4.1) reduce to a 4-dimensional theory of gravity coupled to a scalar field λ(x µ ). We remark that in this theory λ and g µν are viewed as independent degrees of freedom. The D-dimensional Einstein-Hilbert action reduces to where the D-dimensional Newton's constant G D is related to the 4-dimensional one G 4 by the area of the unit D − 4-dimensional sphere The resulting equations of motion are where ∇ µ is the covariant derivative with respect to the 4-metric g µν . The energymomentum tensor is given by

Formulation as time evolution problem 4.3.1. ADM formulation
As we have seen in the previous section, the dimensional reduction on a (D −4)-sphere of D-dimensional vacuum GR results in the 4-dimensional Einstein's equations coupled to a scalar field λ, where the EoMs are given by Eqs. (4.15). We wish to evolve this system numerically and therefore rewrite Eqs. (4.15) as time evolution problem by adopting the 3 + 1-decomposition [103,111], discussed in detail in Sec. 2.1. In order to close the system we introduce the conjugated momentum K λ related to the scalar field λ via analogous to the definition of the extrinsic curvature, Eq. (2.16). Then, the evolution equations for the 3-metric γ ij and scalar field λ are given by The evolution equation for K λ is obtained from the EoM of the scalar field, Eq. (4.15a). The 3 + 1-split yields where D i is the covariant derivative associated to the 3-metric γ ij . Next, we consider the various projections of the Einstein's equations (4.15b), as outlined in Sec. 2.1, with the energy-momentum tensor specified in Eq. (4.16). In particular, the projection of Eq. (4.15b) onto the spatial hypersurface Σ t provides an evolution equation for the extrinsic curvature K ij where R ij is the Ricci tensor associated with the 3-metric γ ij . The Lie derivatives of the dynamical variables along the shift vector β i are given by The physical constraints are obtained in an analogous manner as in the "standard" 4-dimensional case, see Eqs. (2.24) and (2.25). In particular the full contraction of Eq. (4.15b) with the vector n µ normal to the hypersurface Σ t yields the Hamiltonian constraint The mixed projection of Eq. (4.15b) once onto the hypersurface and once perpendicular to it results in the momentum constraint In our framework the complete set of evolution equations for the 3-metric γ ij , the scalar field λ, the extrinsic curvature K ij and the scalar curvature K λ in ADM form are given by Eqs.

Generalized BSSN formulation
In order to obtain a hyperbolic formulation of the ADM evolution Eqs. (4.18), (4.19) and (4.20), we reformulate them in terms of the generalized BSSN system [132][133][134] (cf. Sec. 2.4). Employing the approach discussed in Sec. 2.4, we alter the characteristics of the PDE system by adding the definition-differential constraint G i , given in Eq. (2.63), and the Hamiltonian and momentum constraint, Eqs. (4.22) and (4.23). Then the system of evolution equations is modified according to where [ADMHD] denotes the source terms of the ADM-like Eqs. (4.18), (4.19) and (4.20). Next, we perform the conformal decomposition of the dynamical variables in order to evolve quantities that are regular (everywhere in space). Therefore we consider the conformal variables where the densitization constants n K , n Q parametrize the change of variables 1 . The definition of the conformal variables introduces additional algebraic constraints representing the requirements thatÃ ij remains tracefree and the determinant of the conformal metricγ = 1. Under this change of variables the evolution equations (4.24) become where [gBSSN4D] denotes the 4-dimensional vacuum evolution Eqs. (2.64) andD i denotes the covariant derivative with respect to the conformal metricγ ij . The coupling terms SK, SÃ ij and SΓ j are given by The adopted treatment and the regularization of terms ∼ 1/y or ∼ 1/y 2 , which are potentially troublesome when they reach y = 0 throughout a numerical simulation, are summarized in Appendix C.1.
In order to close the system of evolution equations we have to specify the gauge for the lapse function α and shift vector β i . In an attempt to find gauge conditions well-suited to the numerical evolution of higher dimensional black holes we modify the standard moving puncture gauge [11,12] by terms that couple the lapse and shift to the scalar field and its momenta. For the lapse we choose a modification of the Bona-Massó [126] condition In applications we always choose the 1+log-slicing with µ L = 2/α and, typically, µ λ = 1. Then, the 1 + log-slicing condition translates into for the densitized lapse Q considered as evolved variable in the generalized BSSN scheme.
For the shift we employ the modified Γ-driver condition The free parameters in the gauge conditions that affect the principal part of the system are the scalar functions (µ L , µ λ , ξ Γ , ξ λ ). The damping parameter η β is taken to be a constant in our applications.

4.4.
Hyperbolicity of the evolution system 4.4.1. Well-posedness of the initial value problem and strong hyperbolicity A system of partial differential equations is said to be well-posed if it admits a unique solution that depends continuously, in some appropriate norm, on given data. In physical applications well-posedness is a fundamental requirement both for a system to have predictive power, and for its treatment by numerical methods. Consider the linear, constant coefficient, first order in time, second order in space evolution system where S v , S w denote lower order derivatives. We denote s i as a unit spatial vector and we use the subscript s to denote contractions of the partial derivative with s i , A linear combination u of components of ∂ s v and w, is called a characteristic variable with speed λ in the s i direction if it satisfies ∂ t u = λ∂ s u + transverse derivatives. (4.33) The evolution system (4.32) is called strongly hyperbolic if it admits a complete set of characteristic variables that depend continuously on s i . The evolution system (4.32) admits a well-posed initial value problem if and only if it is strongly hyperbolic. These definitions are extended to non-linear systems with variable coefficients by linearizing them and working in the frozen coefficient approximation.

Hyperbolicity of the dimensionally reduced BSSN system
Strong hyperbolicity for quasi-linear and variable coefficient problems is defined by linearizing around a given solution and working in the frozen coefficient approximation. We perform a 2 + 1 split in space against the unit spatial vector s i by defining the standard projection operator Under this split a convenient choice for the scalar variables is Note thatγ qq andK qq are the natural choice of variables in the scalar sector because they correspond to the quantities that would be obtained by a (D−2)+1 split against s i in the higher dimensional BSSN formulation without the dimensional reduction by symmetry.
In this section we use indices A, B, ... to denote projected quantities, so the vectors are written and the tensors are Up to derivatives transverse to s i the system decomposes into scalar, vector and tensor blocks. We introduce the time derivative ∂ 0 of the scalar block for the metric can be read off from Eqs. (4.24).
and for the scalar field For the gauge conditions we consider a modification of the Bona-Massó condition (4.29) coupled to a parametrized modification of the standard 3 + 1 Γ-driver condition (4.31) 2 .
where we have re-defined the gauge parameter employed in Eqs. (4.29) and (4.31) tō We recover the 4-dimensional puncture gauge and the standard 3 + 1 principal symbol of BSSN whenμ λ =μ ζ = 0. Then, the fields (ln λ) andK λ satisfy a decoupled wave equation in the principal symbol. The vector part of the system is Finally, the tensor part of the system is given by Hyperbolicity of the scalar block: The scalar sector generically admits a complete set of characteristic variables, which are given by respectively, and where we employ the shorthandK = K ss +K qq . In the special case that generically distinct characteristic speeds coincide, the system fails to be strongly hyperbolic. In applications however, the speeds will typically intersect only on sets of measure zero. This scenario can be avoided by choosing the gauge parameters (µ L ,μ S ,μ λ ,μ ζ ) appropriately, such that either the characteristic speeds never coincide, or, if they do, the corresonding variables are decoupled in the principal symbol. For completeness we discuss the special cases in the following. In the scalar sector we recover the standard D = 4 restrictions [107,131], namely that the system is only weakly hyperbolic ifμ S = γ 1/3 ξ Γ = 3/4. This clash must be avoided by an appropriate choice ofμ S . For the new parameters we have the additional restrictions that in case µ L = 1, we must chooseμ λ = 0. In other words, if we choose harmonic slicing, corresponding to µ L = 1, it must be chosen in the full D dimensional spacetime to maintain strong hyperbolicity.
Hyperbolicity of the vector block The vector sector always admits a complete set of characteristic variables. They are given by with characteristic speeds (±1, ± √μ S ).
Hyperbolicity of the tensor block The tensor sector always admits a complete set of characteristic variables with speeds ±1. These characteristic variables correspond to the gravitational wave degrees of freedom in a linear, plane wave approximation to GR.

Initial Data
In this section we will discuss our method to construct initial data for D-dimensional BH spacetimes [5], which are then evolved by the effective 3 + 1-approach described in Secs. 4.2 and 4.3. In particular, we will focus on the preparation of initial configurations describing head-on collisions of two BHs starting either from rest, which is realized by modified Brill-Lindquist type initial data, or with non-zero boost, which is provided by generalized Bowen-York type initial data. We start by considering a (D −1)-dimensional spatial hypersurfaceΣ t with the induced (D − 1)-dimensional metricγ ab and extrinisic curvatureK ab embedded in a D-dimensional spacetime, where we use a, b = 1, ..., (D −1) throughout this section. The D-dimensional spacetime metric g M N is given by the line element Performing a (D − 1) + 1-splitting of the D-dimensional vacuum Einstein's equations (4.1) yields the constraintsH whereR andD a are the Ricci scalar and covariant derivative associated with the (D−1)dimensional spatial metricγ ab . Next, we conformally decompose the spatial metric and extrinsic curvatureγ Furthermore, we assume a conformally flat space, i.e.γ ab = η ab , and impose the maximal slicing conditionK = 0. With these choices, the (D − 1)-dimensional Hamiltonian and momentum constraints in vacuum, Eqs. (4.49), become [250,251] We employ a coordinate system which is adapted to the generalized axial symmetry SO(D − 2) in D = 5 dimensions and SO(D − 3) in D ≥ 6 dimensions, as discussed in Sec. 4.2. Therefore, we first consider the (flat) conformal metric in cylindrical coordinateŝ where dΩ D−4 is the metric on the (D − 4)-sphere. Note, that ϕ ∈ [0, π] is a polar rather than an azimuthal coordinate. Next we introduce "incomplete" Cartesian coordinates as x =ρ cos ϕ , y = ρ sin ϕ , (4.54) where −∞ < x < +∞ and 0 ≤ y < +∞. Then, the (D − 1)-dimensional metric is conformally decomposed according tō The ADM-like variables describing the dynamics of the D-dimensional spacetime with SO(D − 2) or SO(D − 3) isometry are the 3-metric γ ij and extrinisic curvature K ij , the scalar field λ and its conjugate momentum K λ . The goal of this section is the construction of initial data for these quantities. Note, that γ ij , K ij , λ and K λ are related to the (D − 1)-dimensional quantities bȳ where A, B = 4, ..., (D − 1) and h AB is the metric on the (D − 4)-sphere. In the following we summarize the construction of modified Brill-Lindquist and Bowen-York type initial data [5], setting up the initial configuration for head-on collisions of two BHs in higher dimensional spacetimes.

Brill-Lindquist initial data
First, we consider time symmetric initial data, i.e.,K ab = 0, which results in Brill-Lindquist type initial data. Then, the momentum constraint (4.52b) is trivially satisfied, whereas the Hamiltonian constraint (4.52a) reduces to the (D −1)-dimensional flat space For asymptotically flat spacetimes, the conformal factor satisfies the boundary condition lim r→∞ ψ = 1 and a solution ψ 0 to Eq. (4.58) is given by where r (i) is the coordinate location of the i th puncture. The mass parameter µ (i) is related to the horizon radius r S (i) and the ADM mass M (i) of the i th BH by where A D−2 is the area of the unit (D − 2)-sphere and we have set the D dimensional Newton constant to unity. By inserting our assumptions of conformal flatness, maximal slicing and time symmetry into Eqs. (4.56) the 3-dimensional variables describing Brill-Lindquist initial data become where ψ = ψ 0 is given by Eq. (4.59). These closed-form analytic data are the Ddimensional generalization of Brill-Lindquist data [119], adjusted to our approach of the (evolution) formalism, and describe a spacetime containing multiple non-spinning BHs at the moment of time symmetry, i.e., with vanishing linear momentum.

Bowen-York initial data
In order to numerically evolve BHs with non-zero boost, we consider the generalization of Bowen-York initial data to higher dimensional spacetimes. As shown by Yoshino et al. [251], we can write a solution of the momentum constraint (4.52b) describing a spacetime of arbitrary dimensionality D containing N BHs in the form where P a (i) corresponds to the ADM momentum of the i th BH in the limit of large separation from all other BHs. Next, we consider the Hamiltonian constraint (4.52a), which reduces to an elliptic equation for the conformal factor ψ whereÂ ab is prescribed by Eq. (4.62). A solution ψ 0 of the conformal factor is given by where r (i) and µ (i) are the (coordinate) position and mass parameter of the i th BH and u is a regular function. Inserting this ansatz into the Hamiltonian constraint (4.52a) yieldsˆ As in D = 4, the higher dimensional extension of Bowen-York extrinsic curvature data can also accommodate angular momentum of the BHs. Here, however, we shall focus on initial data for non-spinning, boosted BHs only. Without loss of generality, we can always choose coordinates such that the BHs are initially located on the z-axis at z 1 and z 2 and have momenta of equal magnitude in opposite directions P a (1) = −P a (2) . Inserting the momenta into Eq. (4.63) then provides the conformal traceless extrinsic cuvature and the differential equation (4.65) which is solved numerically for u. The class of symmetries covered by the formalism presented in Sec. 4.2 includes head-on and grazing collisions of non-spinning BHs with initial position x a , whereas a non-zero P y is not compatible with the assumed symmetries. However, the x-axis can always be oriented such that the collision takes place in the xz-plane. Our formalism therefore covers general grazing collisions of non-spinning BH binaries in D dimensions. For now, we will focus on the case of BHs with momenta in the z direction, so that P x = 0 and the linear momenta are given by The rescaled trace-free part of the extrinsic curvature for such a configuration iŝ whereÂ (1) ab andÂ (2) ab are given by Eq. (4.63) with the linear momentum (4.66). In a coordinate system adapted to the spacetime symmetry we, thus, obtain where h AB is the metric on the (D − 4)-sphere. The expression forÂ (2) ab is analogous, but with z 2 in place of z 1 and −P z in place of P z in Eq. (4.68).
Employing the relations (4.56) and the 4-metric (4.13) results in the Bowen-York type initial data for the 3-dimensional dynamical quantities prescribed by The elliptic equation for u, Eq. (4.65), writes Finally, Eq. (4.73) will be solved numerical by an extension of the TwoPunctures spectral solver [165] presented in [5].

Wave Extraction
In the previous sections we have described our framework for numerical simulations of higher dimensional BH spacetimes, providing the initial data setup as well as a formulation of the evolution system to model the dynamics of these system. The missing ingredient, namely the formalism that enables us to extract information about gravitational radiation and energy emitted throughout the collision of two (higher dimensional) BHs, will be presented in this section.

Coordinate frames
As discussed in Sec. 4.2, in our approach we perform a dimensional reduction of the D-dimensional spacetime to an effectively 4-dimensional gravity model and we consider different generalizations of "axial symmetries" to higher dimensions: either D ≥ 5 dimensional spacetimes with SO(D − 2) isometry group, or D ≥ 6 dimensional spacetimes with SO(D − 3) isometry group. Here, we focus only on the former case, which allows us to model head-on collisions of non-spinning BHs. Hereafter, we dub these spacetimes as axially symmetric. Although the corresponding symmetry manifold is the (D −3)-sphere S D−3 , the quotient manifold in our dimensional reduction is its submanifold S D−4 . The coordinate frame in which the numerical simulations are performed is where the angles φ 1 , . . . , φ D−4 describe the quotient manifold S D−4 and do not appear explicitly in the simulations. Here, z is the symmetry axis, i.e. the collision line.
In the frame (4.75), the spacetime metric has the form (cf. Eqs. (4.13) and (2.4)) is the scalar field and α, β i are the lapse function and the shift vector, respectively. It is worth noting that, although in D = 4 a general axially symmetric spacetime has non-vanishing mixed components of the metric, such as g tφ , in D ≥ 5 these components vanish in an appropriate coordinate frame. With an appropriate transformation of the four dimensional coordinates x µ , the residual symmetry left after the dimensional reduction on S D−4 can be made manifest: To extract gravitational waves far away from the symmetry axis we employ the Kodama-Ishibashi (KI) formalism [252], which generalizes the Regge-Wheeler-Zerilli [150,151] approach to higher dimensions. We require that the spacetime, far away from the BHs, is approximately spherically symmetric. Note, that spherical symmetry in D dimensions means symmetry with respect to rotations on S D−2 ; this is an approximate symmetry which holds asymptotically far away from the axis and which is manifest in the coordinate frame: Note, that throughout this section x a = (t, r) and that we have introduced polar-like coordinatesθ, θ ∈ [0, π] to "build up" the manifold S D−2 in the background, together with a radial spherical coordinate r, which is the areal coordinate in the background. The coordinate frame (4.78) is defined in such a way that the metric can be expressed as a stationary background (ds (0) ) 2 (i.e., the Tangherlini metric) plus a perturbation (ds (1) ) 2 which decays faster than 1/r D−3 for large r: Here, the Schwarzschild radius r S is related to the ADM mass M by where A D−2 is the area of the (D − 2)-sphere (see Eq. (C.59)). For instance, r S = 2M in D = 4, r 2 S = 8M/(3π) in D = 5 and r 3 S = 3M/(2π) in D = 6. When we define the coordinate frame (4.78), we also require that the coordinate θ in this frame coincides with the coordinate θ appearing in Eq. (4.77a). With this choice, the axial symmetry of the spacetime implies that as in Eq. (4.79b), and λ = sin 2 θg θθ , i.e. Eq. (4.77b). The transformation from Cartesian coordinates x µ = (t, x, y, z), in which the numerical simulation is implemented, to the coordinates (x a ,θ, θ) = (t, r,θ, θ), in which the wave extraction is performed, is given by where R = x 2 + y 2 + z 2 and by the reparametrization of the radial coordinate R =R(r) . We assume that the transformation (4.82), (4.83) yields the "Tangherlini+perturbation" splitting (4.79a), (4.79b) during the entire evolution of the system. This statement can be checked numerically by verifying the relations (see Appendix C.4): where K 0D = π 0 dθ(sinθ) D−3 , together with the axisymmetry conditions (4.77b), (4.81). The preservation of the above identities during the numerical evolution, as shown in Secs. 5.2.2 and 5.3, justifies also the identification of the time coordinate in the numerical and wave extraction frames, and our use of the KI formalism. Finally, Eqs. (4.76), (4.79a), (4.79b) yield the 3 + 1-splitting where xμ = (t, r,θ). With the 3 + 1-split, the axisymmetry conditions (4.77b), (4.81) take the form The variable r can be determined from the angular components of the metric (4.85), by averaging out hθθ, h θθ (see Appendix C.4); its explicit expression is given by (4.87) In Secs. 5.2.2 and 5.3 we will varify that the areal radius r indeed is very close to R throughout the numerical evolution.

Harmonic expansion
In the KI formalism [252] (see also [253]), the background spacetime has the form (4.79a) i.e. is the Tangherlini metric, where the x M coordinates refer to the full spacetime. The spacetime perturbations can be decomposed into spherical harmonics on the (D − 2)sphere S D−2 . They are functions of the (D − 2) angles φī = (θ, θ, φ 1 , . . . , φ D−4 ). We denote the metric of S D−2 by γīj, and with a subscript :ī the covariant derivative with respect to this metric. Finally, we denote the covariant derivative with respect to the metric g ab with a subscript |a . As discussed in [252], there are three types of spherical harmonics: • The scalar harmonics S(φī), which are solutions of with k 2 = l(l + D − 3), l = 0, 1, 2, . . . . The scalar harmonics S depend on the integer l and on other indices; we leave such dependence implicit. We also define Observe, that γījSīj = 0. Each harmonic mode of the metric perturbation δg M N = h M N can be decomposed as where f ab , f a , H L , H T are functions of x a = (t, r). Note, that in each of these expressions there is a sum over the indices of the harmonic.
For l > 1, the metric perturbations can be expressed in terms of the gauge-invariant variables [253] F where we have defined • The vector harmonics Vī(φī), which are solutions of . . . These harmonics satisfy the relation The harmonic expansion of the corresponding metric perturbations is given by Eqs. (4.91b)-(4.91c), with Sī replaced by Vī, Sīj replaced by and H L = 0.
• The tensor harmonics Tīj(φī), which are solutions of . . . These harmonics satisfy, In the D = 4 case they vanish. The harmonic expansion of the corresponding metric perturbations is given by (4.91c), with Sīj replaced by Tīj and H L = 0.

Implementation of axisymmetry
In an axially symmetric spacetime, the metric perturbations are symmetric with respect to S D−3 . Therefore, the harmonics in the expansion of h M N depend only on the angleθ. Furthermore, since there are no off-diagonal terms in the metric, the only non-vanishing g aī components are g aθ ; the only components gīj are either proportional to γīj, or all vanishing but gθθ. This implies that only scalar spherical harmonics can appear in the expansion of the metric perturbations. Indeed, if then Eq. (4.95) gives  [63,201,250]; writing explicitly the index l, they take the form where the normalization K lD is chosen such that and . By computing S lī , S līj from Eqs. (4.90) (using Eq. (4.89)) we find where we have defined Therefore, the metric perturbations are given by where h ab = h ab (t, r,θ), h aθ = h aθ (t, r,θ), hθθ = hθθ(t, r,θ), h θθ = h θθ (t, r,θ) and C cosθ). Using Eqs. (4.92b), (4.93), we obtain the gauge-invariant quantities F , F ab in terms of these functions.
As we have discussed above, this approach has been developed for D > 4, since in D = 4 the off-diagonal terms g tφ , g rφ are not vanishing in general axially symmetric spacetimes. However, we can extend our framework to D = 4 if we restrict ourselves to axially symmetric spacetimes with g tφ = g rφ = 0. In this way, we can test our formalism by comparing our results to the existing literature. For instance, we note that in D = 4 the perturbation functions are related to the expressions in Ref. [254], with the identifications We also remark that in the transverse-traceless (TT) gauge, only H T is non-vanishing, but in a generic gauge (like the one used in the numerical simulations) all these quantities are in principle non-vanishing.

Extracting gravitational waves at infinity
In the KI framework, the emitted gravitational waves are described by the master function Φ. To compute Φ in terms of the gauge-invariant quantities F , F ab one should perform a Fourier transform or a time integration (see [252]). This can be avoided if we compute directly Φ ,t , given by 3 (4.108) where k 2 = l(l + D − 3). In the TT-gauge, the gravitational perturbation is described by H T , which decays as r (D−2)/2 with increasing r, whereas the other perturbation functions have a faster decay (see [63]). In this gauge, the asymptotic behaviour of the master function is and tends to an oscillating function with constant amplitude as r → ∞. Writing the index l explicitly, the energy flux in each l−multipole is [63] The total energy emitted in the process is then

Discussion
In this chapter we have presented a framework that allows the generalisation of the current generation of 3+1 numerical codes to evolve, with relatively minor modifications, spacetimes with SO(D − 2) symmetry in D ≥ 5 dimensions and SO(D − 3) symmetry in D ≥ 6 dimensions. The key idea is a dimensional reduction of the problem along the lines of Geroch's [246] procedure that recasts the D-dimensional vacuum Einstein's equations into an effectively 4-dimensional system coupled to a scalar field. The resulting equations can be transformed straightforwardly into the (generalized) BSSN formulation that has proven to be remarkably successful in numerical evolutions of BH configurations in 3 + 1 spacetimes. We have isolated several issues related to the regularisation of the variables used in our formulation and demonstrated how all difficulties related to the coordinate singularity arising from the use of a "radius-like" coordinate can be successfully addressed in a numerical implementation. We have further illustrated how Brill-Lindquist and Bowen-York type initial data can be adapted straightforwardly to the formalism presented in this chapter. More generally, the class of problems that may be studied with our framework includes head-on collisions in D ≥ 5 and a subset of BH collisions with impact parameter and spin in D ≥ 6. Finally, we have presented our formalism to extract gravitational wave observables. The technique is based on the studies of perturbations of Tangherlini BHs by Kodama & Ishibashi [252] and has been adapted to our framework.
The formalism presented in this chapter has been implemented in the HD-Lean code [2,3], an extension to higher dimensional spacetimes of the original Lean code by Sperhake [161] developed for the modelling of 4-dimensional spacetimes.
Our numerical studies, including thorough checks of the implementation as well as head-on collision of BHs in D = 5 and D = 6 spacetime dimensions, will be discussed in the following Chapter 5.

Black hole collisions in higher
dimensional spacetimes -Results

Introduction
In high-energy physics, BHs are a central piece of the gauge/gravity duality [69,70], and are the generic outcome of particle collisions at center-of-mass energies above the Planck scale [41]. In this regime the particular nature of the particles' structure should become irrelevant, as indicated by Thorne's "hoop" conjecture [39], and "no-hair theorem"-type arguments. These arguments have been supported by numerical studies [41,42] and indicate that, in general, trans-Planckian collisions of particles are well described by collisions of highly boosted BHs. In this context, scenarios such as TeV-gravity are especially interesting, as they lower the fundamental Planck scale to the level at which BHs would be produced in cosmic ray interactions with the atmosphere and particle accelerators [48][49][50][51][52][53][54][55][56][57]. Thus, high-energy BH collisions could be used to look for signatures of extra dimensions and BH production in ground-based experiments in the forthcoming years. Note, however, that so far no BH signatures have been observed at the LHC, working now at 7 − 8T eV [58][59][60].
In the present chapter we focus mainly on evolutions of low energy head-on collisions in D = 5 and D = 6 dimensions. We have performed numerical simulations with the HD-Lean code, which implements the framework and formalism discussed Chapter 4. The code is an extension of the original Lean code [161,206], which is based on the Cactus computational toolkit [208] and uses the Carpet mesh refinement package [209,210], the apparent horizon finder AHFinderDirect [166,167] and an extension to D-dimensional spacetimes of the TwoPunctures spectral initial data solver [5,165].
This chapter is organized as follows: In Sec. 5.2 we perform several tests in order to verify the implementation. In particular, we present simulations of single BHs in D = 5 and D = 6 and compare numerical versus analytic data in case of geodesic slicing of the spacetime. Additionally, we have evolved head-on collisions of non-boosted BHs in D = 4 dimensions in order to verify the employed KI wave extraction. This setup allows for a benchmark test of the KI formalism against the well established Newman-Penrose formalism [152] as well as wave extraction based on the Regge-Wheeler-Zerilli perturbation theory [150,151] in 4-dimensional configurations. In Secs. 5.3 and 5.4 we present head-on collisions of non-boosted, equal-mass BHs in D = 5 and D = 6 dimensions. We discuss the gravitational wave emission as well as the energy outflow radiated throughout the collision. This task becomes more challenging with increasing dimension and the fact that perturbations decrease with ∼ 1 r D−3 has far reaching implications for numerical simulations: As we will see, much finer meshes are necessary to properly resolve the BH region as well as the wave extraction zone. Sec. 5.5 focuses on head-on collision of (non-boosted) unequal-mass BHs in D = 5 dimensions. We analyse the emission of gravitational radiation, its energy and linear momentum content and provide a comparison with point-particle calculations [64].
The results presented in this Chapter are based on the publications [2][3][4][5] as well as on as yet unpublished work in progress [6,7].
Note, that length and time are measured in terms of the Schwarzschild radius r S , which is related to the ADM mass M of the system via Eq. (4.80). The total radiated energy is presented in terms of the ADM mass.

Evolution of a single black hole
In order to check our numerical framework and implementation we have performed a number of tests evolving a single BH in D = 5 and D = 6 spacetime dimensions.
Geodesic slicing First, we have evolved the system by employing the geodesic slicing condition, i.e., the lapse function α and shift vector β a have been set to α =1 , β a = 0 , with a = 0, . . . , (D − 1). Although any numerical simulation using geodesic slicing is doomed to fail because geodesics hit the singularity at finite proper time, it allows for a comparison to analytic data. Thus, this strategy provides an excellent benchmark test for the implementation prior to the break-down of the simulation. We start by studying a single Schwarzschild-Tangherlini BH in D = 5 spacetime dimensions. We express the D-dimensional Tangherlini solution in a coordinate system of type (4.13) with geodesic slicing. This coordinate system may be achieved by setting a congruence of in-falling radial time-like geodesics, each geodesic starting from rest at radial coordinate r 0 , with r 0 spanning the interval [r S , +∞[, and using their proper time τ and r 0 as coordinates (instead of the standard (t, r) Schwarzschild-like coordinates). A detailed construction of the Tangherlini solution in D = 5 in these coordinates is given in Appendix C.2. Then, the (5-dimensional) line element becomes where r 0 (R) is given by Eq. (C.27).
In Fig. 5.1 we plot the metric componentγ xx along the x axis (left) and ζ/χ (right), for various values of τ using both the analytic solution and numerical data. The agreement is excellent forγ xx and good for ζ/χ. The latter shows some deviations very close to the puncture, but we believe that it is not a problem for two reasons: (i) the agreement improves for higher resolution; (ii) the mismatch does not propagate outside of the horizon.
(a)γxx along the x-axis (b) ζ/χ = λ/y 2 along the y-axis It is easy to interpret the behaviour observed forγ xx . The geodesic that starts from r = r 0 (in Schwarzschild-like coordinates) hits the physical singularity of the Tangherlini solution within proper time τ = r 2 0 /r S . Moreover, this happens at The earliest time at which the slicing hits the singularity is τ = r S , which happens at R = r S /2. On the x-axis R = x and indeed one sees in Fig. 5.1 thatγ xx diverges at x = r S /2. The divergence then extends to both larger and smaller values of x, as expected from (5.3).
A quick glance at the evolution equations (4.27) and (4.28) indicates that D = 5 may be a special case. In all these expressions there exist terms which manifestly vanish for D = 5. Therefore we next compare the numerical evolution of a single Tangherlini BH in D = 6 spacetime dimensions with the analytic solution, using geodesic slicing. This comparison has been more difficult than in case of D = 5, because the line element analogous to (5.2) cannot be obtained in a simple analytic form. In Appendix C.2 we demonstrate how a semi-analytic solution can be obtained for the metric. In Fig. 5.2 we compare this expression with the numerical values at times τ = 0.5 r S , 0.7 r S and 0.72 r S . The agreement is excellent and demonstrates that our implementation is correct.
Moving puncture gauge In the second part we evolve a single BH in D = 5 and D = 6 dimensions employing the modified moving puncture gauge, Eqs. (4.30) and (4.31), which is expected to yield long-term stable evolutions. In particular, we perform a convergence analysis of the Hamiltonian and momentum constraint obtained from these simulations. We choose Brill-Lindquist type initial data (cf. Sec We have re-scaled the high resolution run by Q 4 = 1.5 4 , the factor anticipated for fourth order convergence [106]. One can see that there is some noise, but the overall convergence is acceptable. For the Hamiltonian constraint the convergence is essentially 4th order and for the momentum constraint it decreases slightly towards 2nd or 3rd order in patches. From experience in 3 + 1-dimensional numerical relativity this is perfectly acceptable, especially given the fact that prolongation in time is second-order accurate. Next, we have evolved a single puncture in D = 6 spacetime dimensions. Whereas we have been able to obtain long-term stable simulations of single BHs in D = 5 by modifying the moving puncture gauge conditions in a very straight-forward manner, the D = 6 case appeared to be much more difficult. In fact, we could only accomplish this task by extensively studying various choices of the gauge parameter (µ L , µ λ , ξ Γ , ξ λ ) in Eqs. (4.30) and (4.31), involving a large number of numerical experiments.
In Fig. 5.4 we show the Hamiltonian and the y-component of the momentum constraint along the y-axis obtained for a single puncture in D = 6 at t = 8 r S . As for D = 5, the high resolution result is amplified by a factor 1.5 4 expected for fourth order convergence [106]. While the convergence appears to be closer to second order in some patches of the momentum constraint, the results are clearly compatible with the numerical discretization.
The tests presented in this section demonstrate that our implementation of the higher dimensional framework is correct and the code exhibits convergent results.

Benchmark test -head-on collisions from rest in D = 4
In order to test our implementation of the KI formalism, discussed in Sec. 4.6, we have simulated head-on collision of an equal-mass, non-spinning BH binary initially at rest in D = 4 spacetime dimensions. Head-on collisions in 4-dimensional spacetimes have been studied extensively in the literature and provide valuable opportunities to calibrate the wave extraction formalism. In particular, we compare our results derived from the KI formalism against those obtained from both, the Regge-Wheeler-Zerilli wave extraction and the Newman-Penrose framework; see, e.g., Refs. [157][158][159]161,254] for corresponding literature studies. The parameters used in the simulations presented in this section are shown in Table 5.1.
In order to perform these tests, we need to relate our master function Φ of Sec. 4.6.4  to the variables used in traditional four dimensional studies. Specifically, a straightforward calculation shows that the Zerilli wavefunctionΦ adopted in Ref. [254] for l = 2 multipoles and the outgoing Weyl scalar Ψ 4 used in [161] can be expressed in terms of Φ according toΦ Tests on the numerical coordinates The procedure described in Section 4.6 assumes that the numerical spacetime consists of a small deviation from the Schwarzschild-Tangherlini metric. In order to ensure that the gravitational waves are extracted in an appropriate coordinate system we perform a number of checks. First, we test the relations (4.84a), (4.84b) and (4.84c). In Fig. 5.5 we show G tt , i.e., the difference between the numerically calculated component g tt , averaged over the extraction sphere and the corresponding component of the assumed background metric. Here we evaluate the background metric by assuming, as a first approximation, that the Schwarzschild radius of the BH is r S = r S,1 + r S,2 .  Table 5.1.: Grid structure and initial parameters of the head-on collisions starting from rest in D = 4. The grid setup is given in terms of the "radii" of the individual refinement levels, in units of r S , as well as the resolution near the punctures h (see Sec. II E in [161] for details). d is the initial coordinate separation of the two punctures and L denotes the proper initial separation. The deviation of the full 4-metric from the Schwarzschild-Tangherlini background decreases as the extraction radius increases. Indeed, a straightforward calculation shows that a deviation δr S of the Schwarzschild radius from the background value leads to G tt ∼ δr D−3 S /r D−3 , i.e., G tt ∼ δr S /r for D = 4. In the left panel of Fig. 5.5 we therefore show the deviation G tt re-scaled by r. We further apply a time shift to account for the different propagation time of the wave to reach the extraction radii. As shown in the figure, the deviation from the Schwarzschild line element is small and decreases ∼ 1/r in accordance with our expectation. We also note that a deviation δr S represents a monopole perturbation of the background which decouples from the quadrupole wave signal at perturbative order, so that its impact on our results is further reduced.
In summary, we can give an uncertainty estimate for the approximation r S = r S,1 +r S,2 for the Schwarzschild radius of the final BH, which ignores the energy loss through gravitational radiation. As demonstrated by the left panel of Fig. 5.5, at late times |(R ex /r S ) G tt | ∼ 0.01, and, since r R ex (as we discuss below), we obtain the upper bound This crude analysis sets an upper bound of ∼ 1% on the fraction of the center-ofmass energy radiated as gravitational waves. We further note that the close agreement between g tt and its Tangherlini counterpart implies that the time coordinate employed in the numerical simulation and the Tangherlini coordinate time coincide. By analysing G tr and G rr in the same manner, we find that relations (4.84a)-(4.84c) are satisfied with an accuracy of one part in 10 2 throughout the evolution, and one part in 10 3 at late times, when the spacetime consists of a single distorted black hole. In practice, gravitational waves are extracted on spherical shells of constant coordinate radius. The significance of the areal radius associated with such a coordinate sphere in the context of extrapolation of GW signals has been studied in detail in Ref. [255]. For our purposes, the most important question is to what extent gauge effects change the areal radius (4.87) of our extraction spheres. For this purpose, we show its time evolution in the right panel of Fig. 5.5 for different values of R ex . The reassuring result is that the areal radius exceeds its coordinate counterpart by about 1 % at R ex = 50 r S and remains nearly constant in time.
Waveforms As a benchmark for our wave extraction, we compare our results obtained with independent wave extraction tools: (i) the explicitly four dimensional Zerilli formalism and (ii) the Newman-Penrose scalars. For this purpose we have evolved model HD4 f and extracted the Zerilli function according to the procedure described in [254] (see also Eqs. (4.107a)-(4.107d)) and the Newman Penrose scalar Ψ 4 as summarized in [161]. These are compared with the KI wave function Φ ,t and its time derivative Φ ,tt in Fig. 5.6. Except for a small amount of high frequency noise in the junk radiation at t ≈ 25r S , we observe excellent agreement between the different extraction methods.
Next, we consider the dependence of the wave signal on the extraction radius. In Fig. 5.7 we show the l = 2 component of Φ ,t extracted at three different radii and shifted in time by R ex . As is apparent from the figure, the wave function shows little variation with R ex at large distances, in agreement with expectations.
A further test of the wave signal arises from its late-time behaviour which is dominated by the BH ringdown [211], an exponentially damped sinusoid of the form e −iωt , with ω being a characteristic frequency called quasinormal mode (QNM) frequency. Using wellknown methods [211,256,257], we estimate this frequency to be r S ω ∼ 0.746 ± 0.002 − i (0.176 ± 0.002). This can be compared with theoretical predictions from a linearized approach, yielding r S ω = 0.747344 − i 0.177925 [211].
Finally, we consider the numerical convergence of our results. In Fig. 5.8, we plot the differences obtained for Φ ,t extracted at R ex = 30 r S , using the different resolutions of the three models HD4 listed in Table 5.1. The differences thus obtained are consistent with 4 th order convergence. This implies a discretization error in the l = 2 component of Φ ,t of about 4% for the grid resolutions used in this work. Radiated energy Once the KI function Φ ,t is known, the energy flux can be computed from Eq. (4.110). For comparison, we have also determined the flux from the outgoing Newman Penrose scalar Ψ 4 according to Eq. (22) in Ref. [8]. The flux and energy radiated in the l = 2 multipole, obtained with the two methods at R ex = 50 r S is shown in Fig. 5.9 and demonstrates agreement within the numerical uncertainties of about 4 % for either result.
We obtain an integrated energy of 5.5 × 10 −4 M and 5.3 × 10 −4 M , respectively, for the gravitational wave energy radiated in l = 2, where M denotes the center-of-mass  We plot the differences between the low and medium resolution (black solid line) and medium and high resolution (red dashed line) run. The latter is re-scaled by the factor Q 4 = 1.58 expected for 4 th order convergence [8]. energy.
The energy in the l = 2 mode is known to contain more than 99% of the total radiated energy [161]. Our analysis is compatible with this finding; while the energy in the l = 3 mode is zero by symmetry, our result for the energy in the l = 4 mode obtained from the KI master function is three orders of magnitude smaller than that of the l = 2 contribution.

Head-on collisions from rest in D = 5
Having tested the wave extraction formalism in D = 4 dimensions in the previous section, we now turn our attention to the results obtained for head-on collisions of BHs in 5-dimensional spacetimes. As before, we consider nonspinning BH binaries initially at rest with coordinate separation d/r S . Note, that in five spacetime dimensions the Schwarzschild radius is related to the ADM mass M via Eq. (4.80), We therefore define the "total" Schwarzschild radius r S such that r 2 S = r 2 S,1 + r 2 S,2 . By using this definition, r S has the physical dimension of a length and provides a suitable unit for measuring both, results and grid setup.
As summarized in Table 5 : Grid structure and initial parameters of the head-on collisions starting from rest in D = 5. The grid setup is given in terms of the "radii" of the individual refinement levels, in units of r S , as well as the resolution near the punctures h (see Sec. II E in [161] for details). d/r S is the initial coordinate separation of the two punctures and L/r S denotes the proper initial separation.

Tests on the numerical coordinates
In order to verify the assumptions underlying our formalism, we have analysed the coordinate system in analogy to Sec. 5.2.2. First, we have evaluated the averaged areal radius on extraction spheres of constant coordinate radius.
The result shown in the left panel of Fig. 5.10 demonstrates that the coordinate and areal radius agree within about 1 part in 10 4 for R ex ≥ 40 r S . The Tangherlini coordinate r equals by construction the areal radius and our approximation of setting r ≈ R ex in the wave extraction zone is satisfied with high precision.
Second, we evaluate the deviation of the metric components according to Eqs. (4.84a)-(4.84c). From the discussion in Sec. 5.2.2 we expect G tt ∼ r 2 /r 2 S in D = 5. Our results in the right panel of Fig. 5.10 confirm this expectation and demonstrate that our spacetime is indeed perturbatively close to that of a Tangherlini metric at sufficient distances from the BHs; deviations in G tt are well below 1 part in 10 3 at R ex = 60 r S . Furthermore, we can estimate the crudeness of the approximation r 2 S = r 2 S,1 + r 2 S,2 for the Schwarzschild radius of the final BH: as shown in the right panel of Fig. 5.10, at late times |R 2 ex /r 2 S G tt | ∼ 0.01; this value gives an upper bound on the radiated energy.
For the third test, we recall that our higher dimensional implementation does not employ the full isometry group of the S 2 sphere in D = 5 dimensions and axial symmetry manifests itself instead in the conditions (4.86) on the metric components and the scalar field. We find these conditions to be satisfied within 1 part in 10 8 and 1 part in 10 16 , respectively, in our numerical simulations which thus represent axially symmetric configurations with high precision.

Newtonian collision time
An estimate of the time at which the BHs "collide", can be obtained by considering a Newtonian approximation to the kinematics of two point particles in D = 5. In the weak-field regime, Einstein's equations reduce to "Newton's law" a = −∇B(x), with h 00 = −2B(x) = r D−3 S /2r D−3 . The Newtonian time it takes for two point-masses (with Schwarzschild parameters r S,1 and r S,2 ) to collide from rest with initial distance L in D dimensions is then given by .

(5.9)
For D = 4, one recovers the standard result t free-fall = π L 3 /r 3 S r S , whereas for D = 5 we get t free-fall = (L/r S ) 2 r S . (5.10) In general relativity, BH trajectories and merger times are intrinsically observer dependent quantities. For our comparison with Newtonian estimates we have chosen relativistic trajectories as viewed by observers adapted to the numerical coordinate system. While the lack of fundamentally gauge invariant analogues in general relativity prevents us from deriving rigorous conclusions, we believe such a comparison to serve the intuitive interpretation of results obtained within the moving puncture gauge. Bearing in mind these caveats, we plot in Fig. 5.11 the analytical estimate of the Newtonian time of collision, together with the numerically computed time of formation of a common apparent horizon. Also shown in Fig. 5.11 is the time at which the separation between the individual BH's puncture trajectory decreases below the Schwarzschild parameter r S . The remarkable agreement provides yet another example of how well numerically successful gauge conditions appear to be adapted to the BH kinematics. It is beyond the scope of this paper to investigate whether this is coincidental or whether such agreement is necessary or at least helpful for gauge conditions to ensure numerical stability. Suffices it to say at this stage that similar conclusions were reached by Anninos et al. [204] and Lovelace et al. [258] in similar four dimensional scenarios.

Waveforms
We now discuss in detail the gravitational wave signal generated by the head-on collision of two BHs in D = 5 dimensions. For this purpose, we plot in Fig. 5.12 the l = 2 multipole of the KI function Φ ,t for model HD5e f obtained at different extraction radii. Qualitatively, the signal looks similar to that shown in the left panel of Fig. 5.6 for D = 4. A small spurious wavepulse due to the initial data construction is visible at ∆t ≈ 0. This so-called "junk radiation" increases in magnitude if the simulation starts with smaller initial separation of the BHs. We return to this issue further below, when we study the  dependence of the gravitational radiation on the initial BH separation. The physical part of the waveform is dominated by the merger signal around ∆t = 50 r S , followed by the (exponentially damped) ringdown, whereas the infall of the BHs before ∆t = 40 r S does not produce a significant amount of gravitational waves. Comparison of the waveforms extracted at different radii demonstrates excellent agreement, in particular for those extracted at R ex = 40 r S and 60 r S . Extrapolation of the radiated energy to infinite extraction radius yield a relative error of 5 % at R ex = 60 r S , indicating that such radii are adequate for the analysis presented in this work. Due to symmetry, no gravitational waves are emitted in the l = 3 multipole, so that l = 4 represents the second strongest contribution to the wave signal. As demonstrated in the right panel of Fig. 5.6, however, its amplitude is two orders of magnitude below that of the quadrupole.
A convergence analysis also using the lower resolution simulations of models HD5e c and HD5e m is shown in Fig. 5.13 and demonstrates overall convergence of third to fourth order, consistent with the numerical implementation. From this analysis we obtain a conservative estimate of about 4% for the discretization error in the waveform.
In practice, numerical simulations will always start with a finite separation of the two BHs. In order to assess how accurately we are thus able to approximate an infall from infinity, we have varied the initial separation for models HD5a to HD5f as summarized in Table 5.2. For small d we observe two effects which make the physical interpretation of models HD5a−HD5c difficult. First, the amplitude of the spurious initial radiation increases and second, the shorter infall time causes an overlap of this spurious radiation with the merger signal. As demonstrated in Fig. 5.14 for models HD5e and HD5f, however, we can safely neglect the spurious radiation as well as the impact of a final initial separation, provided we use a sufficiently large initial distance d 6 r S of the BH binary. Here, we compare the radiation emitted during the head-on collision of BHs starting from rest with initial separations 6.37 r S and 10.37 r S . The waveforms have been shifted in time by the extraction radius R ex = 60 r S and such that the formation of a common apparent horizon occurs at ∆t = 0. The merger signal starting around ∆t = 0 r S shows excellent agreement for the two configurations and is not affected by the spurious signal visible for HD5e at ∆t ≈ −50 r S .
We conclude this discussion with two aspects of the post-merger part of the gravitational radiation, the ringdown and the possibility of GW tails. After formation of a common horizon, the waveform is dominated by an exponentially damped sinusoid, as the merged BH rings down into a stationary state. By fitting our results with an exponentially damped sinusoid, we obtain the characteristic frequency This value is in excellent agreement with perturbative calculations, which predict a lowest quasinormal frequency r S ω = 0.9477 − i 0.2561 for l = 2 [211,250,259]. A well known feature in gravitational waveforms generated in BH spacetimes with D = 4 as well as D > 4 are the so-called power-law tails [260][261][262][263]. In odd dimensional spacetimes an additional, different kind of late-time power tails arises, which does not depend on the presence of a BH. These are due to a peculiar behavior of the wavepropagation in flat odd dimensional spacetimes because the Green's function has support inside the entire light-cone [263]. We have attempted to identify such power-law tails in our signal at late times, by subtracting a best-fit ringdown waveform. Unfortunately, we cannot, at this stage, report any evidence of such a power-law in our results, most likely because the low amplitude tails are buried in numerical noise.

Radiated energy
Comparison of Figs. 5.7 and 5.14 for the GW quadrupole in D = 4 and D = 5 shows a larger wave amplitude in the five dimensional case and thus indicates that this case may radiate more energy. We now investigate this question quantitatively by calculating the energy flux from the KI master function via Eq. (4.110). The fluxes thus obtained for the l = 2 multipole of models HD5e f and HD5f in Table 5.2, extracted at R ex = 60 r S , are shown in Fig. 5.15. As in the case of the KI master function in Fig. 5.14, we see no significant variation of the flux for the two different initial separations. The flux reaches a maximum value of dE/dt ∼ 3.4 × 10 −4 r S , and is then dominated by the ringdown flux. The energy flux from the l = 4 mode is typically four orders of magnitude smaller; this is consistent with the factor of 100 difference of the corresponding wave multipoles observed in Fig. 5.12, and the quadratic dependence of the flux on the wave amplitude.
The total integrated energy emitted throughout the head-on collision is presented in the left panel of Fig. 5.16. We find that a fraction of E rad /M = (8.9 ± 0.6) × 10 −4 of the center-of-mass energy is emitted in the form of gravitational radiation. We have verified  for these models that the amount of energy contained in the spurious radiation is about three orders of magnitude smaller than in the physical merger signal.
An independent estimate for the radiated energy can be obtained from the apparent horizon area A 4 in the effective four dimensional spacetime by using the spherical symmetry of the post-merger remnant BH. Energy balance then implies that the energy E radiated in the form of GWs is given by where M AH is the apparent horizon mass. The estimate E/M is shown in Fig. 5. 16 and reveals a behavior qualitatively similar to a damped sinusoid with constant offset. Indeed, by using a least squares fit, we obtain a complex frequency r S ω ∼ 0.97 − i 0.29, again similar to the fundamental l = 2 quasinormal mode frequency (see discussion around Eq. (5.11)). At late times, 1−M AH /M asymptotes to 1−M AH /M ∼ (9.3 ± 0.8)× 10 −4 which agrees very well with the GW estimate, within the numerical uncertainties.

Head-on collisions from rest in D = 6
Here, I present evolutions of non-boosted, equal-mass BH binaries colliding head-on in D = 6 spacetime dimensions. This goal proved to be numerically more challenging than previous studies in D = 5 and is subject to detailed investigations. One reason might be the increasing fall-off with spacetime dimension which goes as ∼ 1 r D−3 . Thus, not only do we need higher resolution grids to properly resolve the BH region, but also gradients become steeper and might cause problems. Another reason might be the formulation itself. However, in Secs. 4.2 and 4.3 I have shown that the employed formulation in the continuum limit is strongly hyperbolic and therefore well-posed for generic spacetime dimension D ≥ 5 and most gauge parameters (cf. Sec. 4.4). Nevertheless, this is only a necessary but not sufficient condition for the numerical stability of the system. In practice, I have only accomplished this task by extensively studying a wide range of parameter choices (ξ Γ , ξ λ ) of the Γ-driver shift condition, Eq. (4.31). The procedure has resulted in hundreds of numerical experiments. Finally, the specification ξ Γ = 0.25 and ξ λ = 0.0 has led to succesful, long-term stable simulations. Note, however, that there is no strong mathematical reason why this choice turns out to work better than others and it seems to depend on the spacetime dimension in a non-trivial way. Additionally, the particular setting of these parameters might depend sensitively on the particular configuration. In this section, I present the thus obtained results. The BHs, each with mass parameter r 3 S,i = 0.5, start from rest with initial coordinate distance d = 6.37 r S , d = 9.56 r S or d = 12.74 r S . The specifications of the simulations, such as the grid setup and resolution h of the innermost refinement level, the initial coordinate distance d/r S and the employed damping parameter η β in the Γ-driver shift condition are summarized in Table 5 Table 5.3.: Grid structure and initial parameters of the BH head-on collisions in D = 6. The grid setup is given in terms of the "radii" of the individual refinement levels, in units of r S , as well as the resolution near the punctures h (see Sec. II E in [161] for details). d/r S is the initial coordinate separation of the two punctures. η β is the damping parameter in the Γ-driver shift condition, Eq. (4.31).

Waveforms
First, I discuss the gravitational wave signal emitted throughout the head-on collision of two BHs in D = 6 dimensions. In the left panel of Fig. 5.17 I plot the l = 2 multipole of the KI masterfunction Φ ,t , Eq. (4.108), for model HD6a f , extracted at different radii R ex = 30, 40, 50, 60 r S . Qualitatively, the signal looks similar to those presented in the left panel of Figs. 5.6 and 5.12 for D = 4 and D = 5 dimensions, respectively. A small spurious wavepulse, the so-called "junk radiation" which is due to the initial data construction, would appear at ∆t ≈ 0 r S . However, the amplitude of the junk radiation is significantly less than in the previously studied cases with comparable setups in D = 4 and D = 5 dimensions and, in fact, not visible in Fig. 5.17. The physical part of the waveform is dominated by the merger signal around ∆t = 75 r S , followed by the exponentially damped ringdown. The comparison of the waveforms extracted at different radii demonstrates excellent agreement. I will study the influence of the extraction at finite radii in more detail further below during the analysis of the radiated energy and give an estimate of the thus induced error.
I further check our results by estimating the QNM ringdown frequency of the l = 2 and l = 4 multipoles and compare them to the close-limit approximation [250,264]. The QNM ringdown signal is an exponentially damped sinusoid which dominates the gravitational wave signal after the formation of a common horizon. Therefore, I perform a fit of the numerically computed l = 2 and l = 4 waveforms to this functional form. The resulting characteristic QNM frequencies are summarized in Table 5.4, where I list the dimensionless real and imaginary parts r S ω R,l and r S ω I,l . A perturbative, close-limit analysis presented in [250,264] predicts the frequencies and a complementary, point particle computation for the l = 2 mode [211,264] gives Whereas the values for simulation HD6a f with small initial separation of d = 6.37 r S is in good agreement with the theoretical predictions [211,250,264], Eq. (5.13), within about 6.5% for l = 2 and 9% for l = 4, (i.e., within numerical errors), the frequencies of the wave signal obtained in the case of larger separation d = 9.56, 12.74 r S deviate from the close-limit approximation [211,250,264] by about 20% in both multipoles.
In the right panel of Fig. 5.17 I present the l = 2 multipole of the KI masterfunction Φ ,t considering the initial coordinate distances d = 6.37 r S , d = 9.56 r S and d = 12.74 r S . Whereas the waveform for d = 6.37 r S has a significantly smaller amplitude and higher frequency, the wave signal for the cases with initial distance d = 9.56 r S and d = 12.74 r S basically coincide, indicating that the BHs are initially far enough apart to mimic infinite separation. However, the analysis of the QNM frequencies has revealed a significant deviation from the theoretically predicted value in the two latter cases. Therefore, unless denoted otherwise, I will present results obtained for the high resolution run HD6a f with initial coordinate distance d = 6.37 r S , extracted at R ex = 40 r S .
In the left panel of Fig. 5.18 I present the l = 2 and l = 4 multipole of the KI function of run HD6a f , where the latter has been re-scaled by a factor 50. The l = 4 multipole represents the second strongest contribution to the wave signal, with an amplitude which is about two orders of magnitude smaller than the dominant l = 2 mode, whereas the l = 3 multipole vanishes due to symmetry.
Furthermore, I perform a convergence test of the l = 2 waveforms extracted at R ex = 40 r S , presented in the right panel of Fig. 5.18, by employing model HD6a. Specifically, I have run the configuration with resolutions h c = r S /84, h m = r S /88 and h = r S /92, which I denote as coarse, medium and high resolution. I compare the differences of the amplitudes between the first two simulations with the difference between the medium and high resolution run. The latter has been re-scaled by a factor Q 4 = 1.26 demonstrating fourth order convergence, consistent with the numerical implementation. The thus obtained discretization error in the amplitude of the waveform is about ∆|Φ l=2 ,t |/|Φ l=2 ,t | ≤ 11%. The analogous analysis performed for the phase φ reveals a discretization error up to ∆φ l=2 /φ l=2 ≤ 10%.
Run  Table 5.4.: QNM frequencies r S ω = r S ω R + ı r S ω I of the l = 2 and l = 4 multipoles of the KI wavefunction, resulting from a fit of the numerical waveform to an exponentially damped sinusoid.

Radiated energy
I next investigate the emission of gravitational wave energy in more detail. For this purpose, I first plot the energy flux emitted in the l = 2 multipole of the KI wavefunction as well as as the total radiated energy E/M in Fig. 5.19. I present the radiated energy and flux, extracted at R ex = 40 r S , for all three types of models corresponding to headon collisions of BHs with initial (coordinate) separation d = 6.27, 9.56, 12.74 r S . The signals shown in Fig. 5.19 have been shifted in time by the extraction radius R ex = 40 r S and such that the maxima in the energy flux overlap. As I have previously observed in the analysis of the waveforms, also the radiated energy and energy flux resulting from the HD6a f run on the one hand and the HD6b and HD6c runs on the other differ enormously.
In the left panel of Fig. 5.20 I present the radiated energy E/M as function of the extraction radius R ex = 25, . . . , 60 r S for model HD6a f . Whereas the first five values, which have been extracted on the same refinement level, vary little, the energy measured at R ex = 50, 55, 60 r S , located at the next coarser refinement level, shows significant deviations. I find a similar behaviour also for models HD6b and HD6c. Therefore, I consider only the total energy extracted at R ex = 25, 30, 35, 40, 45 r S and perform a fit of the form The limit R ex /r S → ∞ allows us to extrapolate the energy E ∞ /M measured at infinite radius. I summarize the extrapolated energy E ∞ /M in Table 5.5.
I account for the error arising from the measurment at finite extraction radii by comparing E ∞ /M with the value obtained from a second fit including the quadratic term A 2 /(R ex /r S ) 2 . Then, this error is estimated to be about ∆E/E ∞ ∼ 1%, 1.2%, 6% for models HD6a f , HD6b and HD6c, respectively.
In order to investigate the numerical accuracy I have performed a convergence analysis, presented in the right panel of Fig. 5.20. For this purpose, I have employed model HD6a representing a head-on collision of two BHs with initial distance d = 6.37 r S , simulated with resolutions h c = r S /84, h m = r S /88 and h f = r S /92. The corresponding convergence plot is presented in the right panel of Fig. 5.20. I depict the differences of the total radiated energy E/M , extracted at R ex = 40 r S , between the coarse-medium and medium-high resolutions simulations. The latter difference has been re-scaled by the factor Q 4 = 1.25, thus demonstrating convergence of fourth order consistent with the implementation. The discretization error is estimated to be about ∆E/E ≤ 12%. I have not performed convergence runs for the remaining models HD6b and HD6c and assume that the discretization error is of the same order. I list the error in the total radiated energy E ∞ /M due to discretization and finite extraction in Table 5.5. Note, that the energy appears to depend non-monotonously on the initial distance. However, the expected monotonous behaviour is still covered within the error estimates. Nevertheless, because the extracted total energy E ∞ /M appears to depend non-monotonously on the (initial) distance d/r S , and because I only considered three different initial separations I have not extrapolated the radiated energy to infinite initial separation d/r S → ∞. Instead, I consider model HD6a f which has exhibited waveforms with the correct QNM frequencies. Then, I estimate the total radiated energy emitted throughout the headon collision of two non-boosted, equal-mass BHs in D = 6 spacetime dimensions to be E ∞ /M ∼ (0.104 ± 0.013)%.
Nevertheless, this result has to be taken with a grain of salt until ongoing investigations conclusively answer the question arosen by the observed discrepancies. Possible sources of error include (i) lack of numerical accuracy or (ii) radiation effects. In order to test the first possibility, simulations of model HD6c with higher resolutions are underway. Because the radiated energy is expected to be small, I have assumed the total Schwarzschild radius r 3 S = r 3 S,1 + r 3 S,2 = 1. However, the radiation effects might be non-negligible and are subject to further analysis.

Head-on collsions of unequal mass black holes in D = 5
In the previous sections we reported our results of head-on collisions of equal-mass BHs in D = 5 and D = 6 spacetime dimensions. In the present section, we wish to extend that study to the case of unequal mass BH binaries. This is an interesting extension for several reasons, perhaps the most important of which is the non-trivial comparison with pointparticle (PP) calculations in the linearized regime. We will compare radiated energy, momentum and multipolar dependence of our full nonlinear results with results from linearized Einstein equations. It turns out that the agreement is remarkable, providing an outstanding consistency check on our codes and results. A thorough analysis of the linearized Einstein equations has been done by Berti et al [64].
Here, we have evolved BH binaries, colliding head-on from rest with mass ratios q ≡ M 1 /M 2 = r D−3 S,1 /r D−3 S,2 = 1, 1/2, 1/3, 1/4, where M i is the mass of the i-th BH. The mass parameter r D−3 S,1 /r D−3 S of the smaller BH is given in Table 5.6 and we adapt the value of the second BH accordingly. The initial coordinate separation of the two BHs is set to d/r S = 6.37 which translates to a proper initial separation of L/r S = 6.33. Further details of the setup of the simulations are summarized in Table 5.6 and 5.7. Unless denoted otherwise, our discussion will always refer to the highest resolution runs with h f = r S /84, h f = r S /102.9, h f = r S /118.8 and h f = r S /132.8 for models R5a, R5b, R5c and R5d f in Table 5.7, respectively. The energy flux is computed according to Eq. (4.110) (see Eq. (21) in Ref. [63] for the corresponding expression in Fourier space).  Table 5.7.: Grid structure for the simulations of head-on collisions from rest in D = 5. The grid setup is given in terms of the "radii" of the individual refinement levels, as well as the resolution near the punctures h, in units of r S (see Sec. II E in [161] for details).
The momentum flux can be obtained from with n i a unit radial vector on the sphere at infinity S ∞ . This results in an infinite series coupling different multipoles. Using only the first two terms in the series, we find, for instance, that in D = 5 the momentum flux in the collision direction is given by Here, Φ l ,t is the l−pole component of the KI gauge-invariant wavefunction, Eq. (4.108), [3,252,253]. From the momentum radiated, the recoil velocity of the system can be obtained as

Waveforms
In Fig. 5.21 we show the l = 2, 3, 4 waveforms for different mass ratios, zoomed in around the time of the merger. The waveforms have been shifted in time such that ∆t = (t − R ex − t CAH ) = 0 r S corresponds to the time t CAH at which the common apparent horizon forms and taking into account the propagation time of the waves to the extraction radius R ex = 60, 49, 42.4, 37.9 r S . The waveform is similar to previous four-dimensional results (see, e.g. Ref. [205], a more detailed study has been given by Sperhake et al [1]). Although not shown in Fig. 5.21 we observe a small, spurious signal starting around (t − R ex )/r S ≈ 0, which is an artifact of the initial data. The actual physical part of the waveform is dominated by the merger signal at ∆t/r S ≈ 0 followed by the quasinormal ringdown. We estimate that the different ringdown modes are given by These results agree well, and within uncertainties, with estimates from linearized theory [63,211,259,265], providing a strong consistency check on our results. Finally, we consider numerical convergence of our waveforms. This study is summarized in Fig. 5.22 for the l = 2 mode of the KI wavefunction, and for the most challenging mass ratio, q = 1/4, model R5d in Table 5.7. We have evolved this setup at three different resolutions, namely h c = r S /113.8, h m = r S /123.3 and h f = r S /132.8, which we will refer to as "coarse", "medium" and "high" resolution in the following. We show the difference  Table 5.7. We show the differences between the coarse and medium resolution waveform (black solid line) and the medium and high resolution waveform (red dashed line). The latter has been amplified by the factor Q 4 = 1.47, indicating fourth order convergence. The curves have been shifted in time such that the formation of the common apparent horizon corresponds to ∆t/r S = 0 and taking into account the time that it takes the waves to propagate to the extraction radius.
between the coarse and medium as well as between the medium and high resolution waveforms. The latter has been amplified by the factor Q 4 = 1.47, which indicates fourth order convergence. We obtain the same order of accuracy for the higher modes. The discretization error in the waveforms is estimated to be ≈ 1.5%. Table 5.8 lists some of the most important physical quantities which characterize the head-on collision of BHs in D = 5. In particular, we show the radiated energy in units of total mass M , and the recoil velocity of the final BH in km/s. The maximum amount of energy is emitted in the equal mass case (E rad /M = 0.089% as presented in Sec. 5.3 [3]), and it decreases for smaller mass ratios. We estimate the error in the radiated energy to be about 5 %. These results have been obtained by integrating the energy flux as given by the KI master wavefunction. We have also estimated the radiated energy using properties of the apparent horizon as described in Sec. 5.3. We estimate the discretization error to be about ≈ 10 % when using this method. The apparent horizon estimate for the total radiated energy is shown in parenthesis in Table 5.8, and is consistent with the flux computation within numerical uncertainties.  Table 5.8.: Summary of our results concerning unequal-mass head-on collisions of BHs in D = 5. We show the total radiated energy E/M as measured from the energy flux at R ex , the quantity in parenthesis refers to the estimate obtained using properties of the apparent horizon (see Sec. 5.3 for details). The next three columns show the fraction of energy E l excited in the l-th mode as compared to the total radiated energy. The last column refers to the recoil velocity v recoil in km/s. dependence for the total radiated energy [266], E rad /M ∝ η 2 , where η = q/(1 + q) 2 is the dimensionless reduced mass. For clarity, we show the ratio E rad /(M η 2 ) in the top panel of Fig. 5.23, which can be seen to depend very weakly on η 2 . We can refine this argument by fitting our numerical results to an improved expression of the form

Radiated energy
Moreover, the following expressions for the multipolar content provide a good fit to our numerical data, Linearized, PP calculations presented by Berti et al [64] show that in the limit of zero mass ratio one obtains which agrees with the extrapolation of our numerical results within less than 1%. The multipole contents in the PP limit are still in very good agreement with the extrapolation of our full numerical results to the zero mass ratio limit. In fact, bearing in mind that we are extrapolating from mass ratios of 1/4 down to the zero mass ratio limit, the agreement is impressive. Finally, all these results are consistent with the fact that higher multipoles contribute more to the radiation than in D = 4, where for instance the l = 3 mode contributed roughly 10% of the total energy in the PP limit [191]. Linearized, point-particle calculations show that the trend is consistent and continues in higher dimensions [64], which might mean that accurate wave extraction will become extremely difficult, as higher resolution is necessary to resolve higher-l modes [171].

Radiated momentum
For unequal-mass collisions, the asymmetric emission of radiation along the collision axis causes a net momentum to be carried by gravitational waves. As such, the final BH will "recoil", according to Eq. (5.18). Momentum fluxes and recoil velocity for different mass ratios are shown in Fig. 5.24 and the bottom panel of Fig. 5.23, respectively. We estimate the errors in the recoil velocity to be ≈ 5 %. The general functional form for the dependence of momentum on the mass parameters of the individual holes has been worked out by M. Lemos in generic spacetime dimensions [266] and is the same as in four dimensions By fitting this function to our numerical data, we obtain C = 716 km/s. Observe that v recoil reaches a maximum value at q = 2 − ϕ 0.38, where ϕ is the golden ratio. The quality of the fit can be seen in the bottom panel of Fig. 5.23, where we overplot the numerical data points with the fitting function, Eq. (5.24). This exercise is interesting because we can again extrapolate our results to the PP limit. Berti et al. [64] find v recoil =779q 2 km/s , (5.25) in reasonably good agreement (better than 10%) with our extrapolation. We note that momentum emission is given by a non-trivial interference between different multipoles, so this is a non-trivial agreement.

Discussion
In I have shown that the total energy released in the form of gravitational waves is approximately (0.089 ± 0.006)% in D = 5 and (0.104 ± 0.013)% in D = 6 of the initial center-of-mass energy of the system, for a head-on collision of two BHs starting from rest. As a comparison, the analogous process in D = 4 releases a slightly smaller quantity: (0.055±0.006)%. I summarize the main results for head-on collisions of two BHs starting from rest in D = 4, 5, 6 spacetime dimensions in Table 5.9.
The radiative efficiency E rad /M in Table 5.9 shows that head-on collisions starting from rest in D = 5 and D = 6 dimensions generate about 1.6 and 1.9 times as much GW energy as their four dimensional counterparts. It will be very interesting to investigate to what extent this observation holds for wider classes of BH collisions. I can compare the radiation efficiency with the upper limit derived by Hawking [267] from the requirement that the horizon area must not decrease in the collision. This leads to the area bound Evidently, this bound decreases with dimensionality, while in the present computation it increases with increasing dimensions. As also shown in the table, the generation of GWs in head-on collisions starting from rest is about 3 orders of magnitude below this bound. In four dimensions it has already been demonstrated that there exist more violent processes which release more radiation than the headon collisions considered in this work [214,215,217]. A generalization to more generic configurations such as collisions of boosted BHs is subject to ongoing investigations and will be presented elsewhere [7].
In the context of this work it is particulary interesting to point out results in the pointparticle approximation reported in [63,64]. Although covering a complementary setup to our investigations it is interesting to note that the qualitative behaviour of both methods agrees well. Specifically, the authors of [64] have shown that the gravitational radiation, emitted when a point particle without initial boost falls into a higher dimensional BH, increases with the spacetime dimensionality.
I have, furthermore, achieved succesful evolutions of unequal-mass BH binaries in higher dimensions, by solving Einstein equations in the full nonlinear regime. I have focused on head-on collisions in D = 5 spacetime dimensions. The gravitational waveforms exhibit similar features as their 4-dimensional counterparts [1,205], and I were able to estimate the ringdown frequencies of the lowest multipoles. I find good agreement with published values for the quasinormal frequencies, extracted in a linearized formalism [211]. When extrapolated to the zero-mass-ratio limit, our results agree with linearized calculations [64] at the % level or better for the energy and momentum radiated, as well as for the multipolar dependence. Our findings, supported by linearized analysis, indicate that the higher multipoles become more important for larger D. This will certainly make wave extraction at sufficiently large D a more demanding task, since higher resolutions are necessary to resolve these modes. Additionally, perturbations falloff with a power 1 r D−3 resulting in steeper gradients and shorter length scales and, thus, require finer meshes. Finally, it would be very interesting to perform an exhaustive set of simulations in higher D: our results, together with linearized analysis [64], suggest a qualitative change in radiation emission for D ≤ 12 − 13. In fact, this change is required by the fact that Hawking's area theorem forces the total amount of gravitational radiation to decrease with D, at sufficiently large D [3,64]. Understanding the mechanism at play requires extension of our results to arbitrary spacetime dimensions.
The results presented in this chapter focused on low-energy, i.e., non-boosted head-on collisions of BHs in D = 5 and D = 6 spacetime dimensions.
Of utmost importance for potential applications in high energy collider physics are fully dynamical, non-linear evolution of boosted BH collisions. Results of PP calculations [64] of a particle with initial kinetic energy varying from zero to the UV limit raise the expectation for astonishing results. Counter-intuitively, they find that particles with intermediate large boost radiate less energy than their non-boosted counterparts. Also, while the energy released by the infall of a non-boosted particle increases with spacetime dimension, the trend in the UV limit is reversed. Preliminary fully dynamical studies of head-on collisions of boosted BHs in D = 5 and D = 6 dimensions, which will be presented in detail elsewhere [7], indicate a qualitatively similar trend.
Conclusively, I have started to explore a wide range of interesting phenomena and accomplished the numerical modelling of BH collision higher dimensional spacetimes. This is but the tip of the iceberg and leaves plenty of room for exciting future investigations.
6. Black holes in a box: Toward the numerical evolution of black holes in AdS

Introduction
A powerful tool to explore strongly coupled gauge theories by investigating classical gravity (and vice versa) is provided by the gauge/gravity duality, originally proposed by Maldacena in 1998 [69]. In a nutshell, the gauge/gravity duality is a correspondence between a gravity theory in D-dimensional anti-de Sitter (AdS) spacetime (on the one hand) and a strongly coupled gauge theory "living" on its (D −1)-dimensional boundary. In many of these developments BHs play a crucial role, as may be seen by the following list of examples: (i) the successful microscopic computations of the Bekenstein-Hawking entropy for extremal BHs [268] and Hawking emission rates for near extremal BHs [269] are now seen as applications of the correspondence; (ii) the confinement/deconfinement phase transition in QCD-like theories has been identified [70] with the Hawking-Page phase transition for AdS BHs [71]; (iii) moving away from thermal equilibrium, the quasinormal frequencies of AdS BHs have been identified with the poles of retarded correlators describing the relaxation back to equilibrium of a perturbed dual field theory [270,271]; (iv) in a large class of gauge theories with a gravity dual, a universal behaviour was obtained for the ratio of the strongly coupled medium's viscosity to entropy density, by computing the absorption cross-section of low energy gravitons in the dual BH (or black brane) geometry [272]. The result is in good agreement with experimental results from the Relativistic Heavy Ion Collider (RHIC); (v) critical exponents, of the type found in spherical gravitational collapse by Choptuik in 4-dimensional asymptotically flat spacetime [273], have been conjectured to be dual to the asymptotic value of the parton saturation exponent for high energy scattering in QCD, in the Regge limit, at weak coupling [274]. These and other examples have built expectations that the gauge/gravity correspondence will be a useful laboratory for exploring and understanding very difficult problems of both field theory and gravity, such as the confinement problem, the information loss paradox and the problem of singularities. Given the potential of the correspondence exemplified above, a working framework to solve Einstein's equations exactly for a broad range of initial conditions in AdS spaces would, unquestionably, be very useful for deepening the study of the correspondence, particularly in dynamical situations.
In [76], the authors presented (for the first time) fully numerical evolutions of AdS 5 spacetimes with SO(3) symmetry and studied prompt BH formation. Non-linear evolutions of scalar fields in AdS background [74,275,276] and investigations of perturbations of AdS spacetimes [277] (to higher order in perturbation theory) suggest that AdS is subject to the non-linear turbulent instability. On the other hand, recent simulations of a wider class of real and complex scalars in AdS revealed a threshold for BH formation [278] and the evolution of time-periodic solution hint at stable islands in the phase-space [75].
It has been our long-term goal to assemble a framework for fully dynamical, numerical evolutions in AdS spacetimes and the achievements presented in this chapter serve as the first step in this direction. In order to accomplish this goal, one has to go beyond the standard methods of NR in, at least, two obvious points. First of all, AdS spacetimes are not globally hyperbolic. In asymptotically AdS spaces the boundary plays an "active role" for the bulk evolution. This is easily visualised in the Penrose diagram of AdS, which has a timelike boundary. Physically, null geodesics in AdS reach the boundary for a finite affine parameter. One thus often refers to an asymptotically AdS space as a "box", having in mind that AdS boundary conditions directly affect the bulk physics [279][280][281]. In the gauge/gravity correspondence, the choice of the AdS boundary conditions is dictated by a holographic prescription [282][283][284][285]. Secondly, from the viewpoint of the duality, D-dimensional AdS spacetimes are relevant. Thus, we intend to develop a framework that could be used in AdS D , in particular for D = 5, which is related by the correspondence to 4-dimensional gauge theories. The latter issue has been recently addressed, e.g., in Refs. [2][3][4][5]218,219,221,241,243,244,286] and in recent reviews on the topic [17][18][19][20].
Here we shall focus on the former issue: the active role of boundary conditions. The dynamics of BHs in AdS, and especially the role of spatial infinity ("the box") is poorly understood. In contrast to the asymptotically flat case, interesting new phenomena may occur in AdS backgrounds. For instance, superradiance effects have been shown to make small (as measured by the AdS radius) rotating BHs unstable, through a sequence of reflections at the boundary and amplifications close to the ergoregion [77][78][79]99,[287][288][289][290][291]. The final state of this instability could be a new non-axisymmetric BH configuration, which is also supported by recent gravity/hydrodynamics arguments [228,229]. Notice that non-axisymmetric BHs are strictly forbidden in asymptotically flat spacetimes [292,293], so the boundary does play an important role in the description of BHs.
In order to identify the influence of the boundary on the bulk evolution, we consider here a toy model for AdS. We set the cosmological constant to zero and impose mirrorlike boundary conditions on a box that contains the dynamical system. This mimics the AdS global geometry, keeping the local geometry of vacuum models. We choose the dynamical system to be a BH binary, starting at some given distance, producing either a head-on collision or an inspiralling merger. In the latter situation we consider the initial BHs without intrinsic angular momentum. These are, by now, very well tested systems when purely outgoing boundary conditions are imposed. Thus we will be able to see clearly the modifications due to the non-outgoing boundary conditions in systems with non-trivial dynamics.
For the post-merger dynamics the inspiralling binaries provide a more interesting analysis than the head-on collision case, because the initial center-of-mass energy transferred into gravitational radiation is by more than one order of magnitude larger in the former case as compared to the latter. Immediately after the merger, the system will contain a single (spinning or non-spinning) BH plus gravitational radiation. This radiation will then be (repeatedly) reflected off the boundary and interact with the BH.
The first non-trivial result is that we can follow the numerical evolution for up to two reflections off the boundary of the gravitational radiation produced in the merger. A priori it was not guaranteed that this could be achieved, since it is not known whether the formulation of the Einstein equations that we use provides a well defined initial value boundary problem together with the boundary conditions we impose. Our simple setting actually provides a first attempt to test the (numerical) stability of the initial boundary value problem in a non-globally hyperbolic spacetime. We find that our numerical results are at least second-order convergent for at least two reflections off the wall, after which we gradually lose convergence. A deeper study of these issues is clearly needed, as well as an exploration of how the convergence (and remaining results) change for different boundary conditions. During the window of numerical convergence, we study the properties of both outgoing and ingoing gravitational radiation. The usual studies of BH binaries with outgoing boundary conditions, focus only on the Weyl scalar Ψ 4 , which describes outgoing gravitational waves. However, an equally relevant quantity for the description of gravitational radiation is the scalar Ψ 0 , which encodes ingoing waves, but which is seldom discussed in the literature. Due to our special boundary conditions and setup, we are able to verify certain relations between these two quantities in a numerical evolution for the first time.
This also provides a test on the correctness and meaning of the boundary conditions we have imposed.
By analysing the properties of the apparent horizon of the BH produced in the merger and after each interaction with the gravitational wave packet, we estimate the amount of energy and angular momentum that is transferred from the radiation into the BH per interaction. In case of the inspiralling binary the boxed BH is spinning and we expect superradiant scattering of the waves generated during merger to become important; in fact, the back and forth bouncing of the waves at the reflecting wall and their subsequent amplification by superradiance close to the ergoregion are expected to turn the system into a BH bomb [78,80,99]. Thus, these simulations will be the first attempt at a nonlinear study of the BH bomb. An important open question concerns the understanding of how the evolution proceeds and the end point of the instability. This can only be achieved through non-linear studies. A final statement on this issue will, however, require further analysis than that provided herein. This chapter is organized as follows: In Sec. 6.2 we review the numerical framework with particular focus on the implementation of the "spherical" boundary and the imposed boundary conditions. The numerical results of our simulations are shown in Sec. 6.3. In Sec. 6.4 we close with a discussion of the results and prospects for the future. Appendix D.2 exhibits some snapshots for visualising the evolution of the system we have studied. Appendix D.1 describes the simulations with a cubic, rather than spherical, box. The results presented in this Chapter are based on the publication [8].

Numerical framework
In order to numerically generate a solution to the Einstein field equations, it is most convenient to view the problem as a initial value problem. Therefore we employ the 3+1decomposition of the spacetime and evolve the resulting 3 + 1-Einstein's equation in the (generalized) BSSN form together with the moving puncture approach (see Chapter 2 for a detailed discussion). We evolve the system at hand with the Lean code [161,206]. The key ingredient in which our current numerical framework differs from previous implementations of the Lean code and most other NR codes is the outer boundary condition, which we will discuss in more detail in the remainder of this section.
The vast majority of numerical simulations of BH binaries has been concerned with asymptotically flat spacetimes and consequently employed either of the following boundary treatments: (i) outgoing Sommerfeld conditions on Cartesian grids of finite size, as described for example in [127]; (ii) outgoing radiation [294,295] with multipatch methods, including Cauchy characteristic wave extraction [296,297] and (iii) constraint preserving boundary conditions combined with multidomain methods [298,299].
In contrast we will study the dynamics of BH spacetimes under the influence of a reflective outer boundary. It is natural to use for this purpose an outer boundary of spherical shape. Most importantly, this avoids mixing of different gravitational wave multipoles as would occur in the case of a reflective, cubic outer boundary. This is discussed in more detail in Appendix D.1 where we compare simulations using both types of boundary. Except for this comparison, however, we will exclusively study spherical outer boundaries or, rather, approximate these by using so-called Lego spheres; cf. Sec. 3 in [300]. In the left panel of Fig. 6.1 we sketch the foliation of the spacetime under consideration, suppressing one spatial dimension for simplicity. The numerical domain, i.e., the Lego sphere is visualized by a dark (red) domain on each timeslice Σ t+nδt . Their numerical implementation is illustrated in the right panel of Fig. 6.1 which schematically displays a computational domain using four refinement levels with one or two components each. The individual components are labelled G i m where the indices i and m denote the refinement level and component number. Note that one spatial dimension is suppressed for visualisation purposes. In order to update a grid function at a particular vertex, we require information from neighbouring points because of the discretization of spatial derivatives in the evolution equations. The exact number n of neighbouring points required in each direction depends on the finite difference stencils employed. While n = 3 for the 4 th order accurate stencils used in our simulations, we use n = 1 for simplicity in our illustration in the right panel of Fig. 6.1. Consider first the dark (blue) shaded area inside the inner solid circle of radius R B . Each point in this regular domain can be updated straightforwardly provided we also have valid data on the boundary points marked by × symbols. Points outside the circle of radius R B are not required for updating regular points and are simply ignored in the numerical evolution. The specific boundary condition is then determined by the manner in which we update grid functions on the boundary points marked as × in the figure.
In order to mimic the global structure of an Anti-de Sitter spacetime we effectively enclose the BH binary inside a spherical mirror and set at each boundary point with f denoting any of the BSSN variables. The use of fourthorder stencils adds one complication to this picture: the upgrade of a grid point requires two neighbors, so that points right next to the boundary need special treatment. In practice, we have achieved optimal stability properties by evolving these points with second-order stencils. Our implementation requires one further ingredient in order to handle the spurious radiation inherent to numerically generated initial data of BH binary systems; cf. [301]. In order to avoid contamination of our simulations by such spurious radiation being trapped inside our reflective boundary we employ standard outgoing radiation boundary conditions at early times and only switch on our reflective condition at In order to avoid a discontinuous jump from outgoing to reflective boundary conditions, we gradually switch off the time derivative ∂f /∂t using a weighting factor w(t) which smoothly decreases from 1 to 0 over an interval ∆t = 10 M and ∆t = 20 M for the head-on collision and inspiral, respectively. The duration of the spurious wave pulse ∆t pulse is estimated from previous simulations of similar setups in asymptotically flat spacetimes as for example presented in Refs. [118,161,171]. The spurious radiation is thus given sufficient time to leave the computational domain.  Table 6.1.: Grid structure, as well as initial and final parameters of the simulated black holes. Their initial coordinate distance has been set to d = 6.517 M and the irreducible mass of each BH is M irr,i = 0.483 M . The grid setup is given in terms of the radii of the individual refinement levels as well as the resolution near the punctures h (see Sec. II E in [161] for details). The reflective outer boundary is located at radius R B . The table further shows the Bowen-York [121] parameter for initial linear momentum P i of the individual holes. J AH fin and J rad fin are the spin of the single hole after merger determined from the AH and the merger radiation, respectively. All parameters are given in units of the ADM mass. We did not monitor the AH properties for all runs, therefore the spin of the final BH is not determined (empty cells) for some cases. The Weyl scalars have been extracted at R ex = 35M (IN1, HD1) and R ex = 25M (IN2), respectively. Model VIS is used in Appendix D.2 for visualisation.

Numerical results
Our numerical study focuses on two types of binary BH initial configurations; (i) head-on collisions of non-spinning BHs starting from rest and (ii) quasicircular inspiral of nonspinning holes. In the remainder of this work we label these as HD and IN simulations. The initial parameters of all our simulations as well as the structure of the computational domain and the position of the outer boundary R B are summarised in Table 6.1. In all models the initial coordinate separations between the two BH has been set to d = 6.517 M and their irreducible mass to M irr,i = 0.483 M . Unless denoted otherwise, the results presented refer to the highest resolution available.

Numerical convergence analysis
Before we discuss in detail the physical properties of the BH binary systems, we test the accuracy of our simulations by performing a convergence analysis of model IN2. Numerical simulations based on the finite differencing method typically approximate the continuum solution of differential equations with a leading error that has a polynomial dependence on the resolution, f cont = f num +O(h n ). The order n depends on the specific numerical implementation. Consistency of the code is tested by evolving the same configurations with low, medium and high resolution h c , h m and h f . One straightforwardly shows that the convergence factor is then given by where f h stands for any of the evolved variables obtained for resolution h. For the case of contemporary moving puncture codes, the solution is complicated by the fact that the differential equations are typically discretized with fourth-(or higher) order accurate stencils but prolongation in time between different refinement levels and implementation of outer boundary conditions is only second-order accurate; see Sec. IV in Ref. [302] and our discussion in Sec. 6.2. Here the differences f hm − f h f have been amplified by the convergence factors Q 2 = 1.26 and Q 4 = 1.47 expected for second and fourth-order convergence, respectively. The figure indicates fourth-order accuracy for the first passage of the wave pulse and a gradual deterioration of convergence to second-order accuracy afterwards. We believe this to be a consequence of the different ingredients of the code as discussed above. At early stages, the dominant error is the discretization of derivatives. As the pulse successively passes across mesh refinement boundaries and is reflected off the outer boundary, however, the second-order error in the prolongation operation becomes dominant and reduces the order of convergence. We also note, in this context, that well-posedness of the BSSN evolution system with reflective boundary condition has so far not been demonstrated 1 . We can therefore not rule out adverse effects on the long-term convergence properties due to potential ill-posedness of the continuum system of equations. In the remainder of this discussion we will restrict ourselves to 2-3 passages of the wave pulse as covered in Fig. 6.2 during which the relative uncertainties in Ψ 4 and Ψ 0 are ≤ 5%.

Gravitational wave signal and black hole dynamics
To our knowledge, this work presents the first analysis of gravitational waveforms with both outgoing (Ψ 4 ) and ingoing (Ψ 0 ) contributions for long-term stable numerical simulations of BH binaries. For this reason, we first illustrate the general pattern of the wave signal obtained for model VIS of Table 6 The gravitational wave signal is dominated by the quadrupole contributions which we show in Fig. 6.3. For clarity, the ingoing signal ψ 0 22 has been shifted in time by ∆t = 10 M in order to compensate for the additional propagation time from the extraction radius R ex = 35 M to the boundary R B = 40 M and back after reflection. The reflection introduces an additional phase shift of ∆φ = π which has also been taken into account We show the differences of the l = m = 2 mode between the coarse and medium and the medium and fine resolution run. The latter has been amplified by the factors Q 4 = 1.47 and Q 2 = 1.26. We observe fourth-order convergence in the signal due to the merger whereas the first and second after-merger cycles show only second-order convergence. The first two reflected and ingoing wave pulses show second-order convergence. in the figure. Within numerical errors, we find the resulting outgoing and subsequent ingoing pulses to overlap.
The first outgoing wave pulse, visible in Fig. 6.3 around 150 ≤ t/M ≤ 250, is generated during the inspiral, plunge and merger of the binary and is similar to waveforms obtained for the inspiral of non-spinning BH binaries in asymptotically flat spacetimes (cf. Fig. 1 in Refs. [161,302]). Due to the reflecting boundary, however, this wave pulse does not escape the computational domain. Instead it propagates inwards, interacts with the post-merger remnant hole and eventually manifests itself as a second wave pulse shifted by ∆t ≈ 80 M relative to the first. This process repeats itself many times, with the wave pulse being presumably distorted (by absorption, superradiance and other curved spacetime effects on wave propagation) upon each interaction with the BH. We now investigate in detail these changes of the wave pulse upon interaction with the BH.
Interaction of the wave pulse with the remnant black hole As shown in Fig. 6.3, the outgoing and subsequent ingoing wave pulses overlap within numerical uncertainties. We therefore focus on the outgoing signal in our study of subsequent wave pulses and the gradual changes caused by successive scattering off the BH. Changes in the wave pulse are best illustrated by considering the wave amplitude as shown in Fig. 6.4. Here we superpose the l = 2, m = 0 mode for model HD1 and the l = 2, m = 2 multipoles for models IN1 and IN2 of the first three successive outgoing wave pulses by applying corresponding time shifts to the waveform. Clearly, the wave pulses broaden after each scattering off the BH. We emphasize that this distortion of the pulse is not an artifact of the outer boundary condition as is demonstrated by the good overlap between the ingoing and outgoing pulses in Fig. 6.3.
One possible explanation for this distortion relies on strong field, curved spacetime effects: massless waves in curved geometries do not propagate on the light-cone. Indeed, as shown in the classical work by DeWitt and Brehme, the Green's function for a massless field in a curved spacetime does not generally vanish inside the light cone [260]. This gives 3. This spectrum corresponds solely to the initial wave packet, i.e., the waveform has been truncated immediately before the first reflection off the boundary. The vertical line marks the threshold frequency for superradiance. rise to interesting effects. It is well-known that wave propagation in BH spacetimes shows that the signal can roughly be divided in three parts: (i) the first part is the prompt response, at very early times, whose form depends strongly on the initial conditions. This is the most intuitive phase, being a counterpart of the light cone propagation in flat spacetime; (ii) at intermediate times the signal is dominated by an exponentially decaying ringing phase, and corresponds to the excitation of the BH's characteristic modes of vibration [211]; (iii) a late-time tail, usually a power law falloff of the field [262,263]. Therefore, a variety of possible distortions are possible.
Another possible explanation for the increasing width of the pulse arises in the context of the superradiance mechanism. It is well-known that the scattering of a wave pulse off a rotating Kerr BH will result in superradiant scattering -amplification of the scattered wave packet via extraction of rotational energy from the BH -if the wave pulse satisfies ω <mΩ , (6.4) where Ω ≡ j fin /(2r + ) is the BH angular velocity [88,305,306]. We note, however, that Fig. 6.4 also indicates a broadening of scattered wave pulses in the head-on case where no superradiance is expected. While the observed broadening would be compatible with superradiance, other effects appear to also influence the shape of the pulse and our observations do not conclusively demonstrate its presence. In order to investigate this dispersion in more detail, we plot in Fig. 6.5 the energy spectrum for the dominant l = m = 2 mode. The vertical line in this figure denotes the threshold frequency mΩ ≈ 0.4/M corresponding to a final spin j fin = 0.69 as obtained for the post-merger hole for configuration IN2; cf. Sec. 2.5.2. The figure demonstrates that the l = m = 2 mode does contain contributions which would be subject to superradianceinduced amplification. These results then suggest that the low-frequency component of the pulse is amplified due to superradiance, while the high-frequency component is absorbed. A linear analysis of superradiance in the Kerr geometry [306] shows that superradiant effects are always small, unless the hole is rotating close to the extremal value. Thus, further studies, including larger spins of the post-merger hole, are necessary to comprehensively demonstrate superradiant wave amplification. Black hole dynamics In contrast to the case of asymptotically flat spacetimes, a BH binary under the influence of a reflective shell does not settle down into a stationary configuration soon after merger. This is due to the repeated interaction with the wave pulse passing back and forth across the spatially finite spacetime. The prolonged dynamical state of the system manifests itself prominently in the area A AH of the apparent horizon. In Fig. 6.6 we show the fractional deviation (A AH − A 0 )/A 0 of the apparent horizon area from its value A 0 immediately after merger, i.e. the first instance a common apparent horizon is found. As expected, the horizon area remains nearly constant for the duration of the first passage of the pulse to the reflective shell and back, ∆t ≈ 80 M and 60 M , respectively, for models IN1 and IN2.1 2 . The subsequent increase in A AH demonstrates that some fraction of the gravitational wave pulse energy is absorbed by the hole. It follows another period of approximately constant horizon area, a further increase upon the second scattering of the pulse and so on. Unfortunately, the radiation efficiency is almost 2 orders of magnitude lower for head-on collisions (see Table II in Ref. [161]), so that the increase in horizon area is buried in the numerical uncertainties. The head-on case serves as a useful comparison, however, as it demonstrates that the changes observed for the inspiralling configurations are significant relative to numerical uncertainties.
The BH mass, defined in terms of the equatorial radius of the horizon C e by [170] M = C e 4π , (6.5) shows a similar behaviour as the horizon area. In Fig. 6.7 we plot the fractional deviation (M − M 0 )/M 0 of the mass from its value immediately after merger together with the irreducible mass and the BH spin J of the hole obtained for model IN1. The mass remains approximately constant until the pulse returns after its first reflection, then increases, remains constant during the second passage of the pulse and so on. In contrast, the spin shows a significant increase only during the first scattering of the pulse off the BH. We conclude that in each interaction with the gravitational radiation, the final BH mass increases. It is interesting to compare the increase in the horizon mass with the amount of gravitational wave energy radiated during the last stages of the inspiral, plunge and merger of a corresponding binary system in an asymptotically flat spacetime which is about 3.5 % of the total energy of the system [161,171]. For the IN1 run, we estimate that about 15% of the energy emitted during the merger is absorbed by the central spinning BH per interaction. Our results are consistent with total en- ergy conservation. Moreover, they are not incompatible with superradiant amplification: typically, absorption of high-frequency waves is more effective than superradiant amplification of low-frequency waves, such that the net effect leads typically to absorption by the BH [306]. The prolonged non-stationary character of the post-merger state is also demonstrated by the time dependence of the BH's final spin. Immediately after merger, we obtain j fin = 0.69 from Eq. (2.88) in excellent agreement with corresponding simulations using outgoing radiation boundary conditions [171]. Successive interaction with the reflected wave pulse, however, results in a small but significant increase in the BH's spin as shown in Fig. 6.7. As before, the first increase occurs about ∆t = 80 M after merger, when the pulse has returned to the BH. We estimate the fractional increase in spin resulting from the first scattering at about 5 %. For comparison, the total angular momentum radiated in the case of an asymptotically flat spacetime is reported as J rad /M 2 = 0.246 (28 % of the initial orbital angular momentum of the system) in Table  I of Ref. [171]. It thus appears that a significant amount of angular momentum remains in the form of gravitational waves. Due to numerical uncertainties it is not entirely clear whether later periods of interaction between pulse and hole result in a further transfer of angular momentum from the wave pulse to the hole or vice-versa. Our results indicate, however, that the amount of angular momentum exchanged in subsequent interactions is significantly below 5 %.

Discussion
The dynamics of BHs in generic spacetimes is a fascinating, yet extremely challenging problem. The gauge/gravity duality, however, strongly motivates us to solve dynamical problems with BHs on asymptotically AdS backgrounds. In this work, we have studied a toy model that captures one of the fundamental features of such backgrounds: the active role played by the boundary conditions for the bulk evolution.
We have mimicked the global structure of an AdS background by introducing a reflecting wall at some radius. Within this cavity we evolved an inspiralling BH binary and a BH binary starting from rest at a certain initial distance. Of course, these are very specific initial configurations, and serve merely as tests for future, possibly more complex, situations.
Perhaps the most important conclusion of the present work is that these simulations can be done and represent the first step to a full numerical evolution of BHs in AdS spacetimes. Indeed, as observed in Section 6.3.1, it is not known whether the BSSN evolution scheme together with reflecting boundary conditions is a well-posed initial boundary value problem. Thus, the convergence we have exhibited, which holds up to two reflections off the boundary of the gravitational radiation produced in the merger, is the first of our results. Among the other results presented here we stress the following: 1. For the first time, we were able to numerically study the scalar Ψ 0 , describing ingoing waves, and check in the numerical data the simple relations between Ψ 0 and Ψ 4 (cf. Eqs. (2.75), (2.76), (2.77) and Fig. 6.3).
2. Our results are consistent with the intuitive expectations for a wave packet of radiation (generated during inspiral plus merger) travelling back and forth between the mirrorlike wall and the BH: part of this radiation is absorbed when interacting with the BH (especially high-frequencies). We estimate that about 15% of the wave packet's energy is absorbed by the BH per interaction, at least during the first cycles.
3. The wave packet is clearly distorted upon interaction with the hole, which means some frequencies are absorbed more efficiently than others. In principle, we should observe a small but non-zero superradiance. Unfortunately, we have not obtained incontrovertible proof of non-linear superradiance from our numerical data. Presumably, the system will become unstable after a sufficiently long time, since the radiation should be exponentially amplified, once the high-frequency components had time to be completely absorbed by the BH. In future work we plan to investigate these instability studies further by considering a highly spinning, final BH produced by the inspiral of spinning BHs.
One issue that we have not explored in this paper is the potential influence of the box on the premerger dynamics. For sufficiently small size boxes, it is plausible that the radiation produced in the inspiral may be reflected off the boundary and interact with the binary before the merger. This might produce observable signatures in the premerger dynamics and even in the properties of the remnant black hole. Whereas we seem to observe some hints of this effect in our numerical data, a more exhaustive analysis is required to produce some precise statements.

Introduction
Nowadays, BH physics has grown into a mature and fully developed branch of GR. The theory has been worked out many decades ago, and tremendous progress has been made in actually observing some of the fascinating general relativistic effects. From X-spectra on the inner edge of accretion disks, which probe the innermost stable circular orbit of the geometry, to gravitational wave physics, "precision BH physics" is a new and rapidly developing field [24,25,91,92]. One of the fundamental reasons why precision BH physics is possible at all, are the no-hair and uniqueness theorems: BHs in 4-dimensional, asymptotically flat spacetimes must belong to the Kerr-Newman family and are, thus, fully specified by three parameters only: their mass, angular momentum and electric charge (see e.g. Ref. [40], or Carter's contribution to Ref. [307]). In more colloquial terms, this is commonly expressed by saying that BHs have no hair or, rather, have three hairs only. This simple yet powerful result has far reaching consequences: Given some arbitrary perturbations with the same conserved charges, they must all decay to the same final state, namely one BH with those charges. By now, there are a plethora of studies, at the perturbative and fully non-linear level, investigating the approach to the final stationary state (see, e.g., Ref. [211] and references therein), which we summarize below. The present chapter is part of a larger publication [9] that appeared simultaneously with an accompanying work by Dolan [101], exploring the dynamics of complex, massive scalar fields. Additionally, in Ref. [9] we have extended our study to investigate the time evolution of massive vector fields, going up to very high spins (close to extremality) of the background Kerr BH. A complentary study in the slow-rotation approximation has been presented by Pani et al [308,309]. More recent reviews on the BH bomb mechanism and BH-scalar field systems can be found in Refs. [310,311].
Generic response of a BH spacetime to external perturbations The generic behavior of massless fields around a BH can be divided into three parts: (i) A prompt response, depending on the initial data, at early times, which is the counterpart to light-cone propagation in flat space; (ii) An exponentially decaying "ringdown" phase at intermediate times, where the BH is ringing with its characteristic quasinormal modes (QNMs). This stage typically dominates the signal, and its properties, such as vibration frequency and decay timescale, depend solely on the parameters of the final BH [211]. The long-lasting oscillation of the lowest QNMs is the most important stage in the life of any field around a BH. Its lifetime, or quality factor, depends solely on the BH spin [211]. Specifically, the lifetime tends to increase with growing spin and the decay timescale approaches zero for nearly extremal BHs. Because of the no-hair theorem, the detection of QNMs allows to uniquely determine the BH charges and provides tests of GR [211,256,312]; (iii) At late times, the signal is dominated by a power-law fall-off, dubted "late-time tail" [261,262,313]. Tails are caused by backscattering off spacetime curvature and, as such, also appear in other situations where light propagation is not on the light cone [263].
Superradiant effects Among the most exciting phenomena is the superradiant effect [81,82,314]: In a scattering experiment off a BH the scattered wave will be amplified if its complex frequency ω = ω R + ıω I satisfies the superradiant condition where m is the azimuthal wave quantum number and Ω H the angular velocity of the BH horizon. Then, the excess energy is withdrawn from the object's rotational energy [81,82] and, in a dynamical scenario, the BH would spin down. The effect can be attributed to the existence of negative-energy states in the ergo-region, and dissipation at the event horizon. Superradiance is the chief cause of a number of exciting phenomena in BH physics: (i) Generic perturbations are damped away to infinity and across the event horizon. Because rotating BHs amplify waves that fulfill the superradiant condition (7.1), the amplification factors as well as the quality factor, or lifetime, of these superradiant modes increase with rotation. (ii) Typically, any satellite around a BH spirals inwards as time goes by, due to gravitational wave emission and energy conservation. Emission of radiation at infinity results in a larger binding energy of the particle. Because superradiance implies the extraction of the BH's rotational energy, it is possible that the energy deficit comes entirely from the BH kinetic energy. In this way, satellites around rapidly spinning BHs can in principle orbit on a fixed radius -on so-called floating orbits -for a much longer time, tapping the BH's kinetic energy. In BH binaries, this effect can dominate in the presence of resonances [94,95,315]. Notice, that this phenomenon is analogous to tidal acceleration, e.g., in case of the Earth-moon system [316,317]. (iii) A further interesting effect enters the game if we enclose the spinning BH inside a perfectly reflecting cavity. As was recognized already by Zel'dovich [81,82], this way one would create an unstable system. Any initial perturbation will get successively amplified near the BH and reflected back at the mirror, thus creating an instability, which was termed the "BH bomb" [80,99]. Whereas the setup appears physically artifical at first glance, in Nature the role of the mirror can actually be realized by anti-de Sitter spatial infinity or massive fields. In the former case, the BH bomb translates into a real, physical instability of (small) rotating BHs in asymptotically AdS spacetimes [78,79,289,291]. (iv) Finally, of direct interest for the present study is the fact that massive scalar fields around Kerr BHs are prone to a BH bomb-like instability, because the mass term effectively confines the field [83][84][85][86][87][88]318]. The instability is described by the time dependence of the field, Ψ ∼ e −iωt with complex frequency ω For small coupling M µ S 1 between the BH mass M and scalar field mass µ S , the characteristic frequency giving rise to the instability is [308] where the BH is spinning at angular momentum J = aM . In the opposite limit, i.e., for very large mass couplings M µ S 1, the instability time-scale is [84] 3) The maximum (scalar field) growth rate of M ω I ∼ 1.5 · 10 −7 has been found for the dipole mode and the mass coupling M µ S = 0.42 in the background of a Kerr BH with a/M = 0.99 [87,88]. These calculations have been performed in the linear regime, thus neglecting back-reaction effects such as the BH spin-down or effects due to nonlinear self-interaction of the scalar field. Therefore, the final state of the superradiant instability in the fully dynamical regime is not known, partly because it requires the nonlinear evolution of Einstein's equations for a timescale of order 10 6 M (in geometric units G = c = 1). A plausible evolution scenario consists on an exponentially growing scalar condensate outside the BH, extracting energy and angular momentum from the BH until the superradiant extraction stops, i.e., until condition (7.1) is no longer satisfied. Further interesting new phenomena arise when we consider non-linear interaction terms such as the bosenova-type collapse recently reported in [93,97,98].

Superradiant instability in astrophysical systems
We have argued beforehand that massive fields in the vicinity of BHs give rise to BH bomb-like, superradiant instabilities. However, the effect is very weak for known standard model particles in astrophysical environments: For example, the mass coupling for the lightest known elementary scalar particle, the pion, around a solar mass BH is M µ S ∼ 10 18 , resulting in an instability timescale much larger than the age of the universe. Nevertheless, the superradiant instability around astrophysical BHs might become significant if there exist fields with tiny, but non-vanishing mass. One exciting possibility is provided by axions, ultralight bosonic states emerging from string-theory compactifications, which have not been ruled out by current experiments. In the "axiverse" scenario the existence of an entire landscape of ultra-light pseudo-scalar fields has been proposed, covering a mass range from 10 −33 eV ≤ µ S ≤ 10 −8 eV (see [91][92][93] for recent overviews). In fact, the existence of ultra-light axions leads to a plethora of possible observational implications and signatures, such as modifications of the cosmic-microwave background polarization (for 10 −33 eV ≤ µ S ≤ 10 −28 eV ). They are also anticipated to make up a fraction of dark matter if 10 −28 eV ≤ µ S ≤ 10 −22 eV . Of particular interest in the context of BH physics are axions in the mass range 10 −22 eV ≤ µ S ≤ 10 −10 eV [91][92][93]. Then, the time scales for the superradiant instability becomes astrophysically significant, giving rise to a number of interesting effects: (i) A bosonic cloud bounded in the vicinity of a Kerr BH might create a "gravitational atom", which can be de-excited by the emission of gravitons, thus carrying away BH angular momentum; (ii) If the accretion of bosons from this cloud is efficient, the rotation of the BH can be sustained and it might be turned into a GW pulsar; (iii) If, on the other hand, the accretion from the axionic cloud is not efficient enough, the BH will eventually spin down, thus yielding gaps in the Regge plane (the phase-space spanned by mass and spin parameter of the BH). Further possible effects have been discussed in Refs. [91-98, 310, 311].
Similar superradiant instabilities are expected to occur for massive hidden U (1) vector fields, which are also a generic feature of extensions of the standard model [319][320][321][322]. While superradiant instabilities have been widely studied for massive scalar fields [80,83,87,88,94,99,100], the case of massive vector fields is still uncharted territory, the main reason being that the Proca equations seem to be non-separable in the Kerr background. After the defence of this PhD thesis, investigations in the slow-rotation limit [308,309] and full-blown numerical studies [9] have shed more light on the phase-space and time progression of Proca fields in Kerr. Furthermore, so far most studies on the massive boson instability are done in Fourier space. As far as we are aware there has only been one attempt, besides the most recent studies in [97,101], at studying the massive scalar field instability in the time domain, with generic initial conditions: the work by Strafuss and Khanna [102]. We believe that some of the conclusions reported in that work are wrong; specifically the reported instability growth rate of M ω I ∼ 2 · 10 −5 , which is two orders of magnitude larger than previous analytic results [87,88] and more recent numerical studies [97]. Here, we will attempt a correct explanation for these puzzling results.
In this chapter we will present the evolution of massive scalar fields in Schwarzschild and Kerr BH background in the time-domain, which is part of the publication [9]. We will provide generic Gaussian wave packets as well as quasi-bound states as initial configurations. Investigations of massive vector fields around spinning BHs, which are expected to have a larger amplification factor than massive scalar fields will be reported elsewhere [9]. The exploration of the dynamical coupling of the massive scalar to the BH spacetime as well as non-linear couplings and self-interaction terms are work in progress and subject to future investigations. This chapter is organized as follows: In Sec. 7.2 we present the numerical framework, describing the formulation of the problem as Cauchy problem, the setup of initial configurations and the background spacetime as well as extraction of the scalar field. Sec. 7.3 is devoted to the numerical results of the (massive) scalar field evolution. In particular, we present a number of benchmark tests to verify our implementation before studying more generic setups. We will show that the evolution of a massive scalar has a non-trivial pattern, which can be explained in terms of multi-mode excitation. We believe that this pattern also describes the results reported by Strafuss and Khanna [102]. Finally, we will finish this chapter with a summary of our results and some concluding remarks in Sec. 7.4.

Numerical framework
I intend to explore the stability properties of Kerr BHs against massive scalar field perturbations and track their evolution in this background in the time domain. In the following I will describe the formalism and numerical framework that I have developed in order to achieve this goal.

Action and equations of motion
I focus on the generic action involving massive scalar Ψ and vector fields A µ with mass µ S and µ V , respectively, where k axion parametrizes the axion-like coupling, F µν ≡ ∇ µ A ν − ∇ ν A µ is the Maxwell tensor and * F µν ≡ 1 2 µνρσ F ρσ is its dual. Also µνρσ ≡ 1 √ −g E µνρσ and E µνρσ is the totally anti-symmetric Levi-Civita symbol with E 0123 = 1. 2 The resulting equations of motion are I restrict myselve to vanishing axion-like coupling k axion = 0 and investigate massive scalar fields only. Furthermore, I focus on weak fields throughout this study, whereas non-linear effects and back-reaction on the metric will be studied elsewhere. At linear order in the scalar field amplitude, the equation of motion for the metric, Eq. (7.5c), is identical to Einstein's equations in vacuum. Therefore, I consider a Kerr background and focus on Eqs. (7.5a) in the remainder of the chapter. In Boyer-Lindquist coordinates, the metric is described by where Σ =r 2 + a 2 cos 2 θ, ∆ = r 2 − 2M r + a 2 .
This geometry describes a rotating BH with mass M and angular momentum J = aM 2 . Note, that in order to ensure the regularity of the the spacetime, i.e. the existence of an event horizon, the BH spin is constrained by the Kerr bound a/M ≤ 1. I will set M = 1 throughout our study. Because I intend to solve the equations of motion, Eqs. (7.5a), numerically it is convenient to reformulate them as time evolution problem. For this purpose I employ the 3 + 1-decomposition of the spacetime, discussed in Chapter 2 for the case of dynamical BH spacetimes.
Evolution equations for scalar fields I focus on the setup for a massive scalar field without self-interaction, which implies k axion = V = 0. Its dynamics are determined by the Klein-Gordon equation (7.5a). I write the equation of motions as a time evolution problem, i.e., as a set of first order in time, second order in space partial differential equations (PDEs). Therefore, I introduce the conjugated momentum The definition (7.8) provides evolution equations for the scalar field Ψ where L β Ψ = β k ∂ k Ψ. By applying the 3 + 1-split, the Klein-Gordon equation yields the evolution equations for the momentum where L β Π = β k ∂ k Π and I consider the background spacetime in 3 + 1-form given by Eq. (2.4).

Background in horizon penetrating coordinates
For the purposes of our numerical implementation, I consider the Kerr spacetime in Cartesian Kerr-Schild coordinates (t, x, y, z) as background spacetime. In particular, I employ the spacetime metric (7.6) in 3 + 1-form given by Eq. (2.4). Then, considering a Kerr BH with angular momentum aM in the z direction, the lapse function α, shift vector β i , 3-metric γ ij and extrinsic curvature K ij in Kerr-Schild coordinates (t, x, y, z) are given by The Boyer-Lindquist radial coordinate r BL is indirectly defined via For our computations I will consider the root where I use the notation R 2 = x 2 + y 2 + z 2 − a 2 . Inserting the relations into Eqs.(7.11), the spacetime is described explicitly by

Initial data
I will focus on generic-type initial configurations of Gaussian form, which gives rise to many interesting effects as I will see in Sec. 7.3. Additionally, as a code test and a cleaner way to search for instabilities, I have implemented quasi-bound states as initial data. I describe each of these below.
Gaussian initial data In order to set up a generic Gaussian wave packet I choose initial data of the form where r 0 is the center of the Gaussian and w is its width. 0 Σ(θ, φ) is a superposition of spherical harmonics Y lm (θ, φ). Because our evolutions are carried out on a Cartesian grid, I initialize the spherical harmonics not in spherical coordinates (r, θ, φ) but Cartesian coordinates (x, y, z) x =r sin θ cos φ , y = r sin θ sin φ , z = r cos θ , (7.17) where r 2 = x 2 +y 2 +z 2 is the Kerr-Schild radial coordinate. I list the spherical harmonics up to l = 2 in Cartesian coordinates in Appendix B.
Quasi-bound state initial data I consider a second, useful type of initial setup, which I term quasi-bound state initial data. Perturbations of massive scalar fields around Schwarzschild or Kerr BHs give rise to long-lived modes, or quasi-bound states, which have been investigated in the frequency domain in the past [87,88,211,323]. These modes are interesting for our purposes because they prescribe a clean state, while other spurious effects are absent. For instance, a quasi-bound state consists of almost a single pure frequency, which can be superradiant. Because other modes are absent, no interference and beating of the kind described later will occur. Finally, the evolution of a pure mode serves as an additional benchmark test of the code [97].
There is a powerful and simple method to compute these modes, either by direct numerical integration or via a continued fraction approach. In order to construct the bound states for massive scalars I follow Leaver's method [87,88,211,313]. In Boyer-Lindquist coordinates, denoted by the subscript BL, the field reads 20) and ω c = mΩ H is the critical frequency for superradiance. All the remaining terms in this expression are known in closed form and the characteristic frequency ω can be obtained by solving a three-term recurrence relation [87,211]. The coefficients a n are given by, e.g., Eqs. (35)-(48) of [87]. Because this solution is given in Brill-Lindquist coordinates I perform a coordinate transformation to Kerr-Schild coordinates (see e.g. [106]), for clarity denoted by the subscript KS, Then, the bound state scalar field is given by (note that I drop the subscript KS in the following) The corresponding momenta Π lm are computed from Eq. (7.8).

Wave extraction and output
There are several different quantities that one can evaluate and discuss. For scalar fields, I will measure the fields Ψ and Π on spheres with different discrete radius r ex . Therefore I interpolate these functions on a spherical grid at each time step and decompose them into multipoles, with r = r ex fixed. Specifically,

Numerical implementation
It is our goal to investigate the dynamics of massive scalar fields in the background of Kerr BHs in the time domain. For this purpose I have developed the time evolution code Lin-Lean which is based on the Cactus computational toolkit [207,208] and the Carpet mesh refinement package [209,210]. The evolution will be performed using the 3 + 1 approach as explained in the previous sections and by employing the method of lines (see e.g. [106]). Specifically, the evolution equations are integrated in time using the fourth order Runge-Kutta scheme and spatial derivatives are discretized by second or fourth order centered finite difference (FD) stencils. Advection derivatives are discretized by lopsided FDs. Because I are considering a non-dynamical BH background, I prescribe the spacetime in Kerr-Schild coordinates and excise the region inside the event horizon r S = 2M . In practice, I specify this region by a parameter r exc ≤ 2M , that I typically choose to be r exc = 1.0M . Inside the excised region, i.e., r ≤ r exc , I enforce flat spacetime, i.e., α = 1, β i = 0, γ ij = δ ij and K ij = 0. I realize the excision either by employing second order, one-sided, outwards pointing FD stencils or extrapolate the solution in the excision region. For all practical applications, I choose the latter method and use second order FD stencils on the refinement level containing the excision region, whereas I use fourth order FD stencils on all outer refinement levels.

Evolutions of massive scalar fields
Here, we report on our investigation of massive scalar fields in Schwarzschild and Kerr backgrounds. The features of the evolved scalar field can be understood in terms of the proper oscillation modes of the BH in conjuction with late-time tails arising from curvature backscattering. Whereas the massless scalar field clearly shows this quasinormal ringdown signal followed by the expected power-law tail at late times, massive fields exhibit a much richer structure. The analysis of our results will focus mainly on the dipole mode, i.e., l = m = 1, although other modes are also presented. In Table 7.1 we summarize the frequencies of the two lowest lying dipole modes for a number of mass parameters considering either a Schwarzschild or highly rotating Kerr BH with spin parameters a/M = 0.0 or a/M = 0.99, respectively [211,256]. In case of massless  perturbations, the definition of "lowest lying", or fundamental, mode corresponds to the least damped mode. Massless perturbations are always short-lived, unless the BH is nearly extremal. In contrast, massive perturbation have a smaller decay rate. Additionally, if the mass term is present, another family of modes appears, which lives on the effective potential well, and which we call quasi-bound states [87]. These are long-lived modes and we order them by decreasing imaginary part: the fundamental mode is the shortest-lived mode in this family. At very late times, after the relaxation of the BH, there is still scattering off the background curvature present. This gives rise to late-time power-law tails [261,262,313]. In case of massless perturbations they behave as Instead, tails of massive perturbations exhibit two different stages, which strongly depent on the mass term M µ S [324][325][326] Ψ ∼t −l−3/2 sin(µ S t) , at intermediate times , (7.25a) Ψ ∼t −5/6 sin(µ S t) , at very late times . (7.25b)

Code test -space dependent mass coupling
As an initial code test, we take an unphysical space-dependent mass coupling in the background of a Schwarzschild BH. In particular, we have chosen µ 2 S = −10M 2 /r 4 . This unphysical mass coupling quickly leads to an instability of the system, which provides a unique and fast test of the code. A mode analysis of the Klein-Gordon equation is straightforward and shows the existence of at least one unstable mode, Ψ ∼ e 0.071565t . We have evolved a spherically symmetric Gaussian profile centered at r 0 = 12M and with width w = 2M . The grid setup is given by {(1024, 512, 256, 128, 64, 32, 16, 8, 4, 2), h = M/40}, using the notation of Sec. II E in [161]. The results are shown in Fig. 7.1, for which the scalar field is extracted as a function of time at r ex = 10M . Our numerical results are consistent with an exponential growth, Ψ ∼ e 0.07161t , which agrees to within 0.06% with the mode analysis.

Code test -massless scalar fields
The main features resulting from the time evolution of massless fields around BHs have been known for decades (see e.g. [211] and references therein for an overview). Thus, they provide an excellent benchmark test for our code. For this purpose we have evolved massless scalar fields around a non-rotating and highly spinning BH with a/M = 0.99. The field has been initiated as a Gaussian with width w = 2M located at r 0 = 12M . The specifications of the setup are summarized in Table 7.2. In Fig. 7.2 we display the evolution of the l = m = 0 and l = m = 1 modes of the scalar field in the respective backgrounds. The waveform exhibits all main features discussed in the Introduction: a prompt signal, followed by an exponentially decaying sinusoid and, at very late times, a power-law tail. A fit to the ringdown phase yields numerical QNM frequencies, which are summarized in Table 7.3. They are in excellent agreement within less than 2% with tabulated values [211,256]. The agreement increases for higher resolution. The powerlaw tail is of the form t p , with p = −3.08 (a/M = 0) and p = −3.07 (a/M = 0.99) for the l = m = 0 mode, in agreement with a low frequency expansion of the wave equation which yields p = −3, see Eq. (7.24). Additionally, we have performed a convergence analysis for the more challenging case of a highly rotating BH background with a/M = Run a/M 0 Σ(θ, φ) Grid Setup S mu0 a0 1 0.0 Y 00,1 −1,11 {(384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/100} S mu0 a0 2 0.0 Y 00,1−1,11 { (192, 96, 48, 24, 12, 6, 3, 1.5), h = M/144} S mu0 a99 c 0.99 Y 00,1 −1,11 {(384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/100} S mu0 a99 m 0.99 Y 00,1 −1,11 {(384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/104} S mu0 a99 f 0.99 Y 00,1 −1,11 {(384, 192, 96, 48, 24, 12, 6, 3, 1.5), h = M/108} Table 7.2.: Initial setup for simulations of a massless scalar field with Gaussian initial data located at r 0 = 12 M and with width w = 2 M in Schwarzschild or Kerr background with a/M = 0.99. We denote the modes of the initial pulse 0 Σ(θ, φ) and the specific grid setup, in units of the BH mass M , following the notation of Sec. II E in [161].  Fig. 7.3 we plot the differences of the l = m = 1 mode between the coarse and medium resolution together with the difference between the medium and high resolution runs. The latter has been rescaled by the factor Q 3 = 1.16 demonstrating third order convergence. This analysis allows us to estimate the discretization error in the signal which is about ∆Ψ 11 /Ψ 11 ≤ 8% for the l = m = 1 mode at late times of the evolution and about ∆Ψ 00 /Ψ 00 ≤ 3% for l = m = 0. The performed checks confirm that our implementation works well and yields consistent results.

Massive scalar fields: quasi-bound states
The remainder of this section is devoted to the investigation of massive scalar fields. A mass term introduces a new scale in the problem and, thus, raises the expectation that new features appear in the evolution. One of these features is the existence of long-lived, quasi-bound states, whose form and construction have been discussed in Sec. 7.2.3. We have constructed quasi-bound states for fields with mass coupling M µ S = 0.42, following [87,88,211], and evolved them in Schwarzschild or Kerr background with a/M = 0.99 according to Eqs. 7.8 and 7.22. We show the initial data, specifically |Ψ| 2 along the  x-axis, in Fig. 7.4 for the fundamental and first overtone mode (see Table 7.1). In case of non-rotating BHs the fundamental quasi-bound state is localized around the origin, but moves outwards when rotation is included. Specifically, for a highly spinning BH with a/M = 0.99, the fundamental mode is peaked at around r ∼ 12 M . The first overtone on the other hand has a pronounced minimum or node around r node ∼ 22.5 M and r node ∼ 26.5 M for a/M = 0 and a/M = 0.99, respectively. This information will prove useful when trying to understand the interplay between different modes during the time evolution of massive scalar fields initially prescribed as generic Gaussian wave packet (see Sec. 7.3.4 below).
Additionally, we have evolved the l = m = 1 fundamental bound state mode of a scalar field with M µ S = 0.42 in Kerr background a/M = 0.99. The results are shown in the left panel of Fig. 7.5, where we plot the variation of the absolute value of Ψ 11 relative to its inital value at t = 0 M . By construction, the field should remain localized in the vicinity of the BH and its absolute value Ψ 11 Ψ * 11 ∼ exp(−ıωt) exp(ıωt) ∼ const. should remain almost constant in time, with a small growth rate of M ω I ∼ 1.5 · 10 −7 [87,88]. We have tested these properties numerically by evolving the field in time and extracting its l = m = 1 mode at radii r ex = 10, 20, 30, 40M . The absolute value of the scalar field varies by less than 5% (8% for r ex = 40M ) until t ∼ 200M and by less than 1% at late times. Variations of this order of magnitude are expected for two reasons: first, the dependence of the quasi-bound state was imposed to be a spherical harmonic, see (7.22). Secondly, finite resolution effects due to the numerical scheme play a crucial role and have been identified as a further cause for the slight discrepancy. This becomes evident in the right panel of Fig. 7.5, in which we plot the relative absolute value of the scalar for the low and medium resolution run, h = M/60 and h = M/72, extracted at r ex = 20M . In fact, the two curves agree very well until t ∼ 100 M , but differ by about 2% at late times.
The quasi-bound states shown here should be unstable states with an instability time scale of ∼ 10 7 M , as shown in Table 7.1. Naturally, with such a large instability timescale, such growth is completely buried in numerical noise.

Massive scalar fields: ringdown and tails
After extensively testing the code, we now explore the dynamics of generic, massive scalar fields in the background of a Schwarzschild or Kerr BH with a/M = 0.99. Therefore we initiate the field as a Gaussian wave packet, according to Eqs. (7.16), typically with a width of w = 2M centered around r 0 = 12M . The progression of scalar fields with mass coupling M µ S is carried out by the time evolution system, Eqs. 7.9 and 7.10. We summarize the specific configuration of our simulations, such as the mass coupling M µ S , the modes contributing to the initial Gaussian and the grid setup in Table 7.4. Unless denoted otherwise, we set the excision radius to r exc = 1.0M . In order to read off the results we interpolate the scalar field onto a sphere with fixed radius r ex /M as seen by an observer at infinity and decompose it into modes using spherical harmonics Y lm (θ, φ), according to Eq. (7.23).
In Table 7.5 we summarize the results, including the quasinormal frequency ω and the fall-off of the tail Ψ ∼ t p estimated from the numerical evolution. We compute the real part of the QN frequency ω R by considering various intervals, each about 2 − 5 cycles long, and performing a regression on the thus obtained values. We obtain the imaginary part of the QN frequency ω I , which encodes information about the decay or growth rate τ = 1/ω I , by performing a linear fit of the form ln Ψ lm =A + ω I,lm t . (7.26) This dependence follows from Ψ lm ∼ exp(−ıωt). However, because in some cases the decay rate varies over time and because of resonant or beating effects that we observe in some (other) cases, which we will discuss in more detail in Sec. 7.3.5 below, a clean fit of that form is not always possible. In order to estimate the power-law fall-off of the late-time tails Ψ ∼ t p we employ the fit ln |Ψ| ∼p ln t . (7.27) We start by analysing massive scalar fields in a Schwarzschild BH background considering a number of different mass couplings varying from M µ S = 0.1, . . . , 1.0.
In Fig. 7.6 we present the logarithm of the l = m = 1 mode of the scalar field with mass coupling M µ S = 0.1. We find the QN frequency M ω 11 = 0.293 − ı0.036. In this case, we fit the tail to the oscillatory inverse power-law, Eq. (7.25), expected for small mass couplings at intermediate late times. Numerically we obtain Ψ 11 ∼ t −2.543 sin(0.1t), which agrees with the theoretically expected fall-off within less than 2% [324,325]. The corresponding fit is presented in Fig. 7.6 (red-dashed line) which describes the numerical data well for t ≥ 100 M .
Next, we consider a scalar field with mass coupling M µ S = 0.42. In order to perform a convergence analysis of its l = m = 0 mode, we have evolved the setup (run S m42 a0 in Table 7 We show the corresponding convergence plot in Fig. 7.8. In particular, we present the differences between the coarse -medium and medium -high resolution run. The latter has been amplified by the factor Q 2 = 1.28 demonstrating second order convergence. We estimate the numerical error to be ∆Ψ 00 /Ψ 00 ∼ 3.6% at t ∼ 1000 M which increases to ∆Ψ 00 /Ψ 00 ∼ 6.7% at t ∼ 1500 M . Additionally, we have simulated this field with a more generic wave packet for a longer time by employing higher resolution h = M/60 and a larger outer boundary located at r B = 1536 M . The following results refer to this run. In the left panel of Fig. 7.7 we present the logarithm of the l = m = 0 mode of the scalar field extracted at different radii r ex = 25, 35, 45, 55M . In the right panel of Fig. 7.7 we show the log-log representation of the l = m = 0 mode, extracted at the radii r ex = 25M , together with the fit to the late-time power-law tail, Eq. (7.27). Numerically, we find the exponent p = −0.880 which agrees within 5.6% with the theoretically expected late-time behaviour p = −5/6 [324][325][326].
The l = m = 1 mode of the scalar field with mass coupling M µ S = 0.42, depicted in Fig. 7.13 for extraction radii r ex = 22.5, 30, 40, 50M , exhibits a much richer structure. In particular, we observe an overall modulation of the signal whose specific features depend on the location of the measurement. This behaviour can be understood as a result of beating between different overtone modes whose amplitude depends on the specific location of the measurement. We will discuss this effect in more detail in Sec. 7.3.5 below. Because of this beating effect we did not determine the decay rate or late-time tail of this mode.
These features of the beating effect appear also in the l = m = 1 component of a scalar field with mass coupling M µ S = 0.55, albeit not as pronounced. In Fig. 7.14 we depict the logarithm of the l = m = 1 mode, extracted at r ex = 13, 20, 30, 40M . The first extraction radius corresponds to the node of the first overtone of the bound state of a scalar field with M µ S = 0.55 in Schwarzschild background. The behaviour of the l = m = 0 mode is similar to the case with M µ S = 0.42 and we summarize the (numerically computed) QN frequency and late-time tail in Table 7.5.
Finally, we have studied a scalar field with larger mass parameter M µ S = 1.0 evolving in the background of a Schwarzschild BH. In the left panel of Fig. 7.9 we show the logarithm of the l = m = 1 mode, extracted at radii r ex = 10, 20, 30, 40M . We numerically find the QN frequency M ω 11 = 0.965 − ı0.0046. We show a log-log representation of the l = m = 1 mode, measured at r ex = 20M , together with a fit to the late-time tail in the right panel of Fig. 7.9. We estimate the fall-off exponent to be p = −0.873 which is in good agreement, within less than 4.9%, with the previous results presented in [326].
The numerical results of massive scalar fields in a Schwarzschild BH background, specifically the QNM frequencies and power-law tails, are in good agreement with theoretical predictions and previous numerical results (see e.g. [211,[324][325][326] and references therein). Thus, we have verified that our implementation works well and delivers reliable results.
In the second part of this section we focus on the time evolution of a massive scalar field with mass couplings M µ S = 0.42 or M µ S = 0.55 in the background of a Kerr BH with spin parameter a/M = 0.99. Initially, we choose a Gaussian profile of the scalar field located at r 0 = 12M with width w = 2M . The specific setup is summarized in Table 7.4. The excision radius is set to r exc = 1.0M in order to guarantee that it always lies within the BH event horizon. In order to read off the results we interpolate the scalar field onto a sphere at extraction radii r ex = 20, 26.5, 30, 35, 40, 45, 50, 60M as seen by an observer at infinity and decompose it into modes using spherical harmonics Y lm (θ, φ), according to Eq. (7.23).
First, we analyse the l = m = 0 multipoles of both scalar field configurations, which are depicted in Figs. 7.10 and 7.12. In particular, we present the logarithm of the waveforms measured at r ex = 20, 26.5, 35, 45M in the respective left panels of the figures. We estimate the real part of the frequencies to be M ω R,00 = 0.419 ± 0.01 and M ω R,00 = 0.549 ± 0.01 for the scalar field with mass coupling M µ S = 0.42 and M µ S = 0.55, respectively. A clean fit to the imaginary parts of the frequencies has not been procurable because of variations in time. In the right panels of Figs. 7.10 and 7.12 we show a log-log representation of the l = m = 0 mode, extracted at r ex = 20M , in order to visualize the late-time tail. Specifically, we present the numerical data together with a fit of the form ψ ∼ t p with the power-law fall-off p = −0.861 and p = −0.847 in the case M µ S = 0. 42 Table 7.4.: Initial setup for simulations of a massive scalar field with Gaussian initial data located at r 0 = 12M and with width w = 2M in a Schwarzschild or Kerr background with a/M = 0.99. In particular we denote the mass parameter M µ S , the modes of the initial pulse 0 Σ(θ, φ) and the specific grid setup, measured in units of the BH mass M , following the notation of Sec. II E in [161].
good agreement, within ≤ 3.5% and ≤ 2%, with the prediction p = −5/6 [324][325][326]. Additionally, we have performed a convergence analysis for the evolution of a massive scalar field with coupling M µ S = 0.42 in order to estimate the numerical accuracy of the waveforms. Therefore, we have simulated the setup with three different resolutions h c = M/60, h m = M/72 and h h = M/84 (runs S m42 a99 c , S m42 a99 m and S m42 a99 f in Table 7.4). In Fig. 7.11 we present the corresponding convergence plot for the l = m = 0 (left panel) and l = m = 1 (right panel) modes. The numerically found convergence factor is Q 2 = 1.66, thus demonstrating second order convergence. This procedure allows us to estimate the discretization error to be about ∆ψ/ψ ≤ 1% in both modes.
In summary, the l = m = 0 modes of the considered massive scalar fields with coupling M µ S = 0.42 and M µ S = 0.55 in Kerr background with a/M = 0.99 show the expected ringdown signal followed by a late-time power-law tail. The numerically obtained properties of the waveforms, such as (real) part of the ringdown frequency at late-time fall-off are in good agreement with theoretical predictions and previous numerical results [324][325][326]. However, the corresponding l = m = 1 modes, shown in Figs. 7.15 and 7.16, exhibit interesting features which can be explained by a beating phenomenon. We analyse this behaviour in more detail in the following section.

Mode excitation and beating
We now discuss an exciting effect which we observe in the dipole modes of the scalar fields with mass couplings M µ S = 0. 42      We present the differences between the coarse-medium and medium-high resolution runs, where the latter is amplified by Q 2 = 1.66 indicating second order convergence. that r ex < r node , r ex ∼ r node , r ex > r node and r ex ∼ r max , where r node and r max denote the location of the node and local maximum (r max > r node ) of the first overtone of the corresponding quasi-bound state. Although strictly speaking being a different state, the quasi-bound state serves as excellent guide for the numerical simulations of a scalar field with more generic initial profile. By studying Figs. 7.13, 7.14, 7.15 and 7.16, we observe two effects: (i) a modulation of the l = m = 1 waveform and (ii) a dependence of this modulation on the location of measurement. We further note, that these effects are more pronounced in the case of the scalar field with mass coupling M µ S = 0.42 than for M µ S = 0.55. Therefore, we will restrict ourselves to the analysis of the former case. The two phenomena can be understood by a beating effect between the fundamental and overtone mode and their space dependent excitation. In the following, we will investigate both effects in more detail. The observed modulation or envelope of the waveforms can be understood by considering the beating between the fundamental and n-th overtone mode with similar real parts of the frequencies, i.e., ω 0 = ω R,0 +ıω I,0 and ω n = ω R,n +ıω I,n with ω n,R = ω 0,R +δ and δ 1. The imaginary parts are related via ω I,n = ω I,0 /(n + 1). Then, if we restrict ourselves to only two competing modes, the scalar field can be expressed as where A 0 and A n are the respective amplitudes of the fundamental and overtone mode. In order to estimate the frequencies of the modes, their difference δ and amplitudes we fit the numerical data to the real part of the scalar field which is (approximately) given by (Ψ) ∼A 0 exp(ω I,0 t) + A 1 exp(ω I,0 t/2) cos(δt) . Here, we assume a beating between the fundamental and first overtone mode. The results of this fit for the scalar field with M µ S = 0.42 are summarized in Table 7.6. In some cases, the imaginary part of the frequency could not be fitted accurately because of numerical noise. These cases are denoted by n.a.. We find good agreement between the fitted and theoretically predicted values of the real parts of the frequencies ω R,0 and ω R,1 and their difference δ (see Table 7.1), within less than 7%. Additionally, we plot the fitted curves together with the numerical ones in Figs. 7.13 and 7.15 for specific values of the extraction radii. We have focused on the cases r ex < r node , r ex ∼ r node and r ex ∼ r max > r node . This is a reasonable choice because the generic scalar field eventually approaches the bound state, as becomes evident in the top right panel of Fig. 7.15. The plot corresponds to the l = m = 1 mode measured at r ex ∼ r node , where the amplitude of the (quasi-bound state) overtone mode vanishes. Whereas the waveform shows a small modulation at early times, the signal flattens at late times, indicating that the overtone amplitude decreases, thus indeed approaching the quasi-bound state.
Let's now focus on the scalar field with M µ S = 0.42 in Kerr background, Fig. 7.15. The plots nicely illustrate the interplay between the fundamental and overtone mode and the dependence of their amplitude on the location of measurement. In the case r ex < r node (top left panel of Fig. 7.15), the amplitudes of the fundamental and overtone mode are both present with A 0 > A 1 . The signal exhibits a clear, though not very strong modulation of the waveform. In the case r ex ∼ r node (top right panel of Fig. 7.15) we observe a slight modulation of the signal at early times which flattens at late times. This indicates that the amplitude of the overtone mode decreases and eventually only the fundamental mode will be present. In the case r ex ∼ r max we find that the amplitudes of the fundamental and overtone mode are similar, A 1 ≤ A 0 . This results in a strong modulation of the scalar field as can be seen in the bottom left panel of Fig. 7.15.
We have observed that different modes are excited to different amplitudes depending on initial data and the location where the mode is being measured. In a seminal work, Leaver has established some rigorous results in this regard, by deforming the contour integral in the complex-frequency plane; for further details we refer the reader to the original work [313] and to other comprehensive follow-ups [211,[327][328][329]. The upshot is that each quasinormal mode, which corresponds to a pole in the complex-frequency plane, is excited to a different degree depending on the initial data and on the mode in question. In the frequency-domain, it is useful to represent the scalar field by the Klein-Gordon equation in the form where V is a mass-dependent potential (see for instance [85] for details). The function I(r) is a generic source term, which describes for instance the initial data when Ψ(ω) is a Laplace transform of the time-domain wavefunction. The QNM contribution can be isolated from other features of the signal, such as the late-time tail, using the Green's function technique [313,327,328]. First, one defines a solution of the homogeneous equation having the correct behavior at the horizon (only in-going waves), where q = ± µ 2 S − ω 2 . The excitation amplitude of each mode is characterized by the quantity C n ∝B nΨr + (ω n , r) I(r)Ψ r + (ω n , r)dr . (7.32) Here, the functionΨ r + (ω n , r) ≡ Ψ r + (ω n , r)/A out is the normalized wavefunction, to be evaluated at the quasinormal frequency ω n . Thus, two quantities are crucial to determine the excitation of a initial data, and the excitation factor B n which is equal to Thus, for instance, the relative amplitude between different modes depends strongly on the point where this amplitude is evaluated: if it is close to a node, the mode in question will have a very small amplitude (see Eq. (7.32)): by definition a mode is not excited at its node. Likewise, Eq. (7.32) implies that localized initial data close to the node of the mode does not excite the mode in question, a well-known result for closed systems [328].
We have not attempted a complete description of mode excitation for this work, a preliminary analysis indicates that the excitation factors B n differ appreciable for different quasi-bound states overtones, contrary to the usual quasinormal modes [328]. For instance, for a = 0, the fundamental bound-state with M ω = 0.407524 − 0.001026i has B n ∼ 10 −6 while the first overtone M ω = 0.4146906 − 0.0004053i has B n ∼ 10 −9 . On the other hand N n /A 2 out is roughly of the same order of magnitude for both modes.

Conclusions
In this chapter we have studied in detail the evolution of massive scalar fields in the background of Schwarzschild or Kerr BHs in the time-domain. Therefore, we have presented a time evolution code Lin-Lean, which enables us to model the progression of massive (scalar) fields in the background of BHs. The motivation to revisit the evolution of massive scalar fields in the linear regime, which has been widely studied in the literature [84-88, 90, 97, 102, 324-326, 330], has been twofold: (i) extensively testing and verifying our code in preparation for future applications; (ii) to shed some light onto the puzzling results reported by Strafuss & Khanna [102]. The authors found an instability growth rate M ω I ∼ 2 · 10 −5 for a massive scalar field with M µ S = 0.25 in Kerr background with a/M = 0.9999 -two orders of magnitude larger than the results of frequency domain calculations for the maximum possible instability timescale [87,88].
As initial configurations we have considered quasi-bound states, which are localized in the vicinity of the BH and represent very clean states. The second class of initial data have been Gaussian wave packets, representing more general and therefore possibly more realistic setups.
The first part has been devoted to check our implementation. Therefore we have evolved a spherically symmetric scalar field with unphysical, spatial dependent mass coupling which yields an instability. The numerically computed growth rate agrees with the theoretical one within less than 1%. As a second benchmark test we have evolved massless scalar fields in Schwarzschild and highly spinning Kerr background with a/M = 0.99. The thus estimated QN frequencies and, if applicable, late time tails have been in excellent agreement with theoretical predictions [211]. Notice, that the more challenging case of highly rotating background has required very high resolution of the numerical grid close to the BH as well as in the wave extraction zone.
Our next step included the numerical construction of quasi-bound state initial data of a dipole scalar field with mass coupling M µ S = 0.42. As background we have choosen a Kerr BH with spin parameter a/M = 0.99 which is known to give rise to the strongest instability [87,88]. The thus computed fundamental and overtone modes of the quasibound state have been in excellent agreement with a comparative computation employing Leaver's method [87,88,313]. Additionally, in a time evolution we have verified that the amplitude of the field hardly varied with time. However, we have not been able to observe the growth of the instability because its expected value is tiny, M ω I ∼ 1.5 · 10 −7 , and therefore completely buried in numerical noise.
Finally, we have evolved scalar fields with mass coupling M µ S = 0.1, 0.42, 0.55, 1.0 in Schwarzschild and M µ S = 0.42, 0.55 in highly rotating Kerr background with a/M = 0.99. As initial configuration we have chosen generic Gaussian wave packets. We have first focused on the dipole modes of the field with small and large mass couplings M µ S = 0.1 and M µ S = 1.0 as well as the monopole modes for M µ S = 0.42 and M µ S = 0.55. Their numerically computed QN frequencies and late-time tails have been in good agreement, overall within less than ∼ 5%, with previously reported results [87,88,[324][325][326]. The time evolution of the dipole modes of the scalar field with mass couplings M µ S = 0.42 and M µ S = 0.55 have revealed an interesting phenomena: a modulation of the signal over time, where the amplitude of the envelope strongly depends on the location of the measurement. This behaviour can be understood by the beating effect between fundamental and overtone modes present due to the generic initial data. Indeed, the real parts of the fundamental and overtone frequncies as well as the beating frequncy δ estimated from the numerical data are in remarkable agreement with theoretical predictions summarized in Table 7.1.
A further analysis and comparison to the corresponding quasi-bound states revealed that the excitation of the respective modes is space dependent. For example, the amplitude of the overtone mode at its node vanishes by definition, thus suppressing any modulation. Indeed, our numerical data match these expectations (as can be seen in the top right panel of Fig. 7.15): The dipole mode of the generic scalar field extracted at the overtone's node of the corresponding quasi-bound state shows little variation at early times and relaxes to an almost constant amplitude at late times. Further evidence is given by the signal extracted at the local maximum of the overtone. At this position the amplitudes of fundamental and overtone mode are similar, resulting in a strong modulation of the waveform. Indeed, we observe the expected behaviour in the bottom left panel of Fig. 7.15.
The observed beating effect, supported by the findings in Ref. [101], has the potential to explain the puzzling result reported by Strafuss & Khanna [102]. As in our case, they have evolved a Gaussian wave packet which might excite not only the fundamental but also overtone modes. A rough analysis employing Eq. (7.28) with their setup reveals that beating should be present. This results in a local growth rate of M ω I ∼ 10 −5 as has been observed in [102] and misinterpreted as overall instability growth of the massive scalar field.
By developing and testing our 3 + 1 time evolution code that simulates massive scalar fields in the background of BH spacetimes we have performed first important step towards the exploration of more challenging scenarios: (1) Evolution of Proca fields in BH backgrounds: A further exciting line of research is the investigation of massive vector fields in BH backgrounds. They are expected to give rise to a much stronger superradiant instability than the scalar case. Previous investigations have focused on massive vector fields in Schwarzschild [331][332][333] and Tangherlini [334] backgrounds. Especially, Rosa & Dolan [333] have recently reported on the complete QNM and bound state spectrum of Proca fields in Schwarzschild backgrounds. However, the Kerr case proves the be much more challenging because the equations of motion appear to be non-separable for highly spinning BHs (but see [308,309] for a computation in the slow rotation approximation). Therefore, our code provides the ideal setup, with minor modifications, to address this challenging task. First results will be reported elsewhere [9].
(2) Self-interacting scalar fields: Another class of interesting problems is given by massive scalar fields whose dynamics are described by additional non-linear terms, modelling their self-interaction. This open issue has first been addressed by Yoshino & Kodama [97] who modelled the collapse of a so-called bosenova.
(3) Backreaction effects: As far as we are aware, all studies involving massive scalar fields have been performed in the linear regime, i.e., neglecting backreaction effects of the BH spacetime. Therefore it is of utmost interest to explore the fully non-linear regime, which allows for the investigation of backreaction of the spacetime, such as the spin-down of the BH due to (subsequent) superradiant scattering. This type of studies would enable us to throw a glance at the end-state of the superradiant instability or, possibly, equilibrium configurations.
There is an entire playground of exciting future applications of massive fields in BH spacetimes. Our code is perfectly suited to address these challenging problems, which are subject to ongoing work. γ µ ν denotes the projection operator onto the hypersurface and * denotes the Hodge dual. By using the Gauss-Codazzi equations we express the electromagnetic components in terms of the "3+1" variables [161] E ij =R ij − γ kl (K ij K kl − K ik K jl ) , (A.6a) If we insert Eq. (A.4), the definition of the null-tetrad (A.2) and the expressions (A.6a), (A.6b) into the definition of the Newman-Penrose scalars (2.73a), (2.73b) we obtain In the numerical code we use these relations in order to calculate Ψ 0 and Ψ 4 on the entire Cartesian grid. Then, they are interpolated onto coordinate spheres of various extraction radii r ex . The Newman-Penrose scalars Ψ 0 and Ψ 4 are decomposed into spin-weighted spherical harmonics s Y lm according to where cos(2φ) = cos 2 φ − sin 2 φ, sin(2φ) = 2 cos φ sin φ, cos(3φ) = 4 cos 3 φ − 3 cos φ sin(3φ) = 4 cos 2 φ sin φ − sin φ C. Appendix -Black holes in higher dimensional spacetimes C.1. Analysis of troublesome terms at y = 0 The right-hand sides of Eqs. (4.27) and (4.28) contain various terms which cannot be evaluated directly at y = 0 because they involve explicit division by y. Although these terms are regular by virtue of a corresponding behaviour of the numerators, they need to be explicitly evaluated in the numerical implementation. In this Appendix we outline how the regularity of these terms can be implemented in a simple and efficient manner. For convenience we use a special notation: late latin indices i, j, . . . run from 1 to 3, covering x, y and z, but early latin indices a, b, . . . take values 1 and 3 but not 2, i.e., they cover x and z but not y. We begin this discussion by describing a simple manipulation which underlies most of our regularisation procedure. Consider for this purpose a function h which is linear in y near y = 0, i.e., its Taylor expansion is given by h(y) = h 1 y + O(y 2 ). From this relation we directly obtain This trading of divisions by y for partial derivatives extends to higher orders in a straightforward manner and will be used throughout the following discussion. Next, we consider the right-hand sides of Eqs. (4.27) and (4.28) and summarize the potentially troublesome terms as follows Here f stands for either of the scalars or densities ζ, χ and α.

(C.22)
The benefit in using these expressions is purely numerical: we do not need to store the inverse metric in grid functions which reduces the memory requirements of the simulations.

(C.23)
For a radially in-falling massive particle, starting from rest at r = r 0 , the energy per unit mass is f (r 0 ). The geodesic equation may then be written as This coordinate system encodes a spacetime slicing with zero shift and constant (unit) lapse, i.e., of type (4.76) with α = 1, β a = 0 for all times. To compare it with a numerical evolution we must have the initial data for the spatial metric written in a conformally flat form. Taking the initial hypersurface to be τ = 0 we see that this is achieved by a coordinate transformation r 0 → R with This actually coincides with the standard coordinate transformation from Schwarzschild to isotropic coordinates in five dimensions. The line element finally reads (5.2). At the initial hypersurface τ = 0, ρ 2 + z 2 2 dz 2 + dρ 2 + ρ 2 dθ 2 + ρ 2 sin 2 θdΩ 1 , (C. 28) where we have used the metric on the 3-sphere in the form dΩ 3 = dθ + sin 2θ (dθ 2 + sin 2 θdΩ 1 ) , (C. 29) and performed the coordinate transformation (R,θ) → (ρ, z) defined as ρ = R sinθ , z = R cosθ . (C.30) Using (4.82) we get ds 2 τ =0 = r 0 ( x 2 + y 2 + z 2 ) x 2 + y 2 + z 2 2 dx 2 + dy 2 + dz 2 + y 2 dΩ 1 . (C.31) Thus the coordinate transformation from the spherical coordinates (R,θ, θ) used in Eq. (5.2) to the "incomplete" Cartesian coordinates (x, y, z) used in the numerical evolution is x = R sinθ cos θ , y = R sinθ sin θ , z = R cosθ , (C.32) which resembles the usual coordinate transformation from spherical polar coordinates to Cartesian coordinates in R 3 ; but note thatθ and θ are both polar angles with range [0, π], which is the manifestation of the Cartesian coordinates "incompleteness".
The coordinate change (C.32) brings the five dimensional Tangherlini metric in geodesic slicing to a conformally flat form at τ = 0. This matches the initial data for the numerical evolution. One may ask, however, if the coordinate transformation evolves, in order to compare the analytic form with the numerical evolution. This cannot be the case, since the existence of τ -dependent terms in the coordinate transformation would imply a drift away from geodesic slicing. We are thus guaranteed that the coordinate transformation (C.32) is valid for all values of τ . Then, we can predict the value of the metric components that should be obtained from the numerical evolution; say γ xx should be, at time τ γ xx (τ, x, y, z) = x 2 g RR (τ, R) R 2 + x 2 z 2 gθθ(τ, R) R 4 (x 2 + y 2 ) + y 2 g θθ (τ, R) (x 2 + y 2 ) 2 , (C.33) where R 2 = x 2 + y 2 + z 2 and g RR (τ, R), gθθ(τ, R), g θθ (τ, R) are readily obtained from Eq. (5.2) with (C.29) and (C.32). The result forγ xx along the x-axis is plotted in Fig. 5.1 for various values of τ . For D ≥ 6 the situation is more involved because equations (C.24) can no longer be integrated straightforwardly, but require a numerical treatment. First one notices that the coordinate transformation (t, r) → (τ, r 0 ), with initial conditions t(τ = 0) = 0 and r(τ = 0) = r 0 , brings the D dimensional Tangherlini metric to the form
The areal radius r is related to R by a reparametrization R = R(r), given by Eq. (C.77), which depends on the components gθθ, g θθ only. As shown in Secs. 5.2.2 and 5.3.1, we find that this reparametrization is nearly constant throughout our numerical simulations. Therefore, the quantities g rr , g tr , g rθ , g rθ can be obtained from g RR , g tR , g Rθ , g Rθ by a simple rescaling: because dR dr 1 , (C. 40) we have g rr g RR , and similar relations hold for the other components.

C.4. Harmonic expansion of axisymmetric tensors in D dimensions
As discussed in Section 4.6.3, scalar spherical harmonics in D dimensions S l (θ, θ, φ 1 , . . . , φ D−4 ) are solutions of Eq. where h ab = h ab (t, r,θ), h aθ = h aθ (t, r,θ), hθθ = hθθ(t, r,θ), h θθ = h θθ (t, r,θ), C We also note that the background Tangherlini metric depends on the l = 0 harmonic only; the integral of its components over l ≥ 2 harmonics vanish. Therefore, if we decompose the spacetime metric (see Appendix C.3) as g µν = g r k ∂ r f (r) f t + r k ∂ t H T . (C.73) Conversely, since the perturbations do not depend on the l = 0 harmonic, the background metric g µν can be obtained as follows: Finally, to compute the areal radius r we note that gθθ = r 2 +hθθ and g θθ = r 2 sin 2θ +h θθ . Both the perturbations hθθ and h θθ contain harmonics of different type (S l , S l ,īj ); to extract the background we need the combination in Eq. (C.71f): We have also performed simulations using condition (6.1) on a cubic outer boundary. These simulations represent the inspiral of nonspinning BH binaries with initial separation d/M = 6.514 and initial linear momentum P y i /M = ±0.133. The grid setup for

Run
Grid Setup h f /M R B /M b.c. O1 (192,96,48,24,12,6)  these runs is listed in Table D.1 together with a reference model O1 which describes the inspiral of the same binary in an asymptotically flat spacetime using outgoing radiation boundary conditions. Gravitational waves have been extracted at r ex = 20 M in the form of the Newman-Penrose scalar Ψ 4 . In Fig. D.1, we compare the l = m = 2 mode of Ψ 4 obtained by the evolution of models C1.1 and O1. In Fig. D.1 we present the convergence plot of the l = m = 2 mode of Ψ 4 obtained from evolutions of models C1.1, C1.2 and C1.3. The difference between the medium and fine resolution result has been amplified by the factor Q = 1.58 corresponding to fourth-order convergence. While the overall convergence is about fourth-order as in the case of a spherical shell, the cubical outer boundary introduces a substantial amount of numerical noise which is also demonstrated in Fig. D.2 which displays the l = 2, m = 2 and the l = 4, m = 4 multipoles of Ψ 4 for models C1.3 and IN1. We believe that this is at least partly a consequence of mode mixing in the case of the cubical boundary shell which is not well-suited for the geometric shape of the gravitational wave pulse. For this reason, we have exclusively used a spherical shell in the main part of Chapter 6.

D.2. Snapshots of black holes in a box
In this section we present snapshots of the simulation that evolves BH binaries enclosed by a mirror-like box, presented in Chapter 6. In Fig. D.3 we illustrate the emission of the gravitational wave signal during the inspiral and merger and its evolution in the closed (confined) system containing a central, spinning BH. We display snapshots of the waveforms by superposing (the real part of) Ψ 0 and Ψ 4 as obtained for model VIS of Table 6.1. We show a slice of the orbital plane with x, y = −48M, ..., 48M during an interval t/M = 150, ..., 540. The difference in time between the individual pictures is ∆t/M = 10. The series of snapshots starts in the late inspiral phase shortly before the plunge and we see a strong gravitational wave signal that is emitted throughout the merger (first row and first two columns of the second row). This signal reaches the spherical boundary and is reflected back as can be seen in the final three panels of the second row and first two panels of the third row. Starting with the third snapshot in the third row we see a second pulse going outwards again after it has been scattered off the BH. This process is repeated several times and the series of snapshots ends with the fourth outgoing wave pulse. An animation constructed from the numerical data can be found in [335].