Statistics of thermalization in Bjorken Flow

The apparent early thermalization of quark-gluon plasma produced at RHIC and LHC has motivated a number of studies of strongly coupled ${\mathcal N}=4$ supersymmetric Yang-Mills theory using the AdS/CFT correspondence. Here we present the results of numerical simulations of Bjorken flow aimed at establishing typical features of the dynamics. This is done by evolving a large number of far from equilibrium initial states well into the hydrodynamic regime. The results strongly suggest that early thermalization is generic in this theory, taking place at proper times around $0.6$ in units of inverse local temperature. We also find that the scale which determines the rate of hydrodynamic cooling is linearly correlated with the entropy of initial states defined by the area of the apparent horizon in the dual geometry. Our results also suggest that entropy production during the hydrodynamic stage of evolution is not negligible despite the low value of $\eta/s$.


Introduction
Experimental studies of quark-gluon plasma at RHIC and LHC have led to a number of theoretical challenges. The successful phenomenological description of what is observed rests on the assumption that the system thermalizes on a time scale of order less than inverse local temperature. By thermalization we mean here that a hydrodynamic description becomes valid [1,2]. It has recently been appreciated that at least in some model cases this happens while the system is still very anisotropic [3,4], so instead of using the term thermalization one sometimes speaks of "hydrodynamization".
It is known that in a strongly interacting system the approach to hydrodynamics can be very rapid [4,5] -as rapid as seen in heavy ion collisions. One of the purposes of the present work was to see how generic this behaviour is. This type of question was recently considered in the context of isotropization [6], where no hydrodynamic modes are present. Here we consider a case where we can see the onset of hydrodynamic behaviour (as in [4,5]). We consider a large number of far from equilibrium initial states (which mimic the situation right after heavy nuclei collide) and simulate the subsequent evolution until hydrodynamic behaviour sets in. This way we can see, in particular, the distribution of thermalization times.
While the theory underlying the phenomena studied at RHIC and LHC is firmly believed to be textbook QCD, it has so far not been possible to understand the observations from first principles. Since the post-collision plasma state appears to be strongly coupled, much effort has been devoted to consideration of similar physical issues in the context of strongly coupled N = 4 SYM [7], which has a tractable holographic representation in terms of AdS gravity [8,9]. Even then it has not been possible to describe the entire process analytically, making it necessary to resort to suitable numerical techniques.
Numerical studies of the dynamics of Yang-Mills plasma based on the AdS/CFT correspondence have been an area of significant activity in recent years (for a review

INTRODUCTION
and references see [10]). Holography represents strongly interacting Yang-Mills plasma (in d = 3 + 1 dimensions) by a spacetime geometry in one dimension higher.
Initial states of the plasma are represented as initial geometries which are then evolved numerically using Einstein's equations. These computations are in general rather daunting, so much effort has been devoted to situations where symmetries reduce the complexity to a manageable level. A physically important case is that of Bjorken flow [11]. Despite considerable simplification of the dynamics, this example is rich enough to be of phenomenological significance, even though its usefulness is limited to the central rapidity region.
Studies of Bjorken flow in the context of AdS/CFT were initiated in [12], which showed analytically that hydrodynamic behaviour appears at late times. The dynamics at early times was considered soon after [13], but it quickly became clear that a numerical approach to the problem was required. The first numerical studies were pioneered by Chesler and Yaffe [3], but rather than setting general initial conditions (corresponding to nonequilibrium states of the dual gauge theory plasma) in the bulk the authors considered perturbing the system at the boundary.
Formulating consistent initial conditions requires solving the constraints contained in Einstein equations. This rather nontrivial problem was considered in [13], where 29 solutions were found. These initial conditions were the starting point for the numerical studies [4,5], where a number of important observations were made. Some of the results were however limited by the fact that only a small number of initial conditions were available.
The present work follows [14,3] in using Eddington-Finkelstein coordinates. The parameterization of the metric used here is different however, because the original approach of [3] made it difficult to set the initial conditions at very early proper time. The approach described here overcomes this difficulty. The form of the field equations makes it possible to effectively solve the constraints. Initial conditions can be set by choosing a single independent metric coefficient function which satisfies

INTRODUCTION
some mild conditions and the remaining metric coefficients on the initial time slice are then determined by numerically integrating ordinary differential equations.
Given a well posed initial value problem it is possible to solve the field equations numerically. The complications of AdS asymptotics can be dealt with following the approach of [3] (see also [15]). Once the spacetime geometry is determined (up to some final time) the AdS/CFT prescription provides an explicit recipe to determine the expectation value of the energy-momentum tensor of the dual 4-dimensional gauge theory.
At sufficiently late times the system is governed by the equations of hydrodynamics -this was shown analytically in [12] and numerically in [4,5]. The determination of the "typical" time when hydrodynamic behaviour sets in is one of the main tasks we address. Heller et al. [4,5] analysed 29 initial conditions which indicated the hydrodynamics sets in quite early and noted that the pressure asymmetry at that time is still significant (this was also pointed out in a similar context by [3]).
The approach adopted in this paper makes it possible to generate consistent initial conditions at will, so we are able to look at a large number of solutions. This makes it possible to claim a much more firm estimation. In what follows we describe the distribution of thermalization times for a sample of over 600 initial conditions. An important role in the interpretation of our simulations is played by entropy.
After thermalization the hydrodynamic notion of entropy [16,17,18] is valid, but also at earlier times the apparent horizon 1 area (mapped to the boundary) gives rises to an operational notion of entropy, which asymptotes to hydrodynamic entropy at late times.
The organization of the paper is as follows. In section 2 we briefly recall the implications of boost invariance and the hydrodynamics of Bjorken flow. In section 3 we present the dual, gravitational description and some details of the numerics.
The initial states for the simulation are discussed in section 4. The results for thermalization are presented in section 5 and entropy production is the subject of section 6.

Boost invariant flow
An idealized, but very useful description of a heavy ion collision was formulated by Bjorken [11]. The colliding nuclei are regarded as infinite sheets of mater extending in the plane transverse to the collision axis. The collision energy is also infinitely high, and the physics of the system is assumed to be determined only by the proper time evolution of the energy density. In the case of central, ultrarelativistic collisions this is a useful approximation for the description of observables in central rapidity region.
Boost invariant states in four spacetime dimensions can most simply be characterized by saying that in proper time-rapidity coordinates t = τ cosh y , z = τ sinh y , all observables are independent of the rapidity y. Furthermore, we assume independence of the transverse coordinates.
In particular, in a conformal theory the expectation value of the energy momentum tensor in a boost-invariant state takes the form [12] where so it is determined completely by the energy density (τ ). Furthermore ∼ T 4 , and in N = 4 SYM [20] It is very convenient to introduce the notion of effective temperature, which satisfies eq. (4) in any state; thus it is the temperature of an equilibrium state with the same energy density [21]. In what follows we will denote the effective temperature by T , since this introduces no ambiguity 2 .
In a hydrodynamic state, by definition, the energy momentum tensor can be expressed in terms of fluid velocity u and local temperature T . Specifically, it can be written in the form where the ellipsis denotes terms involving gradients of the hydrodynamic variables.
The energy density and pressure p are understood to be expressed in terms of the temperature using the equations of state. For a conformal field theory one has = 3p ∼ T 4 .
A boost invariant configuration has (u µ ) = (−1, 0, 0, 0), and all the physics is encoded in the dependence of the energy density (or, equivalently, the temperature) on τ . Analytic calculations have been performed up to third order in the gradient expansion [23,21,17], the result being The scale Λ appearing in (6) depends on the initial conditions. Indeed, it is the only trace of the initial conditions to be found in the late time behaviour of the system.
It is expected that starting from any initial state the system should evolve in such a way that at late times it will be described by (6) for some value of the constant Λ.
As we shall see, our simulations suggest that for "low entropy" initial states Λ is proportional to the apparent horizon entropy of the initial state.
To determine whether that the system has thermalized (in the sense of reaching hydrodynamics) it is very convenient to define the function [4] f where w = τ T (τ ) is the proper time in units of inverse effective temperature. In a conformal theory the behaviour of this function at large w is independent of Λ for dimensional reasons. By computing it at late times (using (6)) one finds A given solution, (τ ), starts out far from equilibrium, but the damping of nonhydrodynamic modes ensures that after a sufficiently long time only hydrodynamic modes remain. Given the result of a numerical simulation for T (τ ) one can check whether the hydrodynamic regime has been reached by comparing f (w) calculated from eq. (7) with the hydrodynamic form given by eq. (8).
An important role in the following is played by entropy of the expanding plasma.
The thermodynamic (equilibrium) entropy density is given by In a boost invariant setting it is useful to speak of entropy per unit rapidity and transverse area, given by where the factor of τ reflects the expanding volume. In our holographic calculation the entropy per unit rapidity and transverse area can be computed from the apparent horizon area (due to the Bekenstein-Hawking relation), as discussed in section 6.
As discussed in [5] it is convenient to use the dimensionless ratio where T i is the effective temperature at the initial time. At late times, the ratio (11) approaches the limiting value Henceforth when speaking of entropy, we will mean the quantity given in eq. (11).
The the standard AdS/CFT prescription for calculating the energy-momentum tensor expectation value in a given state is to find the gravity solution corresponding to that state and then extract the expectation value from the near-boundary asymptotics.
The object of this paper is to simulate the approach to equilibrium starting from a randomly chosen, far from equilibrium initial state of Yang-Mills plasma at τ = 0. This implies choosing an initial geometry at random and evolving it using Einstein's equations. The choice of initial geometry must satisfy the constraint equations, whose form depends on the choice of coordinates and parameterization of the metric. The choices made here build on those used in [3,5] and make it possible to effectively solve the constraints at τ = 0 so that an arbitrary number of consistent initial geometries can be generated and evolved [24].
As in the work of Chesler and Yaffe [3,10] (as well as the papers on fluid-gravity duality), we adopt Eddington-Finkelstein coordinates. The metric is parameterized as This parameterization is different from that used in [3], which is awkward if one wishes to set initial conditions at t = 0. With the parameterization (13) it is possible to set initial conditions at any t. Note that on the conformal boundary this time coordinate coincides with the proper time coordinate τ .
Assuming the metric Ansatz (13), the Einstein equations lead to a set of 5 nonlinear PDEs for the metric coefficients a, b, c. Three of these equations can be regarded as evolution equations. Of the remaining two, one is a constraint on initial data, and the other is redundant in the sense that it is automatically satisfied by any consistent initial data evolved using the evolution equations.
Due to the presence of negative cosmological constant the solutions describe locally asymptotically AdS geometries. For the purpose of numerical computations it is very convenient to isolate the asymptotic AdS behaviour by defining It turns out that this form is suitable at t = 0, but at t = 0 one has to use in order to obtain a non-singular set of differential equations. These will be used to make the initial time step at t = 0.
With the definitions (15)-(16) the evolution equations assume a rather special form: even though they are nonlinear, they are linear in the time derivatives.
Specifically, letting u denote the 4-component vector (∂ tb , ∂ tc , ∂ zã ,ã) one finds that the 4 independent Einstein equations can be written as where A is a 4x4 matrix and f is a 4-component vector expressed in terms ofã,b,c and their z derivatives. Solving this linear system of first order ordinary differential equations by integration over z one can obtain the time derivatives ∂ tb , ∂ tc andã itself. These can then be used to evolve the metric. The explicit form of A and f is unilluminating, but straightforward to obtain. As mentioned earlier, for the time step at t = 0 one needs to use (16) instead of (15).
The integration over z can conveniently be carried out using the spectral representation [25], which we use in the form of the collocation method with Chebyshev polynomials. For regular configurations, such as those considered here, it gives exponential convergence, allowing for modest grid sizes.
To solve eq. (17) one must impose the correct boundary conditions. These can be determined by solving the field equations in a near boundary expansion, i.e. an expansion around z = 0. The result of this standard procedure is The above expansion is expressed in terms of the single function a 4 (t), which is not determined by the near-boundary analysis of the equations of motion. This function is determined by the asymptotics of the full numerical solution. Holographic renormalisation [26,27,28] is then invoked to calculate the energy-momentum tensor (by construction this is of the form given in eq. (2)- (3)). In this way from the numerical solution at a given time τ one can read off the energy density using the relation [27] (τ ) = − 3N 2 c 8π 2 a 4 (τ ) .
It can be checked that the field equations imply that at early times the expansion of a 4 (t) contains only even powers of t, which due to (19) implies the same property of the energy density -this was established already in [13] using Fefferman-Graham coordinates.
The evolution in time is computed using the 4th order Adams-Bashforth-Moulton predictor-corrector approach. This turned out to be preferable to using a Runge-

Initial states
To specify a consistent initial condition for the asymptotically AdS geometry one needs to satisfy the single constraint equation among (14). This equation involves only b and c 1 2 We choose the function b and then obtain c by numerically solving the above ordinary differential equation. Finally, the metric coefficient a is obtained by solving the linear system (17).
The choice of b has to be consistent with the near boundary expansion, whose leading terms are given in eq. (18). In practice we chose various families of functions (such as ratios of polynomials) depending on some number of coefficients which were then randomly generated a number of times.
Note that the choice of b also determines the expansion of a 4 (t) in powers of t; this happens, because the near boundary expansion of b is expressed in terms of a 4 (t) and its derivatives. This series representation of a 4 (t) can be used to control the numerical solution at early times.
As discussed in the previous section, the initial energy density when expanded in powers of the proper time around τ = 0 can only contain even powers of τ . It is not however clear whether one should assume that the energy density is a nonzero constant at τ = 0. The color glass condensate approach suggests that (0) = 0 > 0, which was assumed in the present study. It could in principle be that 0 = 0, in which case the leading term would be with 2k = 0 for some k > 0. It has been suggested [29,30] that such initial states (with k = 1) appear in consequence of shock wave collisions in AdS. We hope to return to these kind of initial states in the future, but in the work reported here we assumed that the energy density is nonvanishing at the initial time. With this assumption one can for convenience scale it to any given constant. We chose to fix a 4 (0) = −1. Using the relations (4) and (19) this leads to While the initial geometries generated are very diverse, the states of N = 4 SYM which they map to all have the initial energy density scaled to the same value.
These initial states of the gauge theory are nevertheless distinct, as their subsequent evolution shows.
A very interesting question is how these states should be characterized in the Yang-Mills theory. Since these are highly nonequilibrium states one cannot use thermodynamic notions to describe them. Indeed, it seems that one should not expect any simple description. However, as noted by Heller et al. [4], from a holographic perspective there is one characteristic which seems to be useful and interesting: it is "initial entropy". This quantity can be defined unambiguously on the gravity side in terms of the apparent horizon area of the black brane geometry mapped to the boundary along null geodesics. This quantity is guaranteed to be nondecreasing, and for large times becomes identical to hydrodynamic (and ultimately thermodynamic) entropy. It is not a priori obvious that an apparent horizon should always exist in the chosen time slicing of AdS, but in all states considered here it exists either at the initial time τ = 0 or very shortly thereafter. Using the Bekenstein-Hawking relation to associate entropy with the apparent horizon area one finds, with the assumed normalizations, where a AH is the area of the apparent horizon. For all the numerical solutions considered this indeed coincides with the thermodynamic result (6) at late times.
The range of entropies of all the initial states generated is shown in Fig. 1.

Thermalization
The results presented in this and the following sections follow from picking initial gravity configurations as described above and evolving them until times large enough for hydrodynamics to be valid.
For each initial condition the simulation gives the proper-time dependence of the effective temperature. Similarly to the findings of [5] we find that there are generically three different behaviours. For states of initial entropy below 0.3 typically the temperature either decreases monotonically from the start, or after an initial plateau. For states with higher entropy the effective temperature rises at first and only after some time the system begins to cool. In [5] this phenomenon was refereed as reheating.
Since we have normalized the effective temperature to 1/π at t = 0 for all the initial configurations and for large times the system is cooling in accordance with hydrodynamics, the maximal temperature is a measure of the reheating phenomenon. Fig. 2 shows that while reheating is absent for almost all low entropy initial states, it is typical for entropies of around 0.3 or higher. There seems to be no functional dependence of T max on S i .  [4,5] it was found that hydrodynamics become a good description at some w < 1 in all the 29 cases considered. In the study presented here evolution was carried on until w = 1.2. In all 600 cases examined here hydrodynamics became an accurate description significantly earlier, in agreement with [4,5].
The criterion for thermalization which we used is the same as that proposed in [4].
Specifically, for a given solution T (τ ) we evaluate the corresponding function f (w) defined by (7). We then calculate the difference between this and the asymptotic form given by (8). The thermalization "time" w H is defined as the maximum value The distribution of the values of w H found can be seen in the histogram shown in  The corresponding thermalization times are plotted against initial entropy in Fig. 4. The units are those of initial temperature in order to easily compare with reference [4], where it was observed that for the sample of initial states considered in that work, there appeared to be a correlation between the initial entropy and the thermalization time. In the larger set of solutions examined here this does not appear to be the case.
Typically this is quite large ∆ ∼ 0.7 at thermalization, as seen in Fig. 5. The conclusion is then that hydrodynamics becomes valid rather early, and the large pressure asymmetry predicted by hydro at such early times is in fact reliable. The oscillations visible on these plots as hydrodynamics is approached indicate the presence of quasinormal modes [31].
histogram of the values of Λ obtained in this way is shown in Fig. 6. In QCD this would be of direct physical interest, since the final charged particle multiplicity is proportional to the square of Λ.
The distribution seen in the left panel of Fig. 6 is reminiscent of the histogram of initial entropies. In fact, the simulations reported here suggest that the scale Λ is linearly correlated with the initial entropy. This is seen in the right panel of Fig. 6.
Since in QCD the final charged particle multiplicity is proportional to Λ 2 , this result suggests that this multiplicity is related to the square of the initial entropy.
One can see that the linear correlation between Λ and the initial entropy is valid for a large range of initial entropies, but at the high end the correlation seems to be lost. This suggests that in this regime the initial entropy is not the dominant characteristic of the initial state, and other factors come into play.
Finally, it is instructive to view the hydrodynamization process from a slightly different perspective [32]. Having computed the scale Λ, one can rescale the proper time and effective temperature in such a way that the rescaled data reach hydrodynamics with Λ = 1. Such a rescaling amounts to normalizing the data so all solutions have the same final entropy. If this rescaling is done for all solutions they will start out at various values of the initial effective temperature and will all converge at late times. In Figure 7 we see this for a few sample configurations. The observed viscosity of quark-gluon plasma implies that entropy rises during hydrodynamic evolution, as well as in the pre-hydrodynamic stage -to the extent that it makes sense to speak of entropy in a highly nonequilibrium system. In the case when a holographic dual is available we do have a notion of entropy as defined by the apparent horizon area, so we can quantify entropy production at all times (at least in those cases where an apparent horizon exists and was identified within the computational grid).
Entropy production is an important aspect of the dynamics -it is related to particle production in heavy ion collisions. An interesting observation was made in reference [4], whose authors pointed out that the entropy produced is correlated with the initial entropy. In fact, for Bjorken flow the final entropy is given in terms of the scale Λ by eq. (12), so using eq. (22) one has As explained earlier, the actual value of Λ is obtained by fitting the 3rd order formula (6) to the tail of the numerical solution. This relation, together with our previous observation of a linear correlation between Λ and the initial entropy S i , implies that for small initial entropies there should be a correlation between the entropy produced, S f − S i , and the initial entropy S i , which can be fitted by a quadratic function.
This is indeed seen in the data obtained in our simulations (see Fig. 9) for evolution from initial states with entropy below 0.3. For initial states with higher entropy the correlation is lost. This is of course consistent with our earlier observations concerning the appearance of a "chaotic phase" when the initial entropy exceeds 0.3.
It is also interesting to consider how much entropy is produced during the hydrodynamic stage of evolution. Since hydrodynamics provides an accurate description of the dynamics already at rather early times, it is to be expected that a significant part of the final entropy is due to hydrodynamic evolution -more than often assumed. This can be taken as motivation for efforts to resum the hydrodynamic series [33,34,35,36,37].
In the following, by entropy at thermalization (denoted by S th ) we mean the apparent horizon entropy given by eq. (23) evaluated at thermalization time. It is interesting to estimate the size of gradient corrections to thermodynamic entropy at this stage of evolution. This can be done by comparing the "exact" entropy given in terms of the apparent horizon area by eq. (23) with hydrodynamic entropy defined by the entropy current within the gradient expansion [17]. Specifically, using eq. (9) and (10) one finds that the leading order (thermodynamic) approximation to the entropy defined in eq. (11) is assuming the normalization (22). Gradient corrections to hydrodynamic entropy are then quantified by the ratio (S − S (0) )/S. In the simulations described here we find that at thermalization time this ratio assumes values in the range 0.05 -0.15. This is consistent with the fact that corrections to the leading order term, eq. (27), first appear at 2nd order in the gradient expansion [17]. For this reason, even though the pressure asymmetry is large at thermalization time, the gradient corrections to thermodynamic entropy are relatively small. Therefore, in the situation studied here, thermodynamic entropy is a reasonable approximation as soon as hydrodynamics becomes valid. However, this conclusion need not hold when conformal symmetry is broken, since one then expects corrections to the entropy current already at first order in the gradient expansion.
Entropy production at different stages of evolution of the system is illustrated in Fig. 8. The left panel shows the ratio of entropy produced during the hydrodynamic stage of evolution to the total entropy production. This can be very high, especially for states of low initial entropy, for which often more than half of the entropy produced comes from the hydrodynamic evolution.
Usually one compares the entropy produced during the hydrodynamic stage not to the entropy produced during complete evolution, but rather to total entropy.
This ratio is plotted in the right panel of Fig. 8. This quantity is more in line with expectations, although it can still be quite significant: about 15% − 25% of the total entropy is produced by the dissipative effects during hydrodynamic evolution.  Apart from the motivation coming from heavy ion collisions, the behaviour of entropy as equilibrium is approached is an interesting issue in its own right. It has been speculated in the past that some suitable notion of entropy grows linearly until close to the time when hydrodynamic behaviour sets in. In the simulations described here there does seem to be an interval where this growth is linear (for small t), but this ends well before hydrodynamics becomes valid.

Conclusions
We have analyzed the time evolution of 600 randomly generated far from equilibrium initial states of strongly coupled N = 4 SYM plasma well into the hydrodynamic regime. Having a fairly large sample makes it possible to look for generic features of the dynamics. While the selection of initial states does not have any built-in bias that we are aware of, we cannot claim that the set we used is actually representative.
It is also important that in the approach presented here the number of initial states is not limited by anything else than patience, and more solutions could easily be generated.
The main motivation for this study was to see how generic is the early thermalization observed in the limited sample of initial states considered in [5]. Our results indicate that early thermalization is the norm. Hydrodynamic behaviour typically sets in at a time of the order of 0.6 over the effective temperature, when the pressure asymmetry is still significant (with ∆ of the order of 0.7). At the onset of thermodynamics one can estimate the entropy in field theory using the leading order hydrodynamic approximation, which we have verified to be reasonably accurate, as expected [17]. This serves as a reference point to quantify the amount of entropy produced during the hydrodynamic stage of evolution -the result is that the increase in entropy is typically about 20% of the final entropy, which amounts to about 40% of the total entropy produced during evolution.
In the previous study by Heller et al. [4] it was pointed out that in the context of holography one could try to characterize far from equilibrium initial states by their entropy measured by the area of the apparent horizon appearing in the dual geometry. Our numerical experiments confirm that for a range of initial entropies this is indeed a useful concept. One of the surprises that have emerged was that the parameter Λ which sets the scale for the hydrodynamic tail is linearly correlated with this initial entropy. However, for large entropies (larger than 0.3 in our units) this relation is lost and some "chaotic" behaviour develops. This indicates that for such states there are other characteristics which are important. Even for initial entropies which are not so high, it is not true that they fully characterize the initial state, because thermalization times are not correlated with them at all.
The fact that initial entropy seems to be a useful concept raises the question of its field theory meaning. In other words, how can one define it without appealing to the holographic picture? Perhaps this could be addressed in the context of theoretical models of the initial post-collision state based on QCD, such as models implementing the idea of a color glass condensate [38]. While a priori these far from equilibrium initial states may have no straightforward characterization, the fact that in this approach the initial states are characterized by the single effective gluon saturation scale Q s offers hope that a simple picture may exist.