Asymmetric shockwave collisions in $\text{AdS}_{5}$

Collisions of asymmetric planar shocks in maximally supersymmetric Yang-Mills theory are studied via their dual gravitational formulation in asymptotically anti-de Sitter spacetime. The post-collision hydrodynamic flow is found to be very well described by appropriate means of the results of symmetric shock collisions. This study extends, to asymmetric collisions, previous work of Chesler, Kilbertus, and van der Schee examining the special case of symmetric collisions. Given the universal description of hydrodynamic flow produced by asymmetric planar collisions one can model, quantitatively, non-planar, non-central collisions of highly Lorentz contracted projectiles without the need for computing, holographically, collisions of finite size projectiles with very large aspect ratios. This paper also contains a pedagogical description of the computational methods and software used to compute shockwave collisions using pseudo-spectral methods, supplementing the earlier overview of Chesler and Yaffe.


Introduction and summary
Despite the fact that QCD is not conformal, supersymmetric, or infinitely strongly coupled, and has only a small number (N = 3) of colors, the comparison of heavy ion phenomenology with predictions based on AdS/CFT duality (of "holography") has turned out to be quite fruitful [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. At temperatures above the QCD phase transition the lack of supersymmetry is of minor importance and effects caused by the other differences can be described perturbatively, either on the QCD or gravity side of the duality. For example, corrections due to large but finite values of the 't Hooft coupling λ = g 2 YM N relevant for QCD can be calculated perturbatively on the gravity side, while the effects of non-conformality can be studied within QCD either perturbatively or using lattice gauge theory. Hence, it has been possible to identify which results from holographic modeling of heavy ion collisions should be more, or less, applicable to real QCD. Examples of observables with relatively modest corrections due to finite coupling and non-conformality effects include the viscosity to entropy density ratio [17], 4πη/s = 1+15 ζ(3) λ −3/2 ≈ 1.4 for λ ≈ 12, and the short hydrodynamization time predicted by AdS/CFT duality based on calculations of the lowest quasinormal mode (QNM) frequency [18]. For the latter quantity, finite coupling corrections are larger than for η/s, but not so much as to change the picture qualitatively.
(1. 6) In this frame the two incoming shocks will have widths w + and w − , and physical results may now depend on two independent dimensionless combinations which we take to be µw + and µw − . Over a substantial range of incoming shock widths {w + , w − } ranging from 0.35/µ down to 0.075/µ, we find that the spacetime region in which hydrodynamics is applicable has little or no dependence on the shock widths, or their asymmetry, and is sensitive only to the initial energy scale µ. Using the same definition of a hydrodynamic residual and the 15% figure of merit chosen in Ref. [1], we find that the boundary of the hydrodynamic region of validity remains at µ t hydro ≈ 2 , (1.7) even for highly asymmetric collisions. Similarly, the fluid 4-velocity resulting from asymmetric collisions remains very close to ideal boost invariant flow (1.1), while the post-collision proper energy density remains well-described by a Gaussian. However, the amplitude A, meanξ, and width σ of the Gaussian rapidity dependence are now functions of both incoming shock widths, (ξ, τ init ) = µ 4 A(µw + , µw − ) e − 1 2 (ξ−ξ(µw + ,µw − )) 2 /σ(µw + ,µw − ) 2 . (1.8) For asymmetric collisions, the outgoing energy density peaks at a non-zero mean rapidityξ which is well-described bȳ where the coefficient Ξ is constant for τ > 2 (as shown below in Fig. 6) and has the value Ξ ≈ 7 × 10 −2 . We find that the amplitude A is well-described by the geometric mean of the symmetric collision results, (1. 10) In fact, after shifting the rapidity byξ, we find that the geometric mean of the full symmetric collision rapidity distributions provides a good approximation to the asymmetric collision results. For the width of the rapidity distribution, this implies that For asymmetric collisions, the fit to the data provided by the this Gaussian model is good, as may be seen below in Fig. 7, but is not quite as perfect as for symmetric collisions. A more elaborate model, discussed in section 4.2, involves a weighted geometric mean of the symmetric collision profiles and provides an even better description, valid over a wider range of rapidity.
Given the above extension of the "universal" flow resulting from planar shock collisions to the asymmetric case, we now have the ingredients needed to predict initial conditions for the hydrodynamic flow resulting from collisions of bounded projectiles with finite transverse extent, provided the transverse size of the incident projectiles is large compared to their (Lorentz contracted) longitudinal widths, so that spatial gradients in transverse directions are small compared to longitudinal gradients. The following algorithm provides the leading term in an expansion in transverse gradients: • Regard the colliding system as composed of independent subregions in the transverse plane, or "pixels", with each pixel having a size δ ≡ 1/Q s which is small compared to the transverse extent of the projectiles, but large compared to their longitudinal widths.
• Let j label independent transverse-plane pixels, with p ± z (j) the portion of the longitudinal momentum of each incident projectile residing within pixel j.
• For each pixel j, transform to the CM frame in which the total longitudinal momentum within the pixel vanishes, and evaluate the resulting energy scale µ(j) and incident projectile widths w ± (j) for this pixel. Explicitly, µ(j) 6 = 4 p + z (j) p − z (j)/δ 4 .
• Transform each pixel's stress-energy tensor T µν (j) from its CM frame back to the original (lab) frame.
The result is a representation of the full system's stress-energy tensor on the τ init initial surface, with transverse variation on the pixel scale δ, suitable for use as initial data for further hydrodynamic evolution. This procedure uses strongly coupled holographic dynamics to map energy density profiles of the initial projectiles, which may include initial state fluctuations and have non-vanishing impact parameter, into hydrodynamic initial data, without the need to perform full 5D numerical relativity calculations which are very challenging [8]. As noted above this procedure, based on planar shock results, should be viewed as the first term in an expansion in (small) transverse gradients. It would, of course, be interesting to derive, systematically, subsequent terms in this expansion. Pixels near the periphery of the overlap region of colliding nuclei, illustrated in Fig. 1, will have decreasing CM frame transverse energy density µ 3 due to the rapid fall-off of the transverse energy density of the colliding nuclei near their periphery. Given the fact that the hydrodynamization time scales inversely with µ (1.7), this implies that pixels near the periphery of the overlap region (shown in orange) will enter the hydrodynamic regime much later than pixels in the middle of the overlap region. 1 How this impacts an appropriate choice of the initial Cauchy surface used in hydrodynamic modeling, and the resulting uncertainties in estimates of, for example, the elliptic flow parameter v 2 , is deserving of further study.
The remainder of this paper is organized as follows. In section 2 we review the characteristic formulation of general relativity in asymptotically anti-de Sitter spacetimes and the initial data for planar shock collisions, largely following Ref. [2]. Section 3 describes the numerical procedure and software used to compute shock collisions, highlighting several issues in greater detail than in Ref. [2]. Results are presented in section 4, followed by a brief final discussion in section 5. Readers primarily interested in results should feel free to turn directly to section 4. Additional computational details are presented in the appendix.
2 Planar shock collisions in asymptotically AdS spacetime

Characteristic formulation
As shown in Refs. [2-4, 8, 16], the characteristic formulation of general relativity, originally developed by Bondi and Sachs [19,20], provides a computationally effective method for handling the diffeomorphism invariance of general relativity when studying collisions dynamics in asymptotically AdS spacetimes. The characteristic formulation is based on a null slicing of the geometry in which coordinates are directly tied to a congruence of null geodesics. We will use X ≡ (x, r) to denote 5D coordinates, with x = (x 0 , x i ) ≡ (t, x i ) representing ordinary Minkowski coordinates on the boundary of the AdS spacetime. Requiring that t = const. surfaces be null hypersurfaces implies that the one-form k = ∇t is null, 0 = k A k A = g AB k A k B , which means that g tt = 0. Requiring the spatial coordinates x i to be constant along the null rays tangent to k A implies that 0 = k A ∂ A x i = g AB (∂ A t)(∂ B x i ), which means that g ti = 0. These conditions on the contravariant components of the metric then imply that g rr = g ri = 0. Hence, under these assumptions the most general line element may be written in the generalized infalling (or Eddington-Finkelstein) form, (2.1) It will be convenient to factor the spatial metric G ij into a scale factor Σ times a unimodular matrix g, with det( g) ≡ 1. One may fix one further condition, controlling the parameterization of the null geodesics tangent to k A . Bondi and Sachs [19,20] chose to fix the scale factor Σ(X) = r, convenient for problems with spherical symmetry. We instead follow Chesler and Yaffe [2] and choose to set This condition leaves a residual reparametrization invariance in the metric (2.1) consisting of radial shifts, r →r = r + δλ(x) , (2.4) with the shift δλ depending in an arbitrary fashion on the boundary coordinates x. Under such a shift, the metric coefficient functions transform as From these transformations of A and F i it is apparent that they may be regarded as temporal and spatial components of a gauge field representing radial shifts. It is possible to write the Einstein equations in a manner which is manifestly covariant under radial shifts. To do so, it is convenient to define modified temporal and spatial derivatives, Given these definitions, the Einstein equations, acquire a nested structure with the schematic form, Each equation is a first or second order linear radial differential equation for the indicated metric component(s) or their modified time derivatives. The square brackets of each coefficient or source function indicates on which fields the term depends. Explicit form of these equations, for the case of planar shocks, are given in appendix A. Given the rescaled spatial metric g on any time slice, plus suitable boundary conditions, each radial differential equation may be integrated in turn, thereby determining both the other metric coefficients and the time derivative of g on that time slice. The required boundary conditions may be inferred from the near-boundary behavior which can be obtained by solving equations (2.8a-2.8g) order by order in r. One finds [2], The coefficients a (4) , f (4) i and g (4) ij cannot be determined by a local near-boundary analysis. Note that g (4) ij is necessarily traceless (because g has unit determinant). These coefficients are mapped, via gauge/gravity duality, to the stress-energy tensor of the dual field theory. In our infalling coordinates this relation is given by [2] 0i ≡ −f (4) i , and h (4) ij ≡ g (4) ij . Here N c is the number of colors in the dual field theory, and η = diag(−1, +1, +1, +1) is the Minkowski metric tensor. Explicitly, The radial shift parameter λ(x) is completely undetermined in expansion (2.9) and may be chosen arbitrarily. As in previous work [2-4, 8, 16], we use this freedom to set the radial position r h (x) of the apparent horizon equal to a fixed value, It is sufficient to solve for the spacetime geometry in the region between the horizon and the boundary because information hidden behind the horizon cannot propagate outward and influence boundary observables. Thus, the choice (2.12) results in a convenient rectangular computational domain. With our metric ansatz (2.1), demanding a fixed radial position of the apparent horizon leads to a condition on d + Σ [2]. To derive this condition, one may write the tangents to a radial infalling null congruence in the form k A (X) = µ(X) ∇ A φ(X) for some scalar functions φ and µ. Demanding that the one-form k be null allows one to reexpress the time derivative of φ in terms of spatial derivatives. Requiring that the congruence satisfy the (affinely parameterized) geodesic equation k A k B;A = 0 allows one to reexpress the time derivative of the multiplier function µ in terms of its spatial derivatives. Given these time derivatives, one may then compute the expansion θ = ∇ · k on the time slice of interest. Demanding that the expansion vanish on a surface φ(X) = const. implies that this surface is an apparent horizon. Applying this procedure to the metric ansatz (2.1) and specializing to the case φ(X) = r leads to the desired condition [2], This condition must hold at all times if the radial position of the horizon is to remain fixed at some given value r h . Consequently, on every time slice the condition is also required to hold. When combined with the Einstein equation (2.8g), this final condition leads to an elliptic differential equation for the value of the metric function A on the (apparent) horizon. Explicit forms of the horizon equation (2.13) and the horizon stationarity condition (2.14) may be found in appendix A.

Solution strategy
To solve the nested form ( are needed as boundary conditions for the Σ, F i , and d + Σ equations and serve to fix the coefficient of a homogeneous solution to the corresponding differential equation. The asymptotic coefficients a (4) and f (4) i , proportional to the boundary energy and momentum density, are dynamical degrees of freedom (in addition to the metric g ij ) and are determined by integrating the stress-energy continuity equation as discussed below. The radial shift λ(x) will also be treated as a dynamical degree of freedom, as described below, and adjusted in a manner which ensures that the apparent horizon remains at a fixed radial position.
Given this boundary data, together with the value of g on some given time slice, the radial differential equations (2.8a)-(2.8d) may each be integrated in turn, at every spatial location x i , leading to a determination of d + g ij on the time slice. Two boundary conditions are needed to integrate the second order equation (2.8e) for the metric function A. As seen in table 1, the value of the radial shift λ supplies one condition. The second boundary condition is supplied by the value of A at the apparent horizon, which is determined by solving the horizon stationarity condition (2.14). field homogeneous solution(s) near-boundary behavior Table 1. Near-boundary asymptotic behavior of the homogeneous solutions to the radial differential equations (2.8a)-(2.8e) for the indicated fields, together with the desired asymptotic behavior of physical solutions. The asymptotic coefficients a (4) , f i , and g (4) ij determine respectively the energy density, momentum density, and traceless stress tensor of the dual field theory. The leading terms in the near-boundary behavior of all fields except Σ are driven by the inhomogeneous source terms in the various equations and do not correspond to homogeneous solutions.
Having determined both d + g and A, the actual time derivative for the rescaled spatial metric g is then reconstructed as (2.15) Knowing d + Σ and A (on a given time slice), the near boundary expansion (2.9) shows that the time derivative of the the radial shift λ(x) may be extracted as Similarly, the asymptotic coefficient g (4) ij determining the traceless stress tensor is extracted from the boundary limit of either r 4 ( g ij − δ ij ) or − 1 2 r 3 d + g ij . This information then allows one to determine the time derivatives of a (4) and f (4) i using the boundary stress-energy continuity equation, ∇ µ T µν = 0, which is an automatic consequence of the Einstein equations. Explicitly, ij . (2.17) The above procedure, involving integration of a sequence of linear ordinary differential equations in the radial direction plus one spatial elliptic equation on the apparent horizon, determines the time derivatives of the dynamical data { g ij , λ, a (4) , f (4) i } given initial values of this data on some time slice. These time derivatives are then input into a conventional time integrator, such as fourth order Runge-Kutta, to advance to the next time slice where the entire process repeats.
Overall, this characteristic formulation transforms the highly non-linear coupled Einstein equations into a set of nested linear ordinary differential equations and first order time evolution equations. We solve the radial differential equations, and the horizon stationarity equation, using spectral methods as described in some detail in section 3 and appendix C.

Planar shocks
By "planar shock" we mean an asymptotically anti-de Sitter solution of the vacuum Einstein equations whose boundary stress-energy tensor describes a "sheet" of energy density which moves at the speed of light in some longitudinal direction and is translationally invariant in the other two transverse spatial dimensions. For regular solutions, such a sheet of moving energy density will have some smooth longitudinal profile and non-zero characteristic thickness.
Let {x i } ≡ (x ⊥ , z) denote spatial coordinates separated into transverse and longitudinal components, and consider shocks moving in the ±z direction. To specialize the general infalling metric ansatz (2.1) to the case of planar shock spacetimes, we impose translation invariance in transverse directions plus rotation invariance in the transverse plane, which implies that all metric components are functions of only of r and x ∓ ≡ t ∓ z, that F i only has a longitudinal component, and that the (rescaled) spatial metric has the form [2], g = diag(e B , e B , e −2B ) . (2.18) Consequently, The boundary asymptotics (2.9) implies that the "anisotropy" function B behaves as For later computational convenience, let denote an inverted radial coordinate, so that the spacetime boundary lies at u = 0. In general it does not seem possible to find analytic forms of planar shock solutions using the infalling Eddington-Finkelstein (EF) coordinates (2.19). But analytic solutions are available in Fefferman-Graham (FG) coordinates [3,8,23]. Using {x µ ,ρ} ≡ {t,x ⊥ ,z,ρ} as our FG coordinates, withx ± ≡t ±z andρ an inverted bulk radial coordinate, the metric is a planar shock solution describing a shock moving in the ±z direction with arbitrary longitudinal energy density profile h(z). In the calculations described below, we use simple Gaussian profiles with width w and longitudinally integrated energy density µ 3 , The associated boundary stress-energy tensor is just with all other components vanishing.
Focusing, for ease of presentation, on shocks moving in the +z direction, the translational symmetries imply that the EF and FG coordinates will be related by a transformation of the form [2], As discussed above, the required initial data for the characteristic evolution scheme consists of the anisotropy function B plus the boundary data {a (4) , f (4) z } and the radial shift λ. Inserting a transformation of the form (2.25) into the FG metric (2.22), a short exercise [2] shows that while the boundary data is given by To solve for the transformation functions {α, β, γ}, one approach, used in Refs. [2,3], is to insert the transformation (2.25) into the FG metric (2.22) and demand that the result have the EF form (2.19). 2 To simplify the resulting equations, it is helpful 2 An alternative approach, used in Ref. [8] for more general metrics, is based on observing that the curve defined by fixed values of the EF boundary coordinates and all values of r, X A (r) = (t 0 , x i 0 , r), is a null geodesic of the EF metric (2.1) with r an affine parameter. Therefore the same path in FG coordinates, Y (X(r)), will satisfy the geodesic equation denoting the FG coordinate Christoffel symbols. Explicit forms of the resulting equations can be found in appendix B.
to redefine the transformation functions α and β via (2.28) One finds [2] that the functions ζ and δ satisfy a pair of coupled differential equations, . The desired solutions have the nearboundary behavior Integrating equations (2.29) with boundary conditions ensuring the behavior (2.29c), and inserting the resulting transformation functions into Eqs. (2.26) and (2.27), yields the anisotropy function B and associated boundary data describing of a single shock. To construct initial data for colliding shocks, we superpose counter-propagating single shock data at an initial time t 0 when the two shocks are sufficiently widely separated that their overlap is negligible, However, unlike for the other functions, the overlap of the radial shifts λ ± of the left and right moving shocks in the region close to z = 0 is significant. Since we choose the shocks on the first time slice to be well separated, we may regard the geometry in between as deviating negligibly from pure AdS. This justifies modifying the initial shift function λ in the neighborhood of z = 0, without changing the physical data , we adjust the initial radial shift by setting In practice, we slightly modify the above superposition procedure. Following Refs. [2,3], we replace Eq. (2.30b) with From the form (2.11) of the stress-energy tensor, one sees that 0 is a constant additive shift in T 00 . In other words, 0 is an (artificial) uniform background energy density. Adding a small background energy density helps alleviate numerical problems, as discussed below, and physically means that the colliding shocks will be propagating through a background thermal medium. If the background energy density 0 is sufficiently small compared to the energy densities in the colliding shocks, then the background will effectively be very cold (compared to the energy scale µ of the shocks) and there will be little dissipation to the medium. This modification is done purely for numerical convenience and we will be interested in results extrapolated to vanishing background energy density.

Computational methods and software construction
The aim of this section is to describe the construction of a planar shockwave collision code in sufficient detail so that an interested reader could create their own version with relatively modest effort. Those primarily interested in results should skip to the next section.

Transformation to infalling coordinates
As explained in Ref. [2] and the previous section, the transformation from Fefferman-Graham to infalling coordinates may be computed by first solving for the congruence of infalling geodesics in FG coordinates. Or, in the special case of planar shock geometries, one can directly solve the simplified transformation equations (2.29). We implemented both approaches, and found them to have comparable numerical efficiency. Here, we focus on the direct approach of solving Eqs. (2.29) for the case of a right moving shock. Henceforth, for convenience, we also set µ = 1. Appropriate factors of µ can always be reinserted via dimensional analysis. We solve the coordinate transformation equations (2.29) in the rectangular region u ∈ [0, u end ], z ∈ [−L z /2, L z /2] using Newton-Raphson iteration (i.e., linearizing each equation in the deviation of the solution from the current approximation), and solving the resulting linear equations using spectral methods with domain decomposition. 3 Periodic boundary conditions are imposed in the longitudinal direction and functions of z are approximated as truncated Fourier series. This is exactly equivalent to characterizing any function f (z) by a list of its values, {f l ≡ f (z l )}, on an evenly spaced Fourier grid composed of N z points, for k = 0, · · ·, N z −1. Derivatives with respect to z turn into the application of a Fourier grid differentiation matrix D z = (D z ) kl applied to the vector of function values, Explicit expressions for the Fourier grid differentiation matrix components (D z ) kl are given in appendix C. A rather fine longitudinal grid is required to accurately represent thin shocks within a large longitudinal box. We used Fourier grids with N z = 960 for L z = 12 and shock widths down to 0.075.
To represent the dependence of functions on the radial coordinate u we first decompose the domain [0, u end ] into M equally sized subdomains, and then use a Chebyshev-Gauss-Lobatto grid with N u points within each subdomain. This amounts to using a radial grid composed of the points for j = 1, · · ·, M and k = 0, · · ·, N u −1. The radial dependence of some function g(u) is represented by the list of M × N u function values on this grid, {g jk ≡ g(u jk )}, and derivatives with respect to u turn into the application of a (block diagonal) Chebyshev differentiation matrix D u applied to this list of function values, Explicit expressions for the components of the Chebyshev differentiation matrix D u are given in Eq. (C.12). As discussed in Ref. [2], using domain decomposition (i.e., M > 1) helps to avoid excessive precision loss in the numerical evaluation of equations near the u = 0 boundary, and allows the use of a larger time step without running afoul of CFL instabilities. To integrate radial equations down to u end = 2, we used radial grids with up to M = 22 domains and N u = 12 points within each subdomain. The product of these 1D grids defines our 2D spectral grid. Any function f (u, z) becomes a set of N tot ≡ M × N u × N z values on these grid points, For each set of equations, linearization around some initial, or current, guess for a solution leads to a set of linear equations of the generic form M f = −S, where f is the unknown vector of function deviations from the current guess, S is the vector of residuals, and M is the spectral approximation to the linear operator which results from the linearization of the differential equation(s) at some given value of z.
At this point, these linear equations are singular. First, u = 0 is a regular singular point of the differential equations (2.29a) and (2.29b); one cannot simply evaluate, numerically, these equations at u = 0. Moreover, solutions to these differential equations are, of course, non-unique. One must complement the differential equations with suitable boundary conditions to specify a unique solution. With spectral methods, fixing one of these problems fixes the other. Prior to linearization, one simply replaces the (ill-defined) evaluation of the equations at u = 0 by constraints encoding required boundary conditions.
Examining equations (2.29a) and (2.29b), one sees that the most general nearboundary behavior is for arbitrary values of the coefficients ζ −1 , λ, γ 0 and δ 0 . We want to set the leading coefficients ζ −1 , γ 0 and δ 0 to zero. To implement this Dirichlet condition for γ and δ in a manner which avoids unnecessary precision loss when computing derivatives of these functions at the boundary, it is convenient first to redefine and then reexpress equations (2.29) in terms ofγ andδ. Unwanted solutions with non-zero boundary values for γ or δ are then simply not representable when using our spectral representation forγ orδ. Similarly, using our spectral representation for ζ automatically eliminates unwanted solutions where ζ has singular 1/u behavior. The continuum differential equations imply thatγ andδ both vanish, and have vanishing first derivatives, at the boundary. To deal with the u = 0 regular singular point in the discretized equations forγ andδ it is sufficient to replace the equations at u = 0 with constraints settingγ andδ to zero. If we wished to fix the radial shift λ by simply specifying its value, we could similarly redefine ζ = λ + uζ and requireζ to vanish at the boundary. However, we found it more convenient to fix λ indirectly by demanding that ζ vanish at our chosen value of u end . Referring to Eqs. (2.25) and (2.28), one sees that this condition will make the u = u end surface coincide with a surface of constant FG radial coordinate,ρ = u end . In other words, with this condition the FG computational domainρ ∈ [0,ρ end ] is the same as the EF domain u ∈ [0, u end ].
The net effect of the above procedure, in the discretized equations for ζ,δ andγ at longitudinal position z l , is to replace the the (degenerate) equations at u = 0 by the respective constraints 4 ζ M,Nu−1,l = 0 ,δ 1,0,l = 0 ,γ 1,0,l = 0 . (3.8) In addition to applying boundary conditions at u = 0, when using domain decomposition one must also impose continuity conditions at subdomain boundaries. Our set (3.3) of radial grid points redundantly duplicates the interior endpoints of each subdomain, u j,Nu−1 = u j+1,0 for j = 1, · · ·, M −1, and hence two different rows of the linear equation Mf = −S represent the differential equation evaluated at the same physical point. One could deal with this by eliminating the duplication of subdomain endpoints and suitably redefining the differentiation matrix D u . But it is even easier to fix the problem by simply replacing one of the rows representing an interior subdomain endpoint with a constraint equation enforcing the equality of duplicated function values at this point, f j,Nu−1,l − f j+1,0,l = 0. 5 After these row replacements, the modified linear system is reasonably well conditioned and, with a sufficiently good initial guess, Newton iteration rapidly converges quadratically. To generate an initial guess, it is natural to work sequentially in z. If the shock is propagating in the +z direction with the profile function h(z) having its maximum at z = 0, then at the furthest point behind the shock, z 0 = −L z /2, the geometry deviates negligibly from pure AdS and ζ =γ =δ = 0 is a fine initial guess. Thereafter, we use the converged solution at each z i as an initial guess for the solution at z i+1 . This provides a good initial guess provided the longitudinal grid spacing is sufficiently fine.
The above procedure for solving the transformation equations (2.29) using spectral methods works well as long as the radial depth u end to which one integrates is not too large. The key advantage of this approach is that the precision of the obtained solutions do not degrade near the boundary, even through u = 0 is a singular point of the differential equations. That is to say, spectral methods are excellent for finding well-behaved solutions of equations having regular singular points. However, as u end increases the linear operators one inverts in this Newton iteration scheme become increasingly ill-conditioned. Unfortunately, the depth to which one must integrate in order to locate the apparent horizon (discussed next) after superposing shocks grows with increasing separation of the initial shocks. We used two strategies to cope with this difficulty.
First, following Refs. [2,3], we added a small artificial background energy density 0 when superposing shocks as described above. Increasing the background energy density decreases the depth at which an apparent horizon forms. Second, after using the above approach to find the transformation functions for u < u end , we integrate further into the bulk by switching to an adaptive 4th order Runge-Kutta integrator, with the spectral solution at u end providing initial data. (A description of this standard integrator is given in appendix E.) For simplicity, we choose to integrate to a fixed value u = u max , instead of a fixed value ofρ.
For our chosen range of shock parameters, with widths down to w = 0.075, using a spectral grid down to u end = 2 worked well. With a longitudinal box size L z = 12 and background energy densities in the range of 1-5% of the peak energy density, it turned out that only a modest further integration with the adaptive integrator down to u max = 2.11 was sufficient to reach the apparent horizon throughout the longitudinal box. 6 Having transformed a right-moving single shock solution to infalling coordinates, and extracted the resulting initial data {B + , a z+ , λ + } for evolution using Eqs. (2.26)-(2.27), a simple reflection generates corresponding data for a left-moving shock, We construct initial data for counter-propagating shocks by combining single shock solutions as described earlier in Eqs. (2.30)-(2.32). We chose the initial time t 0 for this superposition so that the initial separation between the shocks, ∆z 0 = −2t 0 , is large compared to the shock widths. We used ∆z 0 = 4 for symmetric collisions of broad shocks, ∆z 0 = 2 for symmetric collisions of thin shocks, and ∆z 0 = 3 for asymmetric collisions of shocks.
For thin shockwave collisions with small background energy density, avoiding numerical instabilities associated with short wavelength perturbations is challenging. As discussed in Ref. [2], it is helpful to damp discretization induced perturbations using appropriate filtering. We constructed and applied smoothing filters to the initial data in both longitudinal and radial directions. Details of these filters are presented in appendix D.2.

Horizon finding
After transforming chosen single shock solutions to infalling coordinates, as just discussed, and combining two counter-propagating shocks as shown in Eqs. (2.30)-(2.32)), the final step in the construction of initial data is locating the apparent horizon which serves as an IR cutoff in the bulk. 7 In our planar shock geometries, the apparent horizon condition (2.13) becomes (3.10) A radial shift, r =r + δλ, corresponds in our inverted radial coordinates to u =ū 1 +ū δλ . (3.11) Ifū ∈ [0, u max ] represents the radial coordinate used in the transformation to infalling coordinates, then we wish to determine the value of a further shift δλ(z) such that condition (3.10) holds at some value of u h ≡ 1/r h which may, for convenience, be chosen to equal the same value u max from the coordinate transformation. With this choice, δλ must be negative for the sought-after apparent horizon to lie within the coordinate transformation domain. Equation (3.10) is a nonlinear but ordinary differential equation for the shift function δλ(z). To solve it, we use spectral methods (with the same Fourier grid in z) combined with a root finding routine. Linearizing equation (3.10) in δλ allows us to apply Newton iteration. Each iteration step starts with a trial value of the radial shift, δλ (m) in iteration m, and computes the residual (i.e., the right-hand side of Eq. (3.10)) and its variation with respect to δλ, and solves the linearized equation to find an improved value δλ (m+1) of the shift.
To evaluate the residual and its variation, we first integrate Eqs. (2.8a)-(2.8c), using the current value of B(z, u) and λ(z), to find the auxiliary functions Σ, F and d + Σ. 8 After each step we convert the spectral representation of B(z, u) to a new radial grid with grid points shifted according to Eq. (3.11). To do so, we perform off-grid interpolation using a sum of Chebyshev cardinal functions [21] with coefficients given by the on-grid values of B(z, u).
For our settings of longitudinal box size and shock parameters, we found it advantageous to choose the initial guess δλ (0) to be 0.1. It was also helpful to start with a relatively large background energy density 0 of about 10% of the peak shock energy density, and then gradually decrease 0 during each iteration step until it reached the desired final value before Newton iteration convergence.
During time evolution, described next, solving the horizon stationarity condition (2.14) on each time step yields the time derivative of the radial shift thereby providing the information needed to integrate λ forward in time. (The explicit form of Eq. (2.14) for our planar shock geometries is given in Eq. (A.4).) To prevent discretization errors from driving long term drift away from the desired horizon condition (3.10), we also directly recomputed the apparent horizon position every 10-100 time steps using the above iterative procedure.

Time evolution
As described above in section (2.2), the data on any time slice needed to integrate forward in time consists of {B(z, u), a (4) (z), f (4) (z), λ(z)}. To compute the time derivative of this data, we successively solve Eqs. (2.8a)-(2.8e) as discussed earlier. Explicit forms of these equations are given in Eqs. (A.2a)-(A.2e) of appendix A. We use the same multi-domain spectral methods described above in section 3.1. These methods presume that functions being represented by their values on the spectral grid are well behaved throughout the computational domain. 9 Our functions Σ and A have divergent near-boundary behavior, as shown in Table 1, so for computational purposes we use redefined functions in which the leading near-boundary behavior is subtracted. For most functions, we also choose redefinitions such that the new functions either have known non-zero boundary values or vanish linearly at the boundary. Specifically, we use the following redefinitions, We use factors of u/(1 + uλ) = (r + λ) −1 in these redefinitions, instead of pure powers of u, so that the new functions transform simply under radial shifts. This is natural as it preserves manifest radial shift covariance in the equations for the new functions, but is not essential. In relations (3.12d) and (3.12e), and henceforth, d + σ and d + b are simply names for redefined functions encoding d + Σ and d + B, respectively, and are not themselves modified d + time derivatives applied to σ or b.
Referring to Table 1 and Eq. (2.9), one sees that the new functions b and d + b vanish linearly as u → 0, while f and d + σ have non-zero boundary values of f (4) z and a (4) , respectively. The new function a has a boundary value of −∂ t λ which is an output, not an input, of the radial integration determining a.
Arranging to have constant or linear near-boundary behavior of redefined functions minimizes the precision loss which can occur when evaluating derivatives very near the boundary. In particular, extracting the third power of u/(1 + uλ) in the definition (3.12a) of b is essential for the numerical stability.
After inserting the redefinitions (3.12) into the relevant radial equations (A.2a)-(A.2e), it is crucial to simplify the resulting equations, prior to numerical implementation, in such a way that cancellations of terms with the most divergent near-boundary behavior are performed exactly, analytically. When each radial differential equation is written in canonical form (with a unit coefficient of the highest order u-derivative), no term in the inhomogeneous source term of the equation should be more singular than 1/u for first order and 1/u 2 for second order equations, otherwise unnecessary precision loss will occur during the numerical evaluation of the equation. 10 In solving the successive radial equations Eqs. (2.8a)-(2.8e) [or (A.2a)-(A.2e)], we implement the following boundary conditions at u = 0 using the row replacement technique described in section 3.1, Finally, for the A equation (2.8e) [or (A.2e)], after conversion to an equation for the new function a the non-normalizable homogeneous solution is automatically excluded by the spectral representation for a. One boundary condition is needed to fix the coefficient of the normalizable homogeneous solution. Referring to table 1, specifying the boundary value of a is the same as fixing the time derivative of the radial shift, 10 Such analytic simplification, eliminating what would otherwise be huge cancellations near the boundary, is essential when performing calculations using machine precision (64 bit) arithmetic. If one instead uses arbitrary precision arithmetic (in, for example, Mathematica), one might think such careful simplification prior to programming is unnecessary. However, failure to properly simplify expressions will then require the use of extraordinarily high precision arithmetic with concomitant poor performance. a(0, z) = −∂ t λ(z). But we do not wish to input some arbitrary choice for this time derivative. Instead, prior to solving the A equation (2.8e) we first solve the horizon stationarity condition (2.14) which determines the value of A on the horizon. This is an inhomogeneous differential equation involving A and its first and second order longitudinal derivatives, evaluated on the horizon. The explicit form is given in (A.4) of appendix A. Then, to solve the radial equation (2.8e) [or (A.2e)], converted to an equation for a, we replace the u = 0 row in the spectral discretization of this equation with a row fixing the value of a at the horizon, i.e., equating a(u max , z) to the value determined by the horizon stationarity condition.
To recap, every time step begins with the sequential solution of equations (A.2a)-(A.2e), plus the horizon stationarity condition, using the same spectral methods and Chebyshev grid employed in the preparation of initial data. This yields Σ, F , d + Σ, d + B and A, from which the ordinary time derivatives of B, a (4) , f (4) z and λ are extracted using relations (2.15)-(2.17). This is the information needed to integrate forward in time.
To perform time integration we use a discrete approximation with non-zero time step δt. We specifically choose the well known fourth order Runge-Kutta method (RK4), which uses four "substeps" each involving the evaluation of time derivatives, performed as described above for each point on our longitudinal grid. (The relevant RK4 formulas are shown in appendix E.) A time step δt = 0.002 was used in all integrations, which was sufficient to deliver stable evolution for all shock widths considered. For broader shocks (w > 0.3) a lower order time integrator would have sufficed, but for shocks with width w < 0.1 we found using at least RK4 to be essential, with our time step, to achieve accurate results. After each time step of the evolution we filter the final results for the propagating data {B, a (4) , f (4) , λ} in the longitudinal direction using a low-pass filter, as detailed in D, which damps the upper third of the spectral bandwidth. The filtering is applied to the final RK4 outputs, not during the RK4 substeps. This damps short wavelength discretization-dependent fluctuations; such filtering should be viewed as a part of the spectral discretization prescription. We do not apply filtering to any interim results while solving the radial equations (2.8a)-(2.8e).

Calculated collisions
Using the above described techniques and associated software, planar shock collisions were computed for various combinations of incoming shock widths. All initial shocks  Table 2. Physical and computational parameters of specific computed collisions. Shown are the incoming shock widths w ± , number of longitudinal grid points N z , and background energy densities 0 . Shock widths w ± are measured in units of µ −1 . The background energy density 0 is in units of the peak energy density of the narrower shock, or µ 3 w −1 + / √ 2π. Computed results at the two listed values of 0 were used to extrapolate to vanishing background energy density.
had Gaussian profiles (2.23) and identical transverse energy density µ 3 . In units in which µ ≡ 1, shock widths ranged between 0.075 and 0.35. For technical reasons involving the damping of numerical artifacts, as discussed above, an artificial background energy density was added whose size ranged from 5.5% down to 1.2% of the peak energy density of the narrower shock. Periodic boundary conditions were applied in the longitudinal direction, with this dimension then discretized with a uniformly spaced (Fourier) grid having of up to N z = 720 points. The longitudinal period L z was set to 10, 11, or 12 for collisions of narrow, asymmetric, or broad shocks, respectively. In the radial direction, domain decomposition with M = 22 subdomains of uniform size in the inverted radial coordinate u = 1/r was used, with a Chebyshev-Gauss-Lobatto grid of N u = 13 points within each subdomain. Time evolution used RK4 time-stepping with a step size δt = 0.002 and total time duration ranging from t = 5/µ to t = 6/µ. Table  2 lists the parameters of specific calculations. Figure 2 shows the energy density T 00 (t, z), in units of µ 4 , for three representative collisions. The top row shows symmetric collisions of shocks with widths w ± = 0.35 (upper left) and w ± = 0.075 (upper right), while the lower row displays results from the corresponding asymmetric collision with (w + , w − ) = (0.075, 0.35).
Local maxima in the energy density are present on the forward lightcone, as clearly seen in Fig. 2. These local maxima lie outside the hydrodynamic region (discussed below). In asymmetric collisions, the width of a given postcollision local maxima largely reflects the width of the corresponding incoming projectile. As shown in Fig. 3, the amplitude of these local maxima decay with the same power-law time dependence seen in symmetric collisions.

Hydrodynamic flow
At every spacetime event inside the forward lightcone of a collision, the timelike eigenvector and corresponding eigenvalue of the holographically computed stress-energy tensor determine the fluid 4-velocity u µ and proper energy density , 11 with normalization u µ u µ = −1 and u 0 > 0. Given the flow velocity and energy density, we use the first order hydrodynamic constitutive relation to construct the hydrodynamic approximation to the stress-energy tensor, where the viscous stress (to first order in gradients) is given by For the conformal fluid of N = 4 Yang-Mills theory, the pressure p = /3 and the shear viscosity η = ( /3) 3/4 / √ 2. 12 Following Ref. [1], we define the spacetime region R in which hydrodynamics provides a good description as the largest connected region within the future lightcone in which the normalized residual, measuring the difference between the holographically computed stress energy tensor and its hydrodynamic approximation, is smaller than 0.15. For all collisions studied, symmetric and asymmetric, with various combinations of incoming shock widths ranging from 0.35 down to 0.075, we found that the boundaries of the region R differ very little from one another, as illustrated in Fig. 4. 13 At z = 0 we find that time at which hydrodynamics first becomes valid (i.e., ∆ < 0.15) to be essentially the same for asymmetric and symmetric collisions and given by t hydro ≈ 2 . (4.5) In the symmetric case this confirms earlier results found in Refs. [1,2]. For symmetric collisions, we reproduced the key results of Ref. [1]: boost invariant flow (1.1) within the hydrodynamic region to within a precision of O(10 −3 ), Gaussian 13 By suitably adjusting the filtering of discretization induced artifacts, as discussed in the appendix D, we could decrease the background energy density in our computations of asymmetric collisions to about 1% of the peak value of the energy density of the narrower shock. For asymmetric collisions it turned out to be quite challenging to achieve high precision and numerical stability with significantly smaller background energy densities. In this and subsequent figures, we perform a linear extrapolation to vanishing background energy density using calculated results at the non-zero background energy densities shown in table 2. At sufficiently late times, this linear extrapolation ceases to be a reliable approximation to the limit of vanishing background energy density. A simple linear extrapolation, with our values of 0 , is adequate in the t ≤ 4 interval displayed in Fig. 4, which coincides with the time interval shown in Ref. [1] of the hydrodynamic region R. Turning to asymmetric collisions of shocks with differing widths, we again find that flow within the hydrodynamic region R is very close to ideal boost invariant flow (1.1), as illustrated in Fig. 5 for rapidity ξ ∈ [−1, 1]. Moreover, the rapidity distribution of the proper energy density on a surface of constant proper time τ τ hydro continues to be well approximated by a Gaussian but now with a peak which is shifted away from vanishing rapidity: Our results for the rapidity shiftξ(w + , w − ; τ ) are shown in Fig. (6) for three examples of asymmetric collisions. To a good approximation, the width dependence of the rapidity shift has a simple factorized form for τ 2, with a coefficient Ξ ≈ 0.07 that is essentially constant for τ > 2.
for the same three cases.
We find that the rapidity distribution of of the proper energy density for the asymmetric collisions is well approximated by the shifted geometric mean of the corresponding symmetric collision results, (4.8) The efficacy of this relation is illustrated in Fig. 7, which shows the proper energy density as a function of the rapidity ξ at proper times τ = 2 (top) and 3 (bottom) for the case of (w + , w − ) = (0.075, 0.35) (left) and (w + , w − ) = (0.075, 0.25) (right). In each plot the solid blue line shows the asymmetric collision result while the red dashed curve shows the shifted geometric mean of the corresponding symmetric collision results. For |ξ| < 1 this model fits almost perfectly, while for |ξ| > 1 small deviations from this simple description begin to show.
To motivate a more elaborate model which captures these deviations from the simple model (4.8), let denote the generalized mean with power p of some quantity X which is observable in symmetric collisions of shocks with widths w + and w − , and then define p[X] as the power for which the generalized mean of symmetric collision results gives the result X(w + , w − ) of this observable in an asymmetric collision with shock widths (w + , w − ). In other words, p[X] is the solution to the equation (4.10) Recall that the geometric mean is the p → 0 limit of the generalized mean (4.9).  directly compares the amplitude A(τ ) for asymmetric collisions with the geometric mean of the corresponding symmetric collision results. For times τ > 2, the difference is negligible.
To construct an improved model, let denote the Gaussian of a symmetric collision rapidity distribution (without the corresponding amplitude). Then replace the geometric mean of the simple model (4.8) by a biased mean of symmetric collision Gaussians, where, once again, A(w ± , τ ) is the amplitude of the rapidity distribution for symmetric collisions of width w ± , and the rapidity shiftξ is given in Eq. (4.7). If the bias a(w + , w − ; ξ, τ ) vanishes, then this form reduces to the previous simple model (4.8). Fitting the improved model (4.12) to our numerical results, we find that the resulting bias function a(w + , w − ; ξ, τ ) is remarkably insensitive to the widths (w + , w − ) and is also constant for τ > 2 to quite good accuracy. Our results for a are displayed in Fig. 10 for the the cases (w + , w − ) = (0.075, 0.35), (w + , w − ) = (0.075, 0.25), and (w + , w − ) = (0.1, 0.25) which, as shown, differ negligibly from each other. The resulting bias function a(ξ) is well described by the simple universal function a(w + , w − ; ξ, τ ) ≈ a(ξ) ≡ − 1 4 tanh ξ .  To show the efficacy of the improved model (4.12) and the improvement as compared with the simple model (4.8), we again compare in Fig. 11 the proper energy density rapidity distributions from asymmetric collisions along with the predictions of the above improved model (4.12) with bias function (4.13), for the same cases shown earlier in Fig. 7. As one sees, the curves are now essentially indistinguishable.  Figure 11. The proper energy density as a function of rapidity ξ at constant proper time τ = 2 (first row) and τ = 3 (second row) for asymmetric collisions with (w + , w − ) = (0.075, 0.35) (left) and (w + , w − ) = (0.075, 0.25) (right), displayed by the solid blue curve. The overlaid red dashed curve shows the result obtained from the improved model (4.12), with bias function a(ξ) = − 1 4 tanh ξ, and the respective Gaussian distributions for the corresponding symmetric collisions.

Discussion
The goals of this work were twofold: On the one hand by studying and quantitatively modeling asymmetric planar shock collisions via holography, we aim to help bridge the gap between descriptions of very early states of a quark gluon plasma formed during heavy ion collisions and the later hydrodynamic regime to which the system evolves. On the other hand, we also hope that a relatively didactic and detailed description of the computational techniques and software construction will be useful to others.
By studying asymmetric collisions of planar shockwaves in AdS 5 , we found that the simple "universal flow" description of symmetric shock collisions, found in Ref. [1], generalizes very naturally to asymmetric shock collisions. Within the hydrodynamic regime, the fluid flow is extremely close to ideal boost invariant flow, while the proper energy density has a Gaussian rapidity dependence. Characterizing the dependence on the amplitude and widths of the initial shocks enabled the construction of a simple model for mapping initial state energy density distributions to hydrodynamic initial data, valid to leading order in transverse gradients and having potential applicability to non-central collisions of highly relativistic nuclei.
The hydrodynamization time was confirmed to be very insensitive to the widths of the colliding shocks, and dependent only on the CM frame energy density. Viewing asymmetric collisions of planar shockwaves as models for "pixels" within non-central collisions of finite sized projectiles with large aspect ratios, this result implies that the hydrodynamization time, measured in the lab frame, increases towards the fringes of the almond-shaped overlap region that forms the post-collision quark-gluon plasma. Suitably modeling the initial state transverse energy density as a function of the distance to the center of the Lorentz-contracted nuclei allows one to estimate the hydrodynamization time of different layers of the quark-gluon plasma.
Possible topics for future work include the analysis of non-local observables and entropy production during asymmetric collisions of planar shocks, explicit comparison of holographic results for localized shock collisions with our model for hydrodynamic initial data, and systematic incorporation of higher terms in an expansion in transverse gradients into this model.

A Einstein equations for planar shocks
In this section we write down explicit forms for the Einstein equations (2.8a)-(2.8g) for planar shocks. We parametrize the rescaled spatial metricĝ aŝ with a single anisotropy function B(u, t, z). (Recall that u ≡ 1/r.) The time-space metric components F x and F y vanish due to rotational invariance in the transverse plane and, for brevity, we write just F below in place of F z for the remaining timespace component. The resulting Einstein equations in our infalling coordinates have the schematic form: which specialize the general infalling form (2.8) to the case of planar shocks. Denoting radial derivatives with primes, f ≡ ∂f /∂r, and f ,z ≡ ∂f /∂z for longitudinal derivatives, the explicit form of the various coefficient and source functions are as follows: The condition (2.13) that the apparent horizon lie at a fixed radial position r h has the explicit form (3.10). For planar shocks the horizon stationarity condition (2.14) becomes:

B Transformation to infalling coordinates
The metric withx ± ≡t ±z, describes a single shock moving in the +z direction using Fefferman-Graham (FG) coordinates. It gives a solution to the Einstein equations for any longitudinal profile function h(x + ), To construct initial data decsribing two counterpropagating shocks, we first transform a single shock solution to the infalling Eddington-Finkelstein (EF) form, with the metric functions A, F , Σ, and B depending only on t−z and the inverted radial coordinate u ≡ 1/r. In other words, the components g uA all vanish except for g ut = −u −2 . We relate the FG and EF coordinates according tõ Demanding that this change of coordinates yields a metric of the desired form, i.e., leads to the following equations for the transformation functions, arising from the specified values of (g EF ) uu , (g EF ) uz , and (g EF ) ut + (g EF ) uz , respectively.
Here primes denote radial derivatives ∂/∂u, and H ≡ h (t − z + u + α + γ). The dependence of the functions H, α, β and γ on their two arguments of t−z and u is suppressed for brevity. The desired solutions to Eqs. (B.5) have the near-boundary behavior Following Ref. [2], it is helpful to redefine α and β in terms of two new functions δ and ζ via Inserting these expressions into equations (B.5) and taking appropriate linear combinations of the results leads to a pair of coupled equations for δ and ζ, plus a single decoupled equation for γ, Alternatively, starting from the infalling form (B.2), it is easy to show that curves along which r ≡ 1/u varies with all other coordinates held fixed are null geodesics (with r as an affine parameter). Since coordinate transformations are isometries, the same curves must satisfy the geodesic equation expressed in FG coordinates, i.e.
where Γ are the Christoffel symbols evaluated in FG coordinates. The solution Y A (r) to this geodesic equation which begins at boundary coordinates x µ = (t, x ⊥ , z) with null tangent d dr Y A (r) = (δ A t + δ A ρ ) on the boundary directly gives the FG coordinates corresponding to the event with EF coordinates of x M = (t, x ⊥ , z, 1/r). Parametrizing the resulting coordinate transformation using Eq. (B.3), the non-trivial t, z and u components of the geodesic equation (B.10) lead to second order equations for the transformation functions, (B.11c)

C Pseudo-spectral methods
Pseudo-spectral methods are a class of mean weighted residual approximation techniques. These methods provide highly efficient techniques for constructing accurate numerical approximations to linear differential equations of the form with L being a linear differential operator. One approximates the solution f by a linear combination of a finite set of of basis functions, f (N ) = N −1 m=0 c m φ m , and defines the residual Given some scalar product (·, ·) for the function space in which the basis functions φ m reside, and a chosen sequence {ξ m } of test functions, one solves for the coefficients {c m } of the spectral approximation f (N ) by demanding that the residual vanish on these test functions, for m = 0, · · ·, N −1. Different weighted residual methods are distinguished by the choice of the test functions. So-called "pseudo-spectral" or "collocation" methods are a subclass of mean weighted residual algorithms in which one chooses the test functions to have point support. In, for example, one dimensional problems one chooses for some selected set of points {x m }. In other words, in pseudo-spectral approximations, one demands that the residual vanish identically on some discrete set of grid points distributed across the computational domain. For a given basis set {φ m }, m = 0, · · ·, N −1, there is a corresponding optimal choice of grid {x m }, m = 0, · · ·, N −1, namely the abcissas of a Gaussian quadrature integration scheme associated with this basis set [21]. The original spectral approximation f (N ) = N −1 m=0 c m φ m is then exactly equivalent to a linear combination of cardinal functions, in which each coefficient is the value of the function approximation on a given grid point, For one dimensional problems on a finite interval, the most commonly used basis functions are Chebyshev polynomials. There are actually two corresponding sets of optimal spectral grids differing in whether the endpoints of the interval are themselves gridpoints. It is easiest to deal with boundary conditions when endpoints are included in the spectral grid in which case, for the interval [−1, 1], the appropriate N +1 point grid consists of the points x m = cos(mπ/N ) , m = 0, · · ·, N .
This is sometimes referred to as a Chebyshev-Gauss-Lobatto grid. For problems on a periodic interval, a truncated Fourier series provides the most useful spectral approximation. In this case, an appropriate 2N point spectral grid consists of 2N evenly spaced points around the periodic interval. So, for the interval [0, 2π], one may use x m = π m/N , m = 0, · · ·, 2N −1 . (C.8) If the differential equation of interest (C.1) involves an M -th order differential operator, L = M k=0 p k (x) d k dx k , then computing the values of the residual R on all grid points, using the cardinal representation (C.6), requires the evaluation of up to M -th order derivatives of each cardinal function at every point on the grid. This computation need only be performed once, and defines a set of "spectral differentiation matrices" with components Given these matrices, the application of the differential operator L to the spectral approximation of some function reduces to the application of the finite matrix L (N ) ≡ L (N ) to the vector of function values on the spectral grid, (Lf (N ) )(x m ) = n L (N ) mn f n . Solving for (the spectral approximation to) the solution of the differential equation (C.1) then reduces to the standard algebraic problem of solving a finite system of linear equations.

C.1 Explicit expressions
Analytic expressions for cardinal functions and differential matrix components, for many different sets of basis functions, may be found in appendix F of Ref. [21]. For a Chebyshev basis and the Chebyshev-Gauss-Lobatto grid (C.7), cardinal functions satisfying (C.5) are given by where T k (x) denote Chebyshev polynomials of the first kind and c j ≡ 1 for 0 < j < N while c 0 = c N ≡ 2. The interior grid points lie at extrema of T N (x). Derivatives of these cardinal functions, evaluated on the Gauss-Lobatto grid, can be evaluated explicitly.
For the first derivative one finds [21] (D ) k . For the Fourier grid with endpoint (C.8), cardinal functions can be expressed as 13) and the first two differentiation matrices are given by A linear transformation, y = ax + b, may be used to convert the above expressions into forms appropriate for arbitrary finite intervals.

C.2 Domain decomposition
The error in an N -term spectral approximation to some function u decreases exponentially with increasing N , provided u satisfies appropriate analyticity conditions [21]. In practice, this desirable behavior only holds as long as there is negligible round-off error from finite precision numerical arithmetic. Unfortunately, differentiation matrices become increasingly ill-conditioned as N increases and this leads to progressively worsening numerical errors in the eventual solution of the linear system. Moreover, with Chebyshev grids, the spacing between grid points is non-uniform and near the endpoints of the interval the grid spacing decreases as 1/N 2 . This rapid decrease of grid spacing can lead to short wavelength (so-called 'CFL') instabilities in time evolution problems. These difficulties can be alleviated by partitioning the computational domain into multiple subdomains, inside each of which one constructs an independent spectral approximation. This is known as domain decomposition. In effect, one solves the differential equation of interest independently in each subdomain with boundary conditions which enforce appropriate continuity conditions connecting adjoining subdomains. Differentiation matrices for the entire domain become block-diagonal.
For a one dimensional problem, if one partitions the full domain into M subdomains, and uses an N -point Chebyshev grid containing endpoints within each subdomain, then each interior subdomain boundary will appear twice in the resulting complete list of grid points (3.3). For a second order differential equation, before solving the resulting linear system, L (M N ) f (M N ) = g (M N ) , one simply replaces each pair of rows which represent the same interior subdomain boundary by a near pair of linear equations which encode continuity of the function, For further detail refer to Ref. [21].

D.1 Longitudinal filter
Numerically filtering the propagating data, namely the functions {B, a, f, λ}, to remove small amplitude noise, specifically cutoff-scale rapid variations in the longitudinal direction, is essential to achieve stable time evolution with low background energy density, especially for narrow shock collisions where a very fine longitudinal grid is required. The reasons behind this, involving spectral blocking in non-linear equations, are discussed in Refs. [2,21]. Such filtering must be applied carefully. To maintain consistency of the solution of the nested set of Einstein equations (2.8), we only filter at the end of each time step, not within RK4 substeps and not in between solutions of the nested radial differential equations.
There are many ways to implement a low-pass filter. For our periodic functions of z, we use a smooth multiplicative filter in k-space. A cardinal function representation using the uniform grid (C.8) and Fourier cardinal functions (C.13), φ(z) = 2N −1 j=0 φ j C j (z), is exactly equivalent to a truncated Fourier series, φ(z) =  We suppress the amplitude of modes with large |k| by multiplying the Fourier coefficients {φ k } by the filter function The parameter Λ is the fractional bandwidth of the filter while N δ/(2π) is the characteristic width in wavevectors of the filter roll-off. We chose to use Λ = 2/3 and δ = 1/2. Transforming back to real space produces the smoothed function

D.2 Radial filter
After transforming single shock solutions to infalling coordinates, as described in Appendix B, we found that moderately high derivatives of the resulting anisotropy function, such as ∂ 3 u ∂ 3 z b(u, z, t 0 ), after the longitudinal filtering as just described, would still show visible noise with rapid radial and longitudinal variation. Such noise grows and becomes problematic upon time evolution. To suppress such artifacts in the initial data, we perform radial filtering on the initial anisotropy function in a matter designed to suppress radial noise near the boundary while simultaneously ensuring the correct near boundary asymptotic behavior.
We do this by first constructing, analytically, the near boundary expansion of the transformation functions solving Eqs. (2.28)-(2.29c), and thence the resulting metric anisotropy function B(u, z, t 0 ) via Eq. (2.26) [or equivalently the rescaled function b(u, z, t 0 ) defined in Eq. (3.12a)]. Explicit expressions for these near-boundary expansions appeared in Appendix B.1. Let b (K) (u, z, t 0 ) denote the K-term partial sum of the near-boundary expansion for the rescaled anisotropy function b(u, z, t 0 ). We define a correction function ∆ (K) (u, z) by the condition that for m = 0, · · ·, K−1, while simultaneously requiring that ∆ (K) (u, z), evaluated on the radial grid, is only non-zero on the first K radial grid points closest to the boundary. On the left side of condition (D.6), the radial derivatives are evaluated using the spectral derivative matrix D u applied to the list of values of b and ∆ (K) on the radial grid. These conditions uniquely determine the correction function ∆ (K) (represented on the spectral grid). The corrected function b improved ≡ b + ∆ (K) coincides with the input function b away from the boundary (by more than K grid points), while having corrected values of radial derivatives up through order K−1 at the boundary. Choosing K = 7, we find that this procedure is effective in suppressing numerical noise in initial data up to quite high orders of derivatives in both radial and longitudinal directions. Unlike a conventional filter, the effect of this procedure is restricted to a small region near the boundary.

E Runge-Kutta methods
Given a first order differential equation for some R k -valued function Φ(t), with initial condition Φ(t 0 ) = Φ 0 , the standard fourth order Runge-Kutta (RK4) algorithm iteratively constructs an approximate solutionΦ at times t n ≡ t n−1 + δt, via the recursion relationΦ (t n+1 ; δt) ≡Φ(t n ) + δt where K j (t n ) ≡ F t n + α j δt,Φ(t n ) + α j δt K j−1 (t n ) , (E.3) withΦ(t 0 ; δt) = Φ 0 . The coefficient vectors defining the RK4 "substeps" are given by To convert this method to an adaptive stepsize integration method, one needs a local error estimation, i.e., some estimate of the difference betweenΦ(t n ) and the desired solution Φ(t n ), assuming thatΦ(t n−1 ) is correct, together with an algorighm for decreasing or increasing the time step δt based on this error estimate. The easiest way to achieve this is to compare the results of performing a single RK4 step with timestep δt versus two RK4 steps with timestep δt/2. The latter (more time consuming) calcuation will suffer from less error due to timestep discretization and, if δt is sufficiently small, this difference will be a decent approximation to the deviation from the true solution. We define err(t+δt) = Φ (t+δt; δt) −Φ(t+δt; δt/2) , (E.5) given a common starting value at time t. For the choice of norm, we use an L ∞ norm, or the maximum over all components ofΦ. If the goal of the numerical calculation is to achieve a relative precision of 10 −a , then we adjust the time step according to δt n+1 = δt n 10 −a err 1/4 .

(E.6)
The "learning rate" of this adaptive algorithm is governed by the exponent 1/4 in this rule. This value reflects the fact that in the basic RK4 method, the error scales as (δt) 4 for sufficiently small timestep δt. In our code we did not impose minimum or maximum step sizes. And in our specific application of transformation to infalling coordinates, when starting with an initial step size of δu = 0.0001 it turned out to be sufficient to update the step size using Eq. (E.6) and always advance directly to the next slice without further adjustments. More generally, it can be necessary to reject a trial step and repeat the the calculation with a smaller step size if the initial error exceeds the desired limit.