Prepared for submission to JHEP Holography and off-center collisions of localized shock waves

Using numerical holography, we study the collision, at non-zero impact parameter, of bounded, localized distributions of energy density chosen to mimic relativistic heavy ion collisions, in strongly coupled N = 4 supersymmetric Yang-Mills theory. Both longitudinal and transverse dynamics in the dual field theory are properly described. Using the gravitational description, we solve 5D Einstein equations with no dimensionality reducing symmetry restrictions to find the asymptotically antide Sitter spacetime geometry. Implications of our results on the understanding of early stages of heavy ion collisions, including the development of transverse radial flow, are discussed.


Introduction
The recognition that the quark-gluon plasma (QGP) produced in relativistic heavy ion collisions is strongly coupled [1] has prompted much work using gauge/gravity duality (or "holography") to study aspects of non-equilibrium dynamics in strongly coupled N = 4 supersymmetric Yang-Mills theory (SYM), which may be viewed as a theoretically tractable toy model for real QGP. Work to date has involved various idealizations (boost invariance, planar shocks, imploding shells) which reduce the dimensionality of the computational problem, at the cost of eliminating significant aspects of the real collision dynamics [2][3][4][5][6][7]. In this paper, we report the results of the first calculation of this type which does not impose unrealistic dimensionality reducing restrictions. We study the collision, at non-zero impact parameter, of bounded, localized distributions of energy density which mimic colliding relativistic Lorentz-contracted nuclei and form stable incoming projectiles in strongly coupled SYM.
In the dual gravitational description our initial state consists of two localized incoming gravitational waves, in asymptotically anti-de Sitter (AdS) spacetime, which are arranged to collide at a non-zero impact parameter. The precollision geometry contains a trapped surface and the collision results in the formation of a black brane. We numerically solve the full 5D Einstein equations for the geometry during and after the collision and report on the evolution of the SYM stress-energy tensor T µν .

Gravitational formulation 2.1 Single shocks
We construct initial data for Einstein's equations by combining the metrics describing gravitational shock waves moving at the speed of light in opposite directions. In Fefferman-Graham (FG) coordinates, the metric of a single shock moving in the ±z direction is given by [8,9] . This is a sourceless solution to Einstein's equations with cosmological constant Λ ≡ −6/L 2 , provided the function h ± satisfies the linear differential equation Solutions to this equation which vanish at the boundary, s = 0, may be written in the form where H ± is an arbitrary function of the 2D transverse wavevector k ⊥ and the longitudinal variable z ∓ (and I 2 is a modified Bessel function). The geometry described by eqs. (2.1) and (2.3) represents a state in the dual SYM theory with a stress-energy tensor expectation value given by [10].
with all other components vanishing. Here, H ± is the 2D transverse Fourier transform of H ± , and the constant κ ≡ L 3 /(4πG N ) = N 2 c /(2π 2 ) (with N c the gauge group rank of the SU (N c ) SYM theory). In other words, the function H ± specifies the energy density (and longitudinal stress and momentum density) of a shock wave moving, nondispersively, at the speed of light in the ±z direction. Given any choice of this energy density profile, eqs. (2.1) and (2.3) give an explicit form for the unique dual gravitational geometry which describes this shock wave.
For simplicity, we consider Gaussian energy density profiles, with longitudinal width w, transverse width R, and a transverse offset ±b/2. Hence, b will be the impact parameter when these oppositely directed shock waves collide. The amplitude A gives the maximum value of the longitudinally integrated energy density (divided by κ).
For the explicit computations presented below, we adopt units in which the amplitude A equals unity and choose longitudinal and transverse widths w = 1 2 and R = 4, respectively, and impact parameter b = 3 4 Rx. Consequently, the stress tensor (2.4) describes localized lumps of energy centered about x = ±b/2, y = 0, and z = ±t. The AdS curvature scale L is not a physical scale in the dual QFT and may independently be set to unity.

Infalling coordinates
Our time evolution scheme [11] for asymptotically-AdS gravitational dynamics uses infalling Eddington-Finkelstein (EF) coordinates in which the spacetime metric has the form with Greek indices denoting spacetime boundary coordinates, x µ ≡ (t, x, y, z), and r the bulk radial coordinate. For the asymptotically-AdS geometries of interest, the metric coefficients g µν have the near-boundary behavior g µν = η µν + g (4) µν /r 4 + O(1/r 5 ) as r → ∞, with η µν ≡ diag(−1, +1, +1, +1) the usual Minkowski metric tensor. The sub-leading coefficients g (4) µν determine the SYM stress tensor expectation value. In these coordinates, 00 . (2.7) Although the single shock solution to Einstein's equations has the nice analytic form (2.4) in Fefferman-Graham coordinates, this geometry does not have a simple analytic form in infalling EF coordinates; the transformation to infalling coordinates must be performed numerically. One may easily show that, in infalling coordinates, curves along which r varies, with the other coordinates held fixed, are infalling null geodesics (with r an affine parameter). To transform the geometry (2.1) to the infalling form (2.6) one must locate the same congruence of infalling radial null geodesics in FG coordinates. Let Y ≡ {y µ , s} denote the FG coordinates of some event, and let X ≡ {x µ , r} denote the EF coordinates of the event at affine parameter r along the radial infalling geodesic which begins at boundary coordinates x µ . Then the solution Y (X) to the geodesic equation (in FG coordinates) for the same null geodesic which begins at boundary coordinates x µ provides the required mapping between EF and FG coordinates. Given this mapping, the required transformation of the FG metric components G M N (Y ) to the metric components G M N (X) in our infalling coordinates To compute the congruence Y (X), we first periodically compactified spatial directions (to obtain a finite computational domain) with transverse size L x = L y = 32 and longitudinal length L z = 12. We employed spectral methods [11], and used a rectangular grid built from single domain Fourier grids with 32 points in the transverse directions, a 256 point Fourier grid in the longitudinal direction, and three Chebyshev domains of 32 points each in the radial direction. We used a Newton iterative procedure to solve the non-linear geodesic equation starting at each grid point on the boundary (modulo the C 4v transverse cubic symmetry). 1 We integrated the geodesic equations to a depth of s = 5 and fixed the residual radial shift reparameterization freedom [11] by demanding that the surfaces u ≡ 1/r = 5 in infalling coordinates, and s = 5 in FG coordinates, coincide.

Colliding shocks: initial data
For early times, t −w, the Gaussian profiles H ± of the oppositely directed incoming shocks have negligible overlap and the precollision geometry can be constructed by replacing the last term in the single shock metric (2.1) with the sum of corresponding terms from left and right moving shocks. The resulting metric satisfies Einstein's equations, at early times, up to exponentially small errors.
Initial data for the subsequent time evolution consists of the values of the spatial metric coefficients, rescaled to have unit determinant,ĝ ij ≡ g ij /(det g ij ) 1/3 , on every point of the initial time slice, together with the values of the energy and momentum density (or equivalently the asymptotic coefficients g (4) 0ν ) at the initial time. The initial time slice was chosen to lie at t = −2, We slightly modified the initial data by adding a small uniform background energy density, equal to 3.7% of the peak energy density of the incoming shocks. This modestly displaced the location of the apparent horizon toward the boundary, improving numerical stability and allowing use of a coarser grid, thereby reducing memory requirements. 2

Colliding shocks: time evolution
Time development of the geometry was calculated using the procedure described in detail in ref. [11]. Time evolution was performed using a spectral grid of size N x = N y = 39, N z = 145, and N r = 40. A fourth-order Runge-Kutta time integration algorithm was used with a time step of 0.005. Integration continued to a final time of t = 4. Even though our gravitational dynamics is five dimensional, with no symmetry restrictions imposed, we were able to perform the time evolution on a 6 core desktop computer with a 14 day runtime.

Results
In Fig. 1 we plot the energy density T 00 (top) and energy flux T 0i (bottom) in the plane y = 0 at several values of time. 3 Note that the color scaling varies from plot to plot. The color scaling in the lower plots denotes |T 0i | and the flow lines indicate the direction of T 0i . At time t = −2 the system consists of two well separated lumps of energy moving towards each other at the speed of light. The non-zero impact parameter in the x-direction is apparent. At time t = 0, when both shocks are centered at z = 0, the energy density and flux differ by only 8% and 15%, respectively, from a linear superposition of two unmodified shocks. Nevertheless, the subsequent evolution is very different from two unmodified outgoing shocks: the remnants of the initial shocks, which remain close to the lightcone, z = ±t, are significantly attenuated in amplitude with the extracted energy deposited in the interior region.
At sufficiently late times, in accord with fluid/gravity duality [12,13], the evolution of the stress tensor should be governed by hydrodynamics. To compare to hydrodynamics, we define the fluid velocity to be the normalized time-like (u µ u µ = −1) future directed eigenvector of the stress tensor, with the proper energy density. Given the flow field u µ and energy density , we construct the hydrodynamic approximation to the stress tensor T µν hydro using the constitutive relations of first order viscous hydrodynamics.
In Fig. 2 we plot the stress tensor components T xx and T zz , and their hydrodynamic approximations T xx hydro and T zz hydro , at the spatial origin x = y = z = 0, as a function of time. At this point the stress tensor is diagonal, the flow velocity u = 0, and T xx and T zz are simply the pressures in the x and z directions. As shown in the figure, the  pressures increase dramatically during the collision, reflecting a system which is highly anisotropic and far from equilibrium. After a time t ≈ 1.25 the pressures are well described by viscous hydrodynamics. Remarkably, at this time the transverse pressure T xx is nearly ten times larger than T zz . (This latter phenomena has also been seen in 1 + 1 dimensional flow [2][3][4].) To quantify the domain in which hydrodynamics is applicable, we define a residual measure ∆ ≡ (1/p) ∆T µν ∆T µν , with ∆T µν ≡ T µν − T µν hydro andp ≡ /3 the average pressure in the local rest frame. The quantity ∆ is frame-independent but, when evaluated in the local fluid rest frame, reduces to the relative difference between the spatial stress in T µν and T µν hydro . Regions with ∆ 1 are evolving hydrodynamically. In Fig. 3 we plot ∆ in the transverse plane at proper times τ = 1, 1.25, and 2, and rapidities ξ = 0 and 1. The color scaling is the same in all plots. Focusing first on ξ = 0 (top row), at τ = 1 one sees that ∆ 0.5 in the central region (x, y ≈ 0), and hydrodynamics is not a good description. However, by τ = 1.25 a fluid droplet with ∆ 0.15 and transverse radius x ⊥ ≡ |x ⊥ | 5.3 has formed, with subsequent evolution well described by hydrodynamics. At τ = 2 the transverse size of the droplet has increased and ∆ < 0.15 for x ⊥ 8.6. Turning now to the behavior at rapidity ξ = 1 (bottom row), one sees that for small x ⊥ the system is already evolving hydrodynamically at τ = 1. Moreover, the onset of hydrodynamics occurs earlier for x < 0 than for x > 0. This feature reflects the fact that the receding maxima remain far from equilibrium and non-hydrodynamic, and (as seen in Fig. 1), the maxima with ξ > 0 lies at x > 0. Interestingly, the inclusion of transverse dynamics seems to hasten the approach to local equilibrium: the equilibration time t hydro ∼ 1.25 is about 30% smaller than was the case in our previous studies [3,11] of planar shock collisions. Recent work [14,15] has found that equilibration time scales of far-from-equilibrium states can be understood, at least semi-quantitatively, in terms of the spectrum of quasinormal modes. Postcollision, a distribution of quasinormal modes will be excited. With the inclusion of transverse dynamics, the dominantly excited quasinormal modes will be ones with non-zero transverse wavevectors of order 1/R. Faster decay rates of quasinormal modes with increasing |k ⊥ | should lead to faster apparent relaxation of the entire sum (i.e., the metric perturbation) independent of location in the transverse plane. This will be the case provided one considers sufficiently coarse measures of relaxation, such as our ∆ < 0.15 criterion, which can become satisfied when numerous quasinormal modes are comparably excited. With a much more stringent local equilibration criterion, we  would expect to see less of a difference relative to the planar case.
A striking feature of the the post-collision evolution in Fig. 1 is the appearance of flow in the transverse plane at early times. The early-time acceleration imparted on the transverse flow can have a significant impact on the subsequent transverse expansion. In Fig. 4 we plot the fluid 3-velocity v ≡ u/u 0 in the z−x, z−y, and x−y planes at time t = 4. The color scaling, which indicates |v|, is the same in each plot. The flow lines show the direction of v. Regions in which ∆ > 0.15, and the system is not behaving hydrodynamically, have been excised. Already at time t = 4 and radius x ⊥ ≈ 5 the transverse fluid velocity in the x−y plane has magnitude 0.3. In contrast, the maximum of the longitudinal velocity, which occurs in the neighborhood of the receding maxima, is 0.64.
One sees from Fig. 4 that the fluid velocity in the x−y plane is nearly radial: we see no strong signatures of elliptic flow. At vanishing rapidity (or z = 0), the traceless part of the spatial stress, as well as the hydrodynamic residual shown in the upper row of fig. 3, deviate significantly from rotational symmetry in the x−y plane (at a 5-6% level at t = 4). However, the energy density and the resulting fluid velocity deviate only slightly from azimuthal symmetry; at t = 4 their transverse plane anisotropy is about 1%. Our computed evolution does not extend through the entire hydrodynamic phase of the plasma and continued hydrodynamic expansion may generate more elliptic flow. However, the very modest elliptic flow, despite the substantial impact parameter, may reflect an unphysical aspect of our choice of Gaussian initial energy density profiles. Gaussian profiles differ, of course, from more realistic models of nuclear density (such as Wood-Saxon profiles). But a Gaussian profile in the transverse plane is computationally convenient. The very rapid decrease of its Fourier transform allowed us to use a rather modest transverse grid. However, a special feature of Gaussian transverse plane profiles is that the overlap function (i.e., the product of the energy densities of the two incident projectiles) is rotationally symmetric for any impact parameter: With Gaussian profiles (and equal size incident projectiles), there simply is no "almond" in the overlap function. The actual dynamics depends, of course, on the full spacetime dependence of the initial data, not just on this overlap function and, as noted above, the stress tensor and the residual ∆ depart significantly from transverse rotational symmetry. Nevertheless, a simple picture in which the early transverse flow reflects, in large part, this initial overlap function appears plausible. A notable feature in our results is that the fluid flow, in the z−x plane, is not symmetric about the z axis and the longitudinal flow does not vanish at z = 0. The latter observation is a direct violation of the simplified model of boost invariant flow, in which v z = z/t and vanishes at z = 0. Nevertheless, at t = 4 and in the region of space where > 0.6 max( ), the longitudinal flow is roughly described by boost invariant flow at the 20% level or better, with larger deviations appearing at larger rapidities. 4 For planar shock collisions, the deviation of the longitudinal fluid velocity from boost invariant flow decreases as the shock thickness decreases [4,11]. 5 It will be interesting to see if this also holds when transverse dynamics is included.
We conclude by discussing the early-time transverse flow predicted in ref. [16]. There, using assumptions of boost invariance and transverse plane rotational symmetry, it is argued that at early times the transverse energy flux is proportional to the gradient of the energy density and grows linearly with time, In Fig. 5 we plot the angular averaged radial flow T 0⊥ ≡ x i ⊥ T 0i , together with the approximation (3.4), at proper times τ = 1.25 and 2, and rapidities ξ = 0 and 1. The approximation (3.4) works remarkably well at both times and rapidities, although the agreement is not quite as good at ξ = 1 where the assumption of boost invariance is more strongly violated. It would be interesting to see if the agreement with the approximation (3.4) improves when the shock thickness decreases.

Final remarks
We have presented results from the first calculation, using numerical holography, of the evolution in strongly coupled N = 4 SYM theory of an initial state which resembles two colliding heavy nuclei, Correctly treating the dynamics, both longitudinal and transverse, requires solving a five dimensional gravitational initial value problem. This was challenging, but feasible, using our characteristic formulation and (good) desktop-scale computing resources. The results presented above are a proof-of-principle. Clearly, it will be interesting to explore how results change as the impact parameter, longitudinal thickness, and transverse size of each projectile are varied. It should also be possible to explore the effects produced by using more realistic initial energy density profiles, possibly including shorter distance fluctuations of the type which have been suggested to be relevant in real heavy ion collisions. (See, for example, ref. [17][18][19][20] and references therein.) We hope future work will shed light on some of these topics.