Flow correlations from a hydrodynamics model with dynamical freeze-out and initial conditions based on perturbative QCD and saturation

We extend the applicability of the hydrodynamics, perturbative QCD and saturation -based EKRT (Eskola-Kajantie-Ruuskanen-Tuominen) framework for ultrarelativistic heavy-ion collisions to peripheral collisions by introducing dynamical freeze-out conditions. As a new ingredient compared to the previous EKRT computations we also introduce a non-zero bulk viscosity. We compute various hadronic observables and flow correlations, including normalized symmetric cumulants, mixed harmonic cumulants and flow-transverse momentum correlations, and compare them against measurements from the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC). We demonstrate that the inclusion of the dynamical freeze-out and bulk viscosity allows a better description of the measured flow coefficients in peripheral collisions and enables the use of an extended centrality range when constraining the properties of QCD matter in the future.

We extend the applicability of the hydrodynamics, perturbative QCD and saturation -based EKRT (Eskola-Kajantie-Ruuskanen-Tuominen) framework for ultrarelativistic heavy-ion collisions to peripheral collisions by introducing dynamical freeze-out conditions. As a new ingredient compared to the previous EKRT computations we also introduce a non-zero bulk viscosity. We compute various hadronic observables and flow correlations, including normalized symmetric cumulants, mixed harmonic cumulants and flow-transverse momentum correlations, and compare them against measurements from the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider (LHC). We demonstrate that the inclusion of the dynamical freeze-out and bulk viscosity allows a better description of the measured flow coefficients in peripheral collisions and enables the use of an extended centrality range when constraining the properties of QCD matter in the future.

I. INTRODUCTION
Heavy-ion collisions at ultrarelativistic energies provide the means to produce and investigate experimentally quark-gluon plasma (QGP), a strongly interacting fluid of quarks and gluons. In recent years the two main collider experiments that have investigated QGP properties are the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL), and Large Hadron Collider (LHC) at CERN. In these experiments a small, short-lived, fluid-like behaving droplet of strongly interacting matter is created at nearly zero net-baryon density. The matter properties of QGP such as its equation of state (EoS) and transport coefficients are reflected in the detailed behavior of various experimental observables; see, e.g., Refs. [1][2][3][4][5][6][7].
The equation of state of strongly interacting matter at zero net-baryon density is currently well known from lattice-QCD computations, and the expected transition temperature T c ≈ 150 − 160 MeV [8][9][10][11] from hadronic matter to QGP is well within the reach of the LHC and RHIC experiments. Currently there are some experimental constrains on the equation of state [12][13][14], but even the lattice-QCD data allows some freedom in the EoS parametrizations [15]. The best knowledge about the transport properties of QCD matter is coming from the global fits of fluid dynamical computations to the available low-p T data from RHIC and LHC [15][16][17][18][19][20][21][22]. Currently, at least within the given models, the shear viscosity at temperatures near the QCD transition temperature is quite well constrained. However, the same cannot be said about the bulk viscosity. Even if the different analyses are based on very similar underlying models, the final constraints on the bulk viscosity can differ quite significantly depending on the details of the selected data and fine details of the models.
The experimental information about the collective dynamics and the spatial structure of the initial conditions is primarily encoded in the flow measurements.
The most basic quantities are the Fourier components of the azimuthal hadron spectra, usually called the flow coefficients v n . The measured flow coefficients reflect the collective fluid dynamical behavior of the system, as they are generated during the evolution of the system when the initial spatial inhomogeneities are converted into momentum-space anisotropies. In the fluid dynamical limit the driving force for this conversion is the inhomogeneus pressure gradients, and the effectiveness of the conversion is dictated by the EoS and the transport properties of QCD matter.
In the actual collisions the flow coefficients fluctuate strongly from event to event, and the fluctuations need to be explicitly considered when modeling the collisions. The presence of the flow fluctuations complicates the modeling, but at the same time they offer also a possibility to probe the initial conditions and the space-time evolution in much greater detail. For example, the relative fluctuation spectra of the elliptic flow coefficient v 2 are practically independent of the QCD matter properties, and reflect mainly the initial density fluctuations, giving thus a way to directly constrain the initial particle production [23] at least at the LHC energies, see the discussion in Ref. [24]. Moreover, the various observables measuring the correlations between the flow coefficients react to the matter properties and initial conditions in a nontrivial way, and offer further constraints on both of them. In particular, the correlations cannot be trivially reproduced just by reproducing the flow coefficients themselves [25].
The aim of this paper is to calculate various measurable flow-correlators by using relativistic second-order fluid dynamics with QCD-based initial conditions. The main ingredients that go into the computation are the matter properties, equation of state and transport coefficients, initial conditions for the fluid dynamical evolution given by the primary production of particles, and finally the conditions when the fluid dynamical evolution ceases and the fluid decouples into free hadrons.

arXiv:2206.15207v3 [hep-ph] 30 Jun 2023
The initial conditions are computed by using the perturbative QCD based Eskola-Kajantie-Ruuskanen-Tuominen (EKRT) saturation model [25,26], where the primary quantity is the minijet transverse energy computed in next-to-leading order perturbative QCD. The low-p T production of the particles is then controlled by a saturation conjecture, detailed in Sec. II. The EKRT saturation model is the main feature that gives a predictive power to our computation. Once the framework is fixed at some collision system, e.g. in central Pb+Pb collisions at the LHC, the collision energy, centrality, and nuclear mass number dependence of hadronic observables are predictions of the model [25,[27][28][29] Once the initial conditions are given, the remaining inputs to the fluid dynamical computation are the matter properties. The EoS is provided by the s95p parametrization of lattice-QCD results [30], and the specific shear viscosity η/s, is parametrized such that it has a minimum around the QCD transition temperature. As a new ingredient compared to the previous EKRT computations we introduce nonzero bulk viscosity, parametrized such that it is peaked close to T c . The main impact of bulk viscosity is to reduce the average p T of hadrons [6]. This allows us to relax our earlier [25,[27][28][29] rather high chemical freeze-out temperature T chem = 175 MeV, in order to better reproduce the measured identified hadron abundances, while still reproducing the measured average transverse momentum of hadrons.
Another new feature in the computation is the dynamical condition to decouple the system into free hadrons. The earlier EKRT results were computed using a constant-temperature decoupling at T dec = 100 MeV. It can be argued that the system decouples when the mean free path of hadrons is larger than the size of the system. The mean free path is a function of temperature, and if the system size is fixed the condition gives a constant temperature. However, the system size actually changes as function of time when the system expands, and moreover the system size varies from collision to collision: Central nuclear collisions produce a much larger system than peripheral ones. In order to account for the differences in the size of the systems, we introduce two conditions for decoupling. The global condition compares the overall size of the system to the mean free path, or here rather to the relaxation time in the second-order fluid dynamics, and the local condition that requires that the Knudsen number Kn, the ratio of microscopic and macroscopic length or time scales, is sufficiently small the fluid dynamics to be applicable [31,32]. We note that this approach, in particular the global condition, is slightly different from the earlier works where dynamical decoupling was developed [33,34].
The main advantage of using dynamical decoupling, besides that it is physically better motivated than the constant-temperature decoupling, is that it allows one to extend the agreement between the fluid computation and the measured flow coefficients towards peripheral nuclear collisions. In particular, the success of fluid dynamics in reproducing the flow coefficients in high-multiplicity proton-nucleus collisions [35][36][37][38][39][40][41][42] suggests that fluid dynamical models should then also describe peripheral nuclear collisions with similar hadron multiplicities. This paper is organized in the following way: In Sec. II we shortly review the EKRT saturation model. In Sec. III we introduce the second-order fluid dynamics, and give the parametrizations of shear and bulk viscosities, and the corresponding corrections to the hadron momentum distributions. In Sec. IV we detail the dynamical freezeout conditions, and in Sec. V we introduce the definitions of the experimental observables. The results from the computations are given in Sec. VI, where we show the new results with bulk viscosity and dynamical decoupling and compare those to the earlier predictions of the EKRT model. Finally the summary and conclusions are given in Sec. VII.

II. INITIAL CONDITIONS
The initial energy density profile is computed by using the EKRT saturation model [25,26,43,44]. It is based on the next-to-leading-order perturbative QCD (pQCD) computation of transverse energy (E T ) production, controlled by the low-p T cutoff scale p 0 determined from the local saturation condition [44], where ∆y is the rapidity interval, b is the impact parameter, K sat quantifies the uncertainty in the onset of saturation, and β quantifies the freedom in the NLO E T definition with low-p T cutoff. The solution p 0 = p sat of the saturation condition then inherits the √ s NN and A dependence from the NLO pQCD computation of E T , and the nuclear geometry enters through the product T A T A of the nuclear thickness functions, The local energy density at the formation time τ s = 1/p sat can then be written using p sat as At each point in the transverse plane the energy density is further evolved into a common initialization time τ 0 = 1/p sat,min ≈ 0.2 fm by using (0 + 1)-dimensional Bjorken expansion, where the minimum saturation scale p sat,min = 1 GeV. Below this scale the computed energy density profile is connected smoothly to the e ∝ T A T A profile. As in the earlier works, we take β = 0.8, and K sat is fixed from the charged particle multiplicity measured in central √ s NN = 2.76 TeV Pb+Pb collisions. For further details and explicit parametrizations of p sat , see Refs. [25,27,28].
The nuclear thickness functions are computed by first randomly sampling the nucleon positions from the Woods-Saxon nucleon density profiles. The Au and Pb nuclei are taken as spherical with a radius R = 6.38(6.7) fm for Au (Pb), and a thickness parameter d = 0.55 fm. As in Ref. [45], in the case of Xe we take into account the deformation by introducing the parameters β 2 = 0.162 and β 4 = −0.003 [46]. The Xe radius is R = 5.49 fm and the thickness parameter d = 0.54 fm.
The nuclear thickness functions are then computed by summing up the individual nucleon thickness functions, where T n is a Gaussian with a width σ = 0.43 fm. The event-by-event fluctuations emerge from the random positions of the nuclei, and impact parameter: The fluctuating T A T A profile leads to a fluctuating energy density profile through the T A T A -dependence of the saturation scale in Eq. (3). A randomly sampled collision event, i.e. the nucleon positions in the nuclei and the impact parameter between the two nuclei, is accepted using a geometric criterion: We require that there is at least one pair of colliding nucleons with a transverse distance less than σ NN /π, where σ NN is the inelastic nucleon-nucleon cross section. Here we take σ NN = 42 mb in √ s NN = 200 GeV Au+Au, mb in √ s NN = 5.023 TeV Pb+Pb, and σ NN = 72 mb in √ s NN = 5.44 TeV Xe+Xe collisions. We emphasize that this criterion is only used as a condition that nuclear collision happens at all, it is not needed in the computation of the initial profile.

III. FLUID DYNAMICAL EVOLUTION AND PARTICLE SPECTRA
After the hot strongly interacting system is produced at τ 0 ∼ 1/p sat , the subsequent spacetime evolution is computed using relativistic dissipative fluid dynamics. The basic equations of fluid dynamics are the local conservation laws of energy, momentum and conserved charges like net-baryon number. These can be expressed in terms of the energy-momentum tensor and charge 4currents as ∂ µ T µν = 0, and ∂ µ N µ i = 0. In what follows we shall neglect the conserved charges so that it is sufficient to consider only the energy-momentum tensor. It can be decomposed with respect to the fluid 4-velocity u µ as where the fluid velocity is defined in the Landau picture, i.e. as a time-like, normalized eigenvector of the energy momentum tensor, T µ ν u ν = eu µ . Here e = T µν u µ u ν is the local energy density, P = − 1 3 ∆ µν T µν is the isotropic pressure, and π µν = T µν is the shear-stress tensor.
The angular brackets denote the projection operator that takes the symmetric and traceless part of the tensor that is orthogonal to the fluid velocity, i.e., A µ = ∆ µν A ν and where ∆ µν = g µν − u µ u ν , and g µν is the metric tensor for which we use the g µν = diag(+, −, −, −) convention. The bulk viscous pressure is defined as Π = P −P 0 , where P is the total isotropic pressure and P 0 is the equilibrium pressure.
The conservation laws are exact, but they do not give sufficient constraints to solve the evolution. The simplest fluid dynamical theory follows by neglecting the dissipative effects completely. In that case the system is always in a strict thermal equilibrium, entropy is conserved, and the equation of state in the form P 0 = P 0 (e) closes the system. The dissipation plays, however, a significant role in the evolution of the system in heavy-ion collisions, and it cannot be readily neglected. The dissipative effects are contained in the shear-stress tensor and in the bulk viscous pressure. Therefore the remaining task is to write evolution equations for them. In the formalism of Israel and Stewart [47] the equations take the form τ π d dτ π µν + π µν = 2ησ µν + 2τ π π µ α ω ν α − δ ππ π µν θ − τ ππ π µ α σ ν α (8) + ϕ 7 π µ α π ν α + λ πΠ Πσ µν , where σ µν = ∇ µ u ν is the strain-rate tensor, ω µν = 1 2 (∇ µ u ν − ∇ ν u µ ) is the vorticity tensor, and θ = ∇ µ u µ is the expansion rate. The shear and bulk relaxation times are denoted by τ π and τ Π respectively, while first-order transport coefficients are the shear viscosity η and the bulk viscosity ζ. The coefficients of the nonlinear terms δ ΠΠ , λ Ππ , δ ππ , τ ππ , ϕ 7 , λ πΠ are second-order transport coefficients. Formally these equations can be derived from kinetic theory [47][48][49][50][51][52][53][54], by expanding around equilibrium and keeping terms up to the first order in gradients (or Knudsen number, a ratio of microscopic and macroscopic time/length scales, such as Kn ∼ τ π ∇ µ u µ , [55]), second order in inverse Reynolds number ∼ π µν /P 0 , and product of Knudsen number and inverse Reynolds number.
In this work the fluid dynamical setup is the same as in our previous works [4,5,25,27,28], i.e. we assume boostinvariant longitudinal expansion, so that it is enough to solve the equations of motion numerically in (2+1) dimensions [56]. The second-order transport coefficients in the Israel-Stewart equations are taken from the 14moment approximation to massless gas [48,49,51] and  bulk-related coefficients are from Ref. [57], i.e., where c 2 s is the speed of sound. The shear and bulk relaxation times are given by The remaining input to the equations of motion are the equation of state and the temperature dependence of the shear and bulk viscosities. The parametrizations of the shear viscosity to entropy density ratio are shown in Fig. 1, where η/s = 0.20 and η/s = param1 are the same as implemented in earlier works [25,27,28]. The new parametrization η/s = dyn has a similar linear QGP part as the previous parametrizations while the hadronic part follows a power law, with power P H , reaching its minimum (η/s) min at temperature T H followed by a constant part with width W min , i.e.
where S H and S Q are the slope parameters below T H and above T Q = T H + W min , respectively. The bulk viscosity is included together with the new η/s = dyn parametrization and its ratio to entropy density is plotted as a function of temperature in Fig. 2. Formally our paramaterization is written in a form: where (ζ/s) max , T ζ/s max , (ζ/s) width and a ζ/s are free parameters. The asymmetry parameter a ζ/s describes the asymmetry of the bulk viscosity peak in such a way that a ζ/s = 0 gives a completely symmetric peak. For the EoS we use the s95p parametrization [30] of the lattice QCD results that includes the chemical freeze-out, implemented as effective chemical potentials in the hadronic part of the EoS [58][59][60]. The earlier η/s = 0.20 and η/s = param1 parametrizations use chemical freeze-out temperature T chem = 175 MeV while the η/s = dyn parametrization uses T chem = 155 MeV.
The transverse momentum spectra of hadrons are obtained by computing the Cooper-Frye freeze-out integrals on the kinetic decoupling surface for the hadrons included in the hadronic part of the EoS. The 2-and 3-body decays of unstable hadrons are accounted for. For the earlier parametrizations η/s = 0.20 and η/s = param1 the kinetic decoupling surface is set to a constant T dec = 100 MeV temperature hypersurface while the η/s = dyn parametrization uses dynamical criteria (see Sec. IV for details) to determine the decoupling surface. The Cornelius algorithm [61] is employed to find the decoupling surface. The viscous correction δf i to each single-particle equilibrium momentum distribution, needed in the Cooper-Frye integrals, is implemented as in Refs. [2,[62][63][64], where k µ is the four-momentum of a given hadron, E k = u µ k µ is the energy of the hadron in the local rest frame, Dynamical freeze-out: (15) Here g i is the degeneracy factor of a given hadron species i, and the sum includes all the species in the EoS.
The fluid dynamical evolution and the transverse momentum spectra are computed for each collision event.
The events are then grouped to the centrality classes according to the final charged particle multiplicities. However, if the experiments report the centrality of the collision by using the number of wounded nucleons, we can compute it by using the geometric collision criterion detailed at the end of Sec. II.
Numerical values of the parameters used here for the η/s = dyn parametrization are shown in Table I. The initial state parameter K sat is tuned to produce the same charged particle multiplicity in 2.76 TeV Pb+Pb collisions as obtained in the ALICE measurements. Parameters of the shear viscosity and the dynamical freeze-out are iteratively adjusted to obtain results that match with ALICE measurements of v n {2} in 2.76 TeV Pb+Pb collisions. Further tuning of the hadronic part of the η/s parametrization is done to also match STAR measurements of v n {2} in central to mid-central 200 GeV Au+Au collisions. The chemical freeze-out temperature is adjusted together with the parameters of bulk viscosity to achieve a good simultaneous agreement of the pion average p T and the proton multiplicity.
We note that the idea here is that bulk viscosity in hadronic evolution is mainly described by chemical freeze-out [65][66][67]. In chemical freeze-out the corresponding bulk relaxation time is formally infinite, or at least much longer than the evolution time of the system, and the dynamics of the bulk pressure related to the non-equilibrium chemistry in this case cannot be readily computed using Israel-Stewart type of theory that assumes that the relaxation times are smaller than the evolution timescale. Instead, the bulk viscosity that is parametrized here should be thought as the residual bulk viscosity that is not included in the partial chemical freeze-out formalism [60]. In practice, the condition that low-temperature bulk viscosity is described mainly by chemical freeze-out is set by adjusting the asymmetry parameter a ζ/s in the parametrization such that bulk viscosity over entropy density becomes very small near and below the chemical freeze-out temperature. We want to emphasize here that this is only one example parametrization which seems to give a good agreement with the LHC and RHIC measurements. To get more detailed estimates of the parameters and their errors and correlations, a global analysis of heavy-ion observables and the parameter space is needed.

IV. DYNAMICAL FREEZE-OUT
When modeling heavy-ion collisions using hydrodynamics the kinetic freeze-out is usually set to take place at a constant-temperature hypersurface. The basic argument is that the fluid decouples into free particles when the temperature dependent mean free path of the particles becomes of the same order as the size of the system R, i.e. λ mfp (T ) ∼ R. If the system size was a constant, this condition would give a constant freeze-out temperature. However, in reality the system size changes as a function of time, and moreover it can differ significantly from collision to collision. In particular, the systems created in central collisions are much larger than the ones created in peripheral collisions.
A typical way to solve this issue is to connect fluid dynamics to a microscopic hadronic afterburner that automatically takes care of the freeze-out. However, a drawback in this approach is that it can easily lead to unphysical discontinuities in the transport coefficients, as at typical temperatures at the switching between fluid dynamics and hadron cascade the η/s values in the fluid evolution are O(0.1), whereas on the hadron cascade side they are O(1) [68][69][70]. Instead of a coupling to hadron cascade, in this work we treat the whole evolution, including the hadronic phase, using fluid dynamics. This has the specific advantage that it allows us to keep all the transport coefficients continuous throughout the whole temperature range realized in the evolution.
In order to account for the nontrivial system size dependence of the freeze-out, we determine the decoupling surface dynamically [33,34] using two different conditions. The applicability of fluid dynamics requires that the local Knudsen number is sufficiently small, and fluid evolution becomes effectively free streaming when Kn 1. In comparisons between kinetic theory and fluid dynamics it was shown that a constant Knudsen number freeze-out in fluid dynamics catches very well the freezeout dynamics of the kinetic evolution [32]. On the other hand, even if the local condition gives that fluid dynamics is applicable, the overall size of the system can still be small compared to the mean free path of the particles. In order to account for this kind of nonlocal freeze-out, we impose a second condition that the fluid element decouples when the mean free path is of the same order as the system size. Hence, our dynamical freeze-out setup is determined by the following two conditions: where C Kn , C R ∼ 1 are some proportionality constants and R is the size of the system. Here we have assumed that the mean free path is proportional to the relaxation time. The additional gamma factor in the second equation takes into account that the size of the system is calculated in the center-of-momentum frame of the nuclear collision, while the relaxation time is calculated in the fluid rest frame. To make sure that we are not in the QGP phase when freeze-out happens we also require that at the freeze-out surface T < 150 MeV. In order to use latter condition (17) we need to have some kind of estimate for the system size which, however, is not uniquely determined. In this work we define the size of the system as where A is the area in the x, y plane in which Kn < C Kn . Additionally we take into account the possibility that the system may consist of multiple separate areas of a fluid and calculate the system size for each of these regions separately. We note that our approximation of the system size is close to the maximum length that a particle must travel from the center to the edge of the system. In practice, however, most of the matter is distributed closer to the edges of the system and most of the particles are moving with the fluid also towards the edge. For this reason the actual size of the system that the particles see can be significantly smaller than R, and as a result the proportionality constant C R can also be significantly smaller than 1.
In summary, here we have on the one hand reduced a possibly complicated non-equilibrium dynamics of the hadronic evolution in the dynamical treatment of kinetic freeze-out, and on the other hand we treat the nontrivial chemistry in the hadronic evolution as a constanttemperature chemical freeze-out. While such an approach may not catch the full microscopic details of the freeze-out dynamics, the purpose is that it would still capture its essential features. A clear advantage is, as mentioned above, that it allows us to keep the transport coefficients of the matter continuous throughout the evolution, and at the same time it also allows us to get constraints for the hadronic part of the transport coefficients. As we can see, the physical picture of the evolution is somewhat different from the typical hybrid hy-dro+cascade models, where the low viscosity QGP evolution is immediately followed by high-viscosity hadronic evolution. In our picture the peripheral collisions decouple practically immediately after the hadronization, but in the central collisions there can still be quite long low-viscosity evolution in the hadronic phase. This is demonstrated in Fig. 3 where the entropy-flux-weighted average freeze-out temperature is plotted as a function of centrality for the η/s = dyn parametrization introduced in Sec. III. We can also notice that the average freeze-out temperature is sensitive to the collision energy and size of the colliding nuclei.

V. FLOW COEFFICIENTS AND CORRELATORS
The fluid dynamical computation gives a singleparticle transverse momentum spectrum of hadrons for each event, and its azimuthal modulation can be expressed by its p T dependent Fourier components v n (p T ) and the phases or event-plane angles Ψ n (p T ), The flow coefficients can be expressed in a convenient way by a complex flow vector V n as where the angular brackets denote an average Similarly, the p T -integrated flow coefficients can be defined as where the average is defined as and the p T -integrated multiplicity dN dy is defined with the same p T integration limits p T,min and p T,max as above. In addition it is possible to use a p T or an energy dependent weight w in the p T integration.
In the following we will write down the expressions of various measurable p T -integrated quantities, but suppress the rapidity, weight, and p T integration limits from the notation. The p T limits will be denoted explicitly when we show our results. Unless otherwise stated, the weight function w = 1.
In the fluid dynamical simulations of heavy-ion collisions we are working directly with continuous particle distributions. In the experiments this is not the case, but each event is measured as a finite number of particles. Therefore, the definitions above are not directly applicable, but the flow coefficients are rather defined through particle correlations. As an example of a two-particle correlation and its continuum limit we can write where N e is the number of hadrons in the event, and dN 2 /dφ 1 dφ 2 is a two-particle distribution function that can be written as a sum of the product of the singleparticle distribution functions and a direct correlation where the direct part emerges, e.g. due to hadron decays. It is a genuine two-particle correlation that is absent if all the correlations between the hadrons are due to the underlying collective flow. If the direct component can be neglected, the two-particle correlation above can be written in the continuum limit as In this limit the two-particle correlator can be written in terms of the flow coefficient. This particular correlator is referred to as the two-particle cumulant, and its average over events gives the two-particle cumulant v n {2}, where · · · ev denotes the average over the events. A similar reasoning leads to a multitude of flow observables. Here we write down only the continuum limit in the absence of direct or nonflow correlations. It should be noted, however, that although the experimental procedures try to suppress the nonflow parts by e.g. requiring a rapidity gap between each pair of hadrons, it is still possible that some of the observables are still plagued by the non-flow. With the current setup we cannot address the non-flow part theoretically, but will assume that the experimental techniques remove them completely.
In a naive picture one may think that the flow coefficients are generated independently as a fluid dynamical response to the corresponding eccentricities of the initial conditions, v n ∝ ε n . In practice, however, this picture holds only for the elliptic flow coefficient v 2 and to a lesser degree for v 3 [23,71], and even then the relation between v 2 and ε 2 ceases to be linear when ε 2 becomes large in noncentral collisions [25]. In general, the flow coefficients are not independent of each other, but both the correlations between the eccentricities in the initial conditions and the nonlinear fluid dynamical evolution generate correlations between them. The degree of the correlation can be measured through various observables that correlate both the magnitudes of the flow, v n , and the event-plane angles Ψ n [71].
A measurable way to quantify the degree of correlation between the flow coefficients is the so called symmetric cumulant [72], defined as where it is important to notice that the event-average is performed with powers of multiplicity as a weight, as denoted in the above equation. An advantage of this definition is that at the particle correlation level the latter term in the definition removes the direct two-particle correlations from the first term, which in turn is a fourparticle correlator at the particle level. Thus the direct two-particle nonflow does not affect the symmetric cumulant. The symmetric cumulant is not a correlator in a sense that it depends not only on the degree of correlation between v n and v m , but also on their absolute magnitudes. On the other hand, the normalized symmetric cumulant, defined as is a measure of only the correlation. The downside of the normalized version is that the normalization can be affected by the direct two-particle nonflow contributions.
The symmetric cumulants measure only correlations involving two second-order flow coefficients. The more general mixed harmonic cumulants (M HC) were introduced in Ref. [73] to give observables that can quantify the correlations between between more than two flow coefficients with higher-order moments of v n 's. Like symmetric cumulants, mixed harmonic cumulants are also constructed in such a way that lower order correlations are removed from multiparticle correlations and the definition of M HC containing two second order flow coefficients is identical to the symmetric cumulants, i.e.
Mixed harmonic cumulants for six-particle correlations involving moments of v 2 and v 3 can be defined as where · · · i = · · · ev,N i . Similarly one can define mixed harmonic cumulants for eight-particle correlations be-tween v 2 and v 3 as and for six-particle correlations between v 2 , v 3 and v 4 as Analogously to normalized symmetric cumulants one defines normalized mixed harmonic cumulants as A complementary observable to the symmetric cumulants, usually referred to as the event-plane correlator, is defined as [74] cos(k 1 Ψ 1 + · · · + nk n Ψ n ) where the k n 's are integers with the property n nk n = 0 so that the correlator is independent of the azimuthal orientation. Despite its name it actually measures a correlation between both the magnitudes of the flow and eventplane angle, and in this sense provides complementary information to the symmetric cumulants above.
These correlations as such provide information that is independent from the flow magnitudes themselves, and give further independent constraints to the initial conditions and transport coefficients. However, it is interesting that the event-plane correlations are closely related to the magnitude of nonlinear response to the initial conditions [75]. The basic idea in quantifying the nonlinear response is that the complex flow vector V n is divided into a linear part V nL that is assumed to correlate only with the corresponding initial state eccentricity ε n , and into a non-linear part that is independent of ε n [71]. If we consider the simplest possible nonlinear contributions, we can write where χ's are the nonlinear response coefficients. Note that the nonlinear parts include only the largest flow vectors V 2 and V 3 that can also, to a reasonable approximation as discussed above, assumed to have only the linear part V 2 = V 2L and V 3 = V 3L . If we further assume that the linear and nonlinear parts are uncorrelated, we may express the response coefficients as and the linear parts of V 4 and V 5 can be written as The connection between the event-plane correlators and the nonlinear response coefficients can be seen by observing, e.g. that so that the two measures differ by a normalization factor that depends on the magnitude of the flow, but not on correlators. A similar connection can also be made between the other χ's. A more complete list of relations can be found from Refs. [75,76]. Even though the nonlinear response coefficients and the correlations between the flow harmonics give information about the initial state eccentricities and their conversion to momentum space anisotropies, they do not directly probe the size of the initial nuclear overlap region which is more sensitive to the average p T fluctuations. Thus, the correlation between the flow coefficients and the average p T is a good probe of the initial state structure [77]. This flow-transverse-momentum correlation is defined by a modified Pearson correlation coefficient [78]  where the event-by-event variance at a fixed multiplicity for some observable O is defined bŷ

VI. RESULTS
In this section we present the results for hadron multiplicities, average p T , flow coefficients and correlations calculated from the EKRT pQCD + hydrodynamics framework with the bulk viscosity and the dynamical freezeout, and compare these against the results from our earlier works [25,[27][28][29] with the constant-temperature freeze-out and without the bulk viscosity. The systems we show here are 200 GeV Au+Au, 2.76 TeV Pb+Pb, 5.023 TeV Pb+Pb, and 5.44 TeV Xe+Xe collisions. As explained in Sec. III, the initial conditions, the transport coefficients and the freeze-out parameters are fixed on the basis of 200 GeV Au+Au and 2.76 TeV data from RHIC and LHC. For both Pb+Pb collision systems we run 40000 event simulations to get better statistics for the symmetric cumulants while for other collision systems we did 20000 event simulations. The statistical errors for different quantities are estimated, as in Ref. [76], via jackknife resampling.

A. Multiplicity, average pT and flow
In Fig. 4 we show the centrality dependence of charged hadron multiplicities for all the above systems compared to the STAR [82], PHENIX [83], and ALICE [79][80][81] data. The essential parameter that controls the multiplicity is K sat in the local saturation criterion. This coefficient is fixed from the multiplicity in 0-5 % 2.76 TeV Pb+Pb collisions. The centrality, √ s NN , and nuclear mass number dependence are predictions of the model. The value of K sat depends on the chosen η/s(T ) and ζ/s(T ) parametrizations due to the different entropy production with different shear and bulk viscosities. However, the final results for the multiplicities are in practice the same for all parametrizations and they agree excellently with the experimental data across all centrality classes and collision energies. The centrality dependences of identified particle multiplicities for 200 GeV Au+Au 2.76 TeV Pb+Pb and 5.023 TeV Pb+Pb collisions are shown in Fig. 5 (left). All of the parametizations manage to produce the same pion multiplicities as the ALICE and PHENIX measurements while the kaon multiplicities differ significantly from the experimental data. The ratio between the proton and pion multiplicities is mostly controlled by the chemical freeze-out temperature. Parametrizations η/s = 0.2 and η/s = param1 use T chem = 175 MeV in order to obtain the same average p T for pions in 2.76 TeV Pb+Pb collisions as the ALICE measurements. However this comes with the drawback that the proton multiplicities differ from the experimental data by a factor of ∼ 2. The addition of the bulk viscosity in the η/s = dyn parametrization enables the possibility to use T chem = 155 MeV which clearly improves the proton multiplicities. However, there is still some discrepancy left that is most visible in the most central collisions at the LHC.
In Fig. 5 (right) we show the average p T of identified particles as a function of centrality for 200 GeV Au+Au, 2.76 TeV Pb+Pb and 5.023 TeV Pb+Pb collisions. Compared to the earlier results, the η/s = dyn parametrization improves the agreement with the experimental data across both collision systems, except for kaons at the LHC energies. In particular, the relative change of the proton p T as a function of centrality is reproduced better. This improvement is not only due to the addition of the bulk viscosity but also the dynamical freeze-out plays a major part by affecting the lifetime of the fluid.
The centrality dependencies of the p T -integrated flow coefficients v 2 {2}, v 3 {2}, and v 4 {2} in all studied systems are shown in Fig. 6. The shear viscosity and the dynamical freeze-out parameters of the η/s = dyn parametrization were tuned to approximately reproduce v 2 {2} in 2.76 TeV Pb+Pb collisions while also reproducing v 2 {2} in central to mid-central 200 GeV Au+Au collisions. The most essential feature of the dynamical freeze-out is that the smaller collision systems freeze out earlier in the hadronic phase. This means that there is less time for the initial state eccentricities to convert to the momentum space anisotropies in peripheral collisions. Indeed, as seen in Fig. 6, all p T -integrated flow coefficients for the η/s = dyn parametrization are significantly smaller in peripheral collisions than the results of the η/s parametrizations from the earlier works that used a constant-temperature decoupling surface. As can be seen from the comparison to measurements, the η/s = dyn parametrization reproduces well the centrality dependence of all flow coefficients in all LHC collision systems and clearly improves the results from the earlier ones in peripheral collisions. The biggest discrepancy with the data and the model calculation is the 40 − 80% -centrality range in 200 GeV Au+Au collisions. In this region especially the predictions for the flow coefficients v 3 {2} and v 4 {2} are well outside of the error bars of the measurements. There are multiple possible reasons for this. First of all, due to the lower multiplicity in the 200 GeV Au+Au collisions it is reasonable to expect significantly larger nonflow effects compared to the LHC systems. Additionally, the δf corrections to the particle spectra are much larger at RHIC than at LHC which adds additional uncertainty to the RHIC results. Lastly, we do not include any nucleon substructure [91], initial flow or nonzero π µν to our initial state model and effects of these modifications are still under investigation. We note that other groups report very similar flow coefficients in peripheral RHIC collisions; see, e.g., Refs. [19,92] The change in the magnitude of the flow coefficients is quite modest from 2.76 to 5.023 TeV Pb+Pb collisions, and a better way to quantify the change is to plot the ratio of the coefficients between the two collision energies. The ratio is also a more robust prediction from fluid dynamics and less sensitive to fine tuning of η/s(T ), for a discussion see Ref. [93]. The predictions for the ratios of v n {2} in Pb+Pb collisions at 2.76 to 5.023 TeV are shown in the upper panel of Fig. 7. The predicted increase ranges from up to 8 % for v 2 to up to 25 % for v 4 . The predictions match well with the ALICE measurements for central to mid-central collisions, only in the most peripheral collisions the η/s = dyn parametrization overestimates the data slightly, especially in the case of v 4 , but there the experimental errors of the ratios are also quite large.
The situation is quite different in the case of Xe+Xe collisions. The ratio of the flow coefficients between the 5.44 TeV Xe+Xe and 5.023 TeV Pb+Pb collisions is shown in the lower panel of Fig. 7. The change in the flow coefficients is significantly larger than in the previous case, even if the collision energy is almost the same in Xe+Xe as in Pb+Pb collisions. The reason is that the system size is quite different when the nuclear mass number changes from A = 208 to A = 129. The most striking feature is the strong increase of v 2 in central Xe+Xe collisions compared to Pb+Pb collisions. A significant factor in the increase is the shape deformation of Xe nuclei. The deformation enhances the initial elliptic eccentricity fluctuations compared to the spherical double magic Pb nuclei. As a result the elliptic flow is 30 % higher in the Xe case. The fact that we correctly predict this increase by taking into account the nuclear deformation is further evidence that the azimuthal asymmetries in the p T spectra are resulting from a fluid dynamical response to the initial geometry. The event-plane correlations, defined in Eq. (35), quantify the correlation between the event-plane angles Ψ n , and also between the flow magnitudes v n . The com-puted event-plane correlations in 2.76 TeV Pb+Pb are shown in Fig. 8. Only a slight separation between the dynamical freeze-out and earlier η/s(T ) parametrizations can be seen and all parametrizations are able to describe the data. The most notable exceptions are the correlations involving the event-plane angle Ψ 6 , which are very sensitive to δf corrections. In these, the η/s = dyn parametrization slightly improves the agreement with the data from the earlier works. This is mostly due the fact that the η/s = dyn parametrization has lower shear viscosity and thus smaller δf corrections. The eventplane correlations have only been measured for 2.76 TeV Pb+Pb collisions which is why we do not show results for other collision systems. The symmetric cumulants, defined through Eq. (28), are complementary to the event-plane correlators in the sense that they depend on the correlation between the flow magnitudes v n like the event-plane correlators, but are independent of the event-plane angles. The symmetric cumulants themselves are not a measure of correlation, but depend explicitly on the magnitude of v n , and not only on the degree of correlation. The corresponding correlation measure is defined through the normalized symmetric cumulants, Eq. (29).
The normalized symmetric cumulants in 2.76 TeV Pb+Pb collisions are shown in Fig. 9 compared to the ALICE data [95]. As in the case of event-plane correlations, there are only small differences between the three η/s parametrizations. The overall agreement between the data and the computations is good, but with a notable exception that in peripheral collisions we underpredict the N SC(2, 4) correlation. The collision energy dependence of the normalized symmetric cumulants is weak, as can be seen in Fig. 10 where we show them in 5.023 TeV Pb+Pb collisions.
In Fig. 11 we show the normalized symmetric cumulants in 200 GeV Au+Au collisions. Note that here the centrality of the collisions is given by the number of participants, as reported by the STAR Collaboration [96]. Compared to Pb+Pb collisions we see much more separation between the dynamical freeze-out and earlier parametrizations for the N SC (3,4), N SC(3, 5) The correlations between higher order moments of two or three flow coefficients can be studied using the mixed harmonic cumulants which provide information that is independent of the normalized symmetric cumulants. The EKRT model predictions for nM HC(v 2 2 , v 2 3 , v 2 4 ) and nM HC(v k 2 , v l 3 ) are compared against the ALICE measurements for 5.023 TeV Pb+Pb collisions in Fig. 12. As can be seen there are only modest differences between the parametrizations and the statistical errors in our simulations are already quite large, especially with nM HC(v 4 2 , v 4 3 ). This is expected, since the correlations between v 2 and v 3 are thought to be more sensitive to the initial state rather than to the dynamics of the system. Our predictions seem to agree quite well with the data except for nM HC(v 4 2 , v 4 3 ), for which we predict a stronger correlation in peripheral collisions than what is measured. Fig. 13  In Fig. 14 we show the higher-order flow coefficients v 4 , v 5 , and v 6 compared to the ALICE data [99] in 2.76 TeV Pb+Pb collisions. As can be seen in the figure, the η/s = dyn parametrization seems to slightly underpredict the higher order flow coefficients in peripheral collisions, while the η/s = 0.2 parametrization manages to reproduce the data quite well. For v 6 we point out that the measured flow is larger in 2.76 TeV than in 5.023 TeV collisions, as can be seen by comparing measurements with Fig. 16, which is in conflict with the behavior of the other flow coefficients. We also note that the difference between the earlier parametrizations η/s = 0.2 and η/s = param1 is more visible here than in the case of lower-order flow coefficients. The corresponding nonlinear response coefficients are shown in Fig. 15. As explained in Sec. V they are closely related to the event-plane correlations, and the good agreement of the calculated response coefficients with the ALICE data is consistent with the good agreement between the calculated and the measured ATLAS eventplane correlations in Fig. 8.

Finally in
The same flow and response coefficients as above, but for 5.023 TeV Pb+Pb collisions, are shown in Figs. 16 and 17, respectively. Together with other higher order flow harmonics we also show v 7 , v 8 , and v 9 , which are only measured for the 5.023 TeV energy. Here we see that the parametrization that uses dynamical freeze-out predicts the higher order flow coefficients quite well while the parametrizations from earlier works are slightly above the measurements.
The response coefficients are not directly proportional to the magnitude of the flow coefficients, or the proportionality is partly canceled by the normalization. That is to say that the agreement in the response coefficients with the ALICE data is similar as at the lower collision energy even though we cannot exactly reproduce the higher order v n 's for both collision energies simultaneously.
The overall agreement with the higher-order flow coefficients with the data is quite similar for both the earlier and current EKRT setup. The improvements due to the dynamical decoupling are not as clear as for v 2 .
However, the differences between the parametrizations are also larger, highlighting the fact that higher-order coefficients, and their √ s NN dependence give important constraints to the determination of shear viscosity.

VII. SUMMARY AND CONCLUSIONS
We have presented the results for the low-p T observables in Pb+Pb, Au+Au, and Xe+Xe collisions at RHIC and LHC energies from the fluid dynamical computations using the NLO pQCD based EKRT model for the initial conditions. Compared to the previous EKRT works in Refs. [25,27,28] we have now added the bulk viscosity together with the dynamical decoupling conditions to improve the validity of our model in peripheral collisions.
The overall agreement of the computed results with the data is very good in particular for the √ s NN , A, and centrality dependence of the charged hadron multiplicity. This is mainly a feature of the EKRT initial conditions. The main uncertainty in the EKRT model is the K sat parameter in the saturation condition, but this can be essentially fixed from one measurement of charged hadron multiplicity. Even if the value of K sat depends on the η/s parametrization through the entropy production during the fluid dynamical evolution, the final results for the √ s NN , A, and centrality dependence are practically independent of the K sat value, making them very robust predictions of the EKRT model. The most significant effect of the dynamical freeze-out can be seen in the absolute magnitude of the flow coef- ficients v n . We have demonstrated that we can reproduce the experimental data for v 2 and v 3 across the centrality range 0 − 80 % in all the collision systems with the exception of peripheral RHIC collisions. This is a significant improvement from the constant-temperature freeze-out which only manages to describe the data up to the 30 − 40 % centrality class. The higher harmonics v 4 , v 5 , and v 6 are quite similarly described by both the earlier computations and the current setup, but the differences between the η/s parametrizations are also more pronounced. On the other hand, the relative increase of the flow coefficients from 2.76 TeV Pb+Pb to 5.023 TeV Pb+Pb and 5.44 TeV Xe+Xe collisions is well described in all the centrality classes shown here. The addition of the dynamical freeze-out together with the bulk viscosity has also made it possible to improve the simultaneous agreement of the identified particle multiplicities and the mean transverse momenta with the measurements. We have also shown the EKRT model predictions for the most recent correlation measurements. Our results for the symmetric cumulants, the mixed harmonic cumulants, the response coefficients, and closely related eventplane correlators are very similar to the earlier EKRT results and the agreement with the data remains reasonably good. The most notable differences are in N SC (2,4) correlators in peripheral collisions, where the predictions are visibly below the experimental data. The effect of the dynamical freeze-out and the bulk viscosity can be seen in the flow-transverse-momentum correlators ρ(v 2 n , [p T ]), where we demonstrated a better quantitative agreement with the experimental measurements in peripheral col-lisions than given by the previous EKRT computations. Especially, we obtained the correct sign in ρ(v 2 4 , [p T ]) correlation in peripheral collisions.
In conclusion, we have introduced dynamical freezeout conditions to model the decoupling of the fluid to free hadrons. In particular, the aim was to capture the essential features of the decoupling that take into account the system size variations at different collision energies and centralities. The clear benefit here is that it allows us to keep the transport coefficients continuous throughout the whole temperature range, without unphysical discontinuities that can appear at a switching between fluid dynamics and hadron cascade. At the same time it is then possible to use the measured data to constrain the QCD matter transport properties also in the hadronic phase.
We emphasize that in spite of the extensive iteration work done, the parametrizations shown here do not necessarily represent the absolute best fit to the data. For that we would need to do a full statistical global Bayesian analysis of the parameter space. This we have left as a future work. However, we have demonstrated that we can reproduce the measured LHC and RHIC low-p T observables reasonably well, and the dynamical decoupling leads to quite a different spacetime picture compared to many hydro+cascade models. Instead of a very viscous hadronic evolution directly after the low-viscosity QGP evolution, in the picture presented here the low-viscosity evolution can extend to quite low temperatures on the hadronic side. The data are from the ATLAS Collaboration [98].

ACKNOWLEDGMENTS
We thank Pasi Huovinen for discussions, and for providing us with the EoS tables. We acknowledge the financial support from the Jenny and Antti Wihuri Foundation, and the Academy of