How (non-) linear is the hydrodynamics of heavy ion collisions?

We provide evidence from full numerical solutions that the hydrodynamical evolution of initial density fluctuations in heavy ion collisions can be understood order-by-order in a perturbative series in deviations from a smooth and azimuthally symmetric background solution. To leading linear order, modes with different azimuthal wave numbers do not mix. Quadratic and higher order corrections are small and can be understood as overtones with corresponding wave numbers.

In recent years, fluid dynamic simulations of relativistic heavy ion collisions have provided strong evidence for a picture according to which the momentum distributions of soft hadrons result from a fluid dynamic evolution of initial density fluctuations, see Refs. [1,2,3,4] for recent reviews. The research focusses now on understanding in detail the mapping from fluctuations in the initial state to experimentally accessible observables in the final state [5,6,7,8,9,10,11,12,13,14,15,16]. The present letter aims at quantifying to what degree this hydrodynamic mapping is linear in the strength of initial fluctuations around some suitably chosen background, and on what scale non-linearities arise. This is of interest since an approximately linear relation (by which we mean a mapping in which non-linearities can be understood as small corrections of a predominantly linear mapping) would provide a particularly simple and thus particular powerful tool for relating experimental observables to the initial conditions of heavy ion collisions and to those properties of matter that govern their fluid dynamic evolution [17].
We consider initial conditions of heavy ion collisions, specified in terms of fluctuating fluid dynamic fields h i on a hyper surface at fixed initial time τ 0 . Here, the index i runs over all independent fields, h i (τ, r, ϕ, η) = w, u r , u φ , u η , π bulk , π ηη , . . . , including e.g. the enthalpy density h 1 = w, three independent fluid velocity components, the bulk viscous tensor, the independent components of the shear viscous tensor, etc. In the following we assume Bjorken boost invariance and drop the rapidity-argument η in the hydrodynamical fields. Following Refs. [17,18], we express h i in terms of a background component h BG i and an appropriately normalized perturbationh i . The background is taken to be a solution of the non-linear hydrodynamic equations initialized at τ 0 with an azimuthally symmetric average over many events. It is evolved with the fluid dynamic solver ECHO-QGP [22]. For any sample of events, this background needs to be determined only once. The time evolution of theh i is viewed as a perturbative series on top of the background fields,  Upper row: the time dependence of the enthalpy density is shown separately for the background w BG (τ, r) and for the perturbationw (2) (τ, r) initialized withw (2) 1 = 0.5. Middle row: the dependence of the perturbationsw (2) on the initial weight for τ = τ 0 + 5fm/c (left) and scaled by the initial weights,w (2) (τ, r)/w (2) 1 (right). This scaling establishes that the fluid dynamic response to perturbations is approximately linear. Lower row: same results as shown in middle row, but for τ = τ 0 + 10fm/c. where r = ∞ 0 dr r, ϕ = 2π 0 dϕ etc. The kernels G ij , H ijk (and corresponding terms for higher orders inh i ) depend on the time-evolved background h BG i only. Due to the azimuthal rotation symmetry of the background, G ij depends on the angles ϕ and ϕ only via the difference ϕ − ϕ and similarly for H ijk . The question we raise in the title can now be made more precise: We ask whether the expansion (2) is possible for a suitably chosen background 1 and whether it is dominated by the first linear term. To address this question, we compare in the following numerical results from a full causal dissipative hydrodynamic evolution to expectations based on the structure and on the symmetries of the perturbative series (2).
For the initial conditions, we make assumptions that are widely spread in the phenomenological literature. The initial transverse velocity components vanish, the longitudinal velocity is Bjorken boost invariant, the shear stress tensor is initialized by its Navier-Stokes value, and the bulk viscous pressure is neglected. Initial fluctuations reside then only in the initial enthalpy density w(τ, r), that we parametrize in terms of an azimuthally averaged background w BG (τ, r) and the weightsw (m) l of the azimuthal (m) and radial (l) wave numbers of a discrete orthonormal Bessel-Fourier decomposition [17] w(τ 0 , r, ϕ) = w BG (τ 0 , r) 1 + Here k is the l-th zero of the modified Bessel function J m and R = 8 fm throughout this work. Sincew(τ, r, ϕ) is real, we havew (m) (τ, r) =w (−m) * (τ, r). In the following, we take the weights with m ≥ 0 as the independent ones and writẽ The corresponding modes with m < 0 are then not independent and are defined by the condition |w We consider first the case for which one single fluctuating basis mode is embedded on top of w BG (τ 0 , r). For example, we specify this mode with the weightw (2) 1 , so that the initial enthalpy density reads For one single mode, we can set without loss of generality ψ (2) 1 = 0. Assuming for simplicity a Bjorkenboost invariant longitudinal dependence, we evolve these initial conditions with the 2+1 dimensional version of the hydrodynamical code ECHO-QGP [22] with a value η/s = 1/4π for the ratio of shear viscosity to entropy density 2 . Following Ref. [19], we use the equation of state s95p-PCE which combines lattice QCD results at high temperatures with a hadron resonance gas at low temperatures. The background w BG used throughout this paper is initialized at τ 0 = 0.6 fm/c with an azimuthally symmetric average of Glauber model initial conditions for Pb+Pb collisions at the LHC, described in Ref. [18]. The time evolution of w BG determined from ECHO-QGP is shown in Fig. 1. The time-evolved fluctuationw (2) (τ, r) is determined from the full hydrodynamic evolution via Fourier analysis. Results are shown in Fig. 1 for different weights w (2) 1 . Fluctuations at time τ 0 are cut-off in the region of very low background density, see e.g.w (2) (τ 0 , r) in Fig. 1 -we have checked that this does not affect our results. The main conclusion from Fig. 1 is that at all relevant times and even for relatively large initial amplitudesw (2) 1 , the fluid dynamic response w BGw (2) (τ, r) to an initial perturbation scales approximately linearly with the weightw (2) 1 . This is the behavior expected from the linear term in eq. (2). We observe this linear dependence with similar accuracy also for other basic modes (data not shown).
The almost exact linear scaling of the hydrodynamic responsew (m) (r, τ ) with the initial weights w (m) l does not imply that non-linearities are absent. To see that, consider the Fourier seriesh i (τ, r, ϕ) = Since the kernels in (2) depend only on the background field, they are invariant under azimuthal rotation and their Fourier expansions read 1 . Right column: Same results but rescaled by the second (third) power of the weightw (2) 1 . This scaling establishes thatw (0) (τ, r) andw (4) (τ, r) (w (6) (τ, r)) can be understood as overtones that are induced by the initial second harmonic perturbation as a perturbative second (third) order correction to (2). The short-range fluctuations in the rescaledw (6) For the case that initial conditions contain only fluctuations of enthalpy density, we haveh with and similarly for H i11 ; l l . According to (7), if one initializes fluctuations with a single mode of weightw (m) l , as done in Fig. 1, then corrections that are quadratic in the fluctuationsh i (τ 0 ) will not appear in the time-evolved harmonicsh (2) shown in Fig. 1 with the leading quadratic (w BGw (0) , w BGw (4) ) and cubic (w BGw (6) ) corrections displayed in Fig. 2. We observe that quadratic (cubic) corrections scale with the square (the cube) of the initial weight w 1 , as expected from (7). Moreover, even for a weightw So far, we have demonstrated with examples that Eq. (8) explains the dominance of linear response and the relative size and ordering of the overtones induced by one basis fluctuation. We have checked extensively that the same equation explains also the structure and symmetries of the hydrodynamic interactions between initial perturbations with different wave numbers. Fig. 3 illustrates this point with a case for which two perturbationsw are embedded on top of the initial background fields. We have checked that the second (third) harmonicsw (2) r)) of the fluid dynamic response scale linearly with the initial weightw 1 ) and that they agree to high accuracy with the response to an initial configuration in which only one modew    1 , as seen in Fig. 3. Also, according to (8), the phases of the first and fifth harmonics are determined by the orientations of the initial perturbations. The comparison with the full numerical results in the middle panel of Fig. 3 shows that this perturbative expectation is realized approximately (strong deviations are seen only for values of the radius r for which either Re w (m) or Im w (m) approach zero and for which the orientation is thus not well defined). To add one level of complication, we consider finally the sixth harmonics w BGw (6) that, according to (8), receives corrections of second order inw (3) 1 and of third order inw (2) 2 . Fig. 3 shows that weighting both contributions with the perturbatively expected information on phases and amplitude provides for a full quantitative understanding of the numerically determined signal w BGw (6) as overtones of the two initial perturbations. Realistic initial conditions for the fluid dynamic evolution of heavy ion collisions are expected to involve fluctuations on many different length scales and with large amplitude. Could it be that the examples discussed so far, although initialized with relatively large amplitudes, are still academic, and that they cannot be extended to deal with the complexity of a realistic heavy ion collision? To lay such concerns at rest, we have embedded simple basis modes in realistic initial conditions with many and large fluctuations. Fig. 4 shows such an initial distribution. It is constructed by subtracting from an arbitrary initial condition generated by a Glauber model the contribution leading to a second harmonic and adding then the perturbation of (5). In this way, we have an analytically controlled initial perturbation on top of a realistically fluctuating background, and we can extract this initial perturbation and the time dependence of its fluid dynamic response via Fourier analysis. The lower panel of Fig. 4 shows that this dynamical response in an event with realistic fluctuations is described to high accuracy by the linear response on top of the smooth background that we had established in Fig. 1.
In summary, the evolution of initial anisotropic density perturbations as determined numerically with the fluid dynamics solver ECHO-QGP seems to follow a simple pattern that can be understood order-by-order in a perturbative expansion for small deviations from an azimuthally symmetric event-averaged background. The leading order is linear and modes with different azimuthal wave numbers do not mix. Quadratic and higher orders can be seen as next-to-leading order corrections. They influence modes with azimuthal wave numbers that can be written as sums (or differences) of the seed wave numbers. Since the non-linear couplings are numerically small, the higher harmonics generated by two-mode or three-mode interactions will often be small in comparison to initially present and linearly evolving perturbations. This ordering may be less pronounced for non-central collisions where the elliptic modes (with m = 2) have a particularly strong weight such that its overtones with m = 4, 6 etc. may dominate over primordial density perturbations  with these wave numbers. We also note that the relative importance of linear and non-linear terms depends significantly on the dissipative properties of the medium. In exploratory studies that lie outside the scope of this letter, we found that for increasing η/s, the linear response is more dominant and the relative weight of higher non-linear orders decreases. Such a more laminar behavior is expected on general grounds. All results shown in this Letter are for a rather minimal value of η/s = 1/4π which maximizes contributions from non-linear corrections.
The perturbative picture expressed by eqs. (2) or (7) provides a useful ordering scheme for a more detailed understanding of the evolution of fluctuations in hydrodynamic fields. The numerical results discussed here provide motivation for a more formal and thorough development of this kind of perturbation theory. It is unclear whether all initial perturbations in hydrodynamic fields follow a similar pattern. In particular, from general fluid dynamic considerations one may expect that non-linear terms are more important for vorticity excitations [24].
Another open question concerns the freeze-out from hydrodynamic fields to particle spectra. To linear order in a mode-by-mode treatment of fluid dynamic perturbations this was discussed recently [21]. However, non-linear corrections arise there as well. For instance, the proportionality v 4 ∝ (v 2 ) 2 was conjectured in [23] as a consequence of this effect. It would be interesting to see how large these non-linearities are in comparison to the ones from the fluid dynamic propagation as discussed here.
It will also be interesting to investigate whether the approach presented here can provide a more detailed dynamical underpinning of results in the recent literature. For instance, the experimentally accessible reaction plane correlations formulated and studied in Refs. [7,25,26,27] are akin of the constraints on the phases of higher order contributions in the perturbative ansatz (8). Similarly, the (non-linear) response to cumulants of the initial density distribution as formulated and studied in Refs. [10,11,27] may be related to some integrated version of eq. (7).