2.1 Viscous hydrodynamics

Viscous hydrodynamics are a widely used tool to describe the expansion of the QGP fireball and to study the soft hadron data for the heavy ion collisions at RHIC and the LHC [18-20, 61-73]. It solves the transport equations of energy momentum tensor and net charge current, which are written as

If the systems are very close to local equilibrium, the energy momentum tensor and the net baryon charge current can be decomposed as: Tμ ν=(e+p) uμuν-pgμν and Nμ=nuμ. Therefore, the fourteen variables in Tμν and Nμ are reduced to six independent unknowns: the energy density, e, pressure, p, and net baryon density, n, and 3 independent components in the four velocity, uμ. The relativistic hydrodynamics are then simplified as ideal hydrodynamics. With an additional input, the equation of state (EoS), p=p(n,e), and the chosen initial and final conditions, the ideal hydrodynamic equations can be numerically solved to simulate the evolution of the bulk matter for the relativistic heavy ion collisions [10].

For a near equilibrium system, one needs to implement the relativistic viscous hydrodynamics (or the so-called relativistic dissipative fluid dynamics). In the Landau frame, Tμν and Nμ are expressed as: Tμ ν=(e+p+∏) u^μu^ν-(p +∏)gμν+πμ ν,$N^{\mu }=n{u}^{\mu }-\frac{n}{e+p}{q}^{\mu }$. Here, πμ ν is the shear stress tensor, ∏ is the bulk pressure, and qμ is the heat flow. From the 2nd law of thermal dynamics or from the kinetic theory, one could obtain the viscous equations of πμ ν, ∏, and qμ, which are expressed as [74, 75]:

where ${\Delta }^{\mu \nu }=g^{\mu \nu }-u^{\mu }{u}^{\nu }$, ${\nabla }^{\langle \mu }{u}^{\nu \rangle }=\frac{1}{2}({\nabla }^{\mu }{u}^{\nu }+{\nabla }^{\nu }{u}^{\mu })-\frac{1}{3}{\Delta }^{\mu \nu }{\partial }_{\alpha }{u}^{\alpha }$, and $\theta =\partial \cdot u$. η is the shear viscosity, ζ is the bulk viscosity, λ is the heat conductivity, and τπ, τ∏, and τq are the corresponding relaxation times.

The above Israel-Stewart formalism can also be obtained from the kinetic theory [76-80] or from the conformal symmetry constraints [77] ⁴. These different derivations give different higher order terms for the 2nd order viscous equations. In general, the contributions of the higher order terms are pretty small or even negligible for a hydrodynamic evolution with small shear and bulk viscosity, which will not significantly influence the final flow observables ⁵.

The equations of state (EoS):

Besides these hydrodynamic equations, one needs to input an EoS to close the system for the numerical simulations or analytical solutions. Currently, many groups use a state-of-the-Art EoS, called s95p-PCE, which combines a parameterized/tablated lattice EoS for the baryon free QGP phase with a hadronic EoS with effective chemical potentials for the partially chemical equilibrium hadronic phase [100, 101]. Ref. [100] also compared the hydrodynamic calculations using various equations of state constructed with a different speed of sound, which found that the spectra and elliptic flow are only slightly influenced by the inputting EoS. The main uncertainties of the hydrodynamic calculations come from the initial conditions, which will be introduced and discussed below.

Initial conditions:

The initial condition is a necessary input for the hydrodynamic simulations. As an open issue related to the creation and thermalization of the QGP, it brings some uncertainties, more or less, for many flow observables in the hydrodynamic calculations. There are many types of initial condition models developed by different groups. The traditional Glauber model assumes zero transverse velocity at the starting time and constructs the initial entropy/energy density profiles from a mixture of the wounded nucleon and binary collision densities [102]. The KLN model treat the degrees of freedom of the initial systems as gluons and calculate the initial density profiles from the kT factorization formula [103]. In the later developed Monte-Carlo versions, called (MC-Glauber and MC-KLN) [104-106], the event-by-event fluctuations are built through the positions fluctuations of individual nucleons inside each colliding nuclei. For the AMPT initial conditions, the initial profiles are constructed from the energy and momentum decompositions of individual partons, which fluctuate in both momentum and position space [107-109]. With an additional Gaussian smearing factor, the fluctuation scales related to the energy decompositions become changeable, which helps to balance the initial eccentricities at different order. As a successful initial condition model, IP-Glasma [110] includes both the nucleon position fluctuations and the color charge fluctuations. It uses the IP-Sat model to generate the wave-functions of high energy nuclei/nucleon and then implements a classical Yang-Mills dynamic to simulate the pre-equilibrium evolution of the early glasma stage. Another successful initial condition model in EKRT [111, 112] combines the PQCD minijet production with the gluon saturation and generates the energy density profiles after a pre-equilibrium Bjorken free streaming. The recently developed TRENTo model [113] is a parametric initial condition model based on the eikonal entropy deposition via a reduced thickness function. With an introduced entropy deposition parameter, the TRENTo model could reproduce the initial eccentricities of various initial condition models that belong to different classes, such as MC-KLN, MC-Glauber, EKRT, IP-Glasma, and etc..

Many initial condition models neglect the initial flow from the pre-equilibrium stage. Recently, the effects of pre-equilibrium evolution have been estimated in Ref [114] through evolving the free streaming particles from MC-Glauber and MC-KLN models, which demonstrated that such pre-equilibrium dynamics significantly increase the initial flow and reduce the initial spacial eccentricities. More sophisticated investigations on pre-equilibrium dynamics can be, in principle, carried on within the framework of dynamical models like EPOS [115], AMPT [107-109], EKRT [111, 112], IP-Glasma [110], URQMD [116, 117], and etc.. After matching the energy-momentum tensor at a switching point, one could principally obtain 3+1-d fluctuating profiles of initial energy density and initial flow for the succeeding hydrodynamic simulations. However, many past studies focus on the initial state fluctuations on the transverse plane, which neglect the fluctuation patterns along the longitudinal direction. The AMPT + ideal hydrodynamic simulations [108] demonstrate that evolving early hot spots in the longitudinal directions could dissipate part of the transverse energy, which leads to a suppression of the final flow anisotropy. Recently, the IP-Glasma model has been extended to three dimension with the explicit small x evolutions of the gluon distributions [118]. Although the related energy momentum tensors can be, in principle, used in the succeeding hydrodynamic simulations, additional works are still required to further extend the distributions to the large rapidity regime with the consideration of large x effects.

Freeze-out / decoupling :

Pure hydrodynamic simulations assume free-streaming hadrons directly emitted from a decoupling surface defined by a constant temperature, energy density, or other kinetic variables[10, 61]. The momentum distributions of various emitted thermal hadrons can be calculated with the Cooper-Frye formula [119] using the freeze-out information on the freeze-out surface (For the details of the Cooper-Frye formula, please refer to [10, 119] as well as the following Section 2.1 for details). With the corresponding decay channels, the unstable hadron resonances delay into stable ones with some momentum distributions that can be further analyzed and compared with the experimental data. In the constant temperature decoupling scenario, the decoupling temperature, T_dec strongly depends on the EoS and other hydrodynamic inputs. For s95p-PCE, T_dec is generally set to 100-120 MeV in order to fit the mean p_T of various hadrons with a sufficient build up of the radial flow [10, 101].

2.2 Hybrid models

A hybrid model matches the hydrodynamic description of the QGP fluid to a hadron cascade simulation for the evolution of the hadron resonance gas at a switching temperature near Tc. The early ideal hydrodynamics + hadron cascade hybrid model simulations have showed that the hadronic matter is highly viscous, which largely suppress the elliptic flow when compared with the pure hydrodynamic calculations with a partially chemical equilibrium EoS [120]. Motivated by this, different groups have extended 2+1-d or 3+1-d viscous hydrodynamic simulations with a hadronic afterburner [121-123]. Such hybrid models give a more realistic description for the hadronic evolution of the QGP fireball, which also naturally imprint the off-equilibrium chemical and thermal freeze-out procedures of various hadron species.

The key component of a hybrid model is the particle event generator that converts the hydrodynamic output on the switching hyper surface into various hadrons for the succeeding hadron cascade simulations. More specifically, such Monte Carlo event generators produce particles with specific momentum and position information according to the differential Cooper-Frye formula [121]:

where fi(x,p) is the distribution function of hadron species, i, gi is the corresponding degeneracy factor, and d³σμ(x) is a surface element on the hyper-surface, ∑, e.g., defined by a constant switching temperature T_sw. Generally, the switching temperature, T_sw, is set to around 160 MeV, which is close to the phase transition temperature of the QCD matter at zero chemical potential [124]. For a viscous QGP fluid, the distribution function f(x,p) includes an ideal part and an off-equilibrium part f=f₀+δf, where δ f generally takes the form: $\delta f=f_0(1\mp f_0)\frac{p^{\mu }{p}^{\nu }{\pi }_{\mu \nu }}{2T^2\left(e+p\right)}$[64-69]⁶.

After converting the fluid into many individual hadrons of various species, the hybrid model implements a hadron cascade model to simulate the microscopic evolution of the hadron resonance gas. The hadron cascade model, for example, the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) [127, 128] solves a large set of Boltzmann equations for various hadron species:

where fi(x,p) is the distribution function and Ci(x,p) is the collision terms for the hadron species, i. With such equations, the hadron cascade model propagates various hadrons with classical trajectories, together with the elastic scatterings, inelastic scatterings, and resonance decays. After all the collisions and decays cease, the system stops evolution and outputs the information of the produced hadron, which can be further analyzed to compared with the corresponding experimental data [127, 128].

Compared with the early hydrodynamic calculations, the hybrid model improves the description of the hadronic evolutions and the decoupling procedure, which leads to a nice description of the flow harmonics of identified hadrons, especially for the mass-splitting between pions and protons [129, 130]. Meanwhile, the imprinted baryon-antibaryon ($B-\overline{B}$) annihilations in the hadronic cascade sector also largely reduce the production of proton and antiproton, which helps to achieve a nice fit of particle yields of various hadron species [129, 131].

2+1-d vs 3+1-d model:

For hydrodynamics or hybrid models, the 2+1-d simulations with a longitudinal boost invariance are more computationally efficient than the full 3+1-d simulations. Before 2010, many developed viscous hydrodynamic codes are (2+1)-dimensional using the Bjorken approximation [64-72]. The published VISHNU code is also basically a (2+1)-d hybrid code since it implements the (2+1)-d viscous hydrodynamic simulations for the evolution of the QGP phase. Although the succeeding UrQMD afterburner are (3+1)-dimensional, the longitudinal boost invariance is still approximately conserved at mid-rapidity after the hadronic evolution [121]. Recently, several groups [73, 132-136] further developed the full (3+1)-d viscous hydrodynamics or hybrid models without a longitudinal boost invariance. Such full (3+1)-d simulations could provide full space-time evolution profiles for the EM and hard probes. They can also be widely used to investigate the longitudinal fluctuations, to study the physics for asymmetric collision systems, such as p+Pb, d+Au, and Cu+Au, etc, and to provide more realistic calculations / predictions for the heavy ion collisions at lower collision energies.

2.3 Event-by-event simulations

As introduced in Sect.2.1, the initial profiles of the created QGP fireball fluctuate event-by-event, which leads to the final state correlations and collective flow for the nucleus-nucleus collisions at RHIC and the LHC [11-13]. For computational efficiency, the early hydrodynamics or hybrid model simulations input smooth initial profiles obtained through averaging a large number of events generated from some specific fluctuating initial conditions and then implement the so-called single-shot simulations. An alternative approach is the event-by-event simulations, which simultaneously run a large number of simulations with the input of individually fluctuating initial profiles. Past research has shown, due to the the approximate linear hydrodynamic response v₂ ∝ε₂ and v₃ ∝ ε₃, the elliptic and triangular flow can be nicely described by the single-shot hydrodynamic simulations with properly chosen initial conditions and well tuned parameter sets. However, the single shot simulation fails to describe other higher order flow harmonics due to the mode coupling effects. Furthermore, some flow observables, such as event-by-event flow harmonics [29, 30], the event-plane correlations [31, 32], and the correlations between different flow harmonics [33-37], etc., can not be directly calculated by the single-shot hydrodynamics or hybrid models, which are required to implement the event-by-event simulations (please also refer to Sect. 5 for details).

Since 2010, many groups have developed event-by-event hydrodynamics / hybrid models to study the initial fluctuations, hydrodynamic response, and the corresponding final state correlations [30, 32, 107-112, 137-141]. In general, such event-by-event simulations is computationally expansive. For example, the iEBE-VISHNU simulations for the correlations between flow harmonics have used 30000 CPU hours in the Tianhe-1A National Supercomputing Center in Tianjin, China. Recently, the OSU-Kent group has developed the massively parallel simulations for 3+1-d viscous hydrodynamics on graphics processing units with CUDA and demonstrated that such GPU simulations are approximately two orders of magnitude faster than the corresponding simulations from CPU [142]. With the development of computer science and the reduced cost of GPU, the GPU-based simulations will become a popular trend for the massive hydrodynamic calculations in the future.

4.1 Semi-quantitative extractions of the QGP shear viscosity

The hydrodynamic calculations from different groups have shown that the flow harmonics are sensitive to the QGP shear viscosity, η/s, which can be used to study the transport properties of the hot QCD matter [18-20, 61-73]. Around 2008, the INT group made an early extraction of the QGP shear viscosity from the integrated and differential elliptic flow data in 200 A GeV Au–Au collisions using the 2+1-d viscous simulations with optical Glauber and KLN initializations [64, 65]. They found these two initial conditions bring large uncertainties for the extracted value of η/s around O(100%). However, it is not reliable to directly read the value of η/s from the direct model to the data comparison since their model calculations neglect the high viscous and even off equilibrium hadronic evolution, which only treat such stages as a pure viscous fluid expansion with both chemical and thermal equilibrium. Ref. [65, 152] further estimated the effects from the late hadronic evolution and concluded that the extracted value of the specific QGP shear viscosity (η/s)_QGP, can not exceed an upper limit of around $5\times \frac{1}{4\pi }$.

For a realistic description of the evolution and decoupling of the hadronic matter, the OSU-LBL group further developed the VISHNU hybrid model [121] which combines the 2+1-d viscous hydrodynamics with a hadron cascade model- UrQMD and then made a semi-qualitative extraction of the QGP shear viscosity from the integrated elliptic flow data in 200 A GeV Au–Au collisions [148, 153]. Figure 1 shows the eccentricity-scaled integrated elliptic flow calculated from VISHNU with MC-Glauber and MC-KLN initial conditions together with a comparison with the corrected experimental data with the non-flow and fluctuation effects removed [149]. From Fig. 1, one finds $\frac{1}{4\pi }<(\eta /s{)}_{\text{QGP}}<2.5\times \frac{1}{4\pi }$, where the main uncertainties of the extracted (η/s)_QGP still come from the undetermined initial conditions. Meanwhile, the corresponding VISHNU simulations with both MC-Glauber and MC-KLN initial conditions could nicely describe the p_T-spectra and differential elliptic flow harmonics, v₂(p_T), for all charged and identified hadrons at various centrality bins in 200 A GeV Au–Au collisions [153]. Compared with the early extractions in Ref. [64], the precision of the extracted value of (η/s)_QGP is largely increased due to a better description of the highly viscous hadronic stage.

Fig. 1.Full-screen View

(Color online) Eccentricity-scaled elliptic flow as a function of final multiplicity per area. The theoretical results are calculated from the VISHNU hybrid model calculations with MC-Glauber (left) and MC-KLN (right) initial conditions [148]. The experimental data are taken from Ref. [149].

In Ref. [154], the VISHNU simulations were further extrapolated to the LHC energies, which systematically investigated the soft hadron data in the 2.76 A TeV Pb–Pb collisions. The related calculations have shown that with the same (η/s)_QGP extracted at top RHIC energies, VISHNU slightly over-predicts the ALICE flow data at the LHC. After slightly increasing (η/s)_QGP (for the MC-KLN initial conditions, (η/s)_QGP increases from ∼0.16 to ∼ 0.20), VISHNU achieves a better description of the elliptic flow of all charged hadrons at various centralities [154].

Many of the early hydrodynamic or hybrid model simulations (includes these 2+1-d hydrodynamic and VISHNU calculations mentioned above) [121, 129, 130, 148, 153-155] belong to the category of single-shot simulations, which input smooth initial energy/entropy profiles from early initial condition models or input some smoothed profiles obtained from averaging millions of events from some specific fluctuating initial condition models. Correspondingly, the effects from initial state fluctuations are neglected. Around 2012, the Mcgill group further developed event-by-event 3+1-d viscous hydrodynamic simulations with the IP-Glasma pre-equilibrium dynamics (MUSIC + IP-Glasma) and calculated the flow harmonics at different orders at the RHIC and the LHC [30]. Figure 2 shows the integrated and differential vn(n= 2...5) of all charged hadrons in the 2.76 A TeV Pb–Pb collisions. Impressively, these different flow harmonics data are nicely described by the MUSIC simulations with η/s=0.2 or a temperature dependent η/s(T) at various centralities. Meanwhile, their simulations also show the averaged QGP viscosity is slightly larger at the LHC than at the RHIC, as found in Ref. [154]. Compared with the VISHNU simulations [148, 153, 154], these MUSIC calculations are purely hydrodynamic, which does not specially treat the hadronic evolution with a hadronic afterburner. However, the main results will not be significantly changed since the flow harmonics at the LHC energies are mainly developed (or even reach saturation) in the QGP phase.

Fig. 2.Full-screen View

(Color online) Root-mean-square anisotropic flow coefficients ${\langle v_n^2\rangle }^{1/2}$ and vn(p_T) in the 2.76 A TeV Pb–Pb collisions. The theoretical curves are calculated from MUSIC with IP-Glasma initial conditions [30]. The experimental data in the left and right panels are measured by the ALICE collaboration [24] and the ATLAS collaboration, respectively.

For the hydrodynamic simulation with IP-Glasma initial conditions, a balanced initial eccentricities at different orders are generated at the beginning, which helps to achieve a simultaneous fit of the elliptic flow, triangular flow, and other higher order harmonics. In contrast, the hydrodynamic calculations with either Mc-Glauber or Mc-KLN initial conditions fail to simultaneously describe all the flow harmonics, vn, at different orders (n= 2... 5), although they can nicely fit the elliptic flow data with a well-tuned QGP shear viscosity. Therefore, these higher-order flow harmonic measurements disfavor these two initial conditions, which also motivated the later developments of other initial condition models. In short, the extracted value of the QGP viscosity may be largely influenced by the initial conditions used in the hydrodynamic calculations. Meanwhile, higher order flow harmonics as well as other flow observables (please also refer to Sect. 5 for details) could put straight constrains for the initial condition models and for the extracted value of the QGP shear viscosity.

Besides the flow data of all charged hadrons, the flow harmonics of identified hadrons could reveal more information on the hadronic evolution of the hot QCD matter and provide additional tests for extracted values of the QGP transport coefficients obtained from the soft hadron data of all charged hadrons. Ref. [129] and [130] have shown, for the extracted constant QGP shear viscosity obtained from the elliptic flow in 2.76 A TeV Pb–Pb collisions, the VISHNU hybrid model could nicely describe the differential elliptic flow data of pions, kaons, and protons [129, 130]. Meanwhile, it could also roughly fit the elliptic flow data of strange and multi-strange hadrons (Λ, Ξ, and Ω) measured at the LHC [155]. Recently, the ALICE collaboration further measured the higher order flow harmonics of identified hadrons in the 2.76 A TeV Pb–Pb collisions, which showed that the triangular and quadratic flow harmonics of pions, kaons, and protons present similar mass ordering as the case for the elliptic flow [156]. In Ref. [109], the PKU group implemented the iEBE-VISHNU hybrid model with the AMPT initial conditions to investigate the flow harmonics of identified hadrons vn(p_T) (n= 2,3,4) at the LHC. After tuning the Guassian smearing factor for initial energy decompositions and the QGP shear viscosity, the differential vn of all charged hadrons can be nicely described by the iEBE-VISHNU simulations. As shown in Fig. 3, iEBE-VISHNU also nicely describes the vn data of pions, kaons, and protons, especially when reproducing correct mass-orderings for these different flow harmonics. Ref. [109] also showed the pure hydrodynamic simulations do not generate enough mass-splittings between the vn of pions and protons. The late hadronic evolution in the iEBE-VISHNU re-distributes the anisotropy flow to various hadron species through the microscopic hadronic scatterings which enhances the vn mass splitting between pions and protons and leads to a nice description of the experimental data [109].

Fig. 3.Full-screen View

(Color online) vn(pT) (n=2,3,4) of pions, kaons, and protons in 2.76 A TeV Pb–Pb collisions, calculated from iEBE-VISHNU with AMPT initial conditions [109]. The experimental data are taken from the ALICE paper [150, 151].

4.2 Quantitative extractions of the QGP shear and bulk viscosity with massive data evaluations

For the flow calculations and predictions at RHIC or at the LHC, most of the hydrodynamics/hybrid model simulations, with different types of initial conditions, input a constant value of the specific QGP shear viscosity and neglect the effects of bulk viscosity. The early model calculations also revealed that the averaged QGP shear viscosity changes with the collision energies, which is slightly larger at the LHC than at RHIC [19, 123, 131, 154]. It is thus very important to extract a temperature-dependent QGP shear viscosity, (η/s)_QGP(T), from the massive soft hadron data in relativistic heavy ion collisions. For the purposes of massive data evaluations, the Livermore group developed the CHIMERA algorithm (a comprehensive heavy ion model evaluation and reporting algorithm) and extracted the initial temperature and the QGP shear viscosity from a simultaneous fit of the p_T spectra, elliptic flow, and HBT radii in 200 A GeV Au + Au collisions [165]. Note that these early massive hydrodynamic simulations around 2012 assume the QGP shear viscosity is a constant value and the bulk viscosity is zero together with an input of the traditional MC-Glauber initial condition, which has been ruled out by some later flow measurements.

To avoid the limitations of simultaneously tuning multiple free parameters in the early work [165], the Duke-OSU group implemented the Bayesian method to the event-by-event hybrid model simulations [166] and then quantitatively estimated the properties of the QGP through a multi-parameter model to data comparison, using the parametric TRENTo initial conditions [167]. Using the newly developed massive data evaluating techniques, the global fitting of the multiplicity, transverse momentum, and flow data at the LHC constrain the free parameters in the TRENTo model, which also give an extracted temperature-dependent specific shear viscosity and bulk viscosity.

Figure 4 shows the estimated temperature dependent shear viscosity, (η/s)(T), from the DUKE-OSU group, obtained from the massive data fitting in 2.76 A TeV Pb+Pb collisions. The blue line is the median with a blue band showing the 90% credible region. Correspondingly, a nonzero bulk viscosity with a peak near the QCD phase transition has also been extracted simultaneously (For details, please refer to [167]). With these extracted QGP transport coefficients and other extracted most probable parameters, the event-by-event hybrid simulations give an excellent overall fit for the multiplicities and mean p_T of all charged and identified hadrons and the integrated vn (n= 2,3,4) of all charged hadrons from the most central collisions to the peripheral collisions in Pb–Pb collisions at the LHC, as shown in Fig. 5.

Figure 4.Full-screen View

(Color online) Estimated temperature dependence of the shear viscosity (η/s)(T) above the QCD phase transition (for Tc> 154 MeV), obtained from a multi-parameter model to data comparison [167].

Figure 5.Full-screen View

(Color online) Multiplicities, mean p_T of all charged and identified hadrons and the integrated vn (n= 2,3,4) of all charged hadrons in 2.76 A TeV Pb–Pb collisions, calculated from event-by-event hybrid model with the high-probability parameters extracted from the massive data fitting [167]. The data are from the ALICE experiment [171, 24].

Note that this extracted η/s(T), within the uncertainty band, is compatible with the well-known KSS bound, $\eta /s<1/4\pi$[168-170], which also supports several early semi-quantitative extractions of the QGP viscosity at RHIC and the LHC. For example, the extracted specific QGP viscosity, $\frac{1}{4\pi }<(\eta /s{)}_{QGP}<2.5\times \frac{1}{4\pi }$, from the VISHNU calculations with MC-Glauber and MC-KLN initial conditions [148, 153] and the implemented η/s = 0.095 (with the same bulk viscosity parametrization) in the MUSIC simulations with the IP-Glasma initialization [164] are both consistent with this quantitative extracted results from the DUKE-OSU collaborations. The early EKRT viscous hydrodynamic calculations for the flow data at RHIC and the LHC also prefer a temperature-dependent η/s(T) with a positive slope [112].

Compared with the early extraction of the QGP viscosity with specific initial condition, Ref. [167] implemented the parametric TRENTo model that could smoothly interpolate among various initial condition schemes through tuning the related parameters. It is thus an ideal initial state model for the massive model-to-data comparison, which helps to make a simultaneous constraint for the initial conditions and the QGP transport coefficients. It was found that the initial entropy deposition from the constrained TRENTo model with fixed parameters is approximately proportional to the geometric mean of the participant nuclear densities, which gives similar scaling to the successful EKRT and IP-Glasma initial conditions.

In Ref. [172], the Bayesian statistical analysis was extended to the massive data fitting in Au–Au collisions at $\sqrt{s_{\text{N}N}}$=19.6, 39 and 62.4 GeV. It was found that the extracted constant QGP specific shear viscosity, η/s, decreases with the increase of collision energy, which shows a similar result obtained from the early hybrid model calculations [123]. In the future, a combined massive data fitting at RHIC (including BES) and the LHC will give more precise temperature-dependent transport coefficients of the QGP.

5.1 Event-by-event vn distribution:

The flow harmonics, vn, are generally measured within a base of the event-average, which mainly reflects the hydrodynamic response to the averaged initial eccentricity coefficients, εn, within some centrality bin. With a large amount of particles produced per event at the LHC, a direct measurement of the event-by-event vn distribution becomes possible, which provides more information on the initial state fluctuation and the underlying probability density function. Around 2012, the ATLAS Collaboration made the first measurement of the event-by-event distributions of vn (n = 2, 3, 4) in Pb–Pb collisions at$\sqrt{s_{\text{N}N}}$ =2.76 TeV [29]. Figure 6 shows the MUSIC hydrodynamic calculations nicely describe the ATLAS data with the IP-Glamsa initial conditions. It also shows, for n=2 and 3, the rescaled vn/< vn > distributions mostly follow the εn/< εn> distributions from the initial state, which are not sensitive to the details of the hydrodynamic evolution [30]. Due to the mode couplings effects for higher flow harmonics, the distributions of v₄ / < v₄> are not necessarily follow ε₄/< ε₄ >. The hydrodynamic evolution balances the distributions of v₄ / < v₄> making a nice description of the experimental data. In Ref. [29], the measured vn distributions were compared with the εn distributions from MC-Glauber and MC-KLN models, which demonstrated certain deviations between model and data for most of the centrality classes. The vn distributions thus provide strong constrains on the initial state models, which do not favor the MC-Glauber and MC-KLN initial conditions.

Figure 6.Full-screen View

(Color online) Scaled event-by-event distributions of vn (n= 2,3,4) from the MUSIC simulations with the IP-Glasma initial conditions [30], together with a comparison with the ATLAS data [29].

The ATLAS measurements can also be used to examine the underlying p.d.f. of the vn distributions. The most popular parameterizations are the Bessel-Gaussian distributions [178]:

where v₀ is the anisotropic flow from the reaction plane, Ψ_RP, and σ is the anisotropic flow fluctuation. It was reported that the Bessel-Gaussian distribution could nicely describe the v₂ distributions for mid-central collisions [178, 179]. Without the constraint of ε₂ 1 for each event it is not expected to work well in peripheral collisions [180]. To fix this problem, a new function, named “Elliptic Power" distribution, was proposed in Ref. [180], which is expressed as:

where α quantifies the fluctuations and v₀ has the same meaning as the Bessel-Gaussian parameterizations. As a promising candidate of underlying p.d.f. of vn distribution, the Elliptic-Power function can nicely describe the event-by-event v₂ and v₃ distributions [180, 181]. However, it can not give an equally nice fitting for these distributions of higher flow harmonics (n≧ 4), which are largely influenced by the non-linear hydrodynamic response. For details, please refer to [180, 181].

5.2 De-correlations of the flow-vector Vn:

Recently, it was realized that the produced particles at different transverse momentums, p_T, and rapidity y do not share a common flow angle or event plane. Such transverse momentum and rapidity dependent flow angles fluctuate event-by-event, which also breaks the factorizations of the flow harmonics [38, 39]. To evaluate the de-correlations of the flow-vector, especially on the transverse momentum dependence, two new observables, vn2/vn[2] <, and the factorization ratio, rn, have been proposed, which are defined as:

where $v_n^a$, ${\Psi }_n^a$(or $v_n^b$, ${\Psi }_n^b$) are the n^th-order flow harmonics and the flow angle at the transverse momentum$p_{\text{T}}^a$(or $p_{\text{T}}^b$). The p_T dependent fluctuations of the flow angle and magnitude make vn{2}/vn[2] < and rn deviated from 1. As shown in Fig. 7, these deviations from unity have already been observed in experiment and qualitatively described by the related hydrodynamic calculations [38], which indicates the existence of the p_T dependent fluctuations of flow angle and magnitude.

Figure 7.Full-screen View

(Color online) The ratio vn2/vn[2] at various centralities in 2.76 A TeV Pb–Pb collisions. The theoretical lines are calculated from VISH2+1 with MC-Glauber and MC-KLN initial conditions [38], the experimental data are measured by the ALICE collaborations [177].

The fluctuations in the longitudinal direction have also been investigated both in experiment and in theory [40, 183-186]. Ref. [184] found that the the final state de-correlations of the anisotropic flows in different pseudo rapidity regime is associated with the spatial longitudinal de-correlations from the initial state. It also predicted a larger longitudinal decorrelation at RHIC than the ones at the LHC, which provide opportunities to further study the longitudinal fluctuation structures of the initial stage.

5.3 Event-plane correlations:

The correlations between different flow vectors could reveal more information on the initial state fluctuations and the hydrodynamic response [187]. In Ref. [31], the ATLAS Collaboration has measured the event-plane correlations among two or three event-plane angles, <cos(cnnΨn + cmmΨm) > and <cos(cnnΨn + cmmΨm + chhΨh) >, in 2.76 A TeV Pb–Pb collisions and observed several different centrality-dependent trends for these correlators. It was also reported that the MC-Glauber model, which only involves the correlations from the initial state, can not reproduce the trends for many of these correlators [31]. Using event-by-event hydrodynamics with MC-Glauber and MC-KLN initial conditions, Qiu and Heinz have systematically calculated the event-plane correlations and demonstrated that the hydrodynamic evolution is essential for an overall qualitative description of various flow angle correlations [32]. Figure 8 presents the model to data comparisons for several selected correlation functions which shows that although correlation strength is sensitive to the initial conditions and the QGP shear viscosity, hydrodynamics successfully reproduces the centrality-dependent trend of these event-plane correlations. In contrast, the correlations of the initial eccentricity plane show large discrepancies with the measured and calculated event correlations of the final produced particles, including magnitudes, qualitative centrality dependence, and even in signs [32]. In Ref. [31, 188], it was found that the AMPT simulations are also able to roughly reproduce the ATLAS data with well tuned parameters. These different model calculations, involving final state interactions [31, 32, 188], demonstrate that the observed event-plane correlations are not solely driven by the initial geometry, but are largely influenced by the complicated evolution of the QGP fireball.

Figure 8.Full-screen View

(Color online) Centrality dependent event-plane correlations calculated from event-by-event VISH2+1 hydrodynamic simulations with MC-Glauber and MC-KLN initial conditions [32]. The data are measured by the ATLAS collaborations [31].

Using a nonlinear response formalism, Ref. [182] calculated the event plane correlations from the initial energy density but expanded with the cumulants method,which roughly reproduces the centrality-dependent trends of several selected correlations. It is also found that the non-linear response of the medium has a strong influence on these related correlators. As shown in Fig. 9, the linear response alone is not able to describe the <cos(4(Ψ₂ - Ψ₄)) > and <cos(2Ψ₂ + 3Ψ₃ - 5Ψ₅) > correlators, while a good description of the data can be achieved after combining the contributions of both linear and non-linear response.

Figure 9.Full-screen View

(Color online) The separate contributions from linear, non-linear, and combined response to the event-plane correlations [182] together with a comparison with the ATLAS data [31].

5.4 Correlations of flow harmonics:

Besides the event-plane correlations, the correlations between different flow harmonics are other important observables closely related to the corrections of the flow vectors, that could further reveal the initial state correlations and the hydrodynamic response. Using the Event-Shape Engineering (ESE) [189], the ATLAS Collaboration firstly measured the correlations between flow harmonics based on the 2-particle correlations and found that v₂ and v₃ are anti-correlated and v₂ and v₄ are correlated [33] ⁸. Recently, a new observable, called Symmetric Cumulants SCv(m, n), was proposed as an alternative approach to measure the correlations between different flow harmonics. It is defined as$S{C}^v(m\mathrm{,}n)=<v_m^2{v}_n^2>-<v_m^2><v_n^2>$ and can be measured by the multi-particle cumulant method. The related Monte-Carlo model simulations imply that SCv(m, n) is insensitive to the non-flow effects [34]. Besides, SCv(m, n) is independent on the symmetry plane correlations by design [146].

Figure 10 (left) shows the centrality dependent symmetric cummulants SCv(4, 2) and SCv(3, 2) in 2.76 A TeV Pb–Pb collisions, measured from ALICE [34] and calculated from the EKRT event-by-event hydrodynamics [112]. The positive values of SCv(4, 2) and negative values of SCv(3, 2) are consistent with the early observation from ATLAS [33], which also illustrates that v₂ is anti-correlated with v₃, but is correlated with v₄. A comparison between the model calculations and the experimental data in Fig. 10 also shows that, although hydrodynamics could successfully reproduce the integrated flow harmonics, vn, it can only qualitatively, but not quantitatively, describe the correlations between these harmonics.

Figure 10.Full-screen View

(Color online) The centrality dependence of symmetric cumulants SC(4,2) and SC(3,2) in 2.76 A TeV Pb–Pb collisions [34].

In Ref. [36], the symmetric cumulants SCv(m, n) and other related observables have been systematically calculated by the event-by-event viscous hydrodynamics VISH2+1 with a focus on investigating the influences from different initial conditions and QGP shear viscosity. Like the case of the early EKRT hydrodynamic simulations, all of these VISH2+1 simulations with MC-Glauber, MC-KLN, and AMPT initial conditions could capture the sign and centrality dependence of SCv(4, 2) and SCv(3, 2), but not be able to archive a simultaneous quantitative description of these two symmetric cumulants for all centrality intervals. Compared with the individual flow harmonics v₂ and v₃, the symmetric cumulants SCv(4,2) and SCv(3,2) are more sensitive to the details of the theoretical calculations. Ref. [36] also predicted other symmetric cumulants SCv(5, 2), SCv(5, 3), and SCv(4, 3) and found that v₂ and v₅ are correlated, v₃ and v₅ are correlated, v₃ and v₄ are anti-correlated for various centralities.

In order to get rid of the influences from individual flow harmonics, it was suggested to normalize SCv(m,n) by dividing the products $<v_m^2><v_n^2>$[34]. Fig. 10 (right) and Fig. 11 (a,b,c,g,h) plot the normalized symmetric cumulants NSCv(n,m) (NSCv(n,m) = SCv(n,m)/$<v_n^2><v_m^2>$) in 2.76 A TeV Pb–Pb collisions. NSCv(4,2) exhibits a clear sensitivity to the initial conditions and the η/s(T) parameterizations, which could provide additional constrains for the initial geometry and the transport coefficients of the hot QCD matter. In contrast, NSCv(3,2) is insensitive to the detailed setting of η/s and the used initial conditions. Fig. 11 also shows that the values of NSCv(3,2) is compatible to the ones of NSCε(3,2) from the initial state due to the linear response of v₂ (v₃) to ε₂ (ε₃). Note that these different NSCv(3, 2) curves in Fig. 11 (g) are almost overlap with each other, which also roughly fit the normalized ALICE data. In contrast, the predicted NSCv(4, 2), NSCv(5, 2), and NSCv(5, 3) are sensitive to both initial conditions and η/s. Due to the nonlinear hydrodynamic response, NSCv(4, 3) does not necessarily follow the sign of NSCε(4, 3) for some certain initial conditions.

Figure 11.Full-screen View

(Color online) The centrality dependence of normalized symmetric cumulants NSC(m,n), and the corresponding normalized symmetric cumulants of the initial eccentricity coefficients NSCε(m,n) in 2.76 A TeV Pb–Pb collisions, calculated from event-by-event VISH2+1 simulations with MC-Glauber, MC-KLN, and AMPT initial conditions [36].

In a recent work [35], the NSCv(m,n) are expressed in terms of the symmetry plane correlations and moments of v₂ and v₃. Considering the relative flow fluctuations of v₃ are stronger than v₂, one expects smaller values for NSCv(5,2) compared to NSCv(5,3), as shown in Fig. 11. On the other hand, it was predicted that NSCv(m,n) involving v₄ and v₅ increases with η/s in the same way as the symmetry plane correlations [190, 191], which qualitatively agrees with the results in Fig. 11 from most central collisions to semi-peripheral collisions.

As discussed above, the low flow harmonics, v₂ or v₃, is mainly determined by a linear response to the initial eccentricity ε₂ or ε₃, while higher flow harmonics(vn with n> 3) not only contains the contributions from the linear response of the corresponding εn, but also has additional contributions from lower order initial anisotropy coefficients. These additional contributions are usually called non-linear response of higher flow harmonics [192, 193]. In Ref. [35], it was proposed that a direct connection between the symmetry plane correlations and the flow harmonic correlations NSCv(m,n) could be built from the nonlinear hydrodynamic response of higher flow harmonics. Besides, the past hydrodynamic calculations have shown that the contributions of nonlinear response can explain the symmetry plane correlations and its centrality dependence [37, 193]. Recently, the proposed nonlinear hydrodynamic coefficient, χmn [193], has been systematically studied and measured [37, 194, 195], which could be used to further constrain the initial conditions and η/s and to provide a better understand of the correlations between different flow harmonics.

6.1 p–Pb collisions at $\sqrt{s_{\text{N}N}}$=5.02 TeV

High energy proton-lead (p-Pb) collisions at the LHC were originally aimed to study the cold nuclear matter effects and provide the corresponding reference data for Pb–Pb collisions. However, lots of unexpected collective phenomena have been observed in the experiments. For example, the measured two particle correlations showed a symmetric double ridge structure on both the near-and away-side in high multiplicity p–Pb collisions at $\sqrt{s_{\text{N}N}}=$ 5.02 TeV [41-44]. Besides, negative 4- and 8-particle cumulants and positive 6-particle cumulants have been observed in the high multiplicity events [43-45]. In particular, all the multi-particle cumulants (including 4-, 6-, and 8-particles cumulants) are compatible to the ones obtained from the all-particle correlations with Lee-Yang Zero’s method, which corresponds to v₂{4} ≈ v₂{6} ≈ v₂{8} ≈ v₂{LYZ} [44]), as shown in Fig. 12 (This observation has also been confirmed by the later ATLAS [43] and ALICE Collaborations [45] measurements). Meanwhile, the obtained v₂ from two or four-particle cumulants are comparable to the ones from Pb–Pb collisions at 2.76 TeV [43, 44, 46, 196]. Recently, the ALICE collaboration has extended the investigation of anisotropic collectivity via azimuthal correlations of identified hadrons [46, 47]. A typical mass-ordering feature among the v₂ of pions, kaons and protons is observed in high multiplicity p-Pb collisions [46]. Similarly, the CMS Collaboration found a v₂ mass-ordering between ${\text{K}}_S^0$and $\Lambda (\overline{\Lambda })$[47].

Fig. 12.Full-screen View

(Color online) Multiplicity dependence of v₂, obtained from the Fourier decomposition of 2-particle azimuthal correlations, from multi-particle cumulants, and via the LYZ method, in Pb–Pb collisions at$\sqrt{s_{\text{N}N}}$=2.76 TeV (left) and p–Pb collisions at$\sqrt{s_{\text{N}N}}$=5.02 TeV (right) [44].

There are many theoretical efforts attempting to provide explanation for the flow-like behavior of the p–Pb collisions. In general, they can be divided into two big categories that don’t involve the final-state evolution of the medium but only account for initial-state effects [197-205], and that includes the final-state interactions, such as the hydrodynamics or kinetic model description [48-54, 206-211]. In this section, we will focus on reviewing the hydrodynamic calculations as well as the kinetic model investigations on the flow-like signals in the small p–Pb systems.

6.1.1 Results from hydrodynamic simulations:

Hydrodynamics is a useful tool to simulate the collective expansion of the created systems and quantitatively study and predict the final flow observable. Recently, the holographic duality calculations have shown that the size of the produced droplet is ∼1/Teff[212, 213], which indicates that hydrodynamics is possibly applicable for the small systems created in the high energy p–Pb and p–p collisions. Using 3+1-d hydrodynamic or hybrid model simulations, different groups have systematically studied the the multiplicities, mean p_T, final state correlations, and related flow data in p–Pb collisions at$\sqrt{s_{\text{N}N}}$=5.02 TeV [48-54]. In general, these hydrodynamic calculations could semi-quantitatively describe these different soft hadron data, which support the observation of collective flow in experiments of high energy p–Pb collisions.

Figure 13 (left) presents the hydrodynamic calculations for flow coefficients v₂ and v₃ of all charged hadrons in high multiplicity p–Pb collisions, which gives a rough fit of the data from the CMS collaborations [50]. It was also found that such fluid evolution also developed in the radial flow, which leads to a flatter transverse momentum spectra for various hadron species. As shown in Ref. [50], the average transverse momentum of the identified hadrons in p–Pb collisions can be consistently fitted by the hydrodynamic simulations. In contrast, the HIJING model, without any collective expansion, fails to describe the data. In the hydrodynamic language, the interaction between radial and elliptic flow re-distribute, the total momentum anisotropy to various hadron species, leading to a mass ordering of the flow harmonics. Figure 13 (right) shows that the hydrodynamic simulations roughly reproduce the v₂ mass-ordering of pions, kaons, and protons. Note that, other hydrodynamic calculations with different initial conditions and transport coefficients also obtained similar results. For details, please refer to Ref. [51-54].

Fig. 13.Full-screen View

(Color online) The hydrodynamic calculations of the elliptic and triangular flow coefficient of all charged particles (left panel) and elliptic flow of identified hadrons (right panel) in p-Pb collisions at $\sqrt{s_{\text{N}N}}$=5.02 TeV [50], together with a comparison with the CMS [196] and ALICE data [46].

Ref. [214] has shown that, in order to reproduce the multiplicity distribution of p–Pb collisions using the hydrodynamic calculations with Glauber initial conditions, the implementation of additional negative binomial fluctuations is necessary. Correspondingly, initial eccentricities are also modified which leads to a simultaneous fit of the v₂{2} and v₂4 data. In contrast, the early IP-glasma initial condition generates the initial energy distributions with an imprinted spherical shape of protons, which yields a very small v₂ for the p–Pb collision systems [54]. This motivates the recent investigations of the proton structure within the saturation framework, which indicates that the shape of the protons also fluctuates event-by-event [215, 216].

Note that the flow-like signals have also been observed in d–Au and ³He–Au collisions at RHIC. Compared to the p–A collisions at the LHC, the d–Au and ³He–Au collisions provide controlled initial geometry deformations, which are less sensitive to the details of initial state models and are helpful to check the hydrodynamic caculations. Recently, the STAR and PHENIX collaboration has measured the elliptic flow v₂ in d–Au collisions at $\sqrt{s_{\text{N}N}}$=200 GeV and the elliptic and triangular flow v₂ and v₃ in ³He–Au collisions at$\sqrt{s_{\text{N}N}}$=200 GeV [217-220]. The hydrodynamic calculations from different groups, using various initial conditions and the QGP shear viscosity, roughly described these extracted flow data. It was also found that v₂ and v₃ follows ε₂ and ε₃ from the initial state, which gives a support for the collective expansion in these small systems created at RHIC [51, 52, 221-223].

Compared with the Pb–Pb collisions, the initial sizes of the created systems in p–Pb collisions are much smaller. The subsequent collective expansion is expected to enlarge the size of the fireball, where the corresponding radii at the freeze-out can be measured by the Hanbury-Brown Twiss (HBT) correlations. In Ref. [224], the ALICE collaboration has measured the three-dimensional pion femtoscopic radii in p–Pb collisions at $\sqrt{s_{\text{N}N}}$=5.02 TeV, which showed that the size of the p–Pb systems is in between the ones obtained from the p–p collisions and peripheral Pb–Pb collisions. In general, the hydrodynamic calculations could roughly describe the HBT measurements, while the quantitative values from different model calculations are sensitive to the initial conditions and the imprinted initial sizes of the created fireball [223, 225, 226].

In Ref. [227], the validity of hydrodynamics for large Pb–Pb and small p–Pb systems at the LHC has been evaluated through tracing the space time evolution of the Knudsen number. It was found for Pb–Pb collisions, hydrodynamic simulations with η/s∼1/4π are always within the validity regime with the Knudsen numbers well below one. However, the related simulations for smaller p–A systems shows that the hydrodynamic descriptions have broken down at the T_dec = 100 MeV freeze-out boundary, even when using a minimum QGP shear viscosity as a input. Although such investigations will not preclude the collective flow and final state interactions, it is worthwhile to explore the physics of the small p–Pb systems within other frameworks beyond hydrodynamics.

6.1.2 Results from other approaches:

Without the final state interactions, the long range rapidity correlations in high energy p–p and p–Pb collisions have been calculated with the framework of Color Glass Condensate (CGC), which shows a good agreement with the di-hadron data from the CMS, ATLAS, and ALICE [197-200]. However the odd harmonics data disfavor these early CGC calculations without the rescattering contributions [201]. Without a proper hadronization procedure, such calculations can also not predict the flow data of the identified hadrons. Recently, it was proposed that the presence of the colored domains inside the proton and the nucleus breaks rotational invariance, which helps to generate elliptic and triangular flow during the scattrings between a dilute projectile of valence quarks and the nucleus [202-205]. An alternative approach is the classical Yang-Mills simulations, which treat both the proton and nucleus as dense QCD objects with high gluon occupancy and are more appropriate to describe the early time evolution of the created p–Pb systems in the high multiplicity events. Within such framework, Schenke and his collaborators have calculated the single and double inclusive gluon distributions and extracted the associated p_T dependent elliptic and triangular flow of gluons in high energy p–A collisions [228]. They found that the final state effects in the classical Yang-Mills evolution build up a non-zero triangular flow but only slightly modify the large elliptic flow of gluons created from the initial state [228]. Although this investigation only focused on the flow of anisotropy of gluons, the obtained large value of v₂ and v₃ indicates that such pre-equilibrium dynamics should be combined with the model calculations of the final state interactions, such as hydrodynamics or the Boltzmann simulations.

The flow signals in the p–Pb collisions have also been investigated within the framework of the multiphase transport model (AMPT) [206-210]. With a tuned cross-sections within the allowed range, σ∼1.5-3 mb, AMPT nicely fit the two particle correlations and the extracted v₂ and v₃ coefficients in high energy p–Pb collisions [206, 207]. Ref. [206, 210] has shown that AMPT generates a mass-ordering of v₂ and v₃ for various hadron species with the coalescence process tuning on. It was also surprisingly observed that the collective behavior in AMPT is built up by a small amount of interactions, where each parton undergoes two collisions on average. The escape mechanism proposed in Ref. [210, 229] seems to be responsible for the anisotropy buildup in AMPT but is dramatically different from the traditional flow development picture of hydrodynamics due to the strong interactions.

With an assumption that the high energy p–Pb collisions do not reach the threshold to create the QGP but only produce pure hadronic systems, Ref. [211] systematically investigated the 2 and 4 particle correlations of all charged and identified hadrons, using the hadron cascade model Ultra-relativistic Quantum Molecular Dynamics (UrQMD) [127, 128, 230]. Figure 14 shows the two and four -particle cumulants c₂{2} and c₂{4} of all charged hadrons, calculated from UrQMD and measured from ALICE. In general, c₂{2} decreases with the increase of the pseudorapidity gap, which is in agreement with the expectation of suppressing the non-flow effects with a large pseudorapidity gap. However, UrQMD still presents a strong centrality dependence of c₂{2} for |Δη|>1.0, which indicates that the remaining non-flow effects are still strong there. In Fig. 14 (right), the c₂{4} from the ALICE exhibits there is a transition from positive to negative values, which indicates the creation of flow-dominated systems for the high multiplicity events. In contrast, c₂4 from the UrQMD simulations keeps positive for all multiplicity classes, which illustrates that the p–Pb systems created by UrQMD are non-flow dominated.

Fig. 14.Full-screen View

(Color online) Centrality dependence of c₂{2} (left) and c₂4 (right) calculated from UrQMD [211] and measured by ALICE [45].

However, the generally believed collective expansion feature, the mass-ordering of v₂(pT), is reproduced in the UrQMD simulations. Fig. 15 shows that these high multiplicity events from UrQMD present a clear v₂ mass-ordering among pions, kaons, and protons, which are qualitatively in agreement with the corresponding ALICE measurement [46]. In UrQMD, the meson-baryon (M-B) cross sections from AQM are about 50% larger than the meson-meson (M-M) ones, which leads to the v₂ splitting between mesons and baryons in the UrQMD simulations. Fig. 15 also shows that, after switching off the M-B and M-M interaction channels, the characteristic feature of the v₂ mass-ordering disappears. Therefore, even without enough flow generation, the hadronic interactions still lead to a v₂ mass-ordering feature for a hadronic p–Pb system.

Fig. 15.Full-screen View

(Color online) v₂(pT) of pions, kaons, and protons in p–Pb collisions at $\sqrt{s_{{}_{\text{N}N}}}$ =5.02 TeV calculated from UrQMD with and without M-M and M-B collisions [211].

In Ref. [231], the created p–Pb systems are described by non-interacting free-streaming particles, followed by a harmonization procedure and a hadronic cascade evolution. Such non-hydrodynamic simulations showed that, although the elliptic flow is under-predicted, the triangular and quadrupolar flow are raised by the free-streaming evolution, which are comparable to the ones obtained from the hydrodynamic simulations. Meanwhile, the vn mass-orderings among pions, kaons, and protons have also been observed in such non-hydrodynamic p–Pb systems due to the hadronic interactions during the late evolution.

6.2 p–p collisions at $\sqrt{s_{\text{N}N}}$ =7 TeV and 13 TeV

Like the case for high energy p–Pb collisions, the long-range two-particle azimuthal correlations with a large pseudo-rapidity separation have also been observed in high-multiplicity p–p collisions at the LHC. This provides new insights for the novel dynamics of the small QCD systems [55-59]. For p–Pb collisions at $\sqrt{s_{\text{N}N}}$=5.02 TeV, the extensive measurements of the 2 particle and multi-particle correlations, extracted flow harmonics for all charged and identified hadrons, as well as the supportive hydrodynamic calculations strongly indicates that collective expansion has been developed in the small p–Pb systems. However, for high-energy p–p collisions at the LHC, the nature of the observed long-range correlation is still an open question (For different theoretical interpretations, please refer to Ref. [51, 54, 60, 197, 200, 232-237]).

Recently, the ATLAS Collaboration has measured the Fourier coefficients, vn, in p–p collisions at $\sqrt{s_{\text{N}N}}$=13 TeV, using the two-particle correlations as a function of the relative azimuthal-angle and pseudo-rapidity [57]. It was found that the extracted v₂ is approximately a constant as a function of multiplicity and its p_T dependence is very similar to the one measured in the p–Pb and Pb–Pb collisions [57]. The CMS collaboration further measured the vn coefficients for all charged hadrons, as well as for$K_S^0$and $\Lambda /\overline{\Lambda }$ in p–p collisions at $\sqrt{s_{\text{N}N}}=$5, 7, and 13 TeV, which observed a clear v₂ mass-ordering among all charged hadrons,$K_S^0$ and $\Lambda /\overline{\Lambda }$[59]. Furthermore, the CMS collaboration has measured the multi-particle cumulants, the key observable to probe the anisotropic collectivity. A negative sign of c₂{4} and a positive sign of c₂{6} appeared in the high multiplicity p–p collisions at$\sqrt{s_{\text{N}N}}=$13 TeV [59], which seems to indicate the development of anisotropic collectivity in high energy p–p collisions. However, the ATLAS Collaboration reported in Hard Probe 2016 conference that the multiplicity fluctuations could significantly bias the measurements of multi-particle cumulants [238], which indicates that non-flow might mimic the flow signal by pushing the c₂{4} to negative values. In order to avoid the bias from multiplicity fluctuations, the so-called “Method1”, which using the same multiplicity selection for the calculations of cumulants and N_trk, is applied. The obtained c₂{4}, which is less affected by multiplicity fluctuations, does not show negative sign for the multiplicity regions where negative values of c₂4 was reported by CMS.

For small systems, it is also very important to address and evaluate the non-flow effects. Generally, the multi-particle cumulants, e.g. c₂{4}, are able to suppress the non-flow of two-particle correlations in traditional Au+Au or Pb+Pb collisions. However, the non-flow contributions to the multi-particle correlations are still remained and might play an non-negligible role in the small p–p collision systems. Recently, the ALICE and ATLAS Collaborations have proposed new 4-particle cumulant methods with |Δη| gap separation, using 2- or 3-subevents [239, 240]. By selecting particles from different regions separated by a |Δη| gap, it is possible to further suppress the non-flow contributions in the multi-particle cumulants. This has been verified in the PYTHIA simulations [241]. The preliminary measurements in p–p collisions at 13 TeV, reported in QM2017 [239, 240], have shown that the non-flow effects are suppressed with these new 4-particle cumulant methods. A negative sign of the 4-particle cumulant was observed by ATLAS collaboration after implementing the 3-subevent method, while ALICE has not confirm the negative sign of c₂{4} with a |Δη| gap separation due to the limited statistics and relatively smaller acceptance.

Besides the multi-particle cumulants for single flow harmonics, the CMS Collaboration also measured the symmetric cumulants SC(m,n) and normalized symmetric cumulants NSC(m,n) in p–p, p–Pb, and Pb–Pb collisions [242]. It was found that the normalized NSC(3,2) are similar in p–Pb and Pb–Pb collisions, indicating that these two systems present similar initial state fluctuation patterns for the correlations between ε₂ and ε₃. While, the normalized NSC(4,2) shows certain orderings for the p–p, p–Pb, and Pb–Pb collision systems, which may associates with the different non-linear response and non-flow effects between the large and small systems.

In short, these recent measurements in p–p collisions $\sqrt{s_{\text{N}N}}$=13 TeV are aimed to evaluate whether or not collective flow has been created in high multiplicity p–p collisions. Future investigations, from both experimental and theoretical sides, are very crucial to further address this question and for a deep understanding of the underline physics in the small collision systems.