Axion Dark Matter around Black Holes

Track: Gravitation and Astroparticle Physics

Master Thesis

Axion Dark Matter around Black Holes

The Initial Value Problem of Superradiant Bound States in Numerical Relativity

by

Christoffel Doorman

August 2020

One of the most profound mysteries of our Universe is the nature of dark matter. A novel dark matter candidate which is predicted by many beyond the Standard Model theories is the axion. Large clouds of axions in the mass range [10−20, 10−10] eV can naturally arise around rapidly spinning astrophysical black holes through a classical wave-amplification process called superradiance. Recent studies based on perturba-tive methods have shown that the presence of a companion black hole can enrich the phenomenological features of these superradiant clouds. Energy level transitions can leave characteristic imprints in the gravitational wave signature of the binary that are in the detectable ranges of near-future detectors such as LISA. In this work we aim to go beyond the limitations of perturbative methods by mimicking the inspiral by a sequence of quasi-circular orbits and numerically solve the nonlinear Einstein field equations on constant-time slices. The main purpose of this thesis is to provide the theoretical framework and numerical setup to study a wide class of aspects that are inaccessible to perturbative methods. To this end, we derive the 3+1 (space+time) decompositions of the Einstein field and Klein-Gordon equations and perform a con-formal transformation of the spatial metric to obtain equations in a form well-suited for numerical implementation. The initial value problem is solved in the so-called XCTS approach with the open source spectral code Kadath. We performed a variety of tests of vacuum black holes and found robust results under changes in the boundary conditions and the numerical resolution. Our studies of the effect of the axion cloud revealed irregular solutions that we attribute to a degeneracy of the Klein-Gordon equation on the black hole horizon.

First and foremost, I would like to sincerely thank my supervisor Samaya Nissanke for making this exciting project possible and for all her guidance and support within and beyond this project. Second, I would like to acknowledge Christoph Weniger as the second reader of this thesis. In particular, I would like to express my deepest appreciation to Tanja Hinderer and Horng Sheng. The many extensive (Zoom) meet-ings and brainstorm sessions were ever constructive and disclosed numerous of the interesting aspects of this project. Tanja’s never ending knowledge and supportive attitude are admirable and have lifted many of the obstacles I thought were immov-able. Horng Sheng’s creative ideas and theoretical expertise have been most inspiring, and I appreciate all the time he has taken to patiently explain the various physical concepts and his valuable comments on this thesis. Furthermore, I am very grateful to Philippe Grandcl´ement, who has been so kind to host me for two weeks at the Laboratory for the Universe and Theory in Meudon to explain the many impressive functionalities of the Kadath library, and who has always come up with clever so-lutions to numerical problems we encountered throughout this project. I gratefully acknowledge the organisers and speakers of the Kavli-RISE Summer School at the University of Cambridge for giving me a broad introduction to numerical relativity. Finally, I wish to thank my fellow students, in particular Dylan, Kosio, Kostas, Sven and Lucas, for making the past two years such an enjoyable endeavour.

1 Introduction 5

2 Axion dark matter 9

2.1 Ultralight dark matter . . . 10

2.2 Superradiance . . . 12

2.3 Gravitational atom . . . 14

3 Gravity and spacetime geometry 18 3.1 General Relativity . . . 19 3.1.1 Spacetime metric . . . 20 3.1.2 Curvature of spacetime . . . 22 3.1.3 Black holes . . . 23 3.2 Geometry of hypersurfaces . . . 25 3.2.1 Hypersurface embedding . . . 25 3.2.2 Curvature of hypersurfaces . . . 28 3.3 Foliations of spacetime . . . 30

3.3.1 Lapse function and shift vector . . . 31

3.3.2 Evolution of spatial metric . . . 33

3.3.3 Projections of the Riemann tensor . . . 34

4 3+1 decomposition of equations of motion 39 4.1 Einstein field equations . . . 40

4.1.1 Constraint equations . . . 41

4.1.2 Evolution equations . . . 42

4.2 Klein-Gordon equation . . . 44

4.2.1 Stationary solutions . . . 44

5 The initial value problem 48

5.1 Conformal decomposition . . . 49

5.1.1 Spatial metric . . . 49

5.1.2 Extrinsic curvature . . . 52

5.1.3 Equations of motion . . . 54

5.2 Solving the initial data . . . 56

5.2.1 Extended conformal thin sandwich method . . . 58

5.2.2 Choosing the free initial data . . . 59

5.3 Gauge invariant quantities . . . 62

6 Numerical relativity results 64 6.1 The Kadath library . . . 65

6.2 Vacuum solutions . . . 68

6.2.1 Isolated black holes . . . 68

6.2.2 Binary black hole . . . 73

6.3 Scalar cloud configurations . . . 78

6.3.1 Isolated black holes . . . 78

6.3.2 Binary black hole . . . 82

6.3.3 Discussion . . . 83

6.4 Future directions . . . 85

7 Conclusion and outlook 87

A Mathematical preliminaries 89

B Equations 91

C Numerical code scripts 96

Introduction

The Standard Model of particle physics is an tremendously successful theory that describes the fundamental building blocks of our Universe. Its development in the 1970s was the culmination of decades of remarkable theoretical and experimental progress. Although the Standard Model provides an elegant description of all known particles and their weak, strong and electromagnetic interactions, we now know that the theory is incomplete. In particular, it describes only a small fraction of the total energy content of the Universe. Through accurate astronomical and cosmological observations, we know that the matter density of the Universe is dominated by a mysterious form of non-luminous matter which has an unresolved nature and does not fit in the present Standard Model [1–4]. This dark matter gives rise to one of the most challenging open problems of contemporary science, fostering numerous theories beyond the Standard Model. To date, however, none of these theories have been empirically proven to be correct and the search for dark matter and fundamental physics beyond the Standard Model is far from complete. At the same time, we have just entered a fascinating epoch in which extraordinary technological and experimental successes enable us to probe the only interaction that dark matter assuredly couples to: gravity. The rapidly evolving field of gravitational waves has the potential to expose areas of physics that have never empirically been explored [5]. In particular, the recent detections of gravitational waves emitted by inspiraling compact objects are – besides resources of valuable astrophysical information – promising probes to search for physics beyond the Standard Model. And while the first gravitational wave detection is still very recent, various new detectors and ambitious future space-based missions such as LISA and DECIGO are getting prepared to reveal the unknown nature of the bulk of our Universe [6, 7]. In addition to the impressive technological innovations of detectors, novel data analysis techniques including machine and deep learning are introduced to find the physics that is hidden in the ripples of spacetime [8].

To guide the gravitational wave experimentalists and engineers it is crucial that consistent theo-retical models are continuously being developed and improved. At present, many feasible models that try to solve the dark matter problem have been proposed, based on different hypotheses and spanning over 80 orders of magnitude in mass. A novel class of models that has recently been re-ceiving considerable attention is ultralight dark matter, consisting of hypothetical particles with masses that could be much smaller than the sum of the neutrino masses [9]. It has been shown that large condensates of these ultralight particles can recover the known large scale properties of dark matter [10, 11]. Due to their low mass, ultralight particles are very weakly coupled to ordinary matter and could have easily been missed by traditional experiments. However, these particles arise naturally in many theories beyond the Standard Model and are promising candidates to solve various open problems in contemporary physics [12].

A notable ultralight particle that currently gains increasing attention is the axion. Axions are light scalar bosonic particles that provide a novel solution to the strong CP problem in quantum chromo dynamics (QCD) and are promising dark matter candidates [13, 14]. While the mass of axions can span the entire sub-eV domain down to the Hubble scale of 10−33 eV, axions in the mass range [10−20, 10−10] eV are of particular interest as these bosons can naturally arise in the vicinity of rapidly rotating astrophysical black holes through a process called superradiance [15]. Superradiance is a classical wave amplification process in which ultralight bosons coherently extract energy and angular momentum from the black hole, forming a large boson ’cloud’ that extends to regions far beyond the event horizon.

An interesting characteristic of these superradiant clouds is that they are in many ways anal-ogous to the hydrogen atom in quantum mechanics [16, 17]. One such analogy is the discrete energy levels in the spectra of the cloud and the fine and hyperfine structure of the hydrogen atom. It has been shown that energy level transitions of these bosonic clouds can produce monochromatic gravitational waves that are in the detectable ranges of Advanced LIGO and LISA [18, 19]. Studying ultralight boson clouds around black holes that are bound in a binary system is therefore particularly interesting as the dynamics of the boson cloud can leave im-prints in the gravitational wave signature of the binary. A recent study based on perturbative methods has shown that the presence of a companion black hole enriches the phenomenologi-cal features of the cloud considerably [20, 21]. However, the analytic approximations used in the large binary separation regime break down when the binary separation decreases and the gravitational perturbation becomes strong and non-perturbative. Moreover, various nonlinear effects such as the backreaction and tidal deformations of the cloud have not yet been found analytically. Understanding these phenomena during late inspiral stage is crucial to probe new physics beyond the Standard Model with gravitational waves.

In this work, we aim to go beyond the limitations of previous studies by numerically solving the nonlinear configurations of a black hole binary where one or both of the black holes carries an ultralight boson cloud. The main novelty of this project is that we can study the effects of the presence of a boson cloud in a black hole binary system in a completely unexplored regime, beyond what is accessible to perturbative methods. However, instead of doing a full numerical relativity study, we approximate the dynamics of the inspiral as a sequence of exactly circu-lar orbits. Solving the Einstein field equations for an equilibrium configuration is substantially faster than for the generic case, and avoids many of the subtleties and difficulties of full nu-merical relativity simulations. Moreover, the quasi-equilibrium setting is more controlled and computationally efficient. Consequently, one can explore more of the parameter space, such as varying mass ratios and binary separations.

By considering quasi-equilibrium orbits, the problem essentially reduces to the initial value problem of numerical relativity [22]. In this formulation, one specifies a set of free data and then solves a subset of the Einstein field equations known as the constraint equations. These equations are time-independent and therefore yield a solution at a single time slice of spacetime. By connecting successive equilibrium configurations of different binary separations, one can mimic the inspiral by extracting information from approximations based on energy balance. From these approximations, one can then compute various physical observables of the system [23].

To solve the Einstein field equations numerically, we have to take several mathematical op-erations into account. First, we split the four dimensional equations into time and spatial components, known as the 3+1 (space + time) formalism of general relativity [24]. In this procedure, the ten independent components of the Einstein field equation are captured in four constraint equations and six dynamical equations. Second, we perform a conformal transforma-tion of the metric to find a set of elliptic second-order differential equatransforma-tions that is numerically more favourable. We then use the source code Kadath to solve our set of equation [25]. Kadath is a relatively new and modern library which uses spectral numerical methods and Newton-Raphson iterations to solve partial differential equations. In contrast to the few other publicly available initial data solvers, Kadath provides various functionalities that are particularly useful when solving the initial value problem of binary systems.

The work done for this thesis is part of an ongoing research project. The results given Chapter 6 represent work that is still in progress and can only be regarded as preliminary. As the numerical computation of initial data has proven to involve various subtleties, final results have not yet been obtained. Possible solutions to current obstacles are discussed and future directions will we presented in the outlook at the end of this thesis.

Outline This thesis is divided into three parts: the first part (Ch. 2-3) provides a broad motivation for this study and reviews the most relevant background; the second part (Ch. 4-5) treats the derivation of the equations that will be solved numerically; the third part (Ch. 6-7) presents some preliminary results and the outlook of the project. More explicitly, the chapters in this thesis are organized as follows. Chapter 2 reviews the dark matter problem and provides the motivation and context of the project. In Chapter 3, we review some of the key concepts of general relativity that are relevant for this work and discuss the 3+1 splitting of spacetime. In Chapter 4, we use these concepts to write the Einstein field equations and Klein-Gordon equations in 3+1 dimensions. In Chapter 5, we perform a conformal transformation of the metric to write the equations of motion in a mathematically more favourable form. In addition, we discuss how we choose the free data that are solved by the constraint equations. The numerical method we use and preliminary results of different black hole systems will be presented in Chapter 6. Furthermore, we discuss here the current obstacles and pivots that we encountered and provide possible next steps. Finally, we review our main conclusions and present the outlook of the project in Chapter 7. This thesis also includes three appendices, which contain additional reference material. In Appendix A, we briefly review Lie derivatives and the conformal Killing operator. In Appendix B, we summarize all the equations that are used in the numerical computations and presents the implementation of complex fields. Finally, we provide some of the numerical codes that are used to compute the results in Appendix C.

Notation and conventions Throughout this thesis we use a Lorentzian spacetime metric with signature (–,+,+,+). Unless otherwise stated, we work in natural units (G = ~ = c = 1). Furthermore, we set the cosmological constant Λ to zero. When it can not be made clear from the context, we denote the dimension of tensors with a prescript, e.g. 3R versus4R. We will use φ to represent the bosonic scalar field, while ϕ is reserved for the azimuthal angle in spherical coordinates (r, θ, ϕ). Moreover, we use Einstein convention to write summations over indexed terms. Spacetime indices will be denoted by Greek letters (µ, ν, ...), while Latin letters will be used to either represent spatial components (i, j, ...), vielbein indices (a, b, ...) or label indices (k, l, ...). Scalar field eigenstates are represented by |n`mi, where the integers {n, `, m} are the principal, orbital and azimuthal angular momentum quantum numbers respectively. We follow the atomic physics conventions where n ≥ ` + 1 and |m| ≤ `.

Axion dark matter

Axions are hypothetical low-mass elementary fields that interact very weakly with the Standard Model particles [26]. They are predicted by many beyond Standard Model (BSM) theories and provide promising solutions to various unsolved problems in contemporary physics [12, 16]. In particular, axions are found to be promising dark matter candidates [27]. In this work we consider axion-like ultralight scalar bosons which have masses that range from 10−20to 10−10eV, far below the sum of the neutrino masses, approximately at 1m eV, which are the lightest particles of the Standard Model [18, 19]. Previous work has found that these ultralight bosons can be coherently amplified around rapidly spinning black holes forming boson condensates through a classical wave-amplification process called superradiance [15, 28]. Such clouds have many interesting features and a structure that is in many ways analogous to the hydrogen atom in quantum mechanics [18]. If axions exist, these superradiant phenomena in the vicinity of black holes are particularly interesting because they can leave imprints in gravitational wave signatures that can be measured with existing or near-future detectors such as LIGO and LISA [19]. We first give a brief introduction to dark matter in general and ultralight dark matter in partic-ular in Section 2.1. Then, we discuss how ultralight bosons can form clouds around black holes and present the basic properties of superradiance in Section 2.2. Finally, we draw the analogy between the energy eigenstates of superradiant scalar clouds and eigenstates of the hydrogen atom in Section 2.3. Although the main purpose of this chapter is to provide the reader with some context and motivation for this project, some of the equations appearing in Section 2.3 will be used in later chapters.

2.1

Ultralight dark matter

One of the most profound mysteries of our Universe is the nature of dark matter. Since the 1930s, when Zwicky observed inexplicable large velocities of galaxies in the Coma cluster [29], scientists have began to question the consensus then that the bulk of the Universe consists of luminous matter.1 Over the last century, a number of evidences has been found that substantiate the large abundance of non-luminous matter in our Universe at different scales. For instance, flat rotation curves at galactic scales imply a uniformly distributed dark matter halo that extends to far regions beyond the galaxy [31]. At galaxy cluster level, discrepancies between the observed virial mass of luminous matter and the virial mass observed by gravitational lensing shows that the majority of matter in galaxy clusters is non-luminous [32]. In addition, evidence of dark matter appears at cosmological scale from the CMB angular power spectrum, in which anisotropies of the CMB are very well explained by the cosmological cold dark matter (ΛCDM, see next paragraph) model [33]. Over the past two decades, various observational studies have achieved remarkable progress in establishing the components of our Universe to high precision [1– 4]. It has been shown that only 4.9% of the energy content of our Universe consists of ordinary matter from the Standard Model [33]. The bulk of the energy content consists of dark energy (68.6%) and dark matter (26.5%). Whereas most physicists agree unanimously on the existence of dark matter, there is no consensus of the component or components dark matter consists of. Before discussing possible dark matter candidates, we first briefly state what we do know about the dark matter component of our Universe.

At present, the ΛCDM model is the most accurate cosmological model that describes our Uni-verse and is based on numerous precise astrophysical and cosmological observations. In the ΛCDM model, the Universe is assumed to be spatially flat and dominated by dark energy and cold dark matter. Dark energy, denoted with the cosmological constant Λ, drives the observed accelerated Hubble expansion of the Universe [34]. On the other hand, the cold dark matter (CDM) component explains to high accuracy the observed structures and dynamics of large-scale formations. The CDM component of the ΛCDM model has several properties. First, one can deduce from observations including the distribution and sizes of galaxies that CDM has to be non-relativistic (hence the name “cold” dark matter). Second, CDM is assumed to be non-baryonic as the predicted baryon density of the ΛCDM model agrees well with Big Bang nucleosynthesis [35]. Third, CDM does not couple to the electromagnetic or strong force, po-tentially couples to the weak force, but does interact with gravity. Last, for CDM to exist it is assumed to be stable in the sense that it has a lifetime comparable to that of the Universe.

Although the properties of CDM are well known on large scales, we do not have a clear under-standing or strong observational evidence of the components CDM consist of on smaller scales. Physicists have proposed numerous feasible models that explain dark matter dynamics on small scales while recovering the large scale properties of ΛCDM. These models are based on various different hypothesis and span over 80 orders of magnitude in mass scales, ranging from new elementary particles to astrophysical scale primordial black holes. A novel class of models in the low-mass regime considers dark matter as a cold boson condensate on galactic scales that resembles the properties of CDM on large scales [10, 11]. These boson condensates consist of non-thermal sub-eV bosonic fields and show impressive agreement with dark matter properties from both a theoretical and phenomenological perspective [9].

A considerably popular bosonic field that is proposed as a dark matter candidate is the axion. Axions are light scalar fields that have been introduced by Peccei in 1977 to solve the strong CP problem [13]. While astrophysics and cosmology provide an upper mass limit of approximately 10−3 eV, these QCD axions can have masses far down the sub-eV domain [14]. However, as the QCD axion mass is directly proportional to its φF ˜F coupling to photons, axions with masses below 10−5 eV are notoriously hard to detect directly using table-top experiments because of their weak interactions [36]. It has been shown that if QCD axions exist and solve the strong CP problem, a plenitude of ultralight axions can arise in the compactification of extra spatial dimensions in string theory, spanning the entire domain of sub-eV down to the Hubble scale 10−33 eV [16]. A subset of these axions in the mass range [10−22, 1] eV are found to be good dark matter candidates that could have easily been missed in experiments due to their low mass and weak interaction with the Standard Model [9].

Since the first direct detection of gravitational waves in 2015 [37], the field of gravitational wave research has experimentally and theoretically matured significantly and has already achieved remarkable successes in astrophysics. Considering the current pace of advancements in the field, the prospects of probing cosmology and fundamental physics with gravitational waves are bright. In particular, as dark matter interacts only through gravity, gravitational waves are the most promising observable to detect dark matter and therefore discover new fundamental physics beyond the Standard Model [5]. Gravitational waves are difficult to detect due to their weak nature and present gravitational wave detection is limited to merging compact objects, which produce relatively strong signals. For this reason, a the large number of studies is consider-ing dark matter models in presence of compact objects. Correspondconsider-ingly, studies on axions or axion-like particles in black hole systems has recently increased significantly. Supposing that axions exist, they can form boson condensates around black holes. In particular, large clouds of ultralight bosons can arise around rapidly spinning black holes by a process called

superradi-ance [15], which will be discuss in more detail in Section 2.2. Oscillations of these bosonic clouds can emit gravitational radiation themselves, possibly altering the gravitational wave signature of the hosting black hole [28]. In optimistic scenarios these gravitational wave imprints might be observable by detectors such as LIGO and LISA [19].

2.2

Superradiance

Superradiance is a wave amplification process that appears in various fields such as quantum mechanics, optics, astrophysics and relativity [15, 38, 39]. In black hole mechanics, superradiant instabilities of bosonic fields can occur when a scalar field is scattered onto a rotating black hole.2 However, these superradiant instabilities can also spontaneously arise from their vacuum quantum fluctuations. If the spin of the black hole is greater than the angular phase velocity of the incoming field, and if the angular velocities are aligned, the black hole can enhance the energy and angular momentum of the field [42]. Explicitly, considering a Kerr black hole of mass M and angular velocity ΩH, and a bosonic field of mass µ, angular frequancy ω and azimuthal

angular momentum m, superradiant instabilities can only grow when ω

m < ΩH = a 2M r+

, (2.1)

where a is the spin parameter of the black hole of spin J ≡ aM and r+ ≡ M +

√

M2_{− a}2 _is

the radius of the black hole’s event horizon [20]. Note that for a co-rotating black hole and field, which is required for superradiance to occur, the azimuthal angular momentum is positive (m > 0). Moreover, ultralight bosonic fields of mass µ that grow superradiant instabilities require Compton wavelengths λC ≡ ~/(µc) similar to or larger than the black hole gravitational

radius rg ≡ GM/c2, as for smaller λC the field would not be scattered but accreted. We now

define the gravitational fine-structure constant or gravitational coupling α as the ratio between the gravitational radius of the black hole and the Compton wavelength of the field

α = rg λC

= M µ , (2.2)

in natural units [20]. The fine-structure constant is the crucial parameter that controls the interaction between the geometry and the ultralight boson cloud. When the requirement α . 1 and (2.1) are met, large bosonic condensates can arise in the vicinity of the black hole, often

2_{Although superradiance instabilities of massive bosons are not restricted to scalar fields and can occur for}

e.g. vector and tensor fields [40, 41], we limit our discussion to superradiant instabilities of massive scalar fields because these are technically simpler to solve (see Section 2.3) and best motivated by most BSM-theories [20].

referenced to as a “cloud”. Although the most interesting regime of α occurs for λC ∼ rg, when

superradiant growth is most efficient, analytical results have only been found for (α  1) and (α  1) [15]. Whereas superradiant instabilities would not grow for α & 1, scalar fields where (α  1) do have the potential to grow clouds and admit a similar eigenvalue problem as the hydrogen-like solutions [18, 43] (see Section 2.3 for more details).

The requirement λC & rg implies bounds on the field’s mass µ induced by the hosting black

hole. It can be shown that the mass of the field ranges from 10−20 eV for supermassive black holes (≤ 108M) to 10−10 eV for stellar mass black holes (≥ 2M).3 Although massless fields

can experience superradiance as well, they cannot grow superradiant instabilities as the field has to have mass to induce a reflecting barrier. This “gravitational mirror” reflects the field back to the black hole, allowing the field to repeatedly extract angular momentum from the black hole [45]. This continuous extraction of angular momentum from the black hole results naturally in a spin-down. The theoretical upper limit of the amount of energy that can be extracted from spinning black holes is found to be 0.29M for extreme Kerr black holes [46]. Therefore, the scalar field can only extract a finite amount of mass from the black hole and grow to some limited extend, although the total amount of energy extracted from the black hole can be significant. Two different types of gravitational radiation can be emitted from three different phenomena in superradiant bosonic clouds. Energy level transitions of bosonic clouds [17] and axion an-nihilation into gravitons [28] can produce monochromatic gravitational waves that are possibly detectable by Advanced LIGO in case of stellar mass black holes and by LISA in case of su-permassive black holes. In addition, when the self-interaction of a superradiant cloud becomes stronger than the gravitational coupling to the black hole, the axion cloud collapses and gener-ates a gravitational wave burst that might be detectable in case of a supermassive black hole host [28]. Moreover, superradiant axion clouds arising around black holes that are bound in binary systems are promising scenarios as they can induce rich phenomenological implications during inspiral [20, 21]. Although gravitational wave emission from cloud oscillations in stellar mass black holes is hard to detect with current ground-based detectors, superradiant clouds of ultralight bosons in the mass range [10−16.5, 10−14] eV have found to be detectable by LISA in case of extreme mass-ratio mergers [47]. To conclude, studying ultralight boson clouds around black holes is notably interesting, especially considering the rapidly evolving field of gravitational waves research and near-future planned detectors and missions.

3

Even axion masses of 10−8eV can form superradiant clouds around primordial black holes of mass 0.01M[44].

2.3

Gravitational atom

The description analogous to the hydrogen atom has given valuable insights into the properties and dynamics of the bound states of superradiant clouds. To illustrate this, we consider a massive complex4 _{scalar field Φ describing the ultralight boson cloud. In a general spacetime,}

Φ satisfies the Klein-Gordon equation

( − µ2)φ = 0 , (2.3)

where  is the d’Alembertian operator of the spacetime metric and µ is the mass of the scalar field in natural units. In the non-relativistic limit and neglecting gravitational wave emission, a superradiant massive scalar cloud surrounding a black hole satisfies the Schr¨odinger equation with a Coulomb-like potential that is regulated by the single parameter α defined in (2.2) [15]. This can be shown by considering the ansatz

φ(t, r) = e

−iµt

√

2µψ(t, r) , (2.4)

where the complex scalar field ψ(t, r) varies on a timescale that is longer than µ−1 [20]. By substituting this ansatz into the Klein-Gordon equation (2.3), one can show that up to leading order of r−1 or α, we find the Schr¨odinger equation with a Coulomb-like potential

i∂ ∂tψ(t, r) =  − 1 2µ∇ 2₋α r  ψ(t, r) . (2.5)

Consequently, the eigenstates of the scalar cloud are hydrogen-like and the stationary eigenstates can be written as the eigenfunctions of the hydrogen-atom [43, 48]. We can therefore characterize the eigenstates of the scalar cloud with the principal, orbital and azimuthal angular momentum quantum numbers {n, `, m} respectively, which satisfy ` ≤ n − 1 and |m| ≤ `. For α  1, it can be shown that the eigenfrequencies of the field are [49]

ωn`m ' µ  1 − α 2 2n2 − α4 8n4 + (2` − 3n + 1)α4 n4_{(` + 1/2)</sub> + 2˜amα5 n3<sub>`(` + 1/2)(` + 1)}  , (2.6) 4

Although this analogy is not restricted to complex fields but applies to real fields as well, we restrict our analysis to the former for practical reasons. Whereas the energy momentum tensor of a complex field is stationary and axisymmetric, that of the real scalar field is dependent and non-axissymmetric. However, these time-dependence and non-axissymmetry are α-suppressed. To find the real solution of the complex scalar field, we just take the real part of the complex solution.

Figure 2.1: A schematic representation of the possible low-energy eigenstates (n ≤ 3) of the scalar field, analogous to that of the hydrogen atom in quantum mechanics. The solid blue lines with labeled eigenstates |n`mi represent growing modes, while the dashed red lines represent decaying modes. Reprinted from [17].

where ˜a = a/M = J/M2. However, in the small α limit, for simplicity we only consider the Bohr term ωn`m ' µ(1 − α2/(2n2)) to be contributing. For higher orders of α, one can find

the complete spectrum of the scalar field [17]. However, from (2.1) we know that not all energy eigenstates of the cloud are growing modes as m > 0 is a necessary requirement for superradiant instabilities to occur. A schematic representation of this spectrum for low-energy states (n ≤ 3) is shown in Figure 2.1.

A major advantage of considering a scalar field instead of e.g. a vector field is that the Klein-Gordon equation of a massive scalar field (2.3) is separable in Boyer-Lindquist coordinates which makes the problem technically simpler to solve [17]. To this end, we write for the scalar field the ansatz [50, 51]

φn`m= e−iωt+imϕRn`(r)P`m(cos θ) , φ∗n`m = eiωt−imϕRn`(r)P`m(cos θ) , (2.7)

where P`m are the associated Legendre polynomials, and Rn` are the Bohr radial functions

Rn`(r) = Ae−µαr/n  2µαr n ` L2`+1<sub>n−`−1</sub> 2µαr n  , (2.8)

Figure 2.2: An illustration of a decaying 1s mode (left) and a growing 2p mode (right) of the gravitational atom. The red arrow indicates the direction of the black hole spin. Figure adapted from [17].

reads [48] A = s 2µα n 3(n − ` − 1)! 2n(n + `)! . (2.9)

From the exponential decay of the radial function, we can estimate the radius rcof the cloud in

terms of the typical Bohr radius rc,n' n2(µα)−1 [17].

In our analysis we will consider two different eigenstates: the 1s state |n`mi = |100i and the 2p state |n`mi = |211i, shown in Figure 2.2. For the 1s state, we find that L2`+1<sub>n−`−1</sub>= P`m = 1 such

that at t = 0, the cloud has a radial profile only that decays exponentially:

φ_(1s)= e−iωtR_(1s)(r) = 2(µα)3/2e−(iωt+µαr). (2.10) Notice that because the total angular momentum number |m| = 0, the 1s profile is actually a decaying mode and does not grow a cloud around a black hole (see Figure 2.1). Nevertheless, the 1s state can be used as an important example to test the robustness of our numeric formulation in subsequent chapters. Because the cloud in this energy state has a radial profile only, cf. (2.10), it is in many ways simpler than higher energy modes as the inner boundary value of the field is constant over the black hole horizon.

For the 2p state, the generalized Laguerre polynomial is constant like in the 1s state. The associated Legrendre polynomial however is not: P_(2p)(cos θ) = −(1 − cos2θ)1/2. Therefore we find φ(2p)= − √ 6 12(µα) 5/2_rei(ϕ−ωt)−µαr/2_{sin θ . (2.11)}

In contrast to the 1s state, the 2p state has |m| 6= 0, which means it has a growing mode (see Figure 2.1). Furthermore, considering the radial profile we can infer that for regions near the black hole horizon the field grows linear with r, while at larger distances the field decays exponentially.

Gravity and spacetime geometry

The aim of this work is to study solutions of ultralight scalar fields around black holes in Ein-stein’s theory of general relativity [52]. This theory provides an elegant mathematical framework to describe the dynamics of and interaction between matter and spacetime. The physics of gen-eral relativity is governed by the famous Einstein field equations, which are a set of coupled and nonlinear second-order partial differential equations. This makes the Einstein field equations notoriously hard to solve and analytical solutions have only been found for highly idealized and linearized systems with high degrees of symmetry [53]. Moreover, exact solutions of the Einstein field equations of binary systems have not been found and can only be approximated during early-stage inspiral (e.g. large separation) by post-Newtonian [54, 55] or general relativistic perturbation theories when the mass ratio is small [56]. Previous work based on perturbative methods found intriguing phenomena of ultralight scalar fields in black hole binaries [20, 21]. However, the validity of the approximation brakes down in the most interesting regime when the binary approaches the merger phase and the system becomes highly nonlinear. To go be-yond these limitations and explore the properties of the ultralight scalar cloud during late-stage inspiral, one has to solve the full nonperturbative Einstein field and Klein-Gordon equations. The full Einstein field equations in the highly nonlinear regimes of a binary can be solved us-ing numerical relativity methods [57]. By direct numerical integration of the full non-linear Einstein field equations one can model a black hole binary without imposing symmetries and simplifications. To this end, one has to reformulate the Einstein field equations by splitting the four-dimensional spacetime into a three-dimensional spatial and a one-dimensional timelike part. With this procedure, called the 3+1 formalism, the Einstein field equations are decomposed into an equivalent system of elliptic and hyperbolic partial differential equations that are mathemat-ically more favourable and numermathemat-ically solvable. The splitting of the Einstein field equations

into 3+1 dimensions essentially reduces the system to a Cauchy problem with constraints [58]. In the 3+1 formalism four-dimensional tensors are decomposed orthogonally with respect to spatial slices of spacetime. These slices, called hypersurfaces, resemble the spatial dimensions of spacetime, on which the time-coordinate is constant throughout the hypersurface. A specific set of these hypersurfaces, called a foliation, resembles the dynamical four-dimensional spacetime. A superficial motivation of this specific splitting is that the tensor manipulations are performed on a three-dimensional spatial space where the scalar product is Riemannian. To project ten-sors either tangent or normal to these hypersurfaces, we require a comprehensive description of (i) a single hypersurfaces embedded in a manifold and (ii) a foliation of three-dimensional hypersurfaces spanning the four-dimensional spacetime.

The main purpose of this chapter is to provide the theoretical concepts and tools to decompose the Einstein field equations and Klein-Gordon equation into 3+1 dimensions. To this end, we first briefly recall some of the key concepts of general relativity in Section 3.1, which are needed to construct the 3+1 formalism. We then define a spatial hypersurface and discuss how hypersurfaces are embedded in spacetime manifolds in Section 3.2. Finally, we show how four-dimensional spacetime can be represented by a foliation of hypersurfaces in Section 3.3. In this chapter and the next, we will regularly map tensors between the spacetime manifold and hypersurfaces. In order to differentiate between tensors on the four-dimensional spacetime manifold and tensors on the three-dimensional hypersurfaces, we specify the dimension of the field as a prescript of the variable when this cannot be made clear from the context. Most of the time, however, the context will be sufficient to understand on which manifold a field is defined. Tensors such as the metric and connection are represented as different variables for different manifolds.

3.1

General Relativity

In this section we recall some of the key theoretical concepts and basic equations of general relativity that are relevant for this work. Note that this section should not be regarded as a complete summary of general relativity. For the reader that seeks a robust treatment of general relativity, we refer to the introductory textbooks by Caroll [59] or Schutz [60], and the advanced textbooks of Wald [58] or Misner, Thorne and Wheeler [61]. Most of the concepts presented in this chapter are taken from the former two, unless otherwise stated.

Arguably one of the most important equations in astrophysics and cosmology is the famous Einstein field equation. One way to derive it is by finding the equations of motion of the

Einstein-Hilbert action SEH = 1 16π Z d4x√−gR , (3.1) which is the standard gravitational action in general relativity. Here, g is the determinant of the spacetime metric g = det gµν and R is the Ricci scalar. One can show that by varying the

Einstein-Hilbert action with respect to the metric yields the Einstein field equations1 Gµν ≡ Rµν−

1

2Rgµν = 8πTµν, (3.2) where Gµν is called the Einstein tensor, Rµν is the Ricci tensor and Tµν the energy tensor. Note

that we have set the cosmological constant Λ to zero. The left-hand side of the Einstein field equations represents the curvature of spacetime which is determined by the metric gµν. The

right-hand side is related to any present matter. The Einstein field equations therefore illustrate how matter curves spacetime, and vice versa how curvature of spacetime dictates the dynamics of its matter content [62]. In particular, the Einstein field equations are essentially a set of differential equations of the spacetime metric and therefore represents how the dynamics of the metric reacts to the presence of energy-momentum.

3.1.1 Spacetime metric

To formulate Einstein’s theory of gravity we use tensor fields that are defined on manifolds. While tensors are objects that reflect physics, manifolds are the fundamental mathematical constructs that serve as the background spacetime on which the physics lives. In particular, the curvature of spacetime is represented by a four-dimensional differential and connected Lorentzian manifold, which is locally isometric to Minkowski spacetime.2 _{The non-degenerate metric that}

is defined on a Lorentzian spacetime manifold M is a symmetric rank-2 tensor that essentially captures all geometric and causal information of spacetime. Given a local coordinate system (xµ) = (t, xi) on M, the metric defines the line element

ds2 = gµνdxµdxν, (3.3)

where dxµ are the one-form gradients of xµ. Note that while the metric is dependent of the choice of coordinates, the line element ds2 is not. The line element represents the invariant

1

For the complete derivation of the Einstein field equations from the action principle we refer to e.g. Caroll [59].

2

distance between two separated points on the spacetime manifold and therefore has a meaningful function in computing lengths, surfaces or volumes. For instance, given some curve C on M that is parameterized by an affine parameter λ, we find for the length s on a finite interval of C(λ) s = Z ds = Z ds dλdλ = Z r gµν dxµ dλ dxν dλ dλ . (3.4)

In a similar fashion, one can use the metric to write surface or volume integrals of finite intervals. For instance, the volume element on a n-dimensional manifold readsp|g|dnx where the square root of the determinant of the metric makes the volume element invariant under transformations. Before we can describe the curvature of manifolds, we have to introduce the covariant derivative that “connects” nearby tangent spaces on a manifold. Covariant derivatives are in some sense the more general form of ordinary derivatives and are needed to describe local changes of tensors with respect to parallel transport. The covariant derivative of a tensor is a tensor itself with components ∇σTµ1...µmν1...νn = ∂σT µ1...µm ν1...νn+ m X i=1 Γµi σλT µ1..λ..µm ν1...νn− n X i=1 Γλ_σνiTµ1...µm ν1..λ..νn, (3.5)

where in the second and third term λ = µ1 and λ = ν1 for i = 1 respectively and λ = µm

and λ = νn for i = m and i = n respectively. Here, the symbols Γσµν represent the Christoffel

symbols that are defined as

Γσ_µν = 1 2g

σρ_(∂

µgνρ+ ∂νgρµ− ∂ρgµν) . (3.6)

Because the metric tensor is symmetric, one can easily notice that the Christoffel symbols are symmetric in its lower two indices. Furthermore, an important property of the invariance of the metric tensor gµν is that the covariant derivative of the metric is vanishing everywhere:

∇ρgµν = 0 . (3.7)

Connections that satisfy (3.7) are said to be metric compatible and such connections commute with the raising and lowering of indices. As a consequence, one can easily show that the covariant derivative of the inverse metric vanishes as well.

3.1.2 Curvature of spacetime

With this notion of the spacetime metric and associated connection, we can return to the left-hand side of the Einstein field equations (3.2). As discussed in the previous section, the curvature of spacetime is encapsulated in the metric tensor gµνthat is defined on a Lorentzian manifold M.

From this metric, one can eventually compute the Ricci tensor Rµν appearing in the Einstein

field equations. The Ricci tensor and Ricci scalar R are both contractions of the Riemann tensor, which is a tensor that measures the extent to which the spacetime metric is not locally isometric to Minkowski spacetime. One way to define the Riemann tensor is as the non-commutativity between covariant derivatives acting on a vector [58]

h ∇µ, ∇ν

i

vρ= Rρλµνvλ. (3.8)

This relation is called the Ricci identity and is not restricted to vectors but holds for tensors of arbitrary rank. From (3.5) and the definition of the Christoffel symbols (3.6), the Riemann tensor is explicitly expressed as a relation between the metric gµν and its gradients:

Rσρµν = ∂µΓσρν− ∂νΓσρµ+ ΓσλµΓλρν− ΓσλνΓλρµ. (3.9)

An arbitrary tensor of rank four defined on a n-dimensional manifold has n4 independent com-ponents. However, the Riemann tensor has several symmetries that reduces the number of independent components. The Riemann tensor is antisymmetric in the first two and last two lower indices and is symmetric in the interchange of these two pairs. Moreover, the antisym-metric part of the last three indices is zero. These symmetries together reduce the independent components of the Riemann tensor to n₁₂2(n2− 1). Therefore, by algebraic arguments, the Rie-mann tensor of a four-dimensional manifold has 20 independent components, while the RieRie-mann tensor of a three-dimensional manifold has six. In a similar fashion, one can show that the num-ber of independent components of the Ricci tensor Rµν = Rλµλν equals n₂(n + 1). In four

dimensions, we find that the number of independent components of the Ricci tensor counts ten. The ten remaining independent components of the Riemann tensor are captured by the Weyl tensor. Whereas the Ricci tensor encodes the part of the Riemann tensor that depends on the energy-momentum tensor, as in the Einstein equation, the Weyl tensor captures the traceless component of the Riemann tensor. In three dimensions, however, we find that the number of independent components of the Riemann tensor and Ricci tensor are equal and consequently the Weyl tensors vanishes.

The Einstein field equation (3.2) is a relation defined on a four-dimensional manifold and there-fore has ten independent components. However, not all of these independent components are

physical degrees of freedom. From the contracted Bianchi identities, one can show that four inde-pendent components are gauge degrees of freedom which arise from diffeomorphism invariance. Therefore, only six degrees of freedom are physical. In vacuum, when Tµν = 0, an additional

gauge freedom exists which reduces the physical degrees of freedom to two.

3.1.3 Black holes

A particular interesting solution of general relativity is the black hole. Black holes are regions of spacetime that are causally disconnected by an event horizon, which is a three-dimensional null hypersurface. The region inside this event horizon contains inevitably a curvature singularity, which makes the numerical integration of the Einstein field equations a non-trivial task [64]. Analytical solutions of the exterior region of the black hole have been found for several cases. The most generic stationary black hole spacetimes can be described by the Kerr-Newman metric, which is determined by only the mass M , angular momentum J and charge Q of the black hole. In this work, we consider two subclasses of the Kerr-Newman metric: the Schwarzschild metric (Q = J = 0) and the Kerr (Q = 0) metric.

Schwarzschild black holes

The Schwarzschild black hole is a static, non-spinning black hole that is spherically symmetric. In the Schwarzschild coordinates {t, r, θ, ϕ}, the metric is given by

ds2 = −1 −2M r  dt2+1 −2M r −1 dr2+ r2dΩ2, (3.10) where dΩ2:= dθ2+ sin2θdϕ2. From (3.10) we can infer various properties. First, Schwarzschild geometry in Schwarzschild coordinates admit two Killing vectors: the azimuthal Killing vector ∂ϕ and the timelike irrotational Killing vector ∂t, which demonstrates the stationary

charac-teristic of the geometry. Second, one can easily show that this metric has two singularities: at r = 0 and at r = 2M . The latter is a coordinate singularity that is located at the event horizon, which can be removed through other choices of coordinates, such as the Eddington-Finkelstein coordinates [65]. In this work we will also consider the Schwarzschild metric in isotropic coordi-nates. The isotropic coordinates is an attractive choice of coordinates because, unlike (3.10), the spatial component of the Schwarzschild solution is Euclidean. This is especially relevant for our discussion in the subsequent chapters, since in numerical relativity we often choose an ansatz in which the spatial metric is conformally flat. To convert (3.10) to isotropic coordinates, one has to rescale the radial coordinate as r = ¯r(1 + m/2¯r)2, where r is the Schwarzschild and ¯r is the

isotropic radial coordinate. We then find for the Schwarzschild metric in isotropic coordinates ds2 = − 2¯r − M 2¯r + M  dt2+  1 +M 2¯r 4 (d¯r2+ ¯r2dΩ2) , (3.11) The Schwarzschild metric in isotropic coordinates describes the spacetime outside the event horizon r ≥ 2M , which in isotropic coordinates is located at ¯r = M/2. Hence, Schwarzschild geomtery in isotropic coordinates has no singularity at the black hole’s horizon. Notice that the Schwarzschild radial coordinate r and isotropic coordinate ¯r asymptotically converge towards each other for distances far from the horizon, such that for large distances the metrics (3.10) and (3.11) become similar.

Kerr black holes

Recall from Chapter 2 that a necessary requirement for superradiance to occur is a rotating black hole. Only when the spin of the black hole is greater than the angular phase velocity of the ultralight boson field, superradiant instabilities can grow. We can describe the spacetime of an isolated and rotating black hole with zero charge with the Kerr metric. The Kerr metric in “Schwarzschild-like” Boyer-Lindquist coordinates (t, r, θ, ϕ) reads

ds2 = −∆ ρ2(dt − a sin 2_θdϕ)2₊ρ2 ∆dr 2_{+ ρ}2_dθ2₊sin2θ ρ2 (adt − (r 2_{+ a}2_)dϕ)2_{, (3.12)}

where we have defined

∆ ≡ r2− 2M r + a2, ρ2 ≡ r2+ a2cos2θ . (3.13) Here, M is the mass and a is the spin of the black hole. From this metric, we can immediately infer that the Kerr solution is stationary and axisymmetric around the z-axis, e.g. ∂tand ∂ϕ are

Killing vector. Moreover, the Kerr metric has similar to the Schwarzschild metric an intrinsic curvature singularity at r = 0. However, the horizon of the Kerr black hole is considerably more involved than in Schwarzschild spacetimes. For a 6= 0, the Kerr black hole has an inner (r+) and

outer (r−) horizon which are located at r±= M ±

√

M2_{− a}2_{. Note that these horizons merge}

when a → M .

For large distances (r  M ) from the horizons, one can show that the Kerr metric in Boyer-Lindquist coordinates shown in (3.12) reduces to the Schwarzschild metric. As analytical solu-tions of the scalar cloud only exists in these far regions, it should be sufficient to assume the Schwarzschild solution when comparing the numerical results.

3.2

Geometry of hypersurfaces

In this section we discuss the geometry and embedding of spacelike hypersurfaces in spacetime manifolds, needed to perform the 3+1 decomposition of the Einstein field equations and Klein-Gordon equation. The embedding of hypersurfaces in some larger manifold is mathematically well developed and is not restricted to Lorentzian spacetimes: hypersurfaces may be spacelike, timelike or even null. However, as the main aim of this work is to solve the initial value problem in a four-dimensional spacetime described by a Lorenzian manifold M, we restrict the discussion to spacelike hypersurfaces.

Before defining a hypersurface, we briefly discuss a general notion of submanifolds. As the name suggests, a submanifold of dimension m is a subset of another n-dimensional manifold and has therefore m ≤ n dimensions. The image Φ(S) of a submanifold S is said to be an embedded submanifold of M if the mapping Φ : S −→ M is one-to-one and has a continuous inverse map Φ−1 : Φ(S) −→ S [59]. In particular, the map Φ induces a push-forward mapping Φ∗ to carry

(p,0)-rank tensors from S to M and a push-backward mapping Φ∗ to carry (0,q)-rank tensors from M to S. An important consequence of the one-to-one property of the embedding Φ is that it forbids self-interactions of Φ(S) in M .

In view of the notion of submanifolds, we define the hypersurface Σ as a (n − 1)-dimensional sub-manifold of a n-dimensional sub-manifold that is embedded with the one-to-one and continuous map Φ. Hypersurfaces embedded in a two-dimensional manifold are plane curves while hypersurfaces embedded in a three-dimensional manifold are surfaces. Spacetime, which is described by a four-dimensional Lorentzian manifold M, embeds hypersurfaces Σ that are three-dimensional. Locally we can define a hypersurface as the set of points for which some scalar function t on M is constant [24]:

∀p ∈ M, p ∈ Σ ⇐⇒ t(p) = const. (3.14)

The hypersurface Σ is then said to be a level surface of t. In the 3+1 formalism of general relativity, the scalar function t is connected to the time-coordinate of the spacetime manifold M, while the spatial coordinates of M are represented by the coordinates of embedded hypersurfaces. For this reason, we limit our discussion in this chapter to spatial hypersurfaces.

3.2.1 Hypersurface embedding

To decompose four-dimensional tensors on (M, gµν) into tangent and normal components with

hy-persurfaces. Given a scalar field t on M, such that a hypersurface Σ is a level surface of t, we can define a vector ζµ _{on M which is orthogonal to Σ by taking the gradient of t:}

ζµ= gµν∇νt . (3.15)

By definition, ζµ is orthogonal to all vectors tangent to the hypersurface it is defined upon. For timelike spacetimes (M, gµν), and therefore timelike normal vectors such as ζµ, hypersurfaces

embedded in M are spacelike and have Riemannian metrics with signature (+,+,+) [59].3 The direction of the normal vector field ζµ is by definition unique and therefore collinear to any normal vector field on the hypersurface. Therefore, we can naturally define a unit normal vector nµ _{by normalizing the general normal vector (3.15):}

nµ= − ζ

µ

| − ζνζν|1/2

. (3.16)

The negative sign arises because nµ is timelike and points in the direction of increasing t. From (3.16), we find that nµnµ = −1.4 The unit normal vector can be regarded as the

four-velocity of a “normal” observer, which we will call an Eulerian oberver. Wordlines of Eulerian observers are always normal to spatial hypersurfaces and from the viewpoint of Eulerian ob-servers, the hypersurface Σ can (at least locally) be interpreted as a set of simultaneous events. We return to the notion of Eulerian observers in Section 3.3 when discussing foliations of hy-persurfaces.

Having defined the unit normal vector on the hypersurface Σ, we can construct the spatial metric γij and inverse metric γij of the hypersurface which are induced by gµν:

γµν = gµν+ nµnν, γµν = gµν+ nµnν. (3.17)

The inverse is obtained by raising both indices of the induced metric with the spacetime metric. From this definition, we can intuitively interpret the spatial metric γij as a tensor that measures

the spacetime distance with the spacetime metric gµν and eliminates the timelike contribution

with nµnν. To see γµνis indeed spatial, we can contract (3.17) with the normal unit vector (3.16):

nµγµν = nµgµν+ nµnµnν = nν− nν = 0 , (3.18)

3

If ζµ is spacelike, the hypersurface would be timelike and if ζµ is null, the hypersurface would be a null surface.

4

In case of a timelike hypersurface, the normal unit vector in (3.16) would lose the minus sign and nµnµ= +1

where we used nµnµ= −1 for spacelike hypersurfaces.

To decompose four-dimensional tensors on M into a spatial part tangent to Σ and a timelike part normal to Σ, we introduce two projection operators. To obtain the purely spatial part on T (Σ) of four-dimensional tensors, we need an coordinate-independent reverse mapping from M to Σ that carries both (p,0)-rank as (0,q)-rank tensors. To this end we construct the orthogonal projector operator as

γµν = gµλγλν = gµν + nµnν = δνµ+ nµnν, (3.19)

where δνµ is the Kronecker delta. An important property of the projection operator is that it is

idempotent. Noting that nλnλ= −1, one can show that

γ_λµγλν = (δ_λµ+ nµnλ)(δλν+ nλnν)

= δµ_ν + nµnν+ nµnν + nµnλnλnν

= γµν.

(3.20)

The projection operator allows us to write four-dimensional tensors on M as a three-dimensional tensor on Σ: 3_Ta1...am b1...bn = γ a1 µ1...γ am µmγ ν1 b1...γ νn bn 4_Tµ1...µm ν1...νn. (3.21)

Another noteworthy property of the projection operator is that it acts like a metric when acting on vectors that are tangent to Σt. Given a vector field vν on T (Σ), we find that

γµ_νvν = δ_νµvν+ nµnνvν = vµ, (3.22)

because nνvν = 0 for vν tangent to Σ, showing the same behaviour as the metric.

To decompose four-dimensional tensors in its timelike component, we simply contract a tensor with −nµnν = δνµ− γµν. With the normal and orthogonal projection operators we can write

tensors in spatial and timelike parts. Moreover, we can apply mixed projections of rank-n tensors of n > 1 by projecting indices both normally and orthogonal to Σ.

As metrics of spacelike and timelike manifolds are non-degenerate, there always exists a unique and torsion-free connection on that manifold which satisfies ∇ρgµν = 0 [61]. Similar to the

spacetime connection, we can define a metric compatible covariant derivative Dk associated to

covariant derivative ∇ρ associated with the spacetime metric gµν and the covariant derivative

Dk associated with the induced metric γij, we project the ordinary covariant derivative ∇ρ of

a (n,m)-rank tensor defined on the spacetime manifold (M, gµν) onto the hypersurface Σ. To

this end, we define the relation between the spacetime and hypersurface covariant derivatives as Di ≡ γµ_i∇µand apply the orthogonal projector (3.19) (n + m) times:

Di3Ta1...amb1...bn = γ a1 µ1...γ am µmγ ν1 b1...γ νn bnγ ρ i∇ρ4Tµ1...µnν1...νm. (3.23)

Hence, contracting the orthogonal projection operator with all indices of the four-dimensional tensor and covariant derivative results in a purely spatial object. Note that from the rela-tion 3.23) it immediately follows that Dkγij = 0. In the next sections, we discuss the Riemann

curvature tensor associated to the connection Dk and how it is related to the four-dimensional

Riemann curvature tensor of spacetime.

3.2.2 Curvature of hypersurfaces

For n-dimensional manifolds of n > 2, the Riemann curvature tensor is not necessarily van-ishing and we may speak of the manifold having curvature. Hypersurfaces embedded in four-dimensional spacetime are three-four-dimensional and can consequently have curvature themselves. The curvature of hypersurfaces is differently described by observers living on the hypersurface itself and observers living on the spacetime manifold M in which it is embedded. Vice versa, the curvature of spacetime is also differently described by observers living on the spacetime manifold and observers living on embedded hypersurfaces. Therefore, we describe the curvature of hypersurfaces with two different constructs: the intrinsic curvature and extrinsic curvature. The intrinsic curvature of a manifold specifies the curvature within the manifold seen from an observer living on that same manifold and is described by the Riemann curvature tensor asso-ciated to the manifold (see Section 3.1.2). Since hypersurfaces embedded in four-dimensional manifolds have dimension n > 2, the Riemann curvature tensor is not necessarily vanishing and we can define the intrinsic curvature of the hypersurfaces by the three-dimensional Riemann curvature tensor. Similar to the four-dimensional spacetime Riemann tensor, the Ricci-identity of the three-dimensional Riemann tensor3Rklij expresses the non-commutativity of two

succes-sive covariant derivatives Dk, cf. (3.8). Given a vector vkon a hypersurface Σ, the Ricci-identity

reads

Although 3Rklij and the corresponding Ricci tensor 3Rij = 3Rkikj and scalar curvature 3R =

γij 3_R

ij contain the intrinsic information of a hypersurface’s curvature, it does not provide any

information on how this intrinsic curvature is embedded in the larger spacetime manifold M. From the Ricci-identity in (3.24), we notice that the intrinsic curvature of a hypersurface is computed by spatial derivatives of the induced metric (which is purely spatial). The spacetime intrinsic curvature, however, is a four-dimensional tensor containing time-derivatives of the spacetime metric (c.f. (3.8)). Therefore, the intrinsic curvature of a spatial hypersurface does not carry the same information as the intrinsic curvature of the spacetime in which it is embedded. This missing information is carried by the extrinsic curvature. The extrinsic curvature can be defined in different contexts. From the viewpoint of an observer living on a single hypersurface Σ, the extrinsic curvature can be defined as the change of the direction of normal vectors on Σ when moving along the hypersurface to neighboring points on Σ. We therefore define the extrinsic curvature Kµν as the negative projection of the gradient of the unit normal vector [57]:5

Kµν = −γρµγσν∇ρnσ. (3.25)

As the unit normal vector nµ_{is rotation-free, one can show that the anti-symmetric part of (3.25)}

vanishes [66]. Hence, by definition, Kµνis a symmetric tensor and completely spatial. To expand

the right hand side of this definition, we substitute the definition of the orthogonal projection operator (3.19) twice, use the identity nµ∇νnµ = 0 and the relation between the spatial and

four-dimensional covariant derivatives from (3.23) to write the extrinsic curvature in terms of the unit normal vector:

Kµν = −∇µnν − nµnρ∇ρnν. (3.26)

From an observer living on the spacetime manifold M, the extrinsic curvature of a hypersurface Σ prescribes how the hypersurface is embedded in the larger space of M. Only by knowing both intrinsic and extrinsic curvature of the hypersurface, one can describe the relation between the four-dimensional and three-dimensional curvature.

By taking the trace of this relation with respect to the spacetime metric, we find a relation between the divergence of the unit normal vector and the trace of the extrinsic curvature K, called the mean curvature

K = gµνKµν = −∇µnµ. (3.27)

5

The negative sign is a common convention in numerical relativity, although some literature define the extrinsic curvature with a positive sign.

The mean curvature, which will appear in the 3+1 decomposition of the Einstein field equations, denotes the local curvature of the hypersurface. From (3.27), we may notice that the divergence of Eulerian observers is controllable when one imposes an appropriate condition on the mean curvature [67]. In Section 3.3.2, we briefly elaborate on this further by presenting an alternative expression of the mean curvature in the context of a foliation of multiple hypersurfaces. In Section 5.2.2 we discuss appropriate choices of the mean curvature in order to simplify our equations.

3.3

Foliations of spacetime

In the previous section, we have defined a hypersurface with associated metric, connection and normal vector. In addition, we described how a single hypersurface (Σ, γij) is embedded in a

four-dimensional manifold (M, gµν) by introducing the extrinsic curvature. To describe the entire

spacetime instead of a single time-slice, we need a set of hypersurfaces that includes the whole domain of M (at least locally). As spacetime respects causality, every spacelike hypersurface should intersect every inextendible causal curve (timelike or null) exactly once [68]. Such a hypersurface is called a Cauchy surface and spacetimes admitting a Cauchy surface are said to be globally hyperbolic [69].

Given a smooth and regular scalar field t on a globally hyperbolic spacetime (M, gµν), we can

formally define a foliation as a set of hypersurfaces of which every hypersurface Σt∈R is a level

surface of t of which no hypersurfaces Σt are intersecting [24]. We call t the time-function.

Except the requirements that the unit normal vector nµ ∝ ∂µt should be timelike and

future-oriented and that t should be a single-valued function of xµ, the time-function is completely arbitrary [66]. Moreover, a spacetime is said to by foliated by hypersurfaces Σtif the whole set

of hypersurfaces covers the spacetime (M, gµν):

M = [

t∈R

Σt. (3.28)

For timelike spacetimes, all hypersurfaces of the foliation are spacelike and therefore have time-like normal vectors and spacetime-like tangent vectors. On each individual hypersurface, a coordinate system (xi) is defined that coincides with the spatial coordinates of M, (xµ) = (t, xi).

First we introduce relations between the coordinates of successive hypersurfaces in a foliation in Section 3.3.1. We then discuss the time evolution of the induced metric γij and present an

alternative definition of the extrinsic curvature in the context of a foliation in Section 3.3.2. Fi-nally, in Section 3.3.3, we decompose the spacetime Riemann tensor, yielding the Gauss-Codazzi

equations and an evolution equation of the extrinsic curvature. These equations constitute the 3+1 decomposition of the Einstein field equations, which will be discussed in Chapter 4.

3.3.1 Lapse function and shift vector

Now that we have a mathematical description of individual hypersurfaces and defined a foliation of hypersurfaces, we can construct fields that describe how coordinates between neighbouring hypersurfaces in a foliation are related. To this end, we introduce two functions: a scalar function called the lapse function to prescribe how quantities evolve in the direction of the normal vector, and a vector function called the shift vector which describes their change in spatial directions. As we will show, we can use the lapse function and shift function to write the line element as

ds2 = −N2dt2+ γij(dxi+ βidt)(dxj+ βjdt) . (3.29)

From this metric, one can clearly see how the lapse function treats the time-evolution, while the shift vector handles the spatial evolution over time. To derive the metric in (3.29, we assume the spacetime manifold (M, gµν) is foliated by the time-function t(xµ) whose level-surfaces

correspond to the slices Σt of the foliation. We then introduce the time vector tµ on (M, gµν)

that is tangent to lines of constant spatial coordinates: tµ= (∂tt)µ= (1, 0, 0, 0). Notice that this

does not necessarily applies for nµ, since the spatial coordinates on each hypersurface are still gauge invariant. Normalizing tµsuch that tµ∇µt = 1, the time-function tµcan be interpreted as

the time flow vector of spacetime. To describe the propagation of coordinates in our 3+1 system, we decompose the time vector tµ _{orthogonal and tangent to the time slices Σ}

t. We now define

the lapse function N as the projection of the time flow vector tµ normal to the hypersurfaces Σt [58]:6

N = −gµνtµnν. (3.30)

The lapse function measures the rate of flow of proper time in the direction of the normal vector nµwith respect to the spacetime time-coordinate t. As the normal vectors of the foliated hypersurfaces are timelike, we may assume that the lapse function is strictly positive and non-vanishing [57]. For flat spacetimes, the lapse function is N = 1 and coincides with the spacetime time-coordinate t.

6_{In most numerical relativity literature the lapse function is denoted by α. However, we reserve α for the}

We define the shift vector βi as the spatial projection of the time flow vector tµ onto to the hypersurface Σt:

βi = γi_µtµ. (3.31) The shift vector measures how the local spatial coordinate system changes with respect to the tangent space of successive hypersurfaces. Notice that the shift vector is a spatial vector and therefore satisfies βµnµ= 0. Now that we have decomposed the time vector tµ, we can write it

in terms of its normal and tangent parts with respect to the hypersurface Σt:

tµ= N nµ+ βµ. (3.32) Choosing a foliation fixes the lapse function N and normal vector nµ[22]. The direction of the time-function tµ is then to be determined by the shift vector. Moreover, tµis tangent to lines of constant spatial coordinates, which are free to be chosen on each hypersurfaces. The shift vector βµ_{therefore encapsulates the spatial coordinate gauge of the hypersurfaces. Furthermore, notice}

that tµtµ= −N2+ βkβk. Hence, for timelike time-functions, we have βkβk < N2. By rewriting

equation (3.32), we can express the normal unit vector in terms of the time vector, the shift vector and the lapse function

nµ= 1 N(t

µ_{− β}µ_{) , (3.33)}

where we often write tµ= (∂t)µ. Substituting this expression into the definition of the (inverse)

spatial metric (3.17), we can write the (inverse) spacetime metric as:

gµν = γµν− nµ_nν _{= γ}µν₋ 1

N2(t

µ_{− β}µ_)(tν _{− β}ν_{) . (3.34)}

As the (inverse) induced metric and shift vector are spatial objects, we can write the spacetime (inverse) metric components explicitly as:

gµν = −N2_{+ β}k_β k βj βi γij ! , gµν = −N −2 _N−2_βj N−2βi γij − N−2βiβj ! . (3.35)

Hence, we have found the metric of the 3+1 formalism displayed in (3.29). Notice that by introducing the lapse function and shift vector, we have constructed a coordinate system that is naturally adapted to the foliation. In this coordinate system, we find that the unit normal vector