Extraction of gravitational-wave energy in higher dimensional numerical relativity using the Weyl tensor

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 qualitatively 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–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 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–19].

The calculation of GW signals in the theoretical modelling of D  =  4 dimensional sources in the framework of general relativity has been increasingly well understood following seminal work by Pirani, Bondi, Sachs and others in the 1950s and 1960s; see e.g. [20–26] and [27] for a review. Applications are now routinely found in numerical and (semi-)analytic calculations [28–37] even though care needs to be taken when applied to numerical simulations on finite domains [38].

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–44] contributed to the above mentioned detection of GW150914.

The wave extraction techniques presently used in numerical simulations of astrophysical GW sources can be classified as follows: perturbative methods based on the formalism developed by Regge, Wheeler, Zerilli and Moncrief [9, 10, 45]; application of the quadrupole formula [46] in matter simulations [47, 48]; a method using the Landau–Lifshitz pseudo-tensor [49, 50]; Cauchy characteristic extraction [51–53]; and, probably the most popular technique, using the Weyl scalars from the Newman–Penrose formalism [24, 54–60]. A comparison of various of these techniques is given in [36] and an extended review in [61].

The calculation of GWs in higher dimensional relativity requires generalisation of these techniques to D  >  4. The extraction of the GW energy flux from the Landau–Lifshitz pseudotensor has been generalised straightforwardly to higher dimensions in [62, 63]. 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 [64, 65] and forms the basis of the wave extraction techniques developed in [66, 67]. The main assumption there is that far away from the strong-field regime, the spacetime is perturbatively close to the Tangherlini [68] 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 [69]. Recent years have also seen considerable progress in the understanding of the peeling properties of the Weyl tensor; see [1, 70] and references therein.

In particular, Godazgar and 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. This analysis provides us with a quantity analogous to the Weyl scalar ${{ \Psi }_{4}}$, from which we can calculate the energy radiated in gravitational waves in a similar fashion to the method in D  =  4. The one qualitative difference between the D  =  4 and D  >  4 cases comes in the availiability of a mode decomposition of the gravitational wave. In D  =  4 we can project the Weyl scalar onto spin weighted spherical harmonics, due to the decoupling of the equations of motion as shown by Teukolsky [71, 72]. In higher dimensions however, a set of conditions identified as sufficient for decoupling are not satisfied in black hole spacetimes [73], and so at present we are unable to project out the angular dependence of the gravitational wave energy. 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\leqslant 3$. This has been achieved in practice by either (i) imposing the symmetries directly on the metric line element [74], (ii) dimensional reduction by isometry of the Einstein field equations [75, 76] or (iii) use of so-called Cartoon methods [13, 77–80]. In our implementation, we use the latter method combined with the Baumgarte–Shapiro–Shibata–Nakamura [81, 82] (BSSN) formulation of the Einstein equations as expanded in detail in [83]. The relevant expressions for the GW computation, however, will be expressed in terms of the Arnowitt–Deser–Misner [84] (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.

The key goal of our work is to extract the GW signal from dynamical, asymptotically flat D dimensional spacetimes, i.e. manifolds $\mathcal{M}$ equipped with a metric g_AB, $A,\,B=0,\,\ldots,\,D-1$, of signature D  −  2 that obey the Einstein equations with vanishing cosmological constant

$$\begin{eqnarray}&&{{G}_{AB}}={{R}_{AB}}-\frac{1}{2}R{{g}_{AB}}=8\pi {{T}_{AB}}.\end{eqnarray} \tag{ 2.1 }$$

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 and Wheeler [85].

The ADM space-time decomposition as reformulated by York [86] writes the metric line element in the form

$$\begin{eqnarray}&&\text{d}{{s}^{2}}={{g}_{AB}}\text{d}{{x}^{A}}\text{d}{{x}^{B}}=\left(-{{\alpha}^{2}}+{{\beta}_{I}}{{\beta}^{I}}\right)\text{d}{{t}^{2}}+2{{\beta}_{I}}\text{d}{{x}^{I}}\text{d}t+{{\gamma}_{IJ}}\text{d}{{x}^{I}}\text{d}{{x}^{J}},\end{eqnarray} \tag{ 2.2 }$$

where $I,\,J=1,\,\ldots,\,D-1$, α and ${{\beta}^{I}}$ denote the lapse function and shift vector, respectively, and ${{\gamma}_{IJ}}$ is the spatial metric that determines the geometry of hypersurfaces $t=\text{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

$$\begin{eqnarray}&&{{\partial}_{t}}{{\gamma}_{IJ}}={{\beta}^{M}}{{\partial}_{M}}{{\gamma}_{IJ}}+{{\gamma}_{MJ}}{{\partial}_{I}}{{\beta}^{M}}+{{\gamma}_{IM}}{{\partial}_{J}}{{\beta}^{M}}-2\alpha {{K}_{IJ}}\,;\end{eqnarray} \tag{ 2.3 }$$

for a detailed review of this decomposition see [87, 88].

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 ${{\tilde{\gamma}}_{IJ}}$, the rescaled traceless extrinsic curvature ${{\tilde{A}}_{IJ}}$, a conformal factor χ, the trace of the extrinsic curvature K and contracted Christoffel symbols ${{\rm{\tilde \Gamma }}^{I}}$ associated with the conformal metric. These are related to the ADM variables by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \chi ={{\gamma}^{-1/(D-1)}},\gamma =\det {{\gamma}_{IJ}},~~~ & K={{\gamma}^{MN}}{{K}_{MN}}, \\ {{{\tilde{\gamma}}}_{IJ}}=\chi {{\gamma}_{IJ}}~~ & \Leftrightarrow ~~{{{\tilde{\gamma}}}^{IJ}}=\frac{1}{\chi}{{\gamma}^{IJ}}, \\ {{{\tilde{A}}}_{IJ}}=\chi \left({{K}_{IJ}}-\frac{1}{D-1}{{\gamma}_{IJ}}K\right)~~ & \Leftrightarrow ~~{{K}_{IJ}}=\frac{1}{\chi}\left({{{\tilde{A}}}_{IJ}}+\frac{1}{D-1}{{{\tilde{\gamma}}}_{IJ}}K\right), \\ {{{\rm{\tilde \Gamma }}}^{I}}={{{\tilde{\gamma}}}^{MN}}\rm{\tilde \Gamma }_{MN}^{I}, & \end{array}\end{eqnarray} \tag{ 2.4 }$$

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 [89–91] or the generalised harmonic formulation [77, 92] the ADM variables are obtained in similar fashion or directly from the spacetime metric components through equations (2.2) and (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 two-fold: (i) the computational cost of simulations in D  >  3  +  1 dimensions massively increases with any dimension added and (ii) the SO(D  −  d) rotational symmetry accommodates 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. The so-called cartoon methods have been developed for this specific purpose. An illustration of the idea is given in figure 1, where rotational symmetry is assumed in every plane spanned by two of the coordinates (z, w^a), where $a=d+1,\,\ldots,\,D-1$. Knowledge of the tensorial components in the d dimensional hypersurface spanned by $\left({{x}^{{\hat{i}}}},\,z\right),\,\hat{i}=1,\,\ldots,\,d-1$ then is sufficient to describe the entire spacetime. In numerical evolutions, however, there remains the question of evaluating derivatives perpendicular to that plane. In the original cartoon version [93], this was achieved through ghost zones off the computational domain which are populated by rotation of the data on the plane and then facilitate evaluation of the derivatives through standard discretisation methods. A modification of this method introduced in [77] (see also [13, 79, 80]), dispenses with the need for ghost zones. Instead, the required derivatives are obtained through analytic relations from derivatives in the on-domain directions. A detailed discussion and some example cases for the derivation of these relations are given in [83].

Zoom In Zoom Out Reset image size

In general, this method reduces the number of effective (spatial) dimensions from D  −  1 to d, where $1\leqslant d\leqslant D-2$. 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 [83], axisymmetry 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\leqslant d\leqslant 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 converted to d  =  3, we shall 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 sections 2 and 3 of [83]. For convenience, we will employ two specific coordinate systems adapted to this situation. The first is a set of Cartesian coordinates

$$\begin{eqnarray}&&{{X}^{A}}=\left(t,\,\underset{(d-1)\times}{{\underbrace{{{{x}^{1}},\,\ldots,{{x}^{d-1}}}}\,}}\,,\,z,\,\underset{(D-d-1)\times}{{\underbrace{{{{w}^{d+1}},\,\ldots,\,{{w}^{D-1}}}}\,}}\,\right)=\left(t,\,{{x}^{{\hat{i}}}},\,z,\,{{w}^{a}}\right)=\left(t,\,{{x}^{i}},\,{{w}^{a}}\right),\end{eqnarray} \tag{ 2.5 }$$

where middle Latin indices without (with) a caret range from $i=1,\,\ldots,\,d$ ($\hat{i}=1,\,\ldots,\,d-1$) and early Latin indices run from $a=d+1,\ldots,D-1$; see figure 1. The second is a system of spherical coordinates

$$\begin{eqnarray}&&{{Y}^{A}}=\left(t,\,r,\,\underset{(D-2)\times}{{\underbrace{{{{\phi}^{2}},\,{{\phi}^{3}},\,\ldots,\,{{\phi}^{D-1}}}}\,}}\,\,\right)=\left(t,\,r,\,{{\phi}^{\alpha}}\right),\end{eqnarray} \tag{ 2.6 }$$

where Greek indices run from $\alpha =2,\ldots,D-1$. We use middle, upper case Latin indices to denote all spatial coordinates of either of these systems, so that ${{X}^{I}}=\left({{x}^{{\hat{i}}}},\,z,\,{{w}^{a}}\right)$ and ${{Y}^{I}}=\left(r,\,{{\phi}^{\alpha}}\right)$ with $I=1,\ldots,D-1$. Rotational symmetry is assumed in all directions ${{\phi}^{a}},\,a=d+1,\,\ldots,\,D-1$. In the special case d  =  3, we use the notation ${{x}^{{\hat{i}}}}\equiv (x,\,y)$, so that equation (2.5) becomes

$$\begin{eqnarray}&&{{X}^{A}}=\left(t,\,x,\,y,\,z,\,{{w}^{4}},\,\ldots,\,{{w}^{D-1}}\right).\end{eqnarray} \tag{ 2.7 }$$

We orient the Cartesian coordinates (2.5) such that they are related to the spherical coordinates (2.6) by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \left({{w}^{1}}\equiv \right)\,{{x}^{1}} & =r\cos {{\phi}^{2}}, \\ \left({{w}^{2}}\equiv \right)\,{{x}^{2}} & =r\sin {{\phi}^{2}}\,\cos {{\phi}^{3}}, \\ \vdots & \\ \left({{w}^{d-1}}\equiv \right)\,{{x}^{d-1}} & =r\sin {{\phi}^{2}}\,\ldots \,\sin {{\phi}^{d-1}}\,\cos {{\phi}^{d}}, \\ \left({{w}^{d}}\equiv \right)\,z & =r\sin {{\phi}^{2}}\,\ldots \,\sin {{\phi}^{d-1}}\,\sin {{\phi}^{d}}\,\cos {{\phi}^{d+1}}, \\ {{w}^{d+1}} & =r\sin {{\phi}^{2}}\,\ldots \,\sin {{\phi}^{d-1}}\,\sin {{\phi}^{d}}\,\sin {{\phi}^{d+1}}\,\cos {{\phi}^{d+2}}, \\ \vdots & {} \\ {{w}^{D-3}} & =r\sin {{\phi}^{2}}\,\ldots \,\sin {{\phi}^{D-3}}\,\cos {{\phi}^{D-2}}, \\ {{w}^{D-2}} & =r\sin {{\phi}^{2}}\,\ldots \,\sin {{\phi}^{D-3}}\,\sin {{\phi}^{D-2}}\,\cos {{\phi}^{D-1}}, \\ {{w}^{D-1}} & =r\sin {{\phi}^{2}}\,\ldots \,\sin {{\phi}^{D-3}}\,\sin {{\phi}^{D-2}}\,\sin {{\phi}^{D-1}}. \end{array}\end{eqnarray} \tag{ 2.8 }$$

Here ${{\phi}^{D-1}}\in [0,2\pi ]$, and all other ${{\phi}^{\alpha}}\in [0,\pi ]$, 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 section 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,\,\ldots$ range over the full spacetime from 0 to D  −  1.
• Upper case middle Latin indices $I,\,J,\,K,\,\ldots$ denote all spatial indices, inside and outside the effective computational domain, running from 1 to D  −  1.
• Lower case middle Latin indices $i,\,j,\,k,\,\ldots$ 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 $\hat{i},\,\hat{j},\,\ldots$ run from 1 to d  −  1 and represent the xⁱ (without caret) excluding z. In d  =  3, we write ${{x}^{{\hat{i}}}}=(x,\,y)$.
• Lower case early Latin indices $a,\,b,\,c,\,\ldots$ denote spatial indices outside the computational domain, ranging from d  +  1 to D  −  1.
• Greek indices $\alpha,\,\beta,\,\ldots$ denote all angular directions and range from 2 to D  −  1.
• Greek indices with a caret $\hat{\alpha},\hat{\beta},\,\ldots$ denote the subset of angular coordinates in the computational domain, i.e. range from $2,\,\ldots,\,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.
• ${{\nabla}_{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 ${{\gamma}_{IJ}}$.
• We denote by R with appropriate indices the Riemann tensor (or Ricci tensor/scalar) of the full D dimensional spacetime, and by $\mathcal{R}$ the spatial Riemann tensor (or Ricci tensor/scalar) calculated from ${{\gamma}_{IJ}}$.

Our wave extraction from numerical BH simulations in D  >  4 dimensions is based on the formalism developed by Godazgar and Reall [1]. In this section, we summarise the key findings and expressions from their work.

The derivation [1] is based on the definition of asymptotic flatness using Bondi coordinates [26] $\left(u,\mathfrak{r},{{\phi}^{\alpha}}\right)$ where u is retarded time, $\mathfrak{r}$ the radius and ${{\phi}^{\alpha}}$ are D  −  2 angular coordinates covering the unit D  −  2 sphere. A spacetime is asymptotically flat at future null infinity [94] if the metric outside a cylindrical world tube of finite radius can be written in terms of functions $\mathcal{A}\left(u,\mathfrak{r},{{\phi}^{\alpha}}\right)$, $\mathcal{B}\left(u,\mathfrak{r},{{\phi}^{\alpha}}\right)$, $\mathcal{C}\left(u,\mathfrak{r},{{\phi}^{\alpha}}\right)$ as

$$\begin{eqnarray}&&\text{d}{{s}^{2}}=-\mathcal{A}{{\text{e}}^{\mathcal{B}}}\text{d}{{u}^{2}}-2{{\text{e}}^{\mathcal{B}}}\text{d}u\,\text{d}\mathfrak{r}+{{\mathfrak{r}}^{2}}{{h}_{\alpha \beta}}\left(\text{d}{{\phi}^{\alpha}}+{{\mathcal{C}}^{\alpha}}\text{d}u\right)\left(\text{d}{{\phi}^{\beta}}+{{\mathcal{C}}^{\beta}}\text{d}u\right),\end{eqnarray} \tag{ 3.1 }$$

with $\det {{h}_{\alpha \beta}}=\det {{\omega}_{\alpha \beta}}$ where ${{\omega}_{\alpha \beta}}$ is the unit metric on the D  −  2 sphere. For an asymptotically flat spacetime ${{h}_{\alpha \beta}}$ can be expanded as [94]

$$\begin{eqnarray}&&{{h}_{\alpha \beta}}={{\omega}_{\alpha \beta}}\left({{\phi}^{\gamma}}\right)+\underset{s\geqslant 0}{\sum}\,\frac{h_{\alpha \beta}^{(s+1)}\left(u,{{\phi}^{\gamma}}\right)}{{{\mathfrak{r}}^{D/2+s-1}}},\end{eqnarray} \tag{ 3.2 }$$

and the Bondi news function is obtained from this expansion as the leading-order correction $h_{\alpha \beta}^{(1)}$.

In analogy with the D  =  4 case, a null frame of vectors is constructed which is asymptotically given by⁴

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} l=-\frac{\partial}{\partial \mathfrak{r}},~\,~k=\frac{\partial}{\partial u}-\frac{1}{2}\frac{\partial}{\partial \mathfrak{r}},~\,~{{m}_{(\alpha )}}=\frac{\partial}{\partial {{\phi}^{\alpha}}}. \end{array}\end{eqnarray} \tag{ 3.3 }$$

Note that all the tetrad vectors are real in contrast to the D  =  4 dimensional case where the two vectors m₍₂₎ and m₍₃₎ 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 ${{ \Omega }^{\prime}}$ and its components are given by

$$\begin{eqnarray} &&\Omega _{(\alpha )(\beta )}^{\prime}\equiv {{C}_{ABCD}}{{k}^{A}}m_{(\alpha )}^{B}{{k}^{C}}m_{(\beta )}^{D}=-\frac{1}{2}\frac{\hat{e}_{(\alpha )}^{\mu}\hat{e}_{(\beta )}^{\nu}\ddot{{h}}\,_{\mu \nu}^{(1)}}{{{\mathfrak{r}}^{D/2-1}}}+\mathcal{O}\left({{\mathfrak{r}}^{-D/2}}\right).\end{eqnarray} \tag{ 3.4 }$$

Here $\hat{e}_{(\alpha )}^{\beta}$ denote a set of vectors forming an orthonormal basis for the unit metric ${{\omega}_{\alpha \beta}}$ 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 ${{ \Psi }_{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 ${{ \Psi }_{4}}$, here $\Omega _{(\alpha )(\beta )}^{\prime}$ is purely real, with the $\alpha,\,\beta$ 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 [94]

$$\begin{eqnarray}&&\overset{\centerdot}{{M}}\,(u)=\frac{1}{32\pi}{{{\int}^{}}_{{{S}^{D-2}}}}\overset{\centerdot}{{h}}\,_{\alpha \beta}^{(1)}{{\overset{\centerdot}{{h}}\,}^{(1)\alpha \beta}}\text{d}\omega.\end{eqnarray} \tag{ 3.5 }$$

By substituting in for $\overset{\centerdot}{{h}}\,_{\alpha \beta}^{(1)}$ from the definition of $\Omega _{(\alpha )(\beta )}^{\prime}$ we obtain an expression for the mass loss.

$$\begin{eqnarray}&&\overset{\centerdot}{{M}}\,(u)=-\underset{\mathfrak{r}\to \infty}{{\lim}}\,\frac{{{\mathfrak{r}}^{D-2}}}{8\pi}{{{\int}^{}}_{{{S}^{D-2}}}}{{\left({\int}_{-\infty}^{u}{{ \Omega }^{\prime}}_{(\alpha )(\beta )}^{{}}\left(\tilde{{u}}\,,\mathfrak{r},{{\phi}^{\gamma}}\right)\text{d}\tilde{{u}}\,\right)}^{2}}\text{d}\omega,\end{eqnarray} \tag{ 3.6 }$$

where the notation ${{(\ldots )}^{2}}$ implies summation over the $(\alpha ),\,(\beta )$ labels inside the parentheses, and $\text{d}\omega$ denotes the area element of the D  −  2 sphere. In practice, we will apply equation (3.6) at constant radius $\mathfrak{r}$, therefore replace retarded time u with 'normal' time t and start the integration at t  =  0 rather than $-\infty$, assuming that GWs generated prior to the start of the simulation can be neglected.

The formalism summarised 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\leqslant d\leqslant D-3$, and will derive the expressions required for applying the GW extraction formalism of section 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 [83]. The key result of these relations for our purposes is that the ADM variables α, ${{\beta}^{I}}$, ${{\gamma}_{IJ}}$, K_IJ for a spacetime with SO(D  −  d) isometry can be expressed completely in terms of their d dimensional components ${{\beta}^{i}}$, ${{\gamma}_{ij}}$ and K_ij as well as two additional functions ${{\gamma}_{ww}}$ and K_ww according to

$$\begin{eqnarray}&&{{\beta}^{I}}=\left({{\beta}^{i}},0\right),\\ {{\gamma}_{IJ}}=\left(\begin{array}{*{35}{l}} {{\gamma}_{ij}} & 0 \\ 0 & {{\delta}_{ab}}{{\gamma}_{ww}} \end{array}\right),\\ {{K}_{IJ}}=\left(\begin{array}{*{35}{l}} {{K}_{ij}} & 0 \\ 0 & {{\delta}_{ab}}{{K}_{ww}} \end{array}\right),\end{eqnarray} \tag{ 4.1 }$$

while the scalar α remains unchanged.

From the viewpoint of numerical applications, the key relations of the procedure reviewed in section 3 are equations (3.4) and (3.6). The first provides $\Omega _{(\alpha )(\beta )}^{\prime}$ in terms of the Weyl tensor and the normal frame, and the second tells us how to calculate the mass loss from $\Omega _{(\alpha )(\beta )}^{\prime}$. 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 wave zone 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.

4.1. The Riemann tensor

4.1.1. The (D  −  1)  +  1 splitting of the Riemann tensor.

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 characterised by a future-pointing, unit normal timelike field n^A. This normal field together with the projection operator

$$\begin{eqnarray}&&{{\bot}^{A}}{{}_{B}}={{\delta}^{A}}{{}_{B}}+{{n}^{A}}{{n}_{B}},\end{eqnarray} \tag{ 4.2 }$$

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 ${{\bot}^{B}}{{}_{A}}$. For a symmetric rank (0,2) tensor, for example, we thus obtain the following three contributions

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{T}_{00}}\equiv {{T}_{AB}}{{n}^{A}}{{n}^{B}},~~~\bot {{T}_{0A}}=\bot {{T}_{A0}}\equiv {{\bot}^{C}}{{}_{A}}{{T}_{CB}}{{n}^{B}},~~~\bot {{T}_{AB}}\equiv {{\bot}^{C}}{{}_{A}}{{\bot}^{D}}{{}_{B}}{{T}_{CD}}. \end{array}\end{eqnarray} \tag{ 4.3 }$$

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. [87])

$$\begin{eqnarray}\begin{array}{*{35}{l}} \bot {{R}_{ABCD}}= & {{\mathcal{R}}_{ABCD}}+{{K}_{AC}}{{K}_{BD}}-{{K}_{AD}}{{K}_{CB}}, \end{array}\end{eqnarray} \tag{ 4.4 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \bot {{R}_{A0CD}}\equiv \bot \left({{R}_{ABCD}}{{n}^{B}}\right)= & -{{D}_{C}}{{K}_{AD}}+{{D}_{D}}{{K}_{AC}}, \end{array}\end{eqnarray} \tag{ 4.5 }$$

$$\begin{eqnarray}\begin{array}{*{35}{rl}} \bot {{R}_{A0C0}}\equiv \bot \left({{R}_{ABCD}}{{n}^{B}}{{n}^{D}}\right)= & \bot {{R}_{AC}}+{{\mathcal{R}}_{AC}}+K{{K}_{AC}}-{{K}_{AE}}{{K}^{E}}{{}_{C}} \\ = & {{\mathcal{R}}_{AC}}+K{{K}_{AC}}-{{K}_{AE}}{{K}^{E}}{{}_{C}}, \end{array}\end{eqnarray} \tag{ 4.6 }$$

where in the last line we used the fact that in vacuum R_AC and, hence, its projection vanishes (note, however, that in general ${{\mathcal{R}}_{AC}}\ne 0$ even in vacuum). Furthermore ${{D}_{C}}{{K}_{AD}}={{\partial}_{C}}{{K}_{AD}}- \Gamma _{CA}^{B}{{K}_{BD}}- \Gamma _{CD}^{B}{{K}_{AB}}$ is the covariant derivative of the extrinsic curvature defined on the spatial hypersurface, with Christoffel symbols calculated from the induced metric ${{\gamma}_{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 equations (4.4)–(4.6) the spacetime indices $A,\,B,\,\ldots$ on the left and right-hand side by spatial indices $I,\,J,\,\ldots$ 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) and (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.

4.1.2. The Riemann tensor in SO(D  −  3) symmetry.

The expressions given in the previous section 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 SO(D  −  d) isometry for $1\leqslant d\leqslant D-3$.

For this purpose we recall the Cartesian coordinate system ${{X}^{I}}=\left({{x}^{{\hat{i}}}},\,\,z,\,{{w}^{a}}\right)$ of equation (2.5), adapted to a spacetime that is symmetric under rotations in any plane spanned by two of the $\left(z,\,{{w}^{a}}\right)$. We discuss in turn how the terms appearing on the right-hand sides of equations (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

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{IJKL}}= & \frac{1}{2}\left({{\partial}_{L}}{{\partial}_{I}}{{\gamma}_{JK}}+{{\partial}_{K}}{{\partial}_{J}}{{\gamma}_{IL}}-{{\partial}_{K}}{{\partial}_{I}}{{\gamma}_{JL}}-{{\partial}_{L}}{{\partial}_{J}}{{\gamma}_{IK}}\right) \\ {} & -{{\gamma}_{MN}} \Gamma _{IK}^{N} \Gamma _{JL}^{M}+{{\gamma}_{MN}} \Gamma _{IL}^{N} \Gamma _{JK}^{M}. \end{array}\end{eqnarray} \tag{ 4.7 }$$

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 section 2.2 and appendix A of [83]. 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

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{ijkl}}= & \frac{1}{2}\left({{\partial}_{l}}{{\partial}_{i}}{{\gamma}_{jk}}+{{\partial}_{k}}{{\partial}_{j}}{{\gamma}_{il}}-{{\partial}_{k}}{{\partial}_{i}}{{\gamma}_{jl}}-{{\partial}_{l}}{{\partial}_{j}}{{\gamma}_{ik}}\right)-{{\gamma}_{mn}} \Gamma _{ik}^{n} \Gamma _{jl}^{m}+{{\gamma}_{mn}} \Gamma _{il}^{n} \Gamma _{jk}^{m}, \end{array}\end{eqnarray} \tag{ 4.8 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{ajkl}}= & 0, \end{array}\end{eqnarray} \tag{ 4.9 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{iajb}}= & {{\delta}_{ab}}{{\mathcal{R}}_{iwjw}}, \end{array}\end{eqnarray} \tag{ 4.10 }$$

$$\begin{eqnarray} \begin{array}{*{35}{l}} {{\mathcal{R}}_{iwjw}}\equiv & \frac{{{\partial}_{(i}}{{\gamma}_{j)z}}-{{\delta}_{z(\,j}}{{\partial}_{i)}}{{\gamma}_{ww}}}{z}-{{\delta}_{z(i}}\frac{{{\gamma}_{j)z}}-{{\delta}_{j)z}}{{\gamma}_{ww}}}{{{z}^{2}}}-\frac{1}{2}{{\partial}_{j}}{{\partial}_{i}}{{\gamma}_{ww}}-{{\gamma}_{mn}} \Gamma _{ij}^{n} \Gamma _{ww}^{m} \\ {} & -\frac{1}{2}\frac{{{\partial}_{z}}{{\gamma}_{ij}}}{z}+\frac{{{\delta}_{z(i}}{{\gamma}_{j)z}}-{{\delta}_{iz}}{{\delta}_{jz}}{{\gamma}_{ww}}}{{{z}^{2}}}+\frac{1}{4}{{\gamma}^{ww}}{{\partial}_{i}}{{\gamma}_{ww}}{{\partial}_{j}}{{\gamma}_{ww}}, \end{array}\end{eqnarray} \tag{ 4.11 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Gamma _{ww}^{m}\equiv & -\frac{1}{2}{{\gamma}^{ml}}{{\partial}_{l}}{{\gamma}_{ww}}+\frac{\delta _{z}^{m}-{{\gamma}^{mz}}{{\gamma}_{ww}}}{z}, \end{array}\end{eqnarray} \tag{ 4.12 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{abcl}}= & 0, \end{array}\end{eqnarray} \tag{ 4.13 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{abcd}}= & \left({{\delta}_{ac}}{{\delta}_{bd}}-{{\delta}_{bc}}{{\delta}_{ad}}\right){{\mathcal{R}}_{wuwu}}, \end{array}\end{eqnarray} \tag{ 4.14 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{wuwu}}\equiv & -\frac{1}{4}{{\gamma}^{mn}}{{\partial}_{m}}{{\gamma}_{ww}}{{\partial}_{n}}{{\gamma}_{ww}}-{{\gamma}_{ww}}\frac{{{\gamma}^{zm}}}{z}{{\partial}_{m}}{{\gamma}_{ww}}+\frac{{{\gamma}_{ww}}-{{\gamma}^{zz}}\gamma _{ww}^{2}}{{{z}^{2}}}. \end{array}\end{eqnarray} \tag{ 4.15 }$$

For the right-hand side of equation (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

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{ij}}= & {{\gamma}^{mn}}{{\mathcal{R}}_{minj}}+(D-d-1){{\gamma}^{ww}}{{\mathcal{R}}_{iwjw}}, \end{array}\end{eqnarray} \tag{ 4.16 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{ab}}= & {{\delta}_{ab}}{{\mathcal{R}}_{ww}}, \end{array}\end{eqnarray} \tag{ 4.17 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\mathcal{R}}_{ww}}\equiv & {{\gamma}^{mn}}{{\mathcal{R}}_{mwnw}}+(D-d-2){{\gamma}^{ww}}{{\mathcal{R}}_{wuwu}}. \end{array}\end{eqnarray} \tag{ 4.18 }$$

Note that the last expression, ${{\gamma}^{ww}}{{\mathcal{R}}_{wuwu}}$, does not involve a summation over w, but merely stands for the product of ${{\gamma}^{ww}}$ with the expression (4.15).

The components of the extrinsic curvature are given by equation (4.1). Its derivative is directly obtained from the expressions (A.1)–(A.12) in appendix A of [83] and can be written as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{D}_{i}}{{K}_{jk}}= & {{\partial}_{i}}{{K}_{jk}}- \Gamma _{ij}^{l}{{K}_{kl}}- \Gamma _{ik}^{l}{{K}_{lj}}, \end{array}\end{eqnarray} \tag{ 4.19 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{D}_{i}}{{K}_{ab}}= & {{\delta}_{ab}}\left({{\partial}_{i}}{{K}_{ww}}-{{K}_{ww}}{{\gamma}^{ww}}{{\partial}_{i}}{{\gamma}_{ww}}\right), \end{array}\end{eqnarray} \tag{ 4.20 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{D}_{a}}{{K}_{bj}}= & {{\delta}_{ab}}\left(\frac{{{K}_{jz}}-{{\delta}_{jz}}{{K}_{ww}}}{z}-\frac{1}{2}{{K}_{ww}}{{\gamma}^{ww}}{{\partial}_{j}}{{\gamma}_{ww}}-{{K}_{ij}} \Gamma _{ww}^{i}\right). \end{array}\end{eqnarray} \tag{ 4.21 }$$

Next, we plug the expressions assembled in equations (4.7)–(4.21) into the Gauss–Codazzi equations (4.4)–(4.6) where, we recall, early Latin indices $A,\,B,\,\ldots$ are now replaced by $I,\,J,\,\ldots$ following our switch to adapted coordinates. Splitting the index range I into $(i,\,a)$ for components inside and outside the computational domain, and recalling that an index 0 denotes contraction with $\boldsymbol{n}$, we can write the resulting components of the spacetime Riemann tensor as

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{ijkl}} & ={{\mathcal{R}}_{ijkl}}+{{K}_{ik}}{{K}_{jl}}-{{K}_{il}}{{K}_{jk}}, \end{array}\end{eqnarray} \tag{ 4.22 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{ibkd}} & ={{\delta}_{bd}}{{R}_{iwkw}}, \end{array}\end{eqnarray} \tag{ 4.23 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{iwkw}} & \equiv {{\mathcal{R}}_{iwkw}}+{{K}_{ik}}{{K}_{ww}}, \end{array}\end{eqnarray} \tag{ 4.24 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{abcd}} & =\left({{\delta}_{ac}}{{\delta}_{bd}}-{{\delta}_{bc}}{{\delta}_{ad}}\right)\left({{\mathcal{R}}_{wuwu}}+K_{ww}^{2}\right), \end{array}\end{eqnarray} \tag{ 4.25 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{ajkl}} & ={{R}_{abcl}}=0, \end{array}\end{eqnarray} \tag{ 4.26 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{i0kl}} & ={{D}_{l}}{{K}_{ik}}-{{D}_{k}}{{K}_{il}}, \end{array}\end{eqnarray} \tag{ 4.27 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{a0ck}} & ={{\delta}_{ac}}{{R}_{w0wk}}, \end{array}\end{eqnarray} \tag{ 4.28 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{w0wk}} & \equiv {{\partial}_{k}}{{K}_{ww}}-\frac{1}{2}{{\gamma}^{ww}}{{K}_{ww}}{{\partial}_{k}}{{\gamma}_{ww}}-\frac{{{K}_{kz}}-{{\delta}_{kz}}{{K}_{ww}}}{z}+ \Gamma _{ww}^{m}{{K}_{mk}}, \end{array}\end{eqnarray} \tag{ 4.29 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{a0cd}} & ={{R}_{i0kd}}={{R}_{a0kl}}=0, \end{array}\end{eqnarray} \tag{ 4.30 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{i0j0}} & ={{\mathcal{R}}_{ij}}+K{{K}_{ij}}-{{K}_{im}}K_{j}^{m}, \end{array}\end{eqnarray} \tag{ 4.31 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} K & ={{\gamma}^{mn}}{{K}_{mn}}+(D-d-1){{\gamma}^{ww}}{{K}_{ww}}, \end{array}\end{eqnarray} \tag{ 4.32 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{a0b0}} & ={{\delta}_{ab}}{{R}_{w0w0}}, \end{array}\end{eqnarray} \tag{ 4.33 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{w0w0}} & \equiv {{\mathcal{R}}_{ww}}+\left(K-{{\gamma}^{ww}}{{K}_{ww}}\right){{K}_{ww}}, \end{array}\end{eqnarray} \tag{ 4.34 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{R}_{a0i0}} & =0. \end{array}\end{eqnarray} \tag{ 4.35 }$$

With these expressions, we are able to calculate all components of the spacetime Riemann tensor directly from the ADM variables ${{\gamma}_{ij}}$, ${{\gamma}_{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.

4.2. The null frame

The null frame we need for the projections of the Weyl tensor consists of D unit vectors as given in equation (3.3): (i) the ingoing null vector k^A, (ii) the outgoing null vector l^A which, however, does not explicitly appear in the scalars (3.4) for the outgoing gravitational radiation, and (iii) the (D  −  2) vectors $m_{(\alpha )}^{A}$ pointing in the angular directions on the sphere.

We begin this construction with the D  −  2 unit basis vectors on the D  −  2 sphere, $m_{(\alpha )}^{A}$, and recall for this purpose equation (2.8) that relates our spherical coordinates $\left(r,{{\phi}^{\alpha}}\right)$ to the Cartesian $\left({{x}^{{\hat{i}}}},\,z,\,{{w}^{a}}\right)$. The set of spatial vectors, although not yet in orthonormalised form, then consists of a radial vector denoted by ${{\tilde{m}}_{(1)}}$ and D  −  2 angular vectors ${{\tilde{m}}_{(\alpha )}}$ whose components in Cartesian coordinates ${{X}^{I}}=\left({{x}^{{\hat{i}}}},\,z,\,{{w}^{a}}\right)$ on the computational domain w^a  =  0 are obtained through the chain rule

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\tilde{m}}}_{(1)}} & =\frac{\partial}{\partial r}=\frac{\partial {{X}^{I}}}{\partial r}\frac{\partial}{\partial {{X}^{I}}}~~\Rightarrow ~~\tilde{m}_{(1)}^{I}=\frac{1}{r}\left({{x}^{1}},\,\ldots,\,{{x}^{d-1}},\,z,\,0,\,\ldots,\,0\right), \end{array}\end{eqnarray} \tag{ 4.36 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\tilde{{m}}\,}_{(\alpha )}} & =\frac{\partial}{\partial {{\phi}^{\alpha}}}=\frac{\partial {{X}^{I}}}{\partial {{\phi}^{\alpha}}}\frac{\partial}{\partial {{X}^{I}}}. \end{array}\end{eqnarray} \tag{ 4.37 }$$

We can ignore time components here, because our coordinates are adapted to the space-time split, so that all spatial vectors have vanishing time components and this feature is preserved under the eventual Gram–Schmidt orthonormalisation. Plugging equation (2.8) into (4.37), we obtain for ${{\tilde{m}}_{(\alpha )}}$ (after rescaling by $r\times \sin {{\phi}^{2}}\times \ldots \times \sin {{\phi}^{\alpha}}$)

$$\begin{eqnarray} &&\underset{=\tilde{m}_{(2)}^{I}}{{\underbrace{{\left(\begin{array}{c} -{\sum}_{s=2}^{D-1}{{\left({{w}^{s}}\right)}^{2}} \\ {{w}^{1}}{{w}^{2}} \\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ {{w}^{1}}{{w}^{D-1}} \end{array}\right)}}\,}}\,,\ldots,\underset{=\tilde{m}_{(\alpha )}^{I}}{{\underbrace{{\left(\begin{array}{c} \left.\begin{array}{c} 0~ \\ \vdots ~ \\ 0~ \end{array}\right\}(\alpha -2)\times \\ -{\sum}_{s=\alpha}^{D-1}{{\left({{w}^{s}}\right)}^{2}} \\ {{w}^{\alpha -1}}{{w}^{\alpha}} \\ \vdots \\ {{w}^{\alpha -1}}{{w}^{D-2}} \\ {{w}^{\alpha -1}}{{w}^{D-1}} \end{array}\right)}}\,}}\,,\ldots,\underset{=\tilde{m}_{(D-2)}^{I}}{{\underbrace{{\left(\begin{array}{c} 0 \\ \vdots \\ \vdots \\ \vdots \\ 0 \\ -{{\left({{w}^{D-2}}\right)}^{2}}-{{\left({{w}^{D-1}}\right)}^{2}} \\ {{w}^{D-3}}{{w}^{D-2}} \\ {{w}^{D-3}}{{w}^{D-1}} \end{array}\right)}}\,}}\,,~\underset{=\tilde{m}_{(D-1)}^{I}}{{\underbrace{{\left(\begin{array}{c} 0 \\ \vdots \\ \vdots \\ \vdots \\ 0 \\ 0 \\ -{{\left({{w}^{D-1}}\right)}^{2}} \\ {{w}^{D-2}}{{w}^{D-1}} \end{array}\right)}}\,}}\,.\end{eqnarray} \tag{ 4.38 }$$

In D  =  4 dimensions, these vectors, together with ${{\tilde{m}}_{(1)}}$ of equation (4.36) would form the starting point for Gram–Schmidt orthonormalisation; see e.g. appendix C in [57]. In $D\geqslant 5$ dimensional spacetimes with SO(D  −  d) symmetry, however, we face an additional difficulty: on the computational domain w^a  =  0, all components of the vectors ${{\tilde{m}}_{(d+1)}},\,\ldots,\,{{\tilde{m}}_{(D-1)}}$ vanish and their normalisation would result in divisions of zero by zero. This difficulty is overcome by rewriting the Cartesian components of the vectors in terms of spherical coordinates and then exploiting the freedom we have in suitably orienting the frame. The details of this procedure are given in appendix B where we derive a manifestly regular set of spatial vectors given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(1)}^{A} & =\left(0~|~{{x}^{1}},\,\ldots,\,{{x}^{d}}~|~0,\ldots,\,0\right), \end{array}\end{eqnarray} \tag{ 4.39 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(2)}^{A} & =\left(0~|~-\rho _{2}^{2},~{{x}^{1}}{{x}^{2}},~{{x}^{1}}{{x}^{3}},~\ldots,~{{x}^{1}}{{x}^{d}}~|~0,~\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.40 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \cdots & \cdots \\ \tilde{m}_{\left(\hat{\alpha}\right)}^{A} & =\left(0,~|~\underset{\left(\hat{\alpha}-2\right)\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,,~-\rho _{{\hat{\alpha}}}^{2},~{{x}^{\hat{\alpha}-1}}{{x}^{{\hat{\alpha}}}},~\ldots,~{{x}^{\hat{\alpha}-1}}{{x}^{d}}~|~0,~\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.41 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \cdots & \cdots \\ \tilde{m}_{(d)}^{A} & =\left(0,~|~\underset{(d-2)\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,,~-{{x}^{d}},~{{x}^{d-1}}~|~0,~\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.42 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(d+1)}^{A} & =\left(0~|~\underset{d\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,~|~1,~0,~\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.43 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \cdots & \cdots \\ \tilde{m}_{(D-1)}^{A} & =\left(0~|~\underset{d\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,~|~0,~\ldots,~0,~1\right), \end{array}\end{eqnarray} \tag{ 4.44 }$$

where ${{\rho}_{I}}={\sum}_{s=I}^{D-1}{{\left({{w}^{s}}\right)}^{2}}$, we have restored, for completeness, the time component and the vertical bars highlight the three component sectors: time, spatial on-domain, and spatial off-domain. Equations (4.43) and (4.44) can, of course, be conveniently written in short-hand notation as $\tilde{m}_{(a)}^{A}={{\delta}^{A}}{{}_{a}}$. For the special case d  =  3, the vectors are given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(1)}^{A} & =\left(0~|~x,~y,~z~|~0,\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.45 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(2)}^{A} & =\left(0~|~-{{y}^{2}}-{{z}^{2}},~xy,~xz~|~0,\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.46 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(3)}^{A} & =\left(0~|~0,~-z,~y~|~0,~\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.47 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(4)}^{A} & =\left(0~|~0,~0,~0~|~1,~0,~\ldots,~0\right), \end{array}\end{eqnarray} \tag{ 4.48 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \cdots & \cdots \end{array}\end{eqnarray} \tag{ 4.49 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \tilde{m}_{(D-1)}^{A} & =\left(0~|~0,~0,~0~|~0,~\ldots,~0,~1\right), \end{array}\end{eqnarray} \tag{ 4.50 }$$

The next step is to orthonormalise these vectors. Clearly the vectors $m_{(a)}^{A}$ with components in the w^a dimensions are normalised by:

$$\begin{eqnarray}&&m_{(a)}^{A}=\frac{1}{\sqrt{{{\gamma}_{ww}}}}{{\delta}^{A}}{{}_{a}}\end{eqnarray} \tag{ 4.51 }$$

For the remaining d vectors given by equations (4.39)–(4.42) or, for d  =  3, the spatial triad consisting of the three vectors (4.45)–(4.47), we use standard Gram–Schmidt orthonormalisation. 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 $\partial /\partial u-\frac{1}{2}\partial /\partial \mathfrak{r}$ asymptotically, we transform out of Bondi coordinates, sending $\left(u,\mathfrak{r}\right)\to (t,r)$ and furthermore use the freedom of rescaling this null vector by applying a constant factor of⁵ $\sqrt{2}$

$$\begin{eqnarray}&&{{k}^{A}}=\frac{1}{\sqrt{2}}\left({{n}^{A}}-m_{(1)}^{A}\right)\end{eqnarray} \tag{ 4.52 }$$

Expressing the timelike unit normal field n^A in terms of our gauge variables $\alpha,{{\beta}^{I}}$ we find

$$\begin{eqnarray}&&{{k}^{A}}=\frac{1}{\sqrt{2}}\left(\frac{1}{\alpha},-\frac{{{\beta}^{I}}}{\alpha}-m_{(1)}^{I}\right),\end{eqnarray} \tag{ 4.53 }$$

where ${{\beta}^{I}}=\left({{\beta}^{i}},\,0,\,\ldots,\,0\right),~m_{(1)}^{I}=\left(m_{(1)}^{i},\,0,\ldots,\,0\right)$. This result provides the ingoing null vector for any choice of d and is the version implemented in the code.

4.3. The projections of the Weyl tensor

Finally, we calculate the projections of the Weyl tensor that encode the outgoing gravitational radiation

$$\begin{eqnarray} &&\Omega _{(\alpha )(\beta )}^{\prime}={{R}_{ABCD}}{{k}^{A}}m_{(\alpha )}^{B}{{k}^{C}}m_{(\beta )}^{D},\end{eqnarray} \tag{ 4.54 }$$

(see equation (3.4)) where k^A is given by equation (4.53) and the normal frame vectors ${{m}_{(2)}},\,\ldots,\,{{m}_{(D-1)}}$ are those obtained from Gram–Schmidt orthonormalising the right-hand sides of equations (4.39)–(4.44).

We first note that $\Omega _{(\alpha )(\beta )}^{\prime}$ is symmetric in $\alpha \leftrightarrow \beta$, so contractions solely with ${{m}_{(2)}},\,\ldots,\,{{m}_{(d)}}$ will result in d(d  −  1)/2 components $\Omega _{\left(\hat{\alpha}\right)\left(\hat{\beta}\right)}^{\prime}$. For the special case d  =  3, we obtain the three components $\Omega _{(2)(2)}^{\prime},\, \Omega _{(2)(3)}^{\prime},\, \Omega _{(3)(3)}^{\prime}$. The null vector k has vanishing w components and from equations (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,\,\ldots$ are zero. The only non-vanishing terms involving the Riemann tensor with off-domain indices $a,\,b,\,\ldots$, therefore, have either four such indices or two and contain a Kronecker delta ${{\delta}_{ab}}$; see equations (4.23), (4.25), (4.28) and (4.33). As a consequence, the mixed components $\Omega _{\left(\hat{\alpha}\right)(a)}^{\prime}=0$ and the purely off-domain components $\Omega _{(a)(b)}^{\prime}\propto {{\delta}_{ab}}$. The list of all non-vanishing components $\Omega _{(\alpha )(\beta )}^{\prime}$ is then given by

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \Omega _{\left(\hat{\alpha}\right)\left(\hat{\beta}\right)}^{\prime}=\frac{1}{4}\left[{{R}_{0k0l}}m_{\left(\hat{\alpha}\right)}^{k}m_{\left(\hat{\beta}\right)}^{l}-{{R}_{mk0l}}m_{(1)}^{m}m_{\left(\hat{\alpha}\right)}^{k}m_{\left(\hat{\beta}\right)}^{l}-{{R}_{0kml}}m_{\left(\hat{\alpha}\right)}^{k}m_{(1)}^{m}m_{\left(\hat{\beta}\right)}^{l} \right. \\ {} \\ \left. \qquad \quad \ +\,\,{{R}_{mknl}}m_{(1)}^{m}m_{\left(\hat{\alpha}\right)}^{k}m_{(1)}^{n}m_{\left(\hat{\beta}\right)}^{l}\right], \end{array}\end{eqnarray} \tag{ 4.55 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Omega _{(a)(b)}^{\prime} & ={{\delta}_{ab}}\, \Omega _{(w)(w)}^{\prime}, \end{array}\end{eqnarray} \tag{ 4.56 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \Omega _{(w)(w)}^{\prime} & =\frac{1}{4{{\gamma}_{ww}}}\left[{{R}_{w0w0}}-{{R}_{w0wk}}m_{(1)}^{k}-{{R}_{w0wl}}m_{(1)}^{l}+{{R}_{wkwl}}m_{(1)}^{k}m_{(1)}^{l}\right], \end{array}\end{eqnarray} \tag{ 4.57 }$$

where $\hat{\alpha},\hat{\beta}=2,\,\ldots,\,d$ and all components of the Riemann tensor on the right-hand sides are listed in the set of equations (4.22)–(4.35). In particular, the components R_w0w0, R_w0wk and R_wkwl, which contain indices in the off-domain directions, are obtained from equations (4.34), (4.29) and (4.24), respectively and thus derived directly from quantities computed in the simulation (the ${{\gamma}_{ij}}$, K_ij, ${{\gamma}_{ww}}$ and K_ww that appear on the right-hand sides of these equations or enter in the calculation of the spatial Riemann tensor). It should be noted here that $\Omega _{(\alpha )(\beta )}^{\prime}$ is trace free, and so $\Omega _{(w)(w)}^{\prime}$ can be calculated from the diagonal terms $\Omega _{(2)(2)}^{\prime},\ldots \Omega _{(d)(d)}^{\prime}$. In a numerical simulation, the components of $\Omega _{(\alpha )(\beta )}^{\prime}$ are calculated as functions of time and then can be integrated according to equation (3.6) to extract the amount of energy radiated in gravitational waves.

4.4. SO(2) symmetry

In the remainder of this work, we will implement the specific version of the wave extraction for 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 [69] who calculate the GW energy using the Kodama–Ishibashi formalism, which enables us to compare our findings with their values.

5.1. Code infrastructure and numerical set-up

We perform evolutions using the LEAN code [57, 97] which is based on CACTUS [98, 99] and uses CARPET [100, 101] 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 [83] 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 '$\Gamma$ driver' conditions for the lapse function and shift vector (see e.g. [102]) according to

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{t}}\alpha = & {{\beta}^{m}}{{\partial}_{m}}\alpha -2\alpha K, \end{array}\end{eqnarray} \tag{ 5.1 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\partial}_{t}}{{\beta}^{i}}= & {{\beta}^{m}}{{\partial}_{m}}{{\beta}^{i}}+\frac{1}{4}{{{\rm{\tilde \Gamma }}}^{i}}-\frac{1}{{{2}^{1/3}}{{R}_{h}}}{{\beta}^{i}}, \end{array}\end{eqnarray} \tag{ 5.2 }$$

with initial values $\alpha =1$, ${{\beta}^{i}}=0$.

The BH initial data is calculated using the higher dimensional generalisation of Brill–Lindquist data [103, 104] given in terms of the ADM variables by

$$\begin{eqnarray}&&{{K}_{IJ}}=0,~~{{\gamma}_{IJ}}={{\psi}^{4/(D-3)}}{{\delta}_{IJ}},~~\psi =1+\underset{\mathcal{N}}{\sum}\,\frac{{{\mu}_{\mathcal{N}}}}{4{{\left[{{\sum}_{K=1}^{D-1}}\,{{\left({{X}^{K}}-X_{\mathcal{N}}^{K}\right)}^{2}}\right]}^{(D-3)/2}}},\end{eqnarray} \tag{ 5.3 }$$

where the summation over $\mathcal{N}$ and K extends over the multiple BHs and spatial coordinates, respectively, and $X_{\mathcal{N}}^{K}$ denotes the position of the $\mathcal{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_{\mathcal{N}}^{1}=\pm {{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

$$\begin{eqnarray}&&\mu =\frac{16\pi M}{(D-2){{\mathcal{A}}_{D-2}}},~~~\mu =R_{\text{h}}^{D-3},~~~{{\mathcal{A}}_{D-2}}=\frac{2{{\pi}^{(D-1)/2}}}{ \Gamma \left(\frac{D-1}{2}\right)},\end{eqnarray} \tag{ 5.4 }$$

where ${{\mathcal{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 characterise 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 $\Omega _{(\alpha )(\beta )}^{\prime}$ 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 equation (3.6). The $\Omega _{(\alpha )(\beta )}^{\prime}$ are scalars and so in our angular coordinate system do not depend on ${{\phi}^{4}},\ldots,{{\phi}^{D-1}}$, so the integral over the sphere in (3.6) can be simplified:

$$\begin{eqnarray}&&\overset{\centerdot}{{M}}\,(u)=-\underset{r\to \infty}{{\lim}}\,\frac{{{r}^{D-2}}}{8\pi}{{\mathcal{A}}_{D-4}}{\int}_{0}^{\pi}{\int}_{0}^{\pi}I\left[{{ \Omega }^{\prime}}^{2}\right]~{{\sin}^{D-3}}\left({{\phi}^{2}}\right)~{{\sin}^{D-4}}\left({{\phi}^{3}}\right)~\text{d}{{\phi}^{3}}\text{d}{{\phi}^{2}},\end{eqnarray} \tag{ 5.5 }$$

where $I\left[{{ \Omega }^{\prime 2}}\right]\equiv {{\left({\int}_{-\infty}^{u} \Omega _{(\alpha )(\beta )}^{\prime}\text{d}\tilde{u}\right)}^{2}}$. A final integration over time of the variable $\overset{\centerdot}{{M}}\,$ then gives the total radiated energy.

5.2. 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=\pm {{x}_{0}}=\pm 4.0~{{R}_{h}}$ using a computational grid of size L  =  181R_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 equation (5.4) with $\mu \equiv {{\mu}_{1}}+{{\mu}_{2}}$, the mass parameters of the initial BHs. The radiated energy as a function of time is shown in the upper panel of figure 2. The radiation is almost exclusively concentrated within a window of $\Delta t\approx 20~{{R}_{h}}$ around merger. During the infall and the post-merger period, in contrast, E_rad remains nearly constant. In comparison with collisions in D  =  4 dimensions, we find the burst of spurious (colloquially referred to as 'junk') radiation significantly weaker, presumably because the Brill–Lindquist data in higher D more closely represent two black holes in isolation due to the higher fall-off rate of the gravitational interaction. 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}}+\mathcal{O}\left({{h}^{n}}\right)$, where n denotes the order of convergence. By evaluating the quotient

$$\begin{eqnarray}&&{{Q}_{n}}=\frac{{{f}_{{{h}_{1}}}}-{{f}_{{{h}_{2}}}}}{{{f}_{{{h}_{2}}}}-{{f}_{{{h}_{3}}}}}=\frac{{{\left({{h}_{1}}/{{h}_{2}}\right)}^{n}}-1}{1-{{\left({{h}_{3}}/{{h}_{2}}\right)}^{n}}},\end{eqnarray} \tag{ 5.6 }$$

we can then plot the two differences ${{f}_{{{h}_{1}}}}-{{f}_{{{h}_{2}}}}$ and ${{f}_{{{h}_{2}}}}-{{f}_{{{h}_{3}}}}$ and test whether their ratio is consistent with a given value n. The results for our study are shown in the lower panel of figure 2 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 a total radiated energy ${{E}_{\text{rad}}}=8.19\times {{10}^{-4}}~{{M}_{\text{ADM}}}$ with a discretisation error of  ∼0.4%, but note that the error in the cumulative energy peaks at a larger value of a few $\%$ during the sharp increase of ${{E}_{\text{rad}}}(t)$ marking the merger phase.

Zoom In Zoom Out Reset image size

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 40R_h to 110R_h and extrapolating these values assuming a functional dependency

$$\begin{eqnarray}&&{{E}_{\text{rad}}}(r)={{E}_{\text{rad}}}(\infty )+\frac{a}{r}+\mathcal{O}\left(\frac{1}{{{r}^{2}}}\right),\end{eqnarray} \tag{ 5.7 }$$

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}_{\text{ex}}}=110~{{\text{R}}_{\text{h}}}$ and 0.4% at ${{R}_{\text{ex}}}=60~{{R}_{h}}$.

An independent check of our results is available in comparing the radiated energy with the predictions of the perturbative extraction method [66] based on the Kodama–Ishibashi formalism. For this purpose, we have calculated using ${{h}_{3}}={{R}_{h}}/76.2$ the gravitational-wave energy radiated in the quadrupole mode as predicted by the Kodama–Ishibashi formalism. Contributions from higher-order multipoles are negligible for this comparison; for odd l they vanish completely by symmetry and for even l up to l  =  8 they are well below the numerical uncertainty budget. This quadrupole energy is compared with the result obtained from the Weyl tensor in figure 3. The difference for the total radiated energy is about 0.3%, though a larger temporary discrepancy for E_rad as a function of time is encountered during the steep increase at merger, up to a few $\%$. This discrepancy is within the error budget of the two extraction methods.

Zoom In Zoom Out Reset image size

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}}=\pm 7.8~{{R}_{h}}$ and ${{x}_{0}}=\pm 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 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

$$\begin{eqnarray}&&{{E}_{\text{rad}}}=(8.19\pm 0.05)\times {{10}^{-4}}~{{M}_{\text{ADM}}},\end{eqnarray} \tag{ 5.8 }$$

in excellent agreement with the value $(8.1\pm 0.4)\times {{10}^{-4}}$ reported in the independent study by [69] using dimensional reduction by isometry and the Kodama–Ishibashi formalism.

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 we have formulated the Weyl tensor based wave extraction technique of Godazgar and Reall [1]—a higher dimensional generalisation of the Newman–Penrose scalars—in a form suitable for numerical simulations of D  >  4 dimensional spacetimes with SO(D  −  d), $1\leqslant d\leqslant D-2$, symmetry employing the modified Cartoon method. The only prerequisite for implementing our formalism is the availability of the ADM variables on each spatial hypersurface of the effective computational domain. These are constructed straightforwardly from all commonly used numerical evolution systems such as BSSN, generalised harmonic or conformal Z4.

The recipe for extracting the GW signal then consists of the following steps.

• (1)  
  Computation of the on and off-domain components of the spatial Riemann tensor (which equals the Weyl tensor in the vacuum extraction region) and the derivative of the extrinsic curvature according to equations (4.8)–(4.21).
• (2)  
  Reconstruction of the components of the spacetime Riemann tensor as well as its contractions with the unit timelike normal from the quantities of the previous step according to equations (4.22)–(4.35).
• (3)  
  Construction of the null-frame vectors through Gram–Schmidt orthonormalising the expressions of equations (4.39)–(4.44) and then using (4.53) for the ingoing null vector.
• (4)  
  Calculation of the projections $\Omega _{(\alpha )(\beta )}^{\prime}$ of the Weyl tensor onto the null frame vectors using equations (4.55)–(4.57).
• (5)  
  Calculation of the energy flux in GWs through equation (3.6) and integration in time of the flux to obtain the total radiated energy.

The most common case of modelling higher dimensional spacetimes with rotational symmetries is the case of d  =  3 effective spatial dimensions which allows for straightforward generalisation 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 head-on collisions of BHs starting from rest in D  =  6 using d  =  3. We find these collisions to radiate a fraction $(8.19\pm 0.05)\times {{10}^{-4}}$ of the ADM mass in GWs, in excellent agreement with a previous study [69] 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 $\Omega _{(\alpha )(\beta )}^{\prime}$ 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\gg 4$, and are therefore beyond the scope of the present study. However, the foundation for analysing in detail the GW energy 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.

We thank Pau Figueras, Mahdi Godazgar, Markus Kunesch, Harvey Reall, Saran Tunyasuvunakool and Helvi Witek for highly fruitful discussions of this topic. This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie SkŁodowska-Curie grant agreement No 690904, from H2020-ERC-2014-CoG Grant No. 'MaGRaTh' 646597, from STFC Consolidator Grant No. ST/L000636/1, the SDSC Comet, PSC-Bridges and TACC Stampede clusters through NSF-XSEDE Award Nos. PHY-090003, the Cambridge High Performance Computing Service Supercomputer Darwin using Strategic Research Infrastructure Funding from the HEFCE and the STFC, and DiRAC's Cosmos Shared Memory system through BIS Grant No. ST/J005673/1 and STFC Grant Nos. ST/H008586/1, ST/K00333X/1. WGC is supported by a STFC studentship.

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 [83], so that in the following we can focus exclusively on the additional terms appearing for $1\leqslant d\leqslant 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 [83]. 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 \tilde{\gamma}=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:

A.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 $\boldsymbol{V}$, we obtain

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{{{V}^{z}}}{z}= & \frac{{{a}_{1}}z+{{a}_{3}}{{z}^{3}}+\ldots}{z}={{a}_{1}}+{{a}_{3}}{{z}^{2}}+\ldots \overset{\ast}{{=}}\,{{a}_{1}}\overset{\ast}{{=}}\,{{\partial}_{z}}{{V}^{z}}, \end{array}\end{eqnarray} \tag{ A.1 }$$

where we have introduced the symbol $\overset{\ast}{{=}}\,$ to denote equality in the limit $z\to 0$.

A.2. 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 $\left({{x}^{{\hat{i}}}},z,{{w}^{d+1}},\ldots,{{w}^{a}}\ldots,{{w}^{D-1}}\right)$ to $\left({{x}^{{\hat{i}}}},\rho,{{w}^{d+1}},\ldots,{{w}^{a-1}},\varphi,{{w}^{a+1}},\ldots,{{w}^{D-1}}\right)$. As no other ${{w}^{b}},~b\ne a$, coordinates will enter into this discussion we shall refer to w^a as w. In these coordinates we have that

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\gamma}_{\rho \rho}}= & \frac{{{z}^{2}}}{{{\rho}^{2}}}{{\gamma}_{zz}}+2\frac{zw}{{{\rho}^{2}}}{{\gamma}_{zw}}+\frac{{{w}^{2}}}{{{\rho}^{2}}}{{\gamma}_{ww}}, \end{array}\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{\gamma}_{\varphi \varphi}}= & {{w}^{2}}{{\gamma}_{zz}}-2wz{{\gamma}_{zw}}+{{z}^{2}}{{\gamma}_{ww}}, \end{array}\end{eqnarray} \tag{ A.3 }$$

and the line element for vanishing $\text{d}{{x}^{{\hat{i}}}}=0$ and $\text{d}{{w}^{b}}=0$, $b\ne a$, is given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} \text{d}{{s}^{2}}= & {{\gamma}_{\rho \rho}}\text{d}{{\rho}^{2}}+{{\gamma}_{\varphi \varphi}}\text{d}{{\varphi}^{2}}. \end{array}\end{eqnarray} \tag{ A.4 }$$

Requiring the circumference to be the radius times $2\pi$, we have that ${{\gamma}_{\varphi \varphi}}={{\rho}^{2}}{{\gamma}_{\rho \rho}}$. Substituting the above expressions and taking the limit $z\to 0$, we obtain

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\gamma}_{zz}}-{{\gamma}_{ww}}\overset{\ast}{{=}}\,\mathcal{O}\left({{z}^{2}}\right). \end{array}\end{eqnarray} \tag{ A.5 }$$

Taking the time derivative of this relation and using the definition of the extrinsic curvature, we find that likewise

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{K}_{zz}}-{{K}_{ww}}\overset{\ast}{{=}}\,\mathcal{O}\left({{z}^{2}}\right). \end{array}\end{eqnarray} \tag{ A.6 }$$

A.3. Inverse metric

Various terms that we need to address contain factors of the inverse metric ${{\gamma}^{IJ}}$. In the practical regularisation procedure, these terms are conveniently handled by expressing ${{\gamma}^{IJ}}$ in terms of the downstairs metric components ${{\gamma}_{ij}}$ and ${{\gamma}_{ww}}$ which are the fields we evolve numerically. We know the metric takes the following form:

$$\begin{eqnarray}{{\gamma}_{IJ}}=\left(\begin{array}{cccc|cccc} {{\gamma}_{{{x}^{1}}{{x}^{1}}}} & \cdots & {{\gamma}_{{{x}^{1}}{{x}^{d-1}}}} & {{\gamma}_{{{x}^{1}}z}} & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \cdots & \vdots \\ {{\gamma}_{{{x}^{d-1}}{{x}^{1}}}} & \cdots & {{\gamma}_{{{x}^{d-1}}{{x}^{d-1}}}} & {{\gamma}_{{{x}^{d-1}}z}} & 0 & 0 & \cdots & 0 \\ {{\gamma}_{z{{x}^{1}}}} & \cdots & {{\gamma}_{z{{x}^{d-1}}}} & {{\gamma}_{zz}} & 0 & 0 & \cdots & 0 \\ \hline 0 & \cdots & 0 & 0 & {{\gamma}_{ww}} & 0 & \ldots & 0 \\ 0 & \cdots & 0 & 0 & 0 & {{\gamma}_{ww}} & \ldots & 0 \\ \vdots & \cdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & {{\gamma}_{ww}} \end{array}\right),\end{eqnarray} \tag{ A.7 }$$

and we shall denote the upper left quadrant by the matrix M_ij. For simplicity, we will use the index $\hat{i}$ to denote ${{x}^{{\hat{i}}}}$ in this section, so e.g. cofactors ${{C}_{12}}={{C}_{{{x}^{1}}{{x}^{2}}}}$ and ${{C}_{1z}}={{C}_{{{x}^{1}}z}}$. Similarly the indices $i,\,j,...$ will stand for the xⁱ, i.e. include the z component.

We can now write the cofactor of an element in the top left quadrant of ${{\gamma}_{IJ}}$ as

$$\begin{eqnarray}&&{{C}_{ij}}={{(-1)}^{i+j}}\gamma _{ww}^{\eta}\det \left({{M}_{kl\left\{k\ne j,l\ne i\right\}}}\right)\end{eqnarray} \tag{ A.8 }$$

where $\eta =D-d-1$ and the notation $\det \left({{M}_{kl\left\{k\ne j,l\ne i\right\}}}\right)$ denotes the determinant of the matrix obtained by crossing out the jth row and ith column of M_kl. 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 use this expression for C_ij and the determinant of the right hand side of equation (A.7),

$$\begin{eqnarray}\begin{array}{*{35}{l}} \det {{\gamma}_{IJ}} & =\gamma _{ww}^{\eta}\det {{\gamma}_{ij}} \\ {} & \overset{\ast}{{=}}\,\gamma _{ww}^{\eta}{{\gamma}_{zz}}\det \left({{M}_{kl\left\{k\ne z,l\ne z\right\}}}\right), \end{array}\end{eqnarray} \tag{ A.9 }$$

in order to obtain expressions for inverse metric components according to

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\gamma}^{ij}}=\frac{{{C}_{ij}}}{\det {{\gamma}_{IJ}}}. \end{array}\end{eqnarray} \tag{ A.10 }$$

For d  =  3, this procedure starts from the spatial metric

$$\begin{eqnarray}{{\gamma}_{IJ}}=\left(\begin{array}{ccc|ccc} {{\gamma}_{xx}} & {{\gamma}_{xy}} & {{\gamma}_{xz}} & 0 & \cdots & 0 \\ {{\gamma}_{yx}} & {{\gamma}_{yy}} & {{\gamma}_{yz}} & 0 & \cdots & 0 \\ {{\gamma}_{zx}} & {{\gamma}_{zy}} & {{\gamma}_{zz}} & 0 & \cdots & 0 \\ \hline 0 & 0 & 0 & {{\gamma}_{ww}} & \ldots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & {{\gamma}_{ww}} \end{array}\right).\end{eqnarray} \tag{ A.11 }$$

The components C_ij of the cofactor matrix (which is symmetric) are given by

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{C}_{xx}}=\gamma _{ww}^{\eta}\left({{\gamma}_{yy}}{{\gamma}_{zz}}-\gamma _{yz}^{2}\right),\, & {{C}_{xy}}=-\gamma _{ww}^{\eta}\left({{\gamma}_{yx}}{{\gamma}_{zz}}-{{\gamma}_{zx}}{{\gamma}_{yz}}\right),\, & {{C}_{xz}}=\gamma _{ww}^{\eta}\left({{\gamma}_{yx}}{{\gamma}_{zy}}-{{\gamma}_{zx}}{{\gamma}_{yy}}\right), \\ \cdots & {{C}_{yy}}=\gamma _{ww}^{\eta}\left({{\gamma}_{xx}}{{\gamma}_{zz}}-\gamma _{zx}^{2}\right), & {{C}_{yz}}=-\gamma _{ww}^{\eta}\left({{\gamma}_{xx}}{{\gamma}_{zy}}-{{\gamma}_{zx}}{{\gamma}_{xy}}\right), \\ \cdots & \cdots & {{C}_{zz}}=\gamma _{ww}^{\eta}\left({{\gamma}_{xx}}{{\gamma}_{yy}}-\gamma _{xy}^{2}\right), \end{array}\end{eqnarray} \tag{ A.12 }$$

the determinant becomes

$$\begin{eqnarray}\begin{array}{*{35}{l}} \det {{\gamma}_{IJ}} & =\gamma _{ww}^{\eta}\left({{\gamma}_{xx}}{{\gamma}_{yy}}{{\gamma}_{zz}}+2{{\gamma}_{xy}}{{\gamma}_{xz}}{{\gamma}_{yz}}-{{\gamma}_{xx}}\gamma _{yz}^{2}-{{\gamma}_{yy}}\gamma _{xz}^{2}-{{\gamma}_{zz}}\gamma _{xy}^{2}\right) \\ {} & \overset{\ast}{{=}}\,\gamma _{ww}^{\eta}{{\gamma}_{zz}}\left({{\gamma}_{xx}}{{\gamma}_{yy}}-\gamma _{xy}^{2}\right), \end{array}\end{eqnarray} \tag{ A.13 }$$

and the inverse metric follows by inserting these into equation (A.10).

Using these techniques, we can regularise all terms in equations (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)  
  $$\begin{eqnarray}&&\frac{\delta _{z}^{i}-{{\gamma}^{zi}}{{\gamma}_{ww}}}{z}\end{eqnarray}$$

  We express ${{\gamma}^{zi}}$ in terms of the metric, and trade divisions by z for derivatives ${{\partial}_{z}}$ and obtain

  $$\begin{eqnarray}\frac{\delta _{z}^{i}-{{\gamma}^{zi}}{{\gamma}_{ww}}}{z}\overset{\ast}{{=}}\,\left\{\begin{array}{*{35}{l}} {\sum}_{\hat{j}=1}^{d-1}{{(-1)}^{\hat{i}+\hat{j}}}\frac{{{\gamma}_{ww}}}{\det {{M}_{mn}}}{{\partial}_{z}}{{\gamma}_{z\hat{j}}}\det \left({{M}_{kl\left\{k\ne \hat{i},k\ne z,l\ne z,l\ne \hat{j}\right\}}}\right)~ & ~\text{if}~i=\hat{i} \\ 0 & ~\text{if}~i=z \end{array}\right. \end{eqnarray} \tag{ A.14 }$$

  For the d  =  3 case this reduces to

  $$\begin{eqnarray}\frac{\delta _{z}^{i}-{{\gamma}^{zi}}{{\gamma}_{ww}}}{z}\overset{\ast}{{=}}\,\left\{\begin{array}{*{35}{l}} \frac{{{\gamma}_{yy}}{{\partial}_{z}}{{\gamma}_{xz}}-{{\gamma}_{xy}}{{\partial}_{z}}{{\gamma}_{yz}}}{{{\gamma}_{xx}}{{\gamma}_{yy}}-\gamma _{xy}^{2}}~~ & ~\text{if}~i=x \\ \frac{{{\gamma}_{xx}}{{\partial}_{z}}{{\gamma}_{yz}}-{{\gamma}_{xy}}{{\partial}_{z}}{{\gamma}_{xz}}}{{{\gamma}_{xx}}{{\gamma}_{yy}}-\gamma _{xy}^{2}}~~ & ~\text{if}~i=y \\ 0 & ~\text{if}~i=z \end{array}.\right. \end{eqnarray} \tag{ A.15 }$$
• (2)  
  $$\begin{eqnarray}&&\frac{{{\partial}_{i}}{{\gamma}_{jz}}-{{\delta}_{jz}}{{\partial}_{i}}{{\gamma}_{ww}}}{z}-{{\delta}_{iz}}\frac{{{\gamma}_{jz}}-{{\delta}_{jz}}{{\gamma}_{ww}}}{{{z}^{2}}}+\frac{{{\partial}_{j}}{{\gamma}_{iz}}-{{\delta}_{iz}}{{\partial}_{j}}{{\gamma}_{ww}}}{z}-{{\delta}_{jz}}\frac{{{\gamma}_{iz}}-{{\delta}_{iz}}{{\gamma}_{ww}}}{{{z}^{2}}}\end{eqnarray}$$

  Here we simply trade divisions by z for ${{\partial}_{z}}$ and obtain

  $$\begin{eqnarray}&&\frac{{{\partial}_{i}}{{\gamma}_{jz}}-{{\delta}_{jz}}{{\partial}_{i}}{{\gamma}_{ww}}}{z}-{{\delta}_{iz}}\frac{{{\gamma}_{jz}}-{{\delta}_{jz}}{{\gamma}_{ww}}}{{{z}^{2}}}+\frac{{{\partial}_{j}}{{\gamma}_{iz}}-{{\delta}_{iz}}{{\partial}_{j}}{{\gamma}_{ww}}}{z}-{{\delta}_{jz}}\frac{{{\gamma}_{iz}}-{{\delta}_{iz}}{{\gamma}_{ww}}}{{{z}^{2}}}\\ ~~\overset{\ast}{{=}}\,\left\{\begin{array}{*{35}{l}} 2{{\partial}_{z}}{{\partial}_{(\hat{i}}}{{\gamma}_{\hat{j})z}}~~ & ~\text{if}~i=\hat{i},~j=\hat{j} \\ 0 & ~\text{if}~(i,j)=\left(\hat{i},z\right)~\text{or}~\left(z,\hat{j}\right) \\ {{\partial}_{z}}{{\partial}_{z}}\left({{\gamma}_{zz}}-{{\gamma}_{ww}}\right) & ~\text{if}~i=j=z \end{array}\right. . \end{eqnarray} \tag{ A.16 }$$
• (3)  
  $$\begin{eqnarray}&&-\frac{1}{2}\frac{{{\partial}_{z}}{{\gamma}_{ij}}}{z}+\frac{{{\delta}_{z(i}}{{\gamma}_{j)z}}-{{\delta}_{iz}}{{\delta}_{jz}}{{\gamma}_{ww}}}{{{z}^{2}}}\end{eqnarray}$$

  We use ${{\gamma}_{zz}}-{{\gamma}_{ww}}\overset{\ast}{{=}}\,\mathcal{O}\left({{z}^{2}}\right)$ and trade a division by z for a z derivative. The result is

  $$\begin{eqnarray}-\frac{1}{2}\frac{{{\partial}_{z}}{{\gamma}_{ij}}}{z}+\frac{{{\delta}_{z(i}}{{\gamma}_{j)z}}-{{\delta}_{iz}}{{\delta}_{jz}}{{\gamma}_{ww}}}{{{z}^{2}}}\overset{\ast}{{=}}\,\left\{\begin{array}{*{35}{l}} -\frac{1}{2}{{\partial}_{z}}{{\partial}_{z}}{{\gamma}_{\hat{i}\hat{j}}}~~~ & ~\text{if}~i=\hat{i},~j=\hat{j} \\ 0 & ~\text{if}~(i,j)=\left(\hat{i},z\right)~\text{or}~\left(z,\hat{j}\right) \\ -\frac{1}{2}{{\partial}_{z}}{{\partial}_{z}}{{\gamma}_{ww}} & ~\text{if}~i=j=z \end{array}.\right. \end{eqnarray} \tag{ A.17 }$$
• (4)  
  $$\begin{eqnarray}&&\frac{{{\gamma}_{ww}}{{\gamma}^{zj}}{{\partial}_{j}}{{\gamma}_{ww}}}{z}\end{eqnarray}$$

  Using equations (A.7)–(A.10), we express the inverse metric components ${{\gamma}^{zj}}$ in terms of the downstairs metric and trade the division by z for a z derivative. We thus obtain

  $$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{{{\gamma}_{ww}}{{\gamma}^{zj}}{{\partial}_{j}}{{\gamma}_{ww}}}{z}\overset{\ast}{{=}}\, & \underset{\hat{j}=1}{\overset{d-1}{\sum}}\,\underset{\hat{m}=1}{\overset{d-1}{\sum}}\,{{(-1)}^{\hat{m}+\hat{j}-1}}\frac{{{\gamma}_{ww}}}{\det \left({{M}_{pq}}\right)}{{\partial}_{{\hat{j}}}}{{\gamma}_{ww}}{{\partial}_{z}}{{\gamma}_{z\hat{m}}}\det \left({{M}_{kl\left\{k\ne \hat{j},k\ne z,l\ne z,l\ne \hat{m}\right\}}}\right) \\ {} & +\frac{{{\gamma}_{ww}}\det \left({{M}_{kl\left\{k\ne z,l\ne z\right\}}}\right)}{\det \left({{M}_{pq}}\right)}{{\partial}_{z}}{{\partial}_{z}}{{\gamma}_{ww}}. \end{array}\end{eqnarray} \tag{ A.18 }$$

  which in the case d  =  3 reduces to

  $$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{{{\gamma}_{ww}}{{\gamma}^{zj}}{{\partial}_{j}}{{\gamma}_{ww}}}{z}\overset{\ast}{{=}}\, & \frac{\left({{\gamma}_{yx}}{{\partial}_{z}}{{\gamma}_{zy}}-{{\gamma}_{yy}}{{\partial}_{z}}{{\gamma}_{zx}}\right){{\partial}_{x}}{{\gamma}_{ww}}+\left({{\gamma}_{yx}}{{\partial}_{z}}{{\gamma}_{zx}}-{{\gamma}_{xx}}{{\partial}_{z}}{{\gamma}_{yz}}\right){{\partial}_{y}}{{\gamma}_{ww}}}{{{\gamma}_{xx}}{{\gamma}_{yy}}-\gamma _{xy}^{2}} \\ {} & +{{\partial}_{z}}{{\partial}_{z}}{{\gamma}_{ww}}. \end{array}\end{eqnarray} \tag{ A.19 }$$
• (5)  
  $$\begin{eqnarray}&&\frac{{{\gamma}^{zz}}\gamma _{ww}^{2}-{{\gamma}_{ww}}}{{{z}^{2}}}\end{eqnarray}$$

  The regularisation of this term proceeds in analogy to that of term (9) in appendix B of [83], except we do not set $\det \gamma =1$. By rewriting $1={{\gamma}^{zz}}/{{\gamma}^{zz}}={{\gamma}^{zz}}\det {{\gamma}_{IJ}}/{{C}_{zz}}$, trading divisions by z for z derivatives and using ${{\gamma}_{zz}}\overset{\ast}{{=}}\,{{\gamma}_{ww}}+\mathcal{O}\left({{z}^{2}}\right)$, we obtain

  $$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{{{\gamma}^{zz}}\gamma _{ww}^{2}-{{\gamma}_{ww}}}{{{z}^{2}}}\overset{\ast}{{=}}\, & \frac{1}{2}{{\partial}_{z}}{{\partial}_{z}}\left({{\gamma}_{ww}}-{{\gamma}_{zz}}\right) \\ {} & -\frac{{{\sum}_{\hat{i}=1}^{d-1}}\,{{\sum}_{\hat{m}=1}^{d-1}}\,{{(-1)}^{\hat{i}+\hat{m}-1}}{{\partial}_{z}}{{\gamma}_{z\hat{i}}}{{\partial}_{z}}{{\gamma}_{z\hat{m}}}\det \left({{M}_{kl\left\{k\ne \hat{i},k\ne z,l\ne z,l\ne \hat{m}\right\}}}\right)}{\det \left({{M}_{kl\left\{k\ne z,l\ne z\right\}}}\right)}. \end{array}\end{eqnarray} \tag{ A.20 }$$

  which in the case d  =  3 reduces to

  $$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{{{\gamma}^{zz}}\gamma _{ww}^{2}-{{\gamma}_{ww}}}{{{z}^{2}}}\overset{\ast}{{=}}\, & \frac{1}{2}{{\partial}_{z}}{{\partial}_{z}}\left({{\gamma}_{ww}}-{{\gamma}_{zz}}\right) \\ {} & +\frac{-2{{\gamma}_{xy}}{{\partial}_{z}}{{\gamma}_{xz}}{{\partial}_{z}}{{\gamma}_{yz}}+{{\gamma}_{xx}}{{\left({{\partial}_{z}}{{\gamma}_{yz}}\right)}^{2}}+{{\gamma}_{yy}}{{\left({{\partial}_{z}}{{\gamma}_{xz}}\right)}^{2}}}{{{\gamma}_{xx}}{{\gamma}_{yy}}-\gamma _{xy}^{2}}. \end{array}\end{eqnarray} \tag{ A.21 }$$
• (6)  
  $$\begin{eqnarray}&&\frac{{{K}_{iz}}-{{\delta}_{iz}}{{K}_{ww}}}{z}\end{eqnarray}$$

  The division by z is again traded for a derivative if $i\ne z$ and for i  =  z, we use ${{K}_{zz}}={{K}_{ww}}+\mathcal{O}\left({{z}^{2}}\right)$, so that

  $$\begin{eqnarray}\frac{{{K}_{iz}}-{{\delta}_{iz}}{{K}_{ww}}}{z}\overset{\ast}{{=}}\,\left\{\begin{array}{*{35}{l}} {{\partial}_{z}}{{K}_{\hat{i}z}}~~~ & ~\text{if}~i=\hat{i} \\ 0 & ~\text{if}~i=z \end{array}.\right. \end{eqnarray} \tag{ A.22 }$$

In this section, we discuss how the set of spatial normal frame vectors, equation (4.38), can be recast in a form suitable for applying Gram–Schmidt orthonormalisation. It turns out to be convenient to first rescale the ${{\tilde{m}}_{(\alpha )}}$ such that they would acquire unit length in a flat spacetime with spatial metric ${{\delta}_{IJ}}$. Denoting these rescaled vectors with a caret, we have

$$\begin{eqnarray}&&{{\hat{m}}_{(\alpha )}}=\frac{1}{\sqrt{\left(\underset{s=\alpha}{\overset{D-1}{\sum}}\,w_{s}^{2}\right)\,\left(\underset{s=\alpha -1}{\overset{D-1}{\sum}}\,w_{s}^{2}\right)}}\left(\begin{array}{c} \left.\begin{array}{c} 0~ \\ \vdots ~ \\ 0~ \end{array}\right\}(\alpha -2)\times \\ -{\sum}_{s=\alpha}^{D-1}{{\left({{w}^{s}}\right)}^{2}} \\ \left.\begin{array}{c} {{w}^{\alpha -1}}{{w}^{\alpha}} \\ \vdots \\ {{w}^{\alpha -1}}{{w}^{D-2}} \\ {{w}^{\alpha -1}}{{w}^{D-1}} \end{array}~\right\}(D-\alpha )\times \end{array}\right),~\alpha =2,\,\ldots,\,D-1.\end{eqnarray} \tag{ B.1 }$$

Recall that we formally set ${{w}^{1}}\equiv {{x}^{1}},\,\ldots,\,{{w}^{d-1}}\equiv {{x}^{d-1}}$, ${{w}^{d}}\equiv z$. As a convenient shorthand, we define

$$\begin{eqnarray}&&\rho _{I}^{2}\equiv \underset{s=I}{\overset{D-1}{\sum}}\,{{\left({{w}^{s}}\right)}^{2}},\end{eqnarray} \tag{ B.2 }$$

so that, for instance, $\rho _{1}^{2}={{r}^{2}}$, $\rho _{4}^{2}={{\left({{w}^{4}}\right)}^{2}}+\ldots +{{\left({{w}^{D-1}}\right)}^{2}}$, ${{\rho}_{D-1}}={{w}^{D-1}}$. This definition allows us to write

$$\begin{eqnarray} &&\hat{m}_{(\alpha )}^{I}=\frac{1}{{{\rho}_{\alpha}}\,{{\rho}_{\alpha -1}}}\left(\begin{array}{c} \left.\begin{array}{c} 0~ \\ \vdots ~ \\ 0~ \end{array}\right\}(\alpha -2)\times \\ -\rho _{\alpha}^{2} \\ \left.\begin{array}{c} {{w}^{\alpha -1}}{{w}^{\alpha}} \\ \vdots \\ {{w}^{\alpha -1}}{{w}^{D-2}} \\ {{w}^{\alpha -1}}{{w}^{D-1}} \end{array}~\right\}(D-\alpha )\times \end{array}\right).\end{eqnarray} \tag{ B.3 }$$

We can now express the angles ${{\phi}^{\alpha}}$ in terms of the radial variables ${{\rho}_{I}}$,

$$\begin{eqnarray}&&\sin {{\phi}^{\alpha}}=\frac{{{\rho}_{\alpha}}}{{{\rho}_{\alpha -1}}},~~~~\cos {{\phi}^{\alpha}}=\frac{{{w}^{\alpha -1}}}{{{\rho}_{\alpha -1}}}.\end{eqnarray} \tag{ B.4 }$$

Using these relations in (B.3), we obtain

$$\begin{eqnarray} &&\hat{m}_{(\alpha )}^{I}=\left(\begin{array}{c} \left.\begin{array}{c} 0~ \\ \vdots ~ \\ 0~ \end{array}\right\}(\alpha -2)\times \\ -\sin {{\phi}^{\alpha}} \\ \begin{array}{c} \cos {{\phi}^{\alpha}}\,\cos {{\phi}^{\alpha +1}} \\ \cos {{\phi}^{\alpha}}\,\sin {{\phi}^{\alpha +1}}\,\cos {{\phi}^{\alpha +2}} \\ \vdots \\ \cos {{\phi}^{\alpha}}\,\sin {{\phi}^{\alpha +1}}\,\ldots \,\sin {{\phi}^{D-2}}\,\cos {{\phi}^{D-1}} \\ \cos {{\phi}^{\alpha}}\,\sin {{\phi}^{\alpha +1}}\,\ldots \,\sin {{\phi}^{D-2}}\,\sin {{\phi}^{D-1}} \end{array}~ \end{array}\right)=\left(\begin{array}{c} \left.\begin{array}{c} 0~ \\ \vdots ~ \\ 0~ \end{array}\right\}(\alpha -2)\times \\ -\sin {{\phi}^{\alpha}} \\ \begin{array}{c} \vdots \\ \cos {{\phi}^{\alpha}}\left({\prod}_{s=\alpha +1}^{\alpha +n-1}\sin {{\phi}^{s}}\right)\cos {{\phi}^{\alpha +n}} \\ \vdots \end{array} \end{array}\right),\end{eqnarray} \tag{ B.5 }$$

where $n=1,\,\ldots,\,D-\alpha$, and we formally set $\cos {{\phi}^{D-1}}\equiv 1$ and ${\prod}_{s=\alpha +1}^{\alpha}\sin {{\phi}^{\alpha}}\equiv 1$.

Now, in our computational domain $\rho _{d+1}^{2}=0$, which, from the definition of our coordinate system in equation (2.8) gives

$$\begin{eqnarray}&&{{r}^{2}}{{\sin}^{2}}{{\phi}^{2}}\,\ldots \,{{\sin}^{2}}{{\phi}^{d+1}}=0\end{eqnarray} \tag{ B.6 }$$

Since ${{\phi}^{2}},\,\ldots,\,{{\phi}^{d}}$ are arbitrary in our computational domain, we must have either ${{\phi}^{d+1}}=0~\text{or}~\pi$. Without loss of generality, we choose ${{\phi}^{d+1}}=0$, which fixes the d  −  1 vectors

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\hat{m}}}_{(2)}}= & \left(-\sin {{\phi}^{2}},~\cos {{\phi}^{2}}\,\cos {{\phi}^{3}},~\ldots,~\cos {{\phi}^{2}}\underset{s=3}{\overset{d}{\prod}}\,\sin \left({{\phi}^{s}}\right),\,\underset{(D-d-1)\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,\right). \end{array}\end{eqnarray} \tag{ B.7 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \vdots & \\ {{{\hat{m}}}_{\left(\hat{\alpha}\right)}}= & \left(\underset{\left(\hat{\alpha}-2\right)\times}{{\underbrace{{0,\,\ldots,\,0}}\,}}\,,\,-\sin {{\phi}^{{\hat{\alpha}}}},\,\cos {{\phi}^{{\hat{\alpha}}}}\cos {{\phi}^{\hat{\alpha}+1}},\,\ldots,\,\cos {{\phi}^{{\hat{\alpha}}}}\underset{s=\hat{\alpha}+1}{\overset{d}{\prod}}\,\left(\sin {{\phi}^{s}}\right),\,\underset{(D-d-1)\times}{{\underbrace{{0,\,\ldots,\,0}}\,}}\,\right) \end{array}\end{eqnarray} \tag{ B.8 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \vdots & \\ {{{\hat{m}}}_{(d)}}= & \left(\underset{(d-2)\times}{{\underbrace{{0,\,\ldots,\,0}}\,}}\,,\,\,-\sin {{\phi}^{d}},\,\cos {{\phi}^{d}},\,\underset{(D-d-1)\times}{{\underbrace{{0,\,\ldots,\,0}}\,}}\,\right), \end{array}\end{eqnarray} \tag{ B.9 }$$

which, up to rescaling by ${{\rho}_{{\hat{\alpha}}}}{{\rho}_{\hat{\alpha}-1}}$, are equal to the vectors in equations (4.40)–(4.42). For the remaining vectors, we can use the rotational freedom in the angles ${{\phi}^{d+2}},\,\ldots,\,{{\phi}^{D-1}}$. Any choice for these values will satisfy ${{w}^{d+1}}=\ldots ={{w}^{D-1}}=0$ as required on our computational domain and we merely need to ensure that we choose these angles such that the resulting set of vectors is orthogonal. This is most conveniently achieved by setting

$$\begin{eqnarray}&&{{\phi}^{d+2}}=\ldots ={{\phi}^{D-1}}=0,\end{eqnarray} \tag{ B.10 }$$

which, inserted into equation (B.5), implies

$$\begin{eqnarray}&&\hat{m}_{(a)}^{I}=\delta _{a}^{I},\,~~~a=d+1,\ldots,D-1.\end{eqnarray} \tag{ B.11 }$$

Combined with equations (B.7)–(B.9) and restoring the tilde in place of the caret on the ${{\tilde{m}}_{(a)}}$, we have recovered equations (4.43) and (4.44) in section 4.2 for the angular vectors. For the case d  =  3 we have just two non-trivial vectors:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\hat{m}}}_{(2)}}= & \left(-\sin {{\phi}^{2}},~\cos {{\phi}^{2}}\,\cos {{\phi}^{3}},~\cos {{\phi}^{2}}\,\sin {{\phi}^{3}},~\underset{(D-4)\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,\right), \end{array}\end{eqnarray} \tag{ B.12 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\hat{m}}}_{(3)}}= & \left(0,~-\sin {{\phi}^{3}},~\cos {{\phi}^{3}},~\underset{(D-4)\times}{{\underbrace{{0,~\ldots,~0}}\,}}\,\right), \end{array}\end{eqnarray} \tag{ B.13 }$$

recovering equations (4.46)–(4.50).