1 Introduction

Hydrodynamic description of particle production in relativistic heavy-ion collisions

Mikołaj Chojnacki

The Henryk Niewodniczański

Institute of Nuclear Physics

Kraków, Poland

Thesis submitted for the Degree of Doctor of Philosophy in Physics

Prepared under the supervision of Prof. Wojciech Florkowski

Kraków, March 2009

STRESZCZENIE

Niniejsza praca prezentuje nowo opracowany model ewolucji hydrodynamicznej, który w połaczeniu z modelem statystycznej hadronizacji THERMINATOR służy nam do opisu zachowania silnie oddziałujacej materii wyprodukowanej w relatywistycznych zderzeniach cieżkich jonów. Nasze oryginalne podejście wykorzystano do wykonania dopasowań dla danych pochodzacych z eksperymentów realizowanych na akceleratorze RHIC (Relativistic Heavy Ion Collider w Brookhaven National Laboratory) przy najwyższej jego energii = 200 GeV, oraz do sformułowania przewidywań teoretycznych dla przyszłych eksperymentów cieżkojonowych przy wyższych energiach ( = 5.5 TeV, dla akceleratora LHC, skrót od Large Hadron Collider w CERN-ie).

Nasze wyniki odnosza sie do obserwabli jedno- i dwu-czastkowych w zakresie miekkiej fizyki ( 2 GeV). Opisujemy widma czastek w pedzie poprzecznym, współczynnik przepływu eliptycznego , oraz promienie HBT dla identycznych pionów (HBT jest skrótem od nazwisk Hanbury-Brown i Twiss). W ramach prac nad rozprawa skonstruowaliśmy nowe równanie stanu dla materii silnie oddziałujacej, które łaczy model gazu hadronowego z wynikami symulacji QCD na siatkach. Całość programów tworzy platforme obliczeniowa w skład której wchodzi kod hydrodynamiczny połaczony z modelem statystycznej hadronizacji THERMINATOR (skrót od THERMal heavy IoN generATOR).

Stosujac standardowy model optyczny Glaubera jako warunek poczatkowy dla ewolucji hydrodynamicznej, osiagneliśmy bardzo dobry opis danych eksperymentalnych uzyskanych na akceleratorze RHIC. W szczególności osiagneliśmy znacznie lepszy, od wcześniej uzyskanych w modelach hydrodynamicznych, opis promieni korelacyjnych. Dla przyszłych eksperymentów cieżkojonowych na LHC uzyskaliśmy przewidywania teoretyczne dotyczace miekkich obserwabli.

Zaproponowaliśmy również sposób rozwiazania tzw. zagadki HBT na RHIC-u. Sugerujemy zmodyfikowanie warunków poczatkowych i wprowadznie Gaussowskiego profilu gestości materii jako warunku poczatkowego dla hydrodynamiki. Taka modyfikacja prowadzi do szybszej formacji poprzecznego przepływu kolektywnego, co warunkuje uzyskanie wyjatkowo dobrej zgodności naszego modelu z danymi eksperymentalnymi.

Jako ostatni punkt, wprowadziliśmy do naszego modelu proces swobodnego strumieniowania czastek, który w połaczeniu z mechanizmem nagłej, chociaż opóźnionej w czasie termalizacji, tworzy nowe warunki poczatkowe dla kodu hydrodynamicznego. Wprowadzenie przedrównowagowej ewolucji pozwoliło nam na opóźnienie startu fazy hydrodynamicznej. Jest to pożadany efekt, który ma na celu unikniecie założenia o bardzo wczesnej termalizacji układu, które wydaje sie bardzo trudne do uzasadnienia na gruncie mikroskopowym. Podkreślmy, iż właczenie swobodnego strumieniowania czastek nie zmienia wysokiej zgodności uzyskanych wyników modelowych z danymi doświadczalnymi.

I know that this defies the law of gravity,

but, you see, I never studied law.

– Bugs Bunny

Physics isn‘t a religion.

If it were, we‘d have a much easier time raising money.

– Leon Lederman

— to my Parents —

Acknowledgments

I would like to express my deepest thanks to my supervisor Prof. Wojciech Florkowski for his invaluable help, guidance and patience during the course of this Thesis. Furthermore I wish to thank Wojciech Broniowski, Adam Kisiel and Piotr Bożek for a chance of working with them and being part of the team and to everybody in the Department of Theory of Structure of Matter (NZ41).

I am very grateful to all of my family and friends for their support throughout my journey into obtaining this Degree. I thank you all.

Research was supported by the Polish Ministry of Science and Higher Education grant N202 153 32/4247 (2007-2009).

## Chapter 1 Introduction

Nowadays, the relativistic hydrodynamics is regarded as the best theoretical framework for description of the spacetime evolution of strongly interacting matter produced in ultra-relativistic heavy-ion collisions [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 20, 21, 22, 23, 24] , for a recent review see [25]. In particular, the soft hadronic one-particle data describing the transverse-momentum spectra and the elliptic flow coefficient , collected in the RHIC experiments (Relativistic Heavy-Ion Collider at the Brookhaven National Laboratory), have been successfully explained in various approaches based on the perfect-fluid hydrodynamics. In fact, the explanation of the large value of the elliptic flow by perfect hydrodynamics is regarded as the evidence of early thermalization and suggests that the quark-gluon plasma created at RHIC is a strongly interacting system [26].

On the other hand, the approaches based on the hydrodynamics cannot reproduce the two-particle observables such as the pion correlation radii. The latter are commonly called the HBT radii – after Hanbury-Brown and Twiss who in 1950s showed that it was possible to determine the angular sizes of astronomical radio sources and stars from the correlations of signal intensities, rather than amplitudes. The difficulty of the consistent description of the one- and two-particle observables picked up the name ”HBT puzzle”. The HBT puzzle and the problem of the microscopic explanation of very fast thermalization of the matter produced at RHIC represent two issues that challenge the hydrodynamic picture.

In this Thesis we present our recently developed hydrodynamic model and use it to describe the RHIC data. We address both the one- and two-particle observables: the transverse momentum spectra, the elliptic flow coefficient , and the pion HBT radii. We suggest how the HBT puzzle as well as the early thermalization problem may be solved. Similarly to other approaches, our model is based on the perfect fluid hydrodynamics and includes the symmetry against Lorentz boosts along the beam axis, the so called boost-invariance. This restriction means that our results may be applied only to the central regions of relativistic heavy-ion collisions. On the other hand, our framework differs from other approaches in several important aspects, in particular, in the use of a different equation of state, modification of the initial conditions, two-body method of the calculation of the correlation functions, and different treatment of the final hadronic stage.

The main achievements of the Thesis are the following:

• The construction of the realistic equation of state for strongly interacting matter which interpolates between the hadron gas model and the results of the QCD lattice simulations. This equation of state describes the crossover phase transition, i.e., the transition where thermodynamic variables such as energy density or entropy density change very rapidly in the narrow range of the temperature, however, no real discontinuities in the behavior of the thermodynamic variables are present. Noticeably, our equation of state has no pronounced soft point where the sound velocity is very small and possibly drops to zero. It is known that the presence of such a soft point leads to the ratio of the HBT radii that is larger than unity – an effect which has not been confirmed by the experimental data.

• The boost-invariant hydrodynamic equations for baryon free matter have been rewritten in the very concise form which reduces the number of the independent equations to two. This may be done in the formal way if the range of the variable (the distance from the collision axis) is extended to the negative values. The applied procedure is a direct generalization of the formalism introduced earlier by Baym et al. in the studies of cylindrically symmetric systems with constant sound velocity. The new form of the hydrodynamic equations allows for the simple and natural inclusion of the boundary conditions at the origin of the system.

• The computational platform has been constructed which combines the hydrodynamic code with the statistical hadronization model THERMINATOR[27]. This is arranged in such a way that the information about the freeze-out hypersurface obtained from the hydrodynamic code is exported and treated as an input for THERMINATOR. We emphasize that our equation of state (in the region below the critical temperature MeV) describes the hadron gas with the same set of the hadronic species as that included in THERMINATOR. Hence, there is a smooth change between the hydrodynamic and statistical description of the produced matter. THERMINATOR is a Monte Carlo program simulating the decays of resonances. The Monte-Carlo method allows for the direct comparison of our model results with the experimental data. In particular, one can easily include various experimental cuts.

• The successful description of the soft hadronic RHIC data has been achieved with the standard initial conditions obtained from the optical limit of the Glauber model. By this we mean here that the use of the new equation of state helped to reduce the discrepancy between the theoretical and experimental ratio of the HBT radii and . In this version of our calculations we find , significantly closer to the experimental values than in earlier hydrodynamic studies.

• Predictions for the future heavy-ion collisions at LHC have been formulated. In the studied by us central region we expect that the transition from the RHIC energy, 200 GeV, to the LHC energy, 5.5 TeV, results essentially in a higher initial central temperature used as the input for the hydrodynamic calculation. Thus, one can make predictions for the collisions at the LHC energies using a set of values for which are higher than those used at RHIC. At RHIC we found = 320 MeV, hence for LHC we studied the cases = 400, 450, and 500 MeV. Our results for LHC indicate a moderate increase of the HBT radii and saturation of the pion elliptic flow (as compared to the RHIC experiments).

• The solution of the RHIC HBT puzzle has been proposed which suggests the use of the modified Gaussian-type initial conditions. We find that the choice of the initial condition in the form of a two-dimensional Gaussian profile for the transverse energy leads to a complete and consistent description of soft observables measured at RHIC. The transverse-momentum spectra, the elliptic-flow, and the HBT correlation radii, including the ratio are very well described.

• The processes of the free streaming of partons followed by the sudden equilibration were incorporated in the model. Those two processes deliver modified initial conditions for the hydrodynamics. In particular, the inclusion of the free-streaming stage allows for the delayed start of the hydrodynamic evolution, which is a desirable effect in the context of the early thermalization problem.

The Thesis is organized as follows. In Chapter 2 we describe the construction of our equation of state. In Chapter 3 we present the hydrodynamic equations and transform them to the form used in the numerical calculations. The initial conditions and the freeze-out prescription are introduced in Chapters 4 and 5, respectively. The fits to the RHIC data obtained with the standard initial conditions are presented in Chapter 6. The predictions for the LHC are given in Chapter 7. The solution of the RHIC HBT puzzle with modified Gaussian initial conditions is presented and discussed in Chapter 8. In that Section we also discuss the inclusion of the parton free-streaming as the pre-hydrodynamics stage. The Summary and four Appendices close the Thesis.

We use everywhere the natural units with . The signature of the metric tensor is .

The results discussed in this Thesis were published in the following articles:

• M. Chojnacki, W. Florkowski, and T. Csörgő,
Formation of Hubble - like flow in little bangs,
Phys. Rev. C71 (2005) 044902, (nucl-th/0410036).

• M. Chojnacki and W. Florkowski,
Characteristic form of boost-invariant and cylindrically asymmetric hydrodynamic equations,
Phys. Rev. C74 (2006) 034905, (nucl-th/0603065).

• M. Chojnacki and W. Florkowski,
Temperature dependence of sound velocity and hydrodynamics of ultra - relativistic heavy-ion collisions,
Acta Phys. Pol. B38 (2007) 3249, (nucl-th/0702030).

• M. Chojnacki, W. Florkowski, W. Broniowski, and A. Kisiel,
Soft heavy-ion physics from hydrodynamics with statistical hadronization: Predictions for collisions at = 5.5 TeV,
Phys. Rev. C78 (2008) 014905, arXiv:0712.0947 [nucl-th].

• W. Broniowski, M. Chojnacki, W. Florkowski, and A. Kisiel,
Uniform Description of Soft Observables in Heavy-Ion Collisions at
= 200 GeV,
Phys. Rev. Lett. 101 (2008) 022301, arXiv:0801.4361 [nucl-th].

• A. Kisiel, W. Broniowski, M. Chojnacki, and W. Florkowski,
Azimuthally sensitive femtoscopy in hydrodynamics with statistical hadronization from the BNL Relativistic Heavy Ion Collider to the CERN Large Hadron Collider,
Phys. Rev. C79 (2009) 014902, arXiv:0808.3363 [nucl-th].

• W. Broniowski, W. Florkowski, M. Chojnacki, A. Kisiel,
Free-streaming approximation in early dynamics of relativistic heavy-ion collisions,
submitted to Phys. Rev. C, arXiv:0812.3393 [nucl-th].

They were also presented during various international conferences including:

• M. Chojnacki,
Hubble-like Flows in Relativistic Heavy-Ion Collisions,
Acta Phys. Hung. A27 (2006) 331, (nucl-th/0510092).
18th International Conference On Ultrarelativistic Nucleus-Nucleus Collisions: Quark Matter 2005 (QM 2005).

• M. Chojnacki,
Cylindrically asymmetric hydrodynamic equations,
Acta Phys. Polon. B37 (2006) 3391, (nucl-th/0609060).
Cracow School Of Theoretical Physics: 46th Course 2006.

• M. Chojnacki,
Temperature-dependent sound velocity in hydrodynamic equations for relativistic heavy-ion collisions,
J. Phys. G35 (2008) 044074, arXiv:0709.1594 [nucl-th].
International Conference On Strangeness In Quark Matter (SQM 2007).

• W. Florkowski, M. Chojnacki, W. Broniowski, A. Kisiel,
Soft-hadronic observables for relativistic heavy-ion collisions at RHIC and LHC,
Acta Phys. Polon. B39 (2008) 1555, arXiv:0804.0974 [nucl-th].
Cracow Epiphany Conference On LHC Physics.

• W. Florkowski, W. Broniowski, M. Chojnacki, A. Kisiel,
Hydrodynamics and perfect fluids: Uniform description of soft observables in Au+Au collisions at RHIC,
arXiv:0811.3761 [nucl-th] and arXiv:0902.0377 [hep-ph].
38th International Symposium On Multiparticle Dynamics ISMD08.

• W. Florkowski, W. Broniowski, M. Chojnacki, A. Kisiel,
Solution of the RHIC HBT puzzle with Gaussian initial conditions,
arXiv:0812.4125 [nucl-th].
International Conference On Strangeness In Quark Matter (SQM 2008).

• W. Broniowski, W. Florkowski, M. Chojnacki, A. Kisiel,
Initial conditions for hydrodynamics: implications for phenomenology,
arXiv:0812.4935 [nucl-th].
IV Workshop on Particle Correlations and Femtoscopy.

• W. Florkowski, W. Broniowski, M. Chojnacki, A. Kisiel,
Consistent hydrodynamic description of one- and two-particle observables in relativistic heavy-ion collisions at RHIC,
arXiv:0901.1251 [nucl-th].
IV Workshop on Particle Correlations and Femtoscopy.

## Chapter 2 Thermodynamics of relativistic baryon-free matter

In our approach we concentrate on the description of the mid-rapidity region of ultra-relativistic heavy-ion collisions. Statistical analysis applied to the highest-energy RHIC data indicates that the baryon chemical potential at the chemical freeze-out is of about 25 MeV in this region [28, 29, 30, 31, 32]. The predictions of the statistical models for LHC give even smaller values, 0.8 MeV [33]. On the other hand, the expected temperature is of about 150 - 170 MeV, hence the ratio is small and in the hydrodynamic equations we can approximately assume that the baryon chemical potential vanishes. In this situation, as discussed in Ref. [34], the whole information about the equation of state is encoded in the temperature-dependent sound velocity . We assume that at low temperatures the sound velocity is given by the hadron-gas model with a complete set of hadronic resonances. In this case the function approaches zero as , which is the characteristics of the pion gas ( is the pion mass). On the other hand, at high temperatures our equation of state coincides with the recent lattice simulations of QCD [35]. The thermodynamic properties of the hadron gas and the quark-gluon plasma are discussed below in more detail in Sects. 2.1 and 2.2, respectively.

In the transition region between the hadron gas and the plasma, whose position is characterized by the critical temperature , different interpolations between the hadron gas result and the lattice result may be considered. In Ref. [36] we showed, however, that the most promising equation of state is based on the simplest interpolation between the hadron-gas model and the lattice data, see Fig. 2.1. This is so because the sound velocity function which does not exhibit a distinct minimum at the critical temperature leads to the relatively short evolution time and this effect helps to describe correctly the HBT data. The effects of different forms of the sound velocity are discussed in Sect. 2.3. Since the simplest interpolation is the best, most of the results presented in this Thesis are obtained with the sound velocity function shown in Fig 2.1.

The knowledge of the function allows us to determine all other thermodynamic properties of our system. This is achieved with the help of the following thermodynamic identities

 ε+P=Ts,dε=Tds,dP=sdT,c2s=dPdε, (2.1)

where is the energy density, is the pressure, is the temperature, and is the entropy density. In Fig. 2.2 we display the entropy and energy densities as functions of , and the pressure and sound velocity as functions of the energy density. These quantities follow directly from the assumed form of the function , shown in Fig. 2.1.

We note that other equations of state and their impact on the physical observables were recently studied in Ref. [37]. The result of that work was that the transverse-momentum spectra are quite insensitive to the assumed form of the equation of state. On the other hand, a noticeable dependence of the elliptic flow on the equation of state was observed. This dependence favored the strong first order phase transition. In Ref. [37] the effects concerning the HBT were not studied. Our recent work indicates that the constraints from the particle interferometry data exclude the strong first order phase transition since it leads to the unrealistically long evolution times.

We stress that by using the lattice results we take into account the non - perturbative aspects of the plasma behavior, which may be regarded as the effective inclusion of the strongly-interacting quark-gluon plasma; large deviations from the ideal-gas behavior directly indicate the non-negligible interactions present in the plasma. In particular, is significantly below the ideal-gas value of 1/3 also at temperatures way above . Note that in agreement with the present knowledge, no real phase transition is present in the system, but a smooth cross-over, therefore does not drop to zero at but remains a smooth function.

The hadron-gas model is based on the assumption that all known hadrons (including hadronic resonances) form a multicomponent perfect gas. In this case, well known formulas for the thermodynamic variables of the relativistic perfect gases [38, 39] may be applied. Of course, the hadron-gas model is applicable in the temperature region below the critical temperature = 170 MeV. With increasing temperature the density of hadrons becomes so large that they start overlapping and the very idea of hadrons breaks down. We note that only in the case of pions the quantum statistics is necessary. For other particles, which are much heavier and less abundant, the Boltzmann classical limit is sufficient [40].

#### 2.1.1 Pion gas

Pions are the most abundant particles produced in heavy-ion collisions. Because of their large multiplicity the statistical quantum effects cannot be neglected and one should use the Bose-Einstein distributions to describe the pion spectra. Additionally, since pions do not carry the baryon number and strangeness, their chemical potential may be assumed to be zero (possible values of the isospin chemical potential are very small, MeV, and are usually neglected). Therefore, in describing the thermodynamic properties of the pion gas () we use the relations for a massive boson gas with zero chemical potential, which are worked out in Appendix A, namely

 sπGi(T) = 12π2m3i∞∑κ=11κK3(miTκ), (2.2) PπGi(T) = 12π2T2m2i∞∑κ=11κ2K2(miTκ), (2.3) επGi(T) = 12π2Tm2i∞∑κ=11κ2[3TK2(miTκ)+miκK1(miTκ)], (2.4)

where the index specifies the isospin (distinguishes between , , ) and is the mass of the appropriate pion. The infinite sum over is the technical way to include the Bose-Einstein statistics. In the numerical calculations it is sufficient to include only the first four terms – the inverse powers of reduce the higher-order terms fast enough.

Knowing the expression for the entropy density of the pion gas, we calculate the sound velocity from the last equation in (2.1). In the approximation that the pion masses are equal () we obtain

 (c2s)πG(T)=2Tmπ∞∑κ=11κK3(mπTκ)∞∑κ=1[K2(mπTκ)+K4(mπTκ)]. (2.5)

It is interesting to observe the behavior of the sound velocity in the limit when temperature approaches zero. The arguments of the Bessel functions tend to infinity and we can use the power series expansion to find that

 (c2s)πG(T→0)=Tmπ∑∞κ=1[κ−32e−mπTκ+⋯]∑∞κ=1[κ−12e−mπTκ+⋯]→Tmπ, (2.6)

see Eq. (D.8) from Appendix D. On the other hand, when temperature is large () we expand the modified Bessel functions according to (D.7) and get the following expression

 (c2s)πG(T→∞)=13∑∞κ=1[κ−4−18(mπTκ)2+18(mπTκ)4+⋯]∞∑κ=1[κ−4−124(mπTκ)2+196(mπTκ)4+⋯]→13. (2.7)

The last limit illustrates the expected result – at very high temperatures the massive pion gas behaves effectively like a massless gas with and its sound velocity squared equals . (See Tables A.4 and A.5, where the complete formulas for pressure and energy density in various limits are given).

#### 2.1.2 Classical gas

Strictly speaking, all particles produced in heavy-ion collisions obey the quantum statistics, however for all particles other than pions the effects of quantum statistics are numerically negligible. Thus, to calculate the thermodynamic functions of the hadrons other than pions we may use the following formulas

 sCGi(T) = gi2π2m3iK3(miT), (2.8) PCGi(T) = gi2π2T2m2iK2(miT), (2.9) εCGi(T) = gi2π2Tm2i[3TK2(miT)+miK1(miT)], (2.10)

where is the degeneration factor which holds the information about the spin degeneration of -th particle. The sound velocity in the classical massive case is defined as follows

 (c2s)CGi(T)=2TmiK3(miT)K2(miT)+K4(miT), (2.11)

and has the same asymptotic features as the pion gas, namely it tends to the value for very high temperatures and becomes proportional to if the temperature tends to zero.

The massive hadron gas (HG) model is the sum of both massive pion gas and massive classical gas of all other hadrons. In our study, the information on the mass and spin of individual particles comes from the input file to the SHARE program [41]. The table particles.data holds parameters for 371 particles consisting from u, d and s quarks. Thus, for all hadrons we may write

 sHG(T) = Missing or unrecognized delimiter for \left (2.12) (c2s)HG(T) = sHGTdTd(sHG), (2.13)

and all other thermodynamic quantities (, and ) are expressed in the analogous way.

### 2.2 Quark-Gluon Plasma

In our approach we use the results of the lattice simulations of QCD at the finite temperature presented in Ref. [35], see Fig. 2.3. They were obtained for physical masses of the light quarks and the strange quark.

The pressure data obtained from the QCD lattice calculation in Ref. [35] have been recently parameterized for the case of in Ref. [42], see Fig. 2.4. The parameterization has the form

 P=cT4σ(Tc/T),σ(g)=1+e−ab1+eg−abe−λg, (2.14)

where the dimensionless fit parameters equal: , , and .

Taking the lattice result at the face value, one expects that the sound velocity significantly drops down in the region . Similar behavior, with reaching zero, is expected in the case of the first order phase transition where the changes of the energy happen at constant pressure. However, the lattice simulations suggest that for three massive quarks with realistic masses we deal with the cross-over rather than with the first order phase transition, hence the sound velocity remains finite, as is consistently shown in Fig. 2.1.

### 2.3 Modeling the crossover phase transition

The exact values of the sound velocity in the region are poorly known. The lattice calculations are not very much reliable for and, at the same time, the use of the hadron gas model with vacuum parameters becomes unrealistic for large densities (temperatures). The authors of Ref. [35] state that in the hadronic phase the lattice spacing is larger than 0.3 fm and the lattice artifacts cannot be controlled in this region. In this situation, it is practical to consider different interpolations between the lattice and hadron-gas results and to analyze the physical effects of a particular choice of the interpolating function. This type of the study was performed in Ref. [36]. Here we shortly discuss the main conclusions of this analysis.

In Ref. [36] we considered three different sound-velocity functions . Below, we refer to these three options as to the cases I, II and III, see Fig. 2.5. In the case I, we use the sound-velocity function which agrees with the ideal hadron gas model of Ref. [34] in the temperature range and with the lattice result in the temperature range 1. In the region close to the critical temperature, , a simple interpolation between the two results is used. We have checked that such a simple interpolation yields directly the entropy density consistent with the lattice result. Namely, the use of the thermodynamic relation

 s(T)=s(Tmin)exp⎡⎢⎣T∫TmindT′T′c2s(T′)⎤⎥⎦, (2.15)

relating the entropy density with the sound velocity for zero baryon chemical potential, gives the function which agrees with the lattice result at high temperatures, at [35].

In the cases II and III, the sound-velocity interpolating functions have a distinct minimum at . Comparing to the case I [with and ], the value of the sound velocity at is reduced by 25 % in the case II [where and ], and by 50% in the case III [where and ]. From Eq. (2.15) one concludes that the decrease of the sound velocity at leads to the increase of entropy density for high temperatures. Hence, in order to have the same value of the entropy density at high temperatures, a decrease of the sound velocity function in the region should be compensated by its increase in a different temperature range. For our interpolating functions in the cases II and III we assume that the values of in the range are slightly higher than in the case I, see Fig. 2.5. Such modifications may be regarded as the parameterization of the repulsive van der Waals forces in the hadron gas. The values of the maxima are chosen in such a way that the entropy densities for three considered cases are consistent with the lattice result, see the upper left panel of Fig. 2.6 where the functions are shown. We stress that in the three considered cases the values of in the temperature range remain significantly below the massless limit . Such a limiting value is implicitly used in many hydrodynamic codes assuming the equation of state of non-interacting massless quarks and gluons for , see for example the extended 3+1 hydrodynamic model of Ref. [23].

The studies of the hydrodynamic spacetime evolution of matter described by the equations of state I, II and III were performed in Ref. [36]. As expected, we found that the drop of at leads to the prolonged time of the evolution, hence leads to the increase of the ratio of the HBT radii. This behavior is in contradiction with the observed data which indicate . Thus, we have decided to exclude the cases II and III from further analysis and to restrict our consideration to the case I.

## Chapter 3 Relativistic hydrodynamics of perfect fluid

In this Chapter we present the main ingredients of our hydrodynamic model. We start with the general formulation of the relativistic hydrodynamics in the case of vanishing baryon number. Next, we implement the idea of boost-invariance that allows us to restrict our considerations to the plane . In the subsequent Sections of this Chapter we show that the hydrodynamic equations may be rewritten in the form where only two equations are independent (at the expense of the formal extension of the variable to negative values) and the boundary conditions at the origin are automatically fulfilled. Such a form, being the direct generalization of the approach introduced in Ref. [45], turned out to be very convenient in the numerical analyses.

### 3.1 Hydrodynamic equations for baryon free matter

As the system reaches local thermodynamic equilibrium its further evolution is governed by the conservation laws for energy and momentum, which can be expressed by the formula

 ∂μTμν=0, (3.1)

where the energy-momentum tensor has the form

 Tμν=(ε+P)uμuν−Pgμν. (3.2)

Here is the energy density, is the pressure, is the metric tensor (we use the convention where ) and is the four-velocity,

 uμ=γ(1,v). (3.3)

In Eq. (3.3) is the local three-velocity of the fluid and is the Lorentz factor

 γ=(1−v2)−12. (3.4)

In general situations, in the relativistic systems the energy density and pressure depend on the temperature and baryon number density. This requires that the baryon number conservation law,

 ∂μjμB=0, (3.5)

should be considered together with (3.1). However, in the case of the central rapidity region of ultra-relativistic heavy-ion collisions the dominating degrees of freedom are mesons (initially gluons), hence the net baryon number is very close to zero. In such a case the local value of baryon chemical potential is negligible, , and we can express the entropy density , pressure , and energy density by temperature alone. This allows us to use the thermodynamic identities (2.1) and write the conservation laws (3.1) in the form.

 ∂μ(Tsuμuν)=∂νP. (3.6)

After performing simple transformations shown explicitly in Appendix B.2 Eq. (3.6) leads to the two formulas

 uμ∂μ(Tuν) = ∂νT, (3.7) ∂μ(suμ) = 0. (3.8)

Eq. (3.7) is the acceleration equation and represents the relativistic analog of the Euler equation known from the classical hydrodynamics. Eq. (3.8) states that the evolution is adiabatic (entropy is conserved). In the non-covariant notation the hydrodynamic equations (3.7) and (3.8) are expressed by the following expressions [45]

 ∂∂t(Tγv)+∇(Tγ) = v×[∇×(Tγv)], (3.9) ∂∂t(sγ)+∇(sγv) = 0. (3.10)

The four equations above can be written in the Cartesian coordinates in the equivalent form as

 (v2−1)∂lnT∂t+dlnTdt+11−v2vdvdt = 0, (3.11) (1−v2)(vy∂lnT∂x−vx∂lnT∂y)+vydvxdt−vxdvydt = 0, (3.12) (1−v2)∂lnT∂z+vzdlnTdt+dvzdt+vz1−v2vdvdt = 0, (3.13) dlnsdt+v1−v2dvdt+∂vx∂x+∂vy∂y+∂vz∂z = 0, (3.14)

where the total time derivative is defined by the equation

 ddt=∂∂t+v⋅∇. (3.15)

The details of the transformations leading from (3.6) to (3.11) - (3.14) can be found in Appendix B.3. Equations (3.11) - (3.14) do not form a closed system of equations since they contain five independent variables , , , and . An additional equation is needed to close them, i.e., the equation of state is required which introduces the relation between and . Alternatively, the equation of state may be included by the use of the temperature dependent sound velocity

 c2s(T)=∂P∂ϵ=sT∂T∂s, (3.16)

The temperature dependent sound velocity for strongly interacting matter that is used in our work was discussed thoroughly in Chapt. 2 and is plotted in Fig. 2.1.

### 3.2 Implementation of boost-invariance

The experimental data collected at RHIC by the BRAHMS Collaboration [46] suggests that in the midrapidity region the particle yields do not vary much with rapidity. Thus, we can assume that number of particles per unit rapidity in the range of is essentially constant and the midrapidity region (central region) is boost-invariant. This symmetry demands that the longitudinal component of the velocity has the form of the Bjorken flow [47],

 vz=zt, (3.17)

and the thermodynamic scalar variables like temperature or entropy density are functions of the longitudinal proper time and the transverse coordinates and . In practice, these properties mean that we can solve the hydrodynamic equations for and by using the appropriate Lorentz transformations we obtain solutions for .

Adopting the procedure outlined above and restricting our considerations to the plane we observe that Eq. (3.13) is automatically fulfilled and we are left with only three independent equations

 (v2⊥−1)∂lnT∂t+dlnTdt+11−v2⊥v⊥dv⊥dt = 0, (3.18) (1−v2⊥)(vy∂lnT∂x−vx∂lnT∂y)+vydvxdt−vxdvydt = 0, (3.19) dlnsdt+v⊥1−v2⊥dv⊥dt+∂vx∂x+∂vy∂y+1t = 0, (3.20)

where is the transverse velocity. The equations describing transverse evolution can be rewritten in the cylindrical coordinates which are convenient for further analysis [48] (details can be found in Appendix B.4.1 and B.5.1)

 (v2⊥−1)∂lnT∂t+dlnTdt+11−v2⊥v⊥dv⊥dt=0, (3.21) (1−v2⊥)(v⊥sinα∂lnT∂r−v⊥cosαr∂lnT∂ϕ)−v2⊥(dαdt+v⊥sinαr)=0, (3.22) v⊥dlnsdt+11−v2⊥dv⊥dt−∂v⊥∂t−v2⊥sinα∂α∂r+v2⊥cosαr(∂α∂ϕ+1)+v⊥t=0.

Here is the distance from the beam axis, , and is the azimuthal angle, . These two coordinates parameterize the plane . The angle is the function describing direction of the flow, , see Fig. 3.1. The differential operator used in Eqs. (3.21) - (LABEL:eqn:hydro_s_Cyl_BI) is defined by the formula

 ddt=∂∂t+v⊥cosα∂∂r+v⊥sinαr∂∂ϕ. (3.24)

### 3.3 Characteristic form of hydrodynamic equations

In order to obtain the form of the hydrodynamic equations which is convenient for numerical studies we introduce new independent variables. In this respect we follow the method originally proposed by Baym et al. in Ref. [45] and the new form of the equations will be called the characteristic form. We introduce the potential defined as

 dΦ=1csdlnT=csdlns, (3.25)

and the transverse rapidity . The new form of the hydrodynamic equations is expressed by the dimensionless auxiliary functions and defined by the equations

 A±=Φ±η⊥. (3.26)

The sum and difference of Eq. (3.21) and (LABEL:eqn:hydro_s_Cyl_BI) together with Eq. (3.22) take the form

 ∂A±∂t + v⊥±cs1±csv⊥[cosα∂A±∂r+sinαr∂A±∂ϕ] (3.27) − cs1±csv⊥[v⊥sinα∂α∂r−v⊥cosαr(∂α∂ϕ+1)−1t]=0, ∂α∂t + v⊥cosα∂α∂r+v⊥sinαr(∂α∂ϕ+1) (3.28) − cs(1−v2⊥)v⊥[sinα∂Φ∂r−cosαr∂Φ∂ϕ]=0,

where the transverse velocity and the potential are expressed through the auxiliary functions as follows

 v⊥=tanh(A+−A−2),Φ=A++A−2. (3.29)

Temperature and all other temperature dependent variables, e.g., the sound velocity, can be also calculated from the functions ,

 T=TΦ(Φ)=TΦ(A++A−2), (3.30) cs(T)=cs[TΦ(A++A−2)]. (3.31)

Here we have introduced the subscripts to make clear what kind of the argument is expected for a given function. For example, the temperature may be considered as a function of entropy density or . In those two cases one should use the functions or , respectively. If the equation of state is known all such functions can be easily calculated and Eqs. (3.27) and (3.28) may be used to determine three unknown functions , , and .

### 3.4 Boundary conditions

Having in mind the heavy-ion collisions at RHIC and at LHC we consider the collisions of identical nuclei, Au+Au or Pb+Pb, which collide moving initially along the -axis. The positions of the centers of nuclei depend on the impact parameter and for non-central collision they may be located in the transverse plane at and at , see Fig. 3.2. The distribution of matter created just after the collision is not cylindrically symmetric and has an ellipsoidal shape (one commonly speaks of an almond shape). With the use of the coordinate system explained above the transverse velocity must vanish at the origin of the system namely

 v⊥(t,r=0,ϕ)=0. (3.32)

Also the gradients with respect to the distance of temperature and entropy density must converge to zero at

 ∂T(t,r,ϕ)∂r\lx@stackrelr→0⟶0,∂s(t,r,ϕ)∂r\lx@stackrelr→0⟶0. (3.33)

The boundary conditions (3.32) and (3.33) can be naturally fulfilled by the following Ansatz

 ⎧⎪⎨⎪⎩A+(t,r,ϕ)=A(t,r,ϕ),A−(t,r,ϕ)=A(t,−r,ϕ),α(t,−r,ϕ)=α(t,r,ϕ).r>0 (3.34)

The domain of the transverse distance has been extended to the negative values of . The two functions and are replaced in this way by a single function . At the same time the function has been generalized to have negative arguments by the condition that it is a symmetric function of , see Fig. 3.3. With the help of the definitions (3.34), Eqs. (3.27) may be reduced to a single equation for the function ,

 ∂A∂t + v⊥+cs1+csv⊥[cosα∂A∂r+sinαr∂A∂ϕ] (3.35) − cs1+csv⊥[v⊥sinα∂α∂r−v⊥cosαr(∂α∂ϕ+1)−1t]=0.

The transverse velocity and potential from Eq. (3.29) are now expressed as

 v⊥(t,r,ϕ) = tanh(A(t,r,ϕ)−A(t,−r,ϕ)2), (3.36) Φ(t,r,ϕ) = A(t,r,ϕ)+A(t,−r,ϕ)2. (3.37)

In addition, the use of cylindrical coordinates implies the periodicity of all functions in angle thus creating another set of periodic boundary conditions

 {A(t,r,ϕ=0)=A(t,r,ϕ=2π)α(t,r,ϕ=0)=α(t,r,ϕ=2π). (3.38)

## Chapter 4 Initial conditions

Along with the formulation of the hydrodynamic equations that govern the evolution of matter we must also specify a set of initial conditions that are required to unambiguously solve such equations – the hydrodynamic equations are first-order partial differential equations. The physical systems studied here are formed in the collisions of two identical gold nuclei (RHIC experiment) or two identical lead nuclei (the future heavy-ion program at LHC). In this Chapter we discuss in more detail the initial conditions used by us to analyze the collisions at RHIC and LHC. At first we discuss the most common form of the initial conditions that is based on the Glauber model. Next, we present the modified initial conditions where the initial energy density has a form of the two-dimensional Gaussian in the transverse plane. Finally, we present our method of including the parton free-streaming as a process which precedes the hydrodynamic evolution. In particular, we show how to match the free-streaming stage with the hydrodynamic stage.

### 4.1 Standard initial conditions

In the calculations presented in Chapter 6 and 7 we assume that the initial entropy density of the particles produced at the transverse position point is proportional to a source profile obtained from the Glauber approach. This is done in the similar way to other hydrodynamic calculations. Specifically, for the particle source profile we use a mixed model [49, 50], with a linear combination of the wounded-nucleon density and the density of binary collisions , namely

 s(x⊥)∝ρ(x⊥)=1−κ2ρWN(x⊥)+κρBC(x⊥). (4.1)

The case corresponds to the standard wounded-nucleon model [51], while would include the binary collisions only. The PHOBOS analysis [50] of the particle multiplicities yields at and at . Since the density profile from the binary collisions is steeper than from the wounded nucleons, increased values of yield steeper density profiles, which in turn result in steeper temperature profiles.

In our hydrodynamic code, the initial conditions are specified for the temperature profile which takes the form

 T(τi,x⊥)=TS[siρ(x⊥)ρ(0)], (4.2)

where is the inverse function to the function , and is the initial entropy at the center of the system. The initial central temperature equals . Throughout this work, the initial time for the start of the hydrodynamic evolution is denoted by .

The wounded-nucleon and the binary-collisions densities in Eq. (4.2) are obtained from the optical limit of the Glauber model, which is a very good approximation for not too peripheral collisions [52]. The standard formulas are [51]

 ρWN(x⊥) = Missing or unrecognized delimiter for \left (4.3) + TA(−b2+x⊥){1−[1−σinATA(b2+x⊥)]A}

and

 ρBC(x⊥)=σinTA(b2+x⊥)TA(−b2+x⊥). (4.4)

In Eqs. (4.3) and (4.4) is the impact vector, is the nucleon-nucleon total inelastic cross section, and is the nucleus thickness function

 TA(x,y)=∫dzρ(x,y,z). (4.5)

For RHIC energies we use the value , while for LHC we take . The function in Eq. (4.5) is the nuclear density profile given by the Woods-Saxon function with the conventional choice of the parameters:

 ρ0 = 0.17fm−3, r0 = (1.12A1/3−0.86A−1/3)fm, a = 0.54fm. (4.6)

The values of the atomic mass are: 197 for RHIC (gold nuclei) and 208 for LHC (lead nuclei). The value of the impact parameter in Eqs. (4.3) and (4.4) depends on the considered centrality class.

Besides the initial temperature profile (4.2) we also specify the initial transverse flow profile,

 v(τi,r,ϕ) = H0r√1+H20r2, α(τi,r,ϕ) = 0. (4.7)

The results presented below are obtained with fm. The very small value of the parameter means that the system is practically at rest at the moment when the hydrodynamic evolution starts. However, nonzero improves the stability of the numerical method.

In our numerical calculations we use the adaptive method of lines in the way as it is implemented in the MATHEMATICA package. The and directions are discretized, and the integration in time is treated as solving of the system of ordinary differential equations. Typically, we use grids with = 0.25 fm in the direction and 6 degrees in the direction. This method fails in the case where the shock waves are formed, however with our regular equation of state and the regular initial conditions such shocks are not present.

We stress that the shape of the initial condition (4.2) is important, as it determines the development of the radial and elliptic flow, thus affecting such observables as the -spectra, , and the femtoscopic features. On qualitative grounds, sharper profiles lead to more rapid expansion. Several effects should be considered here. Firstly, as discussed in Ref. [3], hydrodynamics may start a bit later, when the profile is less eccentric than originally due to initial free-streaming of partons in the pre-hydro phase. On the other hand, statistical fluctuations in the distribution of the Glauber sources (wounded nucleons, binary collisions) [53, 54, 55, 56, 57, 58, 59] lead to a significant enhancement of the eccentricity, especially at low values of the impact parameter. Thus the initial eccentricity may in fact be smaller or larger than what follows from the application of the Glauber model. This contributes to the systematic model uncertainty at the level of about 10-20%. This uncertainty could only be reduced by the employment of a realistic model of the pre-hydrodynamic evolution. With this uncertainty in place, one should not expect or demand a better agreement with the physical observables than at the corresponding level of 10-20%.

The results obtained with the introduced here standard initial conditions are presented and discussed in Chapters 6 and 7.

### 4.2 Gaussian initial conditions

The situation described in the previous Section corresponds to the typical case where the hydrodynamic evolution is initiated from a source profile generated by the Glauber model, with the initial central temperature or energy serving as a free parameter. The initial density and possibly the initial flow profiles may be also provided by the early partonic dynamics, for instance by the Color Glass Condensate (CGC) [60, 61]. In practice, however, the theory of the partonic stage carries some uncertainty in its parameters, which influences our knowledge of the details of early dynamics. Having in mind such uncertainties we decided to try another class of the initial conditions, i.e., in our calculations we also include the case where the initial energy profiles in the transverse plane have the Gaussian shape.

The modified Gaussian parameterization of the initial energy density profile at the initial proper time has the following form

 ε(x⊥)=εiexp(−x22a2−y22b2). (4.8)

The energy density profile (4.8) determines the initial temperature profile that is used in the hydrodynamic code,

 T(τi,x⊥)=Tε[εiexp(−x22a2−y22b2)], (4.9)

The initial central temperature , which may depend on the centrality, is denoted by and is a free parameter of our approach.

The results obtained with such initial conditions are presented and discussed in Sect. 8.1.

### 4.3 Free streaming

The very early start of hydrodynamics means that the matter equilibrates very fast. Such a short thermalization time is difficult to explain on the microscopical grounds and inspires hot discussion about the nature of matter produced at RHIC. To avoid the problem with the early thermalization we consider also the scenario where the hydrodynamic evolution is preceded by the free-streaming stage. In this version of our calculations we assume that partons behave as free particles in the proper time interval  fm. At fm the sudden equilibration takes place which is described with the help of the Landau matching conditions. The global picture is as follows: early phase (CGC) generating partons at time fm, partonic free streaming until fm, hydrodynamic evolution until the freeze-out at temperature , free streaming of hadrons and decay of resonances. We note that similar ideas have been described by Sinyukov et al. in Refs. [62, 63].

In the remaining part of this Section we describe in more detail the free-streaming stage. We assume that massless partons are formed at the initial proper time and move along straight lines at the speed of light until the proper time when free streaming ends, . We introduce the space-time rapidities

 η0=12logt0−z0t0+z0 (4.10)

and

 η∥=12logt−zt+z. (4.11)

Elementary kinematics [62] links the positions of a parton on the initial and final hypersurfaces and its four-momentum

 pμ=(pTcoshy,pTcosϕp,pTsinϕp,pTsinhy), (4.12)

where and are the parton‘s rapidity and transverse momentum:

 τsinh(η∥−y)=τ0sinh(η0−y),x=x0+Δcosϕp,y=y0+Δsinϕp,Δ=τcosh(y−η∥)−√τ20+τ2sinh2(y−η∥). (4.13)

Thus the phase-space density of partons at the proper times and are related

 d6Ndyd2pTdη∥dxdy = ∫dη0dx0dy0d6Ndyd2pTdη0dx0dy0× (4.14) δ(η0−y−arcsinh[ττ0sinh(η∥−y)])× δ(x−x0−Δcosϕp)δ(y−y0−Δsinϕp).

In our approach we restrict ourselves to the boost-invariant systems, hence it is reasonable to assume the following factorized form of the initial distribution of partons,

 d6Ndyd2pTdη0dx0dy0=n(x0,y0)F(y−η0,pT), (4.15)

where is the transverse density of partons obtained again from the GLISSANDO model, namely

 n(x0,y0)=n0exp(−x202a2−y202b2). (4.16)

When the rapidity emission profile is focused near , for instance if we have , with , and if , then the kinematic condition (4.13) effectively transforms it into

 F∼exp⎛⎜ ⎜⎝arcsinh2[ττ0sin(y−η∥)]2δy2⎞