Mode-by-mode fluid dynamics for relativistic heavy ion collisions

We propose to study the fluid dynamic propagation of fluctuations in relativistic heavy ion collisions differentially with respect to their azimuthal, radial and longitudinal wavelength. To this end, we introduce a background-fluctuation splitting and a Bessel-Fourier decomposition of the fluctuating modes. We demonstrate how the fluid dynamic evolution of realistic events can be build up from the propagation of individual modes. We describe the main elements of this mode-by-mode fluid dynamics, and we discuss its use in the fluid dynamic analysis of heavy ion collisions. As a first illustration, we quantify to what extent only fluctuations of sufficiently large radial wave length contribute to harmonic flow coefficients. We find that fluctuations of short wave length are suppressed not only due to larger dissipative effects, but also due to a geometrical averaging over the freeze-out hyper surface. In this way, our study further substantiates the picture that harmonic flow coefficients give access to a coarse-grained version of the initial conditions for heavy ion collisions, only.

In nucleus-nucleus collisions at the LHC and at RHIC, the dependence of soft hadron spectra on transverse momentum, on azimuthal orientation, on centrality and on particle species is understood since recently as fluid dynamic response to fluctuating initial conditions [1,2,3,4,5], for reviews see ref. [6,7]. This success of a fluid dynamic description is significant mainly for two reasons. First, the high sensitivity of fluctuations to dissipative properties of the produced fluid implies that fluctuations are promising tools for constraining the transport properties of dense QCD matter with unprecedented accuracy [8,9]. Second, since minimal dissipation implies maximal transparency to fluctuations, fluctuations in the initial stage of the collision can survive the time evolution. Therefore, the analysis of fluctuations may give access to the initial pre-equilibrium state and its fast evolution towards local equilibrium [10,11].
Motivated by these perspectives, many recent works have explored the dynamical relation between fluctuating initial conditions and experimentally accessible data. One important line of research is to characterize initial conditions for event averages and event ensembles in terms of eccentricities or closely related cumulants of the initial (entropy) density distribution [13,14], and to propagate entire events in viscous fluid dynamic simulations to the hadronic final state [9,12,15,16,17,18,19,20,21,22]. To date, this approach provides the most detailed test for the validity of a fluid dynamic description of heavy ion collisions. Despite this success of a cumulant-based characterization of initial conditions, several reasons motivate to explore alternative ones [23,24,25,26,27,28,29]. First, it is a problem well-known in probability theory that while any positive transverse density distribution can be characterized uniquely in terms of its infinite set of moments or cumulants, it is not possible to find (beyond the cumulants that determine a Gaussian) a positive density configuration corresponding to a finite set of cumulants such that all higher ones vanish. Strictly speaking, this implies that single cumulants cannot be propagated in fluid dynamics. Second, it is unknown how to extend a cumulant expansion to vector and tensor fields, as is needed e.g. if one wants to explore the natural possibility that fluctuations are manifest not only in the initial densities but also in the velocity field and shear viscous tensor. Finally, each cumulant receives typically contributions from fluctuations on various different wavelengths. There are advantages in decomposing initial fluctuations in an orthonormal basis of modes, but such bases have been used so far only in studies that formulate fluid dynamic perturbations on top of simple analytically known background fields with extended symmetries [24,25,26,27,29].
In a compagnon work [23], we have discussed how to characterize initial conditions in an orthornormal basis for scalar densities as well as vector and tensor fluctuations, constructed such that single fluctuating modes define positive densities and can therefore be propagated individually, mode-by-mode. In the present letter, we provide the first application of such a mode-by-mode fluid dynamics to realistic initial conditions, equation of state and transport properties. More specific, we characterize event samples as probability distributions over a basis of modes, we propagate each mode individually, and we build up experimental observables as superpositions of the individually propagated modes. We comment on how this can improve our understanding of why specific fluctuations in the initial state survive or do not survive the dynamic evolution.
Initial conditions. The dynamical evolution is initialized by specifying fluid dynamic fields on some hyper-surface, usally taken at fixed τ = √ t 2 − z 2 = τ 0 . Most generally, a model of the initial conditions is then defined in terms of a (functional) probability distribution p τ0 of the energy density or enthalpy w = + p, the fluid velocity u µ , the shear stress π µν , the bulk viscous pressure π bulk and possibly other fluid dynamic fields at τ 0 To discuss properties of the distribution p τ0 , we focus on one fluid dynamic field only, the enthalpy w. We consider fluctuations around a smooth background field w BG that is boost invariant and azimuthally symmetric. Longitudinal and azimuthal fluctuations in w are characterized by a standard Fourier expansion, For notational simplicity, we assume in what follows longitudinal boost invariance, i.e., we neglect any dependence on k η . If needed, our discussion is extended easily to the case of a non-trivial k η -dependence. The radial dependence is expanded in terms of Bessel functions J m that have appropriate boundary conditions at r = 0, Here, the radial wave vectors k of J m (z) and an overall scale R (R = 8 fm for the results presented here). The main difference between (3) and the Bessel-Fourier expansion proposed first in [28] is that we include the normalization factor w BG (τ 0 , r) on the right hand side. This ensures that the enthalpy density is positive everywhere even when only one or a few of the coefficientsw   identify the background field in (2) with the event average w BG ≡ w . Also at finite b, it can be advantageous to choose an azimuthally symmetric w BG even though this symmetry is not realized statistically; one can define e.g. w BG as the average over azimuthally randomized events. The azimuthal dependence of the event sample is then encoded in the event-averaged Bessel coefficients w (m) l = w (m) l that can take non-zero values for even integers m when b = 0. We have tested in model studies that the ansatz (2), (3) is as accurate for classes of semiperipheral collisions, as for central ones [23].
Since the coefficientsw (m) l characterize single events fully, the functional probability density p τ0 becomes a function of a set of numbersw (m) l . We have established [23] that for currently used models of initial conditions, p τ0 satisfies to good approximation the properties of a Gaussian probability distribution. For b = 0, Thus, p τ0 is fully characterized in terms of w BG ≡ w and the two-point correlators Fig. 2 shows the m = 2 two-point correlators (5) for the Monte Carlo Glauber model described in [12,23]. From these data (for all m < m max ), event samples of initial conditions can be generated easily. Since a modew on different radial length scale decorrelate as they are separated in scale.
Dynamic evolution. The above classification of initial conditions introduces naturally a background-fluctuation split- of all fluid fields. Instead of solving the relativistic fluid dynamic equations for the fields w, u µ etc. event-by-event, we solve for the smooth non-fluctuating background field once and for all, and we propagate the full basis of initial fluctuating modes with wave-numbers (l, m) as perturbations on this background field.
Relativistic viscous fluid dynamic solutions of eventaveraged background fields are well-documented. We follow ref. [15] in using the equation of state s95p-PCE which combines lattice QCD results at high temperatures with a hadron resonance gas at low temperatures. It implements also a chemical freeze-out at T = 165 MeV/k B . The default value of the shear viscosity to entropy ratio is η/s = 0.08 /k B and the corresponding relaxation time τ Shear = 0.23 /(k B T ). The evolution is initialized at τ 0 = 0.6 fm/c with initial flow and shear stress fields corresponding to the Navier-Stokes form of a longitudinal Bjorken expansion. The background enthalpy is initialized as w BG = w at b = 0. The entropy in both background and fluctuations scale with (1−x)N part /2+x N coll , x = 0.118, as in ref. [15]. Fig. 3 shows the freeze-out curves resulting from fluid dynamic evolution of this azimuthally symmetric background field. They are consistent with published benchmarks.
The evolution equations for the fluctuations w F , u µ F etc. depend on the solution for the background fields. They become particularly simple if treated as small perturbations that can be linearized. For a given Fourier mode specified by m and k η , the evolution equation becomes then 1 + 1 dimensional and with the Bessel expansion as in eq. (3) it reduces for all τ to a set of coupled ordinary differential equations which we solve numerically. This set-up extends the strategy of Refs. [24,27] to arbitrary background fields and arbitrary classes of initial fluctuations, including initial fluctuations in the fluid velocity [27] or shear. In fig. 3, we show for the normalized enthalpy densityw = δw/w BG the spatial evolution for three modes of different radial wave number l. One sees that the viscous damping increases significantly for shorter radial wave-length, thus illustrating the importance of studying the effect of fluctuations differentially in l. We also find that the viscous damping seen in fig. 3 increases strongly with η/s. On the other hand, modes with larger l lead to more strongly oscillating patterns on the freeze-out surface and have therefore less impact on particle spectra even for η/s = 0. Freeze-out and particle spectra.. Hydrodynamics ceases to apply when interaction rates become too small to maintain local kinetic equilibrium. We assume that this happens when the background field drops below T fo = 120 MeV. Particle distributions then freeze out. We determine them using the standard Cooper-Frye prescription neglecting resonance decays and hadronic rescatterings after freezeout. (In principle these effects could be incorporated by solving the corresponding kinetic equations for the background and, in linearized form, for the fluctuations.) The occupation numbers on the freeze-out surface are taken to be of ideal gas form with chemical potentials according to the equation of state s95p-PCE. Viscous corrections due to shear stress are approximated by the quadratic ansatz [34]. Due to its azimuthal rotation invariance, the background field contributes to the φ-and y-independent part of the one-particle spectrum S(p T ) = dN/(p T dp T dφdy), only. If fluctuations on top of this background are small enough, their effect on particle spectra can be linearized, ln dN single event p T dp T dφdy = ln S 0 (p T ) + m,lw Here, the functions θ (m) l (p T ) determine how the fluctuating modes of wave-numbers (m, l) contribute to the hadronic spectrum. In general, the θ (m) l depend also on rapidity and particle species. They are calculated as follows: The linearized hydrodynamical evolution equations on top of the background solution are solved for the initial condition corresponding to the mode (m, l) in enthalpy density. All fluid fields resulting from this initialization are determined on the freeze-out surface and the corresponding contribution to the particle spectrum is determined from an appropriate linearization of the Cooper-Frye formula. Dividing finally by the background contribution to the particle spectrum yields θ (m) l (p T ). One-particle spectra of event samples are obtained by averaging (6) over the probability distribution p τ0 ({w In close analogy, the calculation of two-particle correlations requires knowledge about the initial correlations between pairs of modesw (m) l1w (m) * l2 whose contribution to the hadronic two-particle spectra is then determined by the product θ . The double differential harmonic flow coefficient for event samples reads then to lowest order inw (7) The single differential harmonic flow coefficients v m (p T ) can be obtained from (7) as appropriately weighted p T integrals. Note that in close analogy to the experimental procedure of extracting harmonic flow coefficients, (7) does not invoke knowledge of a reaction plane but determines v m from the azimuthal dependence of two-particle correlations that have their dynamic origin in the azimuthal correlations w }) can be studied by simple matrix multiplication, see (6), (7). We note that this formulation assumes a linear relation between fluctuating modes at τ 0 and hadronic distributions at freeze-out. One could test the accuracy of such a linear relation by comparing for selected events to results from full fluid dynamic simulations. This may also allow to identify characteristic signatures of non-linear fluid dynamic behavior in heavy ion collisions which would be interesting in itself.
Another possibility to estimate the effects of non-linear terms in the hydrodynamical evolution and at freeze-out is to treat them order-by-order in a perturbative expansion. The leading order in this perturbation theory for small fluctuations around a smooth but dynamically evolving background is the linear order presented in this paper. At next-to-leading or quadratic order one can study for instance how an m = 2 and an m = 3 mode interact and how this contributes to a signal for v 5 . A more detailed discussion of this kind of perturbative treatment is left for a future publication.
In fig. 4, we compare calculations of hadronic spectra and flow coefficients from mode-by-mode hydrodynamics to data for central Pb+Pb collisions taken by the ALICE Collaboration [31,32]. For the one-particle spectra of pions, kaons and protons, we find that fluctuations make a very small contribution to (6), so that the model results shown in Fig. 4 are very close to the result for an initial condition defined by the smooth background field without fluctuations. For the evolution of a smooth background field, we have checked on the level of the freeze-out hyper surface and on the level of spectra that our simulations are quantitatively consistent with the hydrodynamic benchmarks established by the TECHQM Collaboration. Our study does not take into account resonance decays and hadronic rescatterings after freeze-out (we do not switch to a hadron cascade model), and we made no attempt to improve agreement with data by optimizing input parameters of the fluid dynamic simulation such as η/s or the equilibration time τ 0 . It is therefore no surprise that one sees some differences between data and numerical results for one-particle spectra. However, these one-particle spectra are sufficiently well reproduced to serve as baseline for the present study, whose purpose is to illustrate how in a mode-by-mode fluid dynamic analysis, results on elliptic, triangular, 4-th and 5-th order flow are built up in terms of the contributions of individual fluctuations of characteristically different radial wavelength. With this respect, our main conclusion is that the sum (7) converges quickly, for small radial wave numbers l ≤ l max ≈ 5. This means that only fluctuations of sufficiently large radial wavelength matter for the dynamics of flow coefficients. For the density distribution of the specific event shown in fig. 1, for instance, it is then only the coarse-grained information shown for l max = m max = 5 in the upper right panel of fig. 1 that affects the value of flow harmonics in fig. 4. Since we observe this rapid convergence in l for minimal dissipative effects (η/s = 1/4π), we expect this finding to be more general than the specific model study in which we have established it here.
The precise numerical values for v n (p T ) will depend in general on the weights and correlations w of the different fluctuating modes in the initial conditions, on the input parameters of the fluid dynamic evolution, and on the treatment of rescattering effects and resonance decays after freeze-out. The present study does not optimize input parameters and it does not account for physics effects after freeze-out. Also, it is limited to the one set of fluctuating initial conditions characterized in fig. 2. Within this non-optimized setting, we find a very reasonable agreement with data for v 2 , v 3 , v 4 and v 5 in the range p T ≤ 1 GeV, while experimental data for v 2 in the range 1 ≤ p T ≤ 3 GeV lie significantly below our calculation. A full exploration of this significant, and other smaller discrepancies between data and calculation will require to op-  Particle spectra for pions, kaons and protons S(p T ) = dN/(2πp T dp T dy), calculated from (6) and compared to data on 5% most central events [31]. Following panels: Elliptic (v 2 ), triangular (v 3 ), 4-th and 5-th order flow of charged particles, calculated for contributions from π ± , K ± , p andp from (7) with lmax = 1, 2, 5, 10, 20 and compared to data on 2% most central events. The dashed curve in the plot of v 2 shows results (for lmax = 20) for a modified event sample in which the weight of one particular fluctuating mode,w is decreased by a factor 0.7.
timize the input parameters which lies outside the scope of the present study. We note, however, that a simple mild rescaling of the weight of one single fluctuating mode in the initial conditions,w (2) 1 → 0.7w (2) 1 , can improve agreement between simulation results and data for v 2 over a much increased p T -range. While this curious observation clearly does not replace a full optimization of all input parameters, it illustrates that mode-by-mode hydrodynamics offers the possibility of "backward engineering" of initial conditions: for any given dynamics and a set of data, eqs. (6) and (7) allow to optimize the correlators w , and thus the event distribution p τ0 . Further differential test of mode-by-mode fluid dyanamics will also include the study of particle-identified flow harmonics. Fig. 5 shows corresponding results for pions, kaons and protons. Close inspection shows that the curves are ordered with the particle mass at small p T according to v m (p T ) Protons < v m (p T ) Kaons < v m (p T ) Pions , while for larger p T the ordering is reversed.
The proposed set-up may also be interesting in the context of the recently proposed "event shape engineering" [33]. Namely, it allows easily for the calculation of event distributions in one-and two-particle spectra from p τ0 , and it thus allows to study the relations between cuts on event distributions and cuts on initial conditions p τ0 .
Since it offers such possibilities, we expect that the proposed fully differential treatment of fluctuations will become a helpful tool used to fully exploit the experimental precision in heavy ion physics.  Figure 5: The transverse momentum dependent n-th order flow coefficients vn, as in Fig. 4, but calculated differentially for contributions from π ± , K ± , p andp for central events.