{'fraction_non_alphanumeric': 0.028363713550372788, 'fraction_numerical': 0.015654434981028947, 'mean_word_length': 4.638839018551765, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hydrodynamics at RHIC and LHC: What have we learned?

Fluid dynamical description of elementary particle collisions has a long history dating back to the works of Landau and Fermi. Nevertheless, it is during the last 10-15 years when fluid dynamics has become the standard tool to describe the evolution of matter created in ultrarelativistic heavy ion collisions. In this review I briefly describe the hydrodynamical models, what we have learned when analyzing the RHIC and LHC data using these models, and what the latest developments and challenges are.


Introduction
The goal of the heavy-ion programs at BNL RHIC a and CERN LHC b is to observe strongly interacting matter. 'Strongly interacting' in a sense that the interactions in the system are not mediated by the electromagnetic, but by the strong interaction, and 'matter' in the sense that to describe the system, we do not need to describe every (quasi-)particle individually, but thermodynamic concepts like temperature and pressure are applicable. Thus it is natural to try to describe the expansion stage of the collision by a macroscopic approach like fluid dynamics. We expect that the system formed in the collision is initially so hot and dense that relevant degrees of freedom are partons, not hadrons. The system expands and cools, undergoes a phase transition to hadrons, and when the system is dilute enough, interactions cease and the particles stream freely to detectors. Relativistic fluid dynamics is basically an application of conservation laws for energy, momentum and conserved charges (if any). When written in differential form, where T µν is the energy-momentum tensor, and j µ x the charge 4-current of charge x (baryon or electric charge, strangeness, isospin,. . . ), the conservation laws provide evolution equations for the system. If there are n conserved charges, there are 4 + n equations, which contain 10 + 4n unknowns. To make this set of equations solvable, further constraints for the unknowns must be provided. The simplest approach is to assume that the system is in exact local equilibrium. In that case the energymomentum tensor and charge currents can be expressed as T µν = (ǫ + P )u µ u ν − P g µν , and respectively, where ǫ and P are energy density and pressure in the local rest frame, u µ flow velocity, g µν is the metric tensor, and n x the charge density of charge x in the local rest frame. Now the number of unknowns is reduced to 5 + n, and the system of equations can be closed by providing the equation of state (EoS) of the fluid in a form P = P (ǫ, {n i }). To model the collisions at RHIC and LHC, the description is usually further simplified by assuming that net baryon density and other conserved charges are zero. After specifying the boundary conditions for this set of partial differential equations the evolution is determined, and all the microscopic physics is contained within the EoS c . However, the simple phrase "after specifying the boundary conditions" contains a lot of physics. Hydrodynamics provides neither the initial distribution of matter nor the criterion for the end of the evolution, but they must be supplied by other models. At ultrarelativistic energies the initial state of the hydrodynamical evolution cannot be two colliding nuclei: The initial collision processes are far from equilibrium and produce large amount of entropy. 3 Thus the initial state is assumed to be a distribution of thermalized matter soon after the initial collision. The most common approach is so-called Glauber model, 4,5 which is basically a geometrical constraint on the distributions. The Woods-Saxon distributions of nuclear matter in colliding nuclei are projected on a plane orthogonal to the beam (so called transverse plane), and the resulting densities on this plane, and nucleon-nucleon cross section at the collision energy, are used to calculate the number density of binary collisions and participants on this plane-participant meaning a nucleon which has interacted at least once (for details, see Ref. 4). The initial energy or entropy density profile is taken to be proportional to the profile of collisions or participants, or to a linear combination of them. 6,7 The proportionality constant is a free parameter chosen to reproduce the observed final particle multiplicity.
Another popular approach is the KLN (Kharzeev-Levin-Nardi) model [8][9][10][11][12][13] where the initial entropy density distribution is proportional to the distribution of gluons produced in primary collisions. The gluon production is calculated using the Color Glass Condensate (CGC) framework, [14][15][16][17] where one applies the feature of QCD that at small-x gluon densities are large. These large densities correspond to classical fields permitting calculations using classical techniques. Another approach to calculate the initial particle production from first principles is so-called EKRT saturation model 18-20 based on perturbative QCD + saturation framework. Besides these approaches, one can use event generators like UrQMD, 21-23 AMPT, 24,25 or EPOS 26 to generate the initial state for hydrodynamic expansion. Nevertheless, no model so far describes a dynamical process leading to thermalized matter, but thermalization has to be postulated and imposed by hand.
When the system expands and cools, the mean free paths increase. Ultimately mean free paths become so long that rescatterings cease, and particle distributions no longer evolve. Particle distributions are 'frozen' at that stage, and we say that the particles freeze out, or alternatively, that the particles decouple from each other. In practice this freeze-out should be a gradual process, but since implementing a gradual freeze-out is difficult d , it is approximated to take place suddenly on a surface of zero thickness. In such a case one can use so-called Cooper-Frye prescription to evaluate particle distributions on this surface: 29 where σ is the surface where the distribution is to be evaluated, dσ µ its normal 4-vector, f distribution of particles, and u 4-flow velocity of the fluid. Hydrodynamics does not tell when decoupling should take place, but the criterion is a free parameter, and its value is chosen to reproduce the observed p T -spectra. Usually a constant temperature or energy density is used as the criterion, although a more realistic criterion would be the ratio of the scattering rate of particles to the expansion rate of the system, i.e., the inverse Knudsen number 30, 31

Dissipative fluid dynamics
The ideal fluid assumption is extremely strong, and in nature gradients in the system always indicate deviations from equilibrium and thus dissipation. In nonrelativistic fluid dynamics Navier-Stokes equations are known to describe viscous fluid well, but unfortunately the relativistic generalization of Navier-Stokes equations is unstable and allows acausal solutions with superluminal signal propagation speeds. [32][33][34] This undesired behavior can be avoided if one assumes that the dissipative currents (shear stress π µν , heat flow q µ and bulk pressure Π) are not directly proportional to gradients in the system, but are dynamical variables which relax to their Navier-Stokes values on time scales given by corresponding relaxation times τ π , τ q and τ Π . The evolution equations for the dissipative currents can be derived phenomenologically from entropy current, 35,36 from kinetic theory using the Grad 14-moment ansatz, [37][38][39] or via gradient expansion. 40,41 The derivation from kinetic theory leads to the commonly used Israel-Stewart equations, but so far it has had a problem: Unlike Chapman-Enskog expansion, 42 which leads to relativistic Navier-Stokes equations, it is not a controlled expansion in some small parameter, in which one could do power counting and improve the approximation if necessary. This problem has recently been solved by rederiving viscous hydrodynamics using the method of moments, and ordering the terms in the expansion according to their (generalized) Reynolds and Knudsen numbers. 43 In lowest order in Knudsen and Reynolds numbers this approach leads to equations identical to the original Israel-Stewart approach, but it has been shown that an adequate description of heat flow would require inclusion of some higher order terms. 44 Strictly speaking all these derivations are for a single component system, and the derivation for a multi-component system has not been completed yet. [45][46][47] For the behavior of the fluid this does not matter-the different components are assumed to behave as a single fluid-but when the fluid is converted to various types of hadrons, this causes an additional uncertainty. The Cooper-Frye prescription (Eq. (3)) is equally valid in dissipative as in ideal system, but dissipation causes the distribution function f to deviate slightly from the thermal equilibrium distribution. If there is only shear, no bulk pressure nor heat flow, in a single component system the Grad 14-moment ansatz leads in Boltzmann approximation to and f 0 is the equilibrium distribution function. In a multicomponent system the effect of shear stress on particle distributions should depend on particle properties, but how exactly, is still a work in progress. 48 In practical calculations one has therefore assumed that the dissipative correction to the thermal distribution is given by Eq. (4) for all hadron species. Note that the form of δf in Eq. (4) is only an ansatz. Other forms have been argued for, 49 but if the thermal distribution is expanded differently than in the 14-moment ansatz, the evolution equations and definitions of the dissipative coefficients may change as well. So far the viscous hydrodynamical calculations have concentrated on studying the effects of shear viscosity, characterized by the shear viscosity coefficient η. In no calculation has heat conductivity been included, and the bulk pressure has usually been omitted as well e . When the only dissipative current is shear stress, the energy-Hydrodynamics at RHIC and LHC: What have we learned? 5 momentum tensor becomes where π µν is the shear stress tensor. Its evolution is usually given by where D = u µ ∂ µ , σ µν = ∇ µ u ν , and the angular brackets denote the symmetrized and traceless projection, orthogonal to the fluid four-velocity u µ . Both Israel-Stewart 39, 43 and gradient expansion 40, 41 approaches lead to some additional terms in the evolution equation, but these terms are usually considered numerically insignificant, although a full analysis of the effect of all terms has not yet been done. It is worth noting that the factor 4/3 in the last term of Eq. (6) is strictly valid for massless particles. Whether modifying this term according to the mass of particles would change any results has not been studied yet. As well, in the original Israel-Stewart papers 36, 38 this term was estimated small and omitted in the final equations. In the context of heavy-ion collisions, however, it improves the applicability of viscous hydrodynamics significantly. 55

Hybrid models
The sudden change from interacting fluid to free streaming particles in the Cooper-Frye description is clearly an oversimplification. We also expect hadron gas to be so dissipative 45, 56 that the applicability of dissipative fluid dynamics is questionable. 57 In so-called hybrid models these problems are avoided by using fluid dynamics to describe only the early dense stage of the evolution. The fluid is converted to particles using the Cooper-Frye description, not at freeze-out, but when rescatterings are still abundant, and the individual particles are fed to a hadron cascade model like UrQMD 22,23 or JAM 58, 59 to describe the dilute hadronic stage. These models have the advantage that hadron cascades describe freeze-out dynamically without free parameters, all dissipative processes are included, and they can describe a system arbitrarily far from equilibrium. Unfortunately they do not remove all the arbitrariness from the end of the evolution: Like the results of pure hydrodynamical models depend on the freeze-out criterion, the results of hybrid models depend on when the switch from fluid to cascade is done. 57, 60 For a detailed discussion of these models see Ref. 60.

What we know
3.1. There are rescatterings The particle production in the primary collisions is azimuthally isotropic, but the distribution of observed particles in A+A collisions is not. The anisotropy can be easily explained in terms of rescatterings of the produced particles: In a non-central collision of identical nuclei the collision zone has an elongated shape, see  particle is moving along the long axis of the collision zone, it has a larger probability to scatter and change its direction than a particle moving along the short axis. Thus more particles end up in the direction of the short axis. Or, in the hydrodynamical language, the pressure gradient between the center of the system and the vacuum is larger in the direction of the short axis, the flow velocity is thus larger in that direction, and more particles are emitted along the short than along the long axis.
This anisotropy is quantified in terms of Fourier expansion of the azimuthal distribution. The coefficients of this expansion v n , and the associated participant angles ψ n , are defined as f v n = cos[n(φ − ψ n )] , and ψ n = 1 n arctan p T sin(nφ) p T cos(nφ) .
Of these coefficients v 1 is called directed, and v 2 elliptic flow. Elliptic flow of charged hadrons as a function of centrality was one of the first measurements at RHIC. 62 The result is shown in Fig. 2, and compared to early fluid dynamical calculations. 63 As seen, the elliptic flow is quite large and increases with decreasing centrality, as expected if it has the described geometric origin. Thus there must be rescatterings among the particles formed in the collision, and an A+A collision is not just a sum of independent pp collisions. The observed elliptic flow is also very close to the hydrodynamically calculated one, which is a strong indication of hydrodynamical behavior of the matter. Another signature of hydrodynamical behavior is shown in Fig. 3: It was observed that the heavier the particle, the smaller its p T -differential v 2 at low p T . As explained in Refs. 2 and 68, such a behavior arises if all the particles are emitted from the same expanding thermal source. Thus, if the produced matter is not close to kinetic equilibrium, at least it behaves as if it was! It has often been argued that the large observed v 2 and its hydrodynamical reproduction requires that the system reaches thermal equilibrium very fast, 69 within 1 fm/c after the initial collision. The hydrodynamical models fitting the data do indeed use short initial times, τ 0 = 0.15-1 fm/c, but in Ref. 70 it was shown that τ 0 = 2 fm/c works as well. The crucial distinction is what the shape of the initial state is: The larger the deformation, the larger the final momentum anisotropy. During thermalization the matter will not stay put, but begins to expand. This expansion reduces the spatial anisotropy, and unless momentum anisotropy is built up during thermalization, the final anisotropy is smaller. 71 If thermalization is fast, we may assume that changes in matter distribution and flow field during thermalization are tiny, and geometry (Glauber) and/or initial gluon production (Color Glass, EKRT) give reasonable constraints to the initial state of hydrodynamical evolution. In Ref. 70 the shape of the initial state was assumed to be independent of the thermalization time, but if thermalization time is long, this is no longer a good approximation. Thus, if thermalization is fast, we know that we can reproduce the data using hydrodynamical models, but if thermalization takes long, we do not know, since we do not know a plausible initial state for hydrodynamical evolution.

Equation of State has many degrees of freedom
The equation of state (EoS) of strongly interacting matter is an explicit input to hydrodynamical models. Thus one might expect hydrodynamical modeling of heavy-ion collisions to tell us a lot about the EoS, but unfortunately that is not the case. The collective motion of the system is directly affected by the pressure gradients in the system, and thus by the EoS, but the effects of the EoS on the final particle p T distributions can to very large extent be compensated by changes in the initial state of the evolution, and the final decoupling temperature. This makes constraining the properties of the EoS very difficult. However, what we do know is that the number of degrees of freedom has to be large.
It was already seen when modeling S+Au collisions at the CERN SPS at E lab = 200A GeV energy, that if we use ideal pion gas EoS, it is not possible to fit the pion rapidity and p T -distributions simultaneously. Once the initial state and decoupling temperature are fixed to reproduce the rapidity distribution, the transverse momentum distribution of pions becomes too flat, 72 see Fig. 4. Or if one chooses parameters to fit the p T -spectrum, the rapidity distribution is not reproduced. On the other hand, if we use an EoS containing several hadrons and resonances, the distributions can be fitted.
One might want to use the elliptic flow to constrain the EoS after the initial state and freeze-out temperature are fixed to reproduce the p T distributions. Unfortunately elliptic flow is only very weakly sensitive to the details of the EoS: 75 The only flow observable affected by the EoS seems to be the p T -differential anisotropy of heavy particles, e.g., protons. As shown in Fig. 5, the v 2 (p t ) of pions is unchanged within the experimental errors no matter whether one uses an EoS with (EoS A) or without phase transition (EoS H), or an EoS with a first order phase transition (EoS A) or with a smooth crossover (EoS qp). On the other hand, the proton v 2 (p t ) is sensitive to the EoS, but surprisingly the EoS with the first order phase transition is closest to the data. Consequently distinguishing between different lattice QCD EoS parametrizations is very difficult, see Ref. 76.

Softest point of the Equation of State is quite hard
The two-particle correlations at low relative momentum provide information about the space-time structure of the source, and thus constraints to dynamical models generating the source. For bosons these correlations are called Bose-Einstein correlations, and for identical particles the method for their interpretation is called Hanbury Brown-Twiss (HBT) interferometry, or, as in general case, femtoscopy. [77][78][79] To make the measured/calculated three-dimensional correlator easier to understand, it is often expressed in terms of multi-dimensional Gaussian parametrization, widths of which are called HBT radii. Reproduction of the HBT radii measured at RHIC was surprisingly difficult, and for many years it was not possible to describe simultaneously the particle spectra, their anisotropies, and HBT radii. An inconsistency referred to as "HBT puzzle". 80,81 It turned out that this inconsistency was due to several small effects. Successful reproduction of the data required that, 1) the transverse collective expansion begins very early, 2) the EoS is quite hard, 3) dissipative effects are included, and 4) contribution from resonance decays is included. [80][81][82][83] When these requirements are taken into account, the present calculations provide an acceptable fit to the data. 26,52,84 The early build up of collective expansion supports the notion of early thermalization discussed in Sec. 3.1, but the HBT radii provide as ambiguous support as elliptic flow. If sufficient collective motion is build up during thermalization, HBT radii can be fitted also if thermalization takes long (so-called pre-equilibrium flow). 80,81,83 The sufficient hardness of the EoS is fortunately more solid requirement. In practice it means that a bag model EoS with a mixed phase, where the speed of sound is negligible, is disfavored, and thus the softest point of the EoS can hardly be any softer than that of hadron resonance gas. 85 A notion in line with the lattice QCD EoS calculations. 86 Nevertheless, since it is so difficult to constrain the EoS, in the present calculations the lattice QCD EoS 86-89 is taken as given, and the main interest is in studying the dissipative properties of the system. At low temperatures the lattice EoS is equivalent to the EoS of noninteracting hadron resonance gas (HRG) with all the hadrons and resonances in the Particle Data Book up to ∼ 2-2.5 GeV mass. In calculations one often uses an EoS which is a HRG EoS at low temperature, and connected smoothly to a parametrized lattice EoS at high temperature. 76,90 In the long run a systematic study of collisions at different energies may reveal some sensitivity to EoS and thus allow to test experimentally the lattice QCD prediction. However, since the parameter space and the amount of data to be fitted are vast, such a study requires the use of model emulators to map the parameter space instead of using actual hydrodynamical model to calculate results at all parameter combinations. 91, 92

Shear viscosity over entropy density ratio has very low minimum
Once it became clear that ideal fluid dynamics can describe the particle spectra and their anisotropies quite well, it was reasonable to assume that the shear viscosity coefficient over entropy density ratio η/s of the matter produced in the collision was very low. But how low in particular? To answer that required the development and use of relativistic dissipative hydrodynamics. Of the other dissipative quantities heat conductivity can be ignored since at midrapidity the matter formed in the collisions at RHIC and LHC is almost baryon free. Thus there is no gradient of chemical potential and no force causing heat flow. The bulk viscosity, on the other hand, is expected to peak around the phase transition, but be small above it. If it is small at lower temperatures too, the effect of bulk viscosity has been evaluated to be smaller than the effect of shear viscosity. 51 Large viscosity in hadronic phase can have a sizable effect 50, 53, 54 but since bulk viscosity in hadron gas is not well known, and there is no reliable method to distinguish the effects of bulk from the effects of shear, the former is largely ignored, and the calculations concentrate on studying the effects of shear viscosity and on extracting the η/s ratio from the experimental data. 93 It has been shown that the shear viscosity strongly reduces v 2 . 94 Thus in principle extracting the η/s ratio from the data is easy: One needs to calculate the p T -averaged v 2 of charged hadrons using various values of η/s and choose the value of η/s which reproduces the data. Unfortunately this approach is hampered by our ignorance of the initial state of the evolution. Fig. 6 shows a viscous fluid calculation of v 2 of charged hadrons from the first attempt to extract η/s from the data. As seen, a curve corresponding to a finite value of η/s fits the data best, but the preferred value depends on how the initial state of hydrodynamic evolution is chosen: Whether one uses Glauber 4 or KLN approach 10 a factor two difference in the preferred value (η/s = 0.08-0.16). Furthermore, the approximations in the description of the late hadron gas stage in these calculations caused additional uncertainties, so it was estimated 96 that based on these results η/s < 5/(4π). The description of hadronic stage in Ref. 70 was problematic for two reasons. It assumed chemical equilibrium until the very end of the evolution, but at RHIC, the final particle yields correspond to a chemically equilibrated source in T = 160-165 MeV temperature. 97,98 The decoupling temperature required to fit the p Tdistributions is usually around T dec = 100-140 MeV g . Thus there must be a stage where the fluid is in local kinetic, but not in chemical equilibrium. Such a "chemically frozen" stage turned out to have surprisingly large effect on the collective flow in general and elliptic flow in particular. 99,100 The second problem was that η/s was constant during the entire evolution, but the theoretical expectation is that hadron gas has much larger viscosity than what the obtained value η/s = 0.08-0.16 indicates. 45,56 In subsequent calculations these uncertainties have been addressed. A state-ofthe-art calculation of Ref. 101 employs a hybrid viscous hydro + UrQMD hadron cascade model. Such a model reduces the uncertainty related to the description of hadron gas since it provides a dynamical description of chemistry and dissipative properties-including bulk viscosity-based on the scattering cross sections of various hadron species. Unfortunately some uncertainties still remain: It is uncertain whether hadron cascade is applicable in the vicinity of phase transition where the switch from hydro to cascade is done, the dissipative properties of the system (i.e. η, ζ, τ π and τ Π ) change abruptly at the switch, and the results also depend on  when the switch from hydro to cascade is done. 57 The elliptic flow coefficient v 2 as function of centrality as calculated in Ref. 101 is shown in Fig. 7. In this figure the coefficients v 2 have been scaled by the anisotropy of the initial shape, ǫ, and consequently the resulting v 2 /ǫ is almost independent of the initialization. Unfortunately ǫ is not a measurable, but a model dependent quantity. Thus the data has to be scaled by the same ǫ which was used in the calculation, and the data points in Fig. 7 depend on the model used to initialize the calculation. The result is almost identical to the one shown in Fig. 6-Glauber initialization favors lower value of viscosity, η/s ≈ 0.08, than KLN initialization, η/s ≈ 0.16. Since the uncertainties are smaller, this result was estimated to provide a limit 1/(4π) < η/s < 2.5/(4π) for the effective QGP viscosity. 96,101 For further discussion of evaluating η/s, see Ref. 96.

η/s(T )
Calculations discussed in section 3.4 assumed that the η/s-ratio is constant. We know no fluid where η/s would be temperature independent, and there are theoretical reasons to expect it to depend on temperature with a minimum around the phase transition. 104 Thus the temperature independent η/s is an effective viscosity, and its connection to the physical, temperature dependent, shear viscosity coefficient is unclear. What complicates the determination of the physical shear viscosity coefficient, is that the sensitivity of the anisotropies to dissipation varies during the evolution of the system. As studied in Ref. 105, and illustrated in Fig. 8 ( √ s NN = 200 GeV) v 2 (p T ) is insensitive to the value of η/s above phase transition, but very sensitive to its value in the hadronic phase. At the present LHC energy ( √ s NN = 2.76 TeV) the shear viscosity in the plasma phase does affect the final v 2 (p T ), but not more than the shear viscosity in the hadronic phase. It is only at the full LHC energy, √ s NN = 5.5 TeV where the viscosity in the plasma phase dominates, and dissipation in the hadronic phase has only a minor effect. Note that a change of the minimum value of η/s would clearly change v 2 (p T ) at all energies.
Additional factor complicating the determination of the temperature dependence of η/s is that the effect of viscosity on the anisotropies does not depend only on the ratio η/s, but also on the relaxation time τ π of the shear stress tensor. This is demonstrated in Fig. 9. If the minimum value of η/s is increased by a factor two, v 2 is reduced as expected, but if the relaxation time is also increased by a factor two, the effect of the increase in η/s is almost completely compensated. The interplay of relaxation time and shear viscosity was discussed also in Refs. 57 and 110, where it was found that to reproduce hybrid model results using viscous hydrodynamics only, it is not sufficient to increase η/s in the hadronic phase, but one should also change the relaxation time. So far we have seen that calculations with constant η/s require slightly larger value of η/s at LHC than at RHIC. 111,112 This is in line with the increase of (η/s)(T ) in high temperatures, but as shown in Ref. 112, one cannot uniquely constrain (η/s)(T ) by fitting the spectra and v 2 alone.

Fluctuations
In the average the matter formed in heavy-ion collisions has a smooth shape as indicated in Fig. 1, but in each event that is not true. Nuclei are not smooth  distributions of nuclear matter, but consist of individual nucleons, and thus the collision zone can be expected to have a highly irregular shape which fluctuates event-by-event, see Fig. 10. The initial state models (see sect. 2.1) can be generalized to so-called Monte-Carlo (MC) models to take this into account. In those models the initial Woods-Saxon distributions of nuclear matter are sampled to generate configurations of nucleons in nuclei, and the number of participants/collisions 4 or gluon production 12,13 in an event is calculated based on the positions of individual nuclei. It was shown already a while ago that the final observables differ whether one first averages the initial state, and evolves it hydrodynamically, or one evolves the event-by-event fluctuating initial states individually, and averages the results. 114,115 This became widely recognized only years later when it was realized that because of the irregular shape of each event, not only even, but also odd anisotropy coefficients v n are finite and measurable. 116 The higher coefficients are very helpful for extracting the transport coefficients from the data, since the larger the n, the more sensitive the coefficient v n is to viscosity, 117 see Fig. 11. Thus the study of fluctuations provides a way to distinguish different initializations, and first results for the p T -dependence of v 2 and v 3 (called triangular flow) seem to favor the MC-Glauber initialization. 118 It has been suggested that initial fluctuations provide also a way to circumvent our ignorance of the initial shape. 119 In most central collisions the anisotropies are entirely driven by fluctuations, which are better understood than the average shape in semi-central collisions. Thus evaluating η/s using v n in the most central collisions should be less sensitive to the model used to calculate the initial state. The preliminary result of such an analysis was that 0.07 < η/s < 0.43. 119 This value is again for effective viscosity which does not exclude the possibility that η/s could be even smaller or larger in some temperature region. For a further discussion of flow and viscosity, see Ref. 93.
In event-by-event studies it is not sufficient to reproduce only the average values of v n , but the fluctuations of the flow coefficients should be reproduced as well. The distributions of these fluctuations provide a way to constrain the fluctuation spectrum of initial state models independently of the dissipative properties of the fluid. As shown in Fig. 12, once the average v n has been scaled out, the distributions of these fluctuations, i.e., (v n − v n )/ v n or v n / v n , are almost independent of viscosity. The independence extends to other details of the evolution to such an extent, that the distributions of the fluctuations of initial anisotropies are good approximations of the measured distributions of v n , 120 and thus it is sufficient to compare the fluctuations of initial shape, ǫ n , to the observed fluctuations of v n . Unlike the recently developed IP-Glasma model, 121,122 neither MC-Glauber nor MC-KLN model seems to be able to reproduce the measured fluctuations. 123

Initialization
During last two years there have been major advances in modeling the initial state of hydrodynamical evolution. So called IP-Glasma model 121,122 is based on Color Glass Condensate and employs the IP-Sat (Impact Parameter dependent Saturation) model 124,125 of nucleon wavefunctions to generate fluctuating gluon fields in the initial collision, and the classical Yang-Mills dynamics to evolve these fields. [126][127][128][129][130][131] Unlike most of the initial state models, the IP-Glasma model includes the fluctuations of color charge in colliding nucleons too, not only the fluctuations of nucleon positions. Another advantage is that the model includes some of the pre-equilibrium dynamics of the gluon fields making the model less sensitive on the time when one switches to hydrodynamics. 111 Unfortunately the description is still incomplete: the present IP-Glasma model does not lead to a thermal system, but final thermalization still has to be assumed. Nevertheless, calculations using the IP-Glasma initialization reproduce both the fluctuations and the average values of v 2 , v 3 and v 4 , 111, 132 which makes this approach very promising. For a general overview of IP-Glasma, fluctuations and hydrodynamics, see Ref. 133.

What worries us
As described in the previous sections, fluid dynamics has been very successful in describing heavy ion collisions. However, at the time of this writing there are some data which may cause difficulties for the conventional fluid dynamical picture.

Photons
Unlike hadrons and partons, photons and leptons interact only electromagnetically, and hardly scatter at all after being produced. Thus the observed photon and lepton spectra contain contributions from all stages of a heavy-ion collision, and can be used to probe the early hot and dense stage. 134 The yield and p T -spectrum of photons in heavy-ion collisions is fairly well understood 135-139 -see also Ref. 140 and references therein-but the recent measurements of direct photon v 2 at RHIC 141 and LHC 142 have presented a puzzle: If photons are emitted during all stages of the evolution, and momentum anisotropy increases during the evolution, photon v 2 should be smaller than the hadron v 2 . Such a behavior is seen in the theoretical calculations, 143 but not in the data where the photon v 2 (p T ) is roughly equal to the observed pion v 2 (p T ). 141,142 The further refinements of the calculations, where event-by-event fluctuations 144 and viscosity have been taken into account, 145,146 have not been able to increase the v 2 to the observed level, see Fig. 13. So far the only approach to get close to the data has used parametrized expansion with very strong initial expansion and large photon production rates in the hadronic phase. 147 Thus the reproduction of the data remains a challenge to our understanding of microscopic production rates and/or expansion dynamics.

p + P b collisions
In p+P b collisions at √ s NN = 5.02 TeV, the measured dihadron correlations, 148-151 multiplicity and species dependence of average p T , 152, 153 and elliptic and triangular flows 149-151 all depict features easily explained by using hydrodynamics [154][155][156][157][158][159] (For a short summary see Ref. 160). Especially striking is the mass ordering of the p Tdifferential elliptic flow, 151 see Fig. 14. As discussed in section 3.1, this was taken as a strong indicator of the formation of a thermal system. However, it is questionable whether hydrodynamics is applicable to such a small system. Gradients are so large, that dissipative corrections should be of the same order than equilibrium pressure. As discussed in Ref. 161, large dissipative corrections are a problem even in P b + P b collisions, but in P b + P b collisions corrections are large only for a small fraction of the lifetime of the system, whereas corrections are large for a significant fraction of the lifetime in p + P b collision.
In the Color Glass Condensate framework similar correlations [162][163][164] and the mass hierarchy of average p T 165 arise as a result of gluon saturation in the proton and nuclear wavefunctions. How much of the observed behavior can be explained as such an initial state effect h , and how to differentiate initial state effects from hydrodynamical final state effects 161, 166 is at the time of this writing under intense study. If it turns out that the apparently hydrodynamical behavior in p + P b collisions can be explained as an initial state effect, then we may wonder whether the hydrodynamical behavior in P b + P b collisions could be explained as an initial state effect as well. On the other hand, if this is not possible, and fluid dynamics is the most viable description of p + P b collisions, then the properties of QCD matter are even more surprising than thought so far.

Summary
Fluid dynamics has been very successful in explaining the features of bulk, i.e. low p T , particle production in ultrarelativistic heavy-ion collisions. We have seen that in particular anisotropies of particle production can be explained if the rescatterings among particles are so frequent that the system is approximately thermal, that the equation of state of such a matter has many degrees of freedom and relatively hard, and that the shear viscosity coefficient over entropy ratio of the produced matter has very low value at some temperature. Providing further experimental constraints on the equation of state as well as figuring out the temperature dependence and the precise value of the minimum of η/s will require a lot of work. The present studies in event-by-event fluctuations and recent advances in modeling the pre-equilibrium processes in heavy-ion collisions are very helpful for this goal, but it is not yet clear how the elliptic flow of photons in A + A collisions and the p + P b collisions at LHC fit in the overall picture. After a few years' work we will know.