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

###### Abstract

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 Alver:2010gr (); Mishra:2007tw (); Broniowski:2007ft (); Sorensen:2008zk (); Takahashi:2009na (), for reviews see ref. Heinz:2013th (); Gale:2013da (). 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 accuracyQiu:2011iv (); Schenke:2011bn (). 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 Bhalerao:2011yg (); Schenke:2012wb ().

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 Teaney:2010vd (); Teaney:2012ke (), and to propagate entire events in viscous fluid dynamic simulations to the hadronic final state Schenke:2011bn (); Holopainen:2010gz (); Qiu:2011hf (); Gardim:2011xv (); Qian:2013nba (); Petersen:2012qc (); Deng:2011at (); Bhalerao:2011bp (); Gale:2012rq (); Niemi:2012aj (). 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 Floerchinger:2013vua (); Staig:2010pn (); Staig:2011wj (); Gubser:2010ui (); Florchinger:2011qf (); ColemanSmith:2012ka (); Springer:2012iz (). 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 Staig:2010pn (); Staig:2011wj (); Gubser:2010ui (); Florchinger:2011qf (); Springer:2012iz ().

In a compagnon work Floerchinger:2013vua (), 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 . Most generally, a model of the initial conditions is then defined in terms of a (functional) probability distribution of the energy density or enthalpy , the fluid velocity , the shear stress , the bulk viscous pressure and possibly other fluid dynamic fields at

(1) |

To discuss properties of the distribution , we focus on one fluid dynamic field only, the enthalpy . We consider fluctuations around a smooth background field that is boost invariant and azimuthally symmetric. Longitudinal and azimuthal fluctuations in are characterized by a standard Fourier expansion,

(2) |

For notational simplicity, we assume in what follows longitudinal boost invariance, i.e., we neglect any dependence on . If needed, our discussion is extended easily to the case of a non-trivial -dependence. The radial dependence is expanded in terms of Bessel functions that have appropriate boundary conditions at ,

(3) |

Here, the radial wave vectors are set by the -th zeroes of and an overall scale ( fm for the results presented here). The main difference between (3) and the Bessel-Fourier expansion proposed first in ColemanSmith:2012ka () is that we include the normalization factor on the right hand side. This ensures that the enthalpy density is positive everywhere even when only one or a few of the coefficients are non-vanishing. The azimuthal and radial wavenumber and can be restricted to the ranges and when the spatial resolution is bound. Lemoine’s discrete Bessel transformation provides a CPU-inexpensive method for determining Lemoine (); Floerchinger:2013vua (). Fig. 1 illustrates for a phenomenologically relevant enthalpy density that fluctuations in a single event can be characterized satisfactorily in terms of a small set of Bessel coefficients in (2), (3).

Event samples can have statistical symmetries that are broken by event-wise fluctuations. For instance, a sample at vanishing impact parameter will have statistical azimuthal symmetry. In this case, we choose to identify the background field in (2) with the event average . Also at finite , it can be advantageous to choose an azimuthally symmetric even though this symmetry is not realized statistically; one can define e.g. as the average over azimuthally randomized events. The azimuthal dependence of the event sample is then encoded in the event-averaged Bessel coefficients that can take non-zero values for even integers when . We have tested in model studies that the ansatz (2), (3) is as accurate for classes of semi-peripheral collisions, as for central ones Floerchinger:2013vua ().

Since the coefficients characterize single events fully, the functional probability density becomes a function of a set of numbers . We have established Floerchinger:2013vua () that for currently used models of initial conditions, satisfies to good approximation the properties of a Gaussian probability distribution. For ,

(4) |

Thus, is fully characterized in terms of and the two-point correlators

(5) |

Fig. 2 shows the two-point correlators (5) for the Monte Carlo Glauber model described in Holopainen:2010gz (); Floerchinger:2013vua (). From these data (for all ), event samples of initial conditions can be generated easily. Since a mode corresponds to a radial wavelength that decreases with increasing , Fig. 2 shows how fluctuations 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 splitting , etc. of all fluid fields. Instead of solving the relativistic fluid dynamic equations for the fields , 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 as perturbations on this background field.

Relativistic viscous fluid dynamic solutions of event-averaged background fields are well-documented. We follow ref. Qiu:2011hf () 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 . The default value of the shear viscosity to entropy ratio is and the corresponding relaxation time . The evolution is initialized at with initial flow and shear stress fields corresponding to the Navier-Stokes form of a longitudinal Bjorken expansion. The background enthalpy is initialized as at . The entropy in both background and fluctuations scale with , , as in ref. Qiu:2011hf (). 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 , 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 and , the evolution equation becomes then 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. Staig:2010pn (); Florchinger:2011qf () to arbitrary background fields and arbitrary classes of initial fluctuations, including initial fluctuations in the fluid velocity Florchinger:2011qf () or shear. In fig. 3, we show for the normalized enthalpy density the spatial evolution for three modes of different radial wave number . 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 . We also find that the viscous damping seen in fig. 3 increases strongly with . On the other hand, modes with larger lead to more strongly oscillating patterns on the freeze-out surface and have therefore less impact on particle spectra even for .

## 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 . Particle distributions then freeze out. We determine them using the standard Cooper-Frye prescription neglecting resonance decays and hadronic rescatterings after freeze-out. (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 Teaney:2003kp (). Due to its azimuthal rotation invariance, the background field contributes to the - and -independent part of the one-particle spectrum , only. If fluctuations on top of this background are small enough, their effect on particle spectra can be linearized,

(6) |

Here, the functions determine how the fluctuating modes of wave-numbers contribute to the hadronic spectrum. In general, the 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 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 .

One-particle spectra of event samples are obtained by averaging (6) over the probability distribution . In close analogy, the calculation of two-particle correlations requires knowledge about the initial correlations between pairs of modes 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 in

(7) |

The single differential harmonic flow coefficients can be obtained from (7) as appropriately weighted 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 from the azimuthal dependence of two-particle correlations that have their dynamic origin in the azimuthal correlations between different fluctuating modes. In this way, once the functions are calculated for a given smooth background and the finite set of wave numbers , the fluid dynamic propagation of arbitrary samples of fluctuations can be studied by simple matrix multiplication, see (6), (7). We note that this formulation assumes a linear relation between fluctuating modes at 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 and an mode interact and how this contributes to a signal for . 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 Abelev:2012wca (); ALICE:2011ab (). 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 or the equilibration time . 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 . 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 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 for minimal dissipative effects (), 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 will depend in general on the weights and correlations 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 , , and in the range GeV, while experimental data for in the range GeV lie significantly below our calculation. A full exploration of this significant, and other smaller discrepancies between data and calculation will require to optimize 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, , can improve agreement between simulation results and data for over a much increased -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 , and thus the event distribution . 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 according to , while for larger the ordering is reversed.

The proposed set-up may also be interesting in the context of the recently proposed “event shape engineering” Schukraft:2012ah (). Namely, it allows easily for the calculation of event distributions in one- and two-particle spectra from , and it thus allows to study the relations between cuts on event distributions and cuts on initial conditions . 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.

## Acknowledgements.

S. F. acknowledges financial support by DFG under contract FL 736/1-1. We thank U. Heinz, M. Luzum, J.-Y. Ollitrault, H. Petersen and C. Shen for discussion.

## References

- (1) B. Alver and G. Roland, Phys. Rev. C 81 (2010) 054905 [Erratum-ibid. C 82 (2010) 039903].
- (2) A. P. Mishra, R. K. Mohapatra, P. S. Saumia and A. M. Srivastava, Phys. Rev. C 77 (2008) 064902.
- (3) W. Broniowski, P. Bozek and M. Rybczynski, Phys. Rev. C 76 (2007) 054905.
- (4) P. Sorensen [STAR Collaboration], J. Phys. G 35 (2008) 104102.
- (5) J. Takahashi, B. M. Tavares, W. L. Qian, R. Andrade, F. Grassi, Y. Hama, T. Kodama and N. Xu, Phys. Rev. Lett. 103 (2009) 242301.
- (6) U. W. Heinz and R. Snellings, arXiv:1301.2826 [nucl-th].
- (7) C. Gale, S. Jeon and B. Schenke, Int. J. Mod. Phys. A 28 (2013) 1340011.
- (8) Z. Qiu and U. W. Heinz, Phys. Rev. C 84 (2011) 024911.
- (9) B. Schenke, S. Jeon and C. Gale, Phys. Rev. C 85 (2012) 024901.
- (10) R. S. Bhalerao, M. Luzum and J. -Y. Ollitrault, Phys. Rev. C 84 (2011) 034910.
- (11) B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 108 (2012) 252301.
- (12) H. Holopainen, H. Niemi and K. J. Eskola, Phys. Rev. C 83 (2011) 034901.
- (13) D. Teaney and L. Yan, Phys. Rev. C 83 (2011) 064904.
- (14) D. Teaney and L. Yan, Phys. Rev. C 86 (2012) 044908.
- (15) Z. Qiu, C. Shen and U. Heinz, Phys. Lett. B 707, 151 (2012).
- (16) F. G. Gardim, F. Grassi, M. Luzum and J. -Y. Ollitrault, Phys. Rev. C 85 (2012) 024908.
- (17) W. -L. Qian, P. Mota, R. Andrade, F. Gardim, F. Grassi, Y. Hama and T. Kodama, arXiv:1305.4673 [hep-ph].
- (18) H. Petersen, R. La Placa and S. A. Bass, J. Phys. G 39 (2012) 055102.
- (19) W. -T. Deng, Z. Xu and C. Greiner, Phys. Lett. B 711 (2012) 301.
- (20) R. S. Bhalerao, M. Luzum and J. -Y. Ollitrault, Phys. Rev. C 84 (2011) 054901.
- (21) C. Gale, S. Jeon, B. Schenke, P. Tribedy and R. Venugopalan, Phys. Rev. Lett. 110 (2013) 012302.
- (22) H. Niemi, G. S. Denicol, H. Holopainen and P. Huovinen, arXiv:1212.1008 [nucl-th].
- (23) S. Floerchinger and U. A. Wiedemann, arXiv:1307.7611 [hep-ph].
- (24) P. Staig and E. Shuryak, Phys. Rev. C 84 (2011) 034908.
- (25) P. Staig and E. Shuryak, Phys. Rev. C 84 (2011) 044912.
- (26) S. S. Gubser and A. Yarom, Nucl. Phys. B 846 (2011) 469.
- (27) S. Floerchinger and U. A. Wiedemann, JHEP 1111 (2011) 100.
- (28) C. E. Coleman-Smith, H. Petersen and R. L. Wolpert, arXiv:1204.5774 [hep-ph].
- (29) T. Springer and M. Stephanov, Nucl. Phys. A904-905 2013 (2013) 1027c [arXiv:1210.5179 [nucl-th]].
- (30) D. Lemoine, J. Chem. Phys. 101, 3936 (1994).
- (31) B. Abelev et al. [ALICE Collaboration], Phys. Rev. Lett. 109, 252301 (2012).
- (32) K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 107, 032301 (2011).
- (33) J. Schukraft, A. Timmins and S. A. Voloshin, Phys. Lett. B 719 (2013) 394.
- (34) D. Teaney, Phys. Rev. C 68 (2003) 034913.