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.


1.
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 in heavy ion collisions, such as the ones carried out at the RHIC and LHC 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 QCD at strong coupling. Since this is beyond our reach, alternative theoretical techniques are needed, such as hydrodynamic models [2] or extrapolations from weakly coupled calculations [3]. Predictions are hard to make because every technique has its own limitations.
In this paper, our technique consists of a simple model that is constructed holographically for the N = 4 super Yang-Mills (SYM) theory, using the gauge/gravity duality (in the N c , λ → ∞ limit) [4]. In this context, several authors [5][6][7][8][9][10] have translated the problem into that of a collision of gravitational shock waves in an asymptotically anti-de Sitter (AdS 5 ) spacetime. This 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 [11] These studies have enlightened ostensibly the longitudinal dynamics involved in heavy ion collisions. However, in the actual experiments, the presence of inhomogeneities and the build-up of momentum in the transverse plane could make transverse dynamics important. In fact, a study including radial flow [12] 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.
Indeed, a 3+1 dimensional inclusion of transverse dynamics is non-trivial (but not impossible, see for instance [13] for a hydrodynamical simulation). 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, 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 take longer to thermalize than the background, which indicates that allowing for transverse dynamics in the calculation would lead to a shorter thermalization time.

2.
Gravitational Description. The ansatz for the 5D spacetime metric is a generalization of that in [5], 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 ) coordinates. 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 null geodesics. Also, the ansatz is invariant under arbitrary reparametrizations r → r + ξ, a gauge freedom that will be fixed by placing the Apparent Horizon (AH) in r = 1. The boundary is located at r = ∞.
The ansatz is complemented with h(r, t, y, x 1 ) = h 0 (r, t, y) + e ikx1 δ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 will be linearized around the background solution of [5]. Note that there is no arXiv:1407.5628v1 [hep-th] 21 Jul 2014 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]. But it is important to write them in a fully covariant way using derivatives along outgoing null raysḣ, as defined iṅ 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 as long as the background coefficients satisfy and the inhomogeneities' coefficients satisfy These equations will be used as boundary conditions during the evolution. They can also be derived from the conservation equations of the stress-energy tensor ∇ µ T µν = 0, and 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: After transforming the asymptotic expansions (5) to Fefferman-Graham coordinates, we can use the holographic renormalization prescription to extract the relations where we have omitted the e ik x1 factors in front of every δ term.
3. Numerics Overview. A generic description of the numerical approach that can be applied to solve the dynamics of this problem is found in [14]. 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 (except Σ). 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 F in → S 0 → F 0 → S 0 →Ḃ 0 → A 0 →Ḟ 0 → δS → δF → δG →δ S → δB →δ C →δ D → δA. Note that the dotted functions are solved as if they were unrelated to their undotted counterparts.  Those steps correspond to first the 6 equations for the background, and then the 8 linearized equations for the inhomogeneities. In our approach, we solve for all 14 of these equations at every time slice and afterwards invert (3)(4) to obtain time derivatives and evolve the F in to the next slice.
In addition, there are also 4 constraints: one is of the background and solves for ∂ tṠ , while the other three solve for ∂ t (δ S), ∂ r (δ F ) and ∂ r (δ G). Their asymptotic analysis near the boundary provides the conditions (6-8), which must be imposed as boundary conditions during the evolution. 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 [5] were considered, H(t, y) ≡ µ 3 √ 2πw 2 e −(t∓y) 2 /2w 2 , with w = 0.75/µ. 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).
Several inhomogeneities of the expected stress-energy tensor are plotted in Figs. (1-3). Energy is typically spread out by t ∼ 10/µ for the background (see [5]). However, we find as a common feature that there is still a very uneven profile in these inhomogeneities well after that time, and they take longer to vanish. Also, note that since the equations are linearized, the overall amplitude of these inhomogeneities is completely irrelevant.
We discretized the equations using pseudospectral methods [15] and we chose a background energy density δ = 0.075µ 4 . For many more details, see[1]. 4.
Apparent Horizon. As discussed in [14],the residual gauge freedom r → r + ξ(t, y), is fixed by imposing the AH to lie at a constant r, for instance r = 1.
To find the position of the AH, r AH (y), [5] 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, (11) must be modified by F → F + ∂ y r AH . This is a significant complication, since the problem becomes an intricate non-linear differential equation for r AH (y). Given a time slice t 0 , (11) 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). 5. Discussion. The transverse dynamics in a collision of gravitational shockwaves in AdS has been previously considered in [12], where the longitudinal dynamics were approximated as boost-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 approximation comes from fixing a particular Fourier mode (with momentum k).
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, so no more than pure gravity is needed.
It is well known that at the experiments of the RHIC and LHC, collisions are not exactly central. As a consequence, there are asymmetries in the resulting flow after the collision. In order for the holographic approaches to make contact with this fact, inhomogeneities need to be included.
Our disposal 2 is limited, but also a first step towards a complete calculation.
Based on the behavior manifested by the inhomogeneities (as shown in Figs. (1-3)), we infer that one is to expect a longer thermalization time if inhomogeneous profiles are taken into consideration. Thus, holographic studies may have been underestimating its value.
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. Due to computational limitations, we leave this to future work [16].
An 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 [17,18]. 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).
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.