Cauchy Evolution of Asymptotically Global AdS Spacetimes with No Symmetries

We present the first proof-of-principle Cauchy evolutions of asymptotically global AdS spacetimes with no imposed symmetries, employing a numerical scheme based on the generalized harmonic form of the Einstein equations. In this scheme, the main difficulty in removing all symmetry assumptions can be phrased in terms of finding a set of generalized harmonic source functions that are consistent with AdS boundary conditions. In four spacetime dimensions, we detail an explicit set of source functions that achieves evolution in full generality. A similar prescription should also lead to stable evolution in higher spacetime dimensions, various couplings with matter fields, and on the Poincare patch. We apply this scheme to obtain the first long-time stable 3+1 simulations of four dimensional spacetimes with a negative cosmological constant, using initial data sourced by a massless scalar field. We present preliminary results of gravitational collapse with no symmetry assumptions, and the subsequent quasi-normal mode ringdown to a static black hole in the bulk, which corresponds to evolution towards a homogeneous state on the boundary.


Introduction
In recent years anti-de Sitter (AdS) space has proven to be a particularly exciting theoretical laboratory for studying the strong-field regime of General Relativity (GR). AdS with reflective boundary conditions plays the role of a box that naturally keeps propagating waves confined to its interior, where they are perpetually interacting. Thus, even the smallest perturbations in AdS can enter the strong-field regime, where qualitatively new gravitational phenomena emerge. One of the most important of these is gravitational collapse -the growth of curvatures that eventually leads to the formation of a singularity in spacetime potentially associated with a black hole. Obtaining the details of this fundamental process in full generality in AdS is still an open problem. In asymptotically flat spacetimes, although it has not yet been proven rigorously, this process of gravitational collapse is expected to generically end in a rotating black hole that is characterized by two conserved numbers: total mass and total angular momentum. In asymptotically AdS spacetimes, the endpoint is less clear. Small, rapidly rotating black holes are unstable due to a process known as superradiance -the amplification of waves that scatter off a rotating object. Along with the box-like nature of AdS, this amplification leads to a runaway process whose endpoint is unknown.
In an unprecedented way, the simulation of asymptotically AdS spacetimes has also opened up the field of numerical relativity to the study of phenomena in areas beyond the traditional astrophysical setting. At the heart of this push to understand AdS is a deep connection between gravity in AdS to certain conformal field theories (CFT), now known as the AdS/CFT correspondence [1][2][3]. Through this connection, the study of AdS spacetimes has become immediately relevant to fundamental questions in many areas in physics, such as fluid dynamics [4][5][6], relativistic heavy ion collisions [7][8][9][10], and superconductivity [11][12][13]. See, for example, [14][15][16][17] for excellent reviews. The reason why the study of AdS is crucial for our understanding of these phenomena is that AdS/CFT provides an important -and in most cases the only -window into the real-time dynamics of strongly interacting quantum field theories far from equilibrium. The dynamical far-from-equilibrium strongly interacting regime is precisely the one that is least explored and understood, and the one that has the best chance of making contact with certain experiments.
Our current understanding of gravity in AdS remains limited for several reasons. First, evolution in AdS is notoriously hard, in part because it is an initial-boundary value problem whose systematic study is still in its infancy. Cauchy evolution in AdS requires data to be prescribed not only at an initial spacelike hypersurface, but also at spatial and null infinity which constitute the timelike boundary of an asymptotically AdS spacetime. Second, the most interesting phenomena involve spacetimes that have very little or no symmetry, making these evolutions beyond the reach of most numerical codes. Third, for many of these phenomena, there is a variety of physical scales that must be adequately resolved to correctly capture the relevant physics.
The main purpose of this article is to present the first proof-of-principle Cauchy evolution of asymptotically AdS spacetimes that has been achieved with no symmetry assumptions, and to describe the framework that makes Cauchy evolution in AdS possible in full generality. The results presented here are based on a code with adaptive mesh refinement (AMR) capabilities that solves the Einstein equations in generalized harmonic form for asymptotically AdS spacetimes, subject to reflective (i.e., Dirichlet) boundary conditions. We couple gravity to a massless scalar field, but the latter does not play any fundamental role in our scheme; we introduce it as a convenient mechanism to arrange for initial data whose future Cauchy development contain trapped surfaces.
Ingoing characteristic (e.g., Eddington-Finkelstein) coordinates have been successfully used to simulate dynamical spacetimes containing black branes in asymptotically AdS space-times in Poincaré coordinates in full generality, i.e., no symmetry assumptions. 1 This method has been applied to a variety of settings and by now the literature on the subject is vast and we will not review it here. We refer the reader to [20] for a detailed review. This approach, however, will fail if the ingoing radial null geodesics form caustics within the numerical domain, which can happen whenever there is a strong localized perturbation of the background spacetime. For instance, the dynamical formation of localized black holes in the background of the AdS soliton spacetime [21] or even a localized black hole falling through the Poincaré horizon of AdS are just two possible examples where the ingoing coordinates of [20] are likely to become singular due to the formation of caustics. 2 On the other hand, Cauchy evolution in conjunction with generalized harmonic coordinates is well-known to successfully handle strong, highly dynamical and localized gravitational fields, such as those produced by the individual black holes in a binary. Whilst it is possible that many problems that have been solved using ingoing coordinates in the Poincaré patch of AdS can also be solved with Cauchy evolution, the latter can be applied to situations where ingoing coordinates will almost certainly fail. Furthermore, the use of Cauchy evolution benefits from the infrastructure developed over many years to numerically solve the black hole binary problem in general relativity [24][25][26]. In particular, the code described in the present work has built-in AMR and is designed to run in large supercomputing clusters; both of these features will likely turn out to be crucial in solving certain key open problems in AdS.
A key requirement for obtaining stable evolution in AdS is a gauge choice that is consistent with the conditions imposed at the AdS boundary (see, for example, [27]). In most cases, a gauge choice leading to stable numerical evolution is typically found in spacetimes with a certain degree of symmetry. In the present work, we detail a gauge choice in D = 4 spacetime dimensions that leads to stable evolution in an asymptotically global AdS setting with no symmetry assumptions. This work is a direct precursor to fully general studies of gravitational collapse and black hole formation in AdS. In this context, Cartesian coordinates are suitable as they are regular everywhere, do not contain coordinate singularities, and do not have the well-known limitation suffered by spherical coordinates in the form of severely shorter time steps imposed by the Courant-Friedrichs-Lewy (CFL) condition. In addition, most AMR infrastructures are designed for this type of coordinates. Similar coordinates were used in [28] to study the non-spherically symmetric collapse of a massless scalar field in global AdS 5 with SO(3) symmetry. In anticipation of fully general studies, we choose to write our prescription in terms of global Cartesian coordinates, using second order finite difference derivative stencils to discretize the initial constraint equations and the evolution equations. The framework we present here straightforwardly generalizes to other settings and other discretization schemes.
The rest of this article is organized as follows. In Section 2 we describe the setup, starting with a short review of anti-de Sitter spacetime, and two complementary characterizations of asymptotically AdS boundary conditions. In Section 3 we detail our prescription for obtaining stable Cauchy evolution with no symmetries in Cartesian coordinates. The crucial ingredients for this perscription are reflective Dirichlet boundary conditions imposed on appropriate evolution variables, and a specific choice of generalized harmonic source functions. In Section 4 we define boundary quantities whose evolution describes the physics at the AdS boundary. In Section 5 we outline the generalized harmonic scheme that we use in our simulations. Section 6 contains preliminary results of simulations of gravitational collapse with no symmetry assumptions. We conclude with a discussion in Section 7. We have relegated some technical details to several appendices. In Appendix A we write down the Einstein equations in harmonic coordinates. In Appendix B, we follow our prescription for the interesting case of global AdS in spherical coordinates and we obtain the corresponding stable gauge. In Appendix C, we do the same for the Poincaré patch. Appendix D contains a description of our construction of initial data for the class of spacetimes considered in the paper, while in Appendix E we provide the details of our complete gauge choice, including the bulk. In Appendix F we explain how we carry out the extrapolation to read off the boundary quantities. Some convergence tests are presented in Appendix G. Throughout, we use geometric units where Newton's constant is set to G = 1 and the speed of light is set to c = 1.

Anti-de Sitter Spacetime
The dynamics of gravity with a cosmological constant Λ in four dimensions coupled to a real massless scalar field ϕ can be described by the following action: where R is the Ricci scalar of the metric g αβ with determinant g. The variation of the action (2.1) with respect to g αβ and ϕ gives the equations of motion We then recast (2.2) into generalized harmonic form. See Appendix A for the explicit form of the resulting equations that we evolve, and [29] for more details about the theoretical aspects of the formulation. The numerical solution we obtain is given in terms of the spacetime metric g αβ , the scalar field ϕ and a choice of gauge source functions H α .
The metric of AdS 4 is the maximally symmetric vacuum (i.e., ϕ = 0) solution of (2.2) and (2.3) in four dimensions. In terms of global coordinates that cover the whole spacetime, given by (t, r, θ, φ) ∈ (−∞, +∞) × (0, +∞) × [0, π) × [0, 2π), this metric can be expressed aŝ with a characteristic length scale L that is related to the cosmological constant by Λ = −3/L 2 , and where dΩ 2 2 = dθ 2 + sin 2 θdφ 2 is the metric of the round unit 2-sphere. A crucial feature of this spacetime is the presence of a timelike boundary at r → +∞, which makes stable evolution of initial data possible only if boundary conditions are imposed on the evolved fields. In other words, any Cauchy problem in this setting is an initial-boundary value problem.
To proceed further, we first compactify r = 2ρ/(1 − ρ 2 / 2 ) so that the AdS boundary at r → +∞ is at a finite value of the new radial coordinate, ρ = . We hereafter set = 1 without loss of generality, so that the AdS boundary is at ρ = 1. In this way, we obtain (compactified) spherical coordinates x α = (t, ρ, θ, φ). Defining a convenient function f (ρ) = (1 − ρ 2 ) 2 + 4ρ 2 /L 2 , the metric of AdS 4 in this set of coordinates readŝ Second, we make use of Cartesian coordinates x µ = (t, x, y, z) defined by x = ρ cos θ, y = ρ sin θ cos φ, z = ρ sin θ sin φ, where ρ = x 2 + y 2 + z 2 . This allows us to bypass the severe restriction that would be imposed on the time step size near ρ = 0 on a grid in spherical coordinates. This brings the metric of AdS 4 into its final form, suitable for our numerical simulations: (2.6)

Asymptotically anti-de Sitter Spacetimes
We will be interested in the Cauchy evolution of asymptotically AdS spacetimes. In this section we present a review of two different characterizations of such spacetimes and the relation between them, specializing to the case of D = 4 spacetime dimensions for concreteness. In doing so, we will also be able to write down the boundary conditions for asymptotically AdS spacetimes in terms of these two different characterizations.
Let us start from the original arguments presented in [30]. The authors implicitly considered spacetimes (M, g) that admit a conformal compactification, and thus a definition of conformal boundary ∂M . Then they define asymptotically AdS spacetimes by requiring that the spacetime asymptotically approaches the pure AdS solution. More precisely, for any set of global coordinates x α , the authors required the deviation of the full metric g αβ from the pure AdS metricĝ αβ , given by h αβ = g αβ −ĝ αβ , to satisfy three conditions: (i) It is consistent with the asymptotic decay of the Kerr-AdS metric near ∂M in this set of coordinates.
(ii) Its fall-off near ∂M is invariant under the global AdS symmetry group O(3, 2), i.e., near the boundary ∂M for any generator X of O(3, 2).
(iii) The surface integral charges associated with the generators of O(3, 2) are finite.
In addition, for the purposes of this article, we restrict this definition to spacetimes that satisfy the Einstein equations (2.2).
It is important to recognize that conditions (i), (ii) and (iii) can be condensed into one. Ref. [30] already shows that the explicit fall-off satisfying (i) and (ii) automatically implies (iii). Furthermore, requiring (ii) is sufficient to obtain the fall-off near the boundary that satisfies also (i) and (iii). This can be seen from the results of [31], in which (2.7) is solved in any spacetime dimension and the 4-dimensional case coincides with the fall-offs in [30].
The condition (ii) amounts to a full spacetime metric g αβ that approaches the pure AdS metricĝ αβ near ∂M . This has two consequences for the terminology commonly used in the literature, as well as in this work. First, we can refer to ∂M as the AdS boundary because it has the same conformal structure as the boundary of pure AdS, i.e, R × S 2 topology, and metric given by that of the Einstein Static Universe. Second, we can define certain classes of coordinates in terms of the corresponding fall-offs of the metric components near the boundary as follows. Given a set of coordinates x α in which the pure AdS metric components areĝ αβ , we denote by x α all sets of coordinates in which the full metric components g αβ approach the pure AdS metric components in the formĝ αβ . For example, we will denote any set of coordinates in which the metric g asymptotes toĝ in the form (2.5) by (t, ρ, θ, φ) and we will refer to them as spherical coordinates. Similarly, we will denote any set of coordinates in which g asymptotes toĝ in the form (2.6) by (t, x, y, z) and we will refer to them as Cartesian coordinates. 3 The fall-offs for the metric obtained by [30] can thus be written in the form for arbitrary functions f αβ (t, θ, φ). These are supplemented by the fall-offs for the scalar field, given in [31]. Here we restrict the discussion to a massless scalar field ϕ with a fast fall-off that preserves the asymptotics (2.8), for which for arbitrary f (t, θ, φ). In Cartesian coordinates, these fall-offs read for arbitrary f µν and f , and where ρ = ρ(x, y, z).
The fall-offs of the source functions, involved in the generalized harmonic formulation employed in this study, can be deduced from (2.8) through the definition In spherical coordinates, denoting the pure AdS values byĤ α , (2.8) and (2.12) imply for arbitrary f α . In Cartesian coordinates, denoting the pure AdS values byĤ µ , (2.10) and (2.12) imply for arbitrary f µ and ρ = ρ(x, y, z).
A different characterization of asymptotically AdS spacetimes can be given in terms of the well-known Fefferman-Graham (FG) expansion [32]. In this approach, one starts with the definition of a locally asymptotically AdS spacetime (M, g) as a spacetime that admits a conformal compactification, thus allowing the definition of a conformal boundary ∂M , and that satisfies the Einstein equations (2.2). No assumption is made at this stage on the topology of the boundary. The FG theorem states that one can always find a coordinate system xᾱ = (t,z,θ,φ) in a neighbourhood of the boundary for which the boundary is at z = 0 and the metric can be written in the form Then, the near-boundary (i.e., aboutz = 0) expansion of the Einstein equations completely determines the coefficient g (2)āb in terms of g (0)āb . Therefore the dynamics that makes this spacetime differ from pure AdS appears at orderz 3 in the expansion of gāb. If we make the further requirement that the topology of the boundary is the same as in the pure AdS case, i.e., R×S 2 , the spacetime becomes globally asymptotically AdS and this characterization becomes equivalent to the one obtained from the original arguments in [30]. The FG form (2.15) of the metric immediately provides the near-boundary behaviour and shows that coordinates can be defined so that thezz component of any asymptotically AdS metric goes as 1/z 2 , and thezā components vanish in a neighbourhood of the AdS boundary.
We conclude this section by showing an explicit example of how FG coordinates can be found in the case of general asymptotically AdS 4 spacetimes. In what follows, we set the characteristic length scale to L = 1 for simplicity. We start from the general form for the asymptotically AdS metric in spherical coordinates x α , given by g αβ =ĝ αβ + h αβ . The deviations h αβ from the pure AdS metricĝ αβ have fall-offs that are given by the asymptotically AdS boundary conditions (2.8). Defining z = 2(1 − ρ)/(1 + ρ), we can bring the pure AdS metric (2.5) into the FG form. Since z asymptotes to 1 − ρ near the AdS boundary ρ = 1, we can use (2.8) to immediately write down the metric fall-offs in terms of our new coordinate z. The metric in these coordinates reads where the coefficients f αβ in the expansion above are functions of (t, θ, φ). Notice that the metric in (2.17) is not in the FG form yet because the zz component is not 1/z 2 up to the desired order in z, and the tz, zθ, zφ components do not vanish up to O(z 2 ). Defininḡ which can be inverted near the boundary as we finally obtain the metric in FG form: where now the coefficients f αβ are functions of (t,θ,φ). Notice that f ρρ has been reabsorbed in gtt, gθθ, gφφ. From the form of the metric in (2.26), we could use the holographic renormalization prescription of [33] to read off the boundary CFT stress tensor. See Appendix B.4 for more details.

Boundary Prescription
In this section, we present our prescription to obtain a choice of generalized harmonic gauge source functions that achieves stable evolution. We choose to do so using Cartesian coordinates, as they provide a suitable chart to evolve points near the centre of the grid, which is necessary when analyzing gravitational collapse and black hole formation. This procedure generalizes in a straightforward manner to other asymptotically AdS spacetimes in D ≥ 4 spacetime dimensions, different coupling with matter fields, and coordinates on global AdS or on the Poincaré patch. We consider the application to spherical coordinates in Appendix B, and to the Poincaré patch in Appendix C. We impose asymptotically AdS boundary conditions (2.10), (2.11), (2.14) as reflective Dirichlet boundary conditions on appropriate evolution variables, as explained in the next section. For a discussion in a simpler context with more symmetry, see [27].

Evolution Variables and Boundary Conditions
The boundary conditions on asymptotically AdS spacetimes, discussed in Section 2.2, can be imposed as Dirichlet boundary conditions at the AdS boundary. This requires appropriately defining and evolving a new set of variables, from which the full solution (g µν , ϕ, H µ ) can be subsequently reconstructed. Here, we define evolution variables in the Cartesian coordinates employed by our numerical scheme. Later, in Section 4 we will show expressions for quantities at the AdS boundary in spherical coordinates. In Appendix B, we explicitly show how these spherical variables relate to our Cartesian evolution variables.
The Cartesian metric evolution variables,ḡ µν , are defined by first considering the deviation from pure AdS in Cartesian coordinates, h µν = g µν −ĝ µν , then stripping h µν of as many factors of (1 − ρ 2 ) as needed so that each component falls off linearly in (1 − ρ) near the AdS boundary at ρ = 1. 4 We see from (2.10) that in four dimensions, the metric evolution variablesḡ µν that satisfy these requirements are simplȳ Similarly, the Cartesian boundary condition on the scalar field (2.11) suggests that we use the evolution variableφ Finally, the boundary conditions (2.14) on H µ suggest the use of For evolved variables defined in this way, the boundary conditions (2.10), (2.11), (2.14) can be easily imposed as Dirichlet boundary conditions at the AdS boundary: (3.4)

Gauge Choice for Stability
Coordinates over the entire spacetime are fully determined only once we choose the gauge source functions H µ . In Cartesian coordinates, as can be seen from (2.14), H µ are fixed up to order 1 − ρ by its pure AdS valuesĤ µ in an expansion near the AdS boundary. As we shall see, the choice of H µ at the next order in this expansion, (1 − ρ) 2 , cannot be completely arbitrary if we wish to achieve stable evolution. A specification of generalized harmonic source functions at order (1 − ρ) 2 that provides stable Cauchy evolution can be obtained following the procedure detailed in this section.
The first step involves expanding the evolved variables,ḡ µν ,H µ andφ, in a power series about (1 − ρ) ≡ q = 0. By construction, these evolved variables are linear in q at leading order:ḡ where all the coefficients are functions of the coordinates (t, x, y, z) on the boundary ρ(x, y, z) ≡ x 2 + y 2 + z 2 = 1 (or q(x, y, z) = 0). We now substitute these variables into the evolution equations (A.3), and we expand each component in powers of q. The three lowest orders, q −2 , q −1 , q 0 , are fixed by the pure AdS metricĝ which itself is a solution of (A.3), so these terms vanish trivially. The remaining orders vanish only ifḡ µν ,H µ ,φ are a solution of (A.3).
We are now interested in identifying the order of q at which the second derivatives of g (1)µν with respect to (t, x, y, z) appear. For each component, we denote their combination by˜ ḡ (1)µν , i.e., (1)µν , (3.8) for some functions c t µν , c x µν , c y µν , c z µν of (t, x, y, z) at ρ(x, y, z) = 1. 5 These derivative terms are included in the first piece of (A.3), namely in − 1 2 g ρσḡ µν,ρσ . From this, we can easily find their order of q by recalling that the leading order of the inverse metric is given by its purely AdS piece, g µν = O(ĝ µν ) = O(q 2 ), andḡ (1)µν is multiplied by q in the near-boundary expression ofḡ µν (see eq. (3.5)). Thus,˜ ḡ (1)µν must appear in the coefficient of order q 3 for every component of (A.3). 6 In other words, each component of the expansion of (A.3) near q = 0 can be written in the schematic form: or, rearranging the terms in order to obtain wave-like equations, Similar arguments show that terms involving the massless scalar field with the fast fall-off that we have chosen in (2.9) appear at the next order with respect to the leading order in˜ ḡ (1)µν . This holds in any dimension and in any set of coordinates, and it implies that the details of the matter sector, e.g., the value of the mass of a matter field, do not affect the results of the prescription presented here, since only the lowest order coefficients in the expansion of the Einstein equations are relevant.
We now write the coefficients in (3.10) explicitly, including only the q −2 terms. The near-boundary expansion is most easily obtained by first writing the Cartesian coordinates (x, y, z) in terms of the boundary-adapted spherical coordinates (q, θ, φ), and then expanding near q = 0. We find: where the coordinates (q, θ, φ) should be understood as functions of (x, y, z). All that remains is to write down the generalized harmonic constraints C µ ≡ H µ − x µ = 0 at leading order in the same near-boundary expansion. We get: In the generalized harmonic formulation, choosing a gauge amounts to choosing a set of generalized harmonic source functionsH µ for the entire evolution. Although we expect that many gauge choices are allowed, [27] mentions a few that do not give rise to stable evolutions. We now present a procedure that provides the stable gauge in our Cartesian simulations. We believe that our prescription provides a stable gauge in a variety of settings of physical interest, such as higher spacetime dimensions, various couplings to matter fields, different types of global coordinates or Poincaré coordinates. Thus, it enables numerical Cauchy evolution in AdS in full generality, that is, with no symmetry assumptions. The steps that lead to our stable gauge, in a form that can be easily applied to all previously mentioned cases, are the following.
1. Solve the leading order of the near-boundary generalized harmonic constraints forH (1)µ .
For example, in the Cartesian case, the leading orders of (3.21)-(3.24) vanish for: 2. Let N (1) be the lowest order in q appearing in the near-boundary expansions of all the˜ ḡ (1)µν . Plug the source functions obtained in step 1 into the q N (1) terms of the near-boundary expansions˜ ḡ (1)µν . This gives a number of independent equations that, together with their derivatives, ensure tracelessness and conservation of the boundary stress-energy tensor (see Section 4). 7 Solve these equations for an equal number of metric coefficientsḡ (1)µν and their derivatives. In the Cartesian case, N (1) = −2 and there is only one independent equation given bȳ which we can solve, for instance, in terms ofḡ (1)tt .
3. Plug the solutions to the equations in step 2 into the gauge obtained in step 1. In Cartesian coordinates, using (3.26) to eliminateḡ (1)tt from (3.25), we havē This is the asymptotic gauge condition that we have empirically verified leads to stable 3+1 evolution of asymptotically AdS 4 spacetimes in Cartesian coordinates. Other choices of asymptotic source functions may enjoy similar stability properties. The choice ofH µ in the bulk is still completely arbitrary and the functional form that we implement in our simulations is detailed explicitly in Appendix E.
The rationale for this procedure is as follows. Recall that if C µ = 0 and ∂ t C µ = 0 are satisfied at t = 0, 8 and the boundary conditions are consistent with C µ = 0 being satisfied at the boundary for all time. Then, at the analytical level, the generalized harmonic constraint C µ = 0 remains satisfied in the interior for all time. The addition of constraint damping terms to the Einstein equations, eq. (A.3), helps to ensure that deviations at the level of the discretized equations remain under control. Thus, in solving the expanded system of equations (A.3), we are assured that only the subset of solutions that are also solutions of the Einstein equations are being considered. With this in mind, the near-boundary form of (A.3), given by (3.9), implies that our task in obtaining a solution is to satisfy A (i)µν = 0 for all i, and for some choice of source function variablesH µ . This task is significantly eased by picking a gauge, through a suitable choice ofH µ , that eliminates A (1)µν , i.e., the lowest order of the expansion of the Einstein equations near the AdS boundary. This is precisely what the above set of steps is designed to do, and it is why we did not stop at the gauge obtained in Step 1, (3.25), which would have resulted in a gauge that does not explicitly set A (1)µν = 0.
Finally, it is also important to develop an understanding of the reason why the choice of H µ is not completely free. Although identifying every cause for the instability of a simulation is usually very complicated, one practical reason is clear and can be understood with the following example in Cartesian coordinates. Suppose we choose a gauge in which, after some time t > t 0 ,H (1)t takes the valuē where b t ∈ R is a possibly vanishing constant. According to (3.25), the requirement that C t = 0 now impliesḡ (1)tz = b t . Even though this condition does not violate any of the requirements above, it is an additional Dirichlet boundary condition that must be imposed for t > t 0 if we hope to find a solution for this example. 9 Although imposing boundary conditions that change with time is of interest in certain studies motivated by the AdS/CFT correspondence, for simplicity we do not consider such cases in this article. It should be straightforward to generalize our prescription for time-dependent boundary conditions.

Boundary Stress Tensor
In the simulations we output the holographic stress-energy tensor of the dual CFT. In this section, we obtain the analytic expression for this object in spherical coordinates x α = (t, ρ, θ, φ), as they are adapted to the metric of the AdS boundary in global coordinates. Thus, in order to obtain their numerical values, we will have to convert the evolution variables in Cartesian coordinatesḡ µν provided by our numerical scheme into their counterpartsḡ αβ in spherical coordinates. We do this in Appendix B, through the transformation (B.2).
Let us denote by x a = (t, θ, φ) the coordinates on timelike hypersurfaces ∂M q at fixed ρ (or q). To compute the holographic stress-energy tensor of the boundary CFT, T ab CF T , we first compute the quasi-local stress-energy tensor (q) T αβ at ∂M q as prescribed in [34]. We , S α is the spacelike, outward pointing timelike unit vector normal to ∂M q and G αβ is the Einstein tensor of ∂M q . 10, 11 We will be interested in the value of (q) T αβ for q close to 0, i.e., near the AdS boundary. Restricting the indices corresponding to the coordinates x a , we can compute the boundary stress-energy tensor as From (q) T αβ we also compute the total AdS mass as follows [34]. At each time t of evolution, we take a spacelike two-dimensional surface S in ∂M q , with induced metric σ ab = ω ab + u a u b , where u a = −N (dt) a is the future pointing unit 1-form normal to S in ∂M q , lapse N and shift N a . The total AdS mass is then given by The holographic stress-energy tensor can be expressed in terms of the leading order coefficients of the near-boundary expansion ofḡ αβ . We find: 12 10 Notice the different sign in the last term of (4.1) with respect to to [34]. 11 All these tensors, although defined on the tangent space of the spacetime manifold M , are invariant under projection ω α β = δ α β − S α S β onto ∂Mq. Therefore, they can be identified, under a natural (i.e., basisindependent) isomorphism, with tensors defined on the tangent space of ∂Mq. The components of tensors on ∂Mq in coordinates x a is simply given by taking the components of tensors on M in coordinates x α and disregarding every combination of indices that includes an index ρ. See [35] for more details on this correspondence. 12 The expressions (4.4) have a factor of 1/G that corresponds to the large-N scaling of the expectation value of the stress tensor in the boundary 2+1-dimensional CFT. When quoting numerical results, we keep convention of working in geometric units with G = 1.
where we have set the characteristic length scale to L = 1. Similarly, for the total mass in AdS we find We can now use the metric of the AdS boundary, λ ab dx a dx b = −dt 2 + dθ 2 + sin 2 θdφ 2 , to raise one index of T ab CF T and solve the eigenvalue problem T a b CF T v b = Λ v v a at each point along the AdS boundary. In this way, assuming that T ab CF T satisfies the weak energy condition, 13 we obtain the energy density of the boundary CFT, , as minus the eigenvalue associated to the unique (up to rescaling) timelike eigenvector. Similarly, the boundary anisotropy is given by ∆p ≡ |p 1 −p 2 |, where p 1 and p 2 are the eigenvalues associated with, respectively, the remaining two spacelike eigenvectors.
One useful quantity to compute is the trace of the stress-energy tensor, trT CF T = λ ab T ab CF T . We obtain: If we convert the spherical quantities into their Cartesian counterparts we see that trT CF T depends only on the factorḡ (1)tt −ḡ (1)xx −ḡ (1)yy −ḡ (1)zz . We saw in (3.26) that this factor vanishes. This is an important sanity check: we see that tracelessness of the stress-energy tensor, expected for a CFT in 2+1 dimensions, is ensured by the lowest order in the near boundary expansion of the Einstein equations, provided that the generalized harmonic constraints are satisfied. In other words, tracelessness of the boundary stress tensor is, in our scheme, directly tied to how close our numerical solution is to a solution of the Einstein field equations. We check that we are indeed converging to such a solution in Appendix G. In practice, we monitor trT CF T to estimate truncation error. Another important check that we performed is the conservation of the analytic form of T ab CF T . The simplest way to prove this is by using the near-boundary expansion of the Einstein equations in spherical coordinates, as done in Appendix B.

Numerical Scheme
In this section we consider the core elements of the numerical scheme used in this study. We start by discussing the numerical features on which this scheme relies for solving the initial-boundary value problem in AdS. We then describe our apparent horizon finder and the method with which we excise trapped regions.

Numerics of the Initial-Boundary Value Problem
We solve the Einstein equations in generalized harmonic form (A.3) with constraint damping terms, coupled with the massless Klein-Gordon equation (A.5). We obtain asymptotically AdS spacetimes in Cartesian coordinates x µ = (t, x, y, z). The solution is determined in terms of the metric, scalar field and source function variables (ḡ µν ,φ,H µ ) defined in Section 3.1.
We substitute the definitions of these variables, (3.1)-(3.3), in the equations of motion and analytically remove all the purely AdS terms. The resulting PDEs are discretized with second order finite difference derivative stencils, and then integrated in time using an iterative Newton-Gauss-Seidel relaxation procedure with a three time level hierarchy. The source function variablesH µ near the AdS boundary are set as we have prescribed in (3.27), whilst deep in the bulk they are set to zero. In between, we use smooth transition functions to interpolate between the near boundary and the bulk regions, see Appendix E for the details of our full implementation.
We use the PAMR/AMRD libraries [36] for running these simulations in parallel on Linux computing clusters. Although these libraries have adaptive mesh refinement capabilities, numerical evolution is performed on a grid with fixed refinement. The numerical grid is in The typical grid resolution uses N x = N y = N z = 325 points in each of the Cartesian directions, with equal grid spacings ∆x = ∆y = ∆z ≡ ∆.
The time step of evolution is determined by ∆t = λ∆. Although we do not perform a detailed analysis of the stability of our finite difference scheme, the Courant-Friedrichs-Lewy (CFL) condition for stability is expected to be satisfied as long as the CFL factor λ is set to a value well below 1. Thus, we use λ = 0.3. Notice that the most remarkable advantage of using Cartesian coordinates is that the CFL condition does not severely restrict the CFL factor as it would in spherical coordinates, hence allowing simulations to reach large evolution times with modest computational resources. In contrast, spherical coordinates (t, ρ, θ, φ) with fixed resolution ∆ρ, ∆θ, ∆φ would necessitate ∆t = λ min(∆ρ, ρ min ∆θ, ρ min ∆θ∆φ). At points next to the origin, which must be evolved in studies of gravitational collapse and black hole formation, ρ takes its smallest value ρ min = ∆ρ. Hence, in spherical coordinates, ∆t would become prohibitively small for higher resolutions, i.e., for smaller ∆ρ, ∆θ, ∆φ.
The following components play a fundamental role in the numerical implementation of the initial-boundary value problem. Reflective Dirichlet boundary conditions (3.4) are imposed at the AdS boundary ρ = 1. In general the AdS boundary does not lie on Cartesian grid points, so we set boundary conditions at points at most one grid point away from the boundary via interpolation. Referring to Figure 1, for any given evolution variable, we set its value at grid points with ρ < 1 − ∆/2 (i.e., the green dots inside the blue dotted line in this figure) by first order interpolation between the Dirichlet value at boundary points (red dots) and the value at the adjacent point further into the interior ρ < 1 (purple dots). To identify the latter, we move along the Cartesian direction corresponding to the coordinate of the green dot with the largest absolute value. This direction is represented by light blue arrows. Notice that points with ρ ≥ 1 − ∆/2 are excised to avoid issues with quantities that would diverge at ρ = 1. Finally, to obtain the values of quantities at the boundary, needed to extract the holographic observables, we use third order extrapolation from their bulk point values. The details of the implementation in our numerical simulations can be found in Appendix F.
Last but not least, time-symmetric initial data, sourced by a massless real scalar field, are obtained by solving the conformal decomposition of the Hamiltonian constraint (D. 15). The solution to (D.15) is computed, after second order finite discretization, through a full approximation storage (FAS) multigrid algorithm with v-cycling and Newton-Gauss-Seidel relaxation, built into the PAMR/AMRD libraries. We ensure that initial data satisfies the generalized harmonic constraints. See Appendix D for more details and the complete choice of initial data.

Apparent Horizon Finder and Excision
Once the solution is obtained at a certain time t, we can search for the position R(θ, φ) of an apparent horizon (AH). We use the following flow method in spherical coordinates (ρ, θ, φ), obtained in the usual way from the Cartesian coordinates of the solution. We consider n two-dimensional surfaces at constant, equally spaced, values of ρ within a user-specified range included in (0,1), and we pick the one with smallest L 2 -norm of the outward null expansion. Let ρ 0 be the ρ coordinate on this surface. Starting from the initial guess R(θ, φ) = ρ 0 , for any (θ, φ) ∈ [0, π) × [0, 2π) we find the solution to the equation where Θ(ρ, θ, φ)| ρ=R(θ,φ) is the outward null expansion of the two-dimensional surface given by F (ρ, θ, φ) ≡ ρ − R(θ, φ) = 0. We iterate this process with starting point given by the solution R(θ, φ) to (5.1) found in the previous iteration. Assuming that the initial guess ρ 0 is not too distant from the position of the AH, R(θ, φ) is expected to progressively approach the AH after each iteration. This process stops when either the L 2 -norm of Θ(ρ, θ, φ)| ρ=R(θ,φ) is below some specified tolerance, i.e., R(θ, φ) is sufficiently close to the AH, or the user-specified maximum number of iterations has been reached, i.e., either there is no AH at time t or this method was not able to find it.
This AH finder is based on a (θ, φ) grid with equal grid spacings ∆θ = ∆φ = ∆ AH . 14 The outward null expansion at a given AH finder grid point (θ, φ) is obtained by first order interpolation in three dimensions from the values of the expansion at Cartesian grid points that surround (θ, φ). These values are calculated from the definition of outward null expansion once the spacetime metric at time t is known. We observe that a N θ × N φ = 9 × 17 resolution is enough to find the AH in the simulations considered in Section 6 in less than 10 4 iterations. Since (5.1) is a parabolic equation, the "time" step ∆s must be at least of order ∆ 2 AH for stability. When using n = 10 initial trial surfaces and an initial range of ρ values between 0.1 and 0.5, as we do in our simulations, we find that the AH finder works effectively if ∆s takes much smaller values. Specifically, we set ∆s = 10 −4 .
When an AH is found, we excise Cartesian grid points in an ellipsoid included in the AH and centred at the centre of the AH, in order to avoid the formation of geometric singularities in the computational domain. 15 More specifically, the excision ellipsoid has Cartesian semiaxes, a ex x , a ex y , a ex z , determined by a ex , where x AH is the x-coordinate value of the intersection between the AH and the x-axis, and similarly for a ex y and a ex z . We set the excision buffer to δ ex = 0.4. In our simulations, we assume that the characteristics of the equations of motion in the AH region flow towards the origin, although we do not compute the characteristics explicitly. As a consequence, the solution at points inside the AH only evolves to affect points, at later times, that are further inside the AH. In other words, the information needed to solve the equations of motion on and outside the excision surface at a certain time is entirely contained in the numerical domain at previous times. This allows us to solve the equations of motion at the excision surface by employing one-sided stencils that do not reference points inside the excised region, with no need to impose conditions at the excision boundary. By construction, the excised surface is the same for all three time levels involved in the Newton-Gauss-Seidel relaxation for evolution variables at time t. Therefore, we only need to use the one-sided version of the spatial stencils.
It commonly occurs that the excised surface moves during evolution and previously excised points become unexcised. In this case, we initialize the value of newly unexcised points closest to the previous surface using fourth order extrapolated values from adjacent exterior points along each Cartesian direction. We do so for any variable and at all three time levels of the hierarchy. Finally, Kreiss-Oliger dissipation [37] is essential to damp unphysical highfrequency noise that arise at excision grid boundaries; we use a typical dissipation parameter of KO = 0.35.

Results
As a proof-of-principle, we evolve initial data that undergoes gravitational collapses within one light-crossing time, and follow the subsequent ring-down to the Schwarzschild-AdS solution. The geometry of the initial slice is sourced by a massless real scalar field with a Gaussian 14 The grid on which the AH finder is executed is completely independent of the specifics of the Cartesian evolution grid that was described in Section 5.1. 15 This method is effective in removing singularities if the following common assumptions are valid on the spacetimes that we consider: (i) weak cosmic censorship is not violated, i.e., geometric singularities are contained inside a black hole event horizon; (ii) the AH at any time t is contained in t-constant slices of the event horizon; (iii) the AH at any t provides a sufficiently accurate approximation for t-constant slices of the event horizon.
profile, distorted along each Cartesian direction and centred at x = y = z = 0: The amplitude of the profile is A = 0.55 and the eccentricities are e x = 0.3, e y = 0.2, e z = 0.25, so that the most prominent distortion is on the (x, y)-plane. The width of the Gaussian is ∆ = 0.2. We choose the initial slice to be a moment of time symmetry, and the details of the time-symmetric initial data sourced by this matter field are collected in Appendix D. As we see in that appendix, the momentum constraint is trivially satisfied for this type of data, so only the Hamiltonian constraint has to be solved. We evolve this initial data up to t = 31 in units of the characteristic length scale L = 1 (approximately 20 light crossing times), well after the end of gravitational collapse and the resulting black hole formation. The initial data has zero total angular momentum, and angular momentum conservation [38] ensures that this is zero at all times. Therefore, we can expect the black hole to settle down to the Schwarzschild-AdS solution. However, for generic initial data with non-vanishing total angular momentum, this may not be the final state: Ref. [39] conjectured that Schwarzschild-AdS, or more generally Kerr-AdS, may suffer from a non-linear instability for generic perturbations. We will leave this interesting problem for future work.

Collapse and ringdown
We describe here the evolution in the bulk: this consists of an initial short phase, in which the scalar field collapses and forms a black hole, and a long ringdown stage, in which the spacetime settles down to Schwarzschild-AdS. Figure 2 shows the profile of the scalar field variable,φ, at four representative times on the equatorial plane z = 0 for the highest resolution grid, with N x = N y = N z = 325 grid points along each Cartesian direction. Notice that in all of these snapshotsφ = 0 at the AdS boundary, as required by the Dirichlet boundary conditions. At t = 0, the asymmetry of the initial Gaussian profile is too small to be visible. At the beginning of evolution, we see that the scalar field starts propagating towards the AdS boundary, but a significant portion of it is attracted back towards the origin and collapses to form an AH. This occurs at t = 0.331 in the highest resolution simulation. The rest of the scalar field remains outside the black hole, where it keeps bouncing back and forth the AdS boundary and is gradually absorbed. The asymmetry on the (x, y)-plane is clearly visible at t = 2.6, where the scalar field is stretched along the x-direction and squeezed along the y-direction. The elongation changes its direction multiple times during the evolution, as shown in the next two plots: it is along the y-axis at t = 5.0 and again along the x-axis at t = 7.2. At later times, t 9, the value of the scalar field becomes consistent with zero up to solution error 16 and the spacetime settles down to a Schwarzschild-AdS black hole spacetime with mass M = 0.403. The late-time solution is close to Schwarzschild-AdS, which can be seen explicitly in Figure 3. Here, we compare the numerical solution at the last time slice, i.e., t = 31, to a slice of the Schwarzschild-AdS metric with conserved mass obtained from our highest resolution run (M = 0.403). This comparison is achieved with the following procedure. First, we compute the Riemann cube scalar R 3 = R µνρσ R ρσγδ R γδ µν , and the Kretschmann scalar K = R µνρσ R µνρσ . Second, we compute the corresponding values, R 3 AdS and K AdS , for pure AdS 4 . Figure 2: Snapshots of the scalar field profileφ on the z = 0 slice in (x, y) coordinates. In each plot, x and y are the horizontal and vertical axes, respectively, and the black square denotes the boundary of the numerical grid, i.e., x = ±1 and y = ±1. The external boundary of the coloured part is the AdS boundary. The black ellipse denotes the approximate position of the AH. This is obtained as the z = 0 slice of the ellipsoid with Cartesian semiaxes, x AH , y AH , z AH , where x AH is the x-coordinate value of the intersection between the AH and the x-axis, and similarly for y AH and z AH . The internal boundary of the coloured region is the excision surface: we excise points inside an ellipsoid whose semi-axes, a ex x , a ex y , a ex z , are given by a ex x = x AH (1− δ ex ), and similarly for a ex y and a ex z . We use the value δ ex = 0.4 for the excision buffer. Highest resolution: N x = N y = N z = 325.
We then use all four quantities to represent the relative Riemann scalar (R 3 /R so going to larger values of K K AdS −1 is equivalent to moving towards the centre of the grid, and closer to the singularity. Therefore the black vertical lines give an indication of the position of the AH relative to the AdS boundary. The two panels of Figure 3 indicate that, sufficiently close to the AdS boundary, the curvature invariants of the numerical solution are almost identical to Schwarzschild-AdS. For clarity, this is shown using only values of the Riemann cube and Kretschmann scalars along the x and y axes, but we verified this for values from the entire grid. At any given resolution, the numerical curvature invariants start to differ from their Schwarzschild-AdS values as we get closer to the AH. This is expected since the gradients become larger as we approach the centre of the grid. However, these differences converge away as resolution is increased. Finally, although there is an asymmetry at any given resolution between the x and y axes even at this last time slice, this late-time asymmetry also converges away as resolution is increased.

Boundary scalar field and stress-energy tensor
In this section we consider the evolution of the holographic quantities at the AdS boundary defined in Section 4. These quantities are obtained via third order extrapolation from points in the interior, with the only exception of the t = 0 plot of Figure 4, which is computed analytically from the initial distorted Gaussian profile (6.1). See Appendix F for a detailed explanation of the extrapolation scheme.
We start by noting that the numerical values for the total mass M in AdS, obtained from equation (4.5), are approximately constant during the evolution, as expected by mass conservation [38]. More precisely, a small drift of the total mass is observed numerically, however this becomes smaller as we increase the resolution and it is consistent with zero within our error estimate for boundary quantities that we will discuss shortly. Figure 4 shows four snapshots of the vacuum expectation value of the dual scalar field operator at the boundary,φ (1) , obtained from the near-boundary expansion of the bulk scalar field in (3.5). Unlike the z = 0 slice snapshots of Figure 2, these plots of the boundary S 2 encode the asymmetry in all three Cartesian directions in the bulk, as they appear on the boundary at ρ = x 2 + y 2 + z 2 = 1. In fact, the asymmetry of the initial data is already visible at t = 0, where the different values of eccentricities along the three Cartesian direction (largest along x and smallest along y) are evident in this plot. At this time, the boundary scalar field is overall very small, which is expected since the initialφ, given by (6.1), is localized near ρ = 0. Notice from Figure 4 that the asymmetry changes axes during evolution, but interestingly it is always strongest along x and weakest along y or vice-versa. Furthermore, a direct comparison with Figure 2 shows that the features present at a certain t at the boundary take approximately π/2 1.6 to reach the interior of the bulk, i.e., about a light-crossing time, as expected. At later times, mirroring the evolution in the bulk,φ (1) decays exponentially in time as the bulk spacetime settles down to Schwarzschild-AdS. Figure 5 displays the energy density of the boundary CFT. At t = 0 this is strongly asymmetric along the x-direction, as expected from the shape of the initial scalar field profile (6.1). After that, undergoes a phase of strong evolution with several changes of elongation axes, sampled at t = 5.6 and terminating at approximately t = 7.2. From that time onwards, settles down to a uniform configuration, as appropriate for the Schwarzschild-AdS black hole. Approach to uniformity is emphasized by using colour scales with fixed interval length, centred at the mean value of at the corresponding evolution time.
More information about the energy density of the boundary field theory can be deduced Figure 4: Snapshots of the vacuum expectation value of the dual scalar field operatorφ (1) . The first snapshot is obtained analytically from the initial scalar field profile. The remaining three are obtained by third order extrapolation and subsequent smoothening via a low-pass filter; see Appendix F. Highest resolution: N x = N y = N z = 325.
from Figure 6. The trace trT CF T vanishes for a conformal field theory in 2 +1 dimensions, which is the case for our R × S 2 boundary. In Section 4, we had spelled out how this trace, in our scheme, is tied to how well we are solving the Einstein field equations. We thus use the L 2 -norm of the numerical values of trT CF T (red line) as an error estimate for boundary quantities. We compare this error with the difference between maximum and minimum of (blue line), the L 2 -norm of the difference between and its Schwarzschild-AdS value Schw-AdS = M 4π (green line), with M = M h = 0.403, i.e., the highest resolution value of M , and the L 2 -norm of ∆p (magenta line). We compute these quantities from the data of the highest resolution simulation, but at any resolution the hierarchy is the same, although it appears at different scales. If we exclude very early times, we see that max( ) − min( ) is consistent with zero, which confirms that the energy density becomes uniform in time. We also see that || − Schw-AdS || 2 is consistent with zero and decreasing in time, which shows that the energy density settles down to Schw-AdS , as expected. Finally, ||∆p|| 2 is consistent with zero, as appropriate for the boundary anistropy of the Schwarzschild-AdS black hole.

Discussion
We have presented the first proof-of-principle Cauchy evolution scheme with no symmetry assumptions that solves the Einstein-Klein-Gordon equations for asymptotically AdS spacetimes. Stability of this numerical scheme is achieved through the gauge choice (3.27) near the AdS boundary. We have used this scheme to obtain preliminary results using stationary initial data constructed from completely asymmetric Gaussian initial profiles of a massless scalar field.
We observe the collapse of the scalar field into a black hole and the subsequent ringdown to a Schwarzschild-AdS black hole spacetime, in both bulk and boundary quantities. Deviations from Schwarzschild-AdS at late times are consistent with zero within estimates of the numerical error. At very late times, the spatial profiles of these small deviations appear to cascade towards higher harmonics. Even though these deviations are consistent with our er- ror estimates, they may nevertheless trigger a non-linear instability that can only be revealed by evolving for longer times and with higher spatial resolutions. 17 It will be interesting to conduct a detailed analysis by decomposing the scalar field profile into spherical harmonics and showing that the radial part is non-vanishing near the boundary for a very long time. We leave this for future studies.
In this work we limited ourselves to D = 4 spacetime dimensions, but the calculation outlined in Section 3.2 would be almost identical if we were to study Cartesian evolution of asymptotically AdS spacetimes in any D ≥ 4 dimensions. In particular, the stable gauge found with this method would be the same up to a numerical factor. Interestingly, a comparison between (3.27) and the corresponding result in [28] (see eq. (S10) in that previous work) clearly suggests a trend for the expression of the stable gauge as we relax symmetries, and thus increase the number of spatial coordinates on which the solution depends. If this trend were confirmed, repeating the calculation above would not be necessary when increasing the number of spatial degrees of freedom. See [21] for an example in higher dimensions where this was done explicitly. Furthermore, the scheme presented here can be applied to cases with 17 Schwarzschild-AdS has been shown to be stable under spherically symmetric deformations [39]. different types of matter fields, different types of global coordinates, and to coordinates on the Poincaré patch. For instance, in Appendix B we followed the prescription of Section 3 to obtain the stable gauge also in spherical coordinates. In Appendix C, the same procedure leads to a gauge that stabilizes evolution on a Poincaré patch of AdS 4 . In other words, this framework makes numerical Cauchy evolution in asymptotically AdS spacetimes possible in full generality, with no need to impose symmetries.
We expect to be able to tackle several interesting problems in asymptotically AdS spacetimes using this Cauchy evolution scheme. We want to highlight two of the most important of these here. The first is the study of gravitational collapse in AdS with no symmetry assumptions and with angular momentum. The numerical study of gravitational collapse in AdS was done in D ≥ 4 spacetime dimensions by [40,41] in spherical symmetry. In these papers it was shown that a class of small perturbations of amplitude undergoes gravitational collapse and forms a black hole on a time-scale O( −2 ), due to a turbulent cascade of energies from large to small distances until a horizon forms. Subsequently, [28] considered the same massless scalar field model in AdS 5 in a 2+1 setting, and it was observed that for a certain class of initial data, the subsequent evolution resulted in collapse that happens faster away from spherical symmetry. On the other hand, the authors in [42] used a particular metric ansatz in a 1+1 setting to consider the inclusion of angular momentum, and observed delayed collapse. A promising direction is provided in [43], [44] with a proof of the instability of AdS in spherical symmetry for the Einstein-massless Vlasov system. The scheme described in this article makes it possible for numerical investigations to help settle this question, by incorporating all the relevant physics needed to study gravitational collapse in AdS in full generality.
The second important problem we wish to highlight is the study of the superradiant instability in AdS. Superradiantly unstable (see [45] for a review of superradiance) initial data around a Kerr-AdS black hole spacetime was evolved in [19], without imposing symmetries, up to approximately 290 light-crossing times using the characteristic scheme presented in [20]. This paper showed a transition of the Kerr-AdS black hole to a rotating black hole with one helical Killing field consistent with a black resonator [46]. Since the known black resonators are rapidly rotating black holes with an ergo-region, they are also unstable to superradiance [47,48]. Hence a cascade to smaller and smaller resonators, potentially leading to a violation of the weak cosmic censorship conjecture, was suggested [49]. The authors of [19] see a second transition at late times that could be the beginning of such a cascade but they do not continue the evolution further. Hence, the endpoint of the Kerr-AdS superradiant instability is still unknown. To settle this question, it will be necessary to keep track of progressively smaller spatial scales during the course of the evolution, which typically go hand-in-hand with progressively richer dynamics. On top of the computationally expensive nature of 3+1 simulations, it will be necessary to keep track of the evolution on sufficiently long time scales until the endpoint is reached, requiring commensurately large-scale computational resources. We leave this important study for future work. grant EP/P020232/1, in this research, as part of the HPC Midlands+ consortium. This research also utilised Queen Mary's Apocrita HPC facility, supported by QMUL Research-IT. http://doi.org/10.5281/zenodo.438045.

A Generalized Harmonic Formulation
The generalized harmonic formulation of the Einstein equations is based on coordinates x α that each satisfies a wave equation x α = H α with source functions H α . As long as the constraints C α ≡ H α − x α = 0 are satisfied, we can then write the trace-reversed Einstein equations in D dimensions with cosmological constant Λ, where the choice of H α = g αβ H β fixes the gauge, Γ α βγ are the Christoffel symbols associated with the spacetime metric g αβ , and T αβ is the matter stress-energy tensor. It can be proven that the constraints C α = 0 are satisfied for all times t ≡ x 0 , as long as C α = 0 and ∂ t C α = 0 at t = 0. In an initial-boundary value problem, this still holds if we assume that the boundary conditions are consistent with C α = 0 being satisfied on the boundary for all times. However, for numerical initial data, the constraints and their derivatives with respect to t vanish at t = 0 typically only up to truncation error. Thus, to suppress constraint-violating solutions, we supplement (A.2) with constraint damping terms as introduced in [50], controlled by the parameters κ and P . We thus obtain the final form of our evolution equations: where n α = −∂ α t is the timelike, future-directed unit 1-form normal to slices of constant t. Notice that the principal part of (A.3), − 1 2 g γδ ∂ γ ∂ δ g αβ , is a wave operator acting on metric components. Thus, the well-posedness of the wave equation suggests that the initial-boundary value problem in generalized harmonic form is well-posed, if we make reasonable assumptions on the remaining components of the problem. See, for example, [27,29] for more details on this formulation. In our simulations we use the values κ = −10 and P = −1. 18 In this work, we are interested in the case where matter fields are given by a single massless real scalar field ϕ, hence the stress-energy tensor reads For completeness, we also write the Klein-Gordon equation (2.3) for the scalar field ϕ in terms of partial derivatives with respect to the chosen set of coordinates:

B Boundary Prescription for Spherical Coordinates
Although spherical coordinates x α = (t, ρ, θ, φ) are not suitable for numerically evolving points near the origin (see discussion in Section 5.1), they are convenient to extract the physics of the CFT at the AdS boundary, since they are adapted to the boundary topology R × S 2 . In this section we apply the prescription outlined in Section 3 to the case of asymptotically AdS spacetimes in D = 4 spacetime dimensions in spherical coordinates. Similarly to the Cartesian case, we first define the spherical coordinate version of the evolution variables (ḡ αβ ,φ,H α ). We also write down the transformations between these variables and their Cartesian version, (3.1)-(3.3). Then, we obtain the stable gauge in spherical coordinates by following the steps introduced in Section 3.2. We compare this with a different potentially stable gauge that can be inferred from the one used in [27]. Finally, we show that tracelessness and conservation of the boundary stress-energy tensor T ab CF T , defined in Section 4, is a consequence of the lowest order of the Einstein equations in the near boundary expansion, provided that the leading order of the generalized harmonic constraints is satisfied. 19

B.1 Evolution Variables and Boundary conditions
We remind the reader that new evolution variables are defined in order to apply the boundary conditions found in Section 2.2 as simple Dirichlet conditions at the AdS boundary ρ = 1. In the same way as in the Cartesian coordinate case where we defined metric evolution variables g µν in (3.1), the metric evolution variables in spherical coordinatesḡ αβ are defined by (i) considering the deviation from pure AdS tensor h αβ = g αβ −ĝ αβ in spherical coordinates, and (ii) stripping h αβ of as many factors of (1 − ρ 2 ) as needed so that they fall off linearly in (1 − ρ) near the AdS boundary.
The boundary conditions on h αβ (2.8) tell us that Despite the notation, we emphasize thatḡ αβ andḡ µν are not in general components of the same tensor (as it should be clear from their definition), therefore the usual transformation between tensor components in different sets of coordinates cannot be applied. The correct transformation can be easily deduced from (B.1) and (3.1), remembering that h is indeed a tensor:ḡ Similarly, the boundary conditions on the scalar field (2.9) suggest that we use the evolution variableφ which the same as the one in Cartesian coordinates, as expected for a scalar field. Finally, the boundary conditions (2.13) on H α suggest the use of the evolution variables in spherical coordinates. Neither H α ,Ĥ α ,H α nor H µ ,Ĥ µ ,H µ are components of the same tensor, so there is no simple transformation from one set to the other. The two triplets of quantities can only be obtained from the definition of source functions in terms of the full metric g in the appropriate set of coordinates, e.g., equation (2.13) in spherical coordinates.
In a numerical scheme in spherical coordinates employing the framework presented in this article, reflective Dirichlet boundary conditions can be easily imposed as

B.2 Gauge Choice for Stability
Since the evolution variables in spherical coordinates, (ḡ αβ ,φ,H α ), are linear in q = 1 − ρ by construction, we can borrow the near-boundary expansions (3.5)-(3.7) . We now substitute these into the evolution equations (A.3), and we expand each component in powers of q.
Rewriting the resulting equations in the wave-like form (3.10), we obtaiñ Doing the same for the generalized harmonic constraints 0 = C α ≡ H α − x α , we have We now follow the three steps of Section 3.2 to obtain a stable gauge choice.
1. Solve the leading order of the near-boundary generalized harmonic constraints, (B.16)-(B.19), forH (1)α . We obtain We prove below that these equations ensure tracelessness and conservation of the boundary stress-energy tensor T ab CF T , defined in Section 4.

B.3 Tracelessness and Conservation of Boundary Stress Tensor
We conclude this subsection by showing that tracelessness and conservation of T ab CF T follow from (B.21)-(B.24), i.e., from the lowest order of the Einstein equations, provided that the leading order of the generalized harmonic constraints are satisfied.
With the notation of Section 4, let x a = (t, θ, φ) be the coordinates along the AdS boundary, λ ab dx a dx b = −dt 2 + dθ 2 + sin 2 θdφ 2 be the metric of the AdS boundary, and D be the Levi-Civita connection of λ ab , i.e., D is torsion-free and D a λ bc = 0. Then, trT CF T = λ ab T ab CF T is the trace of the boundary-stress tensor and D a T ab CF T = λ ac D c T ab CF T is its divergence. We want to prove that trT CF T = 0 and D a T ab CF T = 0. The expression of trT CF T in terms of the leading order of the metric variables in spherical coordinates was already written in (4.6). We repeat it here for completeness: The divergence of the boundary stress tensor is given by We immediately see that trT CF T = 0 as a consequence of (B.21). Moreover, by solving the system of 6 equations given by the first derivatives of (B.21) with respect to t, θ, φ and (B.22), (B.23), (B.24) forḡ (1)tt,t ,ḡ (1)tt,θ ,ḡ (1)tt,φ ,ḡ (1)tθ,t ,ḡ (1)tθ,θ ,ḡ (1)tφ,t , and substituting the solution into the right hand side of (B.28)-(B.30), we see that D a T ab CF T = 0.

B.4 Boundary stress tensor from holographic renormalization
We can straightforwardly compute the boundary stress tensor from (2.26) using the holographic renormalization prescription of [33] in spherical coordinates xā = (t,θ,φ) on the AdS boundary. We have ab are the z 3 terms of the metric components in FG form, (2.26). The explicit components of the stress tensor in (B.31) are given by f tφ , On the other hand, in Section 4 we compute the boundary stress-tensor starting from the metric in global spherical coordinates and then using the prescription of [34]. Of course, the expressions (4.4) and (B.32) are equivalent, as we now explain. To obtain (4.4), we have not imposed that the metric components satisfy the Einstein equations. On the other hand, (B.31) gives the correct boundary stress-energy tensor if the bulk metric solves the Einstein equations, in agreement with the assumptions of the FG theorem. It is thus expected that (4.4) and (B.32) agree if we assume the validity of the lowest order of the Einstein equations in the form that takes into account the generalized harmonic constraints, i.e., (B.21)-(B.24). In fact, we only need (B.21). For example, starting from (4.4), imposing (B.21) and using the fact thatt = t,θ = θ,φ = φ at the boundary ρ = 1 together withḡ (1)αβ = f αβ , 20 we find precisely the expressions (B.32).

C Boundary Prescription for the Poincaré Patch
Here we follow the prescription of Section 3 in the case of Poincaré AdS and display a choice of generalized harmonic source functions that stabilizes the evolution in this case.
The metric of the Poincaré patch of AdS 4 can be written aŝ in terms of Poincaré coordinates (t, z, x 1 , x 2 ). To include the Poincaré horizon z → ∞ in our computational domain, we compactify the bulk coordinate z = (1 − ρ 2 )/ρ 2 to have the Poincaré horizon at ρ = 0 and the AdS boundary at ρ = 1. This gives the following form for the metric of AdS 4 :ĝ Let us now consider asymptotically AdS spacetimes. Since (C.1) is in the form given by the leading order of the FG expansion, (2.15)-(2.16), we see that (t, z, x 1 , x 2 ) are FG coordinates. We can thus read off the fall-offs of the metric components from the rest of the FG expansion. The evolved fields consist of the spacetime metric g µν , possibly a scalar field ϕ, and the generalized harmonic source functions H µ . The fall-offs of the metric components g µν read the same as (2.10), with f µν (t, x 1 , x 2 ) coefficients. The scalar field fall-off that preserves the metric asymptotics is given by (2.11), with c(t, x 1 , x 2 ) coefficient. The fall-offs of the source functions can be inferred from the metric fall-offs, which are given by (2.14), with f µ (t, x 1 , x 2 ) coefficients. As a result, the corresponding evolution variables in this Poincaré setting are given exactly by the same expressions as we had written in (3.1)-(3.3).
Using the same steps as in Section 3.2, we obtain the following gauge: We have verified that this gauge leads to stable evolution in asymptotically AdS 4 spacetimes in Poincaré coordinates. We close by noting that [21] obtained a similar stable gauge to evolve dynamical black holes in the background of the AdS soliton.

D Initial Data
The Cauchy problem in GR requires the prescription of initial data on a spacelike hypersurface Σ and a choice of gauge throughout the entire evolution. In an asymptotically AdS spacetime, in addition we have to specify boundary conditions at the boundary of AdS; we have dealt with boundary conditions in Section 3. We pick Cartesian coordinates x µ = (t, x, y, z) such that t = 0 on Σ. The spatial Cartesian coordinates on Σ are denoted by x i = (x, y, z), and the corresponding indices by i, j, k, . . . . With this notation, the data needed for the Cauchy evolution in the generalized harmonic scheme is composed of the initial dataφ| t=0 ,ḡ ij | t=0 , ∂ tφ | t=0 , ∂ tḡij | t=0 and the source functionsH µ at all times. The gauge used in our numerical scheme at t > 0 is discussed in Appendix E. With regard to the gauge at t = 0, we do not setH µ | t=0 explicitly, but we make an equivalent choice forḡ tµ | t=0 , and ∂ tḡtµ | t=0 , and then computeH µ | t=0 from (2.12). In summary, the complete set of initial data that we prescribe isφ| t=0 ,ḡ µν | t=0 , ∂ tφ | t=0 and ∂ tḡµν | t=0 . In this section we explain how this is done in our simulations, taking into account two crucial facts. Firstly, initial data cannot be chosen in a completely arbitrary way, but it must satisfy the constraints of GR. Secondly, the choice of the initial degrees of freedom must be consistent with the desired gauge (3.27) near the AdS boundary.

D.1 Constraints
Here we review the constraints in GR and how they are solved in our numerical scheme. We start by defining the relevant quantities on the initial spacelike hypersurface Σ.
The timelike, future-directed unit 1-form normal to Σ is given by where α = 1/ −g µν (dt) µ (dt) ν is the lapse function. The projection operator onto Σ is defined by (Notice that γ µ ν is idempotent, i.e., γ µ ρ γ ρ ν = γ µ ν , as appropriate for a projector.) This operator can be applied to any tensor at a point p ∈ Σ to obtain the part of that tensor tangent to Σ. For instance, given a vector X at a point p ∈ Σ, X µ || = γ µ ν X ν is the part of X tangent to Σ, i.e., X µ || n µ = 0. Let us now consider a tensor defined on the tangent space of the spacetime manifold M at a point p ∈ Σ. If the tensor is invariant under projection onto Σ, then it can be identified with a tensor defined on the tangent space of Σ at p, under a natural (i.e., basis-independent) isomorphism. For example, γ µν = g µν + n µ n ν at points on Σ can be identified with the Riemannian metric of Σ defined as the pull-back 21 on Σ of the spacetime metric g µν , given by γ ij in spatial Cartesian coordinates. See [35] for more details. Indices of tensors invariant under projection onto Σ can be raised and lowered by γ µν or g µν , equivalently. Indices i, j, k, . . . of tensors on the tangent space of Σ can be raised and lowered by γ ij .
The projection of ∇ µ n ν defines the extrinsic curvature of Σ: 22 The Lie derivative along the normal direction in the second equality suggests that a choice of K µν on Σ is "morally" equivalent to a choice for the time-derivative of the metric components at t = 0. K µν is identified with the tensor on the tangent space of Σ, given by K ij .
As a final ingredient, the covariant derivative on Σ of a tensor field invariant under projection onto Σ is defined as the projection onto Σ of the covariant derivative ∇ of the tensor field, and we denote it by D. For instance, D µ X ν || = γ ρ µ γ ν σ ∇ ρ X σ || and D µ X ν || is identified with the tensor on the tangent space of Σ given by D i X j || = γ ρ i γ j σ ∇ ρ X σ || . D is the Levi-Civita connection of γ ij , i.e., it is torsion-free and D i γ jk = 0.
We can now write the constraints that initial data on Σ must satisfy. The "normalnormal" projection (i.e., contraction with n µ n ν ) of the Einstein equations gives the Hamiltonian constraint where (3) R is the Ricci scalar associated with the connection D, K = γ ij K ij and ρ = T µν n µ n ν is the matter energy density measured by an observer with 4-velocity n µ . The "tangentnormal" projection (i.e., contraction with γ µν n ρ ) of the Einstein equations gives the momentum constraint where j i = −T ρσ n ρ γ σi is the matter momentum density measured by an observer with 4velocity n µ .
We now explain how these constraints are solved for massless real scalar matter, whose energy-momentum tensor is (A.4), in the simplified case of time-symmetric data. Time symmetry in the scalar sector, ∂ tφ t=0 = 0, (D.6) implies j i = 0. Time symmetry in the gravitational sector, together with the initial gauge choiceḡ ti t=0 = 0, (D.8) implies K ij = 0. Thus, we see that the momentum constraint is trivially satisfied. The Hamiltonian constraint, instead, reduces to This can be solved through the conformal approach, initiated in [51], which assumes that the spatial metric γ ij is conformal to the spatial metricγ ij of the t = 0 slice of pure AdS in Cartesian coordinates: where ζ is a smooth positive function on Σ, satisfying the AdS boundary condition ζ| ρ=1 = 1. LetD be the Levi-Civita connection ofγ ij and (3)R the corresponding Ricci scalar. Using (D.10) and its inverse, γ ij = ζ −4γij , we obtain Plugging (D.11) into (D.9) gives Finally, the version of the Hamiltonian constraint that we are going to solve is obtained by writing the matter energy density ρ in terms of ζ. The time-symmetry requirement ∂ t ϕ t=0 = 0 gives ρ = 1 2ζ 4γ ij ∂ i ϕ∂ j ϕ, (D.14) so the Hamiltonian constraint readŝ For any given choice of scalar field ϕ on Σ, (D.15) is an elliptic equation that can be solved for ζ with boundary condition ζ| ρ=1 = 1. In our simulations we pick the initial scalar field profile ϕ| t=0 =φ| t=0 (1 − ρ 2 ) 2 withφ| t=0 specified by (6.1), and we solve (D.15) with a multigrid algorithm, built into the PAMR/AMRD libraries. The initial metric variablesḡ ij | t=0 are then easily reconstructed from (D.10) and γ ij = g ij | t=0 =ĝ ij | t=0 +ḡ ij | t=0 , i.e.,

D.2 Consistency at the boundary
In the previous section we explained how some components of the initial data for our simulations are obtained: (i) we impose time-symmetry, namely ∂ tφ | t=0 = 0 and ∂ tḡij | t=0 = 0; (ii) we make the initial gauge choiceḡ ti | t=0 = 0; (ii) we choose the massless real scalar field profilē ϕ| t=0 given by (6.1); (iii) we determineḡ ij | t=0 through the conformal decomposition of the Hamiltonian constraint. In this section we determine the remaining necessary components for Cauchy evolution based on the generalized harmonic scheme:ḡ tt | t=0 and ∂ tḡtµ | t=0 .
In doing so, the only restriction to consider is the one already obtained in Step 2 of our gauge prescription in Section 3.2: the Einstein equations in a gauge that satisfies the generalized harmonic constraints impose the conditionḡ (1)tt =ḡ (1)xx +ḡ (1)yy +ḡ (1)zz near the boundary. This will hold at all times of the evolution and it must be imposed on initial data. Given that there is no requirement on the value ofḡ tt in the bulk, we make the simplest choice and set that to zero. In order to smoothly transition from the bulk value ofḡ tt to its required boundary value, we use the smooth transition function where R(ρ) = (ρ b − ρ)/(ρ b − ρ a ) and ρ a , ρ b are the values between which the transition takes place, set to ρ a = 0.5, ρ b = 0.9 in our simulations. Thus, our choice ofḡ tt | t=0 is To conclude, the remaining initial variables can be chosen in a completely arbitrary way so we make the simplest choice everywhere on the grid:

E Complete Gauge Choice
In Section 3.2 we discussed the gauge choice of source functions that we impose near the boundary in order to obtain stable evolutions. Furthermore, the gauge at t = 0,H µ | t=0 , is determined from the initial data, detailed in Appendix D, through the definition of source functions (2.12) at t = 0. All that remains is to make a gauge choice ofH µ in the bulk, and smoothly join this with the target boundary values (3.27) on each spatial slice and with the initial valuesH µ | t=0 during evolution. In this section we describe how all this is implemented in our numerical scheme.
We start by choosing a zero value forH µ in the bulk, as this is the simplest choice. Therefore, the values of the source functions on each spatial slice, after the time transition from t = 0, are given by where the spatial transition function f 1 (ρ) is defined as in (D.17) with transition occurring between ρ 1a = 0.05 and ρ 1b = 0.95.
With these ingredients, we can finally write the complete gauge choice made in our simulationsH (E.3) From the properties of g(t, ρ), we see thatH µ =H µ | t=0 at t = 0 andH µ = F µ for t ξ 1 in the interior and t ξ 2 near the boundary. Since the target gauge is crucial for stability and needs to be reached quickly, ξ 2 is typically set to a small value. On the other hand, it is not necessary, and perhaps even troublesome, to deal with a fast transition in the bulk, therefore ξ 1 takes a larger value. In our simulations, we set ρ 0a = 0.0, ρ 0b = 0.95, ρ 1a = 0.05, ρ 1b = 0.95, ξ 1 = 0.1, ξ 2 = 0.0025.

F Boundary Extrapolation
As explained in Section 4, since the AdS boundary generally does not lie on points of the Cartesian grid, we can only obtain the approximated value of any boundary quantity f through extrapolation from the numerical values of f on grid points near the boundary. In this section we describe how extrapolation is implemented in our scheme with the help of Figure 7.
For simplicity, we consider first order extrapolation, i.e., extrapolation from two grid points. The following can be generalized to higher extrapolation orders in a straightforward way. In particular, third order extrapolation is used for the plots in Section 6.2, since this improves the accuracy of the extrapolated numerical values. 23  Figure 7), where n ∆ denotes the degree of the three resolutions used for convergence, n 9h/4 = 0, n 3h/2 = 1, n h = 2 (notice that 3 2 n ∆ ∆ is a constant for all three resolutions). We have empirically found that considering points outside of this region in the next steps leads to unphysical or non-converging values.
2. For any point in the range defined at step 1, identify the coordinate with the largest absolute value, e.g., x, and its sign, say x > 0. If two coordinates have the same absolute value, then we pick x over y and z, and y over z. Each direction identified in this way is represented by a light blue arrow. Among all the points along the identified direction (x in our example) and within the range of step 1, pick the closest point to the boundary. We denote this point by p 1 and its coordinates by (x 1 , y 1 , z 1 ). For each direction identified as above, the corresponding p 1 point is represented as a green dot in Figure 7.
3. Consider the nearest point to p 1 along the identified axis in the direction of the bulk (decreasing x in the example). We denote this point by p 2 and its coordinates by (x 2 , y 2 , z 2 ). For each p 1 point, the corresponding p 2 is represented as a purple dot in Figure 7. In our example x 2 = x 1 − ∆, y 2 = y 1 , z 2 = z 1 .
in cases where the latter are known, e.g. boundary scalar field values at t = 0.
4. Use first order extrapolation on f ∆ (p 1 ), f ∆ (p 2 ) to determine the value of f bdy ∆ (p bdy ) where p bdy is the boundary point along the identified axis in the direction of the boundary. For each pair p 1 , p 2 , the corresponding p bdy is represented by a red dot in Figure 7 and the AdS boundary is represented by a red line. In our example, p bdy is the point with coordinates (x bdy , y bdy , z bdy ) = ( 1 − y 2 1 − z 2 1 , y 1 , z 1 ) and f bdy ∆ (p bdy ) = 5. In order to avoid issues arising from singularities in the definition of spherical coordinates in terms of Cartesian coordinates, we do not extrapolate boundary points with z bdy = 0. Instead, we fill each of these points by copying the mean value of the closest boundary extrapolated points. This ensures continuity at the semi-circle z bdy = 0, y bdy ≥ 0, i.e., points with φ = 0 ∼ 2π. Figure 7 shows that the extrapolated values are not uniformly distributed on the boundary. We aim to improve this in the future by extrapolating the values at points on a uniform (θ, φ) grid with given resolution on the S 2 at the boundary. For now, we fill the empty regions by linearly interpolating boundary values. The data obtained in this way displays high-frequency noise that does not allow for a clear visualisation of physical features. Therefore, we apply a low-pass filter to quantities to be shown on the boundary S 2 . More precisely, we apply the filter on three copies of the boundary sphere joined along the semi-circle z bdy = 0, y bdy ≥ 0 and then we plot the smooth data of the central copy. After re-enforcing continuity at the semi-circle as explained in step 5 above, this strategy provides regular smooth data at the semi-circle if the original raw data is approximately periodic in φ with period 2π, which is expected for data on a sphere.
Notice that, as (F.1) shows, second order convergence of boundary values f bdy ∆ is a direct consequence of second order bulk convergence of f ∆ , which is confirmed by Figure 8 in our simulations. Despite this fact, some modifications must be made to our extrapolation scheme if we wish to perform explicit convergence tests on our boundary data. We now explain the reason for this and the necessary modifications. We assume the validity of the Richardson expansion [52] for f ∆ at any grid point p, where f (p) is the true value of f at p and the rest of the right hand side is the solution error of f ∆ (p). The validity of this expansion is confirmed by bulk convergence of f ∆ to f . Then, from (F.1), we obtain the Richardson expansion for f bdy ∆ at any extrapolated boundary point p bdy : f bdy ∆ (p bdy ) = f (p bdy ) + e extr (p bdy , p 1 , p 2 ) + e ∆ (p 1 , where the f (p bdy ) is the true value of f at p bdy , e extr (p bdy , p 1 , p 2 ) is the error due to the extrapolation approximation. The remaining error terms come from the solution error in f ∆ . The typical convergence test involves the computation of the convergence factor Q(p bdy ) = 1 ln(3/2) ln f 9h/4 (p bdy ) − f 3h/2 (p bdy ) f 3h/2 (p bdy ) − f h (p bdy ) (F.4) at each boundary point p bdy . We clearly see that Q(p bdy ) can be expected to asymptote to 2 as ∆ → 0, thus confirming second order convergence in the continuum limit, only if the points p 1 , p 2 are the same for all 3 resolutions involved. Therefore, our extrapolation scheme must be modified to select pair of bulk points, p 1 and p 2 , for extrapolation that are present in all three grids involved in the convergence test. In practice, we saw that boundary convergence follows the trend of bulk convergence only if, in addition to this modification, we restrict to p 1 points in the range mentioned in step 1 above. The reason for this should be investigated further.
Finally, (F.3) shows that this type of test does not prove convergence to the true value f (p bdy ), but rather to its approximation f (p bdy ) + e extr (p bdy , p 1 , p 2 ). For this reason, the convergence test (G.1) cannot be performed at the boundary for functions with vanishing true value (such as trT CF T ), because their extrapolated value is not just the term linear in ∆ 2 but it also includes the extrapolation error c extr . A more detailed analysis must be made to examine the explicit form e extr (p bdy , p 1 , p 2 ) and be able to find the rate of convergence to f (p bdy ). In our study, we simply make the natural assumption that e extr (p bdy , p 1 , p 2 ) decreases as we increase resolution, so f bdy ∆ (p bdy ) is a sufficiently accurate approximation of f (p bdy ) for sufficiently high resolution (i.e., sufficiently small ∆).

G Convergence of the Independent Residual
To show that the solution is converging to a solution of the Einstein equations, we compute the independent residual that is obtained by taking the numerical solution, substituting it back into a discretized version of the Einstein equations. At each grid point, we then take the maximum value over all components of the Einstein equations, which we denote by Φ ∆ . The independent residual should be purely numerical truncation error, so we can compute a convergence factor for it by using only two resolutions: Q EF E (t, x, y, z) = 1 ln(3/2) ln Φ 3h/2 (t, x, y, z) Φ h (t, x, y, z) . (G.1) Again, with second-order accurate finite difference stencils and with a factor of 3/2 between successive resolutions, we expect Q to approach Q = 2 as ∆ → 0. Figure 8 displays the L 2 -norm of the convergence factor (G.1) for two pairs of resolutions. It clearly shows second order convergence to a solution of the Einstein equations, after an initial transition phase.