Hydrodynamic fluctuations and ultra-central flow puzzle in heavy-ion collisions

One of the long-standing problems in the field of high-energy heavy-ion collisions is that the dynamical models based on viscous hydrodynamics fail to describe the experimental elliptic flow $v_2$ and the triangular flow $v_3$ simultaneously in ultra-central collisions. The problem, known as the"ultra-central flow puzzle", is specifically that hydrodynamics-based models predict the flow ratio of the two-particle cumulant method $v_2\{2\}/v_3\{2\}>1$ while $v_2\{2\}/v_3\{2\} \sim 1$ in the experimental data. In this Letter, we focus on the effects of hydrodynamic fluctuations during the space-time evolution of the QGP fluid on the flow observables in the ultra-central collisions. Using the (3+1)-dimensional integrated dynamical model which includes relativistic fluctuating hydrodynamics, we analyze the anisotropic flow coefficients $v_n\{2\}$ in 0-0.2% central Pb+Pb collisions at $\sqrt{s_\text{NN}}=2.76~\text{TeV}$. We find that the hydrodynamic fluctuations decrease the model overestimate of $v_2\{2\}/v_3\{2\}$ from the experimental data by about 19% within the present setup of $\eta/s = 1/2\pi$. This means that the hydrodynamic fluctuations qualitatively have an effect to improve the situation for the puzzle, but the effect of the hydrodynamic fluctuations alone is quantitatively insufficient to resolve the puzzle. The decrease of the ratio largely depends on the shear viscosity $\eta/s$, which calls for future comprehensive analyses with, for example, a realistic temperature-dependent viscosity.

reproduce the experimental behavior v 2 {2} ∼ v 3 {2}. This contradicts our expectation that models based on hydrodynamics perform better in central collisions because of larger multiplicities and volumes of locally thermalized domains. This implies that the state-of-the-art dynamical models, which are also intensively used in Bayesian analyses, potentially have missing pieces. So far, a number of attempts have been made to resolve the puzzle from different aspects, such as improved descriptions of initial conditions [26][27][28][29][30][31][32][33][34], effects of the transport coefficients [26,28,35,36], and the equations of state [37], but none has yet succeeded in explaining the experimental behavior v 2 {2} ∼ v 3 {2} within the dynamical models based on hydrodynamics.
In ultra-central collisions, an approximate linear mapping from the initial geometrical anisotropies ε n to the anisotropic flow coefficient v n is known [38][39][40], v n = κ n ε n . (1) Although initial anisotropies of the second and the third orders are almost the same ε 2 ∼ ε 3 in the initial stage, the viscosity leads to the ordering κ 2 > κ 3 , and thus v 2 > v 3 . This is because smaller structures, which are associated with higher anisotropies, are smeared out more by the viscosity [32,33]. However, ε n and v n can be de-correlated by additional eventby-event fluctuations in the hydrodynamic and hadronic stages. This happens more easily in ultra-central collisions due to the smaller geometrical origin of the initial anisotropies. The effect of the initial fluctuations has been intensively investigated in existing studies, but the effects of dynamical fluctuations arising in the hydrodynamic and later stages have not been comprehensively studied so far.
In this Letter, we focus on hydrodynamic fluctuations [41,42] arising in the hydrodynamic stage. Hydrodynamic fluctuations are thermal fluctuations related to the hydrodynamic description of the system. In hydrodynamics, microscopic degrees of freedom are integrated out by coarse-graining so that the system is described by only a few slow macroscopic variables. However, macroscopic dynamics cannot be completely separated from microscopic one. The microscopic dynamics induces the thermal fluctuations of the macroscopic variables on an event-by-event basis, which is nothing but the hydrodynamic fluctuations. Since the hydrodynamic fluctuations and the dissipations are mutually related by the fluctuation-dissipation relation (FDR) [43], it is natural and indispensable to consider hydrodynamic fluctuations in the hydrodynamic models for the highly non-equilibrium dynamics of the heavy-ion collisions [44]. Previous studies have shown that causal hydrodynamic fluctuations [45][46][47][48][49] play an important role to explain the centrality dependence of the rapidity decorrelation of anisotropic flows at LHC [50,51]. Since hydrodynamic fluctuations disturb fluid evolution randomly, they are expected to increase fluctuations of v n and affect v n {2}, especially in ultracentral collisions. Also, hydrodynamic fluctuations are more likely to affect higher-order v n {2}, which are related to smaller structures [48,52], and thus could be a key to resolving the puzzle.
We use the natural units ℏ = c = k B = 1 and the sign convention of the metric g µν = diag(+, −, −, −) throughout this Letter.
We employ the (3+1)-dimensional integrated dynamical model [48,50] with hydrodynamic fluctuations to describe the space-time evolution of high-energy heavy-ion collisions. In the integrated dynamical model, as prescribed in Ref. [8], we combine three models corresponding to the three stages of the collision reactions: a Monte-Carlo version of the Glauber model (MC-Glauber) [53] smoothly extended in the longitudinal direction [54] for the initial entropy deposition, a relativistic fluctuating hydrodynamic model rfh [47] which is relativistic dissipative hydrodynamics with hydrodynamic fluctuations for the space-time evolution of locally thermalized matters, and a hadron cascade model JAM [55] for the microscopic transport of hadron gases.
The main dynamical equations of the causal second-order fluctuating hydrodynamics are the conservation law of the energy-momentum tensor T µν of fluids and the constitutive equation for dissipative currents. We neglect the conservation law of the baryon number since we focus on heavy-ion collisions at LHC energies, at which the net baryon number is almost negligible around midrapidity. The energy-momentum tensor T µν can be tensor-decomposed as follows: where e, P, and π µν are the energy density, the pressure, and the shear-stress tensor, respectively. The flow velocity u µ in the Landau frame is defined as T µ ν u ν = eu µ , and ∆ µν g µν − u µ u ν is a projector for four-vectors onto the components transverse to u µ . We do not consider the bulk pressure in the present study. For an equation of state, we adopt s95p-v1.1 [56] which smoothly combines the equation of state from (2+1)-flavor lattice QCD simulations with that from the hadron resonance gas model corresponding to the hadron cascade model JAM. The constitutive equations for π µν in this study are [47,49,57] τ π ∆ µν αβ u λ ∂ λ π αβ + 1 + 4 3 τ π ∂ λ u λ π µν = 2η∆ µν αβ ∂ α u β + ξ µν π , (3) where η and τ π are the shear viscosity and the relaxation time, respectively. The tensor ∆ µν αβ is a projector for second-rank tensors onto the symmetric and traceless components transverse to u µ . The stochastic term ξ µν π represents hydrodynamic fluctuations, whose magnitude is determined by the FDR. In the Milne coordinates x, y , the FDR is written as where T is the temperature, and ⟨· · · ⟩ means the event average. In the actual calculations, a spatial regularization is needed to tame the ultra-violet divergences caused in the non-linear hydrodynamic equations [47]. We employ the smearing of the noise fields by the Gaussian kernel of the widths λ ⊥ and λ η s in the transverse and longitudinal directions, respectively. Smaller spatial cutoff parameters result in larger effects of hydrodynamic fluctuations.
For the initial conditions of hydrodynamic simulations, we generate the event-by-event entropy density distributions s (η s , x ⊥ ; τ 0 ) at the fixed hydrodynamic initial proper time τ 0 using mckln [8]. We parametrize the initial transverse profile from the linear combination of the participant number density of nuclei A and B, ρ A part (x ⊥ ) and ρ B part (x ⊥ ), respectively, and the number density of the binary collisions, ρ coll (x ⊥ ), for a randomly sampled impact parameter b satisfying P(b)db ∝ bdb, as follows, where C and α are the normalization factor and the hard fraction, respectively. The initial transverse profiles are extended to the longitudinal direction using the modified Brodsky-Gunion-Kühn (BGK) model based on the idea of the rapidity trapezoid [54,58,59]. In the present study, we neglect the initial transverse flow, the initial shear-stress tensor, and the initial longitudinal fluctuations. At a switching temperature T sw , we switch the description from the hydrodynamics to the microscopic kinetic theory using the Cooper-Frye formula [60] with a viscous correction [61,62]. The subsequent space-time evolution of hadron gases including hadronic rescatterings and decays of resonances is described by the hadron cascade model JAM [55]. Let us summarize the parameters of the present study. As in the previous calculations [8,48], we set the initial proper time τ 0 = 0.6 fm and the switching temperature T sw = 155 MeV. For the transport properties of QGP, we choose the specific shear viscosity η/s = 1/4π [63] or 1/2π and the relaxation time τ π = 3η/sT [57,64]. We use the same initial parameters C/τ 0 and α as the previous study [50] for each hydrodynamic model to reproduce the centrality dependence of the experimental charged-particle multiplicity measured by the ALICE Collaboration [65]. These parameters for each hydrodynamic model are summarized in Table 1.
Because of the large computational cost of the event-by-event dynamical simulations including the (3+1)-dimensional hydrodynamics and hadronic cascades, it is impractical to perform the minimum-biased simulations and afterward select the 0-0.2% events. The simplest approach to effectively generate the ultra-central events would be to fix the impact parameter b to 0 fm, but collisions with exact b = 0 fm do not occur in reality. Another approach would be to carry out the centrality selection at the initial stage using, e.g., the total entropy of the initial condition. However, even with a fixed initial entropy, the measured centrality can fluctuate due to the non-trivial evolution in the hydrodynamic and hadronic stages, the finite acceptance, etc [66]. In this study, we employ the importance sampling method to reduce the overall simulation cost while determining the centrality in the final multiplicity as in experiments. We introduce a weight function w(b) = exp (−b) into the impact parameter distribution as P(b)db ∝ b · w(b)db and generate initial conditions with small b intensively. To cancel the artifact of the extra weight in the distribution so that the redistributed events reproduce the proper multiplicity distribution, we need reweighting in the statistical analyses as where ⟨ · · · ⟩ ev is an event average, and X i is a physical quantity in the event i with impact parameter b i . Assuming that events with b larger than 5 fm do not contribute to centrality 0-0.2%, we generate only the initial conditions with b smaller than 5 fm. Figure 1 shows that the cross-section with b smaller than 5 fm is about 11.6% of the total cross-section within our version of the MC-Glauber model. Therefore, for the 0-0.2% centrality events, we can select the events that belong to the top 1.72% (= 0.2/11.6) of the reweighted multiplicity distribution. The resulting fraction of the 0-0.2% central- 11.6% ity events in the simulated number of events was about 10%. This means that the total number of simulated events to get the same number of simulated 0-0.2% events is reduced to about (1/10%)/(1/0.2%) = 1/50 compared to the case of minimumbiased simulations 1 . The impact parameter distribution of the centrality 0-0.2% events has a peak at ∼0.8 fm. The maximum impact parameter within about 500 simulated events of centrality 0-0.2% was less than 2.5 fm, which implies that the events of b > 2.5 fm essentially do not contribute to the 0-0.2% centrality class. This confirms that our choice of the initial cut of b = 0-5 fm is sufficiently large, and there is room to further reduce the initial cut below 5 fm.
Using the integrated dynamical model with the above importance sampling, we perform simulations of Pb+Pb collisions at √ s NN = 2.76 TeV. For each parameter set, we generate 5 000 hydrodynamic events and perform 10 independent particlizations and hadronic cascades for each hydrodynamic event, which ends up with 50 000 events in total. We apply the multiplicity cut and obtain about 5 000 events as the 0-0.2% central events.
The anisotropic flow coefficients from the two-particle cumulant method are calculated as where ∆ϕ represents the difference of the azimuthal angles between two charged hadrons, and the inner ⟨ · · · ⟩ represents the average over the particle pairs in each event. To compare the result with the data of the CMS Collaboration [23], we pick charged particles in the pseudo-rapidity range |η p | < 2.4 and the transverse momentum range 0.3 < p T < 3.0 GeV and calculate flows introducing a pseudo-rapidity gap |∆η p | min .
where r and ϕ are the radius and the azimuthal angle, respectively, with the origin being the center-of-mass of the entropy density distribution, and Φ n defines the nth-order participantplane angle. In Fig. 2, we first notice that the linear relation is not perfect on an event-by-event basis but has fluctuations, v n = κ n ε n + δ n . This means that an initial-state analysis based on the linear relation between v n {2} and ε n {2} = ⟨ε 2 n ⟩ ev would suffer from the flow fluctuations as v n {2} 2 ∼ κ 2 n ε n {2} 2 + ⟨δ 2 n ⟩ ev , which calls for the necessity of the event-by-event dynamical calculations. In Fig. 2, we also see that hydrodynamic fluctuations increase the flow fluctuations of v n and result in fatter distribution with smaller values of the Pearson correlation coefficient r, which means that the linear relation v n = κ n ε n becomes even worse. The hydrodynamic fluctuations also increase the intercept of the fitted line in Fig. 2, i.e., the final anisotropy v n is generated by the hydrodynamic fluctuations even with the vanishing initial anisotropy ε n = 0. These results suggest that the hydrodynamic fluctuations have a large impact on the linear relation v n = κ n ε n , which are typically assumed in the flow-puzzle discussions based on the initial anisotropies, and thus we expect that the consideration of the hydrodynamic fluctuations would possibly change the situation. Figure 3 shows the anisotropic flow coefficients v n {2} with |∆η p | > 2.0 in Pb+Pb collisions at √ s NN = 2.76 TeV for 0-0.2% centrality. The overestimation of v n {2} seen in Fig. 3 is the known one in existing studies 2 . In Fig. 3, the hydrodynamic fluctuations are found to increase v n {2} due to the increased flow fluctuations of v n while the shear viscosity reduces the magnitudes of v n {2}. The increase by the hydrodynamic fluctuations is more significant for larger shear viscosity η as naively expected from the FDR (4). We also observe in Fig 3 that the relative increase of v n by the hydrodynamic fluctuations is more significant in higher orders, which can be understood by the nature of the hydrodynamic fluctuations affecting smaller structures more. This behavior qualitatively makes v 3 {2} closer to v 2 {2} and can improve the situation for the flow puzzle. However, the effect is quantitatively too small to make v 3 {2} have the same magnitude as v 2 {2}, and thus this does not quantitatively solve the flow puzzle. To see the qualitative effects of hydrodynamic fluctuations in detail, we analyze the pseudo-rapidity gap |∆η p | min dependence of v n {2} in Fig. 4. The results from the viscous hydro are nearly independent of |∆η p | min because the present model does not contain major physical sources of non-flow correlations such as jets. On the other hand, in the case of fluctuating hydro, the influences of the hydrodynamic fluctuations on v n {2} become smaller as |∆η p | min increases, and v n {2} approaches that of the viscous hydro. This is because the hydrodynamic fluctuations have larger effects in a short range, and taking large |∆η p | min reduces their effects. Nevertheless, even with a large |∆η p | min , there remain sizable effects of the hydrodynamic fluctuations, especially with large shear viscosity. This means that taking the rapidity gap does not totally remove the effect of the hydrodynamic fluctuations. Figure 5 shows the v 2 {2}/v 3 {2} ratios calculated from the results shown in Fig. 4. With the shear viscosity η/s = 1/4π, the ratios with the viscous and fluctuating hydro have mostly the same value. With the increased shear viscosity η/s = 1/2π, the ratio increases slightly with the viscous hydro but decreases slightly with the fluctuating hydro. Also, the ratio approaches the experimental value with smaller |∆η p | min , i.e., with larger effects of the hydrodynamic fluctuations. The ratios of |∆η p | > 2 are also compared to the experimental value v 2 {2, |∆η p | > 2}/v 3 {2, |∆η p | > 2} ∼ 1.02 in Ref. [23]. With the large shear viscosity η/s = 1/2π, the deviation of the ratio ∼ 1.39 − 1.02 = 0.37 in the viscous hydro is reduced to ∼ 1.32 − 1.02 = 0.30 in the fluctuating hydro by about 19%. Likewise, with the smaller shear viscosity η/s = 1/4π, the deviation is reduced by about 3%, which is though not statistically significant. This means that hydrodynamic fluctuations qualitatively contribute to resolving the ultra-central flow puzzle, especially with larger shear viscosity. Nevertheless, it seems hard to quantitatively resolve the puzzle solely by the hydrodynamic fluctuations. In this Letter, we investigated the effects of hydrodynamic fluctuations during the space-time evolution of the QGP fluid for the ultra-central flow puzzle using a (3+1)-dimensional integrated dynamical model with relativistic hydrodynamic model rfh which includes causal hydrodynamic fluctuations and dissipation. We used importance sampling to intensively generate ultra-central collision events. We performed simulations of 0-0.2% central Pb+Pb collisions at √ s NN = 2.76 TeV with and without hydrodynamic fluctuations for comparison. We showed that the hydrodynamic fluctuations worsen the linear relation v n = κ n ε n in the ultra-central collisions on an event-by-event basis so that the dynamical simulations are important for the quantitative analysis of the ultra-central flow puzzle. We found that the hydrodynamic fluctuations increase the anisotropic flow coefficients v n {2} by increasing the flow fluctuations. The hydrodynamic fluctuations also make the v 2 {2}/v 3 {2} ratio closer to the experimental data, specifically by about 19% with a larger shear viscosity η/s = 1/2π. The effects of hydrodynamic fluctuations on v n {2} decrease with increasing pseudo-rapidity gap |∆η p | min , yet there remain sizable effects of the hydrodynamic fluctuations even with a large rapidity gap |∆η p | min . Our analyses show that hydrodynamic fluctuations qualitatively improve the situation for the ultra-central flow puzzle, though their effects alone are too small to solve the puzzle within the present model. In the future, we shall investigate the effect of hydrodynamic fluctuations with temperaturedependent shear viscosity, where we expect a larger effect coming from larger viscosity at high temperatures. The ultra-central flow puzzle remains a challenge for hydrodynamic models even today, but the hydrodynamic fluctuations would certainly contribute to an improvement in the ultra-central flow ratio.