Towards collisions of inhomogeneous shockwaves in AdS

We perform a numerical simulation of the evolution of inhomogeneities with transverse profile in a collision of gravitational shockwaves in asymptotically anti-de Sitter spacetime. This constitutes a step closer towards an accurate holographic description of the thermalization of a strongly coupled plasma, which can model the dynamics of heavy ion collisions. The results indicate that the considered inhomogeneities typically become hydrodynamical earlier or at the same moment when hydrodynamics applies to the background, even though they decay slowly.


Introduction
Quantum Chromodynamics (QCD) is a challenging theory in its non-perturbative regime, especially when out-of-equilibrium situations are involved. In particular, the formation of a Quark-Gluon Plasma (QGP) in heavy ion collisions, such as the ones carried out at the RHIC and LHC particle colliders, still lacks a satisfying theoretical description.
The description of the initial stage of this process (right after the collision) would require a full non-perturbative calculation of strongly coupled QCD. Since this is currently beyond our reach, the necessity and importance of developing strong coupling techniques for heavy ion phenomenology has received a boost in recent years. For instance, much work has been done using viscous hydrodynamics [2] or extrapolations from weakly coupled calculations [3]. Although each technique comes with limitations of its own, some successful predictions include those of the shear viscosity to entropy density ratio [4] and the elliptic flow coefficient v 2 [5].
Here we shall take advantage of the convenience of applying the gauge/gravity duality, which allows to access the non-trivial quantum dynamics of this out-of-equilibrium stage by solving classical gravitational dynamics. In this context, our particular holomodel will be constructed for the simplest theory with a gravity dual, N = 4 super Yang-Mills (SYM) theory (in the N c , λ → ∞ limit) [6]. Using this framework, several authors [7][8][9][10][11][12] have previously translated the aforementioned problem into that of a collision of gravitational shock waves in an asymptotically anti-de Sitter (AdS 5 ) spacetime. Such an approach would be equivalent to studying the collision of infinitely extended, homogeneous and planar layers of matter in SYM, which in turn would simulate the highly Lorentz-contracted colliding ions. For an excellent review on this topic, see [13].
These studies have enlightened ostensibly the longitudinal dynamics involved in heavy ion collisions. However, in the actual experiments, the presence of inhomogeneities and JHEP07(2015)126 the build-up of momentum in the transverse plane could make transverse dynamics important. In fact, a study including radial flow [14] showed that the momentum distribution reaches local equilibrium quickly, after which hydrodynamics applies. Boost-invariance and rotational symmetry were assumed as approximations, in order to keep the numerical calculation effectively 2+1 dimensional. 1 Here, a different formulation will be used to the same effect.
Particularly, as a first approach towards a full calculation including transverse dynamics, we consider the propagation of inhomogeneous perturbations on top of the dynamical gravitational background that encodes the process of thermalization. The dependence on the transverse direction is specified by fixing it to be that of a planar wave. In other words, in each simulation we consider the evolution of a specific Fourier mode.
This problem requires the use of numerical techniques in order to obtain quantitative results. The connection with the dynamics of the plasma is made by extracting the evolution of the stress-energy tensor. Generically, we find that the inhomogeneities decay slowlier than the background thermalizes, which justifies the inclusion of transverse dynamics for a more accurate calculation. However, the thermalization time, defined according to the applicability of hydrodynamics, is not affected, since the inhomogeneities acquire a hydrodynamic behavior soon enough.

Gravitational description
The ansatz for the 5-dimensional spacetime metric is a generalization of that in [7], including additional components to account for the transverse directions x 1 , x 2 in the form where A, B, C, D, Σ, F and G are functions of the bulk radial coordinate (r), time (t) and spatial longitudinal (y) and transverse (x 1 , x 2 ) coordinates. The boundary is located at r = ∞. Note that the determinant of the spatial part of the metric is a power of a single function, Σ. This allows to simplify Einstein's equations, by following the characteristic formulation of General Relativity (GR). For the same reason, we employ generalized ingoing Eddington-Finkelstein coordinates, where paths of varying r, with the other coordinates fixed, are infalling radial null geodesics. As a matter of fact, the metric ansatz is invariant under arbitrary reparametrizations of this parameter, r → r + ξ(t, y, x). This constitutes a gauge freedom that will be fixed by placing the Apparent Horizon (AH) at r = 1. The ansatz is complemented with h(r, t, y, x 1 ) = h 0 (r, t, y) + e ikx 1 δh(r, t, y) where h represents every function in the metric. This is so that the problem simplifies from 3+1 to 2+1 dynamics. The δh terms will be treated as perturbations, and Einstein's equations 2 will be linearized around the background solution of [7]. Note that there is no

JHEP07(2015)126
background counterpart for the functions C, D and G. The value of k must be fixed from the beginning. These linearized equations are too long to be reproduced here but they can be found in [1]. In order to be able to organize the equations in a favorable structure, it is important to write them in a fully covariant way using derivatives along outgoing null raysḣ, as defined inḣ Note that d 3 h is a derivative in the longitudinal direction orthogonal to radial null geodesics. Taking into account the previous decomposition into background and fluctuations, these definitions apply to the fluctuations aṡ The asymptotic analysis of Einstein's equations near the boundary provides a large r expansion, where we include the extra gauge freedom ξ described above, of the form We identify a 4 , δa 4 , b 4 , δb 4 , δc 4 , δd 4 , f 4 , δf 4 , δg 4 as the normalizable modes which are related to the stress-energy tensor of the dual theory. They are functions of (t, y), and they are not completely independent, since the previous expansions solve the equations only as long as the background coefficients satisfy and the inhomogeneities' coefficients satisfy

JHEP07(2015)126
These equations will be used to evolve the boundary conditions forward in time. They can equivalently be derived without making use of Einstein's equations, but from the conservation equations of the stress-energy tensor ∇ µ T µν = 0. They have a physical interpretation in terms of continuity conditions for the transport of energy and momentum. Specifically, the expectation value of the stress-energy tensor of this problem contains energy density, momentum densities, pressures, and shear stress terms, all holographically mapped via gauge/gravity duality: After transforming the asymptotic expansions (2.5) to Fefferman-Graham coordinates, we can use the holographic renormalization prescription to extract the relations where we have omitted the e ik x 1 factors in front of every δ term.

Numerics overview
A generic description of the numerical approach that can be applied to solve the dynamics of this problem is found in [16]. By applying the characteristic formulation within AdS, the set of coupled partial differential equations of GR can be very conveniently written as a nested set of linear ordinary differential equations. It is necessary to specify the spatial part of the metric on the initial time slice, except for the determinant (that is, except for Σ). Thus, one starts with the initial data provided for F in = {B 0 , δB, δC, δD}. From there, the nested structure allows to solve for the other functions step by step, following the sequence Note that the dotted functions are solved as if they were unrelated to their undotted counterparts. These 14 steps correspond to the 6 equations for the background, followed by the 8 linearized equations for the inhomogeneities. The numerical scheme that was implemented was based on pseudospectral methods, solving these equations at every time slice, then inverting (2.3)-(2.4) to obtain time derivatives and evolving the F in forward to the next slice. An Adams-Bashforth method was carried out for this purpose. 3 In addition, there are 4 constraints (equations that do not provide dynamics, but are evaluated to monitor accuracy): one is of the background and includes ∂ tṠ , while the other three include ∂ t (δ S), ∂ r (δ F ) and ∂ r (δ G) respectively. Their asymptotic analysis near the boundary provides the conditions (2.6)-(2.8), which must be imposed to provide boundary conditions as the time evolution goes along. Given their self-fulfilling nature, the constraints can be left out of the calculation, only to be used as a convergence check of the numerics.
The initial data of the evolution is extracted from the metric of two planar shocks moving towards each other. For the calculations presented here, the same shocks as in [7] were considered, H(t, y) ≡ µ 3 √ 2πw 2 e −(t∓y) 2 /2w 2 , with w = 0.75/µ. We also chose a back- ground energy density δ = 0.075µ 4 . This provides initial data for B 0 (t = 0, r, y), and for the coefficients a 4 (t = 0, y) and f 4 (t = 0, y). But now this must be supplemented with initial data for the perturbations. In principle, any initial state can be considered, as long as Einstein's equations are satisfied (the inhomogeneities may take any shape).
In our calculations, we simply chose δa 4 , δf 4 , δg 4 so that the inhomogeneity behaves like a planar wave proportional to the amplitude of the background at each point, that is, a 4 → a 4 (1+ e ik x 1 ), f 4 → f 4 (1+ e ik x 1 ) and δg 4 = 0. And for the bulk profiles of δB, δC, δD, we chose them to be given by the first terms in their respective expansions (2.5). This fixes the radial dependence and the boundary values are determined by (2.8) so that δB(0, r, y) = a 4 (0, y) 4r 4 , δC(0, r, y) = δD(0, r, y) = 0 . For these initial conditions, functions δC, δD and δG acquire non-vanishing profiles spontaneously. It would be interesting to study different choices of initial conditions, in order to check how they affect the stability of the propagation of the shocks before the collision, but we choose to leave such an analysis to future work. Several inhomogeneities of the expected stress-energy tensor are plotted in figures 1-3. Note that since the equations are linearized, the overall amplitude of these inhomogeneities is completely irrelevant. Their sign is also irrelevant, since each of these figures corresponds to an x 1 = ctant. slice, and they oscillate along the transverse direcion. This is due to the factor e ikx 1 they bear in front.

Apparent horizon
As discussed in [16],the residual gauge freedom r → r + ξ(t, y), is fixed by imposing the AH to lie at a constant r, for instance r = 1. This is easily carried out by absorbing any deviation into the chosen gauge, that is, δξ = r AH − 1. The computation behind is a crucial part of the numerical calculation, so it merits to give a further explanation about this.

JHEP07(2015)126
To find the position of the AH, r AH (y), [7] gives where the functions here can be those of the background (the inhomogeneities make a negligible contribution). This equation can be solved by finding the root of that expression. However, it is derived under the assumption that the AH lie at a constant position r = r AH , instead of a trajectory r = r AH (y). Otherwise, (4.1) must be modified by inverting the gauge transformation that would have led us to having it at a constant position in the first place.
As can be seen in its expansion (2.5b), this entails the explicit change F → F − ∂ y ξ, as well as evaluating every function at r + ξ. This is a significant complication, since the problem becomes an intricate non-linear differential equation for r AH (y) (or, equivalently, ξ(y)). Given a time slice t 0 , (4.1) gives the correct AH after performing an iterative procedure to find ξ(t 0 , y). But during the time evolution, it is more efficient to demand the time derivative to vanish, It can be seen [1] that this constitutes a linear 2 nd order differential equation for ∂ t ξ, so its time evolution can be readily computed. But, since we are using the equation that assumes a constant r AH , this approach introduces some error in the calculation, which must be corrected by performing the explicit calculation of ξ(t 0 , y) occasionally during the evolution (every 10-20 timesteps).

Hydrodynamics
Comparing our results with those for the background (see [7]), a common observation can be drawn: the presence of the inhomogeneities persists even after the background has equilibrated. For instance, we can see that the energy density E is spread out by the time t ∼ 10/µ. However, we find that there is still a very uneven profile in δE well after that time, as shown in figure 1. This is an interesting observation, because the presence of such a perturbation in the energy density at later times would give some non-trivial structure to the otherwise flat spatial profile of the local energy density [18]. However, it is important to keep in mind that the program of research on holographic thermalization has created two well established concepts of what should be understood as thermalization. On the one hand, one could refer to the isotropization of the stress tensor in the local rest frame. This can take a very long time. 4 On the other hand, one could refer to the applicability of viscous hydrodynamic constitutive relations. This is sometimes called "hydrodynamization", and the consensus in this respect is that it is surprisingly short. On the following, we shall refer to the latter definition of thermalization. Thus, in order to draw comparisons with previous work, we need to test the validity of hydrodynamics. To do so, we compare the actual pressures from the boundary stressenergy tensor with the pressures that would follow if the viscous hydrodynamic constitutive relations were satisfied [19,20]. Results from such a comparison can be found In figure 4, where we plot the inhomogeneities in the longitudinal δP y and transverse {δP x 1 , δP x 2 } pressures as a function of time, at two specific points, y = y 0 = 5π/µ (this is the point where our shocks actually collide) and y = y 0 +3/µ. The dashed lines show the fluctuations in the pressures δP hydro as predicted by the hydrodynamic equations.
The agreement with hydrodynamics is quite remarkable. In particular, at y = y 0 , the sudden increase in the pressures due to the collision is reflected in a dramatic rise in their inhomogeneities, as expected. During this stage, the system is very anisotropic and far from equilibrium, and as a consequence hydrodynamics is not expected to be applicable (as reflected in [7] for the background). Surprisingly, the hydrodynamic constitutive relations seem to hold almost from the beginning for the fluctuations considered here.
At late times, the pressure inhomogeneities asymptotically approach each other. This process of isotropization has not been completed even at t = 12/µ, which is much larger than the time it takes for the background to become isotropic. Thus, our fluctuations provide a nice example of a system where isotropization time and the hydrodynamization time are completely different.

Discussion
The transverse dynamics in a collision of gravitational shockwaves in AdS has been previously considered in [14], where the longitudinal dynamics were approximated as boost-

JHEP07(2015)126
invariant, and rotational symmetry is assumed. The dynamics we have discussed here do not require boost-invariance and describe the behavior of inhomogeneities in the shockwave which are allowed to propagate in a transverse direction. Our main approximations are the linearization of the equations of motion (2.2), and the simple choice of initial conditions (3.1).
The results apply qualitatively to all strongly coupled 4D conformal gauge theories with a gravitational dual description. Typically, a closer approach to QCD requires introducing quarks, or fundamental matter. However, the gluons are the dominant degrees of freedom at the timescales involved in a heavy ion collision, and this is what motivates the idealization of considering nothing more than pure gravity. Still, tt is not obvious a priori how good an approximation to the QGP dynamics it can provide. As usual, eventual comparison with experimental observations is what can support the assumptions taken in holographic models.
At the experiments of the RHIC and LHC, density perturbations may exist and play a relevant role in the results. In order for the holographic approaches to make contact with this, inhomogeneities need to be included. Our disposal (2.2) is limited, but also a first step towards a complete calculation. A different but also relevant phenomenon is the presence of radial flow, for which an expanding fireball of finite size would need to be considered, as opposed to a wave of infinite transverse extent. In order to model this, one would need to go beyond the fluctuation approximation.
Based on the typical behavior manifested by our fluctuations (as shown in figure 4), we infer that one is to expect inhomogeneous profiles to keep being inhomogeneous for a relatively long time. Quite longer than the thermalization time, which would be unaffected by the inhomogeneities, given that they experience a very fast hydrodynamization. This result supports the use of hydrodynamic approximations to study the propagation of perturbations during the out-of-equilibrium stage of QGP formation. However, one should note that the agreement observed in our results can be related to the linearized treatment we have used to solve the equations.
Furthermore, interesting information could be extracted from a thorough analysis of the results. In particular, considering different initial configurations for the inhomogeneities would allow to establish which perturbations are physically relevant and which are not. It would also be possible to see at which point of the evolution do they become hydrodynamical, or whether they do at all. Additionally, a byproduct of this calculation would be the spectrum of quasinormal modes for finite k of the final black hole. A hydrodynamic gradient expansion of the dynamics could allow to read them off [21]. We leave this to future work.
Another estimulating extension of this project would be to allow for an interaction between modes with different momenta, since in that case turbulent effects may arise. In recent years, it has been discovered that turbulence is an ubiquitous property of Gravity [22,23]. However, if symmetries are forced into the system, turbulence is missed. Intuitively, one can argue that turbulent effects would lead to a shorter thermalization time, due to the cascading behavior towards higher modes (both in the transverse and longitudinal directions).

JHEP07(2015)126
Finally, it is possible that turbulence may arise in the calculation presented here at later times, since the perturbations propagate in a non-trivial time dependent background. The perturbed mode could resonate with its pattern. This could be an appealing analysis too.