Kinetic freeze-out, particle spectra and harmonic flow coefficients from mode-by-mode hydrodynamics

The kinetic freeze-out for the hydrodynamical description of relativistic heavy ion collisions is discussed using a background-fluctuation splitting of the hydrodynamical fields. For a single event, the particle spectrum, or its logarithm, can be written as the sum of background part that is symmetric with respect to azimuthal rotations and longitudinal boosts and a part containing the contribution of fluctuations or deviations from the background. Using a complete orthonormal basis to characterize the initial state allows one to write the double differential harmonic flow coefficients determined by the two-particle correlation method as matrix expressions involving the initial fluid correlations. We discuss the use of these expressions for a mode-by-mode analysis of fluctuating initial conditions in heavy ion collisions.


I. INTRODUCTION
Fluid dynamic models of relativistic heavy ion collisions at RHIC and at the LHC provide an overall good description of transverse momentum spectra and harmonic flow coefficients up to about p T 2 − 3 GeV [1][2][3][4][5][6][7][8][9][10][11][12][13], for reviews see Refs. [14][15][16]. The initialization of these models needs to account for event wise fluctuations [17][18][19][20][21][22], as they may arise either from quantum fluctuations or from the substructure of nuclei in terms of nucleons. Recently, we have developed an approach [23] that allows for a differential study of the effects of such fluctuations in the fluid dynamic fields. It is based on a mode-by-mode decomposition of fluctuations in single events, on a functional characterization of initial conditions for event classes and on a background-fluctuation splitting for the hydrodynamical evolution. Here, we give a technically complete documentation of how kinetic freeze-out of fluctuating modes is handled in this mode-by-mode hydrodynamics. In particular, we provide explicit expressions for the hadronic response to specific fluid dynamic fluctuations. These expressions are of more general interest since they can be used to establish to what extent specific fluid dynamic features are washed out or survive hadronization.
In general, a fluid dynamic description of an expanding system seizes to be valid at kinetic freeze-out, when densities become too low and microscopic interaction times become too long to maintain the system close to local equilibrium. This transition from fluid dynamic behavior to free-streaming particles occurs within a finite time after the collision and within finite volume around the collision point. The Cooper-Frye freeze-out prescription [24] is based on the assumptions that freeze-out occurs sufficiently rapidly to approximate it along a sharp three-dimensional freeze-out hyper surface Σ f , and that the occupation numbers on Σ f are given by free thermal single particle distribution functions 1 . Both assumptions may be questioned. In particular, one often employs in phenomenological applications a smoother freeze-out prescription in which the switch to particle distribution functions at fixed Σ f is followed by an intermediate regime of hadronic scatterings, described for example by Boltzmann's equation, see for example [25] for a discussion how this can be done within viscous hydrodynamics. This procedure is interesting since it can give insight into the relevance of a four-dimensional freeze-out volume and a regime of hadronic scattering. However, the usual choice for initializing this regime is again some Cooper-Frye freeze-out prescription, although it would be more consistent to base the initialization in this case on occupation numbers for interacting particles. At least in this sense, the currently used refinements of the Cooper-Frye freeze-out prescription involve additional model assumptions. In the present work, we shall study freeze-out of fluid dynamic fluctuations into free streaming particles without taking hadronic scattering into account, but we shall discuss finally how these effects could be taken into account and where they may be relevant. This paper is organized as follows: In section II, we first review the Cooper-Frye freeze-out condition for a fluid with finite bulk and shear viscosity. We focus in particular on azimuthally symmetric freeze-out conditions, but we emphasize already here that such highly symmetric freeze-out conditions are also relevant for collisions at finite impact parameter. The reason is that mode-by-mode hydrodynamics employs a background-fluctuation splitting of all fluid dynamic fields. The background part is chosen to be invariant under rotations in the transverse plane and longitudinal boosts, and the freeze-out hyper surface will inherit these symmetries since we take it as the surface of constant background temperature. This is slightly different from standard Cooper-Frye prescriptions where Σ f is defined as the surface of constant temperature, but the two descriptions can be mapped to each other to a given order. At linear order in perturbations of the fluid dynamic fields around the background they are actually identical as we show in appendix B. All azimuthal asymmetries, irrespective of whether they arise from event-wise fluctuations or from collisions at finite impact parameter, are then encoded in fluctuating modes that are frozen out on the azimuthally symmetric hypersurface Σ f . We shall discuss this procedure for the background fields and for the fluctuating fields in section III. In section IV we finally provide explicit expressions for a mode-by-mode decomposition of the fluctuating fields before shortly discussing our results in an outlook and conclusions section. Appendix A compiles some analytic expressions for rapidity and azimuthal integrals and we discuss the difference between a freeze-out at constant background temperature and one at constant total temperature in appendix B.

II. COOPER-FRYE KINETIC FREEZE-OUT AND OCCUPATION NUMBER
Quite generally we assume that at a space-time point x where the system is still close to local thermal equilibrium the occupation number for a specific particle species i takes the form Here, the single particle distribution function f i depends on the four-momentum p µ of the particle, and it depends on the position x through the values that the fluid dynamic fields take at x. In equation (1), we have written the fluid dynamic fields of temperature T , velocity u µ , shear stress π µν and bulk viscous pressure π bulk explicitly. Equivalently, we could have used reparametrizations of these fields, such as energy density, entropy or enthalpy instead of T . And if further fluid dynamic fields would be relevant for characterizing the fluid, the distribution function f i would depend on them as well. In particular, we neglect here and in what follows non-zero conserved baryon number or electric charge as well as electroagnetic fields, and we do not consider the dependence of occupation numbers on spin. The expression in (1) can be used both for (a set of) on-shell particles where E = m 2 i + p 2 , and for resonances with non-zero width. In that case, E = ν i + p 2 and the values of ν i are determined by the spectral density ρ i (ν i ).
Following the prescription of Cooper and Frye the particle spectrum after kinetic freeze-out is given as The integral is here over a three-dimensional hyper-surface Σ f where kinetic equilibrium is just still valid. This is sometimes refered to as the "surface of last scattering". The minus sign accounts for our choice of the metric with signature (−, +, +, +).
There is no precise theory on how the freeze-out surface is determined by physical principles. 2 The general idea is that it corresponds to the point in the evolution of a fluid element where the scattering rate between particles becomes too small to maintain kinetic equilibrium. This is characterized often by the condition that particle densities become smaller than a certain value corresponding to some freeze-out temperature T fo . However, a dynamical freeze-out criterion may be more physical, see for example [29] for a recent attempt in this direction. This is so since freeze-out should depend also on the history of the collision as can be understood when thinking about a hypothetical very slow expansion where thermal equilibrium would be maintained down to very small temperatures.
A. Freeze-out at constant background temperature As recalled in section IV, mode-by-mode hydrodynamics is based on a background fluctuation splitting of all fluid dynamic fields. Instead of the standard condition of freeze-out on a hypersurface of constant temperature, we shall √ t 2 − z 2 and radius r = x 2 1 + x 2 2 for a central Pb-Pb collision at LHC energy. We compare different choices for the shear viscosity to entropy ratio: η/s = 0 (gray, uppermost curve at small r), η/s = 0.08 (blue, middle curve at small r) and η/s = 0.3 (orange, lowermost curve at small r). The arrows indicate the direction of the fluid velocity at freeze out for the case η/s = 0.08.
work in the following with the related but slightly different prescription of using a hyper surface of constant background temperature. For a particular event the actual temperature as determined by the sum of background and fluctuating hydrodynamical fields can then vary on this surface, and the occupation numbers must be corrected accordingly. To any given order in fluctuations, a freeze-out at strictly constant temperature can be mapped to this description. In principle, one has to include correction terms that account for the change in particle spectra between the two surfaces.
To linear order in fluid dynamic fields these correction terms vanish as we show in appendix B. In any case, in the absence of a precise theory of freeze-out we do not see a physics argument that would prefer one of these closely related freeze-out conditions. In the following, we choose to absorb the entire non-trivial dependence of fluid dynamic fields on azimuthal angle and space-time rapidity in the fluctuating part of fluid dynamic fields. The background part of the fluid dynamic fields is then invariant under azimuthal rotations and Bjorken boost transformations. In general, the freeze out surface at constant temperature will depend on the azimuthal angle 0 ≤ ϕ < 2π, the space rapidity −∞ < η < ∞, the proper time τ and the radius r. However, for the symmetric background field considered here, freeze-out on a hyper surface of constant background temperature is characterized for all ϕ and η by the same one-parameter curve in the τ -r-plane, where without loss of generality we take α ∈ [0; 1]. Fig. 1 shows the freeze-out curve in the τ -r-plane at constant background temperature T 0 for a central Pb-Pb collision at LHC energies. The location of the freeze-out surface, and the values of all fluid dynamic fields along it, vary with the value of the shear viscosity over entropy ratio, see Fig. 1.
For the particle four-momentum we use a parameterization in terms of the transverse mass m 2 T = E 2 −p 2 3 , the transverse momentum p T = p 2 1 + p 2 2 , the momentum-space rapidity η P = arctanh(p 3 /E) and the (momentum) azimuthal angle φ = arctan(p 2 /p 1 ). Transverse mass and momentum are related by m 2 T = m 2 i + p 2 T or m 2 T = ν i + p 2 T for particles and resonances, respectively. In these coordinates, the freeze-out surface element that enters the Cooper-Frye formula (2) reads B. Distribution functions at freeze-out We now discuss the distribution function f i in Eq. (1). Close to thermal equilibrium, it is of the form where f i,eq accounts for the equilibrated part and δf i describes some deviation. The equilibrated part does not depend on the shear stress and the bulk viscous pressure. It is a Lorenz scalar and must therefore be of the form f i,eq = f i,eq (p µ u µ , T ). If kinetic freeze-out happens in a regime where chemical equilibrium is no longer maintained, then the equilibrated occupation numbers f i,eq also depend on chemical potentials µ i (T ). In agreement with our approximation to neglect all conserved currents apart from energy and momentum these chemical potentials are not independent in the thermodynamic sense but are functions of the temperature. If deviations δf i from equilibrium are sufficiently small, then they will be linear in the shear stress tensor and the bulk viscous pressure. Using that δf i is a Lorentz scalar and that u µ π µν = 0, one can write where ∆ µν = g µν + u µ u ν . The ansatz chosen here makes the functions g i and h i dimensionless. According to this ansatz, the most general distribution function f i in Eq. (1) can be written in terms of the three functions f i,eq , g i and h i . The explicit form of the distribution f i,eq is of course fixed from the condition of thermodynamic equilibrium. However, the functions g i and h i that characterize deviations will depend on the microscopic dynamics underlying equilibration processes. In principle, these functions can be determined from a proper, microscopic, quantum field theoretic calculation in thermal equilibrium using linear response theory. So far we made no assumptions about the size of interaction effects and the discussion was actually not constrained to the regime where perturbation theory or kinetic theory are valid. In particular one may use the formula derived so far to switch from the hydrodynamic description to a more microscopic scattering theory when a particular hyper surface, for example at a certain temperature, is crossed. In general, the equilibrated occupation numbers f i,eq on such a hyper surface would be those of an interacting system and would differ from the ideal gas form. We now leave these general considerations aside and specialize to a freeze-out happening at a point where interaction effects cease to be important. Then, the equilibrated occupation numbers are given by their ideal gas form, for bosons and fermions respectively. Also, the energy-momentum tensor at a point x is then linear in the occupation numbers for different particle species i, Identifying the equilibrated parts on both sides of this equation, one finds for the energy density and pressure For simplicity we have chosen here a reference frame where u µ = (1, 0, 0, 0). Similarly, projecting the non-equilibrated part of (11) to the transverse and traceless component yields and taking the trace one gets The form of δf in Eq. (9) makes sure that energy density and pressure are not corrected by the shear viscous part ∼ g i . For h i one has one additional constraint, which ensures that energy density is not modified. Comparison to (14) shows that h i = 0 for massless particles. As mentioned above, the specific form of g i will depend on the microscopic dynamics of equilibration processes. To what extent this form is model-dependent has been discussed for instance in Ref. [34]. In what follows, we shall use the so-called quadratic ansatz [33,34] For a single component gas one can check that with this choice (13) is fulfilled within a few percent for bosons, fermions and Boltzmann distributed particles. Technically, it is convenient to expand f i,eq in powers of the Boltzmann factor exp pν u ν +µi T . We will therefore work in the following with the occupation numbers for bosons and fermions, respectively.

III. PARTICLE SPECTRA AND BACKGROUND-FLUCTUATION SPLITTING
In section II, we have shown that the particle spectrum (2) can be written as where integration is over an azimuthally symmetric and Bjorken boost invariant hyper surface Σ f . The fluid dynamic fields u µ , T , π µν and µ i that appear in (18) are evaluated on Σ f . Choosing this freeze-out to occur at constant background temperature T 0 means that we perform a background-fluctuation splitting of the temperature into a position-independent background and a position-dependent fluctuation term, For all other fluid dynamic fields, the background terms can depend on the position along Σ f , but they share the symmetries of Σ f . We therefore write for the background-fluctuation splitting and likewise for other fluid dynamic fields such as π Bulk . The chemical potentials µ i are assumed to depend only on temperature T so that they can be written as in eq. (19). We recall that the ansatz (19)-(21) is valid also for collisions at finite impact parameter, since it allows for the parametrization of arbitrary ϕ-dependencies of the fluid fields. One way to see that the background-fluctuation splitting (19)-(21) is physically meaningful is to consider the background parts of all fields as characterizing the event average over many collisions (with random azimuthal orientation), while the fluctuating parts characterize event-specific fluctuations. For more technical details on this point, see Ref. [35]. The strategy of the following is to use the ansatz (19)- (21) to expand the spectrum (18) in fluctuations around the background fields. We do this by evaluating first the zeroth order background contribution to (18) in subsection III A. We then discuss in subsection III B the contributions that arise from first order in fluctuations. We anticipate that the final result of this excercise will be a set of equations that relate in a very explicit form specific modes of fluctuating hydrodynamic fields to specific features in experimental observables.

A. Zeroth order in fluctuations: the background contribution
To zeroth order in fluctuations, the temperature-dependence of the spectrum (18) reduces to a dependence on T 0 . As for the background fluid velocity u µ 0 , symmetries imply that only the temporal and radial component can take non-vanishing values in Bjorken coordinates τ, r, ϕ, η. The Lorentz scalar p µ u µ 0 in (18) hence takes the form For the background part of the shear stress tensor symmetries imply that only the components π τ τ 0 , π τ r 0 , π rτ 0 , π rr 0 , π ϕϕ 0 and π ηη 0 are non-zero. Due to the transverse and traceless constraints only two of them are actually independent. We choose a parameterization of the form (still in coordinates (τ, r, ϕ, η)) Here and in the following, we denote by a tilde suitably rescaled field components. Entering with these equations the single inclusive hadron spectrum (18) and performing the integrals over space rapidity η and azimuthal angle ϕ using the identities in appendix A, one finds Here, I n (z) and K n (z) are modified Bessel functions of the first and second kind, respectively, and the integration kernels Z 1 , Z 2 , Z 3 and Z 4 are given in the appendix A.  For a given form of the independent hydrodynamic fields T 0 , u r 0 ,π t 0 andπ ηη 0 , eq. (24) determines the corresponding background particle spectrum as an integral over the freeze-out hyper surface at fixed background temperature. It can be evaluated for different particle species by using the appropriate relation between m T and p T and adding the appropriate degeneracy factors for spin and isospin. In Fig. 2 we show the functions τ (α), r(α), u r 0 (α),π t 0 (α) and π ηη 0 (α) for a hydrodynamic calculation corresponding to central Pb-Pb collisions at the LHC. The parameters such as initialization time, initial temperature etc. have been chosen as in Ref. [36] and the freeze-out curve corresponds to Fig. 1.
In Fig. 3 we plot the background contribution S 0 (m T , p T ) to the single inclusive transverse momentum spectrum of pions and protons, neglecting spin and iso-spin degeneracy factors as well as chemical potentials. We compare different values of the shear viscosity to entropy ratio. We also show the result obtained without the term linear in the shear stress, i. e. withπ t 0 =π ηη 0 = 0 in eq. (24). The result changes only slightly: At LHC energies the lifetime of the fireball is long enough for the background part of the shear stress to relax largely to its equilibrium value.

B. Particle spectra and harmonic flow due to fluctuating fields
We discuss now how fluctuations around the background hydrodynamic fields contribute to the particle spectrum at freeze-out. To that end we first introduce a specific parametrization for fluctuations in all scalar, vector and tensor fields. For the enthalpy density w = + p and the fluid velocity on the freeze-out surface, we write to lowest order in the deviations from the background fields Here, the fluctuating part is parametrized such thatw,ũ r ,ũ ϕ andũ η are dimensionless. It is useful to introduce a circular polarization basis for the transverse components of the vector field, For the shear stress tensor, we parametrize fluctuations around the background field to linear order in the form Not all components of this shear stress tensor are independent. There is a constraint from the symmetry π µν = π νµ that impliesπ µν =π νµ . Also, the shear stress tensor is traceless, π µ µ = 0 and henceπ τ τ =π rr +π ϕϕ +π ηη . In addition, there is the orthogonality relation u µ π µν = 0 that implies to linear order in the fluctuation field (u µ − (u 0 ) µ )π µν 0 + (u 0 ) µ (π µν − π µν 0 ) = 0. As a consequence of these constraints, only five components ofπ µν are independent. We choose them to beπ rϕ ,π ϕϕ ,π rη ,π ϕη andπ ηη . Again it is useful to employ a circular polarization basis writing The dependent components of the shear stress tensor can be expressed in terms of the independent ones as well as the velocity field, In our construction the background fields w 0 , u µ 0 , π µν 0 and the freeze-out surface are symmetric with respect to azimuthal rotations and Bjorken boost transformations. It is therefore useful to make a Fourier expansion for the fluctuation field according tõ and similar forũ − ,ũ + andũ η and the fluctuationsπ µν in the shear viscous tensor. On the linear level, one can determine the contribution of each mode to the spectrum separately. For symmetry reasons, this contribution will be proportional to a phase factor e imφ+ikηη P , where φ and k η denote the momentum space azimuthal angle and rapidity, respectively. We introduce the short hand notatioñ etc. On the freeze-out hyper surface, each mode is then characterized by the set of ten functions which arew(α),ũ − (α), u + (α) andũ η (α) and the five independent components ofπ µν (α). From equation (18), we find for the contribution of a each single fluctuating mode to the spectrum up to linear order The contributions from fluctuations in the shear viscous tensor are not written out explicitly in this intermediary result (32), but we present them for the final result (33) below. To arrive at this final result, we integrate analytically over the azimuthal angle ϕ and rapidity η. The final result takes then the form Here and in what follows, we compose the background fields into the structure h BG , with ten independent field components Analogously, we denote fluctuations around these background fields by the shorthandh, h = w,ũ r ,ũ φ ,ũ η ,π bulk ,π −− ,π ++ ,π −η ,π +η ,π ηη , with i-th componenth i . As seen from the structure of equation (33), the integral kernels Y i,a and Y i,b specify the response of the single inclusive hadron spectrum to fluctuations on the freeze-out surface in enthalpy, fluid velocity, bulk viscous pressure and shear stress tensor. These kernels are given explicitly in appendix A. They depend in general on the value h BG (α) of the background field along the freeze-out hyper surface. They also depend on the summation index j as specified in the appendix A.
In summary, the output of any fluid dynamical simulation can be written in terms of a set of azimuthally symmetric background fields h BG (α) along the freeze-out hyper surface, supplemented by a set of fluctuationsh(α). Equation (33) specifies the sensitivity of the single inclusive hadron spectrum up to linear order in all fluctuations. It can be used to describe the particle distribution for a particular event. Also, event averages of (33) can be determined if the the probability distribution of the fluctuationsh for an event class is known.

IV. MODE-BY-MODE DECOMPOSITION
Mode-by-mode hydrodynamics [23] provides a more general strategy for characterizing single fluctuating events, but also event averages, correlations and probability distributions. It is based on decomposing arbitrary fluctuating initial conditions (in all scalar, vector and tensor fieldsh (i) ) in a complete orthonormal set of basis functions. For each of these basis functions ("modes"), the (linearized) hydrodynamic equations can then be solved. Knowing the decomposition of initial conditions in terms of orthonormal modes and the dynamical propagation of each mode, one has then a direct bridge between initial conditions and particle spectra for single events and event averages.
Here, we illustrate this strategy first for cases for which the initial conditions do not show fluctuations in the fluid velocity and shear stress, so that the contributionw to the enthalpy density is the only origin of fluctuations. This case is most often assumed in the phenomenological literature. For simplicity we also neglect first all dependence on rapidity η so that the right hand side of (30) has support for k η = 0 only. Under these assumptions, we discuss the effect of fluctuations on single inclusive hadron spectra in section IV A, and on two-particle correlation functions in section IV B. In section IV C, we discuss then shortly how these assumptions can be relaxed.

A. Single inclusive hadron spectra
Following Ref. [35] we decompose the initial fluctuations in an orthonormal basis, where f The complete particle spectrum for the event in (36) is then given to first order in the fluctuations by the linear superposition of modes with wave numbers (m, l) on top of the background spectrum S 0 (m T , p T ), We note that in this expression, events are characterized by the set of weightsw    in (36). For the phenomenologically relevant case of determining the particle spectrum as an average over many events, one can write then up to linear order in the fluctuations where the measure Dw denotes integration over all modesw  . In this case, the particle spectrum (39) corresponds therefore to the background contribution S 0 plus a possible contributions proportional to w (0) l . When the background is taken to be the event average of the fluid fields, the expectation values w (0) l vanishes as well. Instead of expanding the particle spectrum in the small fluctuations as above it may actually be advantageous to expand its logarithm. The underlying observation is that the contribution of a particular fluid cell to the particle spectrum at freeze-out depends on the hydrodynamic fields in an essentially exponential way. We write where the response of the logarithm of the spectrum to a single basis mode is now given by the dynamical map of f Due to the non-linear structure, the event-averaged spectrum is now in general affected by fluctuations, For the particular case of a Gaussian probability distribution p [w] as discussed in Ref. [35], one can perform the Ratio of the full one-particle spectrum S(p T ) to the background contribution S 0 (p T ). The deviation from 1 is due to event-by-event fluctuations (see (43)). The curves are for pion spectra calculated from ideal hydrodynamics, and viscous hydrodynamics with η/s = 0.08 and η/s = 0.3, respectively. The color coding is as in Fig. 1. The dashed lines correspond to the result obtained when the shear viscous contributions at freeze-out ("quadratic ansatz" in (18) ) are neglected.
averaging in (42) and one finds Taking into account that w (m) l ∝ δ m0 for an azimuthally symmetric background, one sees easily that the spectra (43) and (39) agree up to linear order in the fluctuationsw (m) l . To quadratic order inw (m) l , neither (43) nor (39) is complete. However, since linear fluctuations in the fluid fields enter essentially in an exponential way, we expect that (43) resums relevant contributions to higher order in the fluctuating fields. Therefore, we typically calculate the single inclusive hadron spectrum from (43).
To quantify the effect of event-by-event fluctuations on the one-particle spectrum, we plot in Fig. 4 the ratio S(m T , p T )/S 0 (m T , p T ) as defined by (43) for an ensemble of central collisions with initial conditions taken from a Glauber Monte-Carlo model as discussed in ref. [23]. The deviation from 1 is relatively small, or, in other words, the one-particle spectrum is given to good approximation by the corresponding background part with only a small contribution coming from fluctuation effects. This figure 4 can also be viewed as indicating that the difference between the spectrum (39) (that reduces to S 0 (m T , p T ) for central collisions with w (m) l = 0) and the spectrum S(m T , p T ) defined in (43) is relatively small. We believe that this observation holds also for other models of initial state fluctuations.

B. Two-particle correlations and flow coefficients
Let us now discuss the two-particle spectrum. We use the standard definition for the ratio of the two particle spectrum (with both particles from the same event) to the product of two independently averaged one-particle spectra, C 2 (m T 1 , p T 1 , φ 1 ; m T 2 , p T 2 , φ 2 ) = dN pairs p T 1 dp T 1 dφ1dη P p T 2 dp T 2 dφ2dη P event average dN p T 1 dp T 1 dφ1dη P event average dN p T 2 dp T 2 dφ2dη P event average . (44) Using an expansion similar to (42) for both the numerator and denominator one finds For a Gaussian probability distribution p[w], one can actually perform the functional integrals and one finds (46) The double differential harmonic flow coefficients v m {2}(m T 1 , p T 1 ; m T 2 , p T 2 ) are now defined by the expansion We note that the experimental determination of C 2 is usually made such that v 2 0 {2}(m T 1 , p T 1 ; m T 2 , p T 2 ) or a p Tintegrated variant thereof, is normalized to 1. In contrast, in our definition v 2 0 {2} does not necessarily equal 1 even when integrated over p T . This reflects a dispersion in event-by-event values of multiplicities, more specific the p Tintegrated v 2 0 {2} equals the ratio of N 2 and N 2 . We could adapt our normalization of C 2 to the experimental one by dropping the factor in the first line of (46). Since the argument of the exponential is typically small this leads to a small correction only.
In general, for very large fluctuations, the relation between (46) and (47) is non-trivial and the harmonic flow coefficients have to be calculated numerically by doing the appropriate Fourier transform. For not too large values of v 2 m {2} one can expand the exponential in (46) which yields for m = 0 and for m ≥ 0, respectively. The p T -integrated harmonic flow coefficients can be obtained by integrating the expression (49) over p T 1 and p T 2 , appropriately weighted by the one-particle spectrum. More specific, (50) The single-differential harmonic flow coefficients can be defined as and by construction they are normalized such that We note that in general, the flow coefficients v 2 m {2}(m T 1 , p T 1 ; m T 2 , p T 2 ) that enter the two-particle spectrum (47) do not factorize into products v m {2}(m T 1 , p T 1 ) v m {2}(m T 2 , p T 2 ). Therefore, testing the breakdown of such a factorization experimentally as in Ref. [37] should not be interpreted in general as indicating the limited validity of a fluid dynamic description.
We note that Eqns. (46), (48) and (49) constitute rather simple relations between the initial correlations in the expansion coefficients w (m) l and the observable two-particle spectrum. The relation is given by the functions θ (m) l which can be calculated and tabulated in an appropriate way.

C. Generalizations
So far, we have concentrated our discussion on initial density fluctuations and vanishing rapidity dependence (Bjorken boost invariance). In general, one might expect that a detailed model of the initial state and early nonequilibrium dynamics predicts also the size and shape of initial fluctuations in other hydrodynamic fields, in particular fluid velocity and shear. Moreover, these fluctuations may have non-trivial rapidity dependence. We now generalize the most important expressions for the particle spectrum and harmonic flow coefficients to this case.
We chose the hydrodynamical background fields (34) such that the event-average over the fluctuating fields (35) vanish, h i = 0. The expansion of the hydrodynamical fields at the initial time τ 0 , Eq. (36) can be generalized to where f (m) i,l (r) are appropriate basis functions, see for example [35]. Equation (53) describes the most general fluctuation around the azimuthal rotation and Bjorken boost symmetric background in the hydrodynamical fields at the initial time τ 0 .
The contribution of each individual mode to the particle spectra at freeze-out can be determined in complete analogy to the density modes discussed above. The functions θ (54) Equation (43) which gives the one-particle spectrum for a Gaussian probability distribution of initial fluctuations gets generalized to S(m T , p T ) = dN p T dp T dφdη P event average Note that this does not depend on φ or η p in agreement with the assumed statistical symmetries. In a similar way, Eq. (46) for the two-particle correlation function generalizes for h i = 0 to

V. CONCLUSIONS
We have discussed here the kinetic freeze-out for heavy ion collisions using a background-fluctuation splitting for the hydrodynamical fields. To this end, we have introduced in section II a version of the Cooper-Frye freeze-out prescription according to which particles decouple instantaneously at constant background temperature rather than at constant temperature. We have then discussed in section III how fluctuations in all fluid dynamic fields can be treated as perturbations on this freeze-out hyper surface. In general, one can adopt a setting in which the background solution is invariant with respect to azimuthal rotations and longitudinal boosts. The background contribution of the particle spectrum has these symmetries, as well. Deviations from these symmetries come then from event-byevent fluctuations in the hydrodynamical fields. In section IV, we have decomposed arbitrary fluctuations in the initial conditions into an orthonormal set of basis functions. To lowest (linear) order, these basis functions can be propagated independently in hydrodynamics, and they can then be hadronized independently at freeze-out. The central results of the present paper are (43) and (46) that describe the single inclusive hadron spectrum and two-hadron correlations functions in terms of the hadronic response θ (m) l (m T , p T ) of basis functions at freeze-out and the weights that these basis functions carry in the initial conditions. After having calculated the hadronic response θ (m) l (m T , p T ) for the basis functions only one time, these equations allow one to determine the spectra and two-particle correlations for arbitrary events without further fluid dynamic simulations.
The formalism presented here is accurate to lowest (linear) order in fluctuating fields only, and we did not attempt to go to quadratic and higher orders in the fluctuating hydrodynamical fields in this paper. We plan to study the numerical relevance of non-linear corrections in subsequent works, addressing the question of non-linear corrections in the fluid-dynamic evolution of basis functions, and the question of additional non-linearities at freeze-out. Within the present paper, we note only that the Cooper-Frye freeze-out formula (2) contains actually contributions from all orders in the velocity and temperature fields, since the hydrodynamic fields enter in the exponent of occupation numbers as in Eq. (17). Therefore, although complete only to leading order in fluctuations, the formalism presented here accounts at least for some of the expected non-linear corrections. We also note that the azimuthal and boost invariance of the background that has greatly simplified our discussion to linear order will have important implications for higher order contributions, as well. For example, for a single event, the particle distribution may contain a mode ∼ e imφ which at quadratic order has contributions from modes e im1ϕ and e im2ϕ in the hydrodynamic fields but these modes must fulfill m 1 + m 2 = m. The situation is similar with respect to Bjorken boost invariance: a mode e ikηη P may have contributions at quadratic order from two modes e ikη1 and e ikη2 but there is the constraint k η = k η1 + k η2 .
Finally we emphasize that our current study is based on a sharp kinetic freeze-out and includes no phase of hadronic scatterings and resonance decays between the hydrodynamic regime and free streaming. In principle this can be incorporated into the mode-by-mode formalism by solving the corresponding kinetic equations for the background contribution and, in linearized form, for the fluctuations. In particular the influence of resonance decays can be quite sizable for certain observables, most prominent particle identified spectra and harmonic flow coefficients. The formula presented in the current paper could be used to initialize this hadronic scattering and decay phase although a description based on occupation numbers for interacting particles would be desirable.

Appendix A: Rapidity and azimuthal integrals
In this appendix, we provide details on how to calculate from (18) the expressions for the spectrum of the background field (24) and the contributions (33) to first order in the fluctuating fields. For the background contribution, we calculate from equation (23) the shear viscous term at freeze-out, The integrals over space rapidity that appear in the Cooper-Frye freeze-out of different terms of (18) are then of the form where the integrand f * (η) on the right hand side is either unity (we write * = 0) or involves up to three powers in cosh η and sinh η. We denote each such power by a subscript 'c' or 's' respectively, so that for instance the integral R ccc is obtained for the integrand f * (η) = cosh 3 (η). All relevant integrals can then be expressed explicitly in terms of Bessel functions K n (z) of the second kind, We define also the abbreviations to be used below. Similarly one can perform the integrals over the azimuthal angle ϕ. Introducing the shorthand and considering for the integrand g * (ϕ) up to three powers in cos ϕ and sin ϕ such that for instance A ccs is the integral for g * (ϕ) = cos 2 ϕ sin ϕ, we find A ccc (m, z) = 1 8 (I m−3 (z) + 3I m−1 (z) + 3I m+1 (z) + I m+3 (z)), A css (m, z) = 1 8 (−I m−3 (z) + I m−1 (z) + I m+1 (z) − I m+3 (z)) .
With the help of the analytic expressions compiled in this appendix one can reduce the three-dimensional integration over the freeze-out surface to a one-dimensional integral along a curve in the τ -r-plane that can be easily done numerically.
Using this in eq. (B2) shows that the difference between the freeze-out at constant temperature and the freeze-out at constant background temperature indeed vanishes to linear order in the fluctuating hydrodynamical fields. More generally, to higher orders in fluctuating fields, one can still write the particle spectrum as an integral over the freeze-out curve determined for the background solution although there might be small correction terms accounting for the fact that freeze-out happens at fixed fluctuating temperature. On the other side, in the absence of a precise theory where and how freeze-out really happens, we do not see a physics argument to prefer one of the prescriptions. Therefore, one may simply work with the freeze-out at strictly constant background temperature to all orders in fluctuations. This should be fine also if the hydrodynamic description is supplemented by an additional phase of hadronic scatterings and decays described by kinetic theory.