Hydrodynamic description of particle production in relativistic heavyion collisions
Mikołaj Chojnacki
The Henryk Niewodniczański
Institute of Nuclear Physics
Polish Academy of Sciences
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 CERNie).
Nasze wyniki odnosza sie do obserwabli jedno i dwuczastkowych 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 HanburyBrown 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 RHICu. 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 (20072009).
Contents:
 1 Introduction
 2 Relativistic thermodynamics
 3 Relativistic hydrodynamics of perfect fluid
 4 Initial conditions
 5 Freezeout prescription
 6 Softhadronic observables at RHIC
 7 Predictions for LHC
 8 Uniform description of the RHIC data
 9 Summary
 A Relativistic thermodynamics of perfect gases
 B Properties of hydrodynamic equations
 C Notation
 D Mathematical supplement
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 ultrarelativistic heavyion 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 oneparticle data describing the transversemomentum spectra and the elliptic flow coefficient , collected in the RHIC experiments (Relativistic HeavyIon Collider at the Brookhaven National Laboratory), have been successfully explained in various approaches based on the perfectfluid 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 quarkgluon plasma created at RHIC is a strongly interacting system [26].
On the other hand, the approaches based on the hydrodynamics cannot reproduce the twoparticle observables such as the pion correlation radii. The latter are commonly called the HBT radii – after HanburyBrown 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 twoparticle 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 twoparticle 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 boostinvariance. This restriction means that our results may be applied only to the central regions of relativistic heavyion 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, twobody 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 boostinvariant 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 freezeout 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 MonteCarlo 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 heavyion 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 Gaussiantype initial conditions. We find that the choice of the initial condition in the form of a twodimensional Gaussian profile for the transverse energy leads to a complete and consistent description of soft observables measured at RHIC. The transversemomentum spectra, the ellipticflow, 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 freestreaming 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 freezeout 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 freestreaming as the prehydrodynamics 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, (nuclth/0410036). 
M. Chojnacki and W. Florkowski,
Characteristic form of boostinvariant and cylindrically asymmetric hydrodynamic equations,
Phys. Rev. C74 (2006) 034905, (nuclth/0603065). 
M. Chojnacki and W. Florkowski,
Temperature dependence of sound velocity and hydrodynamics of ultra  relativistic heavyion collisions,
Acta Phys. Pol. B38 (2007) 3249, (nuclth/0702030). 
M. Chojnacki, W. Florkowski, W. Broniowski, and A. Kisiel,
Soft heavyion physics from hydrodynamics with statistical hadronization: Predictions for collisions at = 5.5 TeV,
Phys. Rev. C78 (2008) 014905, arXiv:0712.0947 [nuclth]. 
W. Broniowski, M. Chojnacki, W. Florkowski, and A. Kisiel,
Uniform Description of Soft Observables in HeavyIon Collisions at
= 200 GeV,
Phys. Rev. Lett. 101 (2008) 022301, arXiv:0801.4361 [nuclth]. 
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 [nuclth]. 
W. Broniowski, W. Florkowski, M. Chojnacki, A. Kisiel,
Freestreaming approximation in early dynamics of relativistic heavyion collisions,
submitted to Phys. Rev. C, arXiv:0812.3393 [nuclth].
They were also presented during various international conferences including:

M. Chojnacki,
Hubblelike Flows in Relativistic HeavyIon Collisions,
Acta Phys. Hung. A27 (2006) 331, (nuclth/0510092).
18th International Conference On Ultrarelativistic NucleusNucleus Collisions: Quark Matter 2005 (QM 2005). 
M. Chojnacki,
Cylindrically asymmetric hydrodynamic equations,
Acta Phys. Polon. B37 (2006) 3391, (nuclth/0609060).
Cracow School Of Theoretical Physics: 46th Course 2006. 
M. Chojnacki,
Temperaturedependent sound velocity in hydrodynamic equations for relativistic heavyion collisions,
J. Phys. G35 (2008) 044074, arXiv:0709.1594 [nuclth].
International Conference On Strangeness In Quark Matter (SQM 2007). 
W. Florkowski, M. Chojnacki, W. Broniowski, A. Kisiel,
Softhadronic observables for relativistic heavyion collisions at RHIC and LHC,
Acta Phys. Polon. B39 (2008) 1555, arXiv:0804.0974 [nuclth].
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 [nuclth] and arXiv:0902.0377 [hepph].
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 [nuclth].
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 [nuclth].
IV Workshop on Particle Correlations and Femtoscopy. 
W. Florkowski, W. Broniowski, M. Chojnacki, A. Kisiel,
Consistent hydrodynamic description of one and twoparticle observables in relativistic heavyion collisions at RHIC,
arXiv:0901.1251 [nuclth].
IV Workshop on Particle Correlations and Femtoscopy.
Chapter 2 Thermodynamics of relativistic baryonfree matter
In our approach we concentrate on the description of the midrapidity region of ultrarelativistic heavyion collisions. Statistical analysis applied to the highestenergy RHIC data indicates that the baryon chemical potential at the chemical freezeout 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 temperaturedependent sound velocity . We assume that at low temperatures the sound velocity is given by the hadrongas 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 quarkgluon 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 hadrongas 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
(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 transversemomentum 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 stronglyinteracting quarkgluon plasma; large deviations from the idealgas behavior directly indicate the nonnegligible interactions present in the plasma. In particular, is significantly below the idealgas 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 crossover, therefore does not drop to zero at but remains a smooth function.
2.1 Hadron gas
The hadrongas 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 hadrongas 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 heavyion collisions. Because of their large multiplicity the statistical quantum effects cannot be neglected and one should use the BoseEinstein 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
(2.2)  
(2.3)  
(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 BoseEinstein statistics. In the numerical calculations it is sufficient to include only the first four terms – the inverse powers of reduce the higherorder 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
(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
(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
(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 heavyion 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
(2.8)  
(2.9)  
(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
(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.
2.1.3 Massive hadron gas
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
(2.12)  
(2.13) 
and all other thermodynamic quantities (, and ) are expressed in the analogous way.
2.2 QuarkGluon 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
(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 crossover 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 hadrongas 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 soundvelocity 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 soundvelocity 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
(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 soundvelocity 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 noninteracting 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 boostinvariance 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
(3.1) 
where the energymomentum tensor has the form
(3.2) 
Here is the energy density, is the pressure, is the metric tensor (we use the convention where ) and is the fourvelocity,
(3.3) 
In Eq. (3.3) is the local threevelocity of the fluid and is the Lorentz factor
(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,
(3.5) 
should be considered together with (3.1). However, in the case of the central rapidity region of ultrarelativistic heavyion 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.
(3.6) 
After performing simple transformations shown explicitly in Appendix B.2 Eq. (3.6) leads to the two formulas
(3.7)  
(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 noncovariant notation the hydrodynamic equations (3.7) and (3.8) are expressed by the following expressions [45]
(3.9)  
(3.10) 
The four equations above can be written in the Cartesian coordinates in the equivalent form as
(3.11)  
(3.12)  
(3.13)  
(3.14) 
where the total time derivative is defined by the equation
(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
(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 boostinvariance
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 boostinvariant. This symmetry demands that the longitudinal component of the velocity has the form of the Bjorken flow [47],
(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
(3.18)  
(3.19)  
(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)
(3.21)  
(3.22)  
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
(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
(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
(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
(3.27)  
(3.28)  
where the transverse velocity and the potential are expressed through the auxiliary functions as follows
(3.29) 
Temperature and all other temperature dependent variables, e.g., the sound velocity, can be also calculated from the functions ,
(3.30)  
(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 heavyion 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 noncentral 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
(3.32) 
Also the gradients with respect to the distance of temperature and entropy density must converge to zero at
(3.33) 
The boundary conditions (3.32) and (3.33) can be naturally fulfilled by the following Ansatz
(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 ,
(3.35)  
The transverse velocity and potential from Eq. (3.29) are now expressed as
(3.36)  
(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
(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 firstorder 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 heavyion 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 twodimensional Gaussian in the transverse plane. Finally, we present our method of including the parton freestreaming as a process which precedes the hydrodynamic evolution. In particular, we show how to match the freestreaming 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 woundednucleon density and the density of binary collisions , namely
(4.1) 
The case corresponds to the standard woundednucleon 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
(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 woundednucleon and the binarycollisions 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]
(4.3)  
and
(4.4) 
In Eqs. (4.3) and (4.4) is the impact vector, is the nucleonnucleon total inelastic cross section, and is the nucleus thickness function
(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 WoodsSaxon function with the conventional choice of the parameters:
(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,
(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 freestreaming of partons in the prehydro 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 1020%. This uncertainty could only be reduced by the employment of a realistic model of the prehydrodynamic 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 1020%.
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
(4.8) 
The energy density profile (4.8) determines the initial temperature profile that is used in the hydrodynamic code,
(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 freestreaming 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 freezeout 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 freestreaming 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 spacetime rapidities
(4.10) 
and
(4.11) 
Elementary kinematics [62] links the positions of a parton on the initial and final hypersurfaces and its fourmomentum
(4.12) 
where and are the parton‘s rapidity and transverse momentum:
(4.13) 
Thus the phasespace density of partons at the proper times and are related
(4.14)  
In our approach we restrict ourselves to the boostinvariant systems, hence it is reasonable to assume the following factorized form of the initial distribution of partons,
(4.15) 
where is the transverse density of partons obtained again from the GLISSANDO model, namely
(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