Numerical solution of gravitational dynamics in asymptotically anti-de Sitter spacetimes

A variety of gravitational dynamics problems in asymptotically anti-de Sitter (AdS) spacetime are amenable to efficient numerical solution using a common approach involving a null slicing of spacetime based on infalling geodesics, convenient exploitation of the residual diffeomorphism freedom, and use of spectral methods for discretizing and solving the resulting differential equations. Relevant issues and choices leading to this approach are discussed in detail. Three examples, motivated by applications to non-equilibrium dynamics in strongly coupled gauge theories, are discussed as instructive test cases. These are gravitational descriptions of homogeneous isotropization, collisions of planar shocks, and turbulent fluid flows in two spatial dimensions.


Introduction
The advent of gauge-gravity duality (or "holography") has revolutionized the study of strongly interacting field theories. 1 In gauge theories with gravitational duals, holography maps problems involving non-equilibrium quantum dynamics, in the limit of strong coupling and large gauge group rank, into problems involving classical gravitational dynamics in higher dimensions. Consequently, holography provides unique opportunities to study strongly-coupled non-equilibrium dynamics -provided one can actually solve the associated classical gravitational dynamics. Gauge theories with known gravitational duals are generally conformal field theories (CFTs) or conformal theories deformed by relevant operators; for such theories the dual gravitational description involves a 5D spacetime with asymptotically anti-de Sitter (AdS) geometry. 2 Much work to date has explored near-equilibrium phenomena via holography. Examples include the study of viscosity and related transport coefficients [6][7][8][9], more general aspects of dissipative hydrodynamics [10][11][12], quasi-normal modes and near-equilibrium response [13,14], and dynamics of probes such as heavy quarks [15][16][17][18][19] or light quarks [20,21]. 3 There has been much less work on far-from-equilibrium dynamics, as this requires the solution of gravitational dynamics with non-trivial initial conditions and (except in extremely special cases) such solutions can only be found numerically. 4 Despite the difficulty, notable progress has been made on gravitational initial value problems involving asymptotically AdS geometries. 5 Recent work applying holography to far-from-equilibrium dynamics includes studies of isotropization in spatially homogeneous systems [24][25][26], boost-invariant expansion with transverse homogeneity [27,28] or with radial flow [34], spherically symmetric thermalization [29], collisions of planar shocks [30][31][32], and turbulence in 2D fluids [33]. With the exception of the last example, all these problems have a sufficient degree of symmetry that the 5D Einstein equations reduce to either 1+1 or 2+1 dimensional partial differential equations (PDEs). An obvious goal of current and future work is the solution of initial value problems involving lower degrees of symmetry.
In this paper, we discuss the computational challenges involved in solving, numerically, asymptotically anti-de Sitter gravitational initial value problems. We describe in detail a particular approach which we have found to be effective in a series of progressively more complex applications, three of which will be examined as specific test cases: homogeneous 1 See, for example, refs. [1][2][3][4] and references therein. 2 The Klebanov-Strassler cascading gauge theory [5] is an example of a more complicated theory whose dual geometry is not asymptotically anti-de Sitter. 3 For additional prior work in this area, see the recent review [22], which is focused on applications to heavy ion collisions, and references therein. 4 For a broad perspective on numerical relativity and applications to high energy physics, see ref. [23]. 5 In the numerical relativity community, the phrase "initial value problem" is sometimes viewed as referring, specifically, to dynamical evolution problems based on spacelike Cauchy surfaces. A "characteristic" formulation refers to evolution schemes based on null surfaces. We will treat the phrase "initial value problem" as encompassing dynamical evolution schemes with initial data given on either spacelike or null surfaces.
isotropization, planar shock wave collisions, and turbulent 2D fluid flows. Results on homogeneous isotropization have been previously reported in ref. [24]; the degree of symmetry for this problem reduces the 5D Einstein equations to a set of coupled 1+1 dimensional PDEs. Initial results for planar shock wave collisions were reported in ref. [30]; in this case the 5D Einstein equations reduce to 2+1 dimensional PDEs. Studies of fluid flows in two spatial dimensions, using holography, involve the solution of 4D Einstein equations [33]. With no simplifying symmetry restrictions, this case requires solving coupled 3+1 dimensional PDEs. We present results for each of these test cases and discuss both the associated physics and computational issues such as stability and accuracy. The results presented in this paper extend and complement earlier work. In particular, for colliding planar shock waves, we show that it is possible to perform numerically stable, accurate, calculations without adding any artificial background energy density, as was done in refs. [30][31][32]. We study stocks of different thicknesses, as in ref. [31], but integrate farther in time. In agreement with ref. [31], we find that collisions of relatively thin shocks are not well approximated as boost invariant. However, we show that the resulting hydrodynamic flow may be characterized as locally boost invariant, in a sense which we discuss in section 4.2.
For simplicity of presentation, most discussion in this paper is limited to problems involving solutions to pure Einstein gravity which are asymptotic to the Poincaré patch of anti-de Sitter space [1]. Many interesting extensions are only touched upon or left to future work. These include generalizations of these methods to problems involving non-flat boundary geometry (e.g., global AdS asymptotics, or explicit time-dependent boundary geometries [24,27]), additional compact dimensions (e.g., dynamics of initial states in N = 4 super-Yang-Mills (SYM) theory which are not invariant under the full SU (4) R-symmetry), or additional dynamical fields (dilaton-gravity, Maxwell-Einstein, etc.). 6

Setup and conventions
Gauge/gravity duality relates certain quantum field theories in D spacetime dimensions to gravitational physics in D+1 dimensions. (As noted above, we are not considering problems in which the dynamics of additional compact dimensions play any role in the gravitational description.) We consider quantum field theories (QFTs) in D dimensional flat Minkowski spacetime, and hence will be interested in gravitational solutions describing D+1 dimensional geometries with boundary, for which the boundary geometry is D dimensional Minkowski space. Using Fefferman-Graham coordinates, the resulting asymptotically AdS metric may 6 As 5D Einstein gravity is a consistent truncation of 10D IIB supergravity on AdS5 × S 5 , all the 5D pure Einstein gravity solutions we discuss may be viewed as supergravity solutions which describe the dynamics of N = 4 SYM states invariant under the SU (4)R symmetry.
be represented as [35-38] 7 ds 2 = L 2 ρ 2 g µν (x, ρ) dx µ dx ν + dρ 2 , where ρ is a "bulk" radial coordinate such that the spacetime boundary lies at ρ = 0, with {x µ } denoting the D remaining "boundary" coordinates. We use L as the spacetime curvature scale; it is related to the cosmological constant via The metric functions g µν (x, ρ) have a near-boundary asymptotic expansion in integer powers of ρ, with the leading term equal to the desired Minkowski boundary metric and subleading terms starting at order ρ D , 8 3) It will be convenient to use a rescaled stress-energy tensor 4) where G N is the D+1 dimensional Newton gravitational constant. [For D = 4, Newton's constant is related to the dual SU (N c ) SYM theory via G N = π 2 L 3 /N 2 c , so κ = N 2 c /(2π 2 ).] The coefficient of the first sub-leading term in the near-boundary expansion (2.3) determines the boundary stress-energy tensor, which coincides with the expectation value of the (rescaled) stress-energy tensor in the dual QFT, µν (x) . (2.5) Einstein's equations imply boundary stress-energy conservation and tracelessness [35,37,38], Given a non-vanishing stress-energy expectation value, one may define an associated velocity field u(x) and (rescaled) proper energy density ε(x) as the timelike eigenvector and corresponding eigenvalue of the stress-energy tensor,
(with normalization u(x) 2 ≡ −1), provided T µν (or − T µν ) satisfies the weak energy condition. 9 An observer moving with spacetime velocity u(x) sees an energy density equal to ε(x) and vanishing energy flux. For later convenience, let u = u µ dx µ denote the one-form dual to the vector field u = u µ ∂ µ .

Computational strategy
A basic issue affecting any numerical relativity calculation is the choice of how to deal with the diffeomorphism invariance of general relativity. This lies at the heart of how one converts Einstein's equations into a well-posed initial value problem. One general approach is to choose an ansatz for the metric, whose form will greatly restrict the remaining diffeomorphism freedom. The Fefferman-Graham form (2.1) is one such possibility. The ansatz must allow an arbitrary metric, consistent with the symmetries of the physical problem under consideration, to be transformed into the chosen form by a suitable change of coordinates. Even when this is possible in any local region, a given ansatz for the metric may fail to provide good coordinates covering the entire domain of interest. This is a known problem with the Fefferman-Graham form (2.1). Although convenient and useful for analyzing near-boundary behavior, in solutions describing gravitational infall and horizon equilibration, the Fefferman-Graham metric develops coordinate singularities in the bulk and fails to remain regular across the future event horizon [39,40]. Consequently, despite its utility for other purposes, the Fefferman-Graham ansatz is not a good choice for numerical initial value problems.
A different approach, avoiding the need to commit to some specific global form of coordinates, is provided by the ADM formalism in which Cauchy surfaces are arbitrary spacelike slices of the geometry, and some chosen lapse function and shift vector field relate the coordinates on neighboring spacelike slices foliating the geometry [41,42]. This approach has been widely used in numerical relativity calculations in asymptotically Minkowski space [43][44][45]. However, employing this approach has some practical downsides. Implementing this method (particularly when combined with adaptive mesh refinement) is complex. One must formulate a scheme for dynamically choosing lapse and shift vectors, or make some a-priori choice, in a manner which, one hopes, will allow the foliation to remain regular throughout the spacetime region of interest. Achieving a numerically stable scheme can be problematic [43,[46][47][48][49].

Metric ansatz
We have chosen to employ the first approach involving a metric ansatz, one which is specifically tailored to gravitational infall problems. The metric ansatz is a generalization of tra-ditional ingoing Eddington-Finkelstein coordinates for black holes. 10 It is based on a null slicing of spacetime constructed from infalling null geodesics, and will lead to a characteristic formulation of gravitational dynamics. 11 The general form of the metric is where r is a non-inverted bulk radial coordinate (so the spacetime boundary lies at r = ∞), and {x µ } denote the D remaining boundary coordinates. 12 The boundary one-form w = w µ dx µ appearing in the second term is independent of the radial coordinate r. This one-form is assumed to be timelike and, without loss of generality, may be taken to satisfy w 2 = −1 (using the boundary metric discussed below). A more explicit representation of the metric g µν which describes the geometry on fixed-r slices will be introduced in section 3.4. From the ansatz (3.1), one immediately sees that lines along which r varies while the other coordinates are held fixed are null curves. One may easily check that these curves are infalling null geodesics for which r is an affine parameter. Therefore, the vector ∂ r is a directional derivative along infalling null geodesics. At the boundary (r = ∞), an observer whose Dvelocity components equal w µ would describe these geodesics as representing trajectories of comoving objects at rest in his frame; their tangent vectors are normal to the D−1 spatial basis vectors in the observer's frame. In our coordinates, these geodesics remain purely radial throughout the bulk geometry.
The form of the metric ansatz (3.1) is preserved by two types of residual diffeomorphisms: arbitrary D-dimensional diffeomorphisms (independent of r), and arbitrary shifts in the radial coordinate (depending on x), The diffeomorphism freedom (3.2) may be used to transform the boundary one-form w into a standard form such as This simple choice will be used in the examples presented in section 4. Alternatively, one could choose to require that the boundary one-form w coincide with the flow field u which will (eventually) be extracted from the boundary stress-energy tensor via eq. (2.7). Circumstances in which this may be desirable will be discussed in the next subsection. 10 Previous work [11,40] studying late time behavior of solutions which approach stationary black brane solutions convincingly demonstrates the virtues of using generalized Eddington-Finkelstein coordinates for this class of asymptotically AdS gravitational infall problems. 11 For useful prior discussions of characteristic formulations of relativity, see ref. [50] and references therein. 12 The inverse metric G M N = (L/r) 2 (g µν −w µ w ν /w 2 ) −w µ /w 2 −w ν /w 2 −(r/L) 2 /w 2 , with w µ (x, r) ≡ g µν (x, r) wν (x).

Boundary metric and asymptotic behavior
We are interested in solutions to Einstein's equations for which the boundary geometry is flat Minkowski space. Using the ansatz (3.1), such solutions may be expanded, asymptotically, in inverse powers of r, The leading term h µν (x) is the D dimensional boundary metric. This equals the r → ∞ limit of the induced metric obtained by restricting the D+1 dimensional metric (3.1) to r = const. slices, ds 2 r=const. = (r 2 /L 2 ) g µν (x, r) dx µ dx ν , after rescaling to remove the overall r 2 /L 2 factor. The order-D coefficient g µν in expansion (3.5) cannot be determined solely by a near-boundary analysis; the value of this coefficient (which depends on the solution throughout the bulk) determines the boundary stress-energy tensor in a manner similar to the Fefferman-Graham case. With w α ≡ h αβ w β , one finds The boundary metric h µν may be chosen to equal the standard Minkowski metric, But demanding a flat boundary geometry does not obligate one to use Cartesian Minkowski space coordinates. Use of the boundary metric (3.7) represents a further, arbitrary choice of coordinates on the boundary geometry. Alternatively, one may choose to describe Minkowski space using some set of coordinates {x µ } which are non-trivially (and non-linearly) related to a set of Cartesian Minkowski coordinates {y α }, so that For some problems, the standard choice (3.7) of Minkowski boundary metric is sufficient. This will be the case for the specific examples presented in subsequent sections. For other problems, exploiting the freedom of using non-Cartesian boundary coordinates, with corresponding boundary metric (3.8), is helpful. This is true, for example, in problems involving cylindrical or spherical symmetry in the dual field theory, where it is natural to use boundary coordinates adapted to that symmetry.
We believe that exploiting the freedom to choose non-Cartesian boundary coordinates will also be helpful in problems involving highly relativistic fluid flow with large gradients of flow velocity. For such problems, it will undoubtedly be preferable to choose the congruence of radial null geodesics underlying the ansatz (3.1) to involve geodesics describing infalling matter which is at rest (or nearly at rest) in the local fluid rest frame -not at rest with respect to some globally defined inertial Lorentz frame which is necessarily divorced from any local physics of interest. This implies that one would like to use the D-dimensional diffeomorphism freedom (3.2) to set the boundary one-form w appearing in the ansatz (3.1) equal to the flow field u, as suggested earlier. 13 However, simplifying features in the form of the resulting Einstein equations (discussed next) are easier to exploit, numerically, if one uses the diffeomorphism freedom (3.2) to transform the choice to a formally equivalent description where the boundary one-form w has the standard form (3.4) and the complexity of the actual fluid flow is isolated in non-trivial boundary coordinates, for some choice of y α (x). This amounts to changing from an Eulerian to a Lagrangian description of fluid mechanics. The desired diffeomorphism is one for which y α (x), for fixed values of the spatial coordinates x ≡ {x i } and varying x 0 , gives the worldline of a fluid cell labeled by x. If one chooses x 0 to coincide with proper time along this worldline, then the required diffeomorphism is one for which

Horizons and IR cutoffs
The radial direction in AdS is related, via holography, to the energy scale in the dual conformal field theory. Dynamics arbitrarily deep in the bulk correspond to arbitrarily low energy processes in the quantum field theory. With bounded resources, any numerical calculation can only be accurate over a finite dynamic range. So it is inevitable that some form of high energy (UV) and low energy (IR) cutoff will be necessary in any numerical calculation. An effective UV cutoff is imposed by the discretization used when solving differential equations; this will be discussed below. Here, we focus on IR issues. The Poincaré horizon of AdS D+1 is the locus of events beyond which no signal can reach any boundary observer. Any infalling perturbation will distort the geometry and hence perturb the Poincaré horizon. A perturbation with uniform (boundary) energy density can deform the geometry to AdS-Schwarzschild (AdS-BH) form [51], describing a black brane embedded in asymptotically anti-de Sitter space. This geometry has a non-compact planar event horizon, with an associated temperature which is related to the radial position of the horizon. In more general cases of gravitational infall in asymptotically AdS space, one should expect a time-dependent geometry which, at least at late times, will resemble the AdS-BH solution in a local "tubewise" fashion [11].
The essential point is that a non-compact event horizon, with the topology of a plane, may be regarded as an effective IR cutoff. From a holographic perspective, the energy scale of this IR cutoff is set by the local temperature of the horizon. Events beyond this horizon cannot affect any physics extracted by a boundary observer.
In a numerical calculation of the evolving geometry, one is free to excise the portion of spacetime beyond such an event horizon. However, the location of the event horizon cannot be determined without knowing the entire future spacetime geometry (because gravitational infall arbitrarily far in the future can change which null congruence is picked out as the event horizon). Of more practical utility is the identification of an apparent horizon, or outermost marginally trapped surface which, if it exists, will lie inside the true event horizon. 14 We will require initial data such that, at some initial time t 0 , there exists an apparent horizon at some radial position r = r h (t 0 , x). And we will require that this horizon smoothly evolve into an apparent horizon located at radial position r = r h (t, x) on subsequent time slices. Hence, we are assuming that there exists an apparent horizon which, on every time slice, has a planar topology and whose radial position is a smooth function of x and t. The location of this apparent horizon will function as an IR cutoff and will be the boundary of our computational domain. The modification of initial data needed to create or adjust the location of such an apparent horizon is simple: it corresponds, in the dual field theory, to adding a small background energy density. Explicit examples will be discussed in the context of our test cases below.
We will find that some of the fields in our metric ansatz grow, in a power-law fashion, as one moves deeper into the bulk. This can lead to increasingly large problems with numerical loss of precision, which will be discussed in more detail below. Such precision loss can be ameliorated by increasing the IR cutoff, or in other words, choosing initial data which leads to larger values of the apparent horizon radius.
The bottom line is that excising the geometry inside the horizon is not only allowable, it is necessary to avoid numerical problems. The location of the apparent horizon may be tuned by suitably adjusting what, in the dual field theory, is a small background energy density. This is an appropriate point at which to discuss the limits of applicability of our methods. We require that the metric ansatz (3.14) provide good coordinates throughout the region of spacetime between the boundary and an apparent horizon at some radial position r = r h (t, x). This could potentially fail if: (a) some coordinate singularity develops in the spacetime region outside the apparent horizon, or (b) an apparent horizon of the assumed form does not exist.
Since our coordinates are directly tied to the congruence of infalling radial null geodesics, possibility (a) would mean some event is not uniquely identified by our coordinates (x, r), which label a particular infalling radial geodesic (originating at point x on the boundary), together with an affine position r along this geodesic. This is precisely what happens when there is focusing of the geodesic congruence, leading to intersections between differing geodesics. The boundary of the region where such intersections occur defines a caustic. As illustrated schematically in fig. 1, a localized perturbation will typically lead to geodesic focusing and r x r x Figure 1. Focusing of null infalling radial geodesics and consequent formation of caustics. Only the radial direction and one spatial direction are shown. The grey shaded "blob" represents some perturbation in the geometry causing focusing of infalling geodesics. The shaded area at the bottom of each figure represents events behind the apparent horizon. Left panel: caustic formation outside the apparent horizon. Right panel: caustic hidden behind apparent horizon. consequent formation of caustics. Our method assumes that any such caustics lie outside the computational domain; in other words, they must be hidden behind the apparent horizon.
Possibility (b), or non-existence of a planar topology apparent horizon, can occur if the apparent horizon changes form discontinuously. For example, gravitational infall could lead to the formation of a compact trapped surface which is disconnected from a non-compact apparent horizon lying deeper in the bulk. This is illustrated schematically in fig. 2. Of course, the formation of such a compact apparent horizon will likely also lead to focusing and caustic formation in nearby infalling geodesics, so these two failure modes are interrelated.
In either case, the applicability of our methods should be restored if the value of the IR cutoff is increased, i.e., if the position of the non-compact planar topology horizon is pushed outward by increasing the size of the background energy density in the dual theory, as illustrated in the right panels of figs. 1 and 2. Consequently, for some problems, one should expect there to be a limit on the maximum scale separation achievable between the IR cutoff and the physics of interest.
Despite this limitation, we find that a large range of interesting problems are amenable to solution using our methods. In fact, we have yet to encounter difficulties with caustic formation or horizon topology change in any gravitational infall problem we have studied. The underlying issue is one of relative scales. The above described pathologies are likely to occur if one is studying situations with variations in the geometry (or bulk sources) which are spatially localized on a scale which is small compared to the gravitational infall time associated with the apparent horizon. In the dual field theory, this corresponds to states with spatial structure on scales which are small compared to the length scale (or inverse temperature) τ x r x r Figure 2. Possible forms of apparent horizon evolution induced by gravitational infall. Only the radial and one spatial direction are shown. The solid, dashed, and dotted lines, bounding progressively lighter shaded regions, show the position of the apparent horizon at three times t 0 , t 1 , and t 2 , respectively, with t 0 < t 1 < t 2 . Right panel: planar horizon topology at all times, to which our methods apply. Left panel: non-planar horizon topology (at times t 1 and t 2 ), requiring different computational methods. set by the energy density. In a strongly coupled theory, one expects fine spatial structure on scales small compared to τ to be washed out on the microscopic time scale τ , with negligible influence on the later evolution. So, in practice, sources or initial conditions of most interest are those with spatial size large compared to τ . Caustics generated by such sources should be hidden by an apparent horizon whose infall time (or inverse local temperature) is set by the microscopic scale τ .

Einstein's equations
Turning Einstein's equations into a computable time-evolution scheme necessitates a separation of these equations into those which specify the dynamical evolution of the geometry, and those which impose constraints on acceptable initial data (or boundary data). To exhibit key aspects of the explicit equations which emerge when our metric ansatz (3.1) is inserted into Einstein's equations, we use the diffeomorphism invariance (3.2) of the ansatz to specialize the boundary one-form w to the standard choice (3.4), and rename the metric components in the ansatz, (3.14) 15 We have redefined A by a factor of two, and flipped the sign of Fi, relative to definitions in our earlier works [24,27,30]. This change simplifies the forms (3.17)- (3.19) for the radial shift transformations and associated radially covariant derivatives presented below. Our metric ansatz (3.14) is closely related to the Here, X = (x, r) ≡ (t, x, r) denotes event coordinates in which t ≡ x 0 is a null time coordinate, r remains the AdS radial coordinate, and x ≡ {x i } denotes the remaining D−1 spatial coordinates. 16 For later convenience, let denote the spatial dimensionality of the boundary theory. Spatial (ν dimensional) diffeomorphisms are a residual invariance of the form (3.14) of the metric, and transform the metric functions G ij , F i , and A in the usual manner (as components of a covariant tensor, one-form, and scalar field, respectively). As mentioned earlier, arbitrary radial shifts, r →r = r + δλ(x) , (3.16) also leave the form of the metric invariant. Metric functions transform as . From the transformations (3.17), it is apparent that A and F i function as the temporal and spatial components of a "radial shift" gauge field. In light of the spatial diffeomorphism invariance of the metric ansatz, it must be possible to write explicit forms of the resulting Einstein equations in a manner which is manifestly covariant under spatial diffeomorphisms. In addition, it is possible, and quite helpful, to write expressions in a form which also makes invariance under the radial shift symmetry manifest. To do so, we introduce derivatives which transform covariantly under both radial shifts and spatial diffeomorphisms. For the temporal derivative, this is accomplished by defining (3.18) As noted earlier, ∂ r is a directional derivative along ingoing radial null geodesics. The d + derivative is the corresponding directional derivative along the outgoing null geodesic which passes through some event X in the radial direction. The analogous definition for spatial derivatives, acting on (spatial) scalar functions, is Geometrically, these are derivatives along spacelike directions which are orthogonal (at the event X) to the plane spanned by tangents to ingoing and outgoing radial null geodesics.
null Bondi-Sachs form (see, for example, ref. [50]). The key difference is that our radial coordinate r is an affine parameter along infalling geodesics, whereas the Bondi-Sachs metric uses a non-affine radial coordinatē r, chosen to make the determinant of the spatial metric a prescribed function ofr. 16 Using the symbol v instead of t for the null time coordinate would be traditional, as this is customary in discussions of black hole geometries using Eddington-Finklestein (or Kruskal) coordinates. We choose not to do so, but readers should keep in mind that t = const. surfaces are null, not spacelike. Near the boundary our coordinates (t, r) are related to Fefferman-Graham coordinates (x 0 FG , ρ FG ) via r = 1/ρ FG and t = x 0 FG − ρ FG .
In the derivatives (3.18) and (3.19), A and F i act like gauge field components, with ∂ r the associated "charge" operator. When acting on spatial tensor fields, one must augment the derivative (3.19) with an affine connection, which we denote by Γ i jk , to build a derivative which is also covariant under spatial diffeomorphisms. The required connection is the usual Christoffel connection associated with the spatial metric G ij except that, to maintain radial shift invariance, the spatial derivatives appearing in the definition of the connection must be replaced by d i derivatives. Hence, Here, and henceforth, we use primes to denote radial differentiation. We denote by ∇ the resulting spatial and radially covariant derivative. When displaying indices, we use a vertical bar (|), instead of the usual semicolon, to indicate this modified covariant derivative. So, for example, if v is a spatial vector and ω a spatial one-form, then The modified covariant derivative is both metric compatible, G ij|k = 0, and torsion free, Γ i jk = Γ i kj . Associated with our modified spatial covariant derivative is a modified spatial Riemann curvature tensor, R ijkl , defined by the usual formula, but with our modified derivatives replacing the usual derivatives. 17 With these preliminaries in hand, we now examine the resulting Einstein equations. The D+1 dimensional set of equations (3.12) must decompose into one symmetric rank two spatial tensor equation, two spatial vector equations, and three spatial scalar equations. After tedious work, one finds the following simple results. The three scalar equations may be written as: 18 The modified spatial Ricci tensor and scalar are given by the usual contractions, R jk ≡ R i jik and R ≡ R k k . The modified Riemann tensor is antisymmetric in the last two indices, as usual, but need not be antisymmetric in the first two indices, or symmetric under The two-form Ω, defined in eq. (3.25), is the "magnetic" field strength associated with the radial shift symmetry. The extra piece ∆ R ijkl of the modified spatial Riemann tensor leads to a corresponding term ∆ Rij = 1 4 [G · Ω − Ω · G + Ω (tr G )] in the modified spatial Ricci tensor which is antisymmetric. 18 We use a mixture of index-free notation (for simple factors like tr G , F · F , or ∇ · F ), together with indices on more involved expressions; this makes the results most concise. Factors of the inverse spatial metric G −1 = G ij are implicitly present in raised spatial indices. Be aware that raising of indices does not commute with radial or temporal differentiation.
The dot products appearing here and in subsequent expressions are defined using the spatial metric G ij . The spatial tensors G and F are defined as the radial derivatives of covariant components, (G ) ij ≡ (G ij ) and (F ) i ≡ (F i ) . Likewise for G , d + G, d + F , etc. Hence, In equation (3.24), the last term involves the square of the two-form which is the spatial (or "magnetic") part of the field strength associated with the radial shift symmetry. The penultimate term involves the corresponding time-space (or "electric") part of the radial shift field strength, The two vector equations are: And the symmetric tensor equation is: with R ij the modified spatial Ricci tensor. The trace of this equation separates from the traceless part, and reads with R the modified spatial Ricci scalar. Every term appearing in eqs.

Propagating fields, auxiliary fields, and constraints
To elucidate the structure of equations (3.22)-(3.30) it is helpful to write them in a more schematic form after extracting an overall scale factor Σ from the spatial metric G ij . Let with the rescaled metricĝ ≡ ĝ ij defined to have unit determinant, 19 detĝ(X) = 1 .
(3.32) 19 The spatial scale factor Σ must be non-zero throughout the computational domain, as any zero in Σ implies a coordinate singularity at which the metric degenerates. The determinant of the spatial metric (3.31) coincides (up to a sign) with the determinant of the complete bulk metric (3.14), det Gij = − det gMN = Σ 2ν .
Equations (3.22), (3.27), and (3.23) are linear second order radial ordinary differential equations (ODEs) for Σ, F , and A, respectively, having the forms 20 The trace (3.30) and traceless parts of the tensor equation (3.29), and the second vector equation (3.28), are first order radial ODEs for the (modified) time derivatives d + Σ, d +ĝ , and d + F , respectively, with the schematic forms The final scalar equation (3.24) directly expresses the (modified) second time derivative of Σ in terms of the fieldsĝ, Σ, F , and A, plus the first d + derivatives of Σ andĝ, The coefficient functions appearing in the above linear operators are appearing in the inhomogeneous 20 To convert eq. (3.22) to the form (3.33), note that detĝ = 1 implies that tr (ĝ ) = 0 and tr (ĝ ) = tr (ĝ 2 ).
ODEs (3.34)-(3.39) depend only on the indicated fields (and their radial and spatial derivatives). Explicit forms of these source functions may be easily extracted from eqs. (3.27), (3.30), (3.29), (3.23), (3.28) and (3.24), respectively. The function A is an auxiliary field; no time derivative of A appears in any of the above equations. One must integrate the second order radial ODE (3.35) on every time slice (after determining the fields appearing in the source term for this equation) to find A. 21 The first order radial ODEs (3.36), (3.37) and (3.38) determine the modified time derivatives of Σ,ĝ, and F . One may regard the functions Σ,ĝ, and F as propagating fields, with the second order ODEs (3.33) and (3.34) serving as constraints on initial data for Σ and F . If these constraints hold at one time, then the dynamical equations (3.36)-(3.38) ensure that these constraints will be satisfied at all later times.
Alternatively, one may choose to regard Σ and F as auxiliary fields which are determined on each time slice by integrating the second order ODEs

Residual gauge fixing
The residual reparameterization freedom associated with radial shifts (3.16) is apparent in the asymptotic near-boundary behavior of solutions to Einstein's equations. After choosing the boundary metric (3.7) one finds the asymptotic behavior: 22 where λ = λ(x) is completely undetermined. Here and henceforth we have, for convenience, set the curvature scale L = 1. (Factors of L can be restored using dimensional analysis.) As mentioned earlier, asymptotic analysis also cannot determine the values of the subleading order-D coefficients in the metric which, after the renaming (3.13) of metric functions, are the coefficients a (D) , f (each of which is a function of x). Reexpressing the 21 The specification of appropriate integration constants for this, and all the other, radial ODEs will be discussed in subsection 3.7. 22 These asymptotic expansions hold for D > 2. For D = 2 (three-dimensional gravity), expansions in 1/r terminate. Exact solutions to Einstein's equations (3.12), with a flat boundary metric (3.7), are given by Σ = r + λ, A = 1 2 (r + λ) 2 − ∂0(λ + χ), and F = −∂1(λ + χ), with λ = λ(x 0 , x 1 ) completely arbitrary and χ = χ(x 0 , x 1 ) an arbitrary solution of the free wave equation, ∂ 2 χ = 0. result (3.6) for the stress-energy tensor using our renamed metric functions, we have 23 One must solve Einstein's equations throughout the bulk to determine the coefficients a (D) , f ij ; our procedure for doing so will be discussed in the next subsection. But λ(x) is determined by fiat -one must simply adopt some scheme for fixing λ.
One seemingly natural approach is to demand that λ vanish identically. That is, one could require that Σ(x, r)−r vanish, for all x, as r → ∞. However, this turns out to be a bad choice as it leads to apparent horizons whose radial positions vary rapidly with x. Such variation causes greater difficulty with numerical loss of precision due to cancellations between terms which grow large deep in the bulk. And it can lead to situations where the radial coordinate r decreases to zero and turns negative before the apparent (or Poincaré) horizon is reachedwhich is a nuisance since it makes the inverted radial coordinate u ≡ 1/r, which is otherwise convenient for numerical work, singular within the computational domain of interest.
A much preferable choice is to use the residual reparameterization freedom to put the apparent horizon at a fixed radial position, for all x. This choice makes the computational domain a simple rectangular region. If the surface r =r h is an apparent horizon, then an outgoing null geodesic congruence, normal to the surface and restricted to a t = const. slice, will have vanishing expansion [52,53]. This translates, in our metric ansatz, to a condition on d + Σ at the apparent horizon. 24 One finds: and all fields evaluated at radial positionr h . 25 We want condition (3.44) to hold at all times. It is convenient to regard this as the combination of a constraint on initial data (which is implemented by finding the radial shift 23 Becauseĝ has unit determinant, the sub-leading coefficientĝ is automatically traceless (as well as symmetric). So the full stress-energy tensor (3.42) of the dual field theory is automatically traceless as well. 24 Ref. [54] has a particularly nice treatment of null congruences. The congruence may be defined as kα(x) = µ(x)φ(x),α where, within the time-slice of interest, the surface φ(x) = C for some value of the constant C will define the apparent horizon. Requiring that k be null fixes the time derivative ∂tφ in terms of spatial derivatives of φ. Requiring the congruence to satisfy the (affinely parameterized) geodesic equation k α k β;α = 0 determines the time derivative of the multiplier function µ in terms of its spatial derivatives. Given these time derivatives, one may then compute the expansion via θ = ∇·k. Demanding that the result vanish on the surface φ = C gives the condition that this surface be an apparent horizon. eq. (3.44) is the result of specializing this condition to the case φ = r, so that the surface under consideration lies at a fixed radial position. 25 This expression and the subsequent horizon stationarity condition (3.47) are written using ordinary spatial covariant derivatives, not our modified derivatives (3.21). These gauge fixing conditions are, by necessity, not invariant under radial shifts and do not have simpler forms when written using the modified derivative ∇.
(3.16) needed to satisfy condition (3.44) at the initial time t 0 ), together with the condition that the horizon position be time-independent, ∂r h /∂t = 0, which requires that the time derivative of condition (3.44) hold at all times, Evaluating this horizon stationarity condition [and using eqs. (3.36), (3.37), and (3.39) to simplify] leads to a second order linear elliptic differential equation for A on the horizon. Explicitly, one finds with R (ν) the spatial Ricci scalar.

Integration strategy
The set of equations (3.33)-(3.39) have a remarkably convenient nested structure, which permits a simple and efficient integration strategy. On some given time slice t 0 , eq. (3.33) is a linear (in Σ) second order radial ODE which may be integrated to determine Σ(t 0 , x, r), providedĝ is already known on the time slice t 0 . Linearly independent homogeneous solutions behave as r 1 and r 0 as r → ∞. Consequently, the two needed integration constants may be fixed using the leading and first sub-leading terms in the asymptotic behavior, Σ ∼ r + λ + · · · [c.f. eq. (3.41b)]. However, this implies that λ(t 0 , x) must be known, in addition toĝ ij (t 0 , x, r), to determine Σ(t 0 , x, r).
Once Σ andĝ are known at time t 0 , the set (3.34) of second order radial ODEs can be integrated to determine the D−1 components F i (t 0 , x, r). 26 Linearly independent homogeneous solutions behave as r 2 and r 2−D as r → ∞. Consequently, the needed integration constants may once again be fixed from the leading and first sub-leading terms in the asymptotic behavior, This assigns a vanishing coefficient to the r 2 homogeneous solution, and a specified coefficient f to the other homogeneous solution. Hence, in addition toĝ and λ at time t 0 , one must also know the subleading coefficient f (D) i (t 0 , x) before integrating the F equations; how to accomplish this will be discussed momentarily.
Next up is eq. (3.36), which is a first order linear radial ODE for d + Σ, whose coefficients and source term depend on the already-determined values ofĝ, Σ, and F at time t 0 . Note 26 In addition to the manifest dependence on F in the first terms of equation (3.27), the second and third terms in the equation generate, through the modified covariant derivatives, terms which depend linearly on F . To solve for F , it is convenient to use the equivalent form (C.4) which uses ordinary covariant derivatives. In the absence of bulk sources, one may decouple the equations for different components of F by integrating first to find G 1/2 G ik (F k ) , extracting (F k ) , and then re-integrating to find the contravariant components of F . that, with time derivatives rewritten in terms of d + , this equation has no explicit dependence on A. The homogeneous solution behaves as r 2−D as r → ∞, so the single needed integration constant may be fixed by the coefficient of the sub-leading asymptotic term, Hence, in addition toĝ, λ, and f (D) i at time t 0 , we also require that the subleading coefficient a (D) (t 0 , x) be known before integrating the d + Σ equation; how to accomplish this will also be discussed momentarily. Now consider eq. (3.37). This is, in general, a set of coupled first order linear radial ODEs for the 1 2 D(D−1) − 1 components of the traceless symmetric spatial tensor d +ĝij . The coefficients and source terms of these equations again depend only on the already-determined values ofĝ, Σ, F , and d + Σ at time t 0 . The homogeneous solution to this equation behaves as r (1−D)/2 as r → ∞; the needed integration constant just corresponds to demanding the absence of any such homogeneous piece, so that d +ĝij = O(r 1−D ) as r → ∞.
Next turn to eq. (3.35), which is a trivial second-order linear radial ODE for A, with a source term depending on the already-determined values ofĝ, Σ, F , d + Σ, and d +ĝ . Linearly independent homogeneous solutions are r 1 and r 0 . The asymptotic behavior A ∼ , shows that knowledge of λ and ∂ t λ (at time t 0 ) determines these integration constants. If one fixes the residual reparameterization invariance (3.16) by choosing, a-priori, the value of λ as a function of both t and x, then this choice determines the two constants needed to integrate eq. (3.35) for A.
However, as discussed in section 3.6, it is preferable to adjust λ dynamically so as to fix the radial position of the apparent horizon, which forms the IR boundary of the computational domain. As described above, the horizon position invariance condition, dr h /dt = 0, reduces to the second order linear elliptic differential equation (3.47) for A on the horizon. The functions (evaluated at a given time t 0 and radiusr h ) appearing in the coefficients and source term of this linear elliptic PDE have all been determined in earlier steps of the integration procedure. Solving the linear PDE (3.47) (with appropriate boundary conditions in the spatial directions) will determine the IR boundary value A(t 0 , x,r h ). This provides one of the two integration constants needed to integrate eq. (3.35) and determine A everywhere on the t = t 0 time slice; the second integration constant is fixed by the asymptotic behavior A ∼ 1 2 r 2 + λr + O(1) as r → ∞, showing that λ is the coefficient of the term linear in r.
After the determination of A in this manner, using the horizon-invariance condition, one may extract the time derivative of λ from the subleading asymptotic behavior (3.41a) of A. The needed term may be isolated most conveniently by combining A with d + Σ, as with corrections to the limit vanishing as O(r 2−2D ). The determination of A also allows one to extract t-derivatives from d + derivatives so that, on the t = t 0 time slice, one can now evaluate To recap, having started at time t = t 0 withĝ ij , λ, f , and a (D) , the above procedure allows one to evaluate the time derivatives ofĝ ij and λ. Using a suitable integration method (such as fourth-order Runge-Kutta), these time derivatives provide the information needed to determineĝ ij and λ on the next time slice at t = t 0 + , up to an error vanishing as a power of the time step (e.g., 5 for fourth-order Runge-Kutta). Appropriate choices for time integration methods are discussed below in subsection 3.12.
However, before one can repeat the entire procedure above on the t 0 + time slice, one must also evaluate the time derivatives of asymptotic coefficients f from the first subleading term in the large r asymptotic behavior of ∂ t F . Likewise, ∂ t a (D) could be obtained by integrating the final radial ODE (3.39) to find ∂ t d + Σ, and then extracting ∂ t a (D) from its subleading asymptotic behavior. However, there is a simpler, far more efficient approach: direct use of boundary stress-energy conservation (2.6). As indicated in eq. (3.42), up to a common overall factor, − 2D−2 D a (D) is the energy density (and the trace of the spatial stress tensor), f  This completes the series of steps needed to turn Einstein's equations into an algorithm for evolving information from a given t = t 0 null slice to a subsequent slice at t 0 + . It should be emphasized that although one is solving the highly non-linear Einstein equations, this approach breaks the central time-evolution process down into a sequence of steps which only require solving the linear first or second order radial ODEs (3.33)-(3.37), plus the linear elliptic horizon PDE (3.47). The specific procedure described above is not, however, unique. Instead of treating Σ as an auxiliary field, to be computed anew on each time slice using the Schrodinger-like eq. (3.33), as mentioned earlier one could choose to treat Σ as a dynamical field which is evolved by computing d + Σ and then extracting ∂ t Σ. Likewise, the vector F could be treated as a dynamical field and evolved using eq. (3.38), instead of computing it as an auxiliary field from the second order eq. (3.34). One could fix the integration constant in eq. (3.36) for d + Σ using the planar horizon condition (3.44) directly on every time slice, instead of using (and evolving) a (D) to fix the subleading large r asymptotic behavior of d + Σ. These are just a few of the possibilities.
Different choices, while formally equivalent, have differing sensitivities to discretization effects and lead to algorithms with quite different numerical stability. Our experience is that stability is improved by computing auxiliary fields afresh on each time slice (instead of dynamically evolving these fields), and by using boundary stress-energy conservation to evolve the relevant subleading asymptotic coefficients directly, as described in the above scheme.

Initial data
To start the integration procedure, one must specify the spatial dependence of the asymptotic coefficients a (D) (t 0 , x) and f (D) i (t 0 , x) on some initial t = t 0 time slice. And one must specify the radial and spatial dependence of the rescaled spatial metricĝ ij (t 0 , x, r). The asymptotic behavior ofĝ ij (specifically the coefficientĝ In practice, there are several options for selecting initial data. One can choose to study "incoming" scattering states which, at time t 0 , contain well-separated excitations that, if considered in isolation, would have simple known evolution. Our case study of colliding planar shock waves in section 4.2 is an example of this type. Alternatively, one can start with a known static (or stationary) geometry describing an equilibrium state in the dual theory and then, after the initial time t 0 , drive the system out of equilibrium using time-dependent external sources. This was the approach used in refs. [24,27], where specified time-dependent boundary geometries represent sources coupled to T µν . Finally, one can simply make an arbitrary choice for the radial dependence ofĝ ij on the initial time slice. To a large extent, features inĝ ij deep in the bulk quickly disappear behind the horizon and have little influence on the future geometry; they reflect initial transients.
Given some choice of initial data, before proceeding with the integration strategy outlined above one must first find the location r h (t 0 , x) of the (outermost) apparent horizon on the t = t 0 slice, and adjust the radial shift λ(t 0 , x) so that the apparent horizon on the initial slice lies at a prescribed radial position, r h (t 0 , x) =r h . This requires integrating eqs. (3.33), (3.34) and (3.36) with λ set to zero (or some other arbitrary choice), to obtain provisional solutions for Σ, F i and d + Σ on the initial slice. Using these functions, one can locate the outermost value of r (for each x) at which the apparent horizon condition (3.44) is satisfied, and then adjust λ, at time t 0 , to shift this radial position to the prescribed value. Explicitly, if initial data {ĝ ij (t 0 , x, r), λ(t 0 , x)} leads to an apparent horizon at r h (t 0 , x), then the shifted initial data,g will lead to an apparent horizon at the desired location r =r h .

Finite spatial volume
With finite computational resources, one needs a finite computational domain in all directions, including the D−1 spatial directions. 27 One could make an r-independent change of variables which would map the range of the spatial {x i } coordinates to a finite interval, while preserving the form of the metric ansatz (3.14). Such r-independent transformations are part of the residual diffeomorphism freedom. However, we have not found such remapping to be desirable, as this leads to equations which are singular and ill-behaved at the ends of the spatial interval. A simple alternative which does not degrade numerical accuracy or stability is compactification of the spatial directions. We impose simple cubic periodic boundary conditions in spatial directions, with period L s . This should be viewed as a complementary part of the IR cutoff needed for computation. This spatial compactification also dictates the appropriate boundary conditions to use in solving the horizon invariance condition (3.47), namely spatial periodicity of A h .
Of course, compactification of spatial directions can have undesirable consequences. In scattering problems, as outgoing excitations separate there will be a limited time duration before the evolution is polluted by "wrap-around" effects caused by the compactification. If one is interested in exploring the uncompactified dynamics for some time duration τ , then one will generally need a spatial compactification with size L s ≥ c τ .

Field redefinitions
For numerical work, it is helpful to make a change of variable which maps the unbounded radial coordinate r to a finite interval. We just invert, and define In all the radial ODEs (3.33)-(3.37), the endpoint u = 0 (or r = ∞) is a regular singular point. As shown in eq. (3.41), the metric functions A and Σ, as well as the time derivative d + Σ, diverge as u → 0. For numerical purposes, it is very helpful to define subtracted functions in which the (known) leading pieces which diverge as u → 0 are removed, and to rescale the subtracted functions by appropriate powers of u so that the resulting functions vanish linearly, or approach a constant, as u → 0. This diminishes the substantial loss of precision which can occur due to large cancellations between different terms near u = 0. Altogether, 27 Problems with translation symmetry in one or more spatial directions, such as our first two examples below, are trivial exceptions to this assertion. One needs a finite computation domain in all directions in which solutions of interest have non-trivial variation.
this has lead us to use the following redefined fields in much of our numerical work: 28 Writing Σ 2 , and not just (u −1 +λ) 2 , in the subtraction terms for A and d + Σ is an arbitrary choice which makes no practical difference as Σ coincides with u −1 +λ up to O(u 2D−1 ) terms which are negligible near the boundary. The resulting u → 0 boundary conditions for these redefined fields are:

Discretization
To integrate the radial ODEs (3.33)-(3.37), and the horizon equation (3.47), one must discretize the radial and spatial coordinates, represent functions as finite arrays of function values on some specified set of points, and replace derivatives with suitable finite difference approximations.
Complications arise from the fact that u = 0 is a singular point in all the radial ODEs. Typical numerical ODE integrators (involving short-range finite difference approximations) do not tolerate such a singular point at the endpoint of the computational interval. One must introduce some finite separation scale u min , use truncated (analytically derived) asymptotic expansions to approximate functions in the near-boundary region 0 < u < u min , and only use numerical integration for u > u min . To achieve accurate results one must carefully select u min , and the order of the asymptotic expansion, so that the (in)accuracy of the truncated asymptotic expansion is comparable to that of the numerical integration. As one uses progressively finer discretizations (together with suitably matched improvements in the treatment of the asymptotic region), the gain in accuracy scales, at best, as a power of the radial discretization, error ∼ (∆u) k , with the exponent k depending on the range of the chosen finite difference approximation.
For many differential equations, substantially improved numerical accuracy can be obtained by using spectral methods. 29 This approach entails the use of very long-range approximations to derivatives. In essence, one represents functions as linear combinations of a (truncated) set of basis functions, and then exactly evaluates derivatives of these functions. For functions periodic on an interval of length L s , the natural basis functions are complex 28 If one introduces an explicit parameterization forĝij which solves the unit determinant constraint, as we do below in the examples discussed in section 4, then the redefinitions (3.53) forĝij and d+ĝij are replaced by analogous rescaling of the individual functions parameterizingĝij and their time derivatives. 29 For a good introduction to spectral methods, see ref. [55].
exponentials, e iknx with k n ≡ 2πn/L s (or the equivalent sines and cosines), and the expansion is just a truncated Fourier series, For aperiodic functions on an interval, convenient basis functions are Chebyshev polynomials, T n (z) ≡ cos(n cos −1 z). For functions on the interval 0 < u < 1, the appropriate expansion reads This is nothing but a Fourier cosine series in the variable θ ≡ cos −1 (2u−1). In so-called pseudospectral or collocation approaches, one determines the expansion coefficients {α n } by inserting the truncated expansion (3.55) or (3.56) into the differential equation of interest and demanding that the residual vanish exactly at a selected set of points whose number matches the number of expansion coefficients. For the Fourier series (3.55), these grid points should be equally spaced around the interval, for m = 0, · · ·, M . Again, knowledge of the expansion coefficients {α n } is completely equivalent to knowledge of the function values {g m ≡ g(u m )} on the collocation grid points. In practice, one uses these function values, plus interpolation formula, which together exactly reproduce the truncated basis expansions (3.55) or (3.56). 31 30 The Chebyshev grid points (3.59) are simply the image, under the mapping u = 1 2 (1 + cos θ), of equally spaced points in θ which would be appropriate for a Fourier cosine expansion. This choice of grid points, which include the interval endpoints, is most convenient when dealing with the imposition of boundary conditions. 31 In brief, for each truncated basis expansion, one reexpresses the expansion in the form f (x) = m fm Cm(x) where the "cardinal" function Cm(x) is the unique function which (i) can be represented in terms of the same truncated basis expansion, and (ii) vanishes identically at all collocation grid points except the m'th point, where it equals unity [so that Cm(xn) = δmn]. Cardinal functions are essentially regularized delta functions. See ref. [55] for more discussion including (in appendix E of that reference) explicit formulas for the appropriate cardinal functions for the Fourier expansion (3.55) and the Chebyshev expansion (3.56).
For linear differential equations, spectral methods convert the differential equation into a straightforward linear algebra problem (albeit one with a dense coefficient matrix, not a banded or sparse matrix as would be the case when using short-range finite difference approximations). One key advantage of spectral methods is improved convergence. For sufficiently well-behaved functions, accuracy improves exponentially as the number of basis functions is increased. A second advantage is that one can directly apply spectral methods to differential equations with regular singular points, as long as the specific solution of interest is well-behaved at the singular point. See ref. [55] for further detail.
We have found the use of (pseudo)spectral methods to be quite advantageous. We use the Fourier series form (3.55) to represent functional dependence on periodic spatial coordinates, and the Chebyshev form (3.56) to represent functional dependence in the radial direction (using the inverted radial variable u). 32

Time integrator
As outlined above in Section 3.7, in our evolution scheme we choose to evolve the minimal set of fields Φ ≡ {ĝ ij , a (D) , f (D) i , λ}. Discretizing the geometry with N i grid points in the x i spatial direction and N u points in the radial direction, the fields in Φ constitute a total of The time evolution portion of the spatially discretized Einstein equations then takes the schematic form In other words, after discretizing the spatial and radial directions, Einstein's equations reduce to a large system of simple, first-order ODEs describing the time-evolution of Φ. Evaluating F[Φ] is tantamount to first solving the nested system of radial equations (3.33)-(3.37) to find d +ĝij , then using eq. (3.49) to extract the discretized field velocities ∂ tĝij from d +ĝij , and finally using eqns. (3.48) and (3.50) , and ∂ t λ. The first order system (3.60) of simple ODEs can integrated using a variety of numerical ODE solvers. For simplicity, we limit our discussion to non-adaptive constant time step schemes. 33 We have used both implicit and explicit evolution schemes. When using explicit 32 Convergence of the spectral approximation (3.56) with increasing order M is naturally related to analytic properties of the functions under consideration. For problems involving a flat boundary geometry, all metric functions have expansions about u = 0 in integer powers of u. After applying the field redefinitions discussed above, expansions of our unknown functions only involve non-negative powers of u. As noted in footnote 8, for problems involving a non-flat boundary geometry, and an even dimension D, the near-boundary expansion necessarily includes logarithmic terms. One can still usefully apply spectral methods in this case, provided one subtracts these log terms (to reasonably high order) in the field redefinitions. Convergence of the spectral expansion will be degraded and non-exponential, but the performance of spectral methods can still be superior to traditional short range discretization methods. 33 Employing adaptive time-step schemes is clearly advantageous for some problems. However, all the issues discussed below, involving trade-offs between stability, accuracy, and computational efficiency, remain relevant for more complicated adaptive schemes. For more extensive discussion of numerical methods for solving ODEs see ref. [56], or most any other book on scientific computing. time evolution schemes, stability of the resulting numerical evolution requires that one use a suitably small time step. The Courant-Friedrichs-Lewy (CFL) condition [57], required for stability, imposes an upper limit on the time step. For diffusive equations, the time step ∆t must satisfy D∆t ∆x 2 , where ∆x is the minimum spatial grid spacing and D is the relevant diffusion constant. For wave equations with unit propagation velocity, the time step must satisfy ∆t ∆x. (In general, the relevant condition is that the numerical domain of dependence of new field values must encompass the appropriate physical domain of dependence.) Gravitational evolution in asymptotically AdS spacetime contains both diffusive and propagating modes. Diffusive gravitational modes are holographically related to diffusive modes in the dual quantum field theory which describe the spreading of (transverse) momentum density or other conserved charge densities. In the gravitational description, diffusive modes characterize the behavior of conserved densities near the horizon, as seen in the membrane paradigm [59] for horizon dynamics. Consequently, diffusive behavior of gravitational modes predominantly occurs in the spatial directions, and not in the radial direction. Therefore, one CFL condition for the time step is D∆t ∆x 2 . (Near equilibrium, with some effective temperature T , the diffusion constant D = (2πT ) −1 [58].) Gravitational waves can propagate in both radial and spatial directions, and outward-going radial waves propagate near the boundary with ∂u/∂t 1. Hence, the time step ∆t must also satisfy the propagating wave CFL conditions ∆t ∆x and ∆t ∆u. 34 In various applications, we have obtained good results using a third order Adams-Bashforth method as well as both implicit and explicit fourth order Runge-Kutta methods. Which solver is best depends on available computing resources, desired accuracy, and stability. Adams-Bashforth methods have the advantage that only one evaluation of F[Φ] is needed per time step. However, stability can require a very small time step. Explicit fourth-order Runge-Kutta methods require four evaluations of F[Φ] per time step, but are more stable than Adams-Bashforth methods and allow use of a larger time step. Implicit fourth-order Runge-Kutta methods are much more stable than explicit evolution. Moreover, with implicit evolution the time step need not satisfy the CFL condition. However, as we discuss below, implicit evolution requires many evaluations of F[Φ] per time step, which is costly.
Runge-Kutta methods, either implicit or explicit, require the computation of a set of "field velocities" {k i }, i = 1, 2, · · ·, M , defined by where ∆t is the time step, Φ n ≡ Φ(t n ), and α ij is an M × M matrix which determines the particular Runge-Kutta method. Once the set of M velocities {k i } have been evaluated at time t n , the new fields at time t n+1 ≡ t n + ∆t are given by for a set of coefficients {b i } which again depend on the particular Runge-Kutta method employed.
For explicit evolution, we use the classic fourth order Runge-Kutta (RK4) method for which For implicit Runge-Kutta methods, the matrix α ij is not lower triangular and the set of equations (3.61) implicitly define the different field velocities. When using implicit evolution, we compute the k i iteratively. Specifically, we start with a guess for the k i (e.g., the values of k i at the previous time step) and compute k i ≡ F[Φ n + α ij k j ∆t]. After evaluating an error norm ∆ ≡ |k i − k i |, we set k i = k i , reevaluate k i ≡ F[Φ n + α ij k j ∆t], and repeat the processes until ∆ approaches zero to within a chosen accuracy threshold. For this very simpleminded iterative process, the maximum time step is limited by convergence of the iterative scheme, and not by stability of the actual numerical evolution in time. The particular implicit Runge-Kutta method we employ is a fourth-order method known as Lobatto IIIC for which α ij =    With this third-order method, errors scale as O (∆t) 4 . Since the AB3 method requires knowledge of F[Φ] on three consecutive time slices, one must use some other scheme to compute Φ for the first two steps. This initialization can be performed using the abovedescribed explicit fourth-order Runge-Kutta method.
In general, when using non-adaptive integrators we recommend either implicit RK4, or explicit RK4 with suitably small time step, if computational resources (and patience) allows, and using AB3 if anything better is too slow. For problems where characteristic time scales lengthen as the evolution proceeds, use of an adaptive integrator (such as one which incorporates and compares RK4 and RK5 steps) is a reasonable choice. For more discussion of performance, see section 3.16.

Filtering
In addition to the CFL instabilities discussed above, discretization of non-linear PDEs can create spurious mechanisms, absent in the continuum limit, that cause artificial, unphysical growth in the amplitudes of short wavelength modes. This unphysical excitation of short wavelength modes leads to a progressive loss of accuracy and may eventually cause complete breakdown of the numerical evolution.
This problem is referred to as "aliasing," or "spectral blocking". 35 To understand how short wavelength modes can become artificially excited in discretizations of non-linear equations consider, for example, the product of two functions, f for which p−k is an integer multiple of 2πM . One says that the high momentum mode with |p| > M has been "aliased" to the low momentum mode with k = p − 2mM (for some integer m).
When computing the time evolution of non-linear PDEs, spectral aliasing typically leads to continuing unphysical growth in the amplitudes of modes near the |p| = M UV cutoff. 36 This growth of short wavelength modes due to aliasing is called spectral blocking. The same phenomena occurs when employing a basis of Chebyshev polynomials. Spectral aliasing can cause truncation error to grow unboundedly, and lead to time evolution becoming numerically unstable [55]. To make numerical evolution stable, for many PDEs, it is necessary to introduce 35 For more extensive discussion of spectral blocking see, for example, ref. [55]. 36 The power spectrum of Fourier coefficients of the exact solution will fall with increasing magnitude of the wavenumber for |p| ≥ M , provided the solution is smooth on the scale of ∆x = 2π/(2M +1). Consequently, of the modes which suffer from aliasing, the largest amplitude modes are those just slightly above the UV cutoff at |p| = M , and these modes are aliased to modes lying just slightly below |p| = M . Therefore, aliasing predominantly transfers power which should have appeared in modes above the UV cutoff to modes just below the cutoff -amplitudes of these modes receive the the greatest damage due to aliasing. some form of artificial dissipation which damps short wavelength modes. This can take the form of explicit addition of higher derivative terms ("numerical viscosity") to the equations of motion (as we did in ref. [24]). Or one can just selectively filter high k modes whose amplitudes are badly affected by spectral blocking [55].
For gravity, which is highly non-linear, one might expect significant aliasing and resultant spectral blocking. However, black branes in asymptotically AdS spacetime allow rapid dissipation of short wavelength modes, with an attenuation scale set by the infall time into the black brane's horizon. 37 Moreover, in infalling Eddington-Finkelstein coordinates, where lines of constant time t are infalling null geodesics, short wavelength modes can propagate into the black brane horizon instantaneously in coordinate time t, and hence need not persist and pollute the subsequent numerical evolution. As a result, for characteristic evolution of black brane geometries in asymptotically AdS spacetime, instabilities resulting from spectral blocking are less serious than might be expected.
Nevertheless, to ameliorate spectral blocking we have found it useful and often necessary either to introduce numerical viscosity, or to selectively filter short wavelength modes. For spatial directions, applying a sharp low-pass filter which sets to zero all Fourier components for which 2 3 k max < |k i | ≤ k max , for any spatial direction i, is a simple, computationally efficient choice. 38,39 For controlling spectral blocking in the radial direction, our preferred method is filtering in real-space. Let {u fine } represent the radial spectral grid (3.59) used to solve Einstein's equations, and let {u coarse } represent a spectral grid with two thirds as many points in the radial direction. At each time step we interpolate the geometry from the fine grid to the coarse grid. The interpolation can be done without losing spectral accuracy by employing the spectral representation (3.56) to evaluate a function at off-grid locations. We then reinterpolate the geometry from the coarse grid back to the fine grid. As the coarse grid is also a spectral grid, the interpolation back to the fine grid can also be performed without losing spectral accuracy. The process of interpolation from fine to coarse and back to the fine 37 For example, high momentum quasinormal modes of black branes in asymptotically AdS spacetime decay on a time scale set by the gravitational infall time [13,14]. In this context, "high-momentum" applies to modes with rapid radial and/or spatial variations. 38 This is the "2/3's rule" [55,60]. For equations with only quadratic non-linearity, removing the uppermost third of Fourier components is sufficient to prevent aliasing from corrupting the components which are retained. Einstein's equations have higher order (cubic, quartic, and worse) non-linearities. Nevertheless, our experience is that filtering using the 2/3 rule is very effective in removing spectral blocking artifacts. The sufficiency of 2/3's rule filtering, despite high order non-linearities in the equations, is undoubtedly a reflection of the above-mentioned dissipation of short wavelength modes which is an intrinsic feature of black-brane geometries in asymptotically AdS spacetimes. 39 To implement this low-pass filter, one can use fast Fourier transforms (FFTs) to transform from real space to momentum space and back. Or one can construct the one-dimensional real space filter which is exactly equivalent to the desired momentum cutoff, and apply this filter as a convolution in real space. Asymptotically, for very fine spatial discretizations, using FFTs is most efficient. However, given matrix multiplication and convolution routines which are optimized for modern multi-core processors, the break-even point beyond which FFTs become preferable to real-space convolution can lie at surprisingly large values of the number of points Ni used in the discretization of a given spatial direction.
grid has the effect of filtering short wavelength modes (albeit with a soft cutoff, instead of a sharp momentum space cutoff). Filtering in real-space with the Chebyshev grid (3.59), which includes the boundary u = 0, has the advantage that Dirichlet boundary conditions at u = 0 are completely unaffected by filtering.

Parallelization
The characteristic formulation Einstein's equations presented above is easily amenable to parallelization. Imposition of the fixed horizon condition (3.44) makes the computation domain a simple rectangular box, with the resulting discretized spatial lattice a tensor product grid. For such a grid, spatial and radial derivatives of all functions can easily be computed in parallel. When computing, for example, the radial derivative of a function, with a tensor product grid one can evaluate the radial derivative independently at each point in space. Computation of the radial derivatives at a given spatial point (or set of points) can be performed independently by different processors. Moreover, as discussed above, Einstein's equations in the characteristic formulation take the form of linear ODEs in the radial coordinate. After all needed spatial and radial derivatives have been computed, these radial ODEs can be solved independently at each spatial point. In other words, the radial ODEs can be integrated, in parallel, using independent CPUs for each point in space. Solving the radial ODEs in parallel greatly increases computation speed; specific performance results will be discussed below in section 3.16.

Domain decomposition
Consider a discretization with N u points in the radial direction and N i points in the spatial x i direction (so the total number of grid points on a timeslice is N u N s , with N s ≡ D−1 i=1 N i the number of spatial discretization points). For a sufficiently fine spatial discretization, the rate limiting step in our time-evolution procedure is the solution of the linear elliptic PDE (3.47) for the function A at the apparent horizon. Employing spectral methods, solving eq. (3.47) requires the solution of a linear system with a dense coefficient matrix with order (N s ) 2 elements, which is the discretization of the linear operator appearing in eq. (3.47). Besides requiring extensive computation time, which scales as O(N 3 s ), memory consumption can become problematic for large N s . Fortunately, it is easy to ameliorate these difficulties.
Linear elliptic PDEs such as eq. (3.47) can be efficiently solved using domain decomposition [55]. In this procedure, the spatial interval in each x i direction is broken up into m i separate subintervals, thereby decomposing the spatial computational domain into a total of m s ≡ i m i subdomains. Let = 1, · · ·, m s index these subdomains. The x i dependence of functions within some subdomain can be represented as a sum of n Domain decomposition can also be usefully employed in the radial direction. This entails breaking the radial interval up into m u subdomains and coupling adjacent subdomains via boundary conditions. There are several reasons to employ domain decomposition in the radial direction. First, the global nature of spectral methods implies that if the metric happens to be badly behaved deep in the bulk, the spectral representation of the metric will converge poorly everywhere, including near the boundary. However, for many situations, including the planar shock collisions discussed below in section 4.2, as time progresses features in the geometry deep in the bulk can rapidly fall through the event horizon. In such situations, there is little practical value in knowing the metric with high spectral accuracy in such regions. By employing domain decomposition in the radial direction, spectral convergence in one subdomain is only weakly dependent, via boundary conditions, on spectral convergence in other subdomains. Consequently, domain decomposition helps improve convergence near the boundary when convergence is poor deep in the bulk.
Domain decomposition in the radial direction can also be helpful in controlling the effects of round-off error. Near the boundary, Einstein's equations contain 1/u 2 singularities. Since the grid is clustered around u = 0, round-off error coming from points close to u = 0 can be greatly amplified by the presence of the nearby singularity. Domain decomposition helps with this simply because it allows the grid spacing ∆u near u = 0 to be much larger than it will be if one uses a single domain with Chebyshev grid points.
The implementation of domain decomposition in the radial direction is completely analogous to the spatial decomposition of elliptic PDEs discussed above. The radial domain is split up into m u subdomains with n u Chebyshev grid points in each subdomain. The endpoints of adjacent subdomains coincide. The boundary conditions on the second order equations (3.33)-(3.35) for A, Σ, and F are simply that these functions, and their radial derivatives, are continuous across interfaces. When solving the first order radial equations (3.36) and (3.37) for d + Σ and d +ĝij , one constructs solutions for these functions which are continuous along infalling geodesics. When computing the time derivative ∂ tĝij via eq. (3.49), one requires that ∂ tĝij also be continuous along outgoing geodesics.

Performance
Key parameters controlling performance of a numerical calculation using our approach are, naturally, the number of points used in the spectral grids in the radial and spatial directions, the number of time steps which are taken, plus the speed (and memory capacity) of available computing resources.
With appropriate use of domain decomposition, both memory requirements and computational cost per time step are essentially linear in the total number of grid points (radial times spatial). For many problems, such as our homogeneous isotropization and 2D turbulence examples below, 20-25 points in the radial grid are sufficient. For our colliding shock example, we use up to 80 radial points (partitioned into multiple subdomains). The size and spacing of the spectral grid used for spatial directions inevitably depends on the nature of the chosen problem. For colliding shocks, we have used a Fourier grid with just over 800 points in the longitudinal direction; this allows us to evolve the outgoing remnants of the collisions quite far before wrap-around effects arise. 40 For 2D turbulence, we have used Fourier grids with several hundred points (in each direction). Here, the challenge is to make the spatial domain large enough to contain flows whose Reynolds number is in the turbulent regime.
For the following examples, and our prior work [24,27,30], we implemented the abovedescribed approach and performed calculations using MATLAB. It is possible that somewhat improved performance could be obtained by carefully programming in a lower level language. However, particularly for problems (such as our examples in secs. 4.2 and 4.3) where symmetries at most reduce the problem to 2+1D or 3+1D PDEs, the bulk of the computational time is spent in linear algebra routines which are already highly optimized in MATLAB. Consequently, we expect that any potential gain from coding in a lower level language is quite modest. (Far more important, for problems with non-trivial spatial dependence, is that one implements the approach in a manner which allows easy parallelization and hence benefits from multi-core processors.) Our calculations have been performed with only desktop or laptop scale computing resources. 41 Using these relatively limited computing resources, calculating the homogeneous isotropization example discussed below (where symmetries reduce the problem to 1+1D PDEs) is quick, taking a few seconds. Evolving the geometry in the colliding shock example (where symmetries reduce the problem to 2+1D PDEs), for the more demanding case of narrow shocks, required approximately 12 hours on a laptop computer. Performing the numerical evolution of the geometry in our turbulent fluid example (where one is dealing with 3+1D PDEs) required approximately three weeks of time.
A different aspect of performance concerns achievable accuracy. At a crude but important level, a key indicator of accuracy is the absence of obvious numerical instabilities which prevent continuing evolution to arbitrarily late times. Achieving stable evolution requires sensible choices for the spectral grids and time step. The use of UV filtering to control spectral blocking (as described in section 3.13) is important for many problems.
Once stable numerical evolution is achieved, a more refined, physically important, measure of accuracy involves comparison of numerical results with analytically derived late-time asymptotic forms. This is discussed below in the context of our specific examples. A further check of numerical accuracy can be obtained by monitoring the validity of constraint equations. As the constraint equations were used in deriving the horizon stationarity condition (3.47), which determines the value of A on the apparent horizon, correct numerical evolution of the gauge parameter λ via eqs. (3.47) and (3.48) is intimately connected to the numerical validity of the constraint equations. Hence, one simple test of the constraint equations comes from monitoring how well the gauge parameter obtained from evolving eq. (3.48) agrees with the value obtained by directly solving the horizon fixing condition (3.44). If the gauge parameter evolved from eq. (3.48) drifts too far from the value which correctly solves the horizon condition, it may be necessary to make small periodic readjustments in λ to ensure that the horizon remains at r = 1.
The bottom line is that with appropriate care (in adjusting grids, filtering, etc.), the achievable accuracy is, in our view, remarkably good.

Motivation
Relativistic heavy ion collisions may be regarded as proceeding through a sequence of stages. Initially, the collision of the (overlapping portions of the) highly Lorentz contracted nuclei may be viewed as liberating a very large phase space density of partons from the colliding nucleons. Within the central rapidity region of the event, the initial distribution of partons is highly anisotropic, with typical transverse momenta much larger than longitudinal momenta.
These partons subsequently interact and scatter. After a "thermalization time" (which, more properly, should be called an isotropization time), the gas of interacting partons may be modeled as a relativistic fluid -a quark-gluon plasma -whose stress tensor, in a local fluid rest frame, is nearly isotropic. This plasma expands, cools, and eventually reaches a temperature where hadrons reform, fly outward, and ultimately reach the detector. 42 Hydrodynamic modeling of the results of heavy ion collisions strongly suggests that the isotropization time of the dense parton gas is remarkably short, less than 1 fm/c [62], and that the resulting plasma behaves as a nearly ideal fluid. Understanding the dynamics responsible for such rapid isotropization in a far-from-equilibrium non-Abelian plasma is a challenge. The nearly ideal (i.e., low viscosity) behavior of the produced plasma is an indication that experimentally accessible quark-gluon plasma is strongly coupled [63]. 43 Due to the difficulty of studying real time quantum dynamics in QCD at strong coupling, it is useful to examine far-from-equilibrium behavior in an instructive toy model, namely N = 4 SYM, whose equilibrium behavior at non-zero temperature mimics many features of real QCD plasma. This was the motivation for our earlier study [24] of isotropization in spatially homogeneous but highly anisotropic states of strongly coupled N = 4 SYM, using the dual gravitational description. In that work, we considered initial states which could be produced by the action of time-dependent (but spatially homogeneous) background fields. The background field which naturally couples to the stress-energy tensor of the field theory is the metric of the four-dimensional geometry in which the QFT is formulated. A timedependent 4D metric in the QFT description corresponds, under the holographic mapping, to a time-dependent boundary geometry in the dual gravitational description. Our earlier work [24] solved the resulting gravitational dynamics (numerically), using the approach presented in section 3.7, with the simplifying assumption of spatial homogeneity but with the complication of a time-dependent boundary geometry.
In the present paper we focus, for simplicity, on problems involving a flat Minkowski boundary geometry. To illustrate the application of our methods, we will present results on far-from-equilibrium isotropization in which the operational driving via a time-dependent boundary geometry of ref. [24] is replaced by a simple (and arbitrary) choice of initial data for our characteristic formulation.

Setup
The boundary dimension D = 4. With the imposition of spatial homogeneity, spatial parity invariance, and O(2) rotation invariance, the only non-zero functions in the metric ansatz (3.14) are A, Σ, and the diagonal elements ofĝ ij which we write in terms of a single "anisotropy" function B, ĝ ij = diag(e B , e B , e −2B ) . (4.1) The unknown functions A, B, and Σ depend only on t and r. Eqs. (3.33), (3.35), and (3.36) for Σ, A, and d + Σ, respectively, become while eq. (3.37) for d +ĝ reduces to We replace the field redefinitions (3.53) involving the spatial metric with the redefinitions for the anisotropy function and its (modified) time derivative. The asymptotic behavior (3.41) implies that b vanishes at the AdS boundary whileḃ approaches a finite limit of −2b (4) . The latter boundary condition is imposed when solving eq. (4.3) forḃ. The apparent horizon condition (3.44) is just Since there is no spatial dependence, the horizon stationarity equation (3.47) becomes a simple algebraic condition for the value of A on the apparent horizon, Initial data consists of a choice of the anisotropy function on the initial time slice, B(t 0 , r), plus a value for the single asymptotic coefficient a (4) (t 0 ) which sets the initial energy density. We make the simple but arbitrary choice: with β = 5, u 0 = 0.25, w = 0.15, and α = 1. Using the result (3.42) for the boundary stress-energy tensor and inserting the holographic relation G N = π 2 L 3 /N 2 c appropriate for N = 4 SYM, the corresponding energy density T 00 = 3 8 N 2 c α/π 2 . The energy density of an equilibrium, strongly coupled N = 4 SYM plasma at temperature T is given by T 00 eq = 3 8 N 2 c π 2 T 4 . Hence, our chosen value of α corresponds to an equilibrium temperature T ≡ α 1/4 /π = 1/π.
To evolve the geometry, we use a spectral grid in the radial direction with 25 points, and employ explicit fourth-order Runge-Kutta for the time-integrator with a time step ∆t = 0.01. For this simple 1 + 1 dimensional problem, where all dynamics takes place in the radial direction only, we do not employ any filtering. Indeed, because radial lines are infalling null geodesics, any high frequency numerical noise generated by the numerical evolution tends to get absorbed effectively instantaneously by the horizon.

Results
The resulting evolution of the anisotropy function B is shown in the left panel of fig. 3. The right panel displays a plot of the pressure anisotropy δp ≡ T zz − 1 2 (T xx + T yy ), relative to the equilibrium pressure p eq = 1 8 N 2 c (πT ) 4 , as a function of time. Inserting the diagonal form (4.1) of the spatial metric into the general result (3.42) for the stress-energy tensor, one sees that the pressure anisotropy is simply proportional to the coefficient b (4) of the leading near-boundary behavior of the anisotropy function, B(t, u) ∼ b (4) (t) u 4 + O(u 5 ).
Examining fig. 3 one sees, first and foremost, that the geometry evolves toward an isotropic equilibrium geometry, which is just the static Schwarzschild black-brane solution. This is a basic test of the numerics; no problems with numerical instabilities, potentially preventing evolution to arbitrarily late times, are seen. The approach to equilibrium shows exponentially damped oscillations. With no spatial gradients, there is no excitation of hydrodynamic degrees of freedom, and hence no hydrodynamic regime in the response.
At sufficiently late times, the damped oscillations of the pressure anisotropy reflect the discrete spectrum of complex quasinormal mode frequencies characterizing infinitesimal departures from equilibrium [13,14], specifically those of = 2 metric perturbations whose linearized dynamics around the AdS-Schwarzschild black brane solution coincides with fluctuations of a minimally coupled scalar field. The late time asymptotic response has the form  where the first few quasinormal mode frequencies, at zero wavevector, are given by [13]: (4.9) As a check on the accuracy of the numerics, in fig. 4 we plot e |Re λ 1 |t δp/p eq , as well as a fit to the lowest quasinormal mode. As is evident from the figure, the rescaled amplitude of e |Reλ 1 |t δp/p eq is constant at late times. Indeed, our fit to the lowest quasinormal mode agrees with the numerics at the level of a part in 10 4 , or better, after time t = 10.
In terms of physics, perhaps the most significant result one sees from fig. 3 (and from the results of ref. [24]) is that the characteristic relaxation time is comparable or shorter than 1/T , even when the system is initially quite far from equilibrium with δp/p eq of O(10). The gravitational infall time in the AdS-Schwarzschild geometry is also order 1/T . This naturally suggests that, even far from equilibrium, one should regard the gravitational infall time as characterizing the relaxation time of non-hydrodynamic degrees of freedom.

Motivation
Collisions of infinitely extended planar shock waves in N = 4 SYM may be viewed as instructive caricatures of collisions of large, highly Lorentz-contracted nuclei. In the dual description of strongly coupled (and large N c ) SYM, this becomes a problem of colliding gravitational shock waves in asymptotically AdS 5 spacetime. In this section, we discuss the setup, preparation of initial data, and results for such planar shock collisions.
Numerical construction of a complete colliding planar shock geometry was first performed in ref. [30]. More recently, the authors of ref. [31] examined the sensitivity of the post-collision energy density and pressure distributions to the width of the initial shocks. In both of these previous works, a small background energy density was added to the initial data to help control numerical instabilities deep in the bulk.
Using the filtering approach discussed in sec. 3.13 to suppress spectral blocking, it is possible to compute, accurately, colliding shock geometries, even for very thin shocks, without adding any background energy density. In other words, it is possible to study collisions of shocks which are truly excitations of the vacuum state. Even with a vanishing background energy density (or temperature), we find no problems associated with caustics or non-planar horizon topology.

Initial data
The boundary dimension D = 4. With the imposition of spatial homogeneity in transverse directions, plus 2D rotation and reflection invariance in the transverse plane, the only nonzero functions in the metric ansatz (3.14) are A, Σ, the longitudinal component F z of the spatial vector F , and the diagonal elements of the rescaled spatial metricĝ ij . The latter we write in terms of a single anisotropy function B which distinguishes the transverse and longitudinal directions, .  We replace the field redefinitions (3.53) involving the spatial metric with the following field redefinitions for the anisotropy function and its (modified) time derivative, The asymptotic behavior (3.41) implies that b vanishes at the AdS boundary whileḃ approaches a finite limit of −2b (4) . The latter boundary condition is imposed when solving eq. (3.37) forḃ. We choose initial conditions corresponding to two well separated, smooth, non-singular planar gravitational waves. In Fefferman-Graham coordinates [denoted (t,x ⊥ ,z,ρ)] the precollision metric reads wherex ± ≡t ±z and h(z) is an arbitrary function characterizing the longitudinal profile of the shocks. In what follows we choose a simple Gaussian profile of adjustable width and amplitude, parameterized as (4.14) In the distant past, the geometry both between and far away from the shocks is just AdS 5 , up to negligible, exponentially small corrections. Via eq. (2.5), the initial boundary energy density and longitudinal stress are given by while the momentum density Therefore, the metric (4.13), with shock profile (4.14), describes two localized planar lumps of energy of width w moving toward each other at the speed of light and colliding at timet = 0.
Restoring the overall factor of κ = L 3 /(4πG N ) [c.f. eq. (2.4)] and inserting the holographic relation G N = π 2 L 3 /N 2 c , appropriate for N = 4 SYM, shows that the energy per unit area of each incoming shock is µ 3 (N 2 c /2π 2 ). Without loss of generality we may set µ = 1 and measure all quantities in units of µ to the appropriate power. We will present results for the collisions of "wide" shocks with w = 0.375, and "narrow" shocks with w = 0.075. (For comparison, ref. [30] used w = 0.75 and ref. [31] investigated widths ranging from 1.9 down to 0.05.) In the distant past, when the two functions h(t ±z) have negligible overlap, the metric (4.13) is arbitrarily close to an exact solution to Einstein's equations (3.12). 44 But near the collision timet = 0, when the functions h(t ±z) begin to overlap significantly, the metric (4.13) ceases to be a (near) solution to Einstein's equations, and one must compute the future evolution numerically. To do so we employ our characteristic formulation.
To obtain initial data suitable for our formulation, the initial metric (4.13) must be transformed from Fefferman-Graham coordinates to infalling Eddington-Finkelstein coordinates, in which the metric takes the form (3.14). To do so we compute, numerically, the needed coordinate transform for a single shock moving in the +z direction and thereby determine the set of functions {b + (t−z, u), a z− (t+z), λ − (t+z)} for a left-moving shock. As we discuss in greater detail below, we then superimpose the pre-collision functions and likewise for f z , λ} forward in time by numerically solving Einstein's equations.
The metric of a single shock moving in the +z direction is given by [64]  The coordinate transformation taking this metric to the Eddington-Finkelstein form (3.14) (with u ≡ 1/r) can be expressed as for suitable functions α, β, and γ whose determination will be described momentarily. The functions {b + , a z+ , λ + } providing the required initial data for our characteristic formulation can be expressed in terms of the profile function h and the transformation functions α, β, and γ. A short exercise shows a (4) The equations determining the coordinate transformation functions (which follow from solving for infalling radial null geodesics in the metric (4.18), or equivalently demanding that the transformed metric have the desired form (3.14)) are simplified by redefining In terms of ξ, δ and γ, the equations of the coordinate transformation reduce to a system of coupled radial ODEs for ξ and δ, with H ≡ h t−z + u + δ − u 2 ξ/(1 + uξ) . The function γ satisfies the first order radial ODE As is clear from inspecting eq. (4.22), on each slice with fixed t−z, the functions ξ and δ can be determined by integrating the coupled ODEs (4.22) from u = 0 to u = 1. With ξ and δ known, one can then integrate the radial ODE (4.23) on the same t−z = const. slice to determine γ. In other words, determination of the coordinate transformation functions is local in t−z; one need only integrate radial ODEs. The desired solutions to eqs. (4.22) and (4.23) are specified by boundary conditions at u = 0 and u = 1. The conditions t =t, z =z, and u =ρ near the AdS boundary imply that the fields ξ, δ, and γ have the asymptotic forms (4.24) Defining further rescaled fields ∆ ≡ δ/u 4 and Γ ≡ γ/u 4 , we therefore impose at the AdS boundary the conditions and integrate eq. (4.22) to find ξ and ∆, and eq. (4.23) to find Γ. One additional boundary condition is needed to fully specify a solution to eq. (4.22). At u = 1 we impose the condition for some choice of the functionρ max (t−z). This boundary condition determines how deep into the bulk one determines the transformation of the initial geometry. Via eqs. (4.19) and (4.21), one sees that at u = 1 the Fefferman-Graham coordinateρ coincides withρ max . The boundary condition (4.26) also largely determines the gauge parameter λ + since, away from the shock where h is negligible, one has Controlling how deep into the bulk one solves for the initial geometry (in Eddington-Finkelstein coordinates) is essential. If one integrates too far into the bulk, the metric functions become very large, causing problems with loss of numerical precision. This can already be seen in the single shock Fefferman-Graham metric (4.18), where metric functions grow likẽ ρ 2 for largeρ. However, the apparent horizon of the colliding shock geometry exists prior to the collision at t = 0 [30]. The mapping of the initial geometry into Eddington-Finkelstein coordinates must go sufficiently deep into the bulk so that the apparent horizon lies within the chosen computational domain u ∈ [0, 1]. Selecting an appropriate value forρ max so that the apparent horizon lies in this interval, and the bulk geometry is reasonably well behaved, can require some trial and error. We setρ max = 8, (4.28) independent of t−z, and comment below on more refined choices ofρ max (t−z). We employ domain decomposition in both the radial and longitudinal (t−z) directions when solving eqs. (4.22) and (4.23). We use 20 Chebyshev polynomials in each subdomain in both directions. We employ 350 subdomains in the t−z direction and 35 subdomains in the u direction, and solve the equations in the interval −18 ≤ t−z ≤ 18. Using domain decomposition in each direction is advantageous for several reasons. First, as mentioned above, the coordinate transformation can become badly behaved deep in the bulk. As discussed in sec. 3.15, if the convergence of the spectral series very deep in the bulk becomes poor, the use of domain decomposition serves to reduce the influence of such poor convergence on fields closer to the boundary. Second, the use of domain decomposition -with relatively few points in each subdomain -allows the function b + [defined in eq. (4.20b)], and its near-boundary On the boundary, u = 0, the shock is centered at z = 0 at the time shown, t = 0. However, in Eddington-Finkelstein coordinates the shock increasingly extends into the +z direction as one goes deeper into the bulk. This also manifests itself in the gauge parameter λ, which differs significantly from its background value in front of the shock. In regions where b + = 0 the geometry is that of AdS 5 . asymptotics, to be determined numerically with very good and controllable accuracy. In particular, the use of domain decomposition allows finely spaced grid points to be used for rapidly varying functions, thereby enabling good spectral convergence, while simultaneously avoiding the significant round-off error that can occur when employing a single global domain with a very large number of grid points. Fig. 5 plots the resulting functions b + and λ + , at time t = 0, for a single narrow shock moving in the +z direction. One sees that b + is non-zero for positive values of z (well beyond the width of the shock) deep in the bulk. Likewise, the gauge function λ + differs significantly from its background value far in front of the shock. This behavior is an unavoidable consequence of our use of infalling Eddington-Finkelstein coordinates, combined with the fact that, in Fefferman-Graham coordinates, the perturbation to the geometry due to the shock extends arbitrarily deep into the bulk at any fixed value oft−z lying within the shock profile. Any radially infalling null geodesic which begins at the boundary at t = 0 and some z w eventually intersects the shock which is moving in the +z direction with unit speed. Since all events along such a geodesic have common values of the Eddington-Finkelstein coordinates t and z this shows that, for any z > 0, sufficiently deep in the bulk, metric functions are influenced by the shock.
In the neighborhood of slices with fixed t−z w, on which a (4) + and f (4) z+ are negligible, Einstein's equations imply that the local geometry is AdS 5 as long as b + is also negligible. The geometry only ceases to be AdS 5 deep in the bulk where b + becomes non-negligible. In other words, the geometry corresponding the dark red "background" region in the left panel of fig. 5 is simply that of AdS 5 .
The fact that the functions b + and λ + are non-zero in front of the shock may appear to constitute a problem for computing the initial geometry of two colliding shocks in Eddington-Finkelstein coordinates. Incoming shocks which, near the boundary, have arbitrarily large separation at some initial time are, in Eddington-Finkelstein coordinates, already colliding sufficiently deep in the bulk. In other words, even for shocks which are widely separated on the boundary, there will always be some region deep in the bulk where the simple superposition (4.17) of the functions {b ± , a (4) ± , f (4) z± , λ ± } is not correct. However, as we demonstrate below, when the shocks are well separated on the boundary, the region where the functions b ± overlap significantly, and hence where the shocks are already colliding deep in the bulk, lies inside the apparent horizon and thus is causally disconnected from the above-horizon geometry. Just as seen in the Fefferman-Graham metric (4.13), the initial above-horizon geometry both between and far away from the shocks is simply AdS 5 .
Using b − (t+z, u) = b + (t−z, u), we superimpose b ± (t, z) and define the initial anisotropy function to be evolution of the geometry, we choose to use a Fourier grid in the z direction with N z points, with periodicity enforced at z = ±z max with z max ≡ 10. For narrow shock collisions we use N z = 801 and for wide shock collisions we use N z = 401. We use domain decomposition in the radial direction with 4 domains, each having 20 Chebyshev points. After computing the functions {b, a (4) , f (4) z , λ} on the new grid, we then apply a radial gauge transform to reposition the apparent horizon to radial coordinate u = 1.
Before proceeding, we address two more technical points. First, in the infinite volume (z max → ∞) limit, the apparent horizon asymptotes to the Poincaré horizon at Fefferman-Graham coordinateρ → ∞. In this limit, our choice of constantρ max in eq. (4.28) will not yield the the entire above-horizon geometry in the computational domain 0 ≤ u ≤ 1. This can present a problem since, for any finite choice of z max , one cannot compute the location of the apparent horizon and thereby know how bigρ max should be until the functions {b, a (4) , f (4) z , λ} are computed (which requires a choice ofρ max ). However, the above-horizon pre-collision geometry at large |z| is simply AdS 5 . Because of this, one may freely adjust λ(t 0 , z) at large |z| without changing the initial geometry. In other words, one may make the redefinition with W (z) = 1 in the vicinity of the shocks and W (z) arbitrary at large |z|. This freedom allows one to compute and superpose the single shock profiles usingρ max = const., and then apply the transformation (4.33) with W (z) chosen such that the apparent horizon lies in the computational domain for any choice of z max . With this technique,ρ max need only be chosen large enough such that the horizon lies in the computational domain u ≤ 1 near z = 0. We employ this technique and parameterize W (z) via with K, z 0 , and s adjustable parameters. We choose K = 21 and s = 0.25. For narrow shocks we use z 0 = 3, and for wide shocks we use z 0 = 6. Second, after computing {b, a (4) , f (4) z , λ} on the grid used to solve Einstein's equations, but before gauge transforming to reposition the apparent horizon at u = 1, we have found it advantageous to filter high momentum modes. This helps eliminate numerical noise generated in the numerical calculation of b and λ. We perform the filtering by Fourier transforming b and λ in z and then setting the coefficients of modes with momentum |k| ≤ k max /2 to vanish. Fig. 6 shows the resulting gauge transformed initial anisotropy function b and gauge parameter λ for incoming narrow shocks with width 0.075. The apparent horizon is at u = 1. In between the shocks, the functions b, a (4) , and f (4) differ negligibly from zero. As mentioned above, in the neighborhood of z = const. slices on which a (4) + = f (4) z+ = b + = 0, Einstein's equations imply that the local geometry is AdS 5 ; only deep in the bulk where b becomes significant does the geometry deviate from AdS 5 . Therefore, the geometry corresponding the background dark red region in the figure (everywhere except in the vicinity of the shocks) is simply that of AdS 5 . Exactly the same description holds for the wide shock initial data.

Results
Fig 7 displays E ≡ T 00 , the energy density rescaled by a factor of κ = N 2 c /(2π 2 ), for both wide (top) and narrow (bottom) shock collisions. The shocks approach each other at the speed of light in the ±z direction and collide at z = 0 at time t = 0. For both cases, the debris leaving the collision event appears dramatically different than the initial incoming shocks. Prior to the collision, all the shock energy lies near the lightcone (smeared only by the width of the shock), while long after the collision nearly all the energy lies inside the lightcone. Inspecting fig. 7, one sees qualitative differences between narrow and wide shock collisions. For wide shocks, there is no sign of any distinct remnant of the shock remaining on the forward light cone; the energy density of the post-collision debris is smoothly distributed in the interior of the forward light cone [30]. In contrast, for the narrow shock collisions there are clear remnants of the initial shocks propagating outward on the forward light cone [31]. But, as can easily be seen in fig. 7, immediately after the collision energy density is transported inside the lightcone and the portion remaining very near the lightcone steadily attenuates. On the left side of fig. 8 we plot the amplitude A of the energy density on the lightcone as a function of time for the narrow shock collisions. At late times our results are consistent with the power-law decay A ∼ t −0.9 . By time t = 9, the amplitude of the null maxima has decreased to 13% its pre-collision value. Evidently, for both wide and narrow shocks the collision event results in the subsequent annihilation of the shocks with essentially all energy lying well inside the forward light cone at late times. Aside from the decay of the null peaks in the energy density, there is another qualitative difference between collisions of narrow and wide shocks. On the right side of fig. 8 we plot the energy density for the narrow shock collision at successive times t = 1, 2, 3, 4, 5. As is evident from the figure, there is a brief period of time after the collision when the energy  density just behind the receding null peaks is locally negative [31]. However, by time t = 4 the energy density is everywhere positive, just as it always is for wide shock collisions. Evidently, the presence of negative energy density is a transient effect. Indeed, as shown in fig. 9, aside from the decaying null maxima on the light cone, at late times the distribution of energy density produced by both wide and narrow shock collisions looks quite similar. It is instructive to compare our results with predictions from the fluid/gravity correspondence [12]. In the limit of asymptotically slowly varying fields (compared to the dissipative scale set by the local temperature T of the system) Einstein's equations (3.12) can be solved perturbatively with a gradient expansion M N (x, r) , (4.35) where g (n) M N is of order (∂/∂x µ ) n in boundary spacetime derivatives [11]. Via eq. (3.6), this implies that the boundary stress tensor also admits a gradient expansion. In ν = D−1 spatial dimensions, the resulting gradient expansion of the boundary stress begins where ε is the (rescaled) proper energy density, u the fluid velocity, η the shear viscosity, and is the relativistic shear tensor (which is symmetric, traceless, and orthogonal to the flow velocity u). The fluid velocity and proper energy density satisfy T µν hydro u ν = −κε u µ . Moreover, the fluid/gravity gradient expansion yields expressions for all transport coefficients as functions of the proper energy density. For D = 4, the shear viscosity η = 1 4 κ(πT ) 3 , where the local temperature T is defined by ε = 3 4 κ(πT ) 4 [6]. Eq. (4.36) is precisely the constitutive relation of first order relativistic conformal hydrodynamics.
To compare our numerical results with the asymptotic predictions of the fluid/gravity correspondence, we first extract the fluid velocity u and rescaled proper energy density ε from the numerically computed stress-energy tensor (by finding the timelike eigenvector and associated eigenvalue of T µ ν , as discussed in section 2). With u and ε obtained via eq. (2.7), we then use eq. (4.36) to construct the hydrodynamic approximation to the spatial stress tensor, T ij hydro . Rotational symmetry in the transverse plane implies that all off-diagonal elements of the spatial stress tensor vanish, and that T xx = T yy . Therefore, we define a simple dimensionless residual function, where the average pressure p avg ≡ 2 3 T xx + 1 3 T zz . The residual R gives a measure of the relative deviation of the spatial stress from the prediction of the hydrodynamic constitutive relation (4.36). Fig. 10 plots R for collisions of both wide shocks (top) and narrow shocks (bottom). In each plot we exclude the region where R > 0.15. Specifically, for every value of z, we define t * (z) as the last time for which R(t, z) > 0.15 and exclude from the plot all points (t, z) for which t ≤ t * (z). We will denote by H the region where viscous hydrodynamics works at the 15% level or better (as measured by R). The dashed line in each plot is the curve with ∆t = 0.43 and τ hydro = 1.5 which, as seen in the figure, nicely approximates the boundary of region H. Fig. 10 clearly shows that our planar shock collisions result in the formation of R R t t z z Figure 10. The relative deviation R of the spatial stress tensor from prediction of first order viscous hydrodynamics for the case of wide shocks (top) and narrow shocks (bottom). As detailed in the text, we only display the region H = {(t, z) : R(t, z) ≤ 0.15} where the residual is no more than 0.15. The dashed curve, discussed in the text, is defined by eq. (4.39). For both cases, viscous hydrodynamics becomes a good description near mid rapidity when t 2.
an expanding volume of fluid which is well described by hydrodynamics everywhere except near the light cone, where non-hydrodynamic effects become important. At mid-rapidity, viscous hydrodynamics becomes a good description when t 2 [30]. As was noted in refs. [27,30], even in the region H where viscous hydrodynamics works at the 15% level or better, the first order viscous corrections are not small. The viscous stress tensor −2ησ µν in eq. (4.36) can be just as large as the zeroth order ideal fluid term. One manifestation of this is that in the local rest frame of the fluid (where u µ = δ µ 0 ), the spatial stress T local ij can be highly anisotropic with very different eigenvalues (i.e. pressures) in each direction. In the local fluid rest frame, this anisotropy is solely due to the gradient corrections in eq. (4.36). To illustrate this point, fig. 11 plots, for narrow shocks, the difference ∆p = T xx − T zz in the eigenvalues of the spatial stress at z = 0 (where by z → −z symmetry the fluid is at rest), normalized by the average pressure p avg . As just asserted, ∆p/p avg is O(1). Given the size of the first order gradient corrections, it is quite remarkable that the hydrodynamic constitutive relation works so well.
It is also illuminating to examine how well boost invariant flow approximates our numerical results. As the name suggests, boost invariant flow is defined by the condition that the system be invariant under arbitrary boosts in the longitudinal direction. Our initial conditions corresponding to two colliding shocks with non-zero widths are not boost invariant, and hence neither is the debris produced by the collision. Nevertheless, in a qualified sense which we make precise below, the produced debris does display some characteristics of nearly boost invariant flow. In what follows we focus on the case of narrow shock collisions, and on the dynamics in the region H, shown in fig. 10, where viscous hydrodynamics is applicable at the 15% level.
From the fluid/gravity correspondence, the fluid velocity and proper energy density (rescaled by κ) for boost invariant flow, up to second order in gradients, are given by [40] u µ dx µ = dτ ≡ cosh y dt + sinh y dz , (4.40a) where τ ≡ √ t 2 − z 2 is proper time, y ≡ tanh −1 z t is rapidity, and The energy scale Λ is set by initial conditions and is otherwise arbitrary. Each subsequent gradient correction to the proper energy density is suppressed by an additional power of (Λτ ) −2/3 ; for boost invariant flow, the fluid/gravity gradient expansion is precisely a late time expansion in inverse powers of proper time.
Our first comparison to boost invariant flow is shown in fig. 12, where we plot the longitudinal component u z of the fluid velocity at time t = 9 for the narrow shock collision. Also shown in the plot is the boost invariant flow result u z = sinh y = z/τ . Again, we display u z only in the region H where viscous hydrodynamics works at the 15% level or better. As is evident from the figure, the numerical result agrees quite nicely with this prediction of boost invariant flow. Fig. 13 shows a contour plot of the proper energy density ε extracted from our numerical  results and multiplied by a factor of τ 4/3 . Lines through the origin corresponds to events with fixed rapidity, t = z coth y. Inspecting eq. (4.40b), it is evident that if the flow was truly boost invariant then ετ 4/3 would asymptote to a constant, independent of rapidity, in the τ → ∞ limit. Fig. 13 shows that this is not at all the case; the flow is not globally boost invariant (as was also found in ref. [31]). However, one striking feature of fig. 13 is that contours of ετ 4/3 , at late times, are approximately straight lines through the origin, t ≈ z coth(y). This observation suggests that on each slice of constant rapidity y, the proper energy density is approximately given by eq. (4.40b) but with a rapidity dependent scale parameter, Λ = Λ(y). To test this hypothesis, on each slice of constant t/z = coth y we fit the proper energy density ε to the boost invariant expression (4.40b) allowing Λ to depend on y. In the left panel of fig. 14 we plot ε at y = 0, 0.85, 1.25, and 1.6, and the corresponding fit to eq. (4.40b). The agreement with eq. (4.40b) is remarkable. In the right panel of fig. 14 we plot the resulting scale parameter Λ(y) emerging from this fit to local (in rapidity) boost invariant flow.
It would be interesting to study more carefully the dependence of Λ(y) on the width of the incoming shocks, and to evolve longer in time in order to examine the asymptotic behavior of Λ(y) at large rapidity.

Motivation
Turbulent flows in relativistic boundary conformal field theories with ν spatial dimensions should be dual, via holography, to dynamical black hole solutions in asymptotically AdS ν+2 spacetime. This connection raises many interesting questions in gravitational physics. For example, what distinguishes turbulent black holes from non-turbulent ones? And what is the gravitational origin of the Kolmogorov scaling and energy cascades observed in turbulent fluid flows?
Gravitational dynamics may also provide insight into turbulence, in particular for problems where microscopic physics plays a crucial role in turbulent evolution. For superfluids, whose turbulent evolution is not governed by ordinary hydrodynamics, holography has already yielded insight into two dimensional turbulent flows [69]. In particular, ref. [69] found that turbulence in a two dimensional holographic superfluid exhibits a direct energy cascade into the UV. This stands in stark contrast to turbulence in normal fluids in two spatial dimensions, where enstrophy conservation gives rise to an inverse cascade to the IR. A fully consistent microscopic description of normal turbulence in three dimensions may also prove useful. Turbulence in three spatial dimensions is characterized by a cascade of energy from the IR to the UV, with dissipation occurring at microscopic length scales which may lie outside the hydrodynamic regime governed by the Navier-Stokes equation. Via holography, black hole solutions to Einstein's equations provide a laboratory in which one can study the domain of validity and late-time regularity of turbulent solutions to the Navier-Stokes equation.
In this section, we numerically construct black hole solutions in asymptotically AdS 4 spacetime dual to ν = 2 turbulent flows, where energy flows from the UV to the IR in an inverse cascade. The following discussion summarizes work first presented in ref. [33].

Setup
The boundary dimension D = 3. We choose an explicit parameterization of the rescaled spatial metricĝ ij that manifestly satisfies detĝ = 1, From the series expansions (3.41) we see that B and C have the near-boundary asymptotics We choose to replace the field redefinitions (3.53) involving the spatial metric with the following field redefinitions for the spatial metric functions The asymptotic behavior (4.43) implies that b =ḃ = c =ċ = 0 at the AdS boundary u = 0. Therefore, when solving eq. (3.37) forḃ andċ, we impose the Dirichlet boundary conditionṡ b =ċ = 0 at u = 0. We choose initial conditions corresponding to a locally boosted black brane. With our metric ansatz, ds 2 = r 2 g µν (x, r) dx µ dx ν + 2 dt dr , (4.45) a boosted black brane geometry is described by where u µ (x) is the boost velocity and r h (x) ≡ 4πT (x)/3, with T (x) the local temperature of the brane. (We are using simple Cartesian boundary coordinates for the boundary geometry.) The function R(x, r) satisfies where u 2 ≡ u i u i . For constant values of u µ and T , the metric (4.46) is an exact solution to Einstein's equations. After applying the time-space split (3.13) to g µν , the initial data for integrating Einstein's equations consists of the rescaled spatial metric with unit determinant, together with the asymptotic coefficients describing the energy and momentum density on the initial slice, with all functions evaluated at the initial time t i ≡ 0. We apply the above setup to the specific case of a boost velocity with sinusoidal variations plus small random perturbations (which serve to break the symmetry of the initial conditions), We study evolution in a periodic square spatial box of size L s and choose the wavevector Q = 10π/L s . The small fluctuations δu i are chosen to be a sum of the first four spatial Fourier modes with random coefficients, with the overall amplitude of the fluctuation adjusted to make |δu i (x)| max = 1/5. These initial conditions are unstable and capable of producing subsequent turbulent evolution if the Reynolds number Re is sufficiently large. For our initial conditions Re ∼ L s T . We choose box size L s = 1500 and the initial temperature 4πT /3 = 1.
A linear combination of the first 20 Chebyshev polynomials is used to represent the radial dependence of all functions, while an expansion of 305 plane waves (in each direction) is used to represent the spatial dependence. The discretized geometry was evolved from the initial time t i ≡ 0 to a final time t f ≡ 3001, using AB3 with timestep ∆t = 1/25. Computations were performed on a single six core Intel i7-3960x processor overclocked to 4.25GHz. With this relatively limited computing resource, producing the following results required approximately three weeks of running time.

Results
To illustrate the turbulent flow which emerges from the solution to Einstein's equations, we plot in fig. 15 the boundary vorticity, at six different times. We extract the fluid velocity u µ from the boundary stress tensor T µν via eq. (2.7), just as we did for the shock collisions in Section 4.2.
At time t = 0, when the fluid velocity is given by eq. (4.49), the vorticity is approximately sinusoidal in the x 1 direction and translationally invariant in the x 2 direction. By time t = 752, an instability is visible and the approximate symmetry of the initial conditions is destroyed. By time t = 1248, the instability has generated many small vortices with fluid rotating clockwise (red) and counterclockwise (blue). Subsequently, vortices with the same rotation tend to merge together producing larger and larger vortices, as seen in the evolution snapshots at times t = 1760, 2192, and 3001. As time progresses, the number of vortices decreases while the typical vortex size grows. This is a characteristic signature of an inverse cascade.
It is instructive to compare the gravitational evolution with predictions from the Kolmogorov theory of turbulence. A simple quantity to study is the power spectrum of the fluid velocity, defined as (4.52) A celebrated result of Kolmogorov is that for driven steady-state turbulence the power spectrum P obeys the scaling P(t, k) ∼ k −5/3 , (4.53) within an inertial range k ∈ (Λ − , Λ + ). The lower limit Λ − is determined by the size of the largest eddies in the system, while the upper limit Λ + is set by the scale on which viscous effects damps small eddies. Despite the fact that our system is not driven or in a steady-state configuration, we do see hints of Kolmogorov scaling. In fig. 16 we plot P at time t = 1008. Our numerical results are consistent with the scaling (4.53) in the inertial range k ∈ (0.025, 0.055). As we are not driving the system, evidence of the k −5/3 scaling is transient and is destroyed first in the UV, with the UV knee at k = 0.055 shifting to the IR as time progresses further. Beyond the inertial range the spectrum decreases like P ∼ k −p with p ∼ 5 until k ∼ 0.15 beyond which P decreases exponentially.
The inverse cascade also manifests itself in bulk gravitational quantities. One interesting quantity to consider is the horizon area element √ γ. In our coordinate system, and in the limit Figure 15. The boundary vorticity at six different times. The initial conditions shown at time t = 0 give rise to an instability which produces many vortices as seen in the subsequent evolution at times t ≥ 1248. Vortices colored red (blue) correspond to clockwise (counterclockwise) fluid rotation. As time progresses, vortices of like rotation tend to combine to produce larger and larger vortices.  of large Reynolds number Re 1, the event and apparent horizons approximately coincide at r = 1 and the horizon area element is √ γ ≈ √ −g r=1 . 45 In fig. 17 we plot √ γ for the same sequence of times displayed in fig. 15. The evolution of √ γ closely mirrors the evolution of the vorticity on the boundary shown in fig. 15. At time t = 0 when the fluid velocity is given by eq. can be constructed from the null normal n M to the horizon and an auxiliary null vector M whose normalization is conveniently chosen to satisfy M n M = −1. The extrinsic curvature is then given by where the projection operator Π M N ≡ δ M N + M n N . Since the horizon is at r ≈ 1 we choose n M dx M = dr and M dx M = −dt. In our coordinate system, the horizon curvature satisfies Θ M N Θ N M = Θ i j Θ j i , where i, j run only over the spatial coordinates. For later convenience we define the rescaled traceless horizon curvature θ i j ≡ (γ/κ 2 ) 1/4 Σ i j , where Σ i j ≡ Θ i j − 1 ν Θ k k δ i j is the traceless part of the extrinsic curvature, and κ is the eigenvalue of the geodesic equation, n M ∇ M n Q = λ n Q . (4.55) We define the horizon curvature power spectrum with θ i j ≡ d ν x θ i j e −ik·x , and plot the ratio A(t, k)/P(t, k) in fig. 18. As this figure makes clear, our numerical results are consistent with the simple scaling relation A(t, k) ∼ k 2 P(t, k) . (4.57) Evidently, these horizon and boundary observables are highly correlated. This follows directly from the applicability of the fluid/gravity correspondence.
Both qualitative and quantitative features of our numerical results can be understood in terms of ideal conformal hydrodynamics and the fluid/gravity correspondence. As discussed in sec. 4.2, in the limit of long wavelength spatial fluctuations (compared to 1/T ) Einstein's equations can be solved perturbatively with a gradient expansion. At leading order the metric is precisely the locally boosted black brane (4.46), with the evolution of the fluid velocity u and temperature T governed by relativistic ideal conformal hydrodynamics [11,66]. In other words, to leading order in the gradient expansion, solutions to Einstein's equations can be generated merely by solving the equations of relativistic ideal conformal hydrodynamics and constructing the bulk metric from the resulting fluid velocity and temperature via eq. (4.46).
It was recently demonstrated that turbulent evolution in two dimensional ideal relativistic conformal hydrodynamics gives rise to an inverse cascade and Kolmogorov scaling (4.53) [67]. Since the relativistic hydrodynamic equations reduce to the non-relativistic incompressible Navier-Stokes equation at low velocities [65], this connects directly to classic results on non-relativistic two dimensional turbulence. It is well known that two dimensional non-relativistic incompressible turbulent flows exhibit Kolmogorov scaling and an inverse cascade, with the latter a consequence of conservation of enstrophy (the square of the vorticity). As demonstrated in ref. [67], the equations of two dimensional ideal relativistic conformal hydrodynamics conserve a relativistic generalization of enstrophy.
We find that our numerical metric is surprisingly well approximated by the boosted black brane metric (4.46). To perform the comparison, we extract the flow field u and the proper energy density ε from T µν via eq. (2.7). The proper energy density is converted to a local temperature via the (static AdS 4 black brane) relation T ≡ 3 4π 3 2 ε 1/3 . The flow field u and local temperature T are then used to construct the boosted black brane metric (4.46). Finally, we compute the difference ∆g µν between the numerical metric and the boosted black brane metric and define the error to be max{|∆g µν |} on a given timeslice t. As shown in fig. 19, the boosted black brane metric ansatz (4.46) approximates the complete geometry, even at early times, to better than 1%! Although the accuracy with which the simple boosted black brane ansatz approximates the numerical solution is remarkable, it should not be too surprising that turbulent evolution in two spatial dimensions gives rise to dual geometries which are reasonably well approximated by the locally boosted black brane ansatz. First of all, irrespective of the dimensionality, turbulent flows require large Reynolds number, Re 1, which (in a strongly coupled fluid) is equivalent to small gradients compared to the local temperature T . This is precisely the regime where the fluid/gravity gradient expansion should be well behaved. Second, the inverse cascade of turbulence in two spatial dimensions implies that gradients become smaller and smaller as energy cascades from the UV to the IR. Therefore, the leading term (4.46) should become a better and better approximation to the metric as time progresses and the inverse cascade develops.
At least for ν = 2, the above observation has powerful consequences for studying turbulent black holes. Instead of numerically solving the equations of general relativity, one can simply  Figure 19. Absolute deviation, as a function of time, between the numerically computed spacetime metric and a boosted black brane metric, with fluid flow and local temperature extracted from the numerical solution. The upper (blue) curve shows the maximum size of metric components on the timeslice t, and the much lower (green) curve shows the maximum, on the given timeslice, of the difference between the numerically computed metric and the boosted black brane ansatz. study the equations of hydrodynamics and construct the bulk geometry via the fluid/gravity gradient expansion. This is particularly illuminating in the limit of non-relativistic fluid velocities |u| 1, where the bulk geometry and boundary stress are asymptotically close to equilibrium. As shown in ref. [65], under the rescalings t → t/s 2 , x → x/s, u → s u, and δT → s 2 δT (with δT the variation in the temperature away from equilibrium), as s → 0 the boundary evolution of δT and u implied by the fluid/gravity correspondence reduces to the non-relativistic incompressible Navier-Stokes equation. Indeed, the above rescalings are symmetries of the Navier-Stokes equation. Likewise, in the s → 0 limit the geometry dual to the Navier-Stokes equation can be computed analytically [65]. At least for two spatial dimensions, where is it known that solutions to the Navier-Stokes equation remain regular, it should be possible to (re)derive results from classic studies of turbulence, such as Kolmogorov scaling (4.53), directly from the dual gravitational dynamics. This is discussed in more detail in ref. [33].

Conclusions
We have presented a characteristic formulation of gravitational dynamics which permits accurate and efficient study of a wide variety of gravitational initial value problems in asymptotically anti-de Sitter spacetimes. The requirement of the approach that geometries of interest have an apparent horizon cloaking any caustics in the infalling null congruence has, in practice, not been a limitation. Problems with numerical stability are less severe than is often the case with numerical relativity, due to helpful attributes of our characteristic formulation, the presence of an apparent horizon, and the asymptotic anti-de Sitter geometry. With only modest computing resources, we have shown that problems whose symmetries reduce the dynamics to 1+1 dimensional partial differential equations (homogeneous isotropization), 2+1 dimensional PDEs (colliding planar shocks), or 3+1 dimensional PDEs (turbulence in two space dimensions), are quite manageable. An obvious question concerns the feasibility of solving 4+1 dimensional gravitational dynamics with no simplifying symmetry restrictions. We are optimistic that various problems in this category, such as studying turbulent fluids in three spatial dimensions, or off-center "heavy ion" collisions, will also be feasible.
These combinations all transform as scalars with respect to radial shifts (3.16).
Einstein's equations, now in the presence of bulk sources, decompose into the scalar equations: 0 = tr G − 1 2 G 2 + 2T 00 , (A.8) 0 = A + 1 2 ∇ · F + 1 2 F · F + 1 2 (tr d + G) + 1 4 tr (G d + G) + 2 ν Λ − 1 ν tr s + (1− 2 ν ) κ , (A.9) 0 = tr [d + (d + G) − A (d + G) − 1 2 (d + G) 2 ] + 2 ∇ · E + 1 2 tr (Ω 2 ) + 2τ , (A.10) two vector equations: 12) and the symmetric tensor equation: Recall that Ω ij (3.25) and E i (3.26) are the "magnetic" and "electric" parts of the radial shift field strength. The trace of the last equation separates from the traceless part, as before, and reads 0 = G 1/2 tr (d + G) G −1/2 − R + 2Λ + ∇ · F + 1 2 F · F − 2κ . (A.14) The simple integration strategy described in section 3.7 relies on the nesting of the equations. To remain applicable, the Σ equation (A.8) must only require knowledge ofĝ ij , the F equation (A.11) must only require knowledge ofĝ ij and Σ, the d + Σ equation (A.14) must only depend onĝ ij , Σ and F , and the d +ĝij equation (A.13) must not depend on A. The radial shift invariance of the linear combinations (A.5)-(A.7) guarantees that the explicit factors of A and F i which appear in these expressions cancel when combined with the metric dependence inside T M N . For either an electromagnetic field, or a scalar field φ with arbitrary potential V (φ), one may easily confirm that the bulk source terms do not upset the nesting of equations which underlie the integration strategy. 46 The apparent horizon condition (3.44) is not affected by the addition of a bulk stressenergy tensor, but the horizon stationarity condition (3.47) receives modifications from source terms (due to the use of Einstein's equations in its derivation) and becomes 0 = ∇ 2 A − ∇A · (F − G F )

B Riemann tensor components
For some purposes, such as evaluating curvature invariants, it is desirable to have explicit expressions for the Riemann tensor components generated by our metric ansatz (3.14). Defining components with respect to the frame (A.1) is convenient, as this makes the results transform as scalars with respect to radial shifts. (Moreover, the corresponding components of the metric (A.3) are especially simple.) One finds: with R ijkl defined in footnote 17.

C Spatially covariant expressions
Using our metric ansatz (3.14), explicit forms of Einstein's equations and Riemann curvature components are most compact when written using the modified spatial derivatives (3.21) which are covariant under both spatial diffeomorphisms and radial shifts, as done in section 3.4 and appendices A and B. Neverthess, there may be occasions where it is helpful to have available equivalent expressions written using ordinary spatial covariant derivatives. These are recorded below, using the decomposition (A.4)-(A.7) of any bulk stress-energy tensor. Einstein's equations may be separated into three scalar equations, 0 = tr G − 1 2 G 2 + 2T 00 , (C.1) two vector equations, and the symmetric tensor equation, (C.6) The trace of this last equation separates from the traceless part and reads 47 0 = √ G tr (d + G) + ∇ · F + 1 2 tr (G )F · F / √ G + 1 2 ∇ · (tr (G ) F ) − R (ν) + 2Λ In the above, R ij and R (ν) denote the spatial Ricci tensor and Ricci scalar, respectively, and