Nonlinear evolution of density and flow perturbations on a Bjorken background

Density perturbations and their dynamic evolution from early to late times can be used for an improved understanding of interesting physical phenomena both in cosmology and in the context of heavy-ion collisions. We discuss the spectrum and bispectrum of these perturbations around a longitudinally expanding fireball after a heavy-ion collision. The time-evolution equations couple the spectrum and bispectrum to each other, as well as to higher-order correlation functions through nonlinear terms. A non-trivial bispectrum is thus always generated, even if absent initially. For initial conditions corresponding to a model of independent sources, we discuss the linear and nonlinear evolution is detail. We show that, if the initial conditions are sufficiently smooth for fluid dynamics to be applicable, the nonlinear effects are relatively small.


Introduction
It has been noted repeatedly [1][2][3] that the collective dynamics of the hot and dense QCD matter produced in ultra-relativistic heavy-ion collisions shares many commonalities with the expansion dynamics of the Early Universe. Both systems are initially relativistic, of high energy density and rapidly expanding. Relativistic fluid dynamics describes at least some stages of their evolution, while particle species freeze-out from this fluid at characteristic densities. With the progress in measuring and analyzing fluctuations in the cosmic microwave background (CMB) [4,5], and with the discovery that flow measurements in heavy-ion collisions result from the propagation of "initial" event-by-event fluctuations [6][7][8], it is by now clear that the little bangs studied in heavy-ion collisions and the Big Bang studied in cosmology correspond to the smallest and largest physical systems, respectively, for which fluctuation analyses can provide information about material properties. Here, we recall these commonalities in order to raise the prospect that progress in either field could be made by transferring techniques and insights from the other. There are, of course, basic features with respect to which the two systems in question are quite distinct: for instance, both systems are characterized by vastly different scales and their dynamics is governed by different fundamental forces. On the other hand, there is an underlying unifying principle: the applicability of a relativistic fluid description.
In heavy-ion collisions, the nonlinearities in the evolution are expected to be maximal initially, but they decrease with time due to dissipative effects. In cosmology, the initial density perturbations are very small and linear perturbation theory is a good description at early times. But at late times, the gravitational instability leads to strong growth, resulting in the emergence of large-scale structure (LSS) and the formation of galaxies and galaxy clusters. It is apparent then that the formulations of cosmological perturbation theory most relevant for our purposes are the ones applied to the problem of LSS.
Even though the focus, both in heavy-ion physics and in LSS cosmology, is on propagating perturbations in fluid dynamics, there are marked differences in the way the dynamics of these perturbations is implemented. In heavy-ion physics, one usually initializes the fluid evolution with the help of Glauber Monte Carlo models or their variants based on saturation physics [9][10][11][12], which specify the fluid dynamic fields including their fluctuations at the initial time. Typically, these fields are then propagated numerically through hydro-solvers [13][14][15][16] up to the final state, in an effort to map out the relation between the (in principle wide) model space of conceivable initial conditions and the experimentally accessible data that can constrain them. Only very recently the first efforts were made to abstract from this explicit propagation of many model-dependent initial conditions, for example by characterizing the statistics of initial perturbations in terms of n-point correlation functions [17][18][19]. In cosmology on the other hand, the fluid dynamic equations of motion usually are not solved for individual fluctuating fields. Instead, the objects of interest are correlation functions, such as the two-point correlation function of fields (the spectrum), the three-point (the bispectrum) as well as higher-order correlation functions. The various analytical methods that have been developed in order to go beyond linear perturbation theory rely on resummations of subsets of perturbative diagrams [20][21][22][23][24][25][26][27][28][29][30][31][32] or the use of effective field theory [33][34][35].
A transparent and efficient approach to the problem of LSS in cosmology is the timerenormalization-group (time-RG) formulation [27]. Within this approach, the physical system is described in terms of a hierarchy of coupled evolution equations for the spectrum, bispectrum and higher spectra. The method has been applied to ΛCDM and quintessence cosmologies [27], allowing for a possible coupling of dark energy to dark matter [36], or a variable equation of state [37]. It has also been used for the study of models with neutrinos [38]. The main aim of the present paper is to apply this formulation of relativistic fluid dynamics commonly employed in cosmology to a simplified test system often studied in heavy-ion physics. More specifically, we shall apply the time-RG to the problem of time evolution of two-and three-point correlation functions of density and flow fluctuations on top of a Bjorken background [39].
It must be noted that the issue of initial conditions is more complex in heavy-ion collisions: Immediately after the collision the local deviations from a smooth background are expected to be large, at least according to many current phenomenological models, in which the initial energy distribution is dominated by the "lumpy" substructure of nuclei in terms of their constituents (nucleons) and the spatial distribution of nucleon scattering centers. On the other hand, it is clear that fluid dynamics strictly applies only after a phase of early non-equilibrium evolution during which dissipative effects are very large and smoothen the spatial structures, at least to some extent. An interesting question is whether the fluid dynamics of heavy-ion collisions is essentially linear or strongly nonlinear at the time and length scales of a typical event and for those regimes of momenta (or wavevectors) for which the gradient expansion that underlies fluid dynamics is applicable.
This question may also depend on the kind of perturbations considered. For vorticity fluctuations, two of us have argued some time ago that the time and length scales of a heavy-ion collision might be sufficient in order to see the onset of nonlinear or turbulent behavior [41]. On the other hand, within the current phenomenological picture density or compressional modes (sound waves) are much more important. In the present paper we concentrate entirely on the latter and assume that the fluid velocity is vorticity free.
This paper is organized as follows: we introduce the Bjorken model and the treatment of fluctuations in section 2. In section 3, we derive the corresponding evolution equations for two-point and three-point correlation functions of fluid dynamic fields and we discuss how they can be solved perturbatively. Section 4 introduces a model for the initial conditions. The UV sensitivity of fluid dynamics determines its range of applicability and is discussed in detail in section 5 before we turn to the discussion of numerical results in section 6. The main results of our study are then summarized in the conclusions. All dimensionful parameters are measured in femtometers (fm) when units are not explicitly indicated.

Fluid dynamics of perturbations in the Bjorken model
We start from the relativistic hydrodynamic equations of motion for a fluid without any conserved charges. Keeping the effects of shear viscosity (η), but not of bulk viscosity, within the first order in a gradient expansion, the equations of motion read Here, and p denote the local energy density and pressure, respectively. We work with metrics of signature (−, +, +, +). In terms of the covariant derivative ∇ µ and the projector ∆ µν = g µν + u µ u ν on the subspace orthogonal to the flow vector u µ , the traceless and transverse (i.e. orthogonal to u µ ) tensor σ µν takes the form 3) The solution of eq. (2.2) in full generality is a formidable (numerical) task. Instead, we shall employ in the following a split between background and fluctuations in order to discuss a class of interesting solutions. More precisely, we start from an analytically solvable smooth background solution of eq. (2.2) that was first formulated by Bjorken [39] and captures, despite its simplicity, phenomenologically relevant features. We then discuss the propagation of small perturbations on top of this background solution. This approach is analogous to the one employed in cosmology, where fluctuations propagate on top of the smooth Hubble flow of a Friedmann-Robertson-Walker (FRW) Universe.
Bjorken's model [39] is based on the idea that in nuclear collisions at ultra-relativistic energy particle production is almost flat in rapidity. This motivates the assumption that all fluid dynamic fields are initially independent of space-time rapidity, y = arctanh(x 3 /|x 0 |). Fluid dynamic evolution preserves this symmetry. Using the coordinates (τ, x 1 , x 2 , y), with proper time τ = x 2 0 − x 2 3 and metric g µν = diag(−1, 1, 1, τ 2 ), the solutions are then y-independent and the problem reduces to a (2+1)-dimensional one. A strong further simplification is obtained if one also assumes translational and rotational invariance of the initial conditions in the transverse plane. This yields a one-dimensional toy model, the so-called Bjorken model [39]. In terms of the enthalpy w = + p = sT and the kinematic viscosity ν = η/w, the evolution equation of this one-dimensional model reads For ν/τ 1 and an ideal equation of state = 3p, we obtain the characteristic time dependencies of the Bjorken background for the energy density 5) or, equivalently, the characteristic time-dependence of the Bjorken background temperature, T Bj (τ ) = T Bj (τ 0 ) (τ 0 /τ ) 1/3 . For a time-independent normalized viscosity η/s, the ratio ν Bj (τ )/τ in eq. (2.4) decreases like (2.6) Therefore, replacing the bracket in eq. (2.4) by unity is an approximation that is consistent with the late time behavior. This approximation is also justified at early times (τ 0 1 fm, T Bj (τ 0 ) 1/fm) for small values of the viscosity to entropy ratio, such as those favored by holography [40]: η/s 1/(4π).

Perturbations on a Bjorken background
The propagation of perturbations on top of the Bjorken background solution has been discussed in ref. [41]. Here, we recall the results relevant for the present work, and the approximations on which they are based, but we refer to ref. [41] for a detailed derivation. The background solution is defined by (2.5) and by the rapidity-independent Bjorken flow field u µ Bj = (1, 0, 0, 0). We use an ideal equation of state, = 3 p and allow for small local fluctuations around this solution. For the flow field, local fluctuations are parameterized by the transverse and rapidity components u 1 , u 2 , u y , respectively. The normalization condition u µ u µ = −1 of the local fluid velocity u µ = (u τ , u 1 , u 2 , u y ) implies then (u τ ) 2 = 1 + (u 1 ) 2 + (u 2 ) 2 + τ 2 (u y ) 2 = 1 + u j u j . (2.7) Here and in the following, we use for latin indices a summation convention that runs over the spatial components 1, 2, y, with the spatial part of the metric being g ij = diag(1, 1, τ 2 ). We consider small local perturbations in the sense that u j u j (x) 1, and we neglect terms that are parametrically suppressed due to u j u j (x) 1 or ν/τ 1 compared to other terms with the same combination of derivatives acting on the fields. For every combination of derivatives there is one term of lowest order which is not neglected. It is then straightforward to derive from eqs. (2.1) and (2.2) equations of motion for the fluctuations u 1 , u 2 , u y in flow velocity, and for fluctuations δ in the energy density (see ref. [41]). These equations of motion take a particularly transparent form if they are written in terms of a rescaled time variable 8) and in terms of the logarithm of the perturbed energy density or, equivalently, temperature: Here and in the following, we assume that the perturbation δ is small compared to the energy density Bj of the background. After the rescaling of the fluid dynamic fields, the equation of motion for the rescaled logarithm d τ of the energy density takes the form 11) and the equations of motion for the rescaled velocity components (j = 1, 2, y) read (2.12) In the following, we denote the velocity divergence, sometimes also referred to as expansion scalar, by The two equations (2.12), (2.13) are the starting point of the present work. For the following, we subject them to the following additional approximations: • We neglect fluctuations in the longitudinal y-direction.
Remarkably, in the equations of motion (2.11), (2.12), contributions from the ydependence of fluctuations are suppressed by powers of 1/τ . At late times, the equations of motion thus become effectively (2+1)-dimensional. In the following, we make the assumption that the fluctuating initial conditions are independent of y. As a consequence, we deal with a (2+1)-dimensional system at all times, and we consider equations of motion for v j = v j (t, x 1 , x 2 ), with j = 1, 2, but not with j = y.
• We assume that the fluid is irrotational.
In a perturbative treatment, fluctuations in vorticity couple only nonlinearly to density perturbations. Since a perturbative treatment of fluctuations is expected to be valid for the short time scale ( < ∼ 10 fm) relevant for heavy-ion collisions, we neglect in the following vorticity and thus the effects of fluid turbulence. For an irrotational fluid, the velocity can then be written as the gradient of a scalar: v m = ∂ m f , with m = 1, 2. Taking the divergence of eq. (2.12), one can write an evolution equation for the scalar f , A dot denotes a derivative with respect to t, while repeated indices are summed. We have not used upper indices in order to avoid confusion with the more general, rapidity-dependent case discussed earlier.
• An IR cutoff can retain important aspect of a finite transverse extension. The treatment of fluctuations on top of the Bjorken background is based on assuming their statistical homogeneity and isotropy in the transverse plane. This idealization can be meaningful only for distances smaller than the typical transverse extension of a heavy-ion collision, which we take to be around 8 fm. It implies that the power spectrum depends only on the magnitude of the momentum vector, while the bispectrum depends on the magnitudes of three momentum vectors. It is a severe approximation, since it prevents for instance the development of a transverse background flow, which would appear inevitably for systems of finite transverse extension. Despite this limitation, an important effect of the finiteness of the system can be retained in our treatment. Namely, the initial power spectrum is expected to vanish below the momentum scale that corresponds to distances larger than ∼ 8 fm, as no correlations within a beam are expected at such length scales. In principle, this can be implemented by imposing an effective IR cut-off on the initial fluctuations. In practice, the results presented in the following are insensitive to this kind of IR cut-off.

Correlations and spectra
In the following we discuss fluctuations on top of a Bjorken background in a set-up akin to recent treatments of cosmological fluctuations. The role played by Friedmann's equation for the scale factor in the cosmological model will be played by the evolution equation for the Bjorken-expanding background fluid derived in section 2. In order to discuss fluctuations around the Bjorken solution similarly to cosmology, we first Fourier transform their evolution equations and derive the corresponding equations for correlations and spectra.

Nonlinear evolution equations in Fourier space
As summarized in section 2, equations (2.11), (2.12) are our starting point for an analysis of the nonlinear evolution of transverse fluctuations induced by mode coupling and the leading dissipative effects of viscosity. Here, we consider the case of an irrotational fluid where v m = ∂ m f for m = 1, 2, v y = 0, and d τ , f functions of t, x 1 , x 2 . The primary dynamic fields d τ and f can then be written in Fourier space as They satisfy the equations of motioṅ We define the combined field φ a (t, k), a = 1, 2 with and the matrix (3.8) where we have not indicated explicitly the dependence of φ on time. The nonzero elements of γ abc are Note that eq. (3.9) implies that one can always define γ abc (p, q) such that The invariance of the time evolution equation In addition, one observes from eqs. (3.10), (3.4), (3.5) and (3.6) that γ abc (p, q) = γ abc (q, p) = γ acb (p, q). (3.14)

Definitions of expectation values and spectra
For symmetry reasons only the zero mode with k = 0 can have an expectation value and we define φ a (k) =φ a δ(k), (3.15) Note that, because of the definition v m = ∂ m f , a spatially constant part of f (t, x) has no physical meaning. One can always change it without changing the velocity in order to achievef k=0 = 0 and thereforeφ 2 = 0. On the other hand, a spatially constant part in d τ (t, x) corresponds through eq. (2.9) to a deviation from the Bjorken background solution of eqs. (2.5), (2.6). This is a physical effect that may arise when energy from the dissipation of perturbations heats up the system. The power spectrum P ab (k) is defined to correspond to the connected part of the two-point function +φ aφbφc δ(k)δ(p)δ(q) , (3.17) and the trispectrum Q abcd (k, p, q, r) by  .18) we have not written out all permutations but indicated in square brackets how many there are of each kind. We shall later neglect the effect of the trispectrum on the evolution of the spectrum and bispectrum. This is equivalent to setting Q abcd (k, p, q, r) = 0 in eq. (3.18) and keeping only the disconnected parts of the four-point correlation function. From the definitions it is clear that P ab (k) = P ba (−k) whereas the bispectrum satisfies and similarly for the trispectrum. The δ-functions in the definitions of the spectra are a consequence of our assumption of (statistical) homogeneity of space. Similarly, statistical isotropy implies that the orientation of the vectors k, p, q does not play any role. In particular we have (3.20)

Evolution equations for expectation values and spectra
The evolution equations for the expectation value, spectrum, bispectrum and higher-order spectra result from the differentiation of eqs. (3.16), (3.17) with respect to t. In principle this leads to an infinite hierarchy of equations, with the time evolution of the n-point correlation function involving the (n+1)-point function. This kind of hierarchy is analogous to the BBGKY hierarchy [44] and appears in related form also at other places in theoretical physics, for example in the context of the functional renormalization group [45]. If one uses the definitions in eqs. (3.15)-(3.18) and the evolution equation (3.9) one finds for the expectation value for the spectrum and for the bispectrum We have used in eqs. (3.22) and (3.23) the modified time evolution matrix Note that the set of equations (3.21)-(3.23) is not closed because the last equation depends on the trispectrum. The essential approximation that we make in the following in order to obtain a closed system is to neglect the effect of the trispectrum on the evolution of the bispectrum.

Perturbative solution
From now on, we assumeΩ ab (k) = Ω ab (k). As explained previously, this is realized explicitly in the model discussed here. The formal solution of the above system is given by (we 1 When inspecting eqs. (3.4)-(3.6) in the limit p → 0, q → 0, one may worry about the fact that the value of β(p, q) in eq. (3.6) depends on how this limit is taken. On the other hand, in eq. (3.21), γ222(p, q) = β(p, q) gives the effect ofφ 2 2 on the time evolution ofφ2. We have argued in subsect. 3.2 that one can always setφ2 = 0 without affecting the physically relevant velocity field. Another possibility to understand this issue is to consider the velocity divergence θ = ∂ivi instead of the field f = φ2 as is discussed in subsect. 5.1.
where g ab (k , t, t ) is the linear propagator, which gives the evolution of the field at the The solutions can be expanded in powers of the interaction vertex γ abc in order to establish the connection with perturbation theory. The lowest order, corresponding to linear theory, is obtained by setting γ abc = 0. The linear spectrum P L ab and bispectrum B L abc are given by the first line of each of the above equations. The O(γ) correction to the bispectrum is obtained by inserting P L ab in place of P ab in the r.h.s. of eq. (3.27). Inserting the bispectrum at this order in eq. (3.26) generates the O(γ) and O(γ 2 ) contributions to the power spectrum. Iterating the procedure generates the higher-order corrections. However, differences with perturbation theory arise at higher orders, because of the approximation Q abcd = 0 that we have made in deriving the evolution equations for the power spectrum and bispectrum.

Linear propagator
In order to obtain an analytical expression for the propagator, which is necessary for the iterative solution of eqs.
Finally, the definition χ k = exp(B k t/2) δ k leads tö This equation has simple analytical solutions when either one of the two terms in the parenthesis dominates. The dominance switches from the first to the second term at a time τ /τ 0 (4ν 2 0 k 2 /3) −3/2 . Approximating the kinematic viscosity by its homogeneous value (2.10), and assuming an equation of state = 3p, we find from eq. (2.6) that .
We expect that τ 0 1 fm, T Bj (τ 0 ) ∼ 200 − 500 MeV, while the relevant momenta satisfy k < ∼ 1 fm −1 . For values of η/s ∼ 1/(4π) the first term in (3.32) dominates for the times of interest, which are t < ∼ 10 fm. In the limit that the second term is neglected (B k = 0), the solution of eq. (3.33) is . (3.37) These expressions indicate oscillatory behavior corresponding to propagating sound waves in the plasma. The exponential suppression at late times describes the attenuation of the high-frequency modes through the effect of viscosity. The propagator can be determined easily by fixing the constants c 1 , c 2 through the initial conditionsd k (t 0 ),f k (t 0 ). The four components of g ab are given by the coefficients ofd k (t 0 ),f k (t 0 ) in the resulting expressions. We do not display the explicit form of the propagator because of its length and complicated structure. We have verified numerically that eqs. (3.36), (3.37) provide a very good approximation for all scales and times relevant for our analysis.
In figs. 1 and 2 we display the 11-and 12-components of the propagator for τ 0 = 1 fm and ν 0 = 0.04 fm. These figures document the appearance of waves and the attenuation of the high-frequency modes.
The explicit analytical expressions derived above allow one to deduce certain quantitative features of the solutions, which are also reflected in the evolution of the spectra.  1, 2 occur at values of order 1 for the combination k t 3/4 that appears in the argument of the Bessel functions in eq. (3.36). Another feature is the attenuation of the high-frequency modes induced by the viscosity. The characteristic time for this effect is set by the exponential in eq. (3.36). As B k ∼ ν 0 k 2 , the high-frequency modes quickly die out within a time ∼ 1/(ν 0 k 2 ). This feature is also apparent in figs. 1, 2.

Initial conditions
Before turning to the discussion of solutions of the equations of motion (3.21), (3.22), (3.23) for the expectation value, spectrum and bispectrum, we introduce in this section the initial conditions from which these equations of motions are evolved.

Modelling the initial fluctuations
The evolution equations are initialized at the time τ 0 at which a fluid description becomes approximately valid. The initial values at τ 0 are a result of initial-state physics and of the early non-equilibrium dynamics that drives (local) thermalization. They are to be determined from a theoretical description of that era. Here, we obtain them from a simple model inspired by a Glauber-type description. As in most of the phenomenological literature, we assume that only the density fluctuates initially and we neglect possible initial flow fluctuations. To model the density d of transverse fluctuations, we assume that the deviations from the Bjorken background are small so that where we have assumed a velocity of sound c 2 S = 1/3. The transverse enthalpy density w is assumed to originate from a number of localized sources at positions that we take to be independently and uniformly distributed in the transverse plane. For a single event, the corresponding density field d is given by [18,19] Here, the index j labels different sources or scattering centers, and the strength of the source is given by a smoothened transverse density f (x − x j ) smeared around x j . In the following, we consider the case that all sources have the same normalized strength, d 2 x f (x) = 1. We note that in a Glauber type description, one may associate each source to a nucleon participanting in inelastic processes, or to a pair of participants. In this work, we keep the number of sources N fixed. But the model could be easily extended to include fluctuations in the number of contributing sources or in the strength (i.e. norm) of each source. In that case, it would account for fluctuations in multiplicity. In agreement with the symmetries of the Bjorken model, we take all sources to be uniformly distributed in a transverse area A which corresponds to a probability distribution p(x) = 1/A. The expectation value of the transverse density is independent of the transverse position. The Bjorken background entering (4.1) can be chosen such that by construction

The initial spectra
To calculate correlation functions, we introduce the partition function [19] Z To simplify this partition function, we have made use of the assumption that the random variables x j are independent and that they are uniformly distributed. We have introduced the Fourier transformsJ(k) andf (x) of the source J(x) and the density f (x) of an individual collision center, respectively. Correlation functions can now be calculated as usual, by taking functional derivatives. For example, we find (4.5) Note that although the second term in the last line seems to be singular for k = 0, it is actually finite, because one has (2π) 2 δ (2) (k)/A → 1. For the particular case of the Gaussian source for whichf (k) = 1 (2π) 2 e − σ 2 k 2 2 , (4.6) one finds the initial spectrum of density fluctuations with (4.8) In the limit of point-like sources (σ → 0) the initial spectrum becomes independent of k. It is clear that this limit is somewhat artificial. In reality one expects that the enthalpy due to each scattering center is smooth and non-singular. This is a consequence of the finite size of the nucleon itself, but also of the early dynamics that drives local thermalization. Soon after the collision the dissipative effects are expected to be large, which leads to a fast smoothening of the enthalpy density. If a fluid dynamic description becomes valid at time τ 0 , one expects that spatial structures on much finer scales than cτ 0 have already disappeared. (For concreteness, one might think of dynamics close to free streaming.) Based on these considerations, one expects that σ ≈ τ 0 is reasonable. Higher-order correlation functions can be calculated in analogy to eq. (4.5). For example, the three-point function is given by In the second line we have written explicitly only one of three permutations of k 1 ,k 2 , k 3 . For nonzero momenta and the Gaussian source we have

UV sensitivity and the range of applicability of fluid dynamics
Fluid dynamics is the long wavelength limit of a consistent transport theory, and thus it needs not be well-behaved in the UV. Understanding the sensitivity of the fluid dynamic evolution on UV cut-offs is therefore important for determining its range of applicability. In the model studied here, the nonlinear couplings γ abc (p, q) that enter the equations of motion (3.21), (3.22), (3.23) arise from the spatial-derivative terms in the fluid dynamic equations. They show power-law growth in the UV. In principle, there are two physical mechanisms that result in UV cutoffs. First, density fluctuations shortly after the collision are expected to be limited to wavenumbers of order |k| ∼ < 1/fm and the initial conditions are therefore UV-tame on physical grounds. Second, dissipative processes are expected to provide for a dynamic mechanism that further damps high-momentum modes. However, despite these two UV-regulating mechanisms, one cannot exclude a relatively strong sensitivity of the solutions on the details of the initial conditions in the UV. Therefore, before turning in the next section to the numerical study of the equations of motion (3.21), (3.22), (3.23), we discuss here the UV behavior of the fluid dynamics for density perturbations, as introduced in section 3.

Derivative couplings and UV properties
We emphasize first that the UV sensitivity of the fluid dynamic evolution cannot be attributed solely to nonlinear terms in the iterative solution of eqs. (3.26)- (3.27). This becomes apparent if one rewrites the equation of motion (3.9) for φ a (t, k) = d k ,f k into an equation for the fluid dynamic fieldsφ a (t, k) = d k ,θ k , where the divergence of the velocity field θ withθ k = −k 2f k replaces the field f . (An analogous formulation is used for the analysis of cosmological perturbations in the formalism that we are employing [27].) The equations of motion forφ a (t, k) are of the same form as (3.9), but with couplings defined byα These couplings are dimensionless. For large momenta, they are much better behaved than (3.4) -(3.6), although there are still divergent directions in momentum space. (For example, the coupling (5.1) diverges for large p and fixed q.) On the other hand, the linearlized evolution is now described by the equationṡ This illustrates that it depends on the choice of independent fluid dynamic fields whether the strong momentum dependence is found in the linear or nonlinear terms of the evolution equation. In eq. (3.9) for φ a , the nonlinear couplings grow like power-laws with increasing momenta. In eq. (5.5), this strong UV-dependence has been reshuffled now mainly into the linear evoluton ofθ k , The same equation (5.5) also illustrates clearly that the reason for this strong momentum dependence is physical. The first term on its r.h.s. indicates that strong flows are induced by density inhomogeneities with large momenta. This is of course a very physical effect: the fluid is accelerated by pressure gradients. The second term results in an exponential suppression of the high-momentum fluctuations through the viscosity. The general evolution involves a competition between these two terms.
It is easy to verify that the use ofd k ,θ k , instead ofd k ,f k , leads to identical results. For our numerical solution we employ the fieldsd k ,f k , mainly because of the symmetry of the effective couplings (3.4)-(3.6) under the exchange p ↔ q.

Mode-mode coupling and the influence of UV modes
To obtain a better understanding of the UV sensitivity of fluid dynamic perturbation theory, we consider here a particularly simple example where some analytic insight is possible, namely the nonlinear corrections to the spectrum P ab (k, t) with k → 0. We concentrate on the O(γ) correction to the density-density component P 11 (0, t) in the presence of an initial bispectrum whose only non-vanishing component is B 111 . The spectrum can be expressed as h(q, t , t 0 ) = t 0 t −4q 3 g 11 (q, t , t 0 ) g 21 (q, t , t 0 ) + 16 9 q 5 ν 0 g 21 (q, t , t 0 ) 2 .

(5.7)
For k = 0 the propagator simplifies to 8) and the couplings are γ 112 (−q, q) = −q 2 /2 and γ 122 (−q, q) = 4q 4 ν 0 /9. In general, nonlinear corrections in the evolution of the spectra arise from the coupling between modes with different momenta. This coupling also results in the transfer of energy between modes. Since much of the energy is initially concentrated in the high-momentum modes, one expects that the nonlinear evolution will induce the enhancement of the lowmomentum ones. We now turn to the numerical solution of (5.6) in order to illustrate these expectations. Numerical results will be presented for t 0 = 3/4 fm (corresponding to τ 0 = 1 fm) and ν 0 = 0.04 fm, and for variations around these default values, where indicated.
In fig. 3, we plot the function h(q, t , t 0 ) of eq. (5.7) that determines the size of nonlinear corrections to the spectrum P 11 (0, t). At early times, h(q, t , t 0 ) is seen to grow very fast with increasing momenta q. At later times, dissipative effects due to the finite shear viscosity dampen these high-momentum modes rapidly. However, despite this viscous damping, the momentum integration in the second term of eq. (5.6) would lead to a divergent result for a constant initial bispectrum B 111 (0, −q, q). To see this, we note first that as a consequence of the exponential suppression factor in eq. (3.36), the propagator g ab (q, t , t 0 ) that defines h(q, t , t 0 ) takes very small values when ν 0 q 2 t > ∼ 1. For large q and constant B 111 (0, −q, q), the leading contribution to the integral of eq. (5.6) comes from the thin region with t − t 0 < ∼ 1/(ν 0 q 2 ). On the other hand, the two terms in the function h(q, t , t 0 ) involve powers of q that overcompensate the suppression by the q 2 , so that the integral is UV divergent. The spectra P 12 and P 22 display similar behavior. This illustrates that UV-tame initial conditions are a prerequisite for the validity of hydrodynamic perturbation theory. The physically motivated initial conditions discussed in section 4 satisfy this requirement. In particular, an initial bispectrum of the form (4.11) amounts to an exponential suppression of the large-q modes. In figs. (4) and (5) we illustrate that the resulting integrand h(q, t , t 0 ) exp(−σ 2 q 2 ) entering the q-integration in (5.6) vanishes now at large q. The scale at which this damping occurs is determined by the physics defining initial conditions (which is parametrized by σ = 0.4 and 1 fm, respectively), but for the entire class of initial conditions, the result is UV finite.
The above discussion clarifies how the general expectation that the applicability of fluid dynamics is limited to sufficienlty small momentum modes is realized technically for the particular case of evolving P 11 (0, t). For the general case, however, analytical insight is more difficult to obtain, and it is difficult to check a priori the convergence of the iterative solution of eqs. (3.26), (3.27). The numerical analysis can then provide a consistency check a posteriori.

Range of applicability of fluid dynamics
In the previous subsection, we emphasized that fluid dynamics can apply only if initial conditions are regulated in the UV. We now turn to estimating the scale below which this UV cutoff must lie for initiating a consistent fluid dynamic evolution. To this end, we consider corrections to the relativistic Navier-Stokes equations. In a second-order formalism, the constitutive relation for the shear stress π µν + 2 η P µν αβ ∇ α u β = 0 (5.9) (P µν αβ is the standard projector to the symmetric, transverse and traceless part of a tensor) gets supplemented by additional terms, in particular by Within the approximations made in section 2 one finds that the relaxation term becomes important when τ π |ḟ k | ∼ > |f k |.
As a rough estimate, one may assume that the time dependence off k is dominated by sound propagation with frequency ω = c S |k|. Eq. (5.11) corresponds then to τ π c S |k| ∼ > 1, and the range of applicability of fluid dynamics can be expected to cover momenta In general, both the relaxation time τ π and the velocity of sound c S will depend on the microscopic dynamics of the fluid. For the sound velocity, however, this dependence is relatively mild, and we assume the value c S = 1/ √ 3 for an ideal equation of state for the purpose of the following estimates. The relaxation time τ π is not known from first principles in QCD, but it is known for several non-abelian gauge theories with gravity duals in both the strong and weak coupling limit. It is expected to behave similarly in QCD. In particular, for N = 4 SYM theory, it is lim λ→∞ τ π = (2 − ln 2)/(2πT ) in the limit of infinite t'Hooft coupling λ = g 2 N c [42]. Assuming this value, i.e. assuming a strongly coupled liquid, we find A similar numerical estimate is obtained, if one assumes that the time dependence off k is dominated by the viscous damping rate 4νk 2 /3, where ν = η/(sT ) is the kinematic shear viscosity. In this case, the second order term remains small for and using η/s = 1/(4π), one finds |k| ∼ < 6π 2 2−ln 2 ≈ 6.7 T for the range of applicability of fluid dynamics.
It is interesting to note that the strong coupling limit of τ π used above can be written in the form lim λ→∞ τ π (2.7/T )(η/s), (where lim λ→∞ (η/s) = 1/4π) while the weak coupling limit is known to satisfy lim λ→0 τ π (5.9/T )(η/s) [43]. Such a monotonous relation between relaxation time and sound attenuation length is also expected on physics grounds. Therefore, the above estimates for the range of applicability of fluid dynamics may be extended to finite coupling strength, where they lead to |k| ∼ < cT /(η/s) with the constant c in the range 0.3 < c < 0.6. In the strong coupling limit, in which η/s = 1/4π, this reduces to the above estimates, but the range of validity of hydrodynamics will decrease with increasing dissipative effects. It will be very small in a weakly coupled system since η/s takes parametrically large values η/s ∼ 1/α 2 s in perturbation theory. In summary, in a strongly coupled liquid we expect that higher order derivatives in the fluid dynamic description may be neglected for the description of wavevectors up to |k| ∼ 8 T , or so. Any increase in dissipation (i.e. any increase of η/s) reduces this range of validity of a fluid dynamic description. On the other hand, the inclusion of second-order terms in the fluid dynamic expansion may extend this range of validity of fluid dynamics somewhat, while the expansion scheme is expected to break down always for sufficiently large momenta.

Linear and nonlinear evolution
In this section we solve numerically the time evolution for the spectrum and bispectrum derived in sect. 3, starting from the initial conditions of the model of independent point sources introduced in section 4.

Spectrum
To solve for the spectrum at leading order, which corresponds to a linear time evolution, we can simply use the propagator obtained in sect. 3.5. For the initial conditions specified in eqn. (4.8), we find At next-to-leading order, this expression gets supplemented by a term that involves the initial bispectrum of eq. (4.11) and a vertex term ∼ γ. We obtain Higher orders are suppressed by additional powers of 1/N .  Figure 6. The density-density correlator P 11 (k) at O(γ 0 ) (blue, dashed line) and O(γ 1 ) (purple, solid line). The fluid dynamic evolution was initialized at τ 0 = 1 fm, for an initial kinematic viscosity ν 0 = 0.04 fm corresponding to η/s = 1/(4π) and an initial temperature T 0 400 MeV (see eq. (2.6)). The initial number density of sources was N/A = 0.2 fm −2 . Sources are chosen with a Gaussian profile of width σ = 0.4 fm (see eq. (4.6)), that amounts to cut-off initial fluctuations at momentum scales above 3 fm −1 (see fig. 4). Results are shown for times t = 1.5, 2.5, 3.5, 4.5 fm (curves from top right to bottom left).
The time dependence of the components P 11 (k, t) , −k 2 P 12 (k, t) and k 4 P 22 (k, t) is shown in figs. 6-8. Since the velocity divergence θ and the field f are related through θ k = −k 2f k , these spectra correspond to d-d, d-θ and θ-θ correlations, respectively, that is the self-and cross-correlations between the logarithmic energy density d and the velocity divergence θ. As seen from figs. 7, 8, velocity correlations have already built up at t = 1.5 fm, even though the initial velocity field is assumed to vanish at t 0 = 0.75 fm. The magnitude of all spectra diminishes at late times as a result of the longitudinal expansion of the Bjorken background and the effect of viscosity.
Figs. 6-8 show the development of oscillatory behavior which corresponds to the appearance of acoustic waves induced by the strong pressure gradients in the relativistic fluid. The location of the minima is determined by the structure of the propagator that we discussed already in subsection 3.5 in the context of equations (3.36) and (3.37). In particular, the argument of the Bessel functions in these propagators involves the combination k t 3/4 , and this explains the continuous shift of the minima of the spectra towards smaller values for increasing time. As seen from figs. 6-8, this behaviour is generic for all fields and it continues up to the time when the plasma is so dilute that the fluid description breaks down. Therefore, the spectra shown here are particularly suited for a discussion of how the spectrum evolves with time to a form where it is dominated by longer wavelengths. This is relevant for addressing one of the central topics in the fluid phenomenology of heavy-ion collisions, namely how the spectrum of fluctuations at decoupling (that is experimentally accessible through final state hadrons) is related to the spectrum of initial fluctuations. There is a close analogy between this question and the question of how the power spectrum of acoustic oscillations at the time of recombination is fixed by the dynamics of microscopic interactions in cosmological scenarios. Access to the latter is provided by the spectrum of temperature fluctuations in the CMB and by the baryonic acoustic oscillations (BAOs) seen in the large scale structure of the Universe. Figs. 6-8 also illustrate clearly that nonlinearities contribute only a small correction to the final form of the spectrum for the parameter values favored by the phenomenology of heavy-ion collisions. In the above figures, this conclusion can be drawn from a comparison of leading O(γ 0 ) and next-to-leading O(γ 1 ) terms. The same conclusion is also supported from the analysis of the O(γ 2 ) backreaction of the initial spectrum that we shall discuss in the following. This is numerical evidence that a perturbative expasion of the spectrum and bispectrum in powers of γ has good convergence properties.
We conclude this section by discussing in more detail the dependence of the numerical results on the parameters chosen. Figs. 9 will serve as our reference, showing the density-density correlator for the same parameter choices as in figs. 6, but now evolved up to the later time t = 7.5 fm. We compare fig. 9 to figs. 10-11 to illustrate how the spectrum depends on the maximal momentum scale of fluctuations in the initial conditions that can controlled be in our model through the choice of the Gaussian width σ. For larger σ the linear spectrum is suppressed at nonzero values of k, as the corresponding modes are eliminated from the initial spectrum. Moreover, the nonlinear corrections are also suppressed because the loop corrections receive contributions only from a short range of lowmomentum modes. As a consequence, nonlinear corrections to P 11 (k) become negligible, see fig. 10. For small σ, a larger number of oscillations is visible in the spectrum, with large amplitudes relative to the large-σ case, and with somewhat larger nonlinear corrections (see fig. 11). We observe that while all these dependencies may have been expected on general grounds, their numerical size is sufficiently small to fully support the conclusions reached above.  Figure 11. Same as fig. 9, but for σ = 0.2 fm.
We finally consider variations in kinematic viscosity by comparing the case of fig. 9 corresponding to a minimal shear viscosity, to the case of fig. 12, where shear viscosity is a factor 2.5 higher. Again, as expected, viscous effects result in the suppression of the high-momentum modes, which are continuously depleted with time. As a result, at t = 7.5 fm the second (third) peak of the acoustic oscillations in the spectrum are suppressed (strongly suppressed) for ν = 0.1 fm relative to the ν = 0.04 fm case. If such an analysis could be elevated to more realistic models of heavy-ion collisions, the study of the relative sizes of the secondary peaks of the acoustic oscillations could provide a precision tool for further narrowing uncertainties in the extraction of η/s. As an aside, we note the analogy between this observation and the current phenomenological analysis of CMB fluctuations, in which the location and strength of peaks in the power spectrum are linked to material peroperties of the Universe.

Bispectrum
We next turn to the time evolution of the bispectrum. If the trispectrum is neglected, as in eq. comes from the linear evolution of the initial bispectrum. For the initial condition of eq. (4.11), it is of the form At next-to-leading order, this gets supplemented by a term that is independent of the initial bispectrum and leads to a non-trivial bispectrum at late times, even if B abc = 0 at t = t 0 . The full expression is of the form B abc (k, q, p; t) NLO = B abc (k, q, p; t) LO + 2 t t 0 dt g ad (k, t, t )g b1 (q, t, t 0 )g c1 (p, t, t 0 )γ dgh (q, p)g g1 (q, t , t 0 )g h1 (p, t , t 0 ) The contribution from the nonlinear interactions does not factorize in general. Note that here the leading and next-to-leading terms are of the same order in A/N .
In fig. 13 we depict the O(γ 0 ) and O(γ 1 ) contributions to the 111-component of the bispectrum at t = 7.5 fm for momenta that form an equilateral triangle, so that k = q = p. We observe that the linear evolution results in a spectrum with maxima and minima that reflect the oscillatory form of the propagator at late times. Contrary to the spectrum discussed in the previous subsection, the O(γ 1 ) correction to the bispectrum is significant and overwhelms the linear contribution. The nonlinear effects tend to further enhance the bispectrum in the momentum ranges in which the linear evolution has generated significant deviations from zero. While fig. 13 shows the bispectrum for a particular choice k = q = p of momenta, our finding of a large O(γ 1 ) contribution persists also for other choices of arguments, such as the one shown in fig. 14. These results indicate that, for our model, the spectrum is the dominant quantity determining the form of the bispectrum at late times. The initial condition for the bispectrum plays a secondary role. Since the form of the spectrum affects the bispectrum only at O(γ 1 ), the large O(γ 1 ) corrections seen in fig. 13 and fig. 14 do not necessarily indicate the break-down of a perturbative approach but they rather indicate that a qualitatively novel physics effect enters the bispectrum at O(γ 1 ). However, the consistency of the perturbative approach requires then that higher order corrections remain small compared to the O(γ 0 ) and O(γ 1 ) contributions. It is important, therefore, to examine the O(γ 2 ) correction to the spectrum, to which we turn now.

O(γ 2 ) corrections and the backreaction
In the following numerical exploration of O(γ 2 ) corrections to the spectrum, we limit ourselves again to the particular case P 11 (k, t) at k = 0. This simplifies the numerical task of solving eqs. (3.26), (3.27), and it allows for some analytical insight. To O(γ 2 ), the solution of eq. (3.26) for P ab (k, t) involves a nonlinear backreaction of the spectrum onto itself, and the analysis of P 11 (k = 0, t) allows us to study this effect. In fig. 15, the time evolution of the 11-component of the spectrum for vanishing momentum is shown in three different approximations: linear, O(γ 1 ) and O(γ 2 ). All three lines are very close to each other, which indicates that the nonlinear corrections are small. This feature persists also for other phenomenologically acceptable values of the parameters of the model, such as the ones we considered in the previous subsection. fig. 16 shows the corresponding ratios of P 11 (k = 0, t) evaluated at O(γ 1 ) and O(γ 2 ) respectively and divided by the leading order contribution. This figure allows one to identify better the time scale up to which nonlinearities develop, as well as the relative size of O(γ 1 ) and O(γ 2 ) corrections. In particular, we see that for the parameter choices made for fig.  16, nonlinearities contribute to a growing deviation from the linear order up to a time t 3 − 4 fm, while the ratio to the linear contribution remains almost unchanged for later times. This may be understood qualitatively by observing that the bispectrum contributes significantly to P 11 (k = 0, t) in (5.7) only in the thin integration region t − t 0 ∼ < 1/(ν 0 q 2 )  and that the most significant contributions comes from the highest momenta allowed by the cut-off in the initial conditions, q ∼ 1/σ. The resulting simple estimate t turn σ 2 /ν 0 provides indeed a good order of magnitude estimate for the time t turn at which the curves in fig. 16 stop rising. Numerically, the nonlinear effects amount to roughly 10% of the final value of the spectrum, and they develop up to the time O(2) fm/c after the onset of the hydrodynamic regime. The O(γ 2 ) correction is of comparable size but not larger than the O(γ 1 ) correction.
We discuss finally the development of a nonzero expectation value for the density. This is a qualitatively novel feature that arises due to the backreaction of the spectrum. Translational invariance requires that nonzero Fourier modes of the fields do not develop expectation values. However, the k = 0 mode can become nonzero during the time evolution of the system because of nonlinear effects, even if one sets intiallyφ a (t 0 ) = 0 through an appropriate field definition at the initial time, as we did in section 3.2. As we have argued already in section 3.2, the expectation value of f (t, x) has no physical meaning and can always be set to zero. On the other side, the expectation valueφ 1 (t) =d τ (t) is physical. In particular, eq. (3.25) gives d τ (t) = t t 0 dt d 2 q g ab (0 , t, t ) −q 2 P 12 (q, t ) + 4 9 ν 0 q 4 P 22 (q, t ) . (6.5) The last term, which is positive semi-definite, describes the effect of dissipation, while the previous one corresponds to a nonlinear convection effect. Physically, eq. (6.5) characterizes how the dissipation and convection of modes at finite momentum affects the spatially averaged background field. In fig. 17 we depict this time evolution of the k = 0 mode of the density field. In configuration space, this mode corresponds to the density expectation value. At the initial time, our choice of initial conditions enforces a vanishing value, d τ (t 0 ) = 0. However, the subsequent evolution leads to a quick deviation from zero, developing within a short time scale of ∼ 1 fm, similarly to what we observed in fig. 16. At late times, the density drops slowly because of the expansion of the Bjorken background. As we discussed in subsection 3.2, an appropriate time-dependent field redefinition, that would absorb this effect, is required for a precise statistical analysis of the late-time fluctuations around the shifted background.

Discussion
In this work, we have calculated the time evolution of two-and three-point functions (a.k.a. spectrum and bispectrum) of fluid dynamic fields for a simplified fluid model of ultra-relativistic heavy-ion collisions, the so-called Bjorken model. The analysis of spectra and bispectra is a standard tool in the study of cosmological perturbations, while the fluid dynamic formulation of perturbations is a recent development in the phenomenology of ultra-relativistic heavy-ion collisions. Even though the structural analogies between fluid dynamic applications in cosmology and in heavy-ion physics have been remarked on repeatedly, the present work is, to the best of our knowledge, the first in which results for spectra and bispectra are presented for a model of heavy-ion collisions. Here, we summarize our main findings: After having introduced a background-fluctuation splitting of fluid dynamic fields for the Bjorken model in section 2, we found in section 3 that the evolution equations for the spectrum and bispectrum of fluctuations can be written in a perturbative expansion that is of the same form as the expansion used for the study of cosmological perturbations. In section 4, we provided a simple but sufficiently generic model for the initial spectrum and bispectrum of fluctuations, through which this perturbative dynamics can be initialized. Within this setting, which is technically different from other fluid dynamic formulations of heavy-ion collisions, we analyzed the range of wavelengths up to which a fluid dynamic evolution is applicable, as well as the range of validity of a perturbative expansion for physically relevant model parameters. To this end, we discussed in detail the UV sensitivity of the fluid dynamic evolution. We observed that, depending on the choice of independent fluid dynamic fields, the UV sensitivity seems to arise either from linear or from nonlinear terms in the fluid dynamics (see subsection 5.1). We also found strong numerical evidence in support of the assumption that a perturbative expansion applies to the calculation of the spectra and bispectra of fluid dynamic fluctuations of realistic size (see section 6).
The correlated two-and three-point functions of fluid dynamic fields calculated here are not directly measurable in heavy-ion collisions. The reason is that we did not supplement the fluid dynamic evolution with a hadronization prescription, while the measurable momentum correlations of hadrons provide only convoluted information about spatial correlators. However, the spectra and bispectra calculated here provide arguably the most direct information about how an initial spectrum of fluctuations evolves as a function of time and how perturbations dissipate (i.e. weaken in strength). Any attempt of constraining the scale and nature of fluctuations at the initial time τ 0 via experimental measurements of fluctuations at freeze-out relies on understanding the dynamic evolution of the scale and size of fluctuations. We expect that the analysis of fluid dynamic spectra and bispectra will be useful to this end.
We finally point out that the dynamics and in particular the dissipation of fluctuations results in a backreaction that can be characterized by a growing, spatially independent zeromode on top of the unperturbed background field. We have explained in subsections 3.3 and 6.3 the physics encoded in this term and how it arises technically. While the present paper was motivated by the idea of adopting an analysis technique used in cosmology to a problem in heavy-ion physics, a companion paper discusses the physics of this backraction effect in cosmology [46].