Dimensional reduction in numerical relativity: Modified cartoon formalism and regularization

We present in detail the Einstein equations in the Baumgarte-Shapiro-Shibata-Nakamura formulation for the case of $D$ dimensional spacetimes with $SO(D-d)$ isometry based on a method originally introduced in Ref.1. Regularized expressions are given for a numerical implementation of this method on a vertex centered grid including the origin of the quasi-radial coordinate that covers the extra dimensions with rotational symmetry. Axisymmetry, corresponding to the value $d=D-2$, represents a special case with fewer constraints on the vanishing of tensor components and is conveniently implemented in a variation of the general method. The robustness of the scheme is demonstrated for the case of a black-hole head-on collision in $D=7$ spacetime dimensions with $SO(4)$ symmetry.


Introduction
For most of its history, numerical relativity, i.e., the construction of solutions to Einstein's field equations with numerical methods, was mainly motivated by the modeling of compact objects as sources of gravitational waves 2 for ground based [Laser Interferometer Gravitational-Wave Observatory (LIGO), Virgo] and space based [Laser Interferometer Space Antenna (LISA)] detectors; see for example. 3,4 Following the breakthroughs in black-hole (BH) binary simulations in 2005, 5-7 however, the field rapidly expanded into a variety of new physics frontiers. 8 These applications often involve higher-dimensional spacetimes where BHs are known from (semi-)analytic studies to exhibit a richer phenomenology as for exam-ple through the existence of topologically non-spherical horizons or gravitational instabilities. 9 Despite considerable progress through analytic, perturbative and numerical methods, our understanding of the properties of these higher-dimensional BHs is still a long way from the level of maturity obtained in the four-dimensional case. Yet, applications of numerical relativity to BHs in D > 4 have already revealed a plethora of exciting results.
Critical spin parameters have been identified above which Myers-Perry a BHs become unstable to bar mode perturbations and migrate to more slowly spinning BHs via GW emission. 11,12 The celebrated Gregory-Laflamme instability 13 has been shown to lead to the formation of naked singularities in finite asymptotic time in numerical simulations of black strings in D = 5 dimensions. 14 Most recently, a similar behaviour has been identified in evolutions of thin black rings demonstrating the first violation of cosmic censorship for a generic type of asymptotically flat initial data; 15 see also 16 for a perturbative study. Applications of the gauge-gravity duality often consider D = 5 dimensional BHs in asymptotically Anti de Sitter (AdS) spacetimes such that the dual Conformal Field Theory (CFT) lives on the D = 4 dimensional (conformal) boundary of the spacetime. Applications of this AdS/CFT correspondence include the thermalization of quark-gluon plasma, turbulence or jet quenching in heavy-ion collisions; see [17][18][19][20][21][22] and references therein.
BH collisions provide fertile ground for numerical relativity in higher dimensions. First, we obtain unprecedented insight into the dynamics of general relativity in its most violent, non-linear regime. Furthermore, the so-called TeV gravity scenarios provide solutions to the hierarchy problem in terms of large (∼ sub millimeter) extra dimensions [23][24][25] or infinite extra dimensions with a warp factor 26,27 which are accessible to gravity but no other standard-model interactions. If correct, these theories open up the possibility of BH formation in particle collisions or in cosmic ray showers. [28][29][30] Valuable input for the analysis of experimental data at the Large Hadron Collider includes the scattering cross-section and energy loss in GWs. Numerical relativity has provided us with a rather comprehensive understanding of these collisions in D = 4, 31-39 which the community now starts extending to higher D. [40][41][42][43] For more details on these new areas in numerical relativity research see. 8 Numerical simulations of BH spacetimes in higher dimensions are a challenging task. First and foremost this is simply a consequence of the required computational resources. Simulations in D = 4 require of the order of O(10 2 ) cores and O(10 2 ) Gb of memory. Each extra spatial dimension introduces an additional factor of O(10 2 ) grid points and correspondingly more memory and floating point operations. Even with modern high-performance computing systems, this sets practical limits on the feasibility of accurately evolving higher-dimensional spacetimes. At the same time, many of the outstanding questions can be addressed by imposing symmetry assumptions on the spacetimes in question such as planar symmetry in modeling asymptotically AdS spacetimes, 18 cylindrical symmetry for black strings 14 or dif-paper Dimensional reduction in numerical relativity: Modified cartoon formalism and regularization 3 ferent types of rotational symmetries. 15 This can be achieved in practice by either (i) using a specific form of the line element that directly imposes the symmetry in question (see e.g. 18 ), (ii) starting with a generic line element and applying dimensional reduction through isometry (see e.g. [44][45][46] or (iii) implementing the symmetry through a so-called Cartoon method. 47 Here we are concerned with the latter approach and, more specifically, with a modification thereof originally introduced in 1 (see also 12,48,49 ) which we will henceforth refer to as the modified Cartoon method.
This article is structured as follows. In Sec. 2, we introduce the notation used throughout our work, and illustrate the modified Cartoon implementation of the symmetries for a specific example. In Sec. 3 we introduce the Baumgarte-Shapiro-Shibata-Nakamura 50, 51 (BSSN) evolution system we use for the Einstein equations, and derive their specific form in SO(D − d) symmetry when rotational symmetry is present in ≥ 2 planes which corresponds to d < D − 2. The axisymmetric case d = D − 2 imposes less restrictive conditions on the vanishing of tensor density components and their derivatives and the particulars of its numerical implementation are discussed in Sec. 4. As an example, we present in Sec. 5 numerical simulations of a BH head-on collision in D = 7 dimensions employing SO(4) symmetry. We summarize our findings in Sec. 6 and include in three appendices a list of important relations for the components of tensors and derived quantities as well as the regularization necessary at the origin in the quasi-radial direction. It is instructive to illustrate the method by considering first a simpler scenario: axisymmetry in three spatial dimensions. Let (x, z, w) denote Cartesian coordinates and assume rotational symmetry about the x axis b i.e., there exists a rotational Killing field in the z, w plane. Evidently, the geometry of such a three-dimensional manifold can be constructed straightforwardly provided all tensors (e.g. the metric) are known on the semi infinite plane w = 0, z ≥ 0, x ∈ R. We note the simplification in the computational task: the w coordinate has dropped out and the quasi-radius z takes on only non-negative values, reducing an originally three-dimensional computational domain to a calculation on half of R 2 . This is the case considered in the original papers. 1, 47 The most common applications will likely consider higher-dimensional spacetimes with SO(D − 3) symmetry, but here we present the general application to a Let us then consider a D dimensional spacetime consisting of a manifold M and a metric g AB of signature D − 2 where A, B, . . . = 0, . . . , D − 1. We assume the spacetime to obey SO(D − d) symmetry and introduce Cartesian coordinates by symmetry implies the existence of rotational Killing vectors in each plane spanned by two of the coordinates (z, w a ). In complete analogy with the axisymmetric scenario discussed above, it is now sufficient to provide data on the d-dimensional semi-infinite hyperplane w a = 0, xî ∈ R, z ≥ 0. The components of a tensor at any point in the spacetime can then be obtained by appropriately rotating data from the hyperplane onto the point in question.
In modeling spacetimes with such symmetries, it is therefore entirely sufficient to compute data on the hyperplane which largely solves the problem of increased computational cost mentioned in the introduction. There remains, however, the difficulty that the Einstein equations, irrespective of the specific formulation one chooses, contain derivatives of tensor components in the w a directions which cannot be evaluated numerically in the usual fashion, as for example using finite differences or spectral methods. Furthermore, the number of tensor components present in the Einstein equations still increases rapidly with the dimension parameter D resulting in a substantial increase of memory requirements and floating point operations. Both of these difficulties are overcome by exploiting the conditions imposed on the tensor components and their derivatives by the SO(D − d) symmetry. It is these conditions which we address next. It turns out to be convenient in this discussion to distinguish between (1) the case d = D − 2 corresponding to SO(2) isometry, and (2) all remaining cases d < D − 2. We defer discussion of the special case d = D − 2 to Sec. 4 where we present a numerical treatment specifically designed for conveniently dealing with it. The description of this treatment will be simpler after first handling the class of symmetries with d < D − 2 which we discuss in the remainder of this and in the next section.

Tensor components in
The key ingredient we use in reducing the number of independent tensor components and relating their derivatives are the rotational Killing vectors and the use of coordinates adapted to the integral curves of these Killing vectors. The method is best introduced by considering a concrete example. Let ξ denote the Killing vector field corresponding to the rotational symmetry in the (z, w) plane, where w ≡ w a for some fixed number a ∈ {d + 1, . . . , D − 1}. We introduce a new coordinate system that replaces (z, w) with cylindrical coordinates and leaves all other coordinates paper Dimensional reduction in numerical relativity: Modified cartoon formalism and regularization 5 unchanged,X In these coordinates, the Killing field is ξ = ∂ ϕ and the vanishing of the Lie derivative L ξ g AB = 0 implies ∂ ϕ g AB = 0. Note that quantities constructed from the spacetime metric directly inherit this property. This applies, in particular, to the Arnowitt-Deser-Misner 52 (ADM) -see also 53,54 -and the BSSN variables widely used in numerical relativity. For d < D − 2, one can furthermore show that the ϕ component of a vector field and those components of a tensor field T AB , where exactly one index is ϕ, vanish c . The concrete example we now discuss in more detail concerns a symmetric tensor density T AB of weight W and, in particular, the mixed components T iw , where the index i stands for any one of the (xî, z) coordinates and w stands for one of the w a . We first consider the components Tî w for some fixed value ofî. Transforming the componentTî ϕ to Cartesian coordinates, one gets where D is the Jacobian det(∂X A /∂X B ) = ρ. Using thatTî ϕ = 0 by symmetry, this equation implies Similarly, transformingT ρϕ to Cartesian coordinates and using thatT ρϕ = 0 by symmetry, one straightforwardly gets Recalling that the computational domain is the hyperplane w a = 0, xî ∈ R, z ≥ 0, we conclude from Eqs. (6),(7) that on the computational domain T iw = 0. This argument holds for any specific choice of the coordinate w, so that we conclude To compute the derivatives with respect to w on the w = 0 hyperplane, one can proceed as follows. For the tensor components in the example above, one can simply use (6) and (7) to calculate ∂ w T ia and then set w = 0. Alternatively, writing the Killing field ξ as c Here the case d = D − 2 represents an exception; an axisymmetric, toroidal magnetic field, for example, satisfies SO(2) symmetry, but has a non-vanishing ϕ component.
and imposing the vanishing of the Lie derivative L ξ T ia = 0 on the w = 0 hyperplane, one gets Repeating this process for all components of scalar, vector and rank 2 tensor densities as well as their first and second derivatives, we get the relations summarized in Appendix A.
We have shown the calculation here explicitly for the case of tensor densities. It can be shown that the vectorial expressions thus obtained also apply to the contracted Christoffel symbol Γ A ≡ g MN Γ A MN constructed from the metric, even though it is not a vector density.

Dimensional reduction of the BSSN equations
In this section, we will apply the symmetry relations obtained above to the specific case of the BSSN formulation of the Einstein equations in D spacetime dimensions. We emphasize, however, that the procedure spelled out here for the BSSN system can be applied in similar form to any of the alternative popular formulations used in numerical relativity.

The D dimensional BSSN equations
The starting point for the BSSN formulation is a space-time, or (D − 1) + 1, split where the spacetime is foliated in terms of a one-parameter family of D − 1 dimensional, spatial hypersurfaces. In coordinates adapted to this split, the line element takes on the form where I, J, . . . = 1, . . . , D − 1 and α and β I denote the lapse function and shift vector, respectively. The ADM equations in the form developed by York 53 then result in one Hamiltonian constraint, D − 1 momentum constraints and D(D − 1)/2 second-order evolution equations for the spatial metric components γ IJ . The latter are formulated as a first-order-in-time system by introducing the extrinsic curvature Space-time decompositions of the Einstein equations typically split the energy momentum tensor analogously into time, space and mixed components according to where n A denotes the future pointing, timelike unit normal field on the spatial hypersurfaces. The complete set of the ADM equations, thus obtained, can be found as Eqs. (52)-(55) in. 8 paper Dimensional reduction in numerical relativity: Modified cartoon formalism and regularization 7 The BSSN system is obtained from the ADM equations by applying a conformal transformation to the spatial metric, a trace split of the extrinsic curvature and promotion of the contracted spatial Christoffel symbols to the status of evolution variables. The BSSN variables are defined as where γ = det γ IJ , andΓ I M N are the Christoffel symbols associated with the conformal metricγ IJ . We formulate here the conformal factor in terms of the variable χ, following. 7 Alternative versions of the equations using variables W ≡ √ χ or φ ≡ −(ln χ)/4 can be found in. 55,56 Note that the definition of the BSSN variables in (14) implies two algebraic and one differential constraints given bỹ The D dimensional BSSN equations are then given by the Hamiltonian and momentum constraints and the evolution system Here, R IJ , R are the Ricci tensor and scalar associated with the physical spatial metric γ IJ , Λ is the cosmological constant, the superscript "TF" denotes the tracefree part and we have added a constraint damping term σG I in the last line, following the suggestion by. 57 The above equations are complemented by the following auxiliary relations, The BSSN equations in this form are general and facilitate the numerical construction of D dimensional spacetimes. Next, we will describe in detail how the equations can be reduced to an effective system in d spatial dimensions for spacetimes obeying rotational symmetry with d < D − 2.  We first note that the spatial metric with SO(D − d) symmetry has the form which simplifies the calculation of the inverse metricγ AB ; see Appendix B.
The constraint equations (16), (17) become and the BSSN evolution equations (18)- (22) are now written as These equations contain a number of auxiliary expressions which are given in terms of the fundamental BSSN variables by Eq. (29) as well as paper Dimensional reduction in numerical relativity: Modified cartoon formalism and regularization 11 The BSSN equations in this form can readily be implemented in an existing "d+1" BSSN code with the addition of merely two new field variables,γ ww andÃ ww . While the BSSN equations acquire additional terms, the computational domain remains d-dimensional. Furthermore, the entire set of Eqs. There only remains one further subtlety arising from the explicit division by z in several of the terms present. Some (though not all) numerical codes require evaluation of these expressions at z = 0 which makes regularization of these terms mandatory. As we show explicitly in Appendix B, this can be achieved for all terms, yielding expressions that are exact in the limit z → 0. The results we discuss in Sec. 5 make use of these regularized terms on the plane z = 0 demonstrating that this procedure provides stable and accurate evolutions.
We conclude this section with a brief remark of the matter terms present in (31)- (51) in the form of the projections ρ, j i , S ij and S = χ(γ ij S ij + nγ ww S ww ) of the energy-momentum tensor. The specific form of these terms will depend on the physical system under consideration and will need to be evaluated separately for each case as will the precise form of the matter evolution equations resulting from the conservation law ∇ A T AB . Many applications of higher-dimensional numerical relativity concern BHs and the example application discussed in Sec. 5 will be an asymptotically flat vacuum spacetime where the matter terms and the cosmological constant are zero.

The special case of SO(2) symmetry
We now return to the special situation where d = D −2 which corresponds to SO(2) isometry, i.e. axisymmetric spacetimes. 1 This case is special in that ϕ components of a vector field or those components of a tensor field T AB where exactly one index is ϕ, do not necessarily vanish. The reason for this exceptional property of SO(2) symmetry is that with only one Killing vector ∂ ϕ , the vector f (x i ) ∂ ϕ , for an arbitrary function f , trivially satisfies the symmetry as it commutes with all Killing vectors. Thus, SO(2) symmetry does not, in general, cause any tensor components to vanish, and only a negligible amount of computational cost and memory would be saved by explicitly inserting the modified Cartoon terms, as derived in the previous two sections, into the BSSN equations. We are still able, however, to capitalize on the substantial reduction in memory and floating point operations that arises from the dimensional reduction of the computational domain. This is most conveniently achieved by retaining the BSSN equations in their full D-dimensional form and only using the modified Cartoon method to fill in derivatives that cannot be calculated directly on the computational grid.
Let us illustrate this process for D = 5 with SO(2) isometry. Due to the symme-try, we can model such a spacetime on a three-dimensional grid on which we store all vector and tensor components, including those of the type V w and T Iw , which do not vanish for SO(2) symmetry. In order to evolve the system by one time step, we need to compute derivatives with respect to all coordinates. Derivatives with respect to x i can be calculated on the grid using standard methods. Derivatives with respect to w, on the other hand, can be calculated on our three-dimensional grid using the modified Cartoon method. For example, for a vector, the procedure is Appendix C lists all necessary expressions for derivatives with SO(2) symmetry. Once all derivatives have been calculated, we have all the information required to use the standard D = 5 BSSN equations (18)(19)(20)(21)(22) without the need for any extra terms. Note that this method for handling SO(2) symmetry can straightforwardly be combined with the method described in Secs. 2, 3. Such a procedure can handle, for example, the symmetry of black rings and has been applied in 15 to speed up the exploration of the gauge parameter space in numerical evolutions of black rings in D = 5. Black rings have horizons of topology 58 S 1 × S 2 and rotate along the S 1 . This rotational symmetry requires handling with the special method for d = D − 2 because ϕ components do not vanish in that case. The second symmetry corresponding to the S 2 , however, is amenable to the treatment presented in Secs. 2 and 3. In practice, Refs. 15, 59 first applied the latter reduction and then the special SO(2) reduction sketched in Eq. (52).

Application to a black-hole collision
In this section we present, as a specific example for the efficacy of the formalism, results from the numerical simulation of a head-on collision of two non-spinning BHs in D = 7 dimensions starting from rest. A non-rotating BH in D spacetime dimensions is described by the Tangherlini 60 solution where dΩ D−2 denotes the area element of the (D − 2) sphere and the parameter µ is related to the BH mass M and the horizon radius R h by Here, Ω D−2 is the surface area of the unit (D − 2) sphere. The Tangherlini solution can be written in isotropic coordinates in the form which facilitates construction of analytic data for a snapshot of a spacetime containing multiple BHs at the moment of time symmetry according to the procedure of Brill and Lindquist. 61 These higher-dimensional Brill-Lindquist data are given in terms of the ADM variables by where the summation over A and K extend over the multiple BHs and spatial coordinates, respectively, and X K A denotes the position of the A th BH. We have implemented these initial data in the Lean code, 62, 63 which is based on Cactus 64, 65 and uses Carpet 66, 67 for mesh refinement. The specific case presented here has been obtained using SO(4) symmetry, i.e. D = 7, d = 3, for a collision along the x axis of two equal-mass BHs initially separated by 7.58 R h , where R h is the horizon radius associated with a single BH with µ = µ 1 = µ 2 . The computational domain consists of a set of seven refinement levels, the innermost two centered on the BHs and the five outer ones on the origin. We employ standard moving puncture gauge conditions 68 [note that we use here β a = 0 in accordance with Eq. (A.3)] having initialized lapse and shift to their Minkowski values α = 1, β i = 0. Two simulations have been performed in octant symmetry with a grid spacing ∆x = R h /52 and ∆x = R h /104, respectively, on the innermost level, that increases by a factor of two on each consecutive level further out. Figure 1 shows the trajectories of the two BHs evolving in time from the initial separation through merger into a single hole centered on the origin, obtained from the high resolution simulation with ∆x = R h /104. In order to check the consistency of our numerical formalism, we have also analyzed the constraint equations for this configuration. A snapshot of the Hamiltonian constraint (16) along the collision axis at evolution time t = 80 R h is shown in Fig. 2. In this figure, the result obtained for the high resolution run has been amplified by a factor of four expected for second-order convergence. The overlap of the two curves demonstrates convergence at second order, compatible with the numerical scheme that employs second and fourth-order accurate discretization and interpolation techniques. We have performed the same analysis for the Hamiltonian and momentum constraints at several points in time and observe the same second-order convergence of both constraints throughout infall and merger. Note that only one BH is present on the computational domain (at about x = 2.5 in the figure) because of the octant symmetry. The other BH is represented in this simulation through the symmetric boundary conditions imposed at x = 0.

Conclusions
In the presence of rotational symmetry, the Einstein equations simplify considerably and the generation of numerical solutions to these equations can be implemented with significant improvements in computational cost and the required amount of computer memory. The Cartoon method proposed in 47 was the first technique designed with the particular goal of efficiently modeling axisymmetric spacetimes in 3+1 numerical relativity. A modification, often dubbed the modified Cartoon method 1 used relations between tensor components in place of spatial interpolation operations, which not only eliminates the need of introducing a few extra grid points in the symmetry directions, but also allows for a particularly convenient generalization to an arbitrary number of spacetime dimensions and number of rotational symmetries. 12,48,49 In this work, we have presented in detail the complete set of equations as obtained for the BSSN formulation of the Einstein equations in D spacetime dimensions with SO(D − d) isometry where d ∈ {1, 2, . . . , D − 2}. Furthermore, we explicitly demonstrate the presence of extra terms for the case d = D − 2, where the symmetry condition allows for a wider class of components of tensor densities to remain non-zero. Finally, we have compiled a list of terms involving division by the quasi-radial coordinate (the z direction in our case) and illustrate how all irregularities at the origin z = 0 can be cured through equivalence with manifestly regular expressions. Even though we used the BSSN formulation for our discussion, the recipes detailed here can be applied straightforwardly to other popular formulations of the Einstein equations such as the generalized harmonic gauge 1, 69 or the conformal Z4 70, 71 systems.
As an example, we have presented results from a head-on collision from rest of two equal-mass, non-spinning BHs in D = 7 spacetime dimensions. Following a rather slow acceleration phase, due to the rapid diminishing of the gravitational force with distance, the two BHs merge and we observe second-order convergence of the constraints. This confirms in yet another type of application the remarkable robustness observed for the modified Cartoon method in applications to spinning BHs 12 or high-energy collisions in D = 5. 42 This seemingly superior robustness as compared with the method of reduction by isometry developed in 45 is, at present, empirical but merits further investigation at the analytic level. The presence of z in the denominator of several terms in the system of Eqs. (31)-(51) merely arises from the quasi-radial nature of the coordinate z and can be handled straightforwardly in analogy to the treatment of the origin in spherical or axisymmetry.
We will present the regularized terms needed in the generic SO(D−d) symmetry; however, it should be noted that terms involving the inverse metric become much more complicated for a large d, and so we will also explicitly show these terms for the most common case, d = 3.
We first require that all components expressed in a fully Cartesian set of coordinates are regular. A well known consequence of this assumption is that tensor density components containing an odd (even) number of radial, i.e. z, indices contain only odd (even) powers of z in a series expansion around z = 0. The same holds for quantities derived from tensors and densities such as the BSSN variableΓ i .
Next, we consider the inverse metric which we obtain through inversion of the matrix equation (30). By constructing the cofactor matrix and dividing by the determinant, we obtain, for d = 3 where we introduced the symbol " * = " to denote equality in the limit z → 0. The components for the inverse BSSN metric in Eq. (B.1) simplify accordingly.
For a general d we know that the matrix takes the form given in Eq (30). Then, denoting the cofactor matrix for a given element ofγ IJ as C IJ , the inverse BSSN metric components are (note that the metric is symmetric, so that C IJ = C J I ) · · · · · · · · ·γ zz = Czz detγIJ .

(B.3)
Again, in the BSSN case detγ IJ = 1, and the inverse metric element is simply the cofactor of that element. For simplicity, we will use indicesî in place of xî in the remainder of this section, so that, for example C 12 ≡ C x 1 x 2 , C 1z ≡ C x 1 z etc. When used without a caret, the lower case Latin indices i, j, . . . also include the z component. Here, the notation det(M kl{k =j,l =i} ) denotes the determinant of the matrix 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. The next regularity condition we require our spacetime to satisfy is the absence of a conical singularity at z = 0. In polar coordinates (ρ, ϕ) constructed as in Sec. 2.2, this condition can be expressed asγ ϕϕ = ρ 2γ ρρ which translates into the conditionsγ We thus arrive at the following list of regularized terms valid in the limit z → 0.
(1) By expanding β z = b 1 z + b 3 z 3 + . . ., and likewise forΓ z and ∂ zγww , we obtain and likewise for α or χ in place ofγ ww in the last expression.
(2) We express the inverse metric components through their cofactors, given for arbitrary d by Eq. (B.4), and then apply the same trading of divisions by z for derivatives as done for β z /z, to obtain (B.8) Here, as well in items (5) and (9) below, we formally set det(M jl{j =z,j =m,l =i,l =z} ) = 1 for the case d = 2 where no entries would be left in the matrix after crossing out two rows and columns. For d = 1, the case i =î does not arise which obviates the need to evaluate the determinant. .
(B.17) (9) The term (γ zzγ ww − 1)/z 2 requires slightly more work and we describe its derivation here in a little more detail. We first rewrite this term in the form and express the inverse metric componentγ zz in terms of the corresponding cofactor matrix component and the determinant as