Inhomogeneous holographic thermalization

The sudden injection of energy in a strongly coupled conformal field theory and its subsequent thermalization can be holographically modeled by a shell falling into anti-de Sitter space and forming a black brane. For a homogeneous shell, Bhattacharyya and Minwalla were able to study this process analytically using a weak field approximation. Motivated by event-by-event fluctuations in heavy ion collisions, we include inhomogeneities in this model, obtaining analytic results in a long wavelength expansion. In the early-time window in which our approximations can be trusted, the resulting evolution matches well with that of a simple free streaming model. Near the end of this time window, we find that the stress tensor approaches that of second-order viscous hydrodynamics. We comment on possible lessons for heavy ion phenomenology.


Introduction
The holographic gauge-gravity duality provides a framework to model the behavior of strongly coupled quantum liquids. Such liquids are studied in many experimental settings: quark-gluon matter created in ultrarelativistic heavy ion collisions, strongly correlated electrons in metals, cuprates and heavy-fermion materials, condensates of ultra-cold atoms. Therefore there is ample motivation for using holographic models to gain insight into strongly coupled dynamics. Holographic models are known to have a very low shear viscosity to entropy density ratio in the strong coupling limit of the boundary quantum field theory [1][2][3], which is also close to the values reported for quark-gluon matter produced at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) [4], for a cold atom gas near the unitarity limit [5][6][7], and for the electronic fluid in graphene [8]. The study of transport coefficients has also evolved to a more direct relation between gravity and fluid dynamics: to the fluid/gravity correspondence [9,10] where one finds that a particular long wavelength limit of Einstein's equations turns them to hydrodynamic equations. One can even derive the complete second order hydrodynamic equations for conformal relativistic fluids [9,11] and generalize the standard Müller-Israel-Stewart theory. More recently attention has been given to equilibration and thermalization from far-from-equilibrium initial conditions. A particularly interesting question is how and when a strongly coupled quantum liquid approaches a regime in which hydrodynamics becomes a good approximation. This paper investigates that question in a particular analytically tractable model, with additional motivation coming from the early dynamics of hot quark-gluon matter in heavy ion collisions.
Relativistic heavy ion collisions have the dual advantage of forming almost perfectly isolated microscopic quantum systems, for which the von Neumann entropy should be exactly conserved, and of being studied in great detail experimentally at RHIC and at the LHC. Experiments at these accelerators have shown that already at very early times, at most 1 fm/c, the matter produced in the heavy ion collisions shows collective behavior in agreement with what is expected from viscous hydrodynamics. The validity of a hydrodynamical description generally relies on the fluid being near local thermal equilibrium, but some degree of collectivity may exist even before local thermal equilibrium is reached. This conjecture is supported by recent experiments studying p+Pb collisions at LHC [12,13] and d+Au collisions at RHIC [14], which show indications of the presence of hydrodynamic behavior in events producing a large particle multiplicity. On the other hand, several theoretical approaches to the dynamics of thermalization in QCD, e.g. the perturbative bottom-up thermalization framework [15] and classical gauge theory [16], predict substantially longer equilibration times. Heavy ion experiments thus offer a playground for the study of equilibration and thermalization in gauge theories.
The argument about the presence of hydrodynamic behavior of the quark-gluon plasma created in heavy ion collisions rests primarily on two observations. First, one observes a cos(2φ) correlation between the azimuthal momentum direction of produced hadrons and the collision plane, which is known as "elliptic" flow, see [17][18][19][20] and references therein. This phenomenon can also be deduced from the azimuthal two-particle correlations among emitted hadrons, and is sometimes referred to as the "ridge" or "double-ridge" effect in studies in Pb+Pb collisions at LHC [20][21][22][23]. The second observation is related to eventby-event fluctuations. Experimentally it was found for symmetric heavy ion collisions that odd Fourier coefficients of the flow are not much smaller than even ones [20,23,24]. By a symmetry argument, odd coefficients can only be generated by fluctuations, so one is forced to conclude that fluctuations are large. The parton saturation model for the initial nuclear state suggests also that they are of short range in the plane transverse to the beam axis (of order of the inverse saturation scale 1/Q s ). Detailed simulations [25][26][27][28][29] show that if one assumes also these small size fluctuations to evolve hydrodynamically one gets excellent fits to all flow coefficents.
The holographic description of strongly coupled gauge theories offers a framework for the study of the rate of thermalization and the approach to hydrodynamical flow. In [30][31][32][33], it was found that in a simple model thermalization after the homogeneous, isotropic injection of energy proceeds very fast. The initial state of heavy ion collisions is, however, quite anisotropic and inhomogeneous. It is characterized by a strong asymmetry between longitudinal and transverse pressure (the former may even be negative initially due to the presence of strong longitudinal gauge fields) and by large density fluctuations in the transverse direction. The effect of both can be studied with a refined AdS/CFT treatment. The effect of the pressure anisotropy was studied in [34][35][36] under the assumption of longitudinal boost invariance and transverse homogeneity; it was found that hydrodynamic behavior is reached on time scales of order 0.3 − 0.5 fm/c for many different initial conditions. This "hydroization" is not equivalent to complete thermalization because viscous hydrodynamic behavior at early times in a boost invariant expansion implies a rather large pressure anisotropy and thus strong deviation from local thermal equilibrium.
In the present paper we analyze the second aspect, namely the question how the approach to hydrodynamic behavior is affected by local density fluctuations, which has not been studied in detail so far. A priori, there is no reason to believe that hydrodynamics, which is a low-energy effective description, should describe the short-time evolution we will be studying. However, since hydrodynamics turned out to apply unexpectedly early in homogeneous models, it is interesting to ask to what extent this surprise extends to inhomogeneous setups.
Thermalization in strongly coupled conformal field theories with a gravity dual corresponds to black brane formation in asymptotically anti-de Sitter (AdS) spacetimes. To study gravitational collapse, one generically needs numerical general relativity, but interesting situations exist where analytic computations are possible. Consider a massless scalar minimally coupled to gravity in d + 1 spacetime dimensions with negative cosmological constant. In this setting, the authors of [37] considered the effect of a homogeneous boundary source on the bulk geometry of an asymptotically AdS spacetime. Specifically, they turned on for a finite lapse of time δt a homogeneous source φ 0 (t) for a marginal boundary operator corresponding to the massless scalar field in the bulk. Solving the field equations in a weak field approximation, i.e. in the limit where the amplitude of the scalar source is small in an appropriate sense 1 , in [37] it was shown that this generically leads to black brane formation in the bulk (see also [38][39][40] for numerical analyses). For d + 1 = 4 and to leading order in the amplitude of the scalar source, the geometry is given by the AdS 4 -Vaidya metric with M (v) a non-decreasing function that vanishes for v ≤ 0. This geometry describes a shell of null dust falling in from the boundary of AdS and collapsing into a black brane. A schematic representation of the process is depicted in Fig. 1.
The AdS-Vaidya model has been very useful as a simple, tractable model of holographic thermalization. Many observables in the dual field theory have been identified whose timeevolution can be computed relatively easily from geometric quantities in AdS-Vaidya [30][31][32][33][41][42][43][44][45][46][47][48][49][50]. When this model is used as a very crude analogue of the equilibrating quark-gluon plasma in heavy ion collisions, an encouraging result is that, at least for the range of length Figure 1: A schematic representation of the dynamical collapse process first studied in [37]. In pure AdS spacetime, a scalar source φ 0 on the boundary, corresponding to the asymptotic boundary value of a massless bulk scalar field, is turned on for a lapse of time δt. This induces a non-vanishing profile for the bulk scalar field which backreacts on the AdS geometry as a shock wave propagating in the bulk leading to the formation of a black hole horizon. scales studied in [32,33], thermalization for homogeneous energy density occurs as fast as allowed by causality. 2 This insight suggests that hydrodynamics can already describe early stages of heavy ion collisions with strong coupling dynamics. Massive infalling shells were studied in [51][52][53][54][55][56][57][58][59], charged shells in [60][61][62][63] and shells with angular momentum in [64].
In this paper we generalize the construction of [37] to the case of an inhomogeneous scalar field source at the boundary. As in [37], we solve the equations of motion in a perturbative expansion in the amplitude of the scalar field boundary value, which is here allowed to depend on the spatial coordinates. Before the source is turned on, the solution coincides with the pure AdS background. When introducing a dependence on the transverse spatial coordinates, the situation becomes much more involved and finding an exact analytic solution is complicated. In order to still get an analytic insight, in this paper we shall therefore attack the problem under the simplifying assumption that the scale of spatial variations is large as compared to all other scales, following a strategy suggested in [37].
As in [37], the case of a four-dimensional bulk space-time turns out to be technically simpler than the five-dimensional one. 3 Ideally, one would like to introduce space dependence for at least two directions, corresponding to the plane in which the pancaked nuclei overlap in the collision, but for simplicity we consider the simplest case involving spatial dependence on a single coordinate (in addition to the radial coordinate in the bulk). So we consider an asymptotically AdS 4 geometry with inhomogeneities along a single spatial direction.
In section 2, we construct the bulk solution up to second order in the amplitude of the source that drives the gravitational collapse and up to fourth order in the gradient expansion. We argue that our solutions should be reliable for times short compared to the inverse local temperature (of the black brane to be formed) and wavelengths large compared to the inverse local temperature. In section 3, we extract from the bulk solution the expectation value of the boundary stress-energy tensor, and study its time-evolution after the inhomogeneous energy injection. We compare this evolution to that of a simple freestreaming model as well as to first and second order hydrodynamics in section 4. Section 5 contains a summary and a discussion of the possible relevance for heavy ion phenomenology.
A short account of our main results can be found in the companion paper [103].

AdS 4 weak field inhomogeneous collapse in the gradient expansion
We consider a massless scalar minimally coupled to gravity in four spacetime dimensions with negative cosmological constant, Here Λ = −3 and the AdS radius has been set to one. The equations of motion following from (2.1) read When assuming inhomogeneities along only one of the two boundary spatial directions, say x, the ansatz in Eddington-Finkelstein coordinates consistent with the symmetries of the problem can be written as where v is an ingoing null coordinate, r the AdS radial coordinate and y denotes the other spatial direction.
On the boundary, where the null coordinate v is identified with the gauge theory time t, we turn on an inhomogeneous scalar field source over a time δt. Here we have singled out the explicit factor parametrizing the amplitude of the scalar source for the bulk field, which with the non-canonical normalization (2.1) of the bulk kinetic term is dimensionless. As in [37], we require the forcing function ϕ 0 (v, x) and its first few time derivatives to be everywhere continuous.
The initial condition that the spacetime should be pure AdS for v ≤ 0 corresponds to while the asymptotically AdS boundary conditions (with planar boundary geometry) read and lim The metric (2.3) supplemented with the boundary conditions (2.7) is not completely gauge fixed. The form of the metric is left unchanged under the transformation r → r + s(v, x). We use this residual freedom to choose the subleading behavior of f (v, r, x) to be f (v, r, x) = r(1 + O(1/r 2 )).
With an ansatz of the form (2.3) one has, in addition to the scalar equation of motion, seven equations coming from (2.2). In fact, of the ten components of E µν , three turn out to identically vanish and one of the remaining is a linear combination of the others together with the scalar equation. The following linear combinations turn out to be a convenient choice to work with: where E ci = 0, i = 1, 2 are conservation equations, and one of the equations is implied by the others.
The approach we adopt here is to solve the system (2.9) in a double expansion. Following [37], we perform a weak field expansion in the amplitude of the scalar field source. On top of this, we consider a spatial gradient expansion along the non-homogeneous space direction, i.e. an expansion in spatial derivatives 4 . The metric components and scalar field are then written in a double expansion in the parameters , which keeps track of the order in the amplitude of the source, and µ, which acts as a formal derivative counting parameter, as (2.10) Accordingly, the boundary source is written as (2.11) and we assume that the function ϕ 0 (v, µx) and its derivatives are everywhere continuous. The derivative counting parameter µ shall be set to one in the solution at the end of the computations. For compactness in what follows we refer to the order in the number of derivatives as the order in µ.

General structure of the equations and of the solution
The different coefficients in the expansion (2.10) are determined order by order as a function of the previous orders solution. Implementing the pure AdS initial conditions (2.5) and (2.6), as well as the asymptotically AdS conditions (2.7) and (2.8), this implies that these coefficients involve only the forcing function ϕ for the bulk scalar profile.
There are a few general considerations about the structure of the two expansions (2.10) worth making before explicitly solving the equations. The background solution of the derivative expansion (O(µ 0 )) coincides with the homogeneous solution of [37], with the xdependence added by hand to the source ϕ. At zeroth order in µ and first order in , φ 1,0 is forced by its boundary condition to be non-zero. The background geometry is pure AdS at this order, since there is no O( ) source for the Einstein's tensor coming from the bulk scalar. Only at order O( 2 ) are the metric components h 2,0 and f 2,0 sourced.
When we add the x-dependence to the homogeneous solution, it ceases to be a solution to the equation of motion. At leading order in the equation for φ is just the massless scalar equation in pure AdS 4 . By symmetry of the background, the terms violating the homogeneous scalar equation are of order µ 2 , i.e. second x-derivatives of the source. The background bulk scalar solution thus needs to be corrected with a term of order µ 2 . A similar reasoning applies recursively at higher orders and for the metric which starts at order 2 . Notice that the scalar profile at order receives only even order corrections in µ. Similarly all metric coefficients at order 2 receive only even order corrections in µ, with the exception of k(v, r, x), which has only odd µ contributions.
Solving the equations (2.9), one recognizes a structure similar to that observed in [37]. In the amplitude and gradient expansion, the scalar equation of motion is solved using a 1/r expansion and this automatically leads to metric components expressed as an expansion in 1/r. As illustrated in Appendix A, the initial conditions are satisfied by imposing that the 1/r expansions should terminate at a finite order.
Order by order in and µ, E 1 = 0 gives a differential equation in r for the metric component f n,i , which is determined completely by imposing the asymptotic boundary conditions. E 2 = 0 similarly gives a differential equation in r for h n,i , which is completely fixed imposing the boundary conditions together with the conservation equation E c1 = 0. In the same way, the coefficients k n,i are determined using E 3 = 0 and E c2 = 0. The former gives a differential equation in r for k n,i while the latter is a conservation equation, which together with the boundary conditions completely fixes the form of k n,i . Using the rest of the equations, E 4 = 0 and E 5 = 0 follow from each other. The resulting independent equation is the only equation involving b n,i and suffices to determine its expression.

Solution at first and second order in the source amplitude
Solving the system (2.9) following the strategy outlined in the previous Section, we obtain the following solution up to second order in the amplitude of the source field and up to fourth order in the spatial derivative expansion. The compact form of the solution with the different orders in the and µ expansions grouped according to their radial asymptotic behaviour is given in Appendix A for convenience. Here we present the solution as obtained order by order in the two expansion parameters.
In the following we adopt the notation

First order in
At zeroth order in µ (we set µ = 1 in writing the solution) which, when removing the x-dependence, coincides with the homogeneous solution obtained in [37]. The next non-vanishing contribution is at O (µ 2 ) while at fourth order in µ one obtains the additional contribution (2.14) All the metric coefficients vanish at this order in .

Second order in
The scalar profile in the bulk does not get contributions at this order, while the metric at lowest (zeroth) order in µ is the one of the homogeneous case, with the spatial dependence added. Explicitly, the solution of (2.9) reads and all other coefficients vanish. The coefficient of the cross component dvdx of the metric gets the first non-vanishing contribution at order µ where (2.18) At order µ 2 all other coefficients get sourced in Einstein's equations. Their expressions are (2.20) Going further with the expansion in spatial derivatives, a third order contribution to the coefficient for the cross component dvdx of the metric comes into play Finally, at fourth order in the gradient expansion

(2.24)
A simple consistency check on the different terms entering in the components of the metric and in the bulk scalar comes from their scaling, with derivatives contributing with weight 1 and factors of 1/r and integrations contributing with weight −1.
All solutions reduce exactly to the homogeneous formulae [37] when the spatial dependence of the scalar source is suppressed.

Regime of validity
Before continuing and analyzing the evolution of the stress-energy tensor of the dual field theory, let us pause to comment on the features of the solution, the details of the expansions and on their regime of validity.
Let us start from the metric itself. As stated above, it describes the dynamical process of black hole formation due to a wave of energy triggered by the forcing function at the boundary, which propagates into the bulk. For simplicity, we shall first refer to the homogeneous limit where ϕ(v, x) = φ 0 (v). The metric in this case takes the form (cf. h 2,0 (v, r, x) and f 2,0 (v, r, x) in (2.15)) Away from the energy injection period 0 < v < δt, the metric (2.25) has the AdS-Vaidya form The spacetime is AdS 4 for v ≤ 0 (M (v) = 0); for v ≥ δt, on the other hand, M (v) is a constant M determined by the integral in (2.25) and the geometry approximates a black hole geometry with temperature T ∼ M 1/3 ∼ 2/3 δt . (See [37] for a more detailed analysis including a discussion of the horizon formation.) By analogy, in our case, we identify the coefficient of the 1/r term in h(v, r, x) as a notion of local "temperature" (by which we do not mean to suggest local equilibrium has yet been established). Once the forcing function has been turned off, we have with the explicit form of C 2,i given above in (2.16), (2.20) and (2.24).
We have obtained the solution in a weak field approximation for the boundary source, i.e. assuming the amplitude of ϕ is small (as indicated by the formal expansion parameter ). An extensive discussion of the regime of validity of this approximation has been presented in [37]. The conclusion is that the amplitude expansion is perturbative in tT , 5 and it is therefore reliable as long as t 1/T . In order to capture the evolution for longer times, a resummed perturbation theory is needed, where one expands around AdS-Vaidya rather than AdS. This amounts to working exactly in T and perturbatively in all other appearances of , which is similar to absorbing temperature-dependent masses in propagators in thermal perturbation theory. At late times, observables approach their thermal values exponentially in tT . In naive (non-resummed) perturbation theory in T ∼ 2/3 /δt, the exponential series are truncated to finite order, leading to polynomial expressions that diverge at late times. Notice, however, that in the homogeneous case this effect is not yet present up to second order in the expansion, in the sense that there are no divergent terms at this order. At leading order in , the scalar field vanishes at late times, while the metric has the black hole form (2.26) following from (2.25).
The second expansion we introduced is the derivative expansion along the inhomogeneous direction. At finite temperature this corresponds to asking that the scalar source should be slowly varying in x over distances set by the local inverse temperature [37]. After the injection of energy has taken place, if λ is the typical length scale over which the source ϕ(v, x) varies, the derivative expansion holds for λ 1/T . However, here the situation is more involved. For v > δt, (2.27) incorporates the general structure of the solution, which very schematically we can write up to second order in the derivative expansion as (see (3.26) for the derivation) The correction proportional to B(x) can be considered small compared to the leading order term as long as λ δt. Since 1/T ∼ δt/ 2/3 , the condition λ δt is automatically implied by λ 1/T for 1. Moreover, there is a perturbative structure in v/λ that is similar to the effective perturbative structure in tT of the amplitude expansion. 6 In order to have small corrections coming from the gradient expansion, it is therefore also necessary that v λ, which is also automatically satisfied since 1/T λ and v 1/T . There is an additional subtlety regarding the period 0 < v < δt during which the energy is injected. This is apparent in the analysis of the scalar field solution at order (see (2.12)-(2.14)) The first correction involving spatial derivatives is present only as long as the scalar forcing function is turned on. At this order the background is pure AdS, for which the local temperature is effectively zero and the condition λ 1/T is not satisfied. However, as long as the correction obtained at order λ −2 is small compared to the leading λ 0 term, one can reliably work within the derivative expansion. This translates into the condition 1/r λ; at smaller values of the radius we lose control over the perturbative description. For large λ, the part of the v ≈ 0 region that is uncontrolled is deep into the bulk and can only affect the boundary at late times, so that close to the boundary the perturbative expansion is reliable except at late times.
All in all, we can conclude that our expansions are reliable as long as we ensure that t 1/T λ.
3 Evolution of the boundary stress-energy tensor

Boundary stress-energy tensor
The relation between an asymptotically AdS solution for the action (2.1) and the expectation values of the boundary stress-energy tensor and of the boundary operator associated to the massless bulk scalar has been obtained in Appendix B. For a solution written in Fefferman-Graham coordinates the relation is summarized in equations (B.12) and (B.13) where α and β are indices over the boundary coordinates. For notational simplicity we henceforth omit the expectation value symbol.
To read out the stress-energy tensor corresponding to our solution we thus need to perform a change of coordinates to reach the Fefferman-Graham form. The ansatz we started with has the form 3) and the metric we obtained up to order 2 and µ 4 is schematically of the form while the scalar field solution can be written as By comparison with the explicit form of the solution obtained in the previous Section one can identify the expressions of the various coefficients (see also Appendix C). The relevant ones here are and where we have suppressed the explicit coordinate dependence. Any asymptotically AdS solution can be brought into Fefferman-Graham form close enough to the boundary (r → ∞). We therefore look for a change of coordinates of the form (v, r, x, y) → (t, , χ, y), such that the metric in the new coordinates takes the form where x α = (t, χ, y). Any dependence on y, both in the change of coordinates itself and in the metric, has been excluded a priori for symmetry reasons. Notice that the coordinate is simply related to the coordinate used in Appendix B by = 1/z. We find it more convenient to work with the coordinate here, so that the boundary r → ∞ corresponds to → ∞. In these coordinates, the scalar field has the Fefferman-Graham expansion Let us remark that we do not look for an exact change of coordinates, but we work in a large radius expansion at a sufficiently high order to determine the boundary field theory stress-energy tensor.
For v ≤ 0 the spacetime is pure AdS; the change of coordinates is exact and it reduces to the standard coordinate transformation relating Eddington-Finkelstein and Poincaré coordinates. In practice we shall look for a change of coordinates such that g (0),αβ = η αβ = diag(−1, 1, 1). Working perturbatively in the radial variable one constructs the transformation that brings the scalar field and the metric in the required form. The details are worked out in Appendix C, the net result is and For the scalar field one obtains in the large expansion It has exactly the structure expected from the analysis in Appendix B: there is no 1/ term For what concerns the metric g αβ : g (0),αβ = η αβ = diag(−1, 1, 1), while g (2),αβ can be written as which is the condition obtained in (B.8). After some algebra (see Appendix C) is also satisfied, as one can check by explicit computation. From this analysis the stress-energy tensor and the operator expectation values (3.1) and (3.2) can be readily obtained. Renaming the boundary coordinates in a natural way we have 17) The explicit expression in terms of the scalar source can be obtained from (3.6). For the operator dual to the bulk scalar instead one has . (3.18) In the homogeneous limit (ϕ(t, x) = φ 0 (t)) (3.17) reduces to which is the stress-energy tensor of a perfect conformal fluid in three spacetime dimensions.

Evolution following the energy injection
The explicit expression of the boundary stress-energy tensor in terms of the forcing function ϕ(t, x) can be directly read from (3.17), although it is quite involved. Some of the main features have been already outlined in the previous Section. We refer to Appendix D for the complete expressions of the stress-energy tensor up to second order in the scalar source and up to fourth order in the gradient expansion. As an example, we give here the explicit form of T tt , at second order in the derivative expansion To make progress in the analysis and to simplify the subsequent numerical treatment, we shall now assume that the temporal dependence of the source function is factorized with respect to the space profile and write Substituting this ansatz in (3.20) and performing a number of integrations by parts, the corresponding energy density profile takes the form (3.23) For a source profile that is well localized in time between 0 and δt as in (3.22), the time dependent functions in (3.23) have the following simple behaviour: Correspondingly, after all energy has been injected into the system, T tt has a simple polynomial time dependence (3.26) As already pointed out in Section 2.3, this structure comes from the nested integrals over time and gives an effective expansion in t times the inverse length scale of spatial variations. The structure, which at first seems to lead to a divergent profile at arbitrary late time, has to do with our approximations, which fail for large times. Notice there is no contradiction with energy conservation. In fact time dependent terms in (3.26) involve spatial derivatives and their contribution to the total energy is vanishing. (The same is true at fourth order in the gradient expansion, as one can check from the explicit expression of T tt in (D.1).
After the injection of energy has completed (t ≥ δt), time dependent contributions to T tt correspond to terms involving nested time integrals in (D.1), all of which are total spatial derivatives.) In order to study the stress-energy tensor in more detail and to better understand its evolution, we shall use a Gaussian profile to mimic the time dependence of the source (3.27) Such a replacement retains the main featurlarger thanes of the evolution discussed above: when t − ν is large enough, the scalar source can be considered to vanish for all practical purposes.
Consider first the simple case where the spatial profile u(x) has the form of a Gaussian on top of a homogeneous background (3.28) Here we insert the explicit factor µ ∼ 1/λ in the space dependence of the source in order to implement more easily the slowly varying regime. We have considered a homogeneous component in the source field, which gives a non-zero energy density everywhere. In this way the local "temperature" is everywhere non-vanishing and we can keep under good computational control both the weak field and the gradient expansion.
In Figure 2 (top-left panel) we report a sample plot of the T tt component of the stressenergy tensor for ν = 0.5, σ 2 = 0.1, µ = 0.01 and = 0.005. In order to stop the evolution at a time which is reliably within the range of validity of our approximations, we examine the effects of sub-leading contributions in the gradient expansion. For a fixed spatial interval, the expansion is valid till times at which the sub-leading contributions become comparable to the leading order result. In line with our general discussion of the regime of validity of the amplitude expansion, we also impose that t 1/T , with T the local "temperature" defined in the previous Section by analogy to the equilibrium case. From the point of view of T tt , the evolution seems to point towards homogenization, although we cannot follow the evolution for arbitrarily large times. The same qualitative conclusion seems to hold looking at the evolution of the other diagonal components. In the top-right and bottom-left panels of Fig. 2, the T xx and T yy components of the stress-energy tensor are plotted for the same values of the parameters, while the only non-vanishing offdiagonal component T tx is plotted in the bottom-right panel. As opposed to the diagonal components, T tx grows in time, departing from a homogeneous configuration. However, the off-diagonal component of the boundary stress-energy tensor also gets a contribution from the local boost velocity, and the corresponding profile seems to be qualitatively compatible with this interpretation: Intuitively, a lump of energy with a density profile as in Fig. 2 would start spreading outward and the local velocity with respect to the background is maximal where the gradient of the energy density is large.
To have a better understanding of the situation it is convenient to perform the boost that locally brings the stress-energy tensor in the diagonal form T αβ = diag(ε, p x , p y ). This boost in the x direction with rapidity α relates the components in the two frames as (3.29) The local velocity of the plasma in the stationary frame is (3.30) Using (3.29) and (3.30), we can plot in Fig. 3 the energy density and pressures of the fluid in the local rest frame, as well as its local velocity V . The energy density and pressures tend to flatten out as expected, while the pressures anisotropies still build up in this phase as shown in Fig. 4. In order to ensure that the gradient expansion is reliable, in Fig. 4 we restrict the time range such that third (fourth) order terms in the gradient expansion are at least two orders of magnitude smaller than first (second) order terms. In Figures 2 and 3 we relax the time range such that the corrections coming from the subleading order in the gradient expansion should be one order of magnitude smaller than the leading order result, because otherwise the time evolution of the various quantities would not be visible in the plots. Similar evolution was found for instance by Chesler and Yaffe in [71], where it was observed that the generation of anisotropies in the stress-energy tensor due to an anisotropic source continues even after the end of energy injection. The regime where the anisotropy is built up should then be followed at later times by an evolution towards isotropy. As expected from the discussion in Section 2.3, this relaxation dynamics is not visible in naive perturbation theory and should start over timescales set by the inverse temperature. Also, the local boost velocity is qualitatively consistent with that of a lump of fluid that spreads out. It decreases for large values of |x|, consistently with expectations for a lump of fluid with higher density than the surrounding homogeneous medium.

Comparison with free streaming and hydrodynamics
In order to gain some insight into these results, we next compare them with the evolution of the stress-energy tensor in a free-streaming model as well as in hydrodynamics.

Anisotropic free streaming
To compare the evolution we find using AdS/CFT with a simple model of free-streaming particles, we construct a kinetic model of the injected energy distribution. Our model is obtained by assuming that the distribution is massless noninteracting dust, composed of particles moving at the speed of light. In terms of the phase space distribution f (t, x, k), the components of the stress-energy tensor are given by If we assume the phase space distribution at some time t = δt to have the factorized form f ( x, k) = n( x)F ( k) = n(x)F (k), then at a later time where v = k/k t is the particle velocity, and F (k) only appears in (4.1) in an overall normalization factor. It turns out that the specific form of the function F (k) is irrelevant,  Figure 4: The spatial profile of the pressures anisotropies p x − p y (rescaled by 16πG N ) for the pressures plotted in Fig. 3. The plot is limited to shorter times as compared to those considered in Fig. 3 to ensure that the corrections to p x and p y coming from higher orders in the derivative expansion are negligible compared to the pressure anisotropies at leading order. In the plot time increases from the top curve down.
what matters is that all particles move at the speed of light. When comparing with our AdS/CFT results, we start the free streaming evolution just after the energy injection has ended, namely at t = δt, and choose n(x) and F (k) such that the energy density in the laboratory frame coincides with that computed using AdS/CFT. We can introduce polar coordinates and write k x = k cos φ, k y = k sin φ, v x = cos φ, v y = sin φ. One easily finds that T xy = 0, which means that for this particular choice there is no effective flow developing. This does not have to be the case when the k and x dependence are not factorized. The non-vanishing components of T αβ are: where the common factor is given by the initial average energy per particle: The tracelessness of T αβ is easily seen from these expressions.
Notice that at t = δt, the stress-energy tensor in (4.4) is isotropic, while in (3.17) anisotropies start to develop already during the energy injection. However, their relative amplitude only starts to build up considerably after a time δt, and therefore the difference between the two stress-energy tensors at t = δt is not significant.
In Fig. 5, we plot the rest frame energy density ε, the pressures p x and p y and the local velocity V of the plasma during a free streaming evolution. The curves closely resemble The energy density in the rest frame ε fs = p x,fs + p y,fs , pressures p x,fs and p y,fs (rescaled with 16πG N ) and local boost velocity V fs within free streaming. The match of the free streaming and the AdS/CFT energy densities in the laboratory frame has been performed at t = 0.009/µ, when about 96% of the energy has been injected. The evolution closely resembles that obtained with the AdS/CFT correspondence and plotted in Fig. 3.
those obtained from the holographic computation in Fig. 3, although the maximal amplitude of the energy density and pressures decreases more slowly under anisotropic free streaming. Nevertheless, the pressure anisotropy in the local rest frame matches well that obtained within AdS/CFT, as shown in Fig. 6.

Viscous hydrodynamics
Hydrodynamics is a low-energy description that is not meant to be valid on time scales smaller than the inverse local temperature. Nevertheless, viscous hydrodynamics turns out to give surprisingly good results in several homogeneous holographic thermalization setups [34][35][36]. Motivated by this, we now compare our AdS/CFT results with what one would obtain by naively applying the formulas of viscous hydrodynamics. In fact, in [37] it has been suggested that the anisotropies present after the injection of energy could be interpreted in terms of first order viscous hydrodynamics, for which the stress-energy tensor is completely expressed in terms of the local temperature and of the local velocity of the fluid (a review of the basic ingredients can be found for instance in [10]). Given the local velocity V (t, x) and the energy density in the local rest frame ε(t, x), we can ask what the stress-tensor would be if hydrodynamics were valid at a given time t. Agreement with the stress tensor we computed using AdS/CFT would be a necessary condition for the validity of a hydrodynamical description from time t onwards.

First order viscous hydrodynamics
First order hydrodynamics requires the equation of state and the shear and bulk viscosities, which can be found, for instance, in [10,106] for the three-dimensional CFT under consideration. The first-order hydrodynamical stress tensor reads Here p ideal = ε/2 and u α is the local three-velocity of the fluid (which we determine from V (t, x)). The first order viscous contributions to the stress-energy tensor in flat threedimensional space are given by [10] Π αβ where P αβ ≡ u α u β + η αβ is the projector onto space in the local fluid rest frame, is the fluid shear tensor and θ ≡ ∂ ρ u ρ is the expansion. The shear viscosity η and bulk viscosity ζ are [10,106] where we define a local "temperature" T (t, x) by using the thermodynamic relation that would be valid in equilibrium in the local rest frame [106], (4.10) To compare with the stress tensor we have worked out, we use the fluid velocity V = − tanh α determined in the previous Section for a single Gaussian scalar profile.
Determining the pressure components p x,hydro and p y,hydro in the local rest frame, the first order pressure anisotropy reads p x,hydro − p y,hydro = −2ηθ. (4.11) In Fig. 7, this is compared with the pressure anisotropies computed in the gradient expansion. The two have amplitudes of the same order of magnitude and similar shape, although they do not coincide. In particular the ones computed in the gradient expansion are everywhere smaller in amplitude than the hydrodynamics ones.

Second order hydrodynamic corrections
Since first order viscous hydrodynamics did not give a sufficiently accurate description of the time-evolution of the initial inhomogeneities, we move to second order formalism.
Second order hydrodynamics for relativistic conformal fluids was derived in [9,11], providing a nonlinear generalization of Müller-Israel-Stewart theory, used in modeling heavy ion collisions. The relevant second order contribution to the dissipative part of the stress tensor is (see eqn (3.11) of [11]) where " + · · · " refers to terms that are not relevant here, as they either contain the Riemann/Ricci tensor or terms that vanish for irrotational flow. In the above, D is the directional derivative along the 3-velocity: while denotes the transverse traceless part. The new parameter τ Π is a relaxation timescale, which for a strongly coupled conformal fluid depends on harmonic numbers [106] and in 2+1 dimensions (4.14) The central feature of the second order theory is that now hydrodynamics becomes causal, with the discontinuities/inhomogeneities propagating at finite velocity [9,11] v disc = η τ Π ( + p ideal ) = 1 3 + Harmonic − 1 3 ≈ 0.665 . For the second order term in (4.12) we need to evaluate first and then construct its transverse traceless projection Dσ αβ . The first term in (4.16) is already transverse and traceless. By an explicit calculation one can verify that the second term projects to a null matrix. The second order contribution to the dissipative part of the stress tensor thus becomes (4.18) and in total then Explicitly in our case the pressure anisotropies in the local rest frame read: (4.20) The first and second order viscous contributions to the pressure anisotropies have opposite signs and the latter dominate over the first at early times. The dominant term at early times is in fact 2ητ Π u t ∂ t θ.
As we saw in Figure 7 in the previous section, in first order hydrodynamics the pressure anisotropies p x −p y were qualitatively similar to those found from the AdS/CFT framework, but their magnitude grew faster and thus there was never good quantitative agreement. Including the second order correction induces a term that tends to drive the anisotropy to the opposite direction. The net effect is that at very early times there is not even qualitative agreement, but shortly before t ∼ 0.1/µ (which is the estimated upper time limit on the validity of the gradient expansion) the evolution of inhomogeneities reaches a stage where in addition to the qualitative agreement (spatial profiles again having similar shapes) there is even good quantitative agreement, as can be seen in Fig. 8. If we were to extend the evolution to later times (say up to t ∼ 0.15/µ), then the match with second order viscous hydrodynamics appears to become somewhat less good in the region x 100, while continuing to improve in the region x 100. On the other hand, free streaming appears to remain a very good approximation. Given the very short time window our methods allow us to describe, more work is required to confirm whether the early-time agreement with second order hydrodynamics that we find is accidental or really signals the onset of a hydrodynamic regime.
Our analysis is performed in the infinite coupling limit, so one may ask how finite coupling corrections are expected to affect the results. In particular, for some observables such as the spectral density in photon production in a holographic model of non-equilibrium plasma, it has been reported that finite coupling corrections become significant already at large values of 't Hooft coupling λ [55,57]. Other recent work on finite coupling corrections in the context of thermalization is reported in [96,[107][108][109]. However, in the present case the relevant parameters to consider are the shear viscosity η and the relaxation time τ Π . Finite coupling corrections to the shear viscosity (more precisely, to its ratio with the entropy density) have been computed in [110][111][112][113] and corrections to the relaxation time have been computed in [114]. In both cases, the finite copling corrections become sizable only at small values of 't Hooft coupling, λ ∼ O(1 − 10). Based on that, one may expect our results to be applicable at strong but finite coupling.

Further examples
We have also performed tests with a modified source. Consider first a spatial profile u(x) given by the superposition of two Gaussians with separation d between the center of the two profiles (4.21) In Fig. 9 we have plotted the evolution of the spatial profile of the energy density ε obtained from the AdS/CFT computation for d = 2/µ (left) and d = 6/µ (right). The bottom of Fig. 9 gives a comparison between the evolution of the spatial profile of the pressure anisotropies obtained from AdS/CFT (solid lines), second order hydrodynamics (dotted lines) and free streaming (dashed lines), as in Fig. 8. It confirms the picture emerging from the single spatial Gaussian analysis. This conclusion naturally extends to more general superpositions of spatial Gaussians. Let us for instance consider a scalar source with four Gaussian peaks of different size at different locations, on top of a homogeneous background: The corresponding energy density profile obtained within the AdS/CFT approach is plotted in the left panel of Fig. 10. In the right panel the pressure anisotropies obtained in the AdS/CFT computation are compared with those of free streaming and of second order viscous hydrodynamics.

Summary and conclusions
We have generalized the "naive" (non-resummed) weak-field perturbation theory part of [37] to include inhomogeneities in a long-wavelength expansion. This has allowed us to study  Figure 9: Top: the energy density in the rest frame ε = p x + p y (and rescaled with 16πG N ). Bottom: the spatial profile of the pressures anisotropies p x − p y (rescaled by 16πG N ) obtained in the AdS/CFT computation (solid lines) compared with those resulting from second order hydrodynamics (dotted) and free streaming (dashed). The source profile has the form (4.21), with ν = 0.5, σ 2 = 0.1, µ = 0.01 and = 0.005. The left column corresponds to d = 2/µ, the right one to d = 6/µ. the evolution of the stress tensor of the dual field theory for times short compared to the local inverse "temperature", which in turn should be small compared to the local scale of spatial variation. For specific temporal and spatial profiles of the source that injects energy, we have compared the evolution to that of a simple free-streaming model, and have found very good agreement. In general, hydrodynamics is not meant to describe evolution on (time)scales short compared to the local inverse temperature. However, in various homogeneous models it was found that the stress tensor agreed with that of viscous hydrodynamics much earlier than one had a right to expect [34][35][36]. In our inhomogeneous model, for the short time window we are able to explore analytically, we find poor quantitative agreement with first order viscous hydrodynamics. Remarkably, though, second order hydrodynamics does provide very good agreement near the end of the time window in which we trust our computations. This could be a sign of a rapid cross-over from a free-streaming regime to second order hydrodynamics. However, more work is needed to verify whether this is really the beginning of a hydrodynamic regime, rather than an accidental agreement that fails to hold at later times than those we can probe using our methods. To settle this, one would need either an inhomogeneous extension of the resummed perturbation theory of [37], or a numerical analysis.
In our analysis, we have made several assumptions for technical simplicity: the injected energy only depends on one spatial coordinate, and the field theory lives in a 3d spacetime. In particular, our analysis does not capture the longitudinal expansion of the quark-gluon plasma. Our results are therefore complementary to those of [34][35][36], where the approach to hydrodynamics is studied for a spatially homogeneous, but longitudinally expanding fluid. Taken together, and with the important caveat mentioned in the previous paragraph, these two results suggest that free streaming followed by second order viscous hydrodynamics may provide a valid description of the space-time evolution of the stress energy tensor in strongly coupled conformal gauge theories over the entire history.
The reader may wonder why free streaming can reproduce the behavior of a strongly coupled gauge theory for any length of time. We do not have a complete answer to this question, but it is tempting to speculate that the short-time behavior of the stress-energy tensor is dominated by the singular behavior of the two-point correlation function of components of the stress-energy tensor. In a conformal theory, the correlation functions have a power-law dependence on the invariant space-time distance and diverge on the light cone. This suggests that free streaming of massless "particles" may reproduce the short-time behavior of these correlation functions [115].
We end this paper with a discussion of how our results may help provide a justification for a standard approach to early time evolution of partonic matter in high-energy heavy ion collisions as they are investigated at the CERN Large Hadron Collider (LHC) and the Relativistic Heavy Ion Collider (RHIC).
Due to the large Lorentz γ factors of these collisions in the center-of-mass frame, the quark and gluon structure of the colliding nuclei remains essentially frozen during the collision, while the coherence of the nuclear quantum states is broken by soft color exchanges. A large effort was invested into the theoretical investigation of this initial state within a framework called the Color Glass Condensate (CGC) Model. Its central assumption is that nonlinear processes limit the growth of the gluon density at high energy, resulting in its eventual saturation. The model posits that the gluon distribution in the transverse direction attains a universal form -the color glass condensate -that is characterized by a single parameter, the saturation scale Q s , which only depends on the size of the nucleus and the beam energy. The very large transverse area density of color charges suppresses all transverse color correlation beyond a distance 1/Q s 1 fm/c small enough to justify perturbative expansions in the strong coupling constant α s (Q 2 s ). (Note that the running of α s is actually more involved [116], but the basic idea remains the same.) Independently of the specific nature of the CGC model, there exist general arguments that gluon saturation has to occur for sufficiently large collision energies, as described in [117], but the evidence that this regime is already reached in present day experiments is not yet conclusive. One of its consequences, the geometric scaling of cross sections in lepton-hadron interactions, is well established, but geometric scaling seems to be a very general property of high energy interactions and thus does not constitute an unambiguous proof of saturation [118]. Keeping these caveats in mind, it is a well-justified assumption that the initial gluon distribution in high-energy heavy-ion collisions shows saturation. While this saturated state does not necessarily have to be of the universal form of the CGC, it is a reasonable expectation that its properties should be qualitatively similar to those of the CGC. It thus makes sense to use the CGC as a model of the initial state for the study of thermalization in high-energy heavy ion collisions.
Two of the present authors [119] argued that the fluctuations of the resulting energy density are surprisingly large, suggesting that the time evolution of such fluctuations could be crucial for early thermalization. This calculation was not done strictly within the CGC model, which assumes Gaussian fluctuations in the transverse color charge distribution. In contrast, the authors of [119] assumed Gaussian fluctuations in the resulting gauge fields for the purpose of analytical tractability. The CGC makes precise predictions for the transverse correlation function of gauge fields and it is this input which determines the results. Schenke et al. [120] generated fluctuating initial energy densities of colliding nuclei numerically, using the Gaussian color charge fluctuations of the CGC, but did not calculate the resulting energy density correlation function to compare with the analytical model based on Gaussian field fluctuations [119].
Some of the many models proposed to describe the surprisingly fast "hydroization" are of the type "Free-streaming plus sudden equilibration", see [121] and references given therein. The basic idea of this approach is that after the initial collision/overlap period of the two nucleons, which for the RHIC energies discussed in [121] is of order 0.1-0.2 fm/c, has defined the initial state for the fireball formation the system is still very far from hydrodynamic behaviour. Basically the kinetic energy completely dominates the dynamics of the partons. As long as this is the case free streaming should provide a good description. At some time local hydrodynamics takes over. In these models this transition is assumed to be an abrupt one, mainly for the benefit of computational ease.
Detailed simulations [25][26][27] use phenomenological models (such as the Color Glass Condensate model) to compute the initial energy deposition and to evolve it for a short amount of time (of order 0.4 fm/c) using the classical Yang-Mills equations. These studies give excellent fits to all flow coefficients by assuming that large-amplitude, small-range fluctuations subsequently evolve hydrodynamically. Note that for times short compared to the scale set by the background Yang-Mills field, typically of order 0.2 fm/c, classical Yang-Mills theory is kinetic energy dominated and should therefore be well-approximated by a free-streaming model. Alternate approaches start from fluctuations in nucleon positions and the energy deposition in individual nucleon-nucleon collisions [122,123] and use either freestreaming of particles [124] or hydrodynamics of anisotropic fluids [125,126] to bridge the gap to the onset of viscous hydrodynamics. An important open question, which motivated the present work, is to what extent this "gluing" of a phenomenological model describing the initial evolution to viscous hydrodynamics can be justified. If the agreement between free-streaming and second order hydrodynamics that we find near the end of our early-time interval really signals a cross-over (rather than an accidental temporary agreement with hydrodynamics), it provides support for this standard gluing prescription.
coordinates reads In Fourier space (for v, x, y) this has two solutions, namely Focusing on the case ω 2 < k 2 , the infalling solution is φ + . We can take k y = 0, and expand the result in a series around k x = 0 to obtain the derivative expansion in x. We find We should substitute iω → ∂ v and ik → ∂ x to get the answer in coordinate space. The result at order k 4 therefore becomes From the structure of (A.2), we see that when we expand φ + in powers of k, for each power of k 2l we get a finite series of terms in 1/r that range from 1/r 3 to 1/r l+1 . The other solution is φ − which is the advanced instead of the retarded solution. It is not difficult to see that φ − (k, ω, r) = e 2iω/r φ + (k, −ω, r) (A.5) Thus, the expansion of φ − in powers of k is similar to that of φ + , except that there is an overall prefactor of e 2iω/r which gives rise to an infinite series in 1/r. When converting to coordinate space, iω → ∂ v , so that e 2iω/r has the effect of shifting v to v + 2/r. Overall, at fixed order in spatial derivatives, φ + truncates in the 1/r expansion and φ − does not. However, we can write φ − as a truncated series acting on φ(v + 2/r, x) instead.
Therefore at fourth order in the derivative expansion in x the scalar field solution consistent with causality (obtained by demanding that the series in 1/r should truncate) is then given by The solution for the bulk metric is obtained similarly and in Eddington-Finkelstein coordinates has the form with metric components at fourth order in the derivative expansion (A.9)

B From bulk solution to boundary stress-energy tensor
Given the bulk solution for the metric and for the bulk scalar, finite physical quantities for the boundary theory can be systematically extracted through a consistent renormalization of the bulk on-shell action. Here we follow the procedure of [127]. 7 The action (2.1), in general d + 1 bulk dimensions, reads including the boundary Gibbons-Hawking term necessary to have a well defined variational problem with Dirichlet boundary conditions. γ αβ is the induced metric at the boundary, K is the trace of the extrinsic curvature of the boundary and Λ = − d(d−1)
In Fefferman-Graham coordinates the metric and the scalar solutions of the equations of motion following from (2.1) have the general asymptotic z → 0 structure with Dirichlet boundary conditions When evaluated on-shell, the action (B.1) diverges. Following [127], it can be regularized restricting the bulk integral to z ≥ and evaluating the boundary term at z = . Setting R 2 AdS = 1, the regularized action reads where we have used the relation R − 1 2 (∂φ) 2 = 2 d+1 d−1 Λ valid on-shell, Λ = − d(d−1) 2 and the Fefferman-Graham expansion (B.2). From this one can extract the IR bulk divergences (terms that diverge when 1/ is large) and construct the counterterm action.
We now specialize to the case d = 3 of interest for us here. Using the Fefferman-Graham form of the metric where dots indicate the finite part of the regularized on-shell action in the limit → 0. Following [127], a covariant counterterm action can be obtained in a minimal subtraction scheme writing the coefficients of the divergences in terms of the boundary induced metric Solving order by order in z the equations of motion following from (B.1), with Dirichlet boundary conditions (B.4), one finds where all the contractions and curvature tensors are obtained using g (0),αβ . This allows to perturbatively invert the relation (B.7) between γ αβ and (g (0),αβ , g (2),αβ ), and to rewrite the coefficients of the divergences in (B.6) in a covariant form. At the order we are interested in using which, the covariant counterterms action in the minimal subtraction scheme reads The renormalized action is then defined as The expectation values of the boundary energy momentum tensor and of the operator O associated to the massless bulk scalar are obtained from the on-shell renormalized action as and

C Asymptotic change of coordinates
Here we summarize some of the intermediate results in the change of coordinates from the Eddington-Finkelstein to the Fefferman-Graham form of the metric.
In a compact form, the metric components of the solution can be written as while for the scalar field The explicit expressions of the various coefficients can be identified by comparison with the solution worked out in Section 2.2.

(C.5)
In the new coordinates, the scalar field allows the asymptotic expansion φ(t, , χ) = φ (0) (t, χ) + φ (2) (t, χ) From the expression for the stress energy tensor (3.1) and the operator dual to the bulk scalar field (3.2), it follows that we are interested in the coefficients of metric and of the scalar field only up to g (3),αβ and φ (3) .

D The boundary stress-energy tensor
For completeness, we here give explicitly the stress-energy tensor components (3.17) up to fourth order in the gradient expansion.