{'fraction_non_alphanumeric': 0.04024327460886713, 'fraction_numerical': 0.016500440651640882, 'mean_word_length': 4.48528215191691, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Collective phenomena in the early stages of relativistic heavy-ion collisions

A recently developed framework of highly-anisotropic and strongly-dissipative hydrodynamics -- ADHYDRO -- has been introduced and used to analyze the space-time evolution of matter produced in ultra-relativistic heavy-ion collisions. The main goal of this analysis was to study the effect of initial highly-anisotropic stages on the final soft hadronic observables typically measured in the experiment. The study was done in the context of the heavy-ion measurements performed at RHIC (Relativistic Heavy Ion Collider) in Brookhaven National Laboratory. Starting from the general assumption about the form of the phase-space distribution function, thermodynamic properties of locally anisotropic systems of particles have been studied and the form of the generalized equation of state has been formulated. The dynamic equations determining the evolution of a highly-anisotropic fluid have been introduced. The form of the entropy source related to the mechanisms leading to thermalization of the system has been defined. In the simplest case of purely-longitudinal and boost-invariant expansion, different features of the model have been analyzed. Using the ADHYDRO model in the general (3+1)D and boost-invariant (2+1)D versions, different possible scenarios of early stages of heavy-ion collisions have been analyzed. The results of the ADHYDRO model have been compared to the results obtained from the reference, perfect-fluid hydrodynamic model -- LHYQUID. The results of the hydrodynamic models were coupled to the statistical Monte-Carlo model THERMINATOR. The following conclusions have been drawn: a) All studied observables are almost insensitive to the initial anisotropic stage provided the initial conditions of the evolution are properly readjusted, b) Complete thermalization of matter may take place only at the times of about 1 fm/c. In this way the early thermalization puzzle may be circumvented.


Searching for new states of strongly interacting matter
Relativistic heavy-ion physics is an interdisciplinary research field connecting nuclear physics and elementary particle physics. It emerged in the middle of 1980s, when the relativistic beams of heavy ions became available at the Brookhaven National Laboratory and at CERN. The fast development of this field was triggered by a broad interest in the experimental and theoretical studies of properties of strongly interacting matter in the wide range of thermodynamic parameters such as temperature, T , and baryon chemical potential, µ B . In the framework of Quantum Chromodynamics (QCD), which is the fundamental non-abelian gauge field theory of strong interactions, one expects that the strongly interacting matter undergoes deconfinement and chiral phase transitions. The experimental verification of those phase transitions is the main aim of the experimental heavy-ion programs.
In the deconfinement phase transition, the hadronic matter changes into a system of interacting quarks and gluons. The initial theoretical concepts, based on the phenomenon of the asymptotic freedom (discovered by Gross, Politzer and Wilczek) and color screening, indicated that this system was weakly interacting -the asymptotic freedom predicts that the interaction strength between color charged partons (quarks and gluons) vanishes at short distances.
Starting from those ideas, Collins and Perry formulated a hypothesis of existence of a new quark state of matter in compressed super-dense star interiors. Also at that time, the limiting temperature introduced by Hagedorn was reinterpreted by Cabibbo and Parisi as the temperature defining the phase transition between hadronic and quark matter. In 1975 they drew the first phase diagram of strongly interacting matter. It indicated that at high temperature and/or large baryon density the hadronic matter should pass a phase transition to a new state of matter. This new phase was termed by Shuryak a quark-gluon plasma (QGP).
Besides the deconfinement phase transition, at high temperature and/or large baryon density one deals also with the restoration of chiral symmetry, i.e., with the chiral phase transition. The order parameter of the chiral phase transition is the quark condensate. The standard QCD vacuum is characterized by a non-vanishing value of the quark condensate, which is a consequence of spontaneous breaking of chiral symmetry at T = µ B = 0. If the system passes the chiral phase transition, the quark condensate drops practically to zero. It was quickly realized that the verification of theoretical predictions concerning properties of strongly interacting matter may be done in laboratory conditions by performing collisions of heavy atomic nuclei at ultra-relativistic energies. In The experimental results have significantly broadened our general knowledge of the properties of strongly interacting matter. First of all, the multi-particle systems produced in high-energy collisions of heavy ions turn out to be very much different from the systems obtained in more elementary hadron-hadron or electron-positron interactions. It is also generally accepted that due to extremely high densities and sufficiently large sizes of the systems produced in heavy-ion collisions it is possible to ensure the proper conditions for the formation of a quite long-living new state of strongly interacting matter.
Several interesting physical phenomena have been discovered. In particular, at RHIC and at the LHC one observes the strong collective behavior of matter (quantified by splitting of the slopes of the transverse-momentum spectra of various hadron species and by the large elliptic flow coefficient v 2 ), strangeness enhancement, and strong in-medium jet quenching. These facts indicate that quark and gluon degrees of freedom play the main role in the hot and dense medium. However, the data suggests also that the produced system is not weakly interacting. The great success of relativistic perfect-fluid hydrodynamics in the correct description of the data suggests short equilibration times, and consequently large cross sections for quark-gluon interactions. This gives hints that the quark-gluon plasma should be considered as a strongly interacting fluid rather than as a weakly interacting gas.
The RHIC experiments suggest the existence of a strongly interacting quarkgluon plasma in the high-temperature region of the QCD phase diagram (above the critical temperature and for negligible baryon chemical potential, see Fig. 1.1). In this region, the lattice QCD simulations predict a crossover rather than a genuine phase transition. In the crossover phase transition, the thermodynamic variables change suddenly in a narrow range of temperature, however, all changes are smooth and show no discontinuities characteristic for the real phase transitions. The aim of the other current (STAR, NA61) and planned (FAIR, NICA) experiments is to explore different parts of the QCD phase diagram. The low beam energy scan program at STAR and the NA61 experiment search for signatures of the QCD phase transition at finite values of the baryon chemical potential. They try to identify the position of the critical point that marks the end of a line describing the first order deconfinement phase transition at finite values of µ B . The existence of such a line is suggested by many QCD inspired theoretical models.

Theoretical tools for different stages of relativistic heavy-ion collisions
In order to interpret the multitude of experimental data and formulate conclusions about physical properties of matter, we have to compare the experimental results with theoretical model predictions. Since, in practice, it is impossible to obtain a complete description of such complicated processes within QCD only, we have to develop different effective models that are adequate for different stages of the evolution. Just before the collision, the coherent nuclear wave functions may be viewed (due to very high energies) as the gluon sheets described in the theory of color glass condensate (CGC). The passage of the nuclei through each other at the collision proper time identified with τ = 0 is followed by a subsequent pre-equilibrium stage which includes variety of processes leading to the formation of incoherent distributions of partons.
One argues that the complex configuration of color electric and magnetic fields named glasma is produced just after the collision. Its further dynamics is described by the classical Yang-Mills theory. Finally, the glasma fields decay by radiating gluons. Similar picture is established in the framework of string models where one considers collisions of hadrons, followed by the spanning of color strings or tubes. In this case, the strings decay producing partons or hadrons.
The initial dynamics of the collision may be also viewed as reciprocal penetration of partons described by their distribution functions, which leads to multiple parton scatterings. The subsequent non-equilibrium evolution of the system is then described in the framework of kinetic theory with collision terms treated in the perturbative limit. This approach is known as the parton cascade model (PCM). However, if the system is dense, the particles are off-shell and the use of the quantum transport theory becomes necessary (although this is not easy in practice).
The pre-equilibrium dynamics leads eventually to isotropization of momentum distributions of partons, and subsequently to local thermal equilibration of the system. The processes leading to thermalization are not understood well enough at the moment, thus the very fast thermalization, expected at RHIC and the LHC, remains still an unresolved puzzle. In this Thesis we address this problem and propose how to circumvent it. A more detailed discussion of the early thermalization puzzle is presented in Chapter 2.
Multiplicities of particles produced in relativistic heavy-ion collisions are very large. Thus, the use of considerations based on statistical physics seems to be well justified. If the system reaches the state of local thermal equilibrium, its further evolution may be described by relativistic hydrodynamics. During hydrodynamic evolution the local values of thermodynamic variables change gradually. These changes are connected with the changes of the underlying phase-space distribution functions. If the equilibration rates in the system are not sufficiently fast, the local thermal equilibrium is not maintained and the hydrodynamic picture fails. In this case, the system's evolution should be described in the framework of the kinetic theory.
In a series of publications, it has been shown that the space-time evolution of matter produced at RHIC may be well described by relativistic hydrodynamics of the perfect fluid [1,2,3,4,5,6,7,8,9,10,11,12,13] or (if the deviations from equilibrium are noticeable but small) by dissipative hydrodynamics with a small viscosity to entropy ratio [14,15]. Relativistic hydrodynamics requires the information about the form of equation of state (EOS). The latter can be obtained from the lattice QCD simulations (LQCD). The present results obtained at µ B = 0 predict a smooth crossover phase transition from hadronic matter to the quark-gluon plasma at the critical temperature T c ≈ 170 MeV.
The hydrodynamic expansion of the fireball leads to the fast cooling of matter. The initially hot and dense quark-gluon plasma passes the phase transition back to hadronic matter. The hadronic interactions are less and less effective. At a certain stage, the interactions cease and particles freely stream to detectors. This process is called the kinetic freeze-out. It is expected that the kinetic freeze-out is preceded by the chemical freeze-out -the stage where hadronic abundances are fixed. Nevertheless, in many approaches the identification of the two freeze-outs (the single freeze-out model) turns out to be a reasonable approximation. In this case, all hadronic observables may be obtained with the help of the Cooper-Frye formula [16]. In particular, this approach allows for the successful description of the HBT radii consistently with other soft hadronic observables such as the transverse-momentum spectra and the elliptic flow.
If different hadron species undergo (kinetic) freeze-out at different temperatures, the modeling of freeze-out becomes a much more complex problem. In this context, there is a growing interest in developing hybrid transport approaches, i.e., the models that link hydrodynamics with Monte-Carlo simulations which describe the scattering between hadrons in the final state, for example, this is achieved with ultra-relativistic quantum molecular dynamics (UrQMD) [17,18].

Objectives of the Thesis
In this Thesis we analyze the space-time evolution of matter produced in relativistic heavy-ion collisions. We propose to describe this evolution in the framework of a new, hydrodynamics-like model that inherently includes the possibility of having high anisotropy of pressure at the early stage. We call this new framework the ADHYDRO model (highly-Anisotropic and strongly-Dissipative HYDROdynamics). Using ADHYDRO we try to reproduce the soft hadronic observables measured at RHIC. The presence of the anisotropic phase in our description, for about 1 fm/c at the beginning of the collision, helps us to circumvent the problem of early thermalization -in order to describe the data we do not have to assume that matter reaches local thermal equilibrium within a fraction of a fermi. On the other hand, a possible connection of ADHYDRO to the underlying kinetic theory suggests that the microscopic collision times are of the order of a few hundredths of a fermi, a result suggesting that our picture may be consistent with the notion of a strongly interacting quark-gluon plasma.
The Thesis is organized as follows. In Chapter 2 we address the problem of fast thermalization in the very early stages of relativistic heavy-ion collisions. We also motivate the development of the ADHYDRO model presented in the Thesis. In Chapter 3 we discuss thermodynamic properties of the locally anisotropic systems of particles and formulate generalized equation of state. In Chapter 4 we derive the hydrodynamic equations determining the evolution of matter in ADHYDRO. Subsequently, in Chapter 5 we discuss different aspects of the model in the simple case of purely longitudinal and boost-invariant evolution. The Chapter 6 contains the information about the initial conditions for realistic hydrodynamic simulations. In Chapter 7 we include the freeze-out prescription. The results of the model for the boost-invariant systems are presented in Chapter 8 and the results for the most general non-boost-invariant case are shown in Chapter 9. We close the Thesis with the Summary and Conclusions.
In this Chapter, the early thermalization puzzle in relativistic heavy-ion collisions is discussed. Different theoretical approaches describing early stages of the collisions are reviewed briefly and different results for the thermalization time-scales are presented. Finally, several recent attempts to circumvent the problem of early thermalization are introduced. In particular, we provide motivation for the development of our framework of highly-anisotropic and strongly-dissipative hydrodynamics (ADHYDRO) which will be presented in greater detail in the next Chapters.

Evidence for early thermalization
The results collected in the heavy-ion experiments at RHIC are interpreted as the evidence that matter produced in relativistic heavy-ion collisions reaches local thermal equilibrium in a very short time 1 . The main argument for such early equilibration comes from the observation that the space-time evolution of matter is very well explained by the relativistic hydrodynamics of perfect fluid [1,2,3,4,5,6,7,8,9,10,11,12,13] or by viscous hydrodynamics with a small viscosity to entropy ratio [14,15]. In particular, a very good description of the soft region of the transversemomentum spectra of hadrons has been achieved in those frameworks. The shape of the spectra as well as their azimuthal dependence are interpreted as the evidence for the radial and elliptic flow of matter (the latter is described by the elliptic flow parameter v 2 ). The key ingredient of the successful hydrodynamic calculations is an early starting time assumed for the hydrodynamic evolution, τ 0 < 1 fm/c. Since the initial starting time τ 0 is identified typically with the thermalization time τ th , the question arises whether such incredibly fast thermalization of matter can be explained by microscopic models of the very early stages of heavy-ion collisions.
Fast equilibration and perfect-fluidity are naturally explained by the concept that the produced matter is a strongly coupled quark-gluon plasma (sQGP) [19,20]. Strong interactions between quarks and gluons might be responsible for the fast thermalization rate and small viscosity. Unfortunately, the large value of the coupling constant prohibits the use of the standard perturbative methods in this case. That is why the methods based on the AdS/CFT correspondence (Anti de Sitter/Conformal Field Theory correspondence) have become now a very attractive tool to analyze the early dynamics of heavy-ion collisions. Of course, we have also at our disposal LQCD which allows us to study thermodynamic properties of the QCD matter [21,22]. In fact, deviations of the lattice results from the Stefan-Boltzmann limit in the temperature range T ∼ (1 − 3)T c (where T c ∼ 170 MeV is the critical temperature) may be interpreted as the evidence for a strong coupling.
On the other hand, we should keep in mind that the physics of relativistic heavyion collisions has been triggered by the idea that at extremely high temperatures, T ≫ T c , the asymptotic freedom property of QCD leads to color screening and deconfinement [23]. In this case, QGP should behave as a weakly interacting gas of quasiparticles (wQGP) treated with perturbative QCD (pQCD) techniques. It seems hard to find reasonable explanations for sufficiently fast equilibration of such a weakly interacting system. However, several approaches have been proposed in order to understand the early time behavior of wQGP.
The thermalization of wQGP was broadly studied in the framework of Monte-Carlo parton cascade models (PCM) where the equilibration of the system was understood typically as the effect of multiple hard and semi-hard rescatterings of deconfined partons. Originally, PCM included only lowest-order binary (2 → 2) collisions [24]. They turned out to be insufficient to provide fast equilibration. Moreover, the local equilibration of the system was observed only in the central region in coordinate space. Further improvements of PCM models led to inclusion of more subtle effects such as: dilated emission from excited partons, the balance between radiative emissions and reverse absorption processes, and the interference of soft gluons [25,26,27,28]. They have been proved to significantly affect the thermalization rate. Nevertheless, the estimates for equilibration time-scales have delivered quite long times, for example, about 1.8 fm/c for Au+Au collisions at the RHIC top energy. Recently, the importance of multi-particle processes have been emphasized. The particle production and absorption processes (2 → 3) [29,30] together with the three-particle collisions (3 → 3) [31] have turned out to be significant in speeding up the equilibration. It has been demonstrated that the overall kinetic equilibration may be achieved much earlier but still later than 1 fm/c. At the very early stages of heavy-ion collisions, the momentum distribution of produced partons is expected to be highly anisotropic. Using theoretical methods known from the electron-ion plasma physics in a framework of the quark-gluon transport theory, one can show that the anisotropies of QGP lead to kinetic instabilities with respect to the formation of unstable plasma modes [32,33]. For QGP, only magnetic instabilities known as filamentation or chromo-Weibel instabilities [34] are relevant. These strong instabilities drive the system towards the locally isotropic state much faster [35,36,37] than the ordinary perturbative scattering processes in the weak coupling regime. This, in turn, indirectly helps to achieve the genuine equilibration 2 .
A quite new and original mechanism of thermalization of QGP has been formulated during studies of statistical mechanics of the Yang-Mills classical mechanics. It has been postulated that equilibration may be understood as the effect of chaotic dynamics of the non-Abelian highly excited classical color fields coupled or not to the classical colored particles [38,39,40]. In this case, the thermalization time estimated from the Lyapunov exponent is expected to be in the range 0.4 − 1.0 fm/c.
The microscopic models of the early stages of relativistic heavy-ion collisions use quite often the concepts of color strings or color-flux tubes. The idea of fast thermalization is particularly difficult to reconcile with the results of such models. For example, in the color glass condensate (CGC) approach that is often used to describe initial stages of relativistic collisions of heavy nuclei [41], the momentum distribution functions are far away from the equilibrium ones [42,43].
On the basis of the QCD saturation mechanism incorporated in the CGC model, the bottom-up thermalization scenario [44] was proposed. Since for very large nuclei and/or sufficiently high collision energies the saturation scale Q s is much larger than the QCD scale Λ QCD (at RHIC one expects Q s ∼ 1 GeV), the weak coupling techniques can be used in the calculations [45]. At the very early times, τ ∼ Q −1 s , small x gluons with the typical transverse momentum of the saturation scale Q s are freed from the incoming nuclei. Subsequently, in the time interval Q −1 s < τ < α −5/2 Q −1 s , the soft gluons are emitted. After the time τ ∼ α −5/2 Q −1 s the soft gluons form a thermal bath and draw the energy from hard gluons. The stop of temperature growth identified with full thermalization happens at τ ∼ α −13/5 Q −1 s . At that time the primary hard gluons loose their all energy. Straightforward calculations done in this approach for the most central collisions at RHIC give the equilibration time of 2.6 fm/c or even larger [41,46].
The CGC model itself gives quite interesting picture of the temporal behavior of the momentum anisotropies of the matter produced in heavy-ion collisions. For the very early proper times, τ ≪ Q −1 s , the dominant gluon fields produced in the collisions are classical and can be described by the classical Yang-Mills equations. In this case, one predicts the energy-momentum tensor with negative longitudinal pressure [47,48]. The pressure asymmetry persists also at later times, τ ≫ Q −1 s , where the analytic perturbative approaches [49] and the full numerical simulations [50] yield the energy-momentum tensor with the longitudinal pressure significantly lower than the transverse pressure. This results indicates that the hydrodynamic description of the gluon system as a perfect fluid may start only at later times .
In the color-flux-tube models, the production of partons is described by the Schwinger tunneling mechanism in strong color fields [51,52]. The Schwinger mechanism leads to transverse-momentum distributions which have Gaussian shapes.
However, if the string tension fluctuates, the Gaussian distributions are changed effectively into the exponential distributions [53]. The resulting forms of the spectra resemble the Boltzmann distributions in the transverse-momentum space. Thus, the produced system may be viewed as a system where only transverse degrees of freedom are thermalized. This phenomenon may be considered as apparent transverse thermalization, since thermal-like spectra are caused by the specific mechanism of particle production and not by the parton rescattering processes.
An intriguing thermalization scenario explaining the equilibration of both the transverse and longitudinal degrees of freedom is based on the Hawking-Unruh effect [54,55] that states that an observer moving with an acceleration a experiences the influence of a thermal bath with the effective temperature a/2π. In the Color Glass Condensate theory the typical acceleration in the strong color fields is of the order of Q s , which leads to the thermalization time 2π/Q s ∼ 1 fm/c.
We conclude our discussion of the early thermalization problem with the statement that microscopic models of heavy-ion collisions cannot explain naturally the thermalization time-scales shorter than 1 fm/c. Moreover, several models suggest that thermalization of the transverse degrees of freedom takes place before the complete three-dimensional thermalization is achieved 3 . As a consequence, these models predict significant momentum anisotropies during early stages of the collisions. It means that the standard hydrodynamic picture may be inadequate during very early times of heavy-ion collisions.
These observations have initiated development of effective, hydrodynamics-like models which describe the whole evolution of matter and relax the assumption of the early thermalization. Typically, this is done by imposing existence of a short initial pre-equilibrium phase followed by more standard hydrodynamic evolution. These models are discussed below in Section 2.2.

Relaxing the early thermalization assumption
Significant differences between the results of various theoretical models of the early stages of heavy-ion collisions reflect our lack of precise knowledge concerning the processes of particle production. The early stages are, however, of great importance as they may have an important impact on further evolution of matter. Based on the qualitative hints coming from the microscopic models discussed above in Section 2.1, we may study different effective scenarios for the early stages. By comparing the experimental consequences of such scenarios, we may try to check which hypothesis is the most likely to be realized in Nature.
Of the special interest are the physical scenarios where the assumption of the early thermalization of matter is relaxed. In practice, this means that one does not start the hydrodynamic description at an early time which is a fraction of a fermi but the early stage is described by an effective non-equilibrium model which is matched with the hydrodynamic description at a later time.
In particular, it has been proposed that the early-stage dynamics consists of the partonic free-streaming (FS) phase initiated at τ 0 ≪ 1 fm/c and followed by a sudden equilibration (SE) process. Typically, the sudden thermalization is realized by applying the Landau matching conditions at a certain transition time τ tr > τ 0 . The remaining part of the evolution is determined by standard relativistic hydrodynamics. The first FS+SE approximation was proposed by Kolb, Sollfrank, and Heinz [56]. Recently, it has been reconsidered by several authors [57,58,59]. It has been demonstrated that the final observables are insensitive to the presence of the early free-streaming phase if this phase lasts no longer than 1 fm/c [60].
According to the results coming from the string and color-flux-tube models, the pre-equilibrium stage may be viewed in a different way, namely, as the transversally thermalized parton system. It is conceivable that the initial quark-gluon matter is produced already in a state close to the thermodynamic equilibrium in the transverse direction only, while its longitudinal dynamics is described as free streaming of two dimensional clusters [61]. This observation resulted in the development of the model that is quite similar to the FS+SE framework. In this approach, the initial dynamics of the parton system is described by transverse hydrodynamics (TH) [62]. Similarly to the FS+SE approach, the system that is initially not completely thermalized undergoes sudden equilibration and its subsequent evolution is described by standard hydrodynamics. Interestingly, also in this framework the insensitivity of final results with respect to details of the initial transverse-hydrodynamics stage has been observed.
The two approaches discussed above, which have turned out to be successful in describing the experimental data, show that the models incorporating high initial anisotropy should also account for the fact that matter is eventually fully thermalized. Therefore, there is a need for models which can, in a comprehensive way, describe the gradual transition between early highly-anisotropic stages and later fully thermalized stages governed by the perfect-fluid hydrodynamics. Of course, in this context the use of the Landau matching conditions seems to be a very crude approximation.
Certainly, for small deviations from equilibrium the most appropriate description should be based on the framework of the second order viscous hydrodynamics [63,64,65]. However, at the earliest stages of the collision the anisotropies are very large. In this regime neither perfect-fluid hydrodynamics nor viscous hydrodynamics are formally applicable. In particular, viscous corrections become compatible with the leading terms, which may lead to negative pressure. Of course, these stages can be described by models based on the kinetic theory, but the latter has its own well known limitations.
The idea of developing the hydrodynamics-like approach that is suitable for uniform description of several stages of evolution of matter has motivated our development of the ADHYDRO model. This is an effective model which is able to describe strongly-dissipative dynamics of initially highly-anisotropic systems. ADHYDRO works in three different regimes. At very early stages, it describes the system with ex-tremely high momentum anisotropy. Close to equilibrium, where the anisotropies are small, ADHYDRO agrees with the viscous hydrodynamics 4 . For very large evolution time the system is (locally) isotropic in momentum space and our model gives standard perfect-fluid hydrodynamic description.

Locally anisotropic phase-space distribution functions
In this Thesis we consider locally anisotropic systems of particles. We assume that such systems are characterized by the anisotropy of their momentum distributions analyzed in the local rest frame. The local rest frame at the space-time point x is the reference frame in which a small element of the considered system, located at this space-time point, is at rest. Locally anisotropic systems of particles may be described by the phase-space distribution functions whose dependence on different momentum components in the local rest frame is different. We restrict our considerations to the case where the anisotropy of the system is connected with the difference of pressures in the longitudinal and transverse directions. These directions are always defined with respect to the collision axis 1 . The asymmetry in momentum may be realized by introducing two different spacetime dependent scales, λ ⊥ and λ , that may be interpreted as the transverse and longitudinal temperatures, respectively. The explicitly covariant form of such a nonequilibrium distribution function is where p µ denotes the four-momentum. The four-vector U µ describes the fluid fourvelocity In comparison with standard equilibrium distributions, we introduce here a new four-vector V µ . Its appearance is related to the special role played by the beam axis (z-axis) and it is defined in the following way The four-vectors U µ and V µ satisfy simple normalization conditions In the local rest frame (LRF) of the fluid element the four-vectors U µ and V µ have simple forms 5) and the formula (3.1) reduces to For massless particles m ⊥ = p ⊥ , and we may also write Here, for our convenience, we have introduced two dimensionless parameters For example, we may consider the Boltzmann distribution that is stretched or squeezed in momentum space, where ξ is expressed by the pressure anisotropy parameter x, (3.10) For isotropic systems x = 1 (ξ = 0) and (3.11)

Moments of anisotropic distributions
The knowledge of the distribution function (3.1) allows us to calculate one-particle characteristics of the system under consideration. In particular, the first and second moments of the distribution function yield the particle density, energy density, and the two pressures. Since at the very beginning of the evolution the considered system is expected to consist mainly of gluons, in this Section we restrict our analysis to massless particles by setting m = 0.
Using the standard definitions of the particle number current N µ and the energymomentum tensor T µν , and (3.13) and the definition of the entropy current where we obtain the following decompositions (the method leading to Eqs. (3.16)-(3.18) is explained in Ref. [67]): and Here g 0 is the degeneracy factor related to internal quantum numbers of particles that form the fluid and p 0 = E p = m 2 + p 2 ⊥ + p 2 . For simplicity, in Eq. (3.15) we have introduced ǫ = +1 for bosons and ǫ = −1 for fermions. In the limit ǫ → 0 we obtain the classical Boltzmann definition.
In the isotropic case, the longitudinal and transverse pressures are equal, P ⊥ = P = P eq , and we recover the standard form of the energy-momentum tensor valid for the perfect fluid. In general, when the system is out of equilibrium, the last term in (3.17) is different from zero. Then, in the local rest frame of the fluid element we find 3 Hence, as expected, the structure (3.17) allows for different pressures in the longitudinal and transverse directions. For massless particles the condition that the energy-momentum tensor is traceless leads us to the following expression From Eq. (3.20) we can infer that thermodynamics-like parameters describing the anisotropic system are functions of two independent parameters. For example, the energy density is a function of two pressures, ε = ε(P ⊥ , P ), or equivalently of two temperatures, ε = ε(λ ⊥ , λ ). This feature distinguishes anisotropic systems from equilibrated matter where the energy density for baryon-free matter is a function of isotropic pressure only, ε = ε(P eq ). For similar out-of-equilibrium systems with non-vanishing baryon number, three parameters are needed to specify the local state of the system.

Pressure anisotropy
The internal consistency of the approach based on the distribution function (3.1) allows us to derive the following concise structure of the thermodynamic parameters describing the anisotropic system of particles (for more details see [67]) where the function R(x) is defined by the integral The range of the integration over ξ ⊥ and ξ is always between 0 and infinity. The function R ′ (x) in Eqs. (3.22) and (3.23) is the derivative of R(x) with respect to x.
The pressure relaxation function R(x) describes bulk properties of the system characterized by the anisotropy x. The particular form of R(x) follows directly from the microscopic calculation if the distribution function (3.1) is known. At the more phenomenological level of description, a special form of the function R(x) can be assumed and treated as an ansatz describing the anisotropic system [69]. The anisotropy parameter x and the non-equilibrium entropy density σ are the most convenient independent variables used in our further analysis.
behaves like x −1/3 . In this limit P ⊥ = 0 and ε = P . Similarly, for x ≫ 1 (λ ⊥ ≫ λ ) the function R(x) behaves like x 1/6 , implying that P = 0 and ε = 2P ⊥ . This behavior is expected if we interpret the parameters λ ⊥ and λ as the transverse and longitudinal temperatures, respectively. Since color fields are neglected in this approach both pressures are expected to be positive, P ⊥ ≥ 0 and P ≥ 0. Furthermore, close to equilibrium the longitudinal and transverse pressures are expected to equilibrate, P ⊥ = P . Thus, according to Eqs. (3.22) and (3.23) we obtain another key restriction for R(x), namely, Based on the RHIC experimental data, we assume that the system eventually achieves the perfect fluid limit, hence we may write that lim τ →∞ x = 1.
The particle and entropy densities introduced in Eqs. (3.16) and (3.18) are defined by the following integrals (3.28) Comparison of Eqs. (3.27) and (3.28) indicates that the particle density and the entropy density are proportional to each other, with the proportionality constant depending on the specific choice of the parton distribution function f .

Special forms of anisotropic distribution functions
So far we have not considered any explicit forms of the distribution functions. In this Section, we turn to the discussion of the anisotropic distribution functions satisfying the ansatz (3.1) and having the following form Equation (3.29) may be regarded as a simple generalization of the equilibrium distributions in the case of vanishing baryon number. In a more general case, we should take into account also the net baryon number density n B and introduce the baryon chemical potential µ B in (3.29). However, if we focus on the central rapidity region in the ultra-relativistic limit (which is the case when analyzing the RHIC and LHC experimental data) we may neglect the baryon chemical potential µ B . Hence, the thermodynamic parameters remain the functions of entropy density and anisotropy. The function f (ǫ = 0) may be considered as the Boltzmann equilibrium distribution that has been stretched (or squeezed) in the longitudinal direction. For completeness, it is also interesting to consider the stretched (or squeezed) Bose-Einstein With the definition (3.29), we recover the structure (3.21)-(3.23) [68]. In the case of the classical Boltzmann statistics, the relaxation function (3.24) is given by the form 4 Equations (3.27) and (3.28) relate the entropy and particle densities in the following way σ = 4n. (3.32) In the similar way we may consider the Bose-Einstein and Fermi-Dirac distributions, see summary in Table 3.1 .
The results presented in Table 3.1 may be used in Eqs. (3.21)-(3.23). In this way we find that the numerical factor ε/ε BE (giving the ratio of the classical energy density to the Bose-Einstein energy density for the fixed σ and x) is very close to unity, ε/ε BE ≈ 1.027 . It means that the equations defining the energy density and pressures in terms of σ and x for the classical and quantum systems are practically the same. Therefore, in our further studies we use the classical statistics given by Eq. (3.29). Moreover, in order to account for the internal degrees of freedom of gluons we set g 0 = 16 (spin × color).
We stress that physical restrictions concerning the pressure relaxation function, discussed in Section 3.3, exclude from the consideration certain forms of distribution functions that formally satisfy the ansatz (3.1). For example, the anisotropic 4 Note that for x < 1 the function (arctan Table 3.1: Thermodynamic properties of the anisotropic systems for different statistics. The rows contain the following information: R(x) -pressure relaxation function (3.24), g -constant defined by Eq. (3.25), σ(n) -relation between entropy density and particle density calculated from Eqs. (3.27) and (3.28).
distribution functions of the factorized form In the upper part of Fig. 3.1 we present the rescaled pressure relaxation function r(x) = π 2 R(x)/(3g 0 ) while in the lower part of Fig. 3.1 we have plotted its derivative r ′ (x). We can see that according to (3.26) r ′ (x) vanishes when P ⊥ = P . In Fig. 3.2 we show the ratios: P /P ⊥ (solid red line), P /ε (green dotted line), and P ⊥ /ε (blue dashed line) for the case of R(x) given by Eq. (3.30). The considered ratios are functions of x only. In agreement with previously given remarks we see that ε = P for x = 0 and ε = 2P ⊥ for x → ∞. The equilibration of pressures is reached for x = 1. To a good approximation P /P ⊥ ≈ x −3/4 , thus, the anisotropy parameter x may be treated as the direct measure of the pressure anisotropy.

Generalized equation of state
The equations of perfect-fluid hydrodynamics form a closed system of equations only if they are supplemented with the equation of state that specifies thermodynamic quantities as functions of one parameter 5 . Similarly, our framework requires that all thermodynamic quantities should be expressed by two arbitrarily chosen parameters, e.g., by the non-equilibrium entropy density σ and the anisotropy parameter x. Such relations play a role of the generalized equation of state and allow us to close the system of dynamic equations (the latter will be formulated in Chapter 4).
If we consider a system of partons described by the anisotropic phase-space distribution function given by Eq. (3.29) and obtained by squeezing (or stretching) the classical Boltzmann distribution, we have to use Eqs. The derivative r ′ (x) of the pressure relaxation function (solid red line).
Equation of state (3.33) is very likely to be realized at the very early stages of the collisions, when the produced system is highly anisotropic and described by the distribution functions discussed in this Chapter. On the other hand, as the system expands, the interactions between its constituents thermalize it, so the system's thermodynamic variables approach the realistic EOS and x tends to unity. The realistic EOS for the strongly interacting matter, encoded in the temperature dependent sound velocity, is shown in Fig. 3.3. At low temperatures it is consistent with the hadron-gas model with all experimentally identified resonances. On the other hand, at high temperatures it coincides with the LQCD results.
The two limiting cases discussed above suggest that one may use the generalized equation of state of the form Here, the equation of state of the ideal relativistic gas, defined by the functions ε id (σ) and P id (σ) has been replaced by the formulas characterizing the realistic equation of state for vanishing baryon chemical potential, i.e., by the functions ε qgp (σ) and P qgp (σ) corresponding to the sound velocity function depicted in Fig. 3.3.
We have to stress that there is no microscopic explanation for the form of Eqs. (3.34). The main reason for this is that we cannot provide simple quasi-particle picture for the realistic QCD equation of state. Despite this Eqs. (3.34) have many attractive features. Firstly, as discussed above, they reproduce the two attractive limits and they interpolate between them in the continuous way. Moreover there is no additional assumptions needed (such as, e.g., the Landau matching conditions) in order to describe the change from one thermodynamic regime to another one. Additionally, Eqs. of matter in heavy-ion collisions in a uniform way.

Highly-anisotropic hydrodynamics with strong dissipation
In this Chapter we introduce the dynamic equations that determine the evolution of a highly-anisotropic fluid. Our approach is based on the energy-momentum conservation law and on the postulate that specifies the entropy production in the system. The resulting equations define the ADHYDRO model.

Space-time dynamics of highly-anisotropic systems
In Chapter 3 we argued that the system formed at the very early stages of relativistic heavy-ion collision may be described by the anisotropic phase-space distribution function (3.1). We assume now that the space-time dynamics of such a system is governed by the following equations, Equation (4.1) expresses the energy and momentum conservation law. It has an analogous form to the equation that defines dynamics of the perfect fluid. On the other hand, Eq. (4.2) describes the entropy production determined by the entropy source Σ. In the standard perfect-fluid approach, the right-hand-side of Eq. (4.2) vanishes and Eq. (4.2) becomes the condition of adiabaticity of the flow. We note that in the perfect-fluid hydrodynamics Eq. (4.2) follows directly from Eq. (4.1) and the appropriate form of the energy-momentum tensor. In our approach, Eq. (4.2) should be treated as an ansatz. In Chapter 3 we showed that the energy-momentum tensor T µν and the entropy flux S µ appearing in Eqs. (4.1)-(4.2) have the structure defined by Eqs. (3.17) and (3.18). In this case, the specific form of the entropy source Σ allows us to close the system of equations and determines the dynamics of the anisotropic fluid 1 .

CHAPTER 4. ADHYDRO MODEL
The projections of Eq. (4.1) on the four-vectors U ν and V ν yield

General three-dimensional expansion
In the general case, where matter expands in the longitudinal and transverse directions without any symmetry constraints, we may use the following parametrization of the four-velocity of the fluid U µ and the four-vector V µ In Eqs. (4.5) and (4.6) we have introduced the two components of the transverse four-velocity of the fluid, u x and u y , and the longitudinal fluid rapidity ϑ. Using the normalization condition (3.4) we also find If the parametrizations (4.5) and (4.6) are used, the differential operators appearing in (4.3) and (4.4) take the following form where the two linear differential operators L 1 and L 2 are defined by the expressions 14) In Eqs. (4.14)-(4.15) we have introduced the space-time rapidity η and the proper time τ , The difference between the longitudinal fluid rapidity ϑ and the space-time rapidity η is defined as We use the boldface notation for two-dimensional vectors in the transverse plane: Besides Eqs. (4.3) and (4.4), the two additional equations describing transverse dynamics of the system are needed. They can be chosen, for example, as the two linear combinations: In this case we find Here D = U µ ∂ µ is the total time derivative and θ = ∂ µ U µ is the volume expansion rate. Using this notation, the entropy production equation, see Eq. (4.2), can be put in the compact form Dσ + σθ = Σ. If the generalized equation of state ε(σ, x), see Sect. 3.5, and the entropy source Σ(σ, x), see Sect. 4.3, are specified, the ADHYDRO equations form a closed system of five equations for five unknown functions: two components of the fluid velocity u x and u y , the longitudinal rapidity of the fluid ϑ, the non-equilibrium entropy density σ and the anisotropy parameter x. These functions depend on transverse coordinates x, y, the space-time rapidity η, and the proper time τ . One has to solve the ADHYDRO equations numerically for a given initial condition specified at a certain initial time τ 0 (to be discussed in Chapter 6).

Entropy source
The entropy source Σ describes the entropy growth due to equilibration of pressures in the system. As we have already mentioned above, Eqs. (4.1) and (4.2) may be solved only if the function Σ(σ, x) is specified. The functional form Σ(σ, x) must be delivered as the external input for the anisotropic hydrodynamics.
The simplest ansatz for Σ(σ, x) satisfying general physical assumptions may be proposed in the following form [68,72,73] where τ eq is a time-scale parameter introduced to control the rate of the processes leading to equilibration of the system. The expression on the right-hand-side of Eq. (4.22) has several appealing features. First of all it is positive, as expected on the grounds of the second law of thermodynamics. Moreover, it has a correct dimension, since Σ is proportional to the entropy density σ. In this way Eq. (4. 22) does not destroy the scale invariance of the perfect fluid hydrodynamics, which allows for multiplication of σ in the evolution equations by an arbitrary constant.
In the natural way, the entropy source Σ defined by (4.22) vanishes in equilibrium, where x = 1, thus Eqs. (4.1)-(4.2) can be reduced to the equations of perfect fluid hydrodynamics in the limit x → 1. The expression (4.22) is also symmetric with respect to the interchange of λ ⊥ with λ (consequently Σ does not change if The quadratic dependence of the entropy source displayed in (4.23) is characteristic for the 2nd order viscous hydrodynamics (to be discussed in Chapter 5.3).

Anisotropy evolution
If both the generalized equation of state and the entropy production term are defined, one can derive a compact expression for the time evolution of the anisotropy parameter x. Of course, this expression should be considered together with other dynamic equations. Nevertheless, its form turns out to be useful in the analysis of the time evolution of the system. In the case described by Eqs. (3.34) we use Eq. (4.3) and find where Similarly, in the case described by Eqs. (3.33) we find Equation (4.26) may be treated as the special case of (4.24) where we set ε qgp = 3P qgp .

Purely longitudinal boost-invariant motion
In this Chapter we study a simple case of purely longitudinal, boost-invariant expansion. This situation is easy to investigate and we can obtain useful information about the behavior of anisotropic systems by applying simple analytic methods. This will help us to understand the behavior of more complicated systems which either expand transversally or do not exhibit boost invariance.

Implementation of boost invariance
The measurements done by the BRAHMS Collaboration [74] at RHIC showed that the rapidity distributions of charged particles are almost constant in the midrapidity region (|y| ≤ 1). This observation suggests that all measured quantities, to a good approximation, are boost-invariant in this region. The boost-invariance implies that all thermodynamic variables are functions of the proper time τ and transverse coordinates only, while the longitudinal flow has the form Since t = τ cosh η and z = τ sinh η, Eq. (5.1) is equivalent to the simple condition i.e., Y = 0. With this assumption at hand, it is sufficient to find the solutions of our model for z = 0 (τ = t), since solutions for other values of z may be obtained by simple Lorentz transformations.
To illustrate the main features of our model with simple examples, we neglect now the transverse expansion of the fluid. This is a reasonable approximation for the early stages of the evolution, where the longitudinal expansion dominates, u x = u y = 0 (u 0 = 1). We also assume now that the system may be treated as homogeneous in the transverse plane, hence, we neglect dependence of thermodynamic quantities on the transverse variables x and y. These assumptions are characteristic for the Bjorken model [75], where the behavior of the perfect fluid has been considered. In our approach the perfect-fluid dynamics is replaced by the hydrodynamics of the highly-anisotropic fluid.

Asymptotic limits
As we have discussed above in Sect which with the help of Eq. (4.25) may be also rewritten as First of all, it is useful to consider the isentropic flow where Σ(σ, x) = 0. In this case, from Eq. (5.4) we obtain the Bjorken solution where the parameters σ 0 and τ 0 are arbitrary constants. For Σ(σ, x) = 0 the righthand-side of Eq. (5.6) vanishes. This implies that either The case x = 1 corresponds to the standard, perfect-fluid hydrodynamics. Indeed, for an ideal relativistic gas in equilibrium we have P = P ⊥ = P eq and Eq. (5.3) implies ε(τ ) ∼ τ −4/3 . In the second case we obtain a simple longitudinal freestreaming solution (see our detailed discussion in Section 5.5). This situation is equivalent to the case where the equilibration time is set to infinity, τ eq = ∞. In addition we can consider also the case where P = 0. In this situation Eq. (5.3) describes the evolution of the energy density valid for transverse hydrodynamics [61,62], ε(τ ) ∼ τ −1 .

Consistency with Israel-Stewart theory
For a near-equilibrium system, the distribution function can be decomposed into an equilibrium part plus a small correction part. The methodology of including such corrections in the consistent way leads to the framework of relativistic viscous hydrodynamics of the first order (non-causal Navier-Stokes theory) or the second order (causal Israel-Stewart formalism [63,64,65]). When deviations from equilibrium are small, which in our approach may be quantified by the condition |ξ| ≪ 1, and the motion is purely longitudinal, one can show that our model is consistent with the Israel-Stewart theory.
As before, we consider again the baryon-free systems. The 2nd order dissipative fluid is characterized by the energy-momentum tensor where the pressure tensor P µν is defined by the formula Here ∆ µν = g µν − U µ U ν is the tensor projecting on the space perpendicular to the four-velocity U µ , P eq is the isotropic pressure, Π is the viscous bulk pressure, and π µν is the viscous shear tensor.
The energy-momentum tensor (3.17), which forms the basis of our approach, may be rewritten in the form (5.8) if we make the following identification Let us ignore now the bulk viscosity, since for massless particles or ultra-relativistic particles the bulk viscosity vanishes, Π = 0. If there is no transverse expansion, we may compare Eqs. (5.9) and (5.10) in the Landau local rest frame, and make the identification of the longitudinal and transverse pressures: P ⊥ = P eq + π/2, P = P eq − π, where P eq = ε eq /3 = ε/3 and π ≡ π 33 is the 33 component of the fluid shear tensor. Then, the expansion for small ξ gives [76,77] π ε eq = 2x We can match our model with Israel-Stewart theory by making the following identifications 1 τ eq = 4 15τ π , τ π = 5η π T σ eq , (5.14) where τ π is the shear relaxation time, η π is the shear viscosity, T is the temperature in the isotropic system, and σ eq is the equilibrium entropy density. In this way, Eq. (5.13) takes the form of the 2nd order viscous hydrodynamic equation for π The entropy production in the Israel-Stewart theory is given by the formula and for the longitudinal expansion considered in this Section, π αβ π αβ = 3π 2 /2. Substituting Eqs. (5.11) and (5.14) into Eq. (5.16) we obtain which is consistent with our main ansatz for the entropy production if deviations from equilibrium are small. One should stress that deviations from equilibrium in the Israel-Stewart theory have to be small, so that the dissipative fluxes are small compared to equilibrium variables, √ π µν π µν ≪ P eq . If this condition is not fulfilled, the use of the Israel-Stewart framework cannot be formally justified.

Deriving entropy production from kinetic theory
Our formulation of the highly-anisotropic hydrodynamics in one-dimensional case was followed independently by Martinez and Strickland [77,78], who derived equations analogous to Eqs. (5.3)-(5.4) starting from the relativistic Boltzmann equation Here C[f (x, p)] is the collision kernel which accounts for all microscopic interaction processes among particles. In the Martinez-Strickland model [77,78], the evolution of matter is described by two coupled ordinary differential equations: the first one for the momentum anisotropy ξ, and the second one for the typical hard momentum scale p hard (average momentum in the parton distribution function). These equations are obtained by taking the 0th and 1st moments of the Boltzmann equation (5.18).
In general, the entropy production may be calculated directly from the definition (3.14) with the help of Eq. (5.18). For classical particles, it is given by the expression Martinez and Strickland used the relaxation time approximation for the collision kernel. This allows for a direct calculation of the right-hand-side of Eq. (5.19) and gives the entropy source Σ M S (analogous to our Σ), and R MS (ξ(x)) is an analog of our R(x), R MS (ξ) = 1 2

Longitudinal free-streaming
In the case of anisotropic distribution functions (see Section 3.1) it is interesting to consider the pure longitudinal free transport (alias free-streaming) of particles, which is the regime in which particles do not interact completely with each other. In this case, the distribution functions should satisfy the collisionless Boltzmann kinetic equation We assume that the effects of free-streaming do not change the shape of the anisotropic distribution function (3.1) and the time evolution of the parameters λ ⊥ , λ , and U µ is determined by the conservation laws. We parametrize the four-momentum of a particle in the standard way where we have introduced the transverse mass m ⊥ = m 2 + p 2 ⊥ and the rapidity of a particle For simplicity, in the following analysis we assume that particles are massless, m = 0.
After simple algebra we obtain This equation is satisfied by any form of the function f provided λ ⊥ = λ 0 ⊥ and λ = λ 0 τ 0 /τ , where τ 0 , λ 0 ⊥ and λ 0 are constants. In the considered case the anisotropy parameter x has the form where x 0 = λ 0 ⊥ /λ 0 2 specifies the initial anisotropy of the system. Thus, we see that our approach with τ eq → ∞ includes the boost-invariant free-streaming as a special case.

Results for purely longitudinal boost-invariant motion
We solve Eqs. ii) x 0 = 0.032 2 , where the longitudinal pressure is much larger than the transverse pressure, in this case the momentum shape is prolate, and finally, iii) the initial system is isotropic, x 0 = 1. This case is the closest to the perfectfluid hydrodynamics. In all the cases the system's evolution is studied in the time interval 0.25 ≤ τ ≤ 10 fm.
In thin lines, respectively. The results are shown for different values of the time-scale parameter: τ eq = 0 (solid red lines), τ eq = 0.25 fm (dashed green lines), τ eq = 0.5 fm (dotted blue lines), τ eq = 1.0 fm (dashed-dotted brown lines), τ eq = ∞ (solid black lines). In the cases τ eq = 0.25, 0.5, 1.0 fm we observe fast changes of x(τ ) caused by the fact that H(x) and H MS (x) are very large for x ≫ 1 (H(x) ∼ 6x 3/2 ). For τ ≫ τ eq the function x(τ ) saturates and approaches unity. In addition, we show the result for τ eq = 0, which corresponds to sudden isotropization of the system. The initial highly-anisotropic system becomes isotropic almost instantaneously. Another interesting situation takes place for τ eq = ∞, where the quadratic behavior (5.29) typical for free streaming is observed. In the latter case, the anisotropy grows with time and the system never reaches equilibrium.
In the lower part of Fig. 5.1 we show the corresponding plots of the entropy density normalized by the Bjorken solution (5.7). During the equilibration process, large amount of entropy is produced. For τ eq = 0.25, 0.5, 1 fm, the entropy increases by about 65%, 80% and 100%, respectively. For τ ≫ τ eq the normalized entropy saturates, indicating that the flow attains the form of the Bjorken flow. In the case where τ eq = ∞, we also obtain the Bjorken-like dependence (5.7) where σ ∼ λ 2 ⊥ /τ and the transverse temperature becomes independent of time. This is an expected result, since without the longitudinal pressure no longitudinal work is done. If, in addition, there is no transverse expansion, the system cannot cool down and its temperature remains constant. In the case τ eq = 0 the sudden jump of entropy by about 50% is observed.
In Fig. 5.2 we show our results for the initially isotropic case iii). The notation is the same as in Fig. 5.1. In the upper part of Fig. 5.2, for τ eq = 0.25, 0.5, 1 fm we observe the non-monotonic behaviour of the functions x(τ ). This is caused by the fact that H(x = 0) = 0 and free streaming dominates the early evolution. In this case, the anisotropy grows with time (dx/dτ > 0). This is an expected feature, since even if the starting distribution is isotropic, the system becomes anisotropic if there are no interactions to sustain the equilibrium distribution. With increasing τ , the second therm on the right-hand-side of Eq. (5.5) starts to dominate over the first therm and the anisotropy decreases (dx/dτ < 0). For τ ≫ τ eq the system becomes isotropic and the term 2x/τ becomes negligible. In the lower part of Fig. 5.2 we observe the corresponding entropy production for τ eq = 0.25, 0.5, 1 fm. Again, the longer is the anisotropic stage, the more entropy is produced.
Finally, for completeness, in Fig. 5.3 we show our results for the case ii). For τ eq = 0.25, 0.5, 1.0 fm we again observe the non-monotonic behaviour of x(τ ). Initially, the dissipation tends to isotropize the system. However, as the equilibrium is reached the dissipative processes stop, H(0) = 0 and the free-streaming stage starts to dominate the dynamics. The subsequent evolution proceeds very much similarly to the case iii). In the lower part of Fig. 5.3, for τ eq = 0.25, 0.5, 1 fm, we observe the corresponding entropy production. We can see that when the system is near equilibrium the entropy production stops. Moreover, contrary to the case i), the higher is the value of τ eq , the less entropy is produced.
Comparing our results presented in Figs. 5.1, 5.2, and 5.3, we observe that the time dependence of x is very much similar in the cases (4.22) and (5.20). It suggests that the particular form of the entropy source has small effect on the time evolution. Far from equilibrium, the cases (4.22) and (5.20) differ significantly but the large entropy production leads to fast decrease of x. As soon as x becomes small these two schemes are equivalent, since for small x Eqs. We also note that τ eq = 0.25, 0.5, 1.0 fm correspond to τ MS eq ≈ 0.03, 0.067, 0.13 fm. Hence, our choice of the parameters corresponds to very short relaxation times used in the kinetic theory (in the way done by Martinez and Strickland). This may sound reasonable if we consider the plasma as a strongly interacting system.

Initial conditions
In order to make realistic comparisons with the data, we should include the transverse expansion of matter and, possibly, relax the assumption of boost-invariance. In this Chapter we introduce the initial conditions used in our model for such situations.

Glauber initial conditions for boost-invariant systems
In order to solve the ADHYDRO equations discussed in Chapter 4 in the boost-invariant 2+1 dimensional [(2+1)D] case, it is necessary to specify initial conditions at a given transverse space point x ⊥ for the initial proper time τ = τ 0 . The proper time τ = 0 corresponds to the moment when the centers of the colliding nuclei pass through the z = 0 plane. In the numerical calculations we choose τ 0 = 0.25 fm. Before this time, the parton interactions are expected to produce at least a partially thermalized system whose later behavior may be described by ADHYDRO. For boostinvariant simulations the initial conditions are defined by four functions: σ(τ 0 , x ⊥ ), x(τ 0 , x ⊥ ), u x (τ 0 , x ⊥ ), and u y (τ 0 , x ⊥ ). We assume that the initial energy density in the transverse plane, ε 0 (τ 0 , x ⊥ ), is proportional to the normalized density of sourcesρ(b, x ⊥ ), where b is the impact parameter corresponding to a given centrality class, The parameter ε i is the initial energy density at the center of the system in most central collisions. Its value may be estimated from the standard hydrodynamic calculations. We use the value ε i = 86.76 GeV/fm 3 . If the considered matter were in equilibrium, its temperature would be equal to T = 485 MeV. The normalized density of sources,ρ(b, x ⊥ ), is constructed as a linear combination of the wounded-nucleon density ρ W (b, x ⊥ ) and the density of binary collisions 3) The choice of the parameter κ (0 ≤ κ ≤ 1) controls the mixing between ρ W and ρ B . In the case κ = 0 we obtain the standard wounded-nucleon model [79]. The parameter κ should be chosen in such a way as to reproduce the experimentally observed dependence of the charged particle multiplicity on the collision centrality. Typically one uses κ = 0.14 for RHIC [80] and κ = 0.15 for the LHC [81]. The wounded-nucleon and binary-collision densities may be obtained from the optical limit of the Glauber model. In this approach, for a collision with the centrality specified by the impact vector b (see Fig. 6.1), we use the formula for the transverse density of wounded nucleons being a sum of the contributions from the nuclei with the mass numbers A (projectile) and B (target). Thus, assuming only symmetric collisions where A = B and large nuclei (A, B ≫ 1) we may write and the density of binary collisions is   The role of σ in is more important in the definition (6.4), since it influences the shape of profiles. In (6.5) it enters just as a multiplication factor. Thus, the eccentricity of ρ W (b, x ⊥ ) for RHIC is larger than that for the LHC, whereas for not too peripheral collisions the spatial eccentricity of ρ B (b, x ⊥ ) is in general the same, see The function ρ WS (x ⊥ , z) in (6.6) is the nuclear density function given by the Woods-Saxon spherical profile

GLAUBER INITIAL CONDITIONS FOR BOOST-INVARIANT SYSTEMS47
where R A is the radius of the nucleus with the mass number A calculated from the formula a is the surface diffuseness (a = 0.54 fm for both, gold and lead nuclei) and ρ 0 is the nuclear saturation density chosen in such a way that ρ WS (x ⊥ , z) is normalized to the number of nucleons in the nucleus. For the initial anisotropy parameter x(τ 0 , x ⊥ ), which we take as independent of x ⊥ and denote simply as x 0 , we consider again three different options: x 0 = 100, x 0 = 1, and x 0 = 0.032. The case x 0 = 1 is, of course, the closest to that described by the standard perfect-fluid hydrodynamics. The case x 0 > 1 corresponds to the initial situation where the transverse pressure is much larger than the longitudinal pressure. In this case the momentum shape is oblate (the momentum distribution is stretched in the transverse direction and squeezed in the longitudinal direction with respect to the beam axis). This type of the initial conditions is considered in the Color Glass Condensate (CGC) approach where the distribution functions in the longitudinal direction are described by the Dirac delta function, δ(p ) at z = 0 [42,43], see also Ref. [41]. The values x 0 < 1 correspond to the prolate momentum shape. This type of the initial conditions has been analyzed, for example, in Refs. [58,83].
By fixing both the initial energy density and the initial anisotropy parameter, we determine the initial entropy density profile from Eq. (3.34), namely Here ε −1 qgp (ε) is the inverse function to ε qgp (σ). The b-dependence displayed on the right-hand-side of Eq. (6.9) induces centrality dependence of the initial entropy density profiles.
We note that the entropy described by Eq. (6.9) corresponds to the entropy produced in the proper time interval 0 ≤ τ ≤ τ 0 , just before the ADHYDRO description is turned on. For τ > τ 0 , the ADHYDRO model accounts for the additional entropy production described by the entropy source term Σ. Since the overall entropy produced in a given class of collisions is connected with the final multiplicity, we may consider it as a fixed quantity. Therefore, if the initial anisotropy x 0 is much different from unity, the initial entropy density σ 0 is expected to be small. In this case, the large amount of entropy is produced in the phase described by ADHYDRO. On the other hand, if the initial anisotropy is close to unity, the subsequent entropy production should be negligible.
Our form of the initial conditions is completed by the assumption that there is no transverse flow present at τ = τ 0 , i.e., we set u x (τ 0 , x ⊥ ) = 0 and u y (τ 0 , x ⊥ ) = 0.

Tilted initial conditions for non-boost-invariant systems
The experimental data collected at RHIC show that the boost-invariance is realized only in a relative narrow central range of rapidities. Thus, the boost-invariant (2+1)D approach discussed in Section 6.1 is not sufficient to perform realistic predictions for observables outside the midrapidity region (|y| > 1). In order to investigate in more detail the influence of partial thermalization on the behavior of matter produced in the whole range of rapidities in relativistic heavy-ion collisions, we have to perform full-geometry non-boost-invariant 3+1 dimensional [(3+1)D] simulations based on the differential equations introduced in Chapter 4. For a general non-boost-invariant case, the initial conditions for ADHYDRO are defined by five functions: σ(τ 0 , η, x ⊥ ), x(τ 0 , η, x ⊥ ), u x (τ 0 , η, x ⊥ ), and u y (τ 0 , η, x ⊥ ), and ϑ(τ 0 , η, x ⊥ ). The initial entropy density profile is defined in the similar way as in the (2+1)D calculations, however, in the (3+1)D calculations we choose the parameter ε i separately for different physical situations (such as different initial anisotropy profiles and/or different values of the time-scale parameter τ eq ) in such a way as to reproduce well the final multiplicity.
The normalized density of sourcesρ(b, η, x ⊥ ) from (6.10), is constructed now as a sum of the two terms: Here, the function f (η) is the initial longitudinal profile The first term in (6.11) introduces contributions to the density of sources from the forward-(+) or backward-moving (−) participant nucleons. This formula assumes a preferred emission from the participating nucleons along the direction of their motion. The second term introduces a symmetric contribution from the binary collisions. The half-width of the central plateau, ∆η, and the half-width of Gaussian tails, σ η , are fitted to reproduce the RHIC multiplicity of charged particles in pseudorapidity.
The initial density of sources produced by a single forward-or backward-moving participant nucleon with rapidity y b = ln( √ s NN /m N ), where m N is the nucleon mass, 14) The range of rapidity correlations equals η m = y b − η s ≃ 3.36 and the shift in the rapidity η s = 2 is a phenomenological parameter. With the assumption (6.14), the hydrodynamic description outside the region (−η m , η m ) is not reliable. Similarly to (2+1)D calculations, for the initial anisotropy profile, x(τ 0 , η, x ⊥ ), we choose three different values: x 0 = 100, x 0 = 1, and x 0 = 0.032. Moreover we check the possibility of the spatial dependence of the initial anisotropy parameter, x 0 = x 0 (ρ(η, x ⊥ )). We assume also that there is no transverse flow present initially, therefore, we set u x (τ 0 , η, x ⊥ ) = 0 and u y (τ 0 , η, x ⊥ ) = 0, and the initial longitudinal rapidity of the fluid follows the simple Bjorken scaling ϑ(τ 0 , η, x ⊥ ) = η.

Freeze-out
In this Chapter we present first the method used in ADHYDRO to extract the freeze-out hypersurface from hydrodynamic evolution. Then, the Cooper-Frye prescription in the form used by the Monte-Carlo event generator THERMINATOR [84, 85] to describe particle production on a given freeze-out hypersurface is derived. Finally, calculations of physical observables with the help of THERMINATOR are explained in more detail.

Single-freeze-out approximation
The rapid expansion of dense and hot matter (described in our approach by ADHYDRO) implies that the density of particles decreases and the mean free path of particles grows. When the time between interactions of particles becomes larger than the expansion rate (described by divergence of the four-velocity field) τ coll ≥ τ exp , (7.1) the collision processes cannot maintain the local equilibrium any longer. The strongly interacting system of particles undergoes a transition to the system consisting of essentially free-streaming particles. This complex dynamic process, called freeze-out, is quite often treated approximately as if it took place locally on the three-dimensional freeze-out hypersurface. Typically, the freeze-out condition (7.1) is realized by assuming that hadrons completely decouple at the space-time points x where the entropy density drops below a certain freeze-out entropy density σ f . For systems in local thermal equilibrium described by the equation of state σ(T ) this condition may be formulated in terms of temperature, T (σ(x)) = T f = T (σ f ).

(7.2)
A set of freeze-out points x extracted from the hydrodynamic evolution by the constraint (7.2) forms the freeze-out hypersurface. One can subsequently use it to generate primordial particles employing the Cooper-Frye formalism.
In general, the inelastic processes in the expanding system cease earlier than the elastic ones. Therefore, one can distinguish chemical freeze-out (the stage when inelastic processes stop and hadron abundances are fixed) from subsequent kinetic freeze-out (the stage when all types of interactions cease and the momentum distributions of hadrons are frozen). Of course, the latter may be considered as the true freeze-out process. In the case of high energy collisions at RHIC, which are the subject of our studies, we can use quite efficient simplification called the singlefreeze-out approach [86,87,88]. In this approach, one assumes that the kinetic freeze-out coincides with the chemical freeze-out, T kin = T chem = T f . In our studies we adopt this scenario.

Cooper-Frye formalism
The momentum distribution of primordial hadrons emitted on the freeze-out hypersurface Σ may be calculated from the Cooper-Frye formula [89] dN where dΣ µ is the element of the three-dimensional freeze-out hypersurface submerged in the four-dimensional Minkowski space, f (x, p) is the phase-space distribution function of the emitted particles and g = 2s + 1 is the spin degeneracy factor, where s denotes the particle's spin. For systems in local thermal equilibrium described by the perfect fluid hydrodynamics one assumes equilibrium Fermi-Dirac (ǫ = −1) or Bose-Einstein (ǫ = +1) statistics, i.e., one uses the distribution functions depending on the scalar product of particle four-momentum p µ and the four-velocity U µ of the fluid element. Clearly, if deviations from equilibrium at freeze-out are large, the distributions (7.4) should be replaced by the appropriate non-equilibrium distributions. In our calculations, for a typical choice of parameters, only initial parts of freeze-out hypersurfaces exhibit large anisotropies. Thus, as long as they weakly contribute to total final multiplicities they may be simply neglected in the calculation of the observables. One can determine thermodynamic parameters on the freeze-out hypersurface with the help of thermal models [90,91,92,93,94] which analyze relative hadronic abundances. At RHIC, this type of study [94] indicates the non-zero values of chemical potentials. Thus in the generation of hadrons we introduce the chemical potential µ µ = Bµ B + I 3 µ I 3 + Sµ S + Cµ C , (7.5) where B, I 3 , S, C is the baryon number, third component of isospin, strangeness, and charm of the emitted hadron, respectively. Similarly, µ B , µ I 3 , µ S , and µ C denote the corresponding chemical potentials.
The particle's four-momentum may be parametrized in the standard way p µ = (m T cosh y, p T cos φ p , p T sin φ p , m T sinh y) , (7.6) where m T = m 2 + p 2 T is the transverse mass, p T is the transverse momentum, y is the (longitudinal) rapidity of the emitted particle, and φ p is the azimuthal angle of the transverse momentum. Together with the parametrization (4.5), one may write The form of dΣ µ may be obtained with the help of the formula known from differential geometry [95] dΣ µ = ε µαβγ ∂x α ∂a ∂x β ∂b ∂x γ ∂c dadbdc, (7.8) where ε µαβγ is the antisymmetric Levi-Civita tensor with ε 0123 = +1. The quantity dΣ µ (x) defines a four-vector that is perpendicular to the hypersurface at point x.
Its norm is equal to the volume of the hypersurface element. The variables a, b and c introduce a coordinate system in Minkowski space parameterizing the positions of points on the freeze-out hypersurface. Their ordering is chosen in such a way that dΣ µ points in the direction of decreasing temperature. The freeze-out hypersurface obtained from the most general (3+1)D hydrodynamic evolution (without assuming any symmetries) may be parametrized in the following way [11] (see Fig. 7.1) x = d(ζ, φ, θ) sin θ cos ζ cos φ, y = d(ζ, φ, θ) sin θ cos ζ sin φ, This parametrization leads to simple formulas for the space-time rapidity η, longitudinal proper time τ , and the transverse distance r 1 , The parameters a, b and c in Eq. (7.8) are identified with three angles ζ, φ and θ, respectively, which are restricted to the ranges: The quantity d(ζ, φ, θ) describes the distance between a point on the freeze-out hypersurface and the coordinate system's origin (τ = τ 0 , x = 0, y = 0, η = 0). The parameter τ 0 is the initial proper time for hydrodynamic evolution, at which the initial conditions are specified. The parameter Λ introduces a scale which changes dimensionless space-time rapidity η into a quantity with dimension of fermi. The parametrization (7.9) works for all typical freeze-out hypersurfaces where the distance d is a function of ζ, φ, and θ.
Substitution of (7.9) in (7.8) leads to the following form of the Cooper-Frye integration measure Using formulas (7.7) and (7.12) in (7.3) we obtain obtain a six-dimensional particle distribution which one can use to generate both stable hadrons and unstable resonances on the freeze-out hypersurface.

Extraction of physical observables
After performing hydrodynamic calculations with ADHYDRO, the freeze-out hypersurface is extracted. This procedure yields four functions: d(ζ, φ, θ), u x (ζ, φ, θ), u y (ζ, φ, θ) and ϑ(ζ, φ, θ) which are subsequently used by the THERMINATOR Monte-Carlo generator as an input. THERMINATOR generates physical events using the Cooper-Frye formula (7.3). In the single event all known hadron species (stable and unstable) are generated on the freeze-out hypersurface. Subsequently, the unstable primordial particles decay through strong and electro-weak interactions. The final particle distributions consist of stable primordial particles and contributions from resonance decays. The complete information about all generated particles and their decays is kept in the program. This information allows for the calculation of soft hadronic observables typically measured in heavy-ion experiments.
Probably, the most important information about the evolution of matter created in relativistic heavy-ion collisions is contained in the momentum distributions of identified particles. Experimental data on transverse-momentum distributions are often corrected for contributions from weak decays (the contributions from weak decays to the hadronic yields and spectra are experimentally determined and subtracted). The detailed output of THERMINATOR allows for applying such corrections in the model calculations in the easy way, making comparisons to the data more realistic.
In non-central collisions of two heavy nuclei, the experimentally observed momentum distribution of particles is azimuthally anisotropic and may be expressed as a Fourier series in the momentum azimuthal angle, φ p = arctan(p y /p x ), with respect to reaction plane [96] 2v n (p T , y) cos(nφ p ) . (7.13) In the hydrodynamic calculations where the initial conditions are defined by a smooth density profile taken from the optical Glauber formula (6.2) the sine terms disappear due to symmetry with respect to the reaction plane. In event-by-event calculations where the fluctuations of initial conditions are present this is not the case. The coefficient v n (p T , y) is the n-th harmonic differential flow which quantitatively characterizes the momentum azimuthal anisotropy. The v 1 (p T , y) coefficient is called the directed flow, and the v 2 (p T , y) coefficient is called the elliptic flow. Due to symmetry reasons the v 1 (p T , y) vanishes for y ≈ 0. The v 3 coefficient is called the triangular flow. The higher harmonics are typically negligible. For instance, the v 4 coefficient is predicted by hydrodynamic calculations to be at the level of ∼ 0.1%.
Having complete information about the particles generated by the THERMINATOR code one can perform the femtoscopic analysis. The calculation of the correlation function for the identical pions is performed using the two-particle Monte-Carlo method described in Ref. [97,98] and implemented in THERMINATOR 2 [85] code. In the numerical procedure the ideal Bose-Einstein correlation function is evaluated using the formula (7.14) The calculation of the correlation function (7.14) is performed with the help of the box function δ ∆ defined in the following way (7.15) where one uses the bin resolution ∆ = 5 MeV. In Eq. (7.14) we use the momentum difference of the pions, 16) and the average momentum of the pion pair, We also consider the space-time separation between pions ∆x = (t, r) = (t 1 − t 2 , r 1 − r 2 ) . (7.18) In the HBT analysis we use the standard Bertsch-Pratt system of coordinates [99,100], where the long direction is parallel to beam axis, the out axis coincides with the transverse momentum of the pion pair k T and the side axis is perpendicular to out and long axes. We perform the analysis in the longitudinally co-moving system (LCMS) in which k long = 0. For each pair of pions the wave function in Eq. (7.14), Ψ(q * , r * ) = 1 √ 2 (e iq * r * + e −iq * r * ), (7.19) is calculated in the pair rest frame (PRF). The q * and r * denote the relative momentum and position calculated in PRF. In the calculation of expression (7.14) we do not include Coulomb effects. Using the Bertsch-Pratt coordinates k T , q out , q side , q long the correlation function (7.14) may be approximated by the formula where the HBT correlation radii R out , R side and R long are the widths of the Gaussian source. From the fit we obtain parameters R out , R side and R long as a function of k T . In order to make a comparison of our results with the available data from STAR we consider four k T bins

Boost-invariant systems with transverse expansion
The data from RHIC and the LHC shows that the observables measured in the central rapidity region are approximately invariant with respect to the Lorentz boost transformations along the beam axis. Thus, in this Chapter, we use ADHYDRO assuming the boost-invariance along the beam direction. In this way, we restrict our analysis to the midrapidity region (|y| < 1) of ultra-relativistic heavy-ion collisions.

Central collisions
First we consider central Au+Au collisions at the energy of √ s NN = 200 GeV at RHIC. We choose the centrality class c = 0 − 5%, which corresponds to the impact parameter b = 2.26 fm. We check three scenarios of the initial momentum anisotropy: i) x 0 = 100 (dashed blue lines), which means that the system is initially oblate in the momentum space, the fact predicted by the microscopic models of heavyion collisions, ii) x 0 = 1 (solid black lines), which corresponds to initially isotropic systems in the momentum space, this case gives the closest description to perfect fluid hydrodynamics, and iii) x 0 = 0.032 (dotted green lines), which means that the system is initially prolate in momentum space, which is less-likely to be realized in nature than the case i), however, we consider this case for completeness.  The results of the ADHYDRO model are compared to the reference perfect-fluid calculations performed with the LHYQUID code [101,82]. In our calculations we use the initial conditions described in Chapter 6.1. The value of the initial energy density ε i = 86.76 GeV/fm 3 in Eq. (6.1) is the same for all analyzed cases. For the system in local thermal equilibrium, it corresponds to the initial central temperature T i = 485 MeV. Following Eq. (6.9), the initial central entropy density is σ i = σ(τ 0 , x ⊥ = 0) = 227 fm −3 for the case x 0 = 1, and σ i = 151 fm −3 for the cases x 0 = 100 and x 0 = 0.032. The value of the time-scale parameter τ eq = 0.25 fm used in our calculations guarantees that the system reaches local thermal equilibrium in the whole fireball volume approximately after 1 fm (this value for thermalization time is easily acceptable by most of the microscopic model calculations).
In the left panel of Fig. 8.1 we present the time dependence of the central entropy density σ(τ, x ⊥ = 0) for three considered cases of the initial anisotropy parameter: x 0 = 100 (dashed blue line), x 0 = 1 (solid black line), and x 0 = 0.032 (dotted green line). We compare our results to the perfect-fluid case (dashed-dotted red line). We observe that in the cases x 0 = 100 and x 0 = 0.032 the entropy density grows initially. This behavior is different from the entropy-conserving perfect-fluid case where σ(τ, x ⊥ = 0) decreases like 1/τ (dashed-dotted red line). For τ ≈ 0.3 fm σ(τ, x ⊥ = 0) reaches maximum and for τ > 0.5 fm it starts to behave similarly as in the perfect-fluid case. The case x 0 = 1 gives results similar to the perfect-fluid case. The entropy densities at the end of evolution are similar, however, small differences are noticeable. It means that the rescaling of the initial entropy density σ i according to Eq. (6.9) assures only approximately the same final multiplicity and entropy density in all considered cases. In the right panel of Fig. 8.1 we show the time dependence of the entropy source Σ(τ, x ⊥ = 0) defined by Eq. (4.22). The notation is the same as in the left panel. The function Σ(τ, x ⊥ = 0) contains information complementary to that encoded in σ(τ, x ⊥ = 0). Due to strong initial anisotropy, in the case x 0 = 100 the initial entropy production is very large. As the pressure anisotropy decreases (see the corresponding time dependence of the anisotropy parameter x(τ, x ⊥ = 0) in the left panel of Fig. 8.2) the production of entropy becomes smaller. In the case x 0 = 0.032, the entropy production is also quite large initially. However, it quickly drops to zero when x passes unity (see also Chapter 5.6). Subsequently, it increases again and becomes similar to the values found for the case x 0 = 100 1 . When x 0 = 1, the anisotropy remains always close to unity. Hence, the entropy production is very small in this case. We note that in all the cases the anisotropy parameter tends asymptotically to unity, which means that for large times the entropy source vanishes and the system approaches the perfect-fluid regime.
The initial large anisotropy in the case x 0 = 100 is equivalent to the large P ⊥ /P ratio which implies strong hydrodynamic flow in the transverse direction at the early stages of the evolution (small longitudinal pressure means that no work is done in the longitudinal direction and, consequently, thermal energy is used to build the transverse flow only). This effect is, however, compensated by the smaller value of the initial entropy density obtained from Eq. (6.9). Thus, our results obtained with ADHYDRO are very much similar to those obtained in the standard hydrodynamics. This property is seen in the right panel of Fig. 8.2, where we plot in-plane (φ = 0) radial velocity v ⊥ as a function of proper time τ for different values of the distance r, for the cases i), ii), and iii) 2 . In Fig. 8.3 we present in-and out-of-plane freeze-out hypersurfaces determined from the numerical solutions of the ADHYDRO equations. The freeze-out condition is specified only for the locally equilibrated matter, thus, we exclude the space-time regions where τ < 1 fm. In those regions the matter is anisotropic and consists mainly of gluons -it is not likely that it contributes significantly to the hadron production.
We use the freeze-out temperature value T f = 150 MeV which has turned out very much satisfactory to reproduce the experimental data from RHIC using perfectfluid hydrodynamics. After extracting the freeze-out hypersurface from the hydrodynamic code, we use THERMINATOR 2 [85] in order to generate physical events. In Fig. 8.3 we can see that for x 0 = 100 the system cools down a little bit slower than in the standard perfect-fluid hydrodynamics, so the freeze-out hypersurface is larger.  [102]. Protons are feed-down corrected for Λ 0 . Right panels: The calculated transverse-momentum spectra of hyperons Λ 0 (upper part), Ξ − (middle part), and Ω + + Ω − (lower part) for Au+Au collisions at √ s NN = 200 GeV and the centrality class c = 0 − 5% compared to the STAR data [103]. Λ 0 corrected for feed-down from Ξ and Ω.
At the same time we find that the values of the hydrodynamic flow are slightly larger. This, as will be shown later, leads to a slightly larger multiplicity and harder spectra. In the case x 0 = 0.032, we deal with the opposite situation. Despite these small differences, in all the cases the freeze-out hypersurfaces are similar. This suggests that the initial anisotropic stage of evolution does not modify significantly the freeze-out conditions, provided the initial entropy density is appropriately rescaled in each case.
In the left panels of Fig. 8.4 we show the transverse-momentum spectra of pions (upper part), kaons (middle part), and protons (lower part), while in the right panels of Fig. 8.4 we present the transverse-momentum spectra of Λ 0 (upper part), Ξ − (middle part), and Ω + + Ω − (lower part). One can observe that for all values of the initial anisotropy the spectra are very similar and lie close to the (2+1)D  We close the study of the central collisions with Fig. 8.5 which presents the results for k T -dependent HBT radii. We observe the overall agreement with the experimental data for the perfect-fluid calculation performed with LHYQUID. We also see that the correlation radii show almost no sensitivity to the initial anisotropic stage. All calculations yield similar results. This is understandable, since we obtain similar shapes of the freeze-out hypersurfaces for all considered initial values of x 0 .
The analysis of the most central collisions in the framework of ADHYDRO and LHYQUID suggests that the considered observables are almost insensitive to the initial anisotropic stage as long as the initial entropy density is properly rescaled.

Non-central collisions
In this Section we turn to the analysis of non-central collisions. We take into consideration two centrality bins: c = 20 − 30% (b = 7.16 fm) and c = 20 − 40% (b = 7.84 fm). The same value of the initial central energy density, ε i = 86.76 GeV/fm 3 , is used now, as before for central collisions. The initial entropy density at the center of In the left part of Fig. 8.6 we present the in-plane (φ = 0) and out-of plane (φ = π/2) freeze-out hypersurfaces for the centrality class c = 20 − 30%. We come to similar conclusions as in the case of central collisions. Moreover, we observe that the pressure acts in the same way in the in-plane and out-of-plane directions. This feature is expected since we introduced the anisotropy only between the longitudinal and transverse pressures. We also observe that the transverse spatial eccentricity of the fireball, ε x , decreases faster with time in the x 0 = 100 case than in the x 0 = 0.032 case, see Fig. 8.7. This implies that the transverse momentum anisotropy ε p gets built faster. However, these effects are small if the equilibration time is short.
In the right panel of Fig. 8.7 we show the differential elliptic flow coefficient v 2 (p T ) of pions as a function of the transverse momentum, for the centrality class c = 20 − 30%. The model results describe quite well the data for p T ≤ 1 GeV. For different initial anisotropy parameters we obtain almost the same results. The slightly larger values of v 2 are observed in the case x 0 = 100, which is consistent with the conclusions coming from the analysis of the parameter ε p . However, for small τ eq that is used by us, this effect is negligible. In Fig. 8.8 we show our model results for the transverse-momentum spectra of pions (upper left part), kaons (middle left part), protons (lower left part), Λ 0 (upper right part), Ξ − (middle right part), and Ω + + Ω − (lower right part) obtained for the centrality classes c = 20 − 30% and c = 20 − 40% and compared with the RHIC data from PHENIX [102] and STAR [103]. One can notice the very good agreement of the model results with the experimental data, except for Λ 0 whose normalization is too large. The model describes the experimental data up to p T = 2 GeV, thus it works well in the region of applicability of hydrodynamics. The discrepancies for large p T are more significant than in the central collisions. Similarly to central collisions we observe improvement of the high p T tails for the case of x 0 = 100. Figure 8.9 presents the k T dependence of the HBT radii for the centrality class c = 20 − 30%. Similarly to central collisions the agreement with the data is quite satisfactory. There is almost no effects of the early anisotropic stages.
Summarizing, the systematic study of the observables in the (2+1)D boostinvariant version of ADHYDRO model brings us to the conclusion that the observables are almost insensitive to the initial highly-anisotropic stage, provided the initial entropy density is properly rescaled. Such results are consistent with the universality of the flow predicted in Ref. [106].

Non-boost-invariant description of RHIC data
The RHIC results obtained outside the midrapidity region [107,108,74,109,110,111], in particular, the measured rapidity distributions of charged particles, suggest that any realistic description of the data should relax the assumption of Bjorken boost-invariance. Following this observation, in this Chapter we use the ADHYDRO model in 3+1 dimensions. We do not impose any symmetries to analyze the effects of initial anisotropic stages on the subsequent behaviour of matter. The analysis is done for Au+Au collisions studied at RHIC at the energy √ s NN = 200 GeV.
In order to make a comparison with the experiment, we use the ADHYDRO model coupled to the Monte-Carlo model of statistical hadronization THERMINATOR [84,85] as described in Chapter 7.

Constant initial momentum anisotropy
First we analyze the simplest case where the initial anisotropy x 0 is constant in space. Once again, we perform calculations for three different cases: i) x 0 = 100 (dashed blue lines), ii) x 0 = 1 (solid black lines), and iii) x 0 = 0.032 (dotted green lines). We keep this notation all the way in this Chapter. We also use three values of the time-scale parameter: τ eq = 0.25 fm, τ eq = 0.5 fm, and τ eq = 1.0 fm to check how the results differ with respect to the length of the anisotropic stage. We show the results for the fixed freeze-out entropy density σ f = 1.79/fm 3 which corresponds to the freeze-out temperature T f = 150 MeV. This value of T f correctly describes the slope of p T spectra of different hadron species, thus satisfactory separating the flow component from the thermal component in transverse-momentum spectra. The half-width of the central plateau in the initial profile, ∆η = 1, and the half-width of gaussian tails, σ η = 1.3, are fixed, unless specified otherwise. They are fitted to approximately reproduce the shape of charged particle distributions in pseudorapidity in the case x 0 = 1. The initial energy density in the center of the fireball, ε i , is chosen in such a way as to reproduce the total experimental charged particle multiplicity N exp ch = 5060 ± 250 [107] for the centrality class c = 0 − 6%. For the case  x 0 = 1 we use the energy density ε i = 107.5 GeV/fm 3 which yields N theor ch = 5020. Since entropy grows during the hydrodynamic evolution due to non-zero entropy production source Σ(τ eq , σ, x), see Eq. (4.22), the initial central energy density must be appropriately renormalized for each pair of parameters x and τ eq to keep N theor ch unchanged. In the case x 0 = 100 the total entropy produced during the hydrodynamic evolution grows with increasing τ eq , while in the case x 0 = 0.032 the total produced entropy decreases with increasing τ eq (similarly to the one-dimensional behavior shown in Figs. 5.1 and 5.3 in Chap. 5.6).
The values of ε i are summarized in Table 9.1. The centrality dependence is reproduced by applying the scaling (6.3). The initial conditions are given by the tilted source described in Chapter 6.2.

Particle spectra
We start our discussion of the results obtained with the ADHYDRO model with the analysis of the momentum distributions of charged and identified particles. These distributions include the feeding from unstable resonances whose decays are implemented in THERMINATOR [84,85].

Pseudorapidity distribution
At first, we consider the pseudorapidity distributions of charged particles. The pseudorapidity is defined by the formula where p is the momentum of a particle. In Fig. 9.1 we show dN ch /dη P S for three different values of the initial momentum anisotropy: x 0 = 100 (dashed blue lines), x 0 = 1 (solid black lines), and x 0 = 0.032 (dotted green lines), and three values of the time-scale parameter: τ eq = 0.25 fm, τ eq = 0.5 fm, and τ eq = 1.0 fm. We observe that as long as the initial anisotropic stage is short its influence on pseudorapidity distributions of particles is small in both x 0 = 100 and x 0 = 0.032 cases, as compared to the case x 0 = 1. As the length of the off-equilibrium stage increases, the difference of pressures lasts longer, and the difference between the longitudinal and transverse expansion rates becomes larger. For instance, in the case x 0 = 100 the pseudorapidity distribution becomes a bit steeper in the region of forward and backward rapidities due to reduced longitudinal pressure (see our discussion in Chapter 5.3), but the width of the central plateau is nearly unchanged. In Fig. 9.2 we present the case x 0 = 1 with τ eq = 0.25 fm and the standard parameters of f (η) together with the case x 0 = 100 with τ eq = 1.0 fm. In the latter case the initial energy density was reduced to ε i = 41.8 GeV/fm 3 (from 48.8 GeV/fm 3 , see Table 9.1) and the width of the plateau of the initial profile f (η) was refitted to ∆η = 1.5. In this way we obtain the same total multiplicity. We can clearly see the improvement in the description of data as compared to the case x 0 = 100 in Fig. 9.1 (c). Thus, the shape of pseudorapidity distributions might be better de- In the latter case the initial conditions were modified to ε i = 41.8 GeV/fm 3 (compare Table 9.1) and ∆η = 1.5 to obtain the same final total multiplicity. The results are compared to the experimental PHOBOS [107](red dots) and BRAHMS [108](green squares) data.
scribed if the initial fireball is transversally thermalized and the profile parameters are properly readjusted. However, it should be stressed that the widths and forms of the initial longitudinal profiles (6.12) are not well known. Therefore, no definite conclusions about the length of the pre-equilibrium phase may be drawn from the correct description of the shape of the rapidity distributions alone, as it has been already noticed in [112].

Transverse-momentum spectra
In Fig. 9.3 we present the transverse-momentum spectra of positive pions (left part) and kaons (right part) for the centrality class c = 0 − 5% (b = 2.23 fm) and for different rapidity windows. The model calculations are done with τ eq = 1.0 fm. The model results are compared to the BRAHMS data [74]. We reproduce well the data in two different cases x 0 = 100 and x 0 = 0.032. The correct description of the p T spectra in a broad range of rapidities is usually understood as an indication for a high level of thermalization of the entire fireball [11]. After appropriate rescaling of the initial central energy density, the spectra are almost insensitive to initial values of anisotropy and to the length of the anisotropic stage (note a relatively large value of τ eq in Fig. 9.3). The small discrepancies in normalizations between different rapidity windows are mainly due to discrepancies between theoretical and experimental dN ch /dη P S profiles (compare with the part (c) of Fig. 9.1). In the left part of Fig. 9.4 we present the transverse-momentum spectra of positive pions (top), kaons (middle), and protons (bottom) at midrapidity (|y| < 1) for the time-scale parameter τ eq = 1.0 fm and for centralities c = 0 − 5% (b = 2.26 fm) and c = 20 − 30% (b = 7.16 fm). The model results are compared to the data from PHENIX [102]. We can see that the model reproduces very well the slope of the spectra of K + 's and protons up to p T = 3 GeV. The model pion spectra for c = 20 − 30% are underestimated for p T > 2.5 GeV. This is an expected feature, since hard particles require more rescatterings to thermalize and this region is out of applicability of hydrodynamics.
The correct normalization of spectra for different hadron species is assured by introducing thermodynamic parameters at freeze-out, which were obtained in the framework of a thermal model using a complete treatment of resonances [94]. The reproduction of the convex curvature of the pion spectra at low p T is a consequence of feeding from resonance decays which are implemented in THEMINATOR [84,85].
On the right-hand-side of Fig. 9.4 we present the transverse-momentum spectra of hyperons for centralities c = 0 − 5% (b = 2.26 fm) and c = 20 − 40% (b = 7.84 fm) compared to the data from STAR [103]. Although our freeze-out temperature is a bit lower than that used in thermal models, the slopes of the hyperon spectra are reproduced very well. The normalization of Λ's is too small, the effect that may indicate their higher freeze-out temperature. Fig. 9.4 shows again a small sensitivity of the model predictions to initial anisotropy, which is consistent with conclusions based on Fig. 9.3.
One would naively expect that the increase (decrease) of the pressure ratio P ⊥ /P in the case x 0 = 100 (x 0 = 0.032), compared to the perfect-fluid case, should lead to slower (faster) cooling of the system in the transverse direction. However the renormalization of the initial energy density (according to Table 9.1) in order to obtain the same final multiplicity leads to a reduction (amplification) of the transverse flow produced during expansion in the anisotropic phase [113]. This explicitly confirms the phenomenon of universality of the flow, predicting that the overall growth of the flow is the same regardless of the pressure anisotropy of the energy-momentum tensor [106].

Anisotropic flows
The significant amount of anisotropic flows is built at the very early stages of the evolution of matter, when the gradients of pressures are the largest due to the large spacial anisotropy of the source. The inclusion of early non-equilibrium stages of the collisions may significantly change the momentum anisotropy observed in the final spectra, making v n coefficients sensitive probes for the possible existence of non-equilibrium stages. at midrapidity for the time-scale parameter τ eq = 1.0 fm, and for two values of the initial anisotropy parameter. The results are compared to the PHENIX Collaboration data (red dots) [114]. The presented errors are statistical. Horizontal bars denote the p T bins.

Elliptic flow
In the midrapidity region of non-central Au+Au collisions, the main information about the momentum anisotropy of the produced particles is contained in the v 2 coefficient. In the left part of Fig. 9.5 we show the elliptic flow of π + + K + as a function of the transverse momentum p T . The model results obtained for the centrality c = 20−40 % (b = 7.84 fm) and with the time-scale parameter τ eq = 1.0 are compared to the data from PHENIX [114]. We observe that the model reproduces well the data up to p T = 1 GeV. The saturation for higher p T is not reproduced, since it requires inclusion of shear viscosity effects [115]. Again, the model results weakly depend on the initial value of anisotropy, since the increase (decrease) of transverse pressure in the initial anisotropic stage is compensated by the decrease (increase) of the initial energy density. For completeness, in the right part of Fig. 9.5 we show the p T -dependence of the elliptic flow of protons. The mass-splitting of v 2 is not reproduced. However, the elliptic flow of protons may be better described by introducing bulk viscosity [116].
In Fig. 9.6 we present the p T -integrated elliptic flow as a function of pseudorapidity. The model calculations are done for the centrality class c = 15 − 25%, three time-scale parameters, and three values of the initial anisotropy parameter. The model results are compared to the PHOBOS data [109]. We observe that the model overshoots the data by about 20% in the case of x 0 = 1 (solid black line). The peaked shape visible in the data is not well reproduced.
We observe that the elliptic flow increases when the transverse pressure dominates. In this case, the elliptic flow is also steeper in pseudorapidity. It is easy to understand this behavior, since the reduction of longitudinal pressure leads to slower expansion of the fireball in the longitudinal direction. In the case of dominating longitudinal pressure the behavior is opposite. The strength of this effect depends on  the length of the anisotropic stage. Thus, we conclude that the p T -integrated v 2 is sensitive to the early anisotropic stages.

Directed flow
When one considers non-central collisions of heavy-nuclei, one may observe nonzero directed flow in the region where η P S = 0. This observable is sensitive to both transverse and longitudinal pressures, which makes it a perfect probe for measuring the early local anisotropies in momentum [113]. In Fig. 9.7 we show the directed flow of charged particles for the two centrality bins c = 0−5% (left part) and c = 5−40% (b = 6.79 fm, right part), and for the two time-scale parameters τ eq = 0.25 fm (a) and τ eq = 1.0 fm (b). The results are plotted together with the experimental points from STAR [110] and PHOBOS [111].  and three values of initial anisotropy parameter. The results are compared to experimental data from STAR (red dots) [110] and PHOBOS (green squares) [111]. Errors are statistical only.
For different choices of our model parameters, the results describing v 1 are rather stable and consistent with the data. This leads to the conclusion that all the anisotropic stages described by ADHYDRO in this work are consistent with the v 1 data. Moreover, the agreement with the data suggests the validity of the idea of a tilted initial source.

HBT correlations
The identical particle interferometry is a useful tool used in the analysis of the system sizes at freeze-out. In Fig. 9.8 we present the model results together with the STAR data [104] for the HBT radii: R out (a), R side (b), and R long (c), all shown as functions of the pair total transverse momentum k T . Differences between the   Figure 9.8: The HBT correlation radii: R out (a), R side (b), R long (c) and ratio R out /R side (d) of positive pions as a function of total transverse momentum of the pair for centrality c = 0 − 5% for time-scale parameter τ eq = 1.0 fm and three values of initial anisotropy parameter. The results are compared to experimental data from STAR Collaboration (red dots) [104]. Errors bars reflect statistical errors. data and the model results are smaller than 10%. We observe that the slopes of the radii are quite well reproduced by our model. At large k T , the side radius, R side , is slightly smaller than that observed in the experiment, and the out radius, R out , is slightly larger. This is reflected in the overprediction of the R out /R side ratio by about 20%. The k T slope of R out /R side is not reproduced.
For the all considered values of x 0 , we obtain almost identical results. This is expected since the correct renormalization of the initial energy density results in the similar observed multiplicity and the flow on the freeze-out hypersurface (see our discussion in Sections 9.1.2 and 9.2). The similarities in the description of the correlation radii for all x 0 's mean that the shapes of the freeze-out hypersurfaces are also very much similar. The last point was already explicitly highlighted in Fig. 8.3 in Chapter 8.

Space dependent initial anisotropy
In Ref. [117] it has been suggested that the system may be less equilibrated in forward and backward rapidities (i.e., in the fragmentation regions) than in the midrapidity region. On the other hand, in Ref. [118] it has been argued that the thermalization may depend on the density of participants, since the degree of thermalization depends on the number of collisions between nucleons. Thus, one may expect that, in general, x 0 = x 0 (ρ(η, x ⊥ )), where x 0 ≫ 1 at the edges of the system.
Here we skip the possible dependence on centrality and check two scenarios: i) a constant initial anisotropy profile x 0 = 100 with τ eq = 1.0 fm, ε i = 41.8 GeV/fm 3 , ∆η = 1.5, and σ η = 1.3, see Fig. 9.2 in Chapter 9.2, and ii) a space dependent initial anisotropy profile x 0 = x 0 (η, x ⊥ ) where with τ eq = 1.0 fm, ε i = 73.8 GeV/fm 3 , ∆η = 1, and σ η = 1.3. Fig. 9.9 shows the pseudorapidity dependence of charged particle multiplicity for the case i) (blue dashed line) and ii) (black dotted line). We observe that the distributions are similar despite different initial anisotropy profiles (note that in the case ii) ∆η is much smaller, since in the central rapidity region the matter evolves faster in the longitudinal direction than in the case i)). Moreover, moving to forward rapidities the anisotropy grows due to the ansatz (9.2). Thus, the tails of the pseudorapidity distribution are described more precisely. The results concerning the transverse momentum spectra are very similar in both i) and ii) cases. On the basis of conclusions from Section 9.1.1 this is an expected result. Interesting observations may be done for the v 2 coefficient at midrapidity. In the case i) the imposed initial anisotropy is large and constant in entire fireball. The rescaling of the initial energy density cancels the effect of larger initial pressure in the transverse direction. Moreover the spatial anisotropy of the initially produced entropy in the case i), Σ(x 0 , σ(τ 0 , x ⊥ )), is the same as the spatial anisotropy of the initial entropy density profile. Hence the spatial anisotropy of the fireball transforms into the same momentum anisotropy as in the perfect-fluid case and we obtain similar results for v 2 (see also Section 9.1.2). In the case ii) the initial anisotropy is space dependent. In the center of the system it is almost isotropic due to high density of wounded nucleons. The anisotropy is much larger in the regions of low density of wounded nucleons, thus the spatial anisotropy of the initially produced entropy in the case ii) is lower than that in the case i). It results in a slightly lower spatial anisotropy of the whole fireball (see left part of Fig. 9.10) and in consequence it suppresses the generation of v 2 (see right part of Fig. 9.10). The effect of lowering of  the v 2 is clearly seen in the Fig. 9.11 where we present the pseudorapidity dependence of the p T -integrated v 2 . We observe better agreement of the calculations with the experimental data in the central rapidities in the case ii). However the model still does not reproduce the peaked shape seen in the data. The spatial dependence of the initial anisotropy has a small effect on the HBT radii. In Fig. 9.12 we present the HBT radii calculated for the case i) (blue dashed line) and ii) (black dotted line). We observe a small improvement of R side which results in better agreement of R out /R side in the case ii) as compared to the case i).  In this Thesis, a recently developed framework of highly-anisotropic and stronglydissipative hydrodynamics -ADHYDRO -has been introduced and used to analyze the space-time evolution of matter produced in ultra-relativistic heavy-ion collisions. The main goal of this analysis was to study the effect of initial highly-anisotropic stages on the final soft hadronic observables typically measured in the experiment. The study was done in the context of the heavy-ion measurements performed at RHIC (Relativistic Heavy Ion Collider) w Brookhaven National Laboratory.
Starting from the general assumption about the form of the phase-space distribution function, thermodynamic properties of locally anisotropic systems of particles have been studied and the form of the generalized equation of state has been formulated. The dynamic equations determining the evolution of a highly-anisotropic fluid have been introduced. The form of the entropy source related to the mechanisms leading to thermalization of the system has been defined. In the simplest case of purely-longitudinal and boost-invariant expansion, different features of the model have been analyzed. In particular, for small deviations from equilibrium, the consistency of the ADHYDRO model with the Israel-Stewart theory has been demonstrated.
Using the ADHYDRO model in the general (3+1)D and boost-invariant (2+1)D versions, the three possible scenarios of early stages of heavy-ion collisions have been analyzed: i) the system initially oblate in momentum space, ii) the system initially prolate in momentum space, iii) the initially isotropic system. The results of the ADHYDRO model have been compared to the results obtained from the reference, perfect-fluid hydrodynamic model -LHYQUID. The results of the hydrodynamic models were coupled to the statistical Monte-Carlo model THERMINATOR which has allowed us to perform systematic study of soft-hadronic observables: transversemomentum spectra, directed and elliptic flows, and the HBT radii.
The following final conclusions may be drawn: a) All studied observables are almost insensitive to the initial anisotropic stage provided the initial conditions of the evolution are properly readjusted.

CHAPTER 10. SUMMARY AND CONCLUSIONS
b) Complete thermalization of matter may take place only at the times of about 1 fm/c. In this way the early thermalization puzzle may be circumvented.
Finally, we have studied the systems where the initial anisotropy is spatially dependent. Our results suggest that the system may be highly anisotropic at the edges and simultaneously well thermalized in the center. It should be stressed, however, that the correct description of the data has been achieved only in the situations where the complete thermalization has been reached at about 1 fm/c.

Initial conditions
ρ density of sources, ρ normalized density of sources, ε i initial energy density at the center of the system in most central collisions, ρ W wounded-nucleon density, ρ B binary-collisions density, ρ WS Woods-Saxon nuclear density profile, κ mixing factor between wounded-nucleon, and binary-collisions densities, T thickness function, σ in total inelastic N N cross-section, ρ 0 , R A , a parameters for Woods-Saxon nuclear profile, b impact vector, b = |b| impact parameter, f (η) initial longitudinal profile, ∆η, σ η parameters of the initial longitudinal profile, θ step function, m N nucleon mass, y b rapidity of the participant nucleon, η m range of rapidity correlations,