Bulk Viscosity and Cavitation in BoostInvariant Hydrodynamic Expansion
Abstract:
We solve second order relativistic hydrodynamics equations for a boostinvariant dimensional expanding fluid with an equation of state taken from lattice calculations of the thermodynamics of strongly coupled quarkgluon plasma. We investigate the dependence of the energy density as a function of proper time on the values of the shear viscosity , the bulk viscosity , and second order coefficients, confirming that large changes in the values of the latter have negligible effects. Varying the shear viscosity between zero and a few times , with the entropy density, has significant effects, as expected based on other studies. Introducing a nonzero bulk viscosity also has significant effects. In fact, if the bulk viscosity peaks near the crossover temperature to the degree indicated by recent lattice calculations in QCD without quarks, it can make the fluid cavitate — falling apart into droplets. It is interesting to see a hydrodynamic calculation predicting its own breakdown, via cavitation, at the temperatures where hadronization is thought to occur in ultrarelativistic heavy ion collisions.
1 Introduction
In recent years, the comparison between data from experiments at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory on the transverse expansion of the matter produced in ultrarelativistic nucleusnucleus collisions with nonzero impact parameters [1, 2] and calculations done using secondorder relativistic viscous hydrodynamics [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] have strengthened the case that the quarkgluon plasma in QCD at temperatures above, but not too far above, the crossover from a hadron gas is a strongly coupled liquid. These comparisons indicate a shear viscosity to entropy density ratio that is within a factor of a few of [15, 14], which is the value of this ratio in any gauge theory with a dual gravity description in the limit of infinite coupling and infinitely many colors [16, 17, 18]. The degree of success of the hydrodynamic description of these collisions, which turn out to be creating exploding droplets of a fluid that is closer to the ideal liquid limit than water is by about two orders of magnitude — and water is the liquid that hydrodynamics is named after — have refocused attention on the question of when a hydrodynamic description applies and when it breaks down. The aspect of this question that has drawn most attention is “Before what early time in the collision is a hydrodynamic description invalid?” This is the question of how, and how quickly, approximate local thermal equilibrium is attained. We shall focus, instead, on the complementary question “Assuming that a hydrodynamic description is valid starting at some early time, after what late time does this hydrodynamic description break down?”
Ultimately, heavy ion collisions produce thousands of hadrons flying outwards toward the detector. The question of how, and when, a hydrodynamic description breaks down and how one then ends up with a cloud of hadrons flying apart is the question of how, and when, “freezeout” occurs. In calculations, freezeout is typically handled in one of two ways. One option is to choose (by hand, or via fitting the output of the calculation to data) a freezeout temperature well below the crossover temperature, and at the surface in spacetime where the fluid in the hydrodynamic calculation cools to this temperature apply the CooperFrye prescription [19] for mapping hydrodynamic volume elements directly onto a phase space distribution of noninteracting hadrons. The second option is to choose a temperature just below the crossover at which to map the hydrodynamic description onto a hadronic transport code, and then let this code describe hadronhadron interactions until the hadrons later freeze out [20]. The second option is an improvement on the first, but it would be even better if the hydrodynamic description itself would tell us when it breaks down. This is impossible within the framework of ideal (zero viscosity) hydrodynamics. An ideal hydrodynamic evolution is fully specified by its initial conditions and by a thermodynamic equation of state giving the pressure in terms of the energy density. Given , an ideal hydrodynamics code will blithely let the expanding fluid evolve until , without ever giving any hint that in reality it has gone far beyond the epoch when hydrodynamics is actually a good approximation. Thus, one must either add a freezeout prescription by hand, or choose by hand to switch from a hydrodynamic to a transport description of the expanding fluid. Once viscous effects are included in the hydrodynamic description, however, it is possible for the hydrodynamic equations to tell us “from within” when they break down.
There are many ways in which a hydrodynamic evolution of an expanding fluid may break down, but we shall focus only on one: cavitation. In an ordinary flowing liquid, cavitation is the formation of bubbles of vapor in regions where the pressure of the liquid drops below its vapor pressure [21]. Cavitation is well studied, both experimentally and theoretically. It can occur on the trailing edges of pump and propeller blades; in this context, the challenge is often to design the blade so as to avoid cavitation since when the bubbles produced by cavitation later collapse they make shock waves that can damage the blades. Cavitation is also used in medicine — one of the treatments of kidney stones involves destroying them via cavitation induced by an ultrasonic pulse. Returning to our context, we shall ask whether the pressure in the expanding droplet of quarkgluon plasma can become negative, which means that it goes below the pressure of the vacuum — our analogue of the vapor phase. This is impossible in ideal hydrodynamics: the thermodynamic is positive. But, in an expanding fluid in the presence of shear and bulk viscosity, shear and bulk stresses make an additional contribution to the spacespace components of the stress energy tensor and these shear and bulk stresses can drive the pressure — which we shall denote — negative. In particular, bulk viscosity is, in essence, a drop in the pressure of an expanding fluid relative to the pressure of an equilibrium fluid at the same energy density, and this decrease in the pressure will play a key role in our considerations. We shall only analyze boostinvariant dimensional expansion of a dimensional fluid [22], meaning that the fluid is only expanding in the direction. We shall ask whether, and when, the shear and bulk stresses are sufficient to drive the longitudinal pressure (which we shall denote ) negative. If goes negative in the hydrodynamic equations that describe boostinvariant expansion, what will happen when the fluid expands to the point where ? At this point, instead of continuing to expand, dilute, and cool in a spatially uniform boostinvariant fashion as before, the fluid will break apart into fragments, separated from each other by vacuum. Vacuum regions are the analogue of the vapor bubbles in conventional cavitation, and their pressure is zero. So, when regions of fluid can stably coexist with regions of vacuuum. The fragments of fluid formed via cavitation will then fly apart, separating from each other, and will subsequently hadronize. Our boost invariant calculation will only allow us to gauge whether and when goes negative: once the fluid cavitates, it is no longer boost invariant, and so our calculation will not be able to describe the subsequent dynamics. Describing the size distribution of the fragments that results from cavitation would first of all require including transverse expansion in the hydrodynamic description, and would second of all require estimating the surface tension associated with the interface between fragments of quarkgluon plasma fluid and vacuum. If this surface tension is large, large fragments will result; if it is small, smaller ones will be favored.
We find that as long as is in the vicinity of the shear stress alone is not enough to trigger cavitation.^{1}^{1}1Well below the crossover temperature , in the hadron gas, rises significantly [23, 24]. If the bulk stress does not do so earlier, the shear stress could trigger cavitation at some temperature well below . However, there is now evidence from a variety of directions [25, 26, 27, 28, 29, 30] that the bulk viscosity is large — — in a narrow range of temperatures around . We shall quantify how large this peak in must be if it is to trigger cavitation when the expanding fluid has cooled to a temperature near . We find that as long as this peak is between 1/4 and 4 times as wide as suggested by current lattice calculations (whose uncertainties we shall discuss), cavitation will occur if the peak is higher than a threshold height that lies between 1/2 and 1/4 that suggested by the lattice calculations.
Our paper is organized as follows. In Section 2 we set up the hydrodynamic equations that describe the boost invariant dimensional expansion of a dimensional fluid, working to second order in derivatives of the velocity field [31, 15, 32, 33, 34, 35, 36, 8, 37, 38]. After setting the full problem up, in the remainder of Section 2 we set the bulk viscosity to zero. We describe how we specify the equation of state (which arises at zeroth order), shear viscosity (first order), and the various coefficients that arise at second order, as well as the initial conditions. We then show results and explore their sensitivity to the second order coefficients and to the shear viscosity. We find much greater sensitivity to the shear viscosity, indicating that, as other authors have found previously, we are using the second order equations in a regime in which the second order effects are much smaller than the first order effects. We turn bulk viscosity on in Section 3. We first describe how we parametrize and the one new second order coefficient that at a minimum must be introduced and in so doing mention some of the uncertainties in our current knowledge of both. We then present our results. In Section 4 we speculate about their implications. The possibility that bulk viscosity could cause the expanding fluid to break apart into fragments has been discussed previously by Torrieri, Tomasik and Mishustin [39] using the formalism of Ref. [40]. We close by sketching several facets of the observed phenomenology of heavy ion collisions that could indicate that freezeout is preceded by cavitation, some previously highlighted in Refs. [39] and some not, some coming from recent analyses of data and some of long standing.^{2}^{2}2Explorations of possible consequences of bulk viscosity in the phenomenology of heavy ion collisions that are not related to cavitation can be found in Refs. [37, 41, 42, 13]. Because we are neglecting transverse expansion throughout we will not be able to make quantitative contact with data. But, our results motivate the importance of including the peak in the bulk viscosity near in hydrodynamic calculations that do include transverse expansion, and the importance of looking for cavitation in these more realistic calculations.
2 Second order relativistic hydrodynamics for a boostinvariant dimensional expansion
2.1 Setup
The energy momentum tensor for relativistic hydrodynamics can be written as
(1) 
where and are the fluid energy density and pressure, is the fluid fourvelocity, normalized such that , is the viscous tensor, satisfying , and where the projector
(2) 
is also orthogonal to . Hydrodynamics is the effective theory describing the longwavelength dynamics of the energy density and the fluid velocity. Its evolution equations describe the conservation of energy and momentum, and are given by
(3) 
where is the geometric covariant derivative. We shall only be interested in hydrodynamics in flat spacetime, but it will be convenient to use curvilinear coordinates to describe boost invariant expansion and we therefore keep the geometric notation. With as in (1), the evolution equations take the form
(4) 
where is the comoving time derivative in the fluid rest frame, is the spatial derivative in the fluid rest frame, and the brackets denote the combination that is symmetric, traceless, and orthogonal to the fluid velocity, namely
(5) 
In general, includes the physics of shear viscosity, bulk viscosity and thermal conductivity. Thermal conductivity is only relevant if there is a nonzero density of some species with a conserved particle number, and we shall work at zero baryon chemical potential throughout. In a conformal fluid, the bulk viscosity vanishes and, including terms up to second order in gradients, satisfies [36, 43]
(6)  
where is the fluid vorticity, where the shear viscosity is the only property of the fluid that enters at first order in gradients, and where , , and are the four properties of the fluid that arise at second order. Obtaining a closed set of evolution equations requires specifying the equation of state and specifying , and in terms of . Equivalently, , , , and can all be specified in terms of the temperature . It is sometimes also convenient to introduce the entropy density
(7) 
Conformality determines the equation of state and implies that , , , and , but conformality alone does not determine any of the dimensionless proportionality constants other than the one in the equation of state.
If we relax the assumption of conformality (while continuing to assume throughout that there is no net baryon density) the only new term that arises on the righthand side of (6) that is firstorder in derivatives is , where is the bulk viscosity. At second order, many new terms arise [44]. As we shall discuss below, it is a standard simplifying assumption to write only the term , where is a new second order coefficient whose role we discuss below.
We shall only consider solutions to the dimensional hydrodynamic equations in which no quantity depends on the transverse spatial directions and (which in particular means no vorticity) and in which the expansion in the direction is boost invariant [22]. This makes it convenient to change variables from to where is the proper time and is the spacetime rapidity. These curvilinear coordinates are comoving with the fluid, meaning that and the spatial components of all vanish. And, in these coordinates boost invariance implies significant further simplifications: is diagonal, and therefore so is , and the diagonal components of depend only on , not on . Upon making these simplifications, the stress energy tensor in coordinates takes the form [31, 15, 32, 33, 34, 35, 36, 8, 37, 38]
(8) 
where the trace of — namely — and the traceless part of — namely — denote the nonequilibrium contributions to the pressure coming from the bulk and shear stresses, respectively.^{3}^{3}3The quantity that we denote as has been called in some of the literature and elsewhere in the literature. For a fluid at rest, the pressure is isotropic and is given by , which is related to the energy density by the thermodynamic equation of state . As the fluid is expanding, unless it is an ideal fluid with its pressure is no longer isotropic — the transverse and longitudinal pressure are given by
(9)  
(10) 
Furthermore, upon making these simplifications the second order evolution equations are
(11)  
(12)  
(13) 
At first order, and are given by their NavierStokes values
(14) 
and
(15) 
We see that if we ignore the terms in square brackets in (12), then the second order equations (12) and (13) describe and relaxing towards their NavierStokes behavior (14) and (15) with time constants and , along the lines of the IsraelStewart approach to second order dissipative relativistic hydrodynamics [45]. If we ignore bulk viscosity, setting , then (11) and (12), including in particular the terms in the square brackets in (12), follow from conformality [36, 43]. However, once we turn on bulk viscosity we are breaking conformality, and there can then be further second order terms in both (12) and (13) [44]. These equations are in this sense provisional, but it should be noted that the terms in square brackets in (12) become neglible at large and we expect the same to be the case for the missing nonconformal terms also.
2.2 Signs of cavitation
Since and at first order, see (14) and (15), it is reasonable to expect them to have these signs in solutions to the second order equations also. We then see from (9) and (10) that if either or is large enough, the longitudinal pressure can be driven negative, and if is large enough, the transverse pressure can also be driven negative. We shall see in Section 3 that if rises high enough at temperatures in the vicinity of the crossover from quarkgluon plasma to hadron gas, the resulting bulk stress does drive negative, indicating cavitation.
2.3 Equation of state, shear viscosity, and
We see that in order for the evolution equations (11), (12) and (13) to be closed we need the equation of state and expressions relating , , , and to . For the remainder of Section 2, we shall set and therefore , deferring our analysis of the effects of bulk viscosity to Section 3. We then need , , and only.
We need an equation of state that describes the quarkgluon plasma phase as well as the crossover to a hadron gas. Lattice quantum field theory is wellsuited to the calculation of thermodynamic quantities at zero baryon chemical potential, and so there are many lattice calculations of in QCD that we could employ. We shall take from Ref. [46], both because it is an example of the state of the art and because these authors have provided a parametrization of their results that is easy to use. They parametrize their results for the trace anomaly using the functional form
(16) 
and give values with error bars for the coefficients , , and for calculations done with two different lattice actions, with and without combining these calculations with hadron resonance gas calculations valid at lower temperatures. We shall use the central values of their results obtained from combining lattice calculations done with the p4 action and hadron resonance gas calculations: GeV, GeV, GeV, and GeV. These authors find a crossover between hadron gas and quarkgluon plasma occurring in a temperature regime MeV MeV. In Section 3 when we need to specify a value of the crossover temperature we shall use MeV, in order to be consistent with the equation of state that we employ throughout. The pressure is related to the trace anomaly by
(17) 
and the results of Ref. [46] are obtained by choosing MeV and . We shall only work at MeV, where there is no effect of these choices. Knowing and as functions of , we know as a function of also, as well as the entropy density (7). And, from and we have the equation of state . We shall use the same equation of state throughout this paper, focussing on the effects of varying other quantities.
Next, we turn to the shear viscosity . We shall use
(18) 
as a baseline value, and we shall explore the effects of varying . The relationship (18) holds for the plasma phase of any gauge theory that has a dual gravity description, in the limit of large numbers of colors and infinitely strong coupling [16, 17, 18]. Even though much larger values of (of order 1 and larger) are expected both in the hadron gas found well below [23, 24] and in the weakly coupled quarkgluon plasma found far above [47], the baseline (18) is seen as a reasonable starting point for the analysis of quarkgluon plasma in the regime being explored by RHIC collisions, say around . Lattice QCD calculations, to date in a gluon plasma without quarks, indicate values of only a few times (18) [48, 30]: at Meyer finds with a statistical error of [30]. (At this temperature, is small compared to .) Comparison between second order relativistic viscous hydrodynamic calculations that include transverse expansion [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] and data from RHIC [1, 2] on the azimuthal anisotropy of ultrarelativistic heavy ion collisions that have a significant impact parameter indicate that in these collisions approximate local thermal equilibrium is attained rapidly and is small, apparently and conservatively [15, 14]. Given that we analyze longitudinal expansion only, we can have nothing to say about the extraction of information about from these data. But, we shall confirm that dimensional boost invariant expansion is modified significantly as we vary between 0 and .
Less is known about the values of the secondorder coefficients and . We shall take as a baseline
(19) 
and
(20) 
which are their values in the plasma of supersymmetric YangMills theory [36, 43, 49], the simplest and best studied example of a strongly coupled plasma with a dual gravity description. Meyer finds at in lattice calculations of QCD without quarks [30], within a factor of a few of (19).^{4}^{4}4In kinetic theory, [50, 35, 51], with the prefactor depending on the value of the coupling constant. Kinetic theory is not quantitatively valid if , but applying it anyway gives in agreement with Meyer’s lattice result. In kinetic theory, [50], which with would give a value of within a factor of two of (20). We shall find that the effects of varying and by large factors are small, indicating that the hydrodynamic calculations are being done in a regime in which second order effects are small compared to first order effects.
2.4 Baseline results
In Fig. 1 we show an example of a solution to the evolution equations (11) and (12) with vanishing bulk viscosity. The pressure and temperature are related to the energy density via the lattice calculations of QCD thermodynamics from Ref. [46] that we have described in Section 2.3. The shear viscosity and the secondorder coefficients and have been set to their baseline values (18), (19) and (20).
We have initialized the evolution in Fig. 1 at a time fm, a reasonable choice given that the RHIC data on the anisotropic expansion in heavy ion collisions at nonzero impact parameter can only be understood if a hydrodynamic description is already relevant earlier than 1 fm after the collision [52]. We have also chosen a value for the initial energy density that is reasonable for collisions at the top RHIC energy. At time fm in the evolution of Fig. 1, the energy density is GeVfm, which is consistent with estimates of the energy density at this time that have been made using data on the final state energy and entropy [1, 53]. (Using the lattice results for , this energy density corresponds to a temperature MeV.) So, although it is misplaced precision to specify initial conditions with GeVfm (and MeV) at fm, we shall make this choice throughout this paper as in Fig. 1 since varying these choices within a reasonable range would have no qualitative effects. Note also that if we were doing phenomenology (which we cannot do given our dimensional expansion) we would want to adjust the initial state energy density as we vary other parameters (which we will do below) in order to maintain the same late time energy density. We shall not do this, since our purpose is to explore how solutions to the evolution equations depend on parameters, and tweaking the initial conditions as we varied the parameters would for this purpose be a complication.
In Fig. 1, we have chosen at fm. There is no phenomenological justification for this choice. Instead, we find that this choice does not matter. What we observe in the evolution is that rapidly (over a timescale that is a few tenths of a fm in Fig. 1 and that is controlled by ) increases until it is close to its NavierStokes behavior (14), and then follows (14) closely during the subsequent evolution. If instead of initializing with we choose at fm to be twice its NavierStokes value, we find the same behavior. And, varying the initial value of over this range makes very little difference — it changes by less than half as much as we shall find when we vary in the next section.
2.5 Insensitivity to and
In Fig. 2 we see that both the second order coefficients and can be increased by large factors without having significant effects on the evolution of . Increasing extends the initial period of time in the evolution when the shear stress changes from its initial value 0 to approach its NavierStokes behavior (14). With increased relative to its baseline value (19) by a factor of 8, these earlytime transients last fm. Even in this case, Fig. 2 shows that the modification of is small. Increasing depresses relative to its NavierStokes behavior (14) at all times, but the effect is small even when is increased relative to its baseline value (20) by a factor of 16, and the effect on is even smaller. To see significant effects on , must be increased by factors of , and even then the effect on is only of order 10 percent. We conclude that we are using the evolution equations in a regime in which the effects of the second order terms are small — smaller, we shall see, than the effects of the first order terms. We shall therefore use the baseline values (19) and (20) exclusively in results that we shall quote throughout the following. But, we have checked that varying these parameters simultaneously with the variations of and that we discuss below does not change any interesting conclusions.
2.6 Sensitivity to shear viscosity
In Fig. 3 we see that changing the coefficient of the first order term in the evolution equations, namely the shear viscosity , has much more significant effects than those we found in Section 2.5. Increasing by a factor of two has a % effect on the energy density . The fact that the evolution equations are much more sensitive to variation of than they are to variations of or provides qualitative support to the program [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] of using comparison between hydrodynamic calculations that include anisotropic transverse expansion and data from RHIC to extract information about the value of — this extraction should not be complicated by our lack of knowledge of the values of or . We will revisit this conclusion in Section 4 after considering the effects of bulk viscosity, which we have so far neglected, in Section 3.
Furthermore, we see from the right panel of Fig. 3 that, as Martinez and Strickland have analyzed in detail [38], increasing makes the longitudinal pressure of (10) negative at early times. With the initial conditions that we are using, we find that is negative for some window of early ’s if . Martinez and Strickland have analyzed how this criterion depends on the choice of initial conditions. It is easy to see why must arise at early times for sufficiently large . We have seen that the shear stress rapidly rises to its NavierStokes value (14). Together with this means that, after transient behavior at very early times, . The zeroth order (ideal hydrodynamics; no shear viscosity) solution to the evolution equations for boostinvariant expansion with an equation of state has and , meaning that if then grows faster for than does. This means that at some early time, must exceed making . Our results and the results of Ref. [38] confirm that the conclusions of this simple argument apply. With the specific initial conditions that we have used — namely with initially — negativity of can be avoided at any given by increasing by a large enough factor. This stretches out the initial transient, delaying the at which is reached until late enough that never exceeds . However, this resolution is specific to our (completely arbitrary) choice of the initial value of — we could have chosen from the beginning — and is therefore not actually a resolution. The conclusion we should draw is simply that the hydrodynamic description, premised on local thermal equilibrium, must break down at early times and the negativity of is one sign that tells us before when we cannot go. (Other evidence for the same qualitative and, essentially, quantitative conclusion has been developed in Refs. [35, 51].) For , initializing the hydrodynamic evolution equations at fm as we are doing is appropriate. For , choosing fm is only appropriate for certain initial conditions (including ), while choosing fm is safe for more generic initial values of . For , we find that in order to avoid we must choose fm. It would be interesting to pursue this line of reasoning further in the higherdimensional hydrodynamic calculations of Refs. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]: the same hydrodynamic calculations that yield information about the allowed values of via comparison to data from RHIC should at the same time constrain the earliest time at which a hydrodynamic description can be valid, via checking before what time these hydrodynamic evolutions feature negative pressure in some region of space.
The shear viscosity to entropy ratio becomes large at late times, in the hadron gas phase [23, 24]. In our calculations with and , the shear stress is much less than at late times. But, if we consider increasing at late times to , as appropriate at MeV according to the calculations of Demir and Bass [24], we find coming close to going negative at the very late times corresponding to MeV. This would be worth further investigation, as a possible indicator of when the hydrodynamic description breaks down at late times, if not for the fact that freezeout in heavy ion collisions is expected to occur earlier than this. And, furthermore, we shall see that including the effects of bulk viscosity can result in the breakdown of hydrodynamics also happening earlier, when . To this we now turn.
3 Effects of bulk viscosity
We now wish to turn on a nonzero bulk viscosity and study its effects on solutions of the evolution equations (11), (12) and, now, (13). We shall set , and to their baseline values (18), (19) and (20) throughout this Section. In the case of and we do so because, as in Section 2.5, we do not expect that changing their values would have significant consequences. But, we have seen that varying is consequential. We are setting in order to be conservative, in the following sense. We shall focus on the question of whether the longitudinal pressure of (10) goes negative. Increasing increases , which makes a negative contribution to . So, if we find that the bulk viscosity drives negative with set to , increasing would only make even more negative.
Both at low temperatures in the hadron gas [54, 55, 56, 57] and at very high temperatures where the quarkgluon plasma is weakly coupled [58] the bulk viscosity is much smaller than the shear viscosity. These calculations indicate that rises as one approaches from both below and above. And, if the crossover at were a second order phase transition, would peak at [25, 28, 29]. The general expectation that may be significant near is supported by the analysis of Refs. [26, 28] which uses sum rules to relate the bulk viscosity to (derivatives of) thermodynamic quantities calculated on the lattice (although this relation is subtle [59, 60, 61]) and by the analyses of Ref. [27] in which itself is constrained via lattice calculations, albeit in QCD without quarks. Both these approaches find in a narrow range of temperatures very near . We shall ask whether such a peak in can drive negative, triggering cavitation.
3.1 Choosing and
Determining via lattice calculations of Euclidean correlation functions is challenging, and the results obtained in Ref. [27] should not yet be seen as definitive [59, 62, 63]. To give just one example of a difficulty [59], just as peaks at a second order phase transition, so does the relaxation time — due to the phenomenon of critical slowing down. The Euclidean lattice calculations are sensitive to the ratio , making it hard to disentangle one from the other. Note, however, that attributing a peak in to is conservative in the sense that by neglecting the rise in one underestimates the rise in . Unfortunately, other more technical challenges can go in the other direction [62, 63]. Progress is nevertheless possible [64, 30]. The current state of affairs is that lattice calculations make a robust case for the existence of a peak in at , at least in QCD without quarks [30], but they do not yet reliably determine the height of the peak. We shall therefore parametrize the results of Ref. [27], and vary the parameters over a considerable range.
In Ref. [27], Meyer reports results for at five values of the temperature, 1.02, 1.24, 1.65, 2.22 and 3.22. One way to parametrize his results is to write
(21) 
with , and the parameters. ( is not a parameter in that QCD without quarks has a first order phase transition at a reliably determined . When we employ (21) together with the equation of state for QCD with quarks specified via (16), we shall set MeV.) Meyer’s results at the three higher temperatures, well above , are consistent with . This is not surprising since the trace anomaly , which like is a dimensionless measure of the breaking of conformal invariance, is at high temperature. (For example, see (16).) If we set and fit to Meyer’s central values of at the three higher temperatures, we find
(22) 
Note that at these higher temperatures, Meyer gives the central values for that we have used but his results remain consistent with . When we vary , therefore, we should consider values as small as zero. With chosen as in (22) and with , the curve (21) is far below Meyer’s results at and . That is, there is no way to use simply to fit Meyer’s high temperature results and his results at and — where has its peak. The parameters and can then be chosen such that (21) passes directly through Meyer’s central values of at these two temperatures, yielding
(23) 
The parameter controls the height of the peak in and controls its width, and we shall vary both considerably.
Little is known about the value of . We shall use
(24) 
with given in (19) as the baseline value for . Although there is no strong argument for this choice, it holds in one class of strongly coupled nonconformal fluids [65, 44].^{5}^{5}5In kinetic theory, [51]. On general grounds (i.e. as a manifestation of critical slowing down) and as in the specific strongly coupled nonconformal fluid studied in Ref. [66], is expected to peak where peaks, and we shall therefore check the effects of greater than in (24) by as much as a factor of 40.
3.2 Boost invariant expansion with bulk viscosity
Once one has picked and , the next step is to choose initial conditions. We shall initialize at fm and choose as in Section 2. We shall choose the initial shear stress as in Section 2, and we shall also choose the bulk stress . When we then evolve the equations of motion (11), (12) and (13), we find that quickly evolves to its NavierStokes value (15) during an initial time that is controlled by . We shall always stop the evolution at the time when has dropped to , since our parametrization of in (21) is only valid for . We expect to drop rapidly below , but even less is known about the shape of the peak below than above it, so we simply stop the evolution when and ask whether by that time the longitudinal pressure has gone negative.
In describing our results, we begin with . That is, we begin with no peak in , just with . If we choose as in (22), suitable to describe the lattice results at [27], we find that the bulk viscosity has negligible effects. Introducing the bulk viscosity changes the energy density by about 6%. And, with as in (22), we find that we must increase the shear viscosity to in order to see at early times — whereas we saw in Section 2.6 that with no bulk viscosity this required . If we reduce relative to (22), the effects of the bulk viscosity become even more negligible. Even if we increase by a factor of 2 relative to (22), we still find no qualitative consequences: the energy density changes by about 13% and at early times requires . Note that with greater than (22) by a factor of two there is already a range of temperatures near where . Increasing by another factor of two makes over a wide range of temperatures, which is not supported by the lattice calculations. Henceforth, we fix as in (22), meaning that if were zero the bulk viscosity would not have interesting consequences.
We now investigate the consequences of the peak in near . Let us begin by choosing and as in (23), meaning that the peak in has a height and width as indicated by Meyer’s lattice results [27]. We illustrate the resulting evolution in the topleft panel of Fig. 4. We see that with and as in (23), the rising bulk viscosity drives the longitudinal pressure negative when is still well above . As we have described in Section 3.1, although there is good evidence for a peak in near , its height and width (which we are parametrizing by and ) are not well known. So, we have explored for what values of these parameters is driven negative near . As the topright panel of Fig. 4 shows, upon reducing by a factor of 0.5 while keeping fixed we continue to find . In fact, we find that with as in (23), the largest at which remains positive for all , which we shall denote , is , which is 41% of the value (23) indicated by the lattice calculations of Ref. [27].^{6}^{6}6 We have used the equation of state (16) and (17) obtained from lattice QCD calculations throughout. If instead we use the conformal equation of state with no phase transition, then . It is easy to see why a larger peak in the bulk viscosity is then required in order to drive the longitudinal pressure negative: in the vicinity of the peak in the bulk viscosity, is significantly larger than in (17) because (17) describes a phase transition; because is larger, it takes more bulk stress to drive negative. One measure of the robustness of our results is that even with the much larger thermodynamic pressure , a peak in the bulk viscosity comparable to that indicated by the lattice calculations of Ref. [27] is sufficient to trigger cavitation. How does changing the width modify this result? If we increase by a factor of two relative to (23), drops only slightly, to . If we increase it by another factor of two, making the peak in four times wider than in Ref. [27], drops to 0.22. If we change in the other direction, decreasing it by a factor of two relative to (23), meaning that we make the peak in twice as narrow as in Ref. [27], rises only very slightly, to . If we decrease by another factor of two, rises to ; if we decrease it by yet another factor of two — making the peak in eight times narrower than in Ref. [27] — increases to 0.54. In the bottomleft panel of Fig. 4 we illustrate the case with reduced relative to (23) by a factor of four, and unmodified, as in (23). Comparing this panel to the topleft panel, we see that reducing the width of the peak in the bulk viscosity by a factor of four delays the time at which goes negative, but does not significantly reduce the amount by which it goes negative. This is consistent with the observation that is not much changed. We conclude that is quite insensitive to over a wide range of , ranging from 0.22 if is four times larger than (23) to 0.45 if is four times smaller than (23). Over this wide range, is well below the value 0.90 indicated by the lattice calculations of Ref. [27].
Note that all the results we have just quoted were obtained with . If we increase , it will take a smaller to push the longitudinal pressure negative. For example, with as in (23) if we increase to this decreases the value of from 0.37 to 0.32.
It is also worth checking that the fact that is being driven negative really is due to the peak in the bulk viscosity, not to the component of in (21). To this end, we set (with and as in (23)) and found that eliminating the component of increased , but only from 0.37 to 0.41.
The curves in all the panels in Fig. 4 end at the when in the calculation, although they have ceased to be relevant earlier when reaches zero and cavitation occurs. Note, however, that the at which is between 4.7 fm and 5.2 fm in all four panels, while in Fig. 1 reaches at 4.3 fm. This is a small effect, but it can be made larger. As Kapusta discovered in Ref. [67], if diverges at , namely if for some positive power , then the diverging bulk viscosity acts in the hydrodynamic equations so as to prevent from dropping below : as the equations show approaching from above but never reaching it. The slight delay in the at which is reached in our calculations, as in Fig. 4, is the small residue of this effect when the peak in is finite. Note that the solution with diverging and approaching from above is academic, since the solution has negative. In fact, with a powerlaw divergent the bulk shear is large enough to drive already at rather early times. Cavitation occurs long before the asymptotic solution discovered in Ref. [67] is reached.
Following Ref. [37], we have investigated the entropy produced according to the hydrodynamic equations that we have solved. For example, comparing the calculation illustrated in the topleft panel of Fig. 4, in which is as in (21) with (22) and (23), with the calculation illustrated in Fig. 1, in which , we find that turning on a bulk viscosity with a peak as in Meyer’s lattice results [27] increases the entropy by about 20%. This suggests that entropy production due to the bulk viscosity is a modest effect, as was found previously in Ref. [37] for cases in which the peak in is not high enough to cause cavitation. In our case, however, this result cannot be trusted because there is no way to use our boostinvariant calculation to estimate how much further entropy is produced upon cavitation.
We now consider the effects of increasing , which is expected to be large where peaks. We will not attempt to model a timevarying ; instead, we ask what are the consequences of increasing relative to (24) at all temperatures. If we increase by a factor of 10, keeping as in (23), we find that decreases from 0.37 to 0.33. If we increase by a further factor of two, meaning that it is 20 times greater than in (24), increases to 0.49. If we increase by one further factor of two, making it 40 times greater than in (24), increases to 0.94, comparable to the in (23). In the bottomright panel of Fig. 4 we illustrate the case with . The bulk stress changes more gradually as a function of time, and as a result does not go as far negative as in the topleft panel of Fig. 4, but the change is not dramatic and as a consequence is greater, but not by much.^{7}^{7}7After the first version of this paper appeared, Song and Heinz found that cavitation does not occur anywhere in their dimensional calculation if they use [80], which is larger than (24) by a factor of several hundred in the vicinity of the peak in the bulk viscosity. Our results are consistent with this. This insensitivity to large changes in the value of is one indication that second order effects are still small at the time when cavitation occurs. Another indication of this is the fact that is small, for example less than 10% at the time when cavitation occurs in the topleft panel of Fig. 4.^{8}^{8}8As noted in Section 2.1, once we break conformality by introducing a nonzero bulk viscosity further second order terms can arise in both (12) and (13) [44]. See Ref. [69] for one example. We have confirmed that adding the second order terms considered in Ref. [69] has only negligible effects on our results.
We can summarize these results as follows:

With and as in (23), the peak in the bulk viscosity above drives the hydrodynamic evolution to negative , indicating cavitation.

Stable hydrodynamic evolution all the way down to requires that the peak in near be less than a threshold that is one quarter to one half as high as the peak found in Ref. [27]; the threshold peak height is fairly insensitive to the width of the peak, for widths between one quarter and four times that found in Ref. [27].

The effects of a peak in near can be washed out by increasing , the relaxation time for the bulk stress . However, the increase must be by a very large factor. If is larger than of (19) by a factor of 10 (or 20), results are little affected and the peak in must still be reduced by about a factor of three (or two) relative to (23) in order to obtain stable hydrodynamic evolution all the way down to . If is 40 times greater than , however, the effects of a bulk viscosity peak are washed out sufficiently that remains just barely positive even for a peak whose height and width are as in (23).
We have focused on of (10) rather than of (9) because the nonzero shear stress implies that as the bulk stress becomes increasingly negative, goes negative first, only later. By the time goes negative in the calculation, the calculation has already broken down at the time when went negative, triggering cavitation. But, we see in all the panels in Fig. 4 that is quite small by the time goes negative, meaning that is already close to zero when reaches zero. It will be interesting to see, therefore, which component of the pressure goes negative first at which location in the fluid in a calculation that includes transverse expansion.
4 Implications
In thinking through the implications of our results, the first possibility to consider is that the peak in near in QCD with quarks is in fact not so high that cavitation results. As we discussed in Section 3.1, the current lattice calculations of the height of the bulk viscosity peak in QCD without quarks have various caveats, meaning that the peak in this theory could be smaller (although it could just as well be larger) than is indicated by the results of Ref. [27]. Furthermore, there are various indications that is somewhat lower in QCD with quarks than in QCD without quarks. At very high temperatures where the quarkgluon plasma is weakly coupled, if one compares the two theories at a fixed small value of the QCD coupling, say , one finds that in QCD with three flavors of quarks is about 56% of that in quarkless QCD [58]. At these high temperatures, is smaller than in value, so this comparison can give only rough guidance to how the height of the peak near will change when quarks are introduced, but it does suggest that it will decrease. A second argument is simply that the transition in QCD with quarks is a crossover whereas that in quarkless QCD is first order, and if adding quarks smooths out the transition then it is reasonable to guess that it will also round off the peak in . The magnitude of this effect can be guessed by looking at which, like , is a dimensionless measure of the breaking of conformality. In quarkless QCD, peaks at a value of 0.85 [68] while in the equation of state (16) for QCD with quarks, peaks at 0.53. So, both this argument and the comparison to what happens at very high temperatures suggest that the peak in in QCD with quarks is (very roughly) about half as high as in (21) with as in (23). Taking our results at face value, this would put it just above , meaning that the peak in would trigger cavitation very close to . The uncertainties are large and it could certainly be that the peak height ends up lower than , and no cavitation occurs near . The previous studies of the effects of the peak in the bulk viscosity in Refs. [37, 11, 13] have explored the consequences of peaks that are not high enough to cause cavitation.^{9}^{9}9It is worth noting that in the examples of nonconformal plasmas in which the authors of Refs. [70, 65, 71, 72] have been able to compute via gauge/gravity duality, a peak in is seen but it is not large enough to cause cavitation. For example, in the model of Ref. [65] the ratio is given by , with the speed of sound, and is therefore everywhere less than . In contrast, in the example analyzed in Ref. [73] via gauge/gravity duality can be and rises comparably high to the peak found in Meyer’s lattice calculations [27], more than high enough to trigger cavitation. At present, these calculations taken together therefore support the existence of a peak in but do not provide sufficient guidance regarding its height. If this is the path that nature chooses, then hydrodynamic calculations can be followed down to temperatures below , and it becomes interesting to investigate the possibility of cavitation at a lower temperature, driven by the rising shear viscosity of the hadronic phase at low temperatures [23, 24].
It is more interesting to consider the possibility that the peak in near in QCD with quarks is large enough to cause the expanding fluid produced in heavy ion collisions to cavitate when it cools through , as our dimensional calculations indicate. There are a variety of aspects of the observed phenomenology of heavy ion collisions that give some support to this possibility:

The possibility that the peak in the bulk viscosity near can cause the fluid to fragment, and then freezeout, has been considered previously in Refs. [39]. These authors suggest that data on two particle momentum correlations (the HBT effect) in heavy ion collisions at RHIC can be understood if the hadrons in the final state come from such fragments. They have also suggested other experimental observables in Ref. [74].

One of the workhorses of heavy ion phenomenology is the statistical hadronization model, reviewed in Ref. [75], which has been used successfully to describe the ratios among the yields of many different hadrons using a few parameters including the chemical freezeout temperature and baryon chemical potential. One of the conceptual underpinnings of this model, described in Ref. [75] and going back to the original formulation of Hagedorn [76], is the assumption that high energy collisions give rise to multiple clusters — colorless, extended, massive objects — which then hadronize statistically (meaning that all hadronic final states consistent with conservation laws are equally likely). Most work in this context has focused on the statistical hadronization; the dynamics that results in the generation of the clusters in the first place has received less attention. This dynamics is no doubt complex, and may be quite different in hadronhadron and heavy ion collisions. Our work suggests that in the case of ultrarelativistic heavy ion collisions, in which a hydrodynamic description in terms of an expanding fluid with is appropriate in the early stages of the collision, the clusters required by the statistical hadronization model may arise via cavitation, and this cavitation may be triggered by the peak in in the vicinity of . Perhaps this could be an explanation of why the chemical freezeout temperatures extracted via the use of the statistical hadronization model seem to be in the vicinity of [1].

In Refs. [77], Broniowski et al have explained the eventbyevent fluctuations in the mean transverse momentum of the particles produced in heavy ion collisions at RHIC via the same assumption that underlies the statistical hadronization model, namely that hadronization is preceded by the material produced in a heavy ion collision falling apart into clusters, each of which then yield 6 to 15 hadrons when they hadronize.

In Ref. [78], the PHOBOS collaboration provides evidence (from twoparticle correlations in pseudorapidity and azimuthal angle) that the hadrons in the final state produced in heavy ion collisions at RHIC come from clusters that decay into 3 to 6 charged hadrons, meaning 5 to 9 hadrons in all.
As we discussed in Section 1, it is pleasing to have a means by which a hydrodynamic calculation can predict its own break down. The peak in the bulk viscosity near can provide a simple and elegant means: if this peak is high enough — as we have quantified — it drives the longitudinal pressure to zero at which point the fluid cavitates. The phenomenological evidence in support of the notion that hadronization in ultrarelativistic heavy ion collisions is preceded by cavitation, with the fluid fragmenting into droplets that play the role of the clusters which have long been employed in the statistical hadronization framework, is perhaps not yet overwhelming. But, together with our investigation, it certainly seems sufficient to take this possibility seriously.
There are many avenues open for further investigation:

An investigation of the effects of other possible second order terms that can arise in the hydrodynamic equations for a nonconformal fluid, as well as third order terms, would be desirable. Little is known about the coefficients of such terms. But, although they will have quantitative effects, given the insensitivity of our results to large variations in , and there is no reason to expect qualitative effects.

A more interesting direction to pursue is to repeat our study using a dimensional hydrodynamic code, describing both longitudinal and transverse expansion. This will make the determination of the height of the peak in that is needed in order to trigger cavitation when the fluid cools through more quantitative. And, having a fluid whose energy density varies with tranverse position will raise new questions and open new possibilities. For example, cavitation should occur much earlier at the (cooler) edges of the collision region, since there earlier [80, 79]. Cavitation should start at the edges and move inward, just as hadronization has long been understood to do.

It is important to investigate whether if freezeout and hadronization are triggered by cavitation near this modifies the extraction of via comparison between data and hydrodynamic calculations of the anisotropic expansion of the fluid produced in collisions with a nonzero impact parameter. The effect of the physics of cavitation near on this comparison may prove minimal, since the anisotropic flow is generated early in the collision, when the hot fluid is still azimuthally anisotropic in shape. This means that the anisotropic flow is generated well before cavitation is triggered, when the bulk viscosity is still small compared to the shear viscosity.

From the point of view of understanding the observable, and perhaps observed, phenomenology of cavitation, the most important question is the determination of the size distribution of the droplets formed when the fluid cavitates. Work in this direction can be found in Refs. [39]. A from first principles determination of the size distribution will be challenging: for one, it will require determining the surface tension associated with the interface between expanding quarkgluon plasma with and vacuum. If this surface tension is small, small droplets will be favored. The requirement that the droplets must be color singlets will require rearrangement of color within of the surface as droplets separate during cavitation. Although hard to quantify, this can be thought of as a contribution to the surface tension which sets a limit on the smallness of the droplets that is of order a few times . It is then interesting to note that a spherical droplet with a radius of 1 fm that has the energy density obtained from the equation of state (16) at contains about 6 GeV of energy, which is in the right ballpark to explain the PHOBOS data [78] mentioned above. This suggests that the surface tension is small enough that cavitation yields many small droplets. The picture to have in mind is that as the quarkgluon plasma produced in an ultrarelativistic heavy ion collision expands and cools through , the fluid falls apart into a mist of hundreds of small droplets, each of which later hadronizes as in the statistical hadronization model. Cavitation into a mist of small droplets which then become a hadron gas is not likely to have dramatic observable consequences
Acknowledgments
We acknowledge very helpful conversations with Ulrich Heinz, Harvey Meyer, Tomoi Koide, Gunther Roland, Paul Romatschke, Huichao Song, Misha Stephanov and Derek Teaney. NT is grateful to the Research Science Institute of the Center for Excellence in Education for supporting his research. This research was supported in part by the DOE Office of Nuclear Physics under contract #DEFG0294ER40818.
References
 [1] K. Adcox et al. [PHENIX Collaboration], Nucl. Phys. A 757, 184 (2005) [arXiv:nuclex/0410003]; B. B. Back et al. [PHOBOS Collaboration], Nucl. Phys. A 757, 28 (2005) [arXiv:nuclex/0410022]; I. Arsene et al. [BRAHMS Collaboration]; Nucl. Phys. A 757, 1 (2005) [arXiv:nuclex/0410020]; J. Adams et al. [STAR Collaboration], Nucl. Phys. A 757, 102 (2005) [arXiv:nuclex/0501009];
 [2] B. Alver et al. [PHOBOS Collaboration], Phys. Rev. Lett. 98, 242302 (2007) [arXiv:nuclex/0610037]; B. I. Abelev et al. [STAR Collaboration], Phys. Rev. C 77, 054901 (2008) [arXiv:0801.3466 [nuclex]]; J. Y. Ollitrault, A. M. Poskanzer and S. A. Voloshin, arXiv:0904.2315 [nuclex]; P. Sorensen, arXiv:0905.0174 [nuclex].
 [3] R. Baier and P. Romatschke, Eur. Phys. J. C 51, 677 (2007) [arXiv:nuclth/0610108].
 [4] P. Romatschke and U. Romatschke, Phys. Rev. Lett. 99, 172301 (2007) [arXiv:0706.1522 [nuclth]].
 [5] H. Song and U. W. Heinz, Phys. Lett. B 658, 279 (2008) [arXiv:0709.0742 [nuclth]].
 [6] K. Dusling and D. Teaney, Phys. Rev. C 77, 034905 (2008) [arXiv:0710.5932 [nuclth]].
 [7] H. Song and U. W. Heinz, Phys. Rev. C 77, 064901 (2008) [arXiv:0712.3715 [nuclth]].
 [8] M. Luzum and P. Romatschke, Phys. Rev. C 78, 034915 (2008) [Erratumibid. C 79, 039903 (2009)] [arXiv:0804.4015 [nuclth]].
 [9] H. Song and U. W. Heinz, Phys. Rev. C 78, 024902 (2008) [arXiv:0805.1756 [nuclth]].
 [10] D. Molnar and P. Huovinen, J. Phys. G 35, 104125 (2008) [arXiv:0806.1367 [nuclth]].
 [11] H. Song and U. W. Heinz, arXiv:0812.4274 [nuclth].
 [12] M. Luzum and P. Romatschke, arXiv:0901.4588 [nuclth].
 [13] H. Song and U. W. Heinz, arXiv:0907.2262 [nuclth].
 [14] For reviews, see U. W. Heinz, arXiv:0901.4355 [nuclth], P. Romatschke, arXiv:0902.3663 [hepph] and D. A. Teaney, arXiv:0905.2433 [nuclth].
 [15] D. Teaney, Phys. Rev. C 68, 034913 (2003) [arXiv:nuclth/0301099].
 [16] G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 87, 081601 (2001) [arXiv:hepth/0104066].
 [17] P. Kovtun, D. T. Son and A. O. Starinets, JHEP 0310, 064 (2003) [arXiv:hepth/0309213].
 [18] A. Buchel and J. T. Liu, Phys. Rev. Lett. 93, 090602 (2004) [arXiv:hepth/0311175].
 [19] F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
 [20] For example, S. A. Bass and A. Dumitru, Phys. Rev. C 61, 064909 (2000) [arXiv:nuclth/0001033] and D. Teaney, J. Lauret and E. V. Shuryak, arXiv:nuclth/0110037.
 [21] For example, see C. E. Brennen, Cavitation and Bubble Dynamics, Oxford, 1995; and http://en.wikipedia.org/wiki/Cavitation.
 [22] J. D. Bjorken, Phys. Rev. D 27, 140 (1983).
 [23] M. Prakash, M. Prakash, R. Venugopalan and G. M. Welke, Phys. Rev. Lett. 70, 1228 (1993) [Nucl. Phys. A 566, 403C (1994)]; M. Prakash, M. Prakash, R. Venugopalan and G. Welke, Phys. Rept. 227, 321 (1993); D. Davesne, Phys. Rev. C 53, 3069 (1996); A. Dobado and F. J. LlanesEstrada, Phys. Rev. D 69, 116004 (2004) [arXiv:hepph/0309324]; J. W. Chen and E. Nakano, Phys. Lett. B 647, 371 (2007) [arXiv:hepph/0604138].
 [24] N. Demir and S. A. Bass, Phys. Rev. Lett. 102, 172302 (2009) [arXiv:0812.2422 [nuclth]]; N. Demir and S. A. Bass, arXiv:0907.4333 [nuclth].
 [25] K. Paech and S. Pratt, Phys. Rev. C 74, 014901 (2006) [arXiv:nuclth/0604008].
 [26] D. Kharzeev and K. Tuchin, JHEP 0809, 093 (2008) [arXiv:0705.4280 [hepph]].
 [27] H. B. Meyer, Phys. Rev. Lett. 100, 162001 (2008) [arXiv:0710.3717 [heplat]].
 [28] F. Karsch, D. Kharzeev and K. Tuchin, Phys. Lett. B 663, 217 (2008) [arXiv:0711.0914 [hepph]].
 [29] C. Sasaki and K. Redlich, Phys. Rev. C 79, 055207 (2009) [arXiv:0806.4745 [hepph]].
 [30] H. B. Meyer, arXiv:0907.4095 [heplat].
 [31] A. Muronga, Phys. Rev. Lett. 88, 062302 (2002) [Erratumibid. 89, 159901 (2002)] [arXiv:nuclth/0104064].
 [32] A. Muronga, Phys. Rev. C 69, 034903 (2004) [arXiv:nuclth/0309055].
 [33] A. Muronga and D. H. Rischke, arXiv:nuclth/0407114.
 [34] U. W. Heinz, H. Song and A. K. Chaudhuri, Phys. Rev. C 73, 034904 (2006) [arXiv:nuclth/0510014].
 [35] R. Baier, P. Romatschke and U. A. Wiedemann, Phys. Rev. C 73, 064903 (2006) [arXiv:hepph/0602249].
 [36] R. Baier, P. Romatschke, D. T. Son, A. O. Starinets and M. A. Stephanov, JHEP 0804, 100 (2008) [arXiv:0712.2451 [hepth]].
 [37] R. J. Fries, B. Muller and A. Schafer, Phys. Rev. C 78, 034913 (2008) [arXiv:0807.4333 [nuclth]].
 [38] M. Martinez and M. Strickland, arXiv:0902.3834 [hepph].
 [39] G. Torrieri, B. Tomasik and I. Mishustin, Phys. Rev. C 77, 034903 (2008) [arXiv:0707.4405 [nuclth]]; G. Torrieri, B. Tomasik and I. Mishustin, Acta Phys. Polon. B 39, 1733 (2008) [arXiv:0803.4070 [hepph]]; G. Torrieri and I. Mishustin, Phys. Rev. C 78, 021901 (2008) [arXiv:0805.0442 [hepph]].
 [40] H. Kouno, M. Maruyama, F. Takagi and K. Saito, Phys. Rev. D 41, 2903 (1990).
 [41] G. S. Denicol, T. Kodama, T. Koide and Ph. Mota, arXiv:0903.3595 [hepph]; G. S. Denicol, T. Kodama, T. Koide and Ph. Mota, arXiv:0907.4269 [hepph].
 [42] A. Monnai and T. Hirano, arXiv:0903.4436 [nuclth]; A. Monnai and T. Hirano, arXiv:0907.3078 [nuclth].
 [43] S. Bhattacharyya, V. E. Hubeny, S. Minwalla and M. Rangamani, JHEP 0802, 045 (2008) [arXiv:0712.2456 [hepth]].
 [44] P. Romatschke, arXiv:0906.4787 [hepth].
 [45] W. Israel, Annals Phys. 100, 310 (1976); W. Israel and J. M. Stewart, Annals Phys. 118, 341 (1979).
 [46] A. Bazavov et al., arXiv:0903.4379 [heplat].
 [47] P. Arnold, G. D. Moore and L. G. Yaffe, JHEP 0305, 051 (2003) [arXiv:hepph/0302165].
 [48] H. B. Meyer, Phys. Rev. D 76, 101701 (2007) [arXiv:0704.1801 [heplat]].
 [49] M. Natsuume and T. Okamura, Phys. Rev. D 77, 066014 (2008) [Erratumibid. D 78, 089902 (2008)] [arXiv:0712.2916 [hepth]].
 [50] M. A. York and G. D. Moore, Phys. Rev. D 79, 054011 (2009) [arXiv:0811.0729 [hepph]].
 [51] P. Huovinen and D. Molnar, Phys. Rev. C 79, 014906 (2009) [arXiv:0808.0953 [nuclth]]; D. Molnar and P. Huovinen, arXiv:0907.5014 [nuclth].
 [52] P. F. Kolb and U. W. Heinz, arXiv:nuclth/0305084.
 [53] B. Muller and K. Rajagopal, Eur. Phys. J. C 43, 15 (2005) [arXiv:hepph/0502174].
 [54] J. W. Chen and J. Wang, arXiv:0711.4824 [hepph].
 [55] D. FernandezFraile and A. G. Nicola, Phys. Rev. Lett. 102, 121601 (2009) [arXiv:0809.4663 [hepph]].
 [56] J. NoronhaHostler, J. Noronha and C. Greiner, arXiv:0811.1571 [nuclth].
 [57] A. Wiranata and M. Prakash, arXiv:0906.5592 [nuclth].
 [58] P. Arnold, C. Dogan and G. D. Moore, Phys. Rev. D 74, 085021 (2006) [arXiv:hepph/0608012].
 [59] G. D. Moore and O. Saremi, JHEP 0809, 015 (2008) [arXiv:0805.4201 [hepph]].
 [60] P. Romatschke and D. T. Son, arXiv:0903.3946 [hepph].
 [61] S. CaronHuot, Phys. Rev. D 79, 125009 (2009) [arXiv:0903.3958 [hepph]].
 [62] H. B. Meyer, Prog. Theor. Phys. Suppl. 174, 220 (2008) [arXiv:0805.4567 [heplat]].
 [63] K. Huebner, F. Karsch and C. Pica, Phys. Rev. D 78, 094501 (2008) [arXiv:0808.1127 [heplat]].
 [64] H. B. Meyer, JHEP 0808, 031 (2008) [arXiv:0806.3914 [heplat]].
 [65] I. Kanitscheider and K. Skenderis, JHEP 0904, 062 (2009) [arXiv:0901.1487 [hepth]].
 [66] A. Buchel, arXiv:0908.0108 [hepth].
 [67] J. I. Kapusta, arXiv:0809.3746 [nuclth].
 [68] G. Boyd, J. Engels, F. Karsch, E. Laermann, C. Legeland, M. Lutgemeier and B. Petersson, Nucl. Phys. B 469, 419 (1996) [arXiv:heplat/9602007].
 [69] G. S. Denicol, T. Kodama, T. Koide and Ph. Mota, J. Phys. G 36, 035103 (2009) [arXiv:0808.3170].
 [70] S. S. Gubser, A. Nellore, S. S. Pufu and F. D. Rocha, Phys. Rev. Lett. 101, 131601 (2008) [arXiv:0804.1950 [hepth]]; S. S. Gubser, S. S. Pufu and F. D. Rocha, JHEP 0808, 085 (2008) [arXiv:0806.0407 [hepth]].
 [71] A. Buchel, Nucl. Phys. B 820, 385 (2009) [arXiv:0903.3605 [hepth]].
 [72] U. Gursoy, E. Kiritsis, G. Michalogiorgakis and F. Nitti, arXiv:0906.1890 [hepph].
 [73] A. Buchel and C. Pagnutti, Nucl. Phys. B 816, 62 (2009) [arXiv:0812.3623 [hepth]].
 [74] I. Melo, B. Tomasik, G. Torrieri, S. Vogel, M. Bleicher, S. Korony and M. Gintner, arXiv:0902.1607 [nuclth].
 [75] For a recent review, see F. Becattini, arXiv:0901.3643 [hepph].
 [76] R. Hagedorn, Nuovo Cim. Suppl. 3, 147 (1965).
 [77] W. Broniowski, B. Hiller, W. Florkowski and P. Bozek, Phys. Lett. B 635, 290 (2006) [arXiv:nuclth/0510033]; W. Broniowski, P. Bozek, W. Florkowski and B. Hiller, PoS C FRNC2006, 020 (2006) [arXiv:nuclth/0611069].
 [78] B. Alver et al. [PHOBOS Collaboration], arXiv:0812.1172 [nuclex].
 [79] U. Heinz, private communication and H. Song, Ph.D. thesis, Ohio State University, August 2009.
 [80] H. Song and U. W. Heinz, arXiv:0909.1549 [nuclth].