Statistics of thermalization in Bjorken Flow

# Statistics of thermalization in Bjorken Flow

Jakub Jankowski111Email: jakubj@th.if.uj.edu.pl Institute of Physics, Jagiellonian University, ul. Łojasiewicza 11, 30-348 Kraków, Poland Grzegorz Plewa222Email: g.plewa@ncbj.gov.pl National Center for Nuclear Research, ul. Hoża 69, 00-681 Warsaw, Poland Michał Spaliński333Email: mspal@fuw.edu.pl
###### Abstract

The apparent early thermalization of quark-gluon plasma produced at RHIC and LHC has motivated a number of studies of strongly coupled 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 just over in units of inverse local temperature at that time. 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 .

## 1 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 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 and references see [10]). Holography represents strongly interacting Yang-Mills plasma (in 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 non-equilibrium 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 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 anisotropy 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 horizon444The boundary-covariant apparent horizon consistent with conformal symmetry is unique in fluid-gravity duality[19]. 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. Section 2 recalls 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. 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.

## 2 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=τcoshy ,z=τsinhy , (1)

all observables are independent of the rapidity . 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]

 Tμν=diag(ϵ,pL,pT,pT) , (2)

where

 pL=−ϵ−τ˙ϵ ,pT=ϵ+12τ˙ϵ , (3)

so it is determined completely by the energy density . Furthermore , and for SYM [20] in equilibrium

 ϵ=38π2N2cT4 . (4)

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 , since this introduces no ambiguity555Note that it may not be possible to do this for more general flows [22]..

In a hydrodynamic state, by definition, the energy momentum tensor can be expressed in terms of the local fluid velocity and temperature . Specifically, it can be written in the form

 ⟨Tμν⟩=ϵuμuν+p(uμuν+ημν)+… (5)

where the ellipsis denotes terms involving gradients of the hydrodynamic variables. The energy density and pressure are understood to be expressed in terms of the temperature using the equations of state. For a conformal field theory one has .

A boost invariant configuration has , 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

 T(τ) = Λ(Λτ)1/3{1−16π(Λτ)2/3+−1+log236π2(Λτ)4/3+ (6) + −21+2π2+51log2−24log221944π3(Λτ)2} .

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(w)=τwdwdτ , (7)

where is the proper time in units of inverse effective temperature. In a conformal theory the behaviour of this function at large is independent of for dimensional reasons. By computing it at late times (using (6)) one finds

 fH(w)=23+19πw+1−log227π2w2+15−2π2−45log2+24log22972π3w3+… . (8)

A given solution, , starts out far from equilibrium, but the damping of non-hydrodynamic modes ensures that after a sufficiently long time only hydrodynamic modes remain. Given the result of a numerical simulation for one can check whether the hydrodynamic regime has been reached by comparing 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

 s=12π2N2cT3 . (9)

In a boost invariant setting it is useful to speak of entropy per unit rapidity and transverse area, given by

 dSdydx2⊥=τs , (10)

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

 S=112π2N2c⋅(Ti)2(dSdydx2⊥) , (11)

where is the effective temperature at the initial time. At late times, the ratio (11) approaches the limiting value

 Sf=Λ2(Ti)−2 . (12)

Henceforth when speaking of entropy, we will mean the quantity given in eq. (11).

## 3 The gravity dual

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 generic, non-equilibrium initial states of Yang-Mills plasma at . In the holographic approach pursued here this is implemented by 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 adopted here draw on the experience gained from references [3, 5] and make it possible to effectively solve the constraints at , 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 [25]), we adopt Eddington-Finkelstein coordinates. The metric is parameterized as

 ds2=1z2(−2dzdt−a(t,z)dt2+b(t,z)dy2+c(t,z)dx2⊥) . (13)

This parameterization is different from that used in [3], which is inconvenient if one wishes to set initial conditions at . With the parameterization (13) it is possible to set initial conditions at any . Note that on the conformal boundary this time coordinate coincides with the proper time coordinate .

Assuming the metric Ansatz given in eq. (13), the Einstein equations with negative cosmological constant

 Rab+4gab=0 (14)

reduce to a set of 5 nonlinear partial differential equations for the metric coefficients . Three of these equations can be regarded as evolution equations. Of the remaining two, one is a constraint on initial data (and will be discussed in detail in section 4), 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 a 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

 a(z,t) = 1+z3 ~a(z,t) b(z,t) = (t+z)2+z3 ~b(z,t) c(z,t) = 1+z3 ~c(z,t) . (15)

The evolution equations expressed in terms of are amenable to numerical computations at , but at one has to set

 a(z,t) = 1+z3 ~a(0)(z,t) b(z,t) = (t+z)2+z5 ~b(0)(z,t) c(z,t) = 1+z3 ~c(0)(z,t) (16)

in order to obtain differential equations which can be used to make the initial time step at .

With the definitions (3)-(3) the evolution equations assume a rather special form: even though they are nonlinear, they are linear in the time derivatives666This is of course a consequence of using the Eddington-Finkelstein gauge. . Specifically, letting denote the 4-component vector one finds that the 4 independent Einstein equations can be written as

 ∂zu=Au+f , (17)

where is a 4x4 matrix and is a 4-component vector expressed in terms of and their derivatives. Solving this linear system of first order ordinary differential equations by integration over one can obtain the time derivatives , and itself. These can then be used to evolve the metric. The explicit form of and is unilluminating, but straightforward to obtain. As mentioned earlier, for the time step at one needs to use (3) instead of (3).

The integration over can conveniently be carried out using the spectral representation [26], 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 . The result of this standard procedure is

 a(t,z) = 1+a4(t)z4+… b(t,z) = t2+2tz+z2+z4 t2(a4(t)+34ta′4(t))+… c(t,z) = 1−12z4(a4(t)+34ta′4(t))+… . (18)

The above expansion is expressed in terms of the single function , 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 [27, 28, 29] 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 [28]

 ϵ(τ)=−3N2c8π2a4(τ) . (19)

It can be checked that the field equations imply that at early times the expansion of contains only even powers of , which due to (19) implies an analogous property of the energy density – this was established already in [13] using Fefferman-Graham coordinates.

Time evolution is computed using the 4th order Adams-Bashforth-Moulton predictor-corrector approach. This turned out to be preferable to using a Runge-Kutta integrator, due to the relatively large cost of spacial integrations and solving eq. (17) at each time step.

In practice setting up the numerical solution of Einstein equations is somewhat subtle here, since one must choose a suitable range for the coordinate: ideally, it should cover the region starting at the boundary and reaching just inside of the event horizon. The location of the latter is however unknown until the entire spacetime is determined. This chicken-and-egg problem is resolved in the spirit of what is universally done in the numerical relativity literature: the apparent horizon is located (if present), and the coordinate range is terminated there. The position of the apparent horizon is computed directly from the definition in terms of the expansion scalars (see for example [17]). For our purposes the area of the apparent horizon is also of prime interest, since it is connected to the entropy of the black brane (and ultimately, of the Yang-Mills theory plasma).

The position of the apparent horizon is not static, so one has to dynamically adjust the range of every few time steps. Whenever this is done, one has to translate the solution at that time to the new spacial grid. The heuristics for this regridding are essential for the simulation to be stable long enough for hydrodynamics to be reached. The accuracy of the numerical evolution can in practice be monitored by checking whether the constraint equation (which is satisfied at the initial time) remains satisfied after each time step is taken.

## 4 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 the metric coefficients and

 12(∂zbb)2−∂2zbb+12(∂zcc)2−∂2zcc=0 . (20)

At we choose the function and then obtain the initial condition for by numerically solving the above ordinary differential equation. Finally, the metric coefficient is obtained by solving the linear system (17).

The choice of has to be consistent with the near boundary expansion, whose leading terms are given in eq. (3). In practice we chose various families of functions (such as ratios of polynomials, elementary as well as special functions) depending on some number of arbitrary coefficients which were then randomly generated a number of times.

Note that the choice of at also determines the expansion of in powers of ; this happens, because the near boundary expansion of is expressed in terms of and its derivatives. This series representation of , valid close to , 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 can only contain even powers of . It is not however clear whether one should assume that the energy density is a nonzero constant at . The color glass condensate approach suggests that (see for example [30]), which was assumed in the present study. It could in principle be that , in which case the leading term would be

 ϵ(τ)=ϵ2kτ2k+… (21)

with for some . It has been suggested [31, 32] that such initial states (with ) 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 non-vanishing at the initial time. With this assumption one can for convenience scale it to any given constant. We chose to fix . Using the relations (4) and (19) this leads to

 Ti=1π . (22)

While the initial geometries generated are very diverse, the states of 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 non-equilibrium 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 non-decreasing, 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 or very shortly thereafter. Using the Bekenstein-Hawking relation to associate entropy with the apparent horizon area one finds, with the assumed normalizations,

 S=aAHπ , (23)

where 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.

## 5 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 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 referred as reheating.

Since we have normalized the effective temperature to at 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 or higher. There seems to be no functional dependence of on .

Despite these qualitative differences at early times, all numerical solutions exhibit cooling at intermediate and late times. In the earlier studies of boost invariant flow [4, 5] it was found that hydrodynamics become a good description at some in all the 29 cases considered. In the study presented here evolution was carried on until . 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 we evaluate the corresponding function defined by (7). We then calculate the difference between this and the asymptotic form given by (8). The thermalization “time” is defined as the maximum value of such that for

 ∣∣∣f(w)fH(w)−1∣∣∣<0.005 . (24)

The distribution of the values of found can be seen in the histogram shown in Fig. 3. It is clear from this that the early thermalization in boost invariant flow observed in [4] is really a generic feature of the dynamics.

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.

At late times all the solutions converge to hydrodynamics, as expected. This can be seen on plots showing the numerically computed function , or equivalently by plotting the pressure anisotropy

 Δ≡pT−pLϵ=6f(w)−4 . (25)

Typically this is quite large (on average ) at thermalization, as seen in Fig. 5. The conclusion is then that hydrodynamics becomes valid rather early, and the large pressure anisotropy 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 [33].

After thermalization is attained the hydrodynamic description becomes valid (by definition). As discussed in section 2, the boost-invariant hydrodynamics of a conformal fluid is determined by the single scale . For any numerical solution this value can be found by fitting formula (6) to the tail of the numerical data. A histogram of the values of obtained in this way is shown in Fig. 6.

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 the final entropy (which is here), this result suggests that this multiplicity is related to the square of the initial entropy777Of course, unless one understands the meaning of initial entropy directly in field theory terms, this observation is not very useful..

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 [34]. Having computed the scale , one can rescale the proper time and effective temperature in such a way that the rescaled data reach hydrodynamics with . 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.

## 6 Entropy

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 non-equilibrium 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

 Sf=π2Λ2. (26)

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 , implies that for small initial entropies there should be a correlation between the entropy produced, , and the initial entropy , which can be fitted by a quadratic function. This is indeed seen in the data obtained in our simulations (see Fig. 8) for evolution from initial states with entropy below . 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 .

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 [35, 36, 37, 38, 39].

In the following, by entropy at thermalization (denoted by ) 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

 S(0)=τπ2T3 , (27)

assuming the normalization (22). Gradient corrections to hydrodynamic entropy are then quantified by the ratio . 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 anisotropy 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.

Entropy production at different stages of evolution is illustrated in Fig. 9. The left panel shows the ratio of entropy produced during the hydrodynamic stage 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 however, 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. 9. This quantity is more in line with expectations, although it can still be quite significant: about of the total entropy is produced by the dissipative effects during hydrodynamic evolution.

## 7 Conclusions

We have analyzed the time evolution of 600 randomly generated far from equilibrium initial states of strongly coupled SYM plasma well into the hydrodynamic regime. It is 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.

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. One possible limitation is that some randomly generated initial geometries do not evolve in a stable way for long enough to reach the hydrodynamic regime. This is most likely because the heuristics we use to estimate the position of the event horizon are not always effective, especially in cases where no apparent horizon exists on the initial time slice. This could be the reason why among the initial states found in [4] there are some with lower initial entropy than what we have found.

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 anisotropy is still significant (with of the order of 0.3). 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 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 models of the initial post-collision state. Interesting attempts going in this directions have already been made in [40, 41]. The connection between the notion of dynamical entropy defined in these references and apparent horizon entropy is not clear, but could perhaps be made via an approach such as [42].

Apart from pursuing answers to the questions discussed above, the work reported here could be generalized in various ways. One direction would be to consider an alternative class of initial conditions, with different early time behaviour. Another natural avenue to pursue is the behaviour of non-local observables (such as Wilson lines or entanglement entropy [43]) in the expanding backgrounds described here. We hope to return to these issues in the future.

Acknowledgments: The authors would like to thank David Blaschke, Michał P. Heller, Romuald Janik, Daniel Nowakowski, Przemysław Witaszczyk and Wilke van der Schee for valuable discussions. This work was partially supported by Polish National Science Center research grant 2012/07/B/ST2/03794 (GP and MS) and a post-doctoral internship grant No. DEC-2013/08/S/ST2/00547 (JJ).

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters