Collision of localized shocks in AdS$_5$ as a series expansion in transverse gradients

We introduce a computational framework to more efficiently calculate the collision of localized shocks in five dimensional asymptotically Anti-de Sitter space. We expand the Einstein equations in transverse gradients and find that our numerical results agree well with exact solutions already at first order in the expansion. Moreover, the Einstein equations at first order in transverse gradients can be decoupled into two sets of differential equations. The bulk fields of one of these sets has only a negligible contribution to boundary observables, such that the computation on each time slice can be simplified to the solution of several planar shockwave equations plus four further differential equations for each transverse plane `pixel'. At the cost of errors of $\lesssim 10 \%$ at the hydrodynamization time and for low to mid rapidities, useful numerical solutions can be sped up by roughly one order of magnitude.


Introduction
In the past decade a considerable progress has been achieved with respect to simulating heavy ion collisions via holography [1][2][3][9][10][11][12][13][14][15][16][17]. The difficulty of simulating the first few fractions of 1 fm/c after a heavy ion collision comes from the strongly coupled and far from equilibrium nature of the incipient quark-gluon plasma. The gauge/gravity duality, or the dynamical equivalence of strongly coupled, conformal, supersymmetric field theories in d dimensions to weakly coupled Einstein gravity in d + 1 dimensions, provides a powerful tool for studying the early stages of a heavy ion collision. The duality between gravity in five dimensional asymptotically Anti-de Sitter (AdS) space and N = 4 super-Yang-Mills (SYM) theory is especially useful in this context, as the latter theory provides a reasonable approximation to a QCD plasma at high temperatures.
Holographic simulations of heavy ion collisions have become progressively more intricate and realistic up to the point of the exact treatment of collisions of localized nuclei in [3]. These studies revealed non-trivial and sometimes unexpected properties of the pre-hydrodynamic quark gluon plasma, from a sizable transverse flow [3], to universal behaviors regarding post collision flow and hydrodynamization at near constant proper time [1,9,10], just to name a few. However, there are still many interesting and demanding problems that remain unexplored. In particular the simulation of localized heavy ion collisions with a granular structure needed to explain the observed large event-by-event fluctuations of flow observables. In this work we develop tools to substantially facilitate the treatment of these problems by introducing a computational framework to efficiently calculate localized shockwave collisions in five dimensional asymptotically AdS space.
Heavy ions collided at LHC and RHIC are highly relativistic and Lorentz contracted. Due to this strong contraction, gradients transverse to the beam axis are very small compared to longitudinal gradients. The idea of approximating heavy ion collisions by neglecting transverse gradients is not new and is the justification for earlier studies of planar shockwave collisions in Anti-de Sitter space [1,2]. Planar collisions can be thought of as the zeroth order solutions to Einstein equations that were expanded in transverse gradients.
In this work we study the next-to-leading order results when the Einstein equations are expanded in transverse derivatives. Already at first order in transverse derivatives we find remarkably good agreement with the exact results in [3]. Moreover, we show that at first order the Einstein equations can be decoupled into two sets of differential equations for two disjoint sets of bulk fields. The fields of the larger set contribute negligibly to boundary observables. We explain why this happens and find that we can reproduce the exact results to good accuracy for times up to the hydrodynamization time, even for broad shocks, with an algorithm that only requires the planar shock solution and four additional differential equations to be solved on each time slice and transverse 'pixel'. Just using a Mathematica implementation, we can reproduce, to good accuracy, the results of [3], using 12 cores on a standard desktop machine, in about 36 hours.
The paper is structured as follows: For those readers less interested in computational details, we present key results at the beginning, in section 2. This is followed by a broader discussion of the transverse gradient expansion in section 3 and more detailed discussion of our calculations in section 4. Section 5 contains further analysis of our numerical results, followed by a few concluding remarks in section 6.

Approximating transverse flow and hydrodynamization time in heavy ion collisions
From [3] we know that during the early, pre-hydrodynamic stage after the collision of two relativistic heavy ions, a sizable transverse flow is produced, which simulations at zeroth order in transverse derivatives, e.g. [1][2][3][9][10][11][12][13][14][15][16][17], fail to capture. While there are models, e.g. [8], which roughly estimate the size of transverse momenta from zeroth order data, the simplifications and assumptions required do not yield quantitative control or reliable error estimates. The systematic approximation discussed in this work requires the aspect ratios of projectiles to be large, an assumption which is very well fulfilled for realistic heavy ion collisions. We find that already for localized collisions of projectiles with a comparatively small aspect ratio of 8, the transverse flow of our approximation and the exact results in [3] agree surprisingly well, both at mid and small rapidity, up until the hydrodynamization proper time. Examples of our results for transverse energy flux are depicted in Fig. 1. The error of the first order approximation increases slowly for larger rapidities and later proper times. As suggested by the approximation in [8] and consistent with [3] we find that the maximum value of the transverse energy flux as a function of proper time reaches a plateau as displayed in Fig. 2 and stays approximately constant. For higher rapidities the plateau is reached at later proper times.
One key observable to extract from holographic heavy ion collisions is the hydrodynamization time. To quantify how well the system is described by hydrodynamics one introduces a residual ∆ (defined later on in Eq. (5.4)), which measures the deviation of the stress energy tensor from its hydrodynamic approximation. For a small ∆, typically below 0.15, the hydrodynamic approximation is appropriate. Remarkably we find that our approximate solution and the exact solution of [3] for ∆ match well already at zeroth order in transverse derivatives, especially for low rapidities and in the central region of the collision. We also find approximately the same hydrodynamization time (i.e., the time at which the residual ∆ drops below a threshold of 0.15) in the central region. We display the residual ∆ and its zeroth order in transverse derivatives approximation in Fig. 3 at various proper times and rapidities.
In general, observables evaluated before the hydrodynamization time are well described by either the zeroth or first order transverse derivative approximation.  Figure 1: The averaged, lab-frame transverse momentum density T 0⊥ = (x ⊥ ) i T 0i up to first order in transverse derivatives, as a function of the transverse plane radius x ⊥ = x 2 + y 2 , at proper time τ = 1.25 (left) and τ = 2 (right). Both plots display the averaged transverse energy flux at rapidity ξ = 0 (blue curve) and at ξ = 1 (yellow curve). The red dashed curves represent the exact results of [3]. At τ = 1.25 and low rapidities the first order approximation and the exact results agree very well. This agreement slightly deteriorates for larger rapidities and later proper times, but still remains quite good. We work in units in which the maximum of the longitudinally integrated energy density is set to 1.

Transverse gradient expansion
Let G µν be a bulk solution to Einstein equations in asymptotically AdS space which varies slowly with respect to the transverse boundary coordinates x ⊥ = (x, y) compared to its variation with respect to the boundary coordinate x || = z. We introduce a formal parameter such that, after definingx ⊥ ≡ x ⊥ , gradients alongx ⊥ and x || are of comparable order Since the metric G µν is assumed to only weakly depend on x ⊥ the parameter is a small number and we can treat the Einstein equations written in terms of the rescaled coordinatex ⊥ as a perturbative expansion in . We express the metric in terms of the rescaled coordinate G µν (x 0 , x || , x ⊥ ) =G µν (x 0 , x || ,x ⊥ ) expand the Einstein equations in powers of , truncate at a given order, and then scale back to the original coordinates x ⊥ when solving the equations. This procedure is equivalent to replacing ∂ ⊥ → ∂ ⊥ , expanding the Einstein equations in powers of , solving them order by order and setting = 1 in the end.
To spell this out more explicitly, we write the Einstein equations for a metric G Max <T 0⟂ > ξ=0 ξ=1 Figure 2: The maximum value of the averaged, lab-frame transverse momentum density T 0⊥ = (x ⊥ ) i T 0i up to first order in transverse derivatives, as a function of proper time τ , for two rapidities ξ = 0 (blue curve) and ξ = 1 (yellow curve). The energy flux appears to reach a stable plateau, consistent with the approximation in [8].
schematically as Expanding in transverse derivatives, we have where the differential operator E (i) contains i powers of transverse derivatives. Let G (i) µν denote an approximate solution to the Einstein equations valid to order O( i ) so that We assume that G (i) µν is known up to some maximum order I on some initial Cauchy surface at time t = t 0 . We want to solve for G (i) µν up to this order I for some future time interval t 0 ≤ t ≤ t 1 . At order zero the metric G ) indicating when and where the hydrodynamic approximation is a reasonable description. For ∆ 0.15 hydrodynamics is applicable. We show the results for the zeroth order in transverse derivative approximation (first row) and the exact solutions (second row) from [3], for collisions with otherwise identical parameters. Especially for low rapidities and close to the central region we see quite good agreement between the exact and approximate solution, already at zeroth order in transverse derivatives. In both cases one finds a hydrodynamization time of t ≈ 1.25 in the central region and for zero rapidity. Visible artifacts are a consequence of the sparse grids used to discretize the transverse directions.
systematically corrects this zeroth order approximation by writing and demands that the Einstein equations hold up to the next order. Let ∆ (i) L be the planar Lichnerowicz operator evaluated on G (i) , Then Eq. (3.4) will be satisfied if

Transverse derivative corrections to shock collisions
In the next two subsections we quickly review the general treatment of (localized) shock collisions in AdS 5 [1][2][3], after which we discuss in more detail the solution of the Einstein equations in the transverse gradient expansion up to first order. We address how to construct consistent initial data up to first order in transverse derivatives in section 4.3.

Single localized shocks in Fefferman-Graham coordinates
As in [1,3], a single shock moving in ±z direction in Fefferman-Graham (FG) coordinates may be described analytically, The metric ansatz (4.1) solves the Einstein equations provided the function h ± fulfills the linear differential equation 2) through first order in transverse derivatives. To match previous work, we will use a Gaussian profile With this choice, the single shock metric is a valid zeroth order solution solving (3.5), which we denote as G (0) [h ± ]. Moreover the same metric is a valid solution at first order, so G In general the metric (4.1) corresponds to a state in the dual field theory for which

Integration strategy
For the time evolution of two colliding localized shocks, we transform to infalling Eddington-Finkelstein (EF) coordinates and employ the characteristic formulation of general relativity. In EF coordinates the line element of the metric G µν may be written as Here x ≡ (t, x ⊥ , x || ). In these coordinates the boundary is positioned at r = ∞.
Following [1], we have used diffeomorphism invariance to fix the component of the line element proportional to drdt to equal 2 drdt. The metric (4.6) is also invariant under shifts of the radial coordinate of the form As in [1] we will exploit the radial shift symmetry to fix the radial position of the apparent horizon r h to correspond to the end of our integration domain. Let us write the metric in EF coordinates more specifically as with i = x, y, z denoting spatial derivatives. Let G ij = Σ 2ĝ ij with detĝ ij = 1. Note that, unlike in section 3, henceforth G ij denotes the spatial part of the EF-metric. For convenience we also slightly modify our convention regarding the superscript index (·) (i) . Henceforth, for a function X, let X (i) denote the difference between the i-th and (i − 1)-th order solution. In other words, using the notation of section 3, we abbreviate δX (i) with X (i) and we will henceforth use the notation X (i) only in this context. Through first order in transverse derivative expansion the unimodular condition onĝ implies that detĝ (1) ) ij = 0 for the first order correction in transverse derivativesĝ (1) . The spatial scalar factor Σ has the near boundary form (4.9) and the other metric functions have a near boundary beahvior given by Let a 4 , f 4 i andĝ 4 ij be the O(r −4 ) expansion coefficients of the near boundary expansion of the rescaled functions (r + λ) −2 A, (r + λ) −2 F i and (r + λ) −2ĝ ij . The CFT stress energy tensor in 3 + 1 dimensions is determined by the O(r −4 ) coefficients a 4 , f 4 i and ĝ 4 ij in the near boundary expansions (4.10)-(4.12): As demonstrated in [5,6] and generalized in [1] the Einstein equations expressed in EF coordinates (4.6) can be written as a nested system of ordinary radial differential equations (ODEs) on each time slice. Defining derivatives along outgoing geodesics, The differential equation that determines the shift λ, so that the horizon matches the end of our integration domain, is obtained by requiring that the congruence of outgoing geodesics has vanishing expansion rate there. This results into the non-linear ODE where r h is the radial position of the apparent horizon, primes denote radial derivatives ∂ r and ∇ is the covariant derivative acting on spatial tensor fields using the Christoffel connection of G ij (see (A.3)). Requiring this expansion rate to be constant at the radial position of the horizon leads to a linear elliptic differential equation for ∂ t λ. Having solved for {∂ t a 4 , ∂ t f 4 i , ∂ tĝij , ∂ t λ} on one time slice, we use the fourth order Runge Kutta method to evolve the metric in time.
In the case of exact localized shock collisions [3] the elliptic differential equation for ∂ t λ is difficult to handle and solving it requires a large amount of wall clock time and memory 1 . The treatment of this equation is substantially simplified if we expand in transverse derivatives. The details of this differential equation are not relevant for the following discussion, for the interested reader, the exact differential equation is given in (3.47) of [1]. Let us write this linear horizon fixing equation schematically as where K ij , L i , M and S are functions of the metric and its derivatives. At zeroth order in transverse derivatives Eq. (4.19) simplifies to an ordinary linear differential equation in the longitudinal coordinate z, which can be trivially parallelized over the transverse coordinates x ⊥ , 20) To compute the first order in transverse derivative correction we can reuse the Green's function for D (0) and transverse derivatives only contribute to the source termS (1) , Eq. (4.22) is clearly still parallelizable over every value of x ⊥ . The ability to parallelize and reuse the Green's function for D (0) allows us to circumvent this bottleneck of the exact treatment. The same principle applies to the Einstein equations written as a system of ODEs that we solve on each time slice. Again we can reuse the radial Green's functions computed at zeroth order in transverse gradients. The main advantage of an expansion in transverse gradients is thus that Green's functions for each point in the transverse plain are no more expensive to compute than in the planar case, both for the nested system of radial ODEs and the elliptic differential equation which fixes the horizon position. We solve the system of ODEs and (4.19)-(4.23) using spectral methods [7]. The details of our numerics can be found at the beginning of section 5.

Initial data up to order O(∇ ⊥ )
We now discuss the construction of a numerical solution for the initial data of localized shocks, in five dimensional AdS space, up to first order in transverse gradients ∇ ⊥ .
force approach to solving this differential equation would require the inversion of a matrix of size (N x N y N z ) × (N x N y N z ). For modest spatial grid sizes of, e.g., 32 × 32 × 512 this requires over 2 TB of memory and more wall clock time than solving all other equations on a given time slice. More sophisticated methods involving domain decomposition, require very large amounts of memory and run time as well.
Following [1] we transform the single shocks (4.1) from FG coordinates to infalling Eddington-Finkelstein (EF) coordinates. In EF coordinates the metric takes the form (4.6). In this coordinate system the family of curves γ x µ 0 (r) = (x µ 0 , r) with varying radial coordinate r is a geodesic congruence with affine parameter r. Demanding that the tangential vector field Y µ of the image of this congruence under the coordinate transformation to FG coordinates solves the geodesic equation, Y µ ∇ µ Y ν = 0, provides a system of differential equations whose solution allows us to explicitly express FG coordinates Y (X) in terms of EF coordinates X. As boundary conditions of the geodesic equations we require that the coordinate systems match on the boundary, lim ρ→0 Y = lim r→∞ X, and that the radial coordinates match at a hypersurface in the bulk where ρ max = 1 r min with r min the end of our integration domain, r ∈ [r min , ∞]. We solve the geodesic equations up to first order in transverse derivatives, which contribute via derivatives acting on the metric in the Christoffel symbols 2 . We write the zeroth order ansatz plus first order corrections for the FG coordinates Y (X) in terms of the EF coordinates X µ = (t, x ⊥ , z, u) with inverted radial coordinate u = 1/r as With (4.24)-(4.28) we solve Here and henceforth we use the usual subscript notation for derivatives (·) ,x = ∂ x (·) for a coordinate x . Let G ab denote the spatial part of the metric in EF coordinates and 2 Even though it might appear so at first glance, writing the geodesic equation does not introduce an ambiguity regarding the transverse derivative expansion, since the relation Y µ ∂µY ν = ∂ 2 τ x ν holds up to every order in . Thus, the most convenient way to expand the geodesic equation in transverse derivatives is to expand the equation It is then easy to see why transverse derivatives only explicitly contribute via the Christoffel symbols. g αβ be the metric in FG coordinates. Then For the numerical time evolution we will need the zeroth order terms (ĝ ab ) ± and first order corrections (ĝ ab ) (1) ± of the rescaled spatial metricĝ. These are given by . The fourth order near boundary expansion coefficients of A and F i in (4.8) for a single Gaussian shock moving in ±z direction are Comparing the near boundary expansion of the determinant of the spatial metric given in (4.31), (4.32)-(4.34) expressed by the functions α, β, γ, δ and their derivatives, with the near boundary expansion of Σ, shows that the shift parameter λ, introduced in (4.7), can be obtained from the relation (λ) Finally the initial data for the collision of localized shocks at order O( i ) in EF coordinates at t = t 0 = −2 is given by combining the data for well-separated counterpropagating single shocks, Since each shock individually is an exact solution to Einstein equations, their sum is an approximate solution with exponentially small errors, if they are sufficiently separated, such that they don't overlap and the region between them is pure AdS up to exponentially small errors. The initial time t 0 = −2 is chosen such that left and right moving shocks are spatially well separated and have no overlap in the bulk region between the apparent horizon and the boundary. As in [1][2][3] we work with a small background energy density to deal with irregularities in the integration domain.

Hierarchy among O(∇ ⊥ )-fields
The Einstein equations, expanded in transverse derivatives through first order, separate into two sets of differential equations for two disjoint sets of bulk fields S 1 and S 0 : with corresponding sets of boundary initial data, Remarkably, the first order transverse derivative corrections to the zeroth order solutions of the fields in S 0 contribute negligibly to the expectation value of boundary observables. Their contributions are two orders of magnitude smaller than those of first order corrections to fields contained in set S 1 . As shown in Appendix B, through first order in transverse derivatives, the differential equations for the fields in S 0 and the functions in I 0 do not contain any explicit transverse derivatives. The only first order corrections to these differential equations are a consequence of negligibly small O( )-corrections to the initial conditions {(ĝ ij ) (1) , (a 4 ) (1) , (f 4 z ) (1) } on the first time slice, a result of small transverse gradients appearing in the coordinate transformation from FG to EF coordinates. Thus, up to an error of about 1% we can ignore their first order corrections to the boundary stress energy tensor. Since their equations also decouple from the differential equations for fields in S 1 and temporal differential equations of functions in I 1 , we can substantially simplify the calculation. In order to compute the relevant contributions to the boundary stress energy tensor up to the first order in transverse derivatives, it is thus sufficient to solve planar shockwave collisions plus only four additional differential equations on each time slice. Those four equations are given by the O( )-terms of the x, y components of Eq. (A.1b), and the {x, z} and {y, z} components of Eq. (A.1d) given in Appendix A. For each point in the transverse plain the numerics are, therefore, only slightly slower than planar shockwave collisions. It should be added that the approximation of expanding in transverse gradients improves as one further decreases the longitudinal thickness of the shocks or more generally increases the aspect ratio. The parameters in [3] were chosen to simulate the collision of rather broad blobs of energy; for collisions of more realistic, highly Lorentz-contracted projectiles, the agreement should, thus, further improve.
We have computed collisions with identical initial conditions comparing the full first order in transverse derivatives with the above simplified approach. We verified that the contributions to boundary observables in both treatments were identical up to errors of ≈ 1%.

Easy O(∇ 2
⊥ )-terms At late times, well after the hydrodynamization time, differences between the exact solution and the approximation by transverse derivative expansion, truncated at first order, become substantial (see e.g. Fig 4, 5, 12). At time t ≈ 4 higher order corrections in transverse derivatives can no longer be neglected. However, treating second order corrections exactly would more than double memory usage and substantially increase run time to a point where the advantage over the exact solution might become small. Solving radial differential equations and the elliptic differential equation, to update the shift parameter, continue to be much faster compared with the exact treatment. Computing the source terms of each radial differential equation will become the new run-time bottleneck when computing higher order corrections.
One can separate second order corrections in transverse derivatives to the boundary stress energy tensor into two categories: contributions from bulk fields and contributions from boundary data which is time evolved using the covariant conservation equation of the stress energy tensor. As explained above, calculating contributions to boundary observables coming from spatial components of the bulk metric is costly in memory and run time. On the other hand, one may easily determine the partial second order corrections in transverse derivatives to the energy density (∆a 4 ) (2) and the longitudinal momentum density (∆f 4 z ) (2) arising from the first order terms of the transverse energy flux T tx , T ty , the shear components T zx , T zy and the conservation equation of the boundary stress energy tensor (4.17) using plus the initial condition (∆f 4 z ) (2) = (∆a 4 ) (2) = 0 on the first time slice. Solving for this subset of second order corrections requires negligible additional memory and wall-clock time.
A priori, one might only expect to obtain a qualitative estimate of the size of second order corrections from such an incomplete calculation of second order terms. Hence, it is all the more remarkable that we find that inclusion of second order terms in the boundary conservation equation (4.46) and (4.47), alone, drastically improves the agreement between the transverse derivative expansion and the exact result for most observables (see Fig. 5, 11, 12). Nonetheless, this partial second order treatment should be employed with care, like every partial higher order calculation without proof one has captured a dominating contribution.

Results
We compute the coordinate transformation from FG to EF coordinates using a Chebyshev grid in the radial direction with three domains containing 31 grid points each and Fourier grids in longitudinal and transverse directions. The longitudinal grid consists of N = 256 and the transverse grids consist of N x,y = 36 grid points. For the time evolution in EF coordinates we found a radial Chebyshev grid with only two domains with 18 grid points each to be sufficient. After computing the EF initial data, to perform the time evolution we remap data to a smaller transverse Fourier grid with 12 grid points in each dimension and a coarser longitudinal Fourier grid with 140 grid points. To enable comparison of our approximation with the exact results our physical parameters are chosen to match those in [3]. The physical sizes of the spatial 3-torus in which we placed the shocks are L z = 12 in longitudinal direction and L x,y = 32 in transverse direction. We chose the longitudinal width parameters of the shocks to be w = 0.5 and R = 4 and chose a unit amplitude A = 1. The impact parameter is b = 3 the energy density a 4 , the momentum density f 4 i and the shift parameter λ (all of which we compute via the fourth order Runge-Kutta method), removing one third of all modes corresponding to the highest frequencies, while keeping the remaining two thirds. In addition we apply a radial filter on the longitudinally filtered updates by transforming the radial Chebyshev grid onto a slightly smaller, single domain grid with 35 grid points, and transforming back to the original grid. We do not filter during the individual substeps of the Runge-Kutta algorithm.
We compare the stress energy tensor after the collision with its hydrodynamic approximation, where the constitutive relations are truncated after the first order in derivatives. For this we solve the Eigenvalue equation order by order in transverse derivatives for the fluid velocity u µ and the proper energy density ε. The hydrodynamic approximation where the viscous stress is given by is also expanded up to order O(∇ ⊥ ). Here p is the pressure and η the shear viscosity. In order to evaluate how well the system is described by hydrodynamics we compute the residual ∆ = 3 ε ∆T µν ∆T µν (5.4) with ∆T µν = T µν − T µν hydro . Again let ∆ (i) denote the O( i ) contribution to ∆. The hydrodynamic approximation is generally taken to be valid for ∆ < 0.15. We show the results for the residual at zeroth order, ∆ (0) , and the first order correction ∆ (1) , at proper times τ = 1.25, τ = 2 and rapidities ξ = 0, ξ = 1 in Fig. 6 and Fig. 7, respectively. Evidently the inclusion of first order transverse dynamics only has a small effect on the residual ∆, and already at zeroth order in transverse derivatives we find in the central region (small x ⊥ ) the same hydrodynamization time 3 of t ≈ 1.25 as in [3].
Our results regarding transverse energy flux, shown earlier in Fig. 1 and 2, illustrate good agreement between the O( )-approximation and the exact solution regarding transverse flow. In Fig. 4 we display the xx and zz-components of the stress tensor, T xx and T zz , at the origin x = y = z = 0 as a function of time. Again both stress components agree well with the results in [3]. Only for T xx and at late times, t > 2, do we observe a noticeable disagreement. As shown in Fig. 5, this error is more than halved by including the partial O( 2 )-corrections discussed in section 4.5. This suggests that even for rather broad shocks, as shown in this work, and late times t ≈ 4 the full inclusion of O( 2 )-corrections should provide a decent approximation, whereas for earlier times, t < 2, first order corrections appear sufficient.
We computed the fluid 3-velocity v =ū/u 0 up to order O( ) and display its absolute value at time t = 4 on both a x = 0 slice and a z = 0 slice in Fig. 8. In Fig. 9 we display the absolute energy flux |T 0i | at various times on the y = 0 slice. In Fig.  10 we show the analogous plots for the energy density T 00 . A general comparison of O( )-results and exact results shows that the worst agreement occurs for the energy density T 00 at late times t = 4. Therefore we study in Fig. 11 and 12 the energy density on several z = const, y = 0 slices at times t = 2 and t = 4 respectively. At time t = 2 we still observe quite good agreement with the exact results, with partial O( 2 )-corrections again improving the approximation. At t = 4 the errors appear to be substantial. However, including partial O( 2 )-corrections again substantially improves agreement in the central region z = 0 and for large z. However, at intermediate values for the longitudinal coordinate z ≈ 2 a notable mismatch between approximation and exact results remains.  [3], we find the first order correction to the residual ∆ to be less than 2%. This is a consequence of the fact that only fields in S 0 given in (4.43) contribute to the residual at first order. Therefore, it appears to be a fluke that the (partial) second order corrected solution and the exact solution for vanishing z fit so well in the top left plot.

Conclusion
In this work we introduce the framework of a transverse derivative expansion to allow more efficient calculations of holographic models of heavy ion collisions. The first order correction to collisions of planar shocks allows us to reproduce the exact solutions in [3] remarkably well. The transverse derivative expansion has substantial computational advantages over the very memory and time consuming exact treatment. At a given point in the transverse plane, Green's functions of the differential operators that have to be solved on each time slice are identical to those computed for planar collisions. Transverse derivatives only contribute to source terms. Moreover, it is possible to decouple the Einstein equations at first order in transverse derivatives into two sets.
One of them only contributes negligibly to boundary observables. In section 4.4 and Appendix section B discuss why this happens. As a consequence the equations to be solved on each time slice, and each transverse pixel, up to first order in transverse derivatives reduce to those for planar shock collisions plus only four additional differential equations. This reduces the computation time, compared with the exact treatment of localized shock collisions, roughly by an order of magnitude.
In future works we intend to employ the methods described in this paper to study more realistic models of collisions of heavy ions including a granular structure and transverse energy density fluctuations, which are required to explain the observed large event-by-event fluctuations of flow observables. As found in [3] there is a considerable transverse flow, which develops before the hydrodynamic description becomes valid. The question of how sizable transverse energy density fluctuations of the localized shock's initial conditions affect the transverse dynamics before hydrodynamization remains an interesting and open question.
With these definitions we can write the so far unspecified terms appearing in (A.1a)-(A.1e) whereR ij is the radial diffeomorphism covariant, spatial Ricci tensor.

B Decoupling of transverse derivative expanded Einstein equations
In the following we are going to show that the differential equations of the functions belonging to set S 1 in (4.42) can be decoupled from the equations of the functions belonging to set S 0 in (4.43) through first order in transverse derivatives. Furthermore the differential equations of functions in S 0 do not depend on transverse derivatives of the zeroth order metric. Let the superscript (1) denote the first order transverse derivatives correction to zeroth order in transverse derivatives solutions which we denote by the superscript (0) . Then the following statements hold: 1. For solutions at zeroth order in transverse derivatives we have (ĝ ij ) (0) i =j = 0 as well as (F x ) (0) = 0, (F y ) (0) = 0.

Similarly ((G )
i =j do not depend on (ĝ ij ) (1) i=j and its derivatives nor on (Σ) (1) and its derivatives. Also ((G ) i j ) i=j do not depend on (ĝ ij ) (1) i =j or its derivatives. 5. Again by only using that (ĝ ij ) i =j = 0, it is easy to show that (∇ k (G ) k z ) (1) does not depend on transverse derivatives of zeroth order functions nor on (ĝ ij ) (1) i =j and its derivatives. And (∇ k (G ) k x,y ) (1) does not depend on (ĝ ij ) (1) i=j , (Σ) (1) and their derivatives. 6. After a slightly more tedious calculation it follows again from statement 1 that (R ij ) i =j does not depend on (ĝ ij ) (1) i=j , (Σ) (1) , (F z ) (1) and their derivatives and (R ij ) (1) i=j does not depend on transverse derivatives of zeroth order functions nor on (ĝ ij ) (1) i =j nor on F x , F y . The same is obviously true for∇ j F i and F i F j . 7. Statement 6 implies thatR = G ijR ij ,∇F and F 2 do not depend on transverse derivatives of zeroth order functions, nor on (ĝ ij ) (1) i =j nor on F x , F y .
14. From Eq. (4.18) it is apparent that the first order correction to the shift function λ (1) does not depend on transverse derivatives of the zeroth order metric, nor on (F x ) (1) , (F y ) (1) .