Gravitational wave extraction in higher dimensional numerical relativity using the Weyl tensor

Gravitational waves are one of the most important diagnostic tools in the analysis of strong-gravity dynamics and have been turned into an observational channel with LIGO's detection of GW150914. Aside from their importance in astrophysics, black holes and compact matter distributions have also assumed a central role in many other branches of physics. These applications often involve spacetimes with $D>4$ dimensions where the calculation of gravitational waves is more involved than in the four dimensional case, but has now become possible thanks to substantial progress in the theoretical study of general relativity in $D>4$. Here, we develop a numerical implementation of the formalism by Godazgar and Reall (Ref.[1]) -- based on projections of the Weyl tensor analogous to the Newman-Penrose scalars -- that allows for the calculation of gravitational waves in higher dimensional spacetimes with rotational symmetry. We apply and test this method in black-hole head-on collisions from rest in $D=6$ spacetime dimensions and find that a fraction $(8.19\pm 0.05)\times 10^{-4}$ of the Arnowitt-Deser-Misner mass is radiated away from the system, in excellent agreement with literature results based on the Kodama-Ishibashi perturbation technique. The method presented here complements the perturbative approach by automatically including contributions from all multipoles rather than computing the energy content of individual multipoles.


Introduction
Gravitational waves (GWs) entered the limelight with the recent detection of GW150914 [2] -soon followed by a second detection GW151226 [3] -which not only constitutes the first observation of a black-hole (BH) binary system, but also marks a true milestone in gravitational physics. This breakthrough has opened a quantitatively new path for measuring BH parameters [4,5], testing Einstein's theory of general relativity [6] and probing extreme astrophysical objects and their formation history [7], and substantially broadens the scope of multi-messenger astronomy [8]. GW modelling, however, finds important applications beyond the revolutionary field of GW astronomy. Many fundamental questions in general relativity in D = 4 and D > 4 spacetime dimensions concern the stability of strong-gravity sources (see [9][10][11][12][13][14][15][16] and references therein) in the context of cosmic censorship violation, the solutions' significance as physical objects or expanding our understanding of the strong-field regime of general relativity. For instance, GW extraction from numerical simulations of rapidly spinning Myers-Perry BHs demonstrates how excess angular momentum is shed in order to allow the BH to settle down into a more slowly rotating configuration [13]. GW emission also represents a channel for mass-energy loss in ultra-relativistic arXiv:1609.01292v1 [gr-qc] 5 Sep 2016 collisions that are studied in the context of the so-called TeV gravity scenarios that may explain the hierarchy problem in physics; for reviews see e.g. [17][18][19].
The numerical study of solutions to Einstein's equations has proven incredibly useful for understanding the behaviour of black holes and other compact objects. Most recently, the application of numerical relativity in the generation of gravitational waveform templates for GW data analysis [37,[39][40][41][42][43][44] contributed to the above mentioned detection of GW150914.
The calculation of GWs in higher dimensional relativity requires generalization of these techniques to D > 4. The extraction of the GW energy flux from the Landau-Lifshitz pseudotensor has been generalized straightforwardly to higher dimensions in [61,62]. An extension of the Regge-Wheeler-Zerilli-Moncrief formalism for perturbations of spherically symmetric background spacetimes is available in the form of the Kodama and Ishibashi (KI) formalism [63,64] and forms the basis of the wave extraction techniques developed in [65,66]. The main assumption there is that far away from the strong-field regime, the spacetime is perturbatively close to the Tangherlini [67] spacetime. The deviations from this background facilitate the construction of master functions according to the KI formalism which in turn provide the energy flux in the different (l, m) radiation multipoles. Even though both methods are in practice applied at finite extraction radius, their predictions have been found to agree within a ∼ 1 % error tolerance when applied to BH head-on collisions starting from rest in D = 5 [68]. Recent years have also seen considerable progress in the understanding of the peeling properties of the Weyl tensor; see [1,69] and references therein.
In particular, Godazgar & Reall [1] have performed a decomposition of the Weyl tensor in higher dimensions, and derived a generalisation of the Newman-Penrose formalism for wave extraction to D > 4. The numerical implementation of this formalism and probing the accuracy for a concrete example application is the subject of this paper.
For this purpose, we require a formulation of the higher dimensional Einstein equations suitable for numerical integration. For computational practicality, all higher dimensional time evolutions in numerical relativity have employed symmetry assumptions to reduce the problem to an effective "d+1" dimensional computation where typically d ≤ 3. This has been achieved in practice by either (i) imposing the symmetries directly on the metric line element [70], (ii) dimensional reduction by isometry of the Einstein field equations [71,72] or (iii) use of so-called Cartoon methods [13,[73][74][75][76]. In our implementation, we use the latter method combined with the Baumgarte-Shapiro-Shibata-Nakamura [77,78] (BSSN) formulation of the Einstein equations as expanded in detail in [79]. The relevant expressions for the GW computation, however, will be expressed in terms of the Arnowitt-Deser-Misner [80] (ADM) variables, and the formalism as presented here can be straightforwardly applied in other common evolution systems used in numerical relativity.
The paper is structured as follows. In Section 2 we describe the notation used in the remainder of this work. In Section 3 we recapitulate the key results of [1] which sets up the formalism. In Section 4 we put the formalism into a form compatible with the modified Cartoon dimensional reduction of our simulations. In Section 5 we describe the numerical set up used in our simulations, analyse the energy radiated in BH collisions in D = 6 and compare the predictions with literature results based on alternative wave extraction techniques.

Notation and Indices
The key goal of our work is to extract the GW signal from dynamical, asymptotically flat D dimensional spacetimes, i.e. manifolds M equipped with a metric g AB , A, B = 0, . . . , D − 1, of signature D − 2 that obeys the Einstein equations with vanishing cosmological constant Here, T AB is the energy momentum tensor which we assume to vanish in the wave zone, R AB and R denote the Ricci tensor and scalar, respectively, and we are using units where the gravitational constant and the speed of light G = c = 1. Furthermore, we define the Riemann tensor and Christoffel symbols according to the conventions of Misner, Thorne & Wheeler [81]. The ADM space-time decomposition as reformulated by York [82] writes the metric line element in the form where I, J = 1, . . . , D − 1 and α, β I denote the lapse function and shift vector, respectively, and γ IJ is the spatial metric that determines the geometry of hypersurfaces t = const. For this choice of coordinates and variables, the Einstein equations (2.1) result in one Hamiltonian and D − 1 momentum constraints as well as D(D − 1)/2 evolution equations cast into first-order-in-time form by introducing the extrinsic curvature K IJ through for a detailed review of this decomposition see [83,84].
In the remainder of this work we assume that the ADM variables are available. The BSSN formulation, for example, employs a conformally rescaled spatial metric γ IJ , the rescaled traceless extrinsic curvatureÃ IJ , a conformal factor χ, the trace of the extrinsic curvature K and contracted Christoffel symbolsΓ I associated with the conformal metric. These are related to the ADM variables by so that the ADM variables can be reconstructed straightforwardly at every time in the evolution. For other popular evolution systems such as the conformal Z4 system [85,86] or the generalized harmonic formulation [73,87] the ADM variables are obtained in similar fashion or directly from the spacetime metric components through Eqs. (2.2), (2.3). Many practical applications of higher dimensional numerical relativity employ symmetry assumptions that reduce the computational domain to an effective three dimensional spatial grid using the aforementioned techniques. The reasons are twofold: (i) the computational cost of simulations in D > 3 + 1 dimensions massively increase with any dimension added and (ii) the SO(D − d) rotational symmetry accomodates many applications of interest that can therefore be handled by relatively straightforward extensions of numerical codes originally developed for astrophysical systems in four spacetime dimensions.
In general, the number of effective (spatial) dimensions can range inside 1 ≤ d ≤ D−2, where d = 1 corresponds to spherical symmetry, i.e. SO(D−1) isometry, and the other extreme d = D − 2 corresponds to the axisymmetric case SO (2) with rotational symmetry about one axis. As already detailed in [79], axisymmetry, i.e. d = D − 2, represents a special case that imposes fewer constraints on the components of tensors and their derivatives and is therefore most conveniently handled numerically in a manner similar but not identical to the general case. We will discuss first in detail the case 1 ≤ d ≤ D − 3 and then describe the specific recipe for dealing with d = D − 2.
Probably the most important situation encountered in practical applications is that of an effective d = 3 dimensional spatial computational domain. Whenever the expressions developed in the remainder of this work for general d are not trivially obtained for the special choice d = 3, we shall therefore explicitly write down the d = 3 version in addition to the general case.
Let us start by considering a spacetime with at least two rotational Killing vectors which is the scenario discussed in Secs. 2 and 3 of Ref. [79]. For convenience, we will employ two specific coordinate systems adapted to this situation. The first is a set of Cartesian coordinates where middle Latin indices without (with) a caret range from i = 1, . . . , d (î = 1, . . . , d − 1), early Latin indices run from a = d + 1, . . . , D − 1, and rotational Killing vector fields are presumed to exist in all planes spanned by any two of the coordinates (z, w a ). The second is a system of spherical coordinates where Greek indices run from α = 2, . . . , D − 1. We use middle, upper case Latin indices to denote all spatial coordinates of either of these systems, so that X I = (xî, z, w a ) and Y I = (r, φ α ) with I = 1, . . . , D − 1. In the special case d = 3, we use the notation xî ≡ (x, y), so that Eq. (2.5) becomes We orient the Cartesian coordinates (2.5) such that they are related to the spherical coordinates (2.6) by Here φ D−1 ∈ [0, 2π], and all other φ α ∈ [0, π], and we have formally extended the w coordinates to also include (in parentheses in the equation) w i = x i which turns out to be convenient in the notation below in Sec. 4. Note that for the orientation chosen here, the x axis denotes the reference axis for the first polar angle rather than the z axis which more commonly plays this role for spherical coordinates in three spatial dimensions. For orientation among the different sets of indices, we conclude this section with a glossary listing index variables and their ranges as employed throughout this work.
• Upper case early Latin indices A, B, C, . . . range over the full spacetime from 0 to D − 1.
• Upper case middle Latin indices I, J, K, . . . denote all spatial indices, inside and outside the effective three dimensional computational domain, running from 1 to D − 1.
• Lower case middle Latin indices i, j, k, . . . denote indices in the spatial computational domain and run from 1 to d. For d = 3, we have x i = (x, y, z). • Lower case middle Latin indices with a caretî,ĵ, . . . run from 1 to d − 1 and represent the x i (without caret) excluding z. In d = 3, we write xî = (x, y). • Lower case early Latin indices a, b, c, . . . denote spatial indices outside the computational domain, ranging from d + 1 to D − 1. • Greek indices α, β, . . . denote all angular directions and range from 2 to D − 1.
• Greek indices with a caretα,β, . . . denote the subset of angular coordinates in the computational domain, i.e. range from 2, . . . , d. As before, a caret thus indicates a truncation of the index range. • Put inside parentheses, indices cover the same range but merely denote labels rather than tensor indices. • An index 0 denotes a contraction with the timelike normal to the foliation, rather than the time component, as detailed below in Section 4.1.1. • ∇ A denotes the covariant derivative in the full D dimensional spacetime, whereas D I denotes the covariant derivative on a spatial hypersurface and is calculated from the spatial metric γ IJ . • We denote by R with appropriate indices the Riemann tensor (or Ricci tensor/scalar) of the full D dimensional spacetime, and by R the spatial Riemann tensor (or Ricci tensor/scalar) calculated from γ IJ .

Theoretical Formalism
Our wave extraction from numerical BH simulations in D > 4 dimensions is based on the formalism developed by Godazgar & Reall [1]. In this section, we summarize the key findings and expressions from their work. The derivation [1] is based on the definition of asymptotic flatness using Bondi coordinates [26] (u, r, φ α ) where u is retarded time, r the radius and φ α are D − 2 angular coordinates covering the unit D − 2 sphere. A spacetime is asymptotically flat at future null infinity [88] if the metric outside a cylindrical world tube of finite radius can be written in terms of functions A(u, r, φ α ), B(u, r, φ α ), C(u, r, φ α ) as with det h αβ = det ω αβ where ω αβ is the unit metric on the D − 2 sphere. For an asymptotically flat spacetime h αβ can be expanded as [88] and the Bondi news function is obtained from this expansion as the leading-order correction h (1) αβ . In analogy with the D = 4 case, a null frame of vectors is constructed which is asymptotically given by * * The construction of the exact analog of the Kinnersley [89] tetrad in general spacetimes at finite radius is subject of ongoing research even in D = 4 (see e.g. [34,90]). In practice, the error arising from the use of an asymptotic form of the tetrad at finite extraction radii is mitigated by extracting at various radii and extrapolating to infinity [37] and we pursue this approach, too, in this work.
Note that all the tetrad vectors are real in contrast to the D = 4 dimensional case where the two vectors m (2) and m (3) are often written as two complex null vectors. Next, the components of the Weyl tensor are projected onto the frame (3.3) and the leading order term in the radial coordinate is extracted. Following [1], we denote this quantity by Ω and its components are given by Hereê β (α) denote a set of vectors forming an orthonormal basis for the unit metric ω αβ on the D − 2 sphere. In practice, this basis is constructed using Gram-Schmidt orthonormalisation starting with the radial unit vector.
As with the Newman-Penrose scalar Ψ 4 in the four dimensional case, we note that this is the contraction of the Weyl tensor with the ingoing null vector twice and two spatial vectors. Whereas in D = 4 the two polarisations of the gravitational waves are encoded in the real and imaginary parts of Ψ 4 , here Ω (α)(β) is purely real, with the α, β labels providing the different polarisations.
The final ingredient for extracting the energy radiated in GWs is the rate of change of the Bondi mass given by [88] By substituting in forḣ (1) αβ from the definition of Ω (α)(β) we obtain an expression for the mass loss.Ṁ where the notation (. . .) 2 implies summation over the (α), (β) labels inside the parentheses, and dω denotes the area element of the D − 2 sphere. In practice, we will apply Eq. (3.6) at constant radius r, therefore replace retarded time u with "normal" time t and start the integration at t = 0 rather than −∞, assuming that GWs generated prior to the start of the simulation can be neglected.

Modified Cartoon Implementation
The formalism summarized in the previous section is valid in generic D dimensional spacetimes with or without symmetries. We now assume that the spacetime under consideration obeys SO(D − d) isometry with 1 ≤ d ≤ D − 3, and will derive the expressions required for applying the GW extraction formalism of Sec. 3 to numerical simulations employing the modified Cartoon method. Throughout this derivation, we will make use of the expressions for scalars, vectors and tensors in spacetimes with SO(D − d) symmetry and the regularisation of their components at z = 0 as listed in Appendices A and B of Ref. [79]. The key result of these relations for our purposes is that the ADM variables α, β I , γ IJ , K IJ for a spacetime with SO(D − d) isometry can be expressed completely in terms of their d dimensional components β i , γ ij and K ij as well as two additional functions γ ww and K ww according to while the scalar α remains unchanged. From the viewpoint of numerical applications, the key relations of the procedure reviewed in Sec. 3 are Eqs. (3.4) and (3.6). The first provides Ω (α)(β) in terms of the Weyl tensor and the normal frame, and the second tells us how to calculate the mass loss from Ω (α)(β) . The latter is a straightforward integration conveniently applied as a post processing operation, so that we can focus here on the former equation. For this purpose, we first note that in practice wave extraction is performed in the wavezone far away from the sources. Even if the sources are made up of non-trivial energy matter fields, the GW signal is calculated in vacuum where the Weyl and Riemann tensors are the same. Our task at hand is then twofold: (i) calculate the Riemann tensor from the ADM variables and (ii) to construct a null frame. These two tasks are the subject of the remainder of this section. The ADM formalism is based on a space-time decomposition of the D dimensional spacetime manifold into a one-parameter family of spacelike hypersurfaces which are characterized by a future-pointing, unit normal timelike field n A . This normal field together with the projection operator allows us to split tensor fields into components tangential or orthogonal to the spatial hypersurfaces by contracting each tensor index either with n A or with ⊥ B A . For a symmetric rank (0,2) tensor, for example, we thus obtain the following three contributions The most important projections for our study are those of the Riemann tensor which are given by the Gauss-Codazzi relations used in the standard ADM splitting of the Einstein Equations (see e.g. [83]) where in the last line we used the fact that in vacuum R AC and, hence, its projection vanishes (note, however, that in general R AC = 0 even in vacuum). Furthermore CD K AB is the covariant derivative of the extrinsic curvature defined on the spatial hypersurface, with Christoffel symbols calculated from the induced metric γ AB . Equations (4.4)-(4.6) tell us how to reconstruct the full D dimensional Riemann tensor from D − 1 dimensional quantities defined on the spatial hypersurfaces which foliate our spacetime.
From this point on, we will use coordinates adapted to the (D − 1) + 1 split. In such coordinates, we can replace in Eqs. (4.4)-(4.6) the spacetime indices A, B, . . . on the left and right-hand side by spatial indices I, J, . . . while the time components of the spacetime Riemann tensor are taken into account through the contractions with the unit timelike normal n A and which we denote with the suffix 0 in (4.5), (4.6). Note that more than two contractions of the Riemann tensor with the timelike unit normal n A vanish by symmetry of the Riemann tensor.

The Riemann tensor in SO(D − 3) symmetry
The expressions given in the previous subsection for the components of the Riemann tensor are valid for general spacetimes with or without symmetries. In this section, we will work out the form of the components of the Riemann tensor in spacetimes with For this purpose we recall the Cartesian coordinate system X I = (xî, z, w a ) of Eq. (2.5), adapted to a spacetime that is symmetric under rotations in any plane spanned by two of the (z, w a ). We discuss in turn how the terms appearing on the right-hand sides of Eqs. (4.4)-(4.6) simplify under this symmetry. We begin with the components of the spatial Riemann tensor, given in terms of the spatial metric and Christoffel symbols by The rotational symmetry imposes conditions on the derivatives of the metric, the Christoffel symbols and the components of the Riemann tensor that are obtained in complete analogy to the derivation in Sec. 2.2 and Appendix A of [79]. We thus calculate all components of the Riemann tensor, where its indices can vary over the coordinates inside and outside the computational domain, and obtain For the right-hand side of Eq. (4.6), we also need the spatial Ricci tensor which is obtained from contraction of the Riemann tensor over the first and third index. In SO(D − d) symmetry, its non-vanishing components are Note that the last expression, γ ww R wuwu , does not involve a summation over w, but merely stands for the product of γ ww with the expression (4.15). The components of the extrinsic curvature are given by Eq. (4.1). Its derivative is directly obtained from the expressions (A.1)-(A.12) in Appendix A of [79] and can be written as Next, we plug the expressions assembled in Eqs.
With these expressions, we are able to calculate all components of the spacetime Riemann tensor directly from the ADM variables γ ij , γ ww , K ij and K ww and their spatial derivatives. There remains, however, one subtlety arising from the presence of terms containing explicit division by z. Numerical codes employing vertex centered grids need to evaluate these terms at z = 0. As described in detail in Appendix A, all the above terms involving division by z are indeed regular and can be rewritten in a form where this is manifest with no divisions by zero.

The null Frame
The null frame we need for the projections of the Weyl tensor consists of D unit vectors as given in Eq. We begin this construction with the D − 2 unit basis vectors on the D − 2 sphere, m A (α) , and recall for this purpose Eq. (2.8) that relates our spherical coordinates (r, φ α ) to the Cartesian (xî, z, w a ). The set of spatial vectors, although not yet in orthonormalised form, then consists of a radial vector denoted bym (1) and D − 2 angular vectorsm (α) whose components in Cartesian coordinates X I = (xî, z, w a ) on the computational domain w a = 0 are obtained through chain rulẽ , . . . , , . . . , . · · · · · · m A (α) = (0, | 0, . . . , 0 Note that under this procedure the components outside the computational domain of these vectors remain zero and can therefore be ignored. The final element of the null frame we need is the ingoing null vector, which we call k A . Given in [1] as ∂/∂u − 1 2 ∂/∂r asymptotically, we transform out of Bondi coordinates, sending (u, r) → (t, r) and furthermore use the freedom of rescaling this null vector by applying a constant factor of * √ 2 (4.52) Expressing the timelike unit normal field n A in terms of our gauge variables α, β I we find where β I = (β i , 0, . . . , 0), m I (1) = (m i (1) , 0, . . . , 0). This result provides the ingoing null vector for any choice of d and is the version implemented in the code.
vanishing w components and from Eqs. (4.22)-(4.35) we see that all components of the Riemann tensor where an odd number of indices is in the range a, b, . . . are zero. The only non-vanishing terms involving the Riemann tensor with off-domain indices a, b, . . ., therefore, have either four such indices or two and contain a Kronecker delta δ ab ; cf. Eqs. (4.23), (4.25), (4.28), (4.33). In consequence, the mixed components Ω (α)(a) = 0 and the purely off-domain components Ω (a)(b) ∝ δ ab . The list of all non-vanishing components Ω (α)(β) is then given by ) )

SO(2) symmetry
In the axisymmetric case d = D − 2 there exists only one w direction (off domain). As discussed in Section 4 of [79], we keep all tensor components as we would in the absence of symmetry, and the modified Cartoon method and, thus, the rotational symmetry only enters in the calculation of spatial derivatives in the w direction. For SO(2) symmetry, the extraction of gravitational waves therefore proceeds as follows.
• All components of the ADM metric and extrinsic curvature are extracted on the D − 2 dimensional computational domain. • The spatial Riemann tensor and its contractions are directly evaluated using Eq. (4.7) using the relations of Appendix C in [79] for off-domain derivatives.
Note that with the existence of more components of the Riemann tensor, more projections of the Weyl tensor now exist, specifically cross-terms such as Ω (2)(w) . This can be seen straightforwardly by using SO(2) modified Cartoon terms from appendix C of [79] and the expressions for the full and spatial Riemann tensor given in Eqs. (4.4) and (4.7). For example, we can see that a component such as R wijk is non-zero. This will contribute to terms of the form Ω (α)(w) . As already emphasized in [79], the key gain in employing the modified Cartoon method for simulating axisymmetric spacetimes does not lie in the elimination of tensor components, but in the dimensional reduction of the computational domain.

Numerical simulations
In the remainder of this work, we will implement the specific version of the wave extraction for the case d = 3 and D = 6 and simulate head-on collisions of equal-mass, non-spinning BHs starting from rest. We will calibrate the numerical uncertainties arising from the numerical discretisation of the equations (fourth order in space and time and second order at the outer and refinement boundaries), the use of large but finite extraction radii and also consider the dependency of the results on the initial separation of the BHs. This type of collisions has already been studied by Witek et al. [68] who calculate the GW energy using the Kodama-Ishibashi formalism, which enables us to compare our findings with their values.

Code infrastructure and numerical Set-up
We perform evolutions using the LEAN code [57,91] which is based on CACTUS [92,93] and uses CARPET [94,95] for mesh refinement. The Einstein equations are implemented in the BSSN formulation with the modified Cartoon method employed to reduce computational cost. For the explicit equations under the SO(D − 3) symmetry that we use, see Section 3.2 of [79] with parameter d = 3. Without loss of generality, we perform collisions along the x-axis, such that the centre-of mass is located at the origin of the grid, and impose octant symmetry. We specify the gauge in terms of the "1+log" and "Γ driver" conditions for the lapse function and shift vector (see e.g. [96]) according to with initial values α = 1, β i = 0. The BH initial data is calculated using the higher dimensional generalization of Brill-Lindquist data [97,98] given in terms of the ADM variables by where the summation over N and K extends over the multiple BHs and spatial coordinates, respectively, and X K N denotes the position of the N th BH. As mentioned above, we place the BHs on the x axis in the centre-of-mass frame, so that in the equal-mass case, we have X 1 = ±x 0 . Our initial configuration is therefore completely specified by the initial separation which we measure in units of the horizon radius R h of a single BH. The BH mass and the radius R h are related through the mass parameter µ by where A D−2 is the surface area of the unit (D − 2) sphere. The computational domain used for these simulations consists of a set of eight nested refinement levels which we characterize in terms of the following parameters: (i) the resolution h on the innermost level which gets coarser by a factor of two on each consecutive outer level, (ii) the size L of the domain which describes the distance of the outermost edge from the origin, and (iii) the resolution H on the refinement level where the gravitational waves are extracted.
For each simulation, we calculate the Ω (α)(β) on our three dimensional computational grid and project them onto a two dimensional array representing a spherical grid at fixed coordinate radius. The data thus obtained on the extraction sphere are inserted into Eq. (3.6). The Ω (α)(β) are scalars and so in our angular coordinate system do not depend on φ 4 , . . . , φ D−1 , so the integral over the sphere in (3.6) can be simplified: A final integration over time of the variableṀ then gives the total radiated energy.

Numerical Results
We begin our numerical study with an estimate of the uncertainty in our GW estimates arising from the discretisation of the equations. For this purpose, we have evolved two BHs initially located at at x = ±x 0 = ±4.0 R h using a computational grid of size L = 181 R h and three resolutions h 1 = R h /50.8, h 2 = R h /63.5 and h 3 = R h /76.2 which corresponds to H 1 = R h /2.12, H 2 = R h /2.65 and H 3 = R h /3.17 in the extraction zone.
We measure the radiated energy in units of the total ADM mass of the spacetime, which for Brill-Lindquist data is given by Eq. (5.4) with µ ≡ µ 1 + µ 2 , the mass parameters of the initial BHs. The radiated energy as a function of time is shown in the upper panel of Fig. 5.1. The radiation is almost exclusively concentrated within a window of ∆t ≈ 20 R h around merger. During the infall and the post-merger period, in contrast, E rad remains nearly constant. By comparing the high-resolution result with that obtained for the coarser grids, we can test the order of convergence. To leading order, the numerical result f h for some variable obtained at finite resolution h is related to the continuum limit solution f by f = f h + O(h n ), where n denotes the order of convergence. By evaluating the quotient we can then plot the two differences f h1 − f h2 and f h2 − f h3 and test whether their ratio is consistent with a given value n. The results for our study are shown in the lower panel of Fig. 5.1 which demonstrates that our numerical results converge at fourth order. The discretisation error of the total radiated energy is then obtained as the difference between the finite resolution result and that predicted by Richardson extrapolation (see upper panel in the figure). We obtain for the high-resolution case E rad = 8.19 × 10 −4 M ADM with a discretisation error of ∼ 0.4 %.
A second source of error arises from the extraction at finite radius. Following standard practice (see e.g. [37]), we estimate this uncertainty by extracting the GW energy at a set of seven or eight finite radii in the range 40 R h to 110 R h and extrapolating these values assuming a functional dependency where a is a coefficient obtained through the fitting of the numerical data. By applying this procedure, we estimate the uncertainty due to the extraction radius at 0.2 % at R ex = 110 R h and 0.4 % at R ex = 60 R h . Finally, we have measured the dependency of the total radiated energy on the initial separation of the BHs. In addition to the simulations discussed so far, we have performed high-resolution simulations placing the BHs at x 0 = ±7.8 R h and at ±12.8 R h . We have found very small variations at a level of 0.1 % in the radiated energy for these cases, well below the combined error budget of 0.6 % obtained above. Compared with collisions in D = 4 dimensions (see e.g. Table II in [57]), E rad shows significantly weaker variation with initial separation in D = 6. We attribute this to the more rapid fall-off of the force of gravity in higher dimensions leading to a prolonged but dynamically slow infall phase which generates barely any GWs.
In summary, we find the total energy radiated in gravitational waves in a head-on collision of two equal-mass, non-spinning BHs to be E rad = (8.19 ± 0.05) × 10 −4 M ADM , (5.8) in excellent agreement with the value (8.1 ± 0.4) × 10 −4 reported in [68] using the Kodama-Ishibashi formalism.

Conclusions
The extraction of gravitational waves from numerical simulations is one of the most important diagnostic tools in studying the strong-field dynamics of compact objects in four as well as higher dimensional spacetimes. In this work   The most common case of modeling higher dimensional spacetimes with rotational symmetries is the case of d = 3 effective spatial dimensions which allows for straightforward generalization of existing codes (typically developed for 3+1 spacetimes) and also accommodates sufficiently complex dynamics to cover most of the important applications of higher dimensional numerical relativity. We have, for this purpose, explicitly given the specific expressions of some of our relations for d = 3 where these are not trivially derived from their general counterparts. For testing the efficacy and accuracy of this method, we have applied the wave extraction to the study of equal-mass, non-spinning headon collisions of BHs starting from rest in D = 6 using d = 3. We find these collisions to radiate a fraction (8.19 ± 0.05) × 10 −4 of the ADM mass in GWs, in excellent agreement with a previous study [68] employing a perturbative extraction technique based on the Kodama-Ishibashi formalism. We find this energy to be essentially independent of the initial separation which we have varied from 8.0 to 15.6 and 25.6 times the horizon radius of a single BH. We attribute this result to the higher fall-off rate of the gravitational attraction in higher dimensions and the correspondingly slow dynamics during the infall stage.
We finally note that the Weyl tensor based wave extraction ideally complements the perturbative extraction technique of the Kodama-Ishibashi formalism. The latter provides the energy contained in individual (l, m) radiation multipoles but inevitably requires cutoff at some finite l. In contrast, the Ω (α)(β) facilitate calculation of the total radiation, but without multipolar decomposition. It is by putting both extraction techniques together, that we obtain a comprehensive description of the entire wave signal. Future applications include the stability of highly spinning BHs and their transition from unstable to stable configurations, the wave emission in evolutions of black rings and an extended study of higher dimensional BH collisions over a wider range of dimensionality D, initial boosts and with non-zero impact parameter. These studies require particularly high resolution to accurately model the rapid fall-off of gravity, especially for D 4, and are therefore beyond the scope of the present study. However, the foundation for analysing in detail the GW emission in these and many more scenarios is now available in as convenient a form as in the more traditional 3+1 explorations of numerical relativity. For the axisymmetric case d = D − 2, we only need to regularise terms appearing in the calculation of derivatives in the off-domain w direction. All these terms are given explicitly in Appendix C of [79], so that in the following we can focus exclusively on the additional terms appearing for 1 ≤ d ≤ D − 3, i.e. for spacetimes admitting two or more rotational Killing vector fields.
The treatment of these terms proceeds in close analogy to that of the BSSN equations in the modified Cartoon approach as described in detail in Appendix B of [79]. In contrast to that work, however, we will not be using the conformally rescaled metric of the BSSN equations, which satisfies the simplifying condition detγ = 1, and so certain regularised terms involving the inversion of the metric will differ from the expressions obtained for the BSSN equations.
We start with a brief summary of the techniques and the main assumptions we will use to regularise expressions: 1. Regularity: We require all tensor components and their derivatives to be regular when expressed in Cartesian coordinates. Under transformation to spherical coordinates this implies that tensors containing an odd (even) number of radial indices, i.e. z indices in our notation, contain exclusively odd (even) powers of z in a series expansion around z = 0. Using such a series expansion enables us to trade divisions by z for derivatives with respect to z. For example, for the z component of a vector field V , we obtain where we have introduced the symbol * = to denote equality in the limit z → 0.

Absence of conical singularities:
We require that the spacetime contain no conical singularity at the origin z = 0. For the implications of this condition, we consider the coordinate transformation from (xî, z, w d+1 , . . . , w a . . . , w D−1 ) to (xî, ρ, w d+1 , . . . , w a−1 , ϕ, w a+1 , . . . , w D−1 ). As no other w b , b = a, coordinates will enter into this discussion we shall refer to w a as w. In these coordinates we have that and the line element for vanishing dxî = 0 and dw b = 0, b = a, is given by Requiring the circumference to be the radius times 2π, we have that γ ϕϕ = ρ 2 γ ρρ . Substituting the above expressions and taking the limit z → 0, we obtain Taking the time derivative of this relation and using the definition of the extrinsic curvature, we find that likewise 3. Inverse metric: Various terms that we need to address contain factors of the inverse metric γ IJ . In the practical regularisation procedure, these terms are conveniently handled by writing expressing γ IJ in terms of the downstairs metric components γ ij and γ ww which are the fields we evolve numerically. We know the metric takes the following form: and we shall denote the upper left quadrant by the matrix M ij . For simplicity, we will use the indexî to denote xî in this section, so e.g. cofactors C 12 = C x 1 x 2 and C 1z = C x 1 z , and similarly indices i, j, ... will denote the same range, but including the z component.
We can now denote the cofactor of an element in the top left quadrant of γ IJ as where η = D − d − 1 and the notation det(M kl{k =j,l =i} ) denotes the determinant of the matrix M kl obtained by crossing out the j th row and i th column. Likewise, we may add further inequalities inside the braces to denote matrices obtained by crossing out more than one row and column. We can then insert this expression for C ij and the determinant of the right hand side of Eq. (A.7), to obtain expressions for inverse metric components according to For d = 3, this procedure starts from the spatial metric The components C ij of the cofactor matrix (which is symmetric) are given by · · · · · · Czz = γ n ww (γxxγyy − γ 2 xy ) , and the inverse metric follows by inserting these into Eq. (A.10).
Using these techniques, we can regularise all terms in Eqs. (4.11), (4.12), (4.15), (4.21) and (4.29) that contain divisions by z. It turns out to be convenient to combine the individual terms into the following six expressions.
(1) δ i z − γ zi γ ww z We express γ zi in terms of the metric, and trade divisions by z for derivatives ∂ z and obtain For the d = 3 case this reduces to Here we simply trade divisions by z for ∂ z and obtain We use γ zz − γ ww * = O(z 2 ) and trade a division by z for a z derivative. The result is Using Eqs. (A.7)-(A.10), we express the inverse metric components γ zj in terms of the downstairs metric and trade the division by z for a z derivative. We thus obtain which in the case d = 3 reduces to γ ww γ zj ∂ j γ ww z * = (γ yx ∂ z γ zy − γ yy ∂ z γ zx )∂ x γ ww + (γ yx ∂ z γ zx − γ xx ∂ z γ yz )∂ y γ ww γ xx γ yy − γ 2 xy + ∂ z ∂ z γ ww . (A.19) The regularisation of this term proceeds in analogy to that of term (9) in Appendix B of [79], except we do not set det γ = 1. By rewriting 1 = γ zz /γ zz = γ zz det γ IJ /C zz , trading divisions by z for z derivatives and using γ zz * = γ ww + O(z 2 ), we obtain d−1 m=1 (−1)î +m−1 ∂ z γ zî ∂ z γ zm det(M kl{k =î,k =z,l =z,l =m} ) det(M kl{k =z,l =z} ) .
(A. 20) which in the case d = 3 reduces to The division by z is again traded for a derivative if i = z and for i = z, we use K zz = K ww + O(z 2 ), so that Appendix B. Normalisation of the spatial normal frame vectors In this section, we discuss how the set of spatial normal frame vectors, Eq. (4.38), can be recast in a form suitable for applying Gram-Schmidt orthonormalisation. It turns out to be convenient to first rescale them (α) such that they would acquire unit length in a flat spacetime with spatial metric δ IJ . Denoting these rescaled vectors with a caret, we havê , α = 2, . . . , D − 1 .
(B.1) Recall that we formally set w 1 ≡ x 1 , . . . , w d−1 ≡ x d−1 , w d ≡ z. As a convenient shorthand, we define We can now express the angles φ α in terms of the radial variables ρ I , Using these relations in (B.3), we obtain