Fluctuations around Bjorken Flow and the onset of turbulent phenomena

We study how fluctuations in fluid dynamic fields can be dissipated or amplified within the characteristic spatio-temporal structure of a heavy ion collision. The initial conditions for a fluid dynamic evolution of heavy ion collisions may contain significant fluctuations in all fluid dynamical fields, including the velocity field and its vorticity components. We formulate and analyze the theory of local fluctuations around average fluid fields described by Bjorken's model. For conditions of laminar flow, when a linearized treatment of the dynamic evolution applies, we discuss explicitly how fluctuations of large wave number get dissipated while modes of sufficiently long wave-length pass almost unattenuated or can even be amplified. In the opposite case of large Reynold's numbers (which is inverse to viscosity), we establish that (after suitable coordinate transformations) the dynamics is governed by an evolution equation of non-relativistic Navier-Stokes type that becomes essentially two-dimensional at late times. One can then use the theory of Kolmogorov and Kraichnan for an explicit characterization of turbulent phenomena in terms of the wave-mode dependence of correlations of fluid dynamic fields. We note in particular that fluid dynamic correlations introduce characteristic power-law dependences in two-particle correlation functions.


Introduction
It is a long-standing idea, first articulated by Landau in the 1950s, that the evolution of matter compressed in nuclear collisions lends itself to a fluid dynamical description [1]. On the experimental side, characteristic correlations of particle production with the event plane had been interpreted as qualitative support for fluid dynamic behavior since the very first relativistic heavy ion collision experiments at the BEVALAC in the 1980s [2]. Early qualitative predictions, based on fluid dynamics, include notably the argument [3] that the second harmonics of the azimuthal particle distribution (elliptic flow v 2 ) changes at mid-rapidity from out-of-plane to in-plane emission at higher center of mass energy. This was confirmed experimentally in heavy ion collisions at the Brookhaven National Laboratory (BNL) Alternating Gradient Synchrotron ( √ s NN < 5 GeV), and at the CERN Super Proton Synchrotron ( √ s NN < 20 GeV), for a review see e.g. Refs. [4]. Also, fluid dynamic arguments provided early on some qualitative understanding of the dependence of elliptic flow on particle species, and on the energy and rapidity dependence of the collective sidewards displacement of particle production at projectile rapidity (sidewards flow v 1 ) [4]. However, conclusions remained largely limited to the qualitative statement that the observed flow in semi-central collisions "retains some signature of the pressure in the high density region created during the initial collision" [5]. This changed soon after the start of the Relativistic Heavy Ion Collider (RHIC) in the year 2000, when several groups [6][7][8][9] noted that fluid dynamic simulations of Au+Au collisions at √ s NN < 200 GeV can account quantitatively for the main manifestations of collectivity at RHIC, including the dominant elliptic flow signal at mid rapidity and its dependencies on transverse momentum, centrality and particle species. These studies were based on simplified 2+1-dimensional simulations, following Bjorken's argument that the initial conditions for fluid dynamic fields are close to longitudinally boost-invariant, and that this boost-invariance is preserved by the fluid dynamics [10]. Moreover, early comparisons to RHIC data relied on ideal fluid dynamic equations of motion without dissipative effects. Soon afterwards, Teaney [11] observed in 2003 that even very small values of the ratio η/s of shear viscosity over entropy density induce dissipative effects that result in a sizable reduction of the elliptic flow signal. Therefore, to the extent to which uncertainties in the comparison of fluid dynamic simulations with data can be controlled quantitatively, measurements of collective flow in heavy ion collisions provide a sensitive tool for constraining transport properties of QCD matter. This is one of the main motivations for the development of more and more detailed fluid dynamic simulations of relativistic heavy ion collisions in recent years, see e.g. [12][13][14][15] for recent reviews. Ideal fluid dynamics is determined fully by the equation of state and conservation laws. For causal viscous fluid dynamics, transport properties and relaxation times enter in addition. But the equation of state, transport properties and relaxation times are in principle calculable from first principles of a given quantum field theory. Therefore, a comparison of fluid dynamic simulations to data of heavy ion collisions has the potential of constraining properties of QCD matter that are fundamental in the sense that they are most directly related to the QCD lagrangian. In practice, the bottleneck for such a program is the limited control over the initial data that are evolved fluid dynamically 1 . Early comparisons of fluid dynamic simulations with RHIC data [6][7][8][9] employed a set of smooth, event-averaged initial conditions that were specified via the average collision geometry in terms of a transverse energy (or entropy) distribution with vanishing flow at initial times. Even within this limited set of initial conditions, one observed that differences in the initial transverse profile of phenomenologically motivated models could result in variations of the dominant elliptic flow signal by up to 30% [16,17].
Within recent years, there has been a growing realization of the importance of event-by-event fluctuations in constraining the initial conditions of fluid dynamic evolution in heavy ion collisions. In particular, event-averaged initial conditions reflect the symmetries of the almond-shaped nuclear overlap region of finite impact parameter collisions and can therefore give rise only to dipole, quadrupole and higher even moments in the initial density distribution. In marked contrast, the azimuthal momentum distributions measured by all experiments at the LHC [18][19][20][21] and at RHIC [22,23] show a prominent third harmonic moment v 3 , as well as non-vanishing moments v 1 and v 5 in addition to the expected even ones. These structures had been attributed previously to other speculative effects ("Mach-cone", "ridge"), but as pointed out first by Alver and Roland [24] (see also Sorensen [25] for a related earlier suggestion) they emerge most naturally from the fluid dynamical evolution of initial density inhomogeneities. In addition, fluctuations increase the spatial eccentricity of initial transverse density distributions, and this accounts naturally for the fact that elliptic flow values remain sizable in the most central collisions and for smaller colliding systems [26]. There is by now compelling evidence that the dynamical evolution of fluctuating initial conditions is a prerequisite for a detailed quantitative understanding of flow in heavy ion collisions [27,28]. And since the various flow moments v n depend differently on the event-averaged initial state and its event-by-event fluctuations, analyzing the dynamical evolution of these initial fluctuations provides a novel tool for constraining the main uncertainty in fluid dynamical simulations of heavy ion collisions.
There has been a significant effort recently in studying fluctuating initial conditions in heavy ion collisions [29][30][31][32][33][34] and studying their propagation in full fluid dynamic simulations or transport models [35][36][37][38][39][40][41][42][43]. Precursors of these developments include e.g. Refs. [44][45][46]. The recent efforts focussed mainly on initial density inhomogeneities. But more general fluctuating initial conditions are conceivable. For instance, (non fluid-dynamical) initial fluctuations in the flow field u µ may be expected to accompany a fluctuating initial energy density profile [37]. Even if fluctuations in the initial spatial distribution may turn out to be sufficient to account for the measured flow components, it is clearly important to constrain such other conceivable sources of initial fluctuations since these may confound any quantitative interpretation of flow phenomena aimed at an extraction of η/s and other fundamental properties of QCD matter. This argues for treating fluctuations in all fluid dynamical fields democratically.
The present paper aims at supplementing the current discussion of fluctuating initial conditions with a model study in which the propagation of initial fluctuations can be followed in a very explicit, partly analytical way. To this end, we formulate the fluid dynamic evolution of fluctuations in all fluid dynamic fields around an event-averaged Bjorken flow profile. The inclusion of fluctuations in all fields will provide access to qualitatively novel features such as the dynamical evolution of vorticity. It will also allow us to discuss anew how one of the most characteristic manifestations of fluid dynamics, namely turbulence, can emerge in the specific expanding geometry of a Bjorken-like flow profile. The issue here is not whether heavy ion collisions can display fully developed turbulence: It has been pointed out previously (see e.g. Ref. [47]) that the relevant Reynolds numbers are typically larger than unity, since they are proportional to the inverse of the normalized viscosity η/s. However, the length and time-scales in heavy ion collisions are so small that Re < O(100) which is well below the conditions under which fully developed turbulence is expected. Rather, what is at stake is the suggestion first made by Mishra et al. [48] and further discussed by Mocsy and Sorensen [33] that the measurement of the harmonic flow coefficients v n for all values of n may provide information about the initial state similar to the power spectrum extracted from Cosmic microwave background (CMB) radiation. CMB analysis tools exploit the fact that Hubble expansion dampens vorticity fluctuations, so that the fluid dynamical evolution stays at all time scales in a linear non-turbulent regime. It is a priori unclear whether the same situation persists for small fluctuations in a Bjorken-type expansion, or whether small fluctuations can become seeds of turbulent behavior. Our discussion shall address this question and characterize the limitations of a linear treatment of fluid dynamic fluctuations in heavy ion collisions, thus gaining some insight into the conditions for onset of turbulent behavior 2 .
Our work is organized as follows. In section 2, we show that mild extension of models of fluctuating initial conditions give rise to fluctuations in all fluid fields. We note in particular that fluctuations in the velocity field can have in general a vorticity component as well as a divergent component, and that both components may be of comparable size. Furnished with this example that fluctuations in all fluid fields may be relevant, we formulate then the equations of motion for fluctuations around a Bjorken flow field in section 3, and we solve them in a linearized approximation of the evolution in section 4. We then turn in section 5 to the case of turbulent fluctuations, when non-linear contributions to the equations of motion matter. In particular, we provide a parametric argument that a large class of fluctuating initial conditions around Bjorken flow evolves at late times towards an effectively two-dimensional, turbulent system. Motivated by this observation, we recall in section 6 pertinent features of turbulence in terms of correlation functions of fluid fields. In section 7, we finally relate these remarks to heavy ion phenomenology by showing how correlations of fluid dynamic fields enter the one-and identical two-particle correlation functions in a blast wave model supplemented with fluctuations. In the conclusion, we finally summarize our main findings and provide a short outlook.

Fluctuating initial conditions and vorticity
We start our discussion with the prototype of an initial density inhomogeneity implemented in current event-by-event fluid-dynamical simulations, as used e.g. in Ref. [39]. Fluctuations in the initial spatial distribution are described at some initial time τ 0 and close to mid-rapidity y = 0 by a two-dimensional transverse energy density profile of the form Here, the coordinates x i denote for one specific heavy ion collision the positions of wounded nucleons in the transverse plane, as obtained from a Monte Carlo Glauber simulation, see e.g. [53]. A class of events corresponds then to a class of independently simulated distributions {x i }, each defining an energy density (x) with different, event-specific fluctuations and each setting the initial data for an individual fluid dynamic evolution. The normalization K in (2.1) can be constrained by data on the total transverse energy produced in the collision per unit rapidity [39]. The smearing parameter σ is a model-dependent input that sets the scale of spatial inhomogeneities. Fig. 1 illustrates that this model accounts for significant fluctuations in the transverse distribution of wounded nucleons and their corresponding energy density (2.1). In Fig. 1, we have chosen a smearing parameter σ = 0.4 fm, consistent with previous simulations of event-by-event fluid dynamics [39].  One possibility is to initialize fluid dynamical fields at initial time τ 0 with event-wise fluctuations in the transverse energy density (x), but with an exactly vanishing non-fluctuating flow field in the two transverse directions 1 and 2 at initial time τ 0 , and in the rapidity component of the flow vector However, a larger class of initial conditions is conceivable, since the initial state may also display fluctuations in the initial flow fields, and since both energy density (2.1) and the flow fields could depend on rapidity y. Fluctuations in the initial velocity fields have been discussed recently e.g. within a model of the early (Non-Equilibrium) dynamics based on free streaming [38]. For more discussions of initial flow and its influence on HBT radii see also Refs. [54][55][56][57][58]. Some fluctuations in the flow field will be generated by the fluid dynamic response to initial density fluctuations and may therefore be regarded as being implicitly included in the ansatz (2.1) -(2.3). However, such fluid dynamically generated flow fluctuations will be constraint in scale and size to the fluctuations in energy density, and they will lack by construction some qualitative features of the most general fluctuating flow field, such as vorticity. Vorticity 3 characterizes the solenoidal part of a general three-dimensional flow field u j (x), We note that the fluid dynamical evolution of fluctuations in vorticity and energy density decouples as long as these fluctuations are small and can be treated in a linearized evolution (see section 4). This is in marked contrast to fluctuations in the irrotational part of u j (sound modes) that can be driven by fluctuations in energy density. This indicates that one cannot expect to generate sizable values of vorticity by evolving initial conditions of the form (2.1)-(2.3). However, if fluctuations in vorticity are part of the initial conditions, then they will propagate and may display particularly interesting dynamical features, as discussed in sections 4 and 6, respectively. To set the stage for our discussion in later section, we demonstrate now that relatively mild extensions of (2.1) can lead naturally to fluctuations in velocity, including a non-vanishing solenoidal component (2.4). It is one arguably mild extension of (2.1) to associate the transverse region around a single wounded nucleon not with an energy density, but with an energy-momentum tensor T µν w , such that the initial energy-momentum tensor of the entire nucleus-nucleus collision takes the form In general, equation (2.5) can account for non-vanishing fluctuating initial conditions in both energy density (x) and flow u j (x). For instance, neglecting for simplicity non-ideal, shear viscous contributions to (2.5) and assuming an ideal equation of state (x) = 3 p(x), one can write the initial conditions for and u j at some fixed rapidity y and initial time τ 0 to linear order in u j in the form The transverse energy density associated to a single wounded nucleon is given by w (x) = T 00 w (x). Therefore, the 0-component of (2.6) defines an equation of the type (2.1) for the energy density, but it has also a fluctuating initial flow field defined by the spatial components of (2.6), . (2.8) Here u j w is defined by writing equation (2.6) for the energy momentum tensor associated to a single wounded nucleon. The size of initial fluctuations in the velocity field (2.8) depends on how the initial transverse motion associated to the the single wounded nucleon is modeled. Taking guidance from blast wave models, one may choose for u j w e.g. an azimuthally symmetric radial flow field with some radial dependence w, u j w (x) = for j = 1, 2 and u y = 0, say. For such an ansatz, one checks easily that (2.8) defines in general a flow field of non-vanishing transverse divergence and non-vanishing longitudinal vorticity, (2.9) A more general ansatz may be based on the observation that in general, the transverse energy deposited by a single wounded nucleon in a finite window of rapidity recoils against transverse momentum outside this rapidity window. This argues for a net transverse velocity component v i associated to the contribution of each wounded nucleon in (2.5). To illustrate this effect, we assume that each wounded nucleon in the sample shown in Fig. 1 is associated with a non-relativistically small random transverse velocity component v i , drawn from a Gaussian distribution of width |v| = 0.1c 4 . The resulting initial transverse flow field (2.8) is shown in Fig. 2. By comparing to Fig. 14 of Ref. [38], we note that a seemingly comparable flow field can be obtained in a model of early non-equilibrium dynamics based on free streaming where the contribution of every particle is taken to be delocalized in position space over a small volume. In Fig. 3, we plot the absolute value of the longitudinal vorticity |ω 3 | = |∂ 1 u 2 − ∂ 2 u 1 | and the transverse divergence |∂ 1 u 1 + ∂ 2 u 2 | for the flow field shown in Fig. 2. Inspection of this figure shows that both components fluctuate with a similar magnitude and over similar transverse length scales.
For illustrative purposes, we have chosen in Fig. 3 velocity fluctuations of an average strength |v| = 0.1. We emphasize, however, that none of our conclusions in this section or in the following sections depends on the precise numerical choice for |v| . In particular, repeating the analysis of Fig. 3 for much smaller values of |v| would equally well support the only conclusion that we draw from it, namely that the irrotational and solenoidal components of the velocity may be of comparable size. We observe at this point that there is no general model-independent argument for the relative size fluctuations in u j and . This motivates us in the following sections to treat all conceivable sources of event-by-event fluctuations on an equal footing.  . The left hand side shows the absolute value of vorticity |∂ 1 u 2 − ∂ 2 u 1 | for the velocity field shown in Fig. 2. Similarly, the right hand side shows the absolute value of the divergence of the fluid velocity |∂ 1 u 1 + ∂ 2 u 2 | for the same velocity field. The color coding is the same on both sides.
effective fluid dynamical description is large and comprises bulk hadron production up to a few GeV in transverse momentum [12][13][14]. In this section, we recall first shortly the fluid dynamic equations of motion. We focus then on the Bjorken model that defines a particularly simple expanding geometry and that encodes important features of a relativistic heavy ion collision. Regarding the Bjorken model as defining the average fluid dynamical field, we then discuss how fluctuations in energy density and flow propagate in the expanding geometry of a relativistic heavy ion collision. The relativistic hydrodynamic equations for a fluid without any conserved charges read [12][13][14] D + ( + p)∂ µ u µ + π µν ∆ α µ ∂ α u ν = 0, Here is the energy density and p is the pressure in the fluid rest frame, π µν is the viscous part of the energy-momentum tensor in the Landau frame, u µ π µν = 0. The partial derivative ∂ α must be replaced by the covariant derivative ∇ α if one works with coordinates other than cartesian. We work with a cartesian metric of signature g µν = diag(−1, 1, 1, 1). The matrix ∆ µν projects to the subspace orthogonal to the fluid velocity, ∆ µν = g µν + u µ u ν . The derivative in the direction of the fluid motion is D = u µ ∂ µ . The viscous part of the energy-momentum tensor can be expanded in a derivative expansion. To lowest order it vanishes, leading to ideal hydrodynamics. The first order contains shear and bulk viscocity terms is transverse (orthogonal to u µ ) and traceless. To second order in the gradient expansion, the fluid dynamic equations of motion contain various relaxation time corrections. It is a peculiar feature of a second order approximation that the evolution equations are hyperbolic and propagation is limited to the forward light cone even for perturbations of large wave-vector k. For this reason, second order fluid dynamics is often referred to as causal viscous fluid dynamics. However, this wanted feature of causality is not guaranteed to persist in higher orders of the gradient expansion. More generally, fluid dynamics is a long distance effective theory that by its very construction cannot be expected to be reliable for large wave-vectors. For the propagation of fluctuations with small gradients (i.e. small wave-vectors k), second order fluid dynamics will make only small corrections to a first order treatment. For this reasons and to keep the formalism simple, we restrict the discussion in the present paper to first order fluid dynamics. For the case of vanishing bulk viscocity ζ, the corresponding equations of motion read There is evidence that dissipation in a heavy ion collision is mainly due to shear viscosity, and we therefore neglect bulk viscosity for the remainder of this paper.

The Bjorken model
We are interested in studying fluid dynamic fluctuations in the expanding geometry of a relativistic heavy ion collision. The Bjorken model is arguably the simplest formulation of a corresponding expanding geometry. Motivated by the idea that in nuclear collisions at ultra-relativistic energy, particle production is almost flat in rapidity, Bjorken [10] proposed to formulate initial conditions for fluid dynamic fields that are independent of space-time rapidity y = arctanh( If this condition is satisfied at some initial proper time τ = τ 0 , then it persists for all proper times τ = x 2 0 − x 2 3 throughout the evolution. This renders the longitudinal evolution trivial, and the numerical task simplifies to the solution of a (2+1)-dimensional problem.
Under the further simplifying assumption that the initial transverse flow field vanishes at initial times and that transverse gradients in energy density are absent, the Bjorken model reduces to an effectively (1+1)-dimensional toy model that allows for an explicit analytical treatment. In this case, the evolution equation for the energy density becomes where energy density and pressure are related by the equation of state = (p). For what follows, it will be useful to rewrite this equation in terms of the enthalpy and the kinematic viscosity ν = η/w .
Throughout this work, we shall neglect terms that are parametrically suppressed by powers of ν/τ compared to some other term of the same structure. With this approximation, equation (3.9) becomes independent of shear viscosity. As will become clear in the following, however, the dominant effect of shear viscosity on fluctuations can be retained within this approximation.
With the approximation ν/τ 1 and for an ideal equation of state = 3p, one then finds the characteristic time dependencies of the Bjorken model for energy density and temperature For a time-independent normalized viscosity η/s, the ratio decreases. Therefore, replacing the bracket in Eq. (3.9) by unity is an approximation that is consistent with the late time behavior. We note at this point that in situations with strong (non-Gaussian) fluctuations, the evolution equation for averaged fields such as energy density (τ ) = (τ ) gets modified by additional terms, see the discussion in Sect. 6.1. For the present paper we assume that these modifications are small and can be neglected.

Fluctuations on top of a Bjorken background field
In this section, we formulate the theory of the dynamics of fluctuations on top of a Bjorken background field without transverse gradients. That means that the hydrodynamical fields u µ , when averaged over many events follow a Bjorken type solution. However, locally and for a particular event we expect deviations which we want to investigate in more detail. We have chosen a Bjorken background field for our study mainly for two reasons. First, the analytical simplicity of this background will allow for a particularly explicit discussion. Second, the Bjorken model contains essential features of realistic expansion scenarios of relativistic heavy ion collisions. 5 We denote fluctuations on top of the Bjorken flow u µ Bj = (1, 0, 0, 0) by relaxing the constraints (2.2) and (2.3) and allowing for local fluctuations in the transverse and rapidity components, Here and in what follows, we work in light-cone coordinates τ, x 1 , x 2 , y with metric g µν = diag(−1, 1, 1, 1/τ 2 ). The latin index j is summed over 1, 2, y and the corresponding three-dimensional metric reads g ij = diag(1, 1, 1/τ 2 ). We consider small local fluctuations in the sense that u j u j (x) 1. In the following we neglect terms that are parametrically suppressed due to u j u j (x) 1 or due to ν/τ 1 compared to other terms with the same combination of derivatives of the velocity and pressure fields. We note that for every combination of derivatives there is one term of lowest order which is not neglected and that the main physical effects of viscosity -the damping of velocity fluctuations and the dissipation of kinetic energy to heat -are correctly taken into account. With this approximation scheme we find from Eq. (3.4) and (3.5) the following equations governing the velocities in the transverse plane (j = 1, 2) and in rapidity direction (j = y) (3.14) 5 It has been pointed out repeatedly that despite its simplicity, the Bjorken model without transverse gradients retains important features of the early time dynamics of heavy ion collisions. The argument is based on the observation that event-averaged initial energy density distributions show typically only small transverse gradients in the central region of the transverse plane; the central region may indeed be approximated by the ansatz u = 0, (x) = .
The transverse evolution of this initial condition may then be thought of qualitatively as being dominated by a rarefaction wave that moves from the outside (vacuum) to more and more central positions in the transverse plane at late times. At a given position in the transverse plane, the dynamical evolution may be viewed as being characterized by the effectively (1+1)-dimensional Bjorken model up to the later time at which the rarefaction wave reaches the corresponding transverse position. Based on such considerations, the Bjorken estimates for the time-dependence of energy-density (3.10) and temperature (3.11) are used regularily in simple phenomenological estimates.
Here, the first two terms describe the change in the velocity along the direction of the fluid motion. They result from writing Du j = u µ ∂ µ u j for small deviations from the Bjorken background. The terms in the first square bracket account for two effects. One is the acceleration of the fluid due to pressure gradients in the transverse direction. The second term proportional to u j is dominated by the decrease of pressure for increasing τ . This dilution of the fluid leads to an acceleration in the direction of u j . Finally, there are effects of viscosity that are similar to the corresponding term in the (non-relativistic) Navier-Stokes equation.
In addition to eq.(3.14), one finds under the same assumptions for small local fluctuations around the Bjorken background the equation of motion for the internal energy density Here, the first two terms describe the change along the fluid direction of motion. The first square bracket describes dilution effects; the first term ∼ 1 τ is due to the expansion of the Bjorkenbackground in the longitudinal direction while the second term measures the effect of a possible dilution (or compression) in the transverse and rapidity directions. Viscous correction that are parametrically suppressed due to η/(wτ ) 1 have been dropped, and the remaining dissipative contribution to the evolution of internal energy are given in the last bracket of (3.15). They describe how kinetic energy is transferred from the macroscopic motion of the fluid to internal energy.
It will turn out to be useful to rewrite Eqs. (3.14) and (3.15) in terms of rescaled fluctuations in velocity, and for a rescaled time variable (Of course, this rescaled time t is not the time variable x 0 in the laboratory frame.) In what follows, we also absorb deviations from the τ -dependence of Bjorken's energy density (3.10) in terms of the quantity .
Finally, we assume that the shear viscosity ν follows the Bjorken behavior (3.12). That means, we neglect local fluctuations in the kinematic viscosity since they are expected to have only a minor effect. The kinematic viscosity ν 0 can then be written as For an ideal equation of state = 3p, using 1 w dp = 1 sT ∂p ∂T dT = 1 T dT , and neglecting a dissipation term ∼ ν that is of higher power in the velocity field, one can show that Eqs. (3.14) -(3.15) lead to the equation for the (rescaled) velocity (j = 1, 2, y) and for the quantity d τ Here and in what follows, we denote the expansion scalar ∂ j v j of the rescaled velocity fields by In the following we will use both the representation of the equations of motion in Eq. (3.14), (3.15) and the one in (3.20), (3.21). In particular the discussion of linear fluctuations in section 4 will be largely based on (3.14) and (3.15), while for the discussion of non-linear fluctuations in section 5 the representation (3.20), (3.21) will be more appropriate.

Linear fluctuations
In this section we discussion the evolution of fluid dynamical fluctuations that are small enough to neglect non-linear terms in (3.14) and (3.15). In addition, it is assumed that the deviation of the temperature field from the homogeneous background is small, d 1. The resulting linearized equations describe laminar flow. They apply to systems with sufficiently small Reynolds number, as we shall discuss in section 5.
For a fixed time τ and given spatial boundary conditions, one can divide the velocity field uniquely into a solenoidal part with vanishing divergence, and an irrotational part with vanishing curl, where Div u = ∂ j u j and Curl u is defined in Eq. (2.4). We recall that the derivative operators in (4.1) introduce an explicit τ -dependence since they involve the three-dimensional metric g ij = diag(1, 1, τ 2 ). Therefore, the splitting of u j into a solenoidal and an irrotational part does not commute with the τ -derivative.
The field u I j can be represented by the expansion scalar (3.22) It will be convenient to write this expansion scalar as a sum of a transverse and a longitudinal contribution, θ = θ T + θ y . From (3.14), (3.15) we obtain the linearized equations Here, we recall that the quantity d denotes the logarithmic temperature (3.18). The symbol h j in (4.6) takes the values h 1 = h 2 = −2 and h 3 = 1. Interestingly, the vorticity modes ω j decouple from the velocity divergence θ, the logarithmic temperature field d and from each other. Using Fourier decomposition with respect to the spatial argument, their diffusion-type equation of motion can be directly solved We assume here a constant ν 0 as defined in (3.19) and we use (3.17). One sees from equation (4.8) that vorticity modes for essentially all wave vectors are dominated at late times by an exponentially decaying function with a decay time set by the product of kinematic viscosity and the square of the wave vector. A somewhat unusal case is the time evolution of modes with k y = 0 where the exponential damping term is modified by a term ∝ 1/ √ t. In particular, for k 2 1 + k 2 2 = 0 but k y = 0 the vorticities do not decay exponentially for τ → ∞. In addition, the exponential decay is modified by a term that decreases algebraically for the transverse components ω 1 and ω 2 and that increases in the longitudinal components ω 3 . For finite times, the algebraic increase of ω 3 can overcome the exponential decay with viscosity.
In Fig. 4, we plot the solution (4.8) of the linearized fluid dynamic equations of motion for phenomenologically motivated input values, namely a small normalized shear viscosity η/s = 1/4π and a temperature of 500 MeV at initial time τ 0 = 1 fm/c. This translates into an initial kinematic viscosity ν 0 0.03 fm. Most generally, Fig. 4 illustrates the interplay between an exponential decay In addition to the evolution equations for vorticity, there are equations for θ T , θ y and d that we discuss now. These describe sound waves and are best solved in Fourier space. We concentrate first on a wave traveling in the transverse direction x 1 , corresponding to k 1 = 0, k 2 = k y = 0. In this case, eq.(4.4) decouples from the others and can be integrated immediately, Equations (4.3) and (4.5) are coupled, 11) and depend also on the solution for θ y . Concentrating for simplicity on the case θ y = 0 and eliminating d, one finds the second order differential equation (We have dropped a term suppressed due to ν/τ 1 in the second bracket.) For vanishing viscosity, this equation can be solved in terms of Bessel functions. The two linear independent solutions are (4.13) One finds an oscillating behavior with the amplitude increasing algebraically with time proportional to (τ /τ 0 ) 1/6 . For non-vanishing viscosity ν there is also an exponential decay which is larger for large wavevectors k 1 . For small values of k 1 τ √ 3 the two independent solutions of (4.12) are proportional to τ and τ 1/3 , respectively. Eq. (4.11) implies that the temperature field d grows according to τ 2 and τ 4/3 for these two solutions. For k 1 = 0 and late enough times τ , the solutions in (4.13) always have an oscillating behavior, however. For large wave vector k 1 one can neglect the terms ∼ 1 τ and ∼ 1 τ 2 in (4.12). Up to viscous damping, the solution corresponds to a perturbation that propagates with the velocity of sound c S = 1/3 into the transverse direction.
We show the time evolution governed by (4.12) on the left hand side of Fig. 5. In close similarity to the case of vorticity, we observe that also fluctuations in the transverse velocity divergence can persist over time scales relevant in heavy ion collisions. For sound waves in the transverse direction, the fluid acts like an efficient low-pass filter, allowing for the unattenuated (or even slightly amplified) passage of fluctuations on sufficiently large scales 1/k 1 ≥ 1 fm, while filtering out fluctuations on smaller length scales 1/k 1 < 0.25 fm.
We finally turn to sound waves traveling in the rapidity direction y, i. e. k y = 0, k 1 = k 2 = 0. Now, equation (4.3) decouples ∂ τ θ T − 1 3τ θ T + ν 1 τ 2 k 2 y θ T = 0 (4.14) and can be integrated immediately, The equations for θ y and d are coupled and depend on θ T , Concentrating on θ T = 0, this yields the following second order differential equation for θ y (we drop again a term suppressed due to ν/τ 1), The two linear independent solutions to this equation for ν = 0 are For k 2 y 16/3 this becomes τ and τ −5/3 . Equation (4.16) implies that this corresponds to perturbations in the temperature field d that grow like τ 2 or decay like τ −2/3 , respectively. For k 2 y > 16/3 the solutions (4.17) correspond to an oscillation with a period that increases as a function of time and an additional decrease in the amplitude. On the right hand side of Fig. 5, we plot this behaviour of θ y for small non-vanishing viscosity, when the algebraic increase is modified by an exponential decay. We observe again that fluctuations of modes with wave vectors k y ∼ O(1) can persist or can be amplified over time scales commensurate with the expected expansion duration of heavy ion collisions.
In the limit of large k 2 y , one can translate (4.18) back to position space and one finds that it corresponds to an excitation that propagates in the rapidity direction according to In general, for a sound mode propagating fastly into an arbitrary direction described by a large wavevector k = (k 1 , k 2 , k y ), we expect the propagation velocity ( ∂x 1 ∂τ , ∂x 2 ∂τ , ∂y ∂τ ) to satisfy This condition is satisfied by the sound modes in θ y and θ T discussed here.

Turbulent fluctuations
In full generality, the non-linear equations (3.14) and (3.15) or, equivalently, (3.20) and (3.21) are difficult to analyze. One can always use a splitting of the velocity into a solenoidal part that carries vorticity ω j and an irrotational part described by the divergence θ. The nonlinear terms in the equation of motion will lead to couplings between these fields and to the logarithmic temperature field d, however. Intuitively, one expects that perturbations in the fluid divergence propagate quickly in the medium with the characteristic velocity given by the velocity of sound. The fluid velocity (u 1 , u 2 , u y ) can be small compared to sound propagation. This is often characterized by a small Mach number For the description of the part of the fluid velocity that carries vorticity one can often assume in this case a vanishing divergence, θ = 0, since the fast sound modes can be viewed as decoupled from the slower modes dominating the solenoidal part of the fluid velocity. In the present section we will study the equations of motion in this situation and show that there are some interesting parallels to non-relativistic, incompressible fluids.
For our discussion in this section, it will be useful to work with the rescaled velocities v i introduced in section 3. For θ = ϑ = 0, Eq. (3.20) simplifies then to Due to the the solenoidal constraint the temperature field d is not independent of the velocity field. More specific, by taking the divergence of (5.2) one derives This is an instant of the Poisson equation for d. For given boundary conditions it can be inverted to yield d as a (non-local) functional of the velocity field. Equation (5.2) has some interesting features 6 . It takes the form of a two-dimensional Navier-Stokes equation in situations where v y = 0 and where the dependence of v 1 , v 2 on y can be neglected. Moreover, for a large class of initial conditions at time τ = τ 0 , the evolution becomes effectively two-dimensional for late times τ /τ 0 1. Indeed, both the non-linear velocity term that couples v y to v 1 and v 2 and the damping term involving the derivative with respect to rapidity contain factors that decrease as 1/τ 2 . Similarly, the solenoidal constraint (5.3) assumes its two-dimensional form in that limit.
Motivated by the observation that for the particular expansion geometry of the Bjorken model, fluctuations are governed by an evolution equation that reduces at late times to a non-relativistic Navier-Stokes equation (in rescaled time coordinates), we now discuss the conditions for a non-linear turbulent evolution of fluctuations by the very concepts that have proven useful in the classification of solutions of the Navier-Stokes equation. To this end, we consider situations where the velocities change notably on distances of order l in the transverse direction or for rapidity differences ∆y. 6 We mention as an aside that in addition to the obvious rotation symmetry in transverse direction and the translational symmetries in transverse and rapidity directions, Eqs. (5.2) and (5.3) have also the following space-time symmetry for constant velocity V j . Indeed, contributions from the time derivative and from the non-linear advection terms cancel. For V y = 0 this corresponds to Galilean symmetry in the transverse plane while the situation is more complicated for V y = 0. However, it is not so clear whether this invariance is very useful since there is no translational symmetry with respect to time.
The damping term in Eq. (5.2) leads then to a damping rate (inverse relaxation time) that can be characterized in terms of the dimensionless number κ = l 2 τ 2 ∆y 2 . (5.6) This damping rate is of order ν 0 /l 2 for κ 1, and it is of order ν 0 /(τ 2 ∆y 2 ) for κ 1. For characteristic velocities v T in transverse, respectively v y in rapidity direction, the flow can then be characterized in terms of the Reynolds numbers and Obviously, these definitions can be extended to intermediate κ, as well.
If both Reynolds numbers are small, Re (T ) 1, Re (y) 1, the resulting flow pattern is expected to be laminar. The viscous damping term dominates then over the nonlinear terms in Eq. (5.2) and velocities are expected to follow a regular behavior dominated by an exponential decay in time. More specifically, if the non-linear term in Eq. (5.2) can be neglected one can use Fourier decomposition with respect to the spatial arguments and one finds the solution This resembles closely the behavior of vorticity in Eq. (4.8).
In contrast, if both Reynolds numbers are large, Re (T ) 1, Re (y) 1, one expects a turbulent regime. The nonlinear terms dominate now over the viscous damping term in Eq. (5.2) and the velocities will change rather irregularly from point to point both in the transverse plane and for different values of the rapidity variable y.
We consider next the case Re (T ) 1 and Re (y) 1. For κ 1 the derivatives with respect to x 1 and x 2 effectively drop out from Eq. (5.2). One might expect a sort of one-dimensional turbulent behavior. However, the unusual explicit time dependence of the non-linear and damping terms could spoil this conclusion. Also for κ 1, we are unable to predict consequences of (5.2) without additional explicit calculations. However, because of the late time behavior of Re (y) 1, we do not expect that this case is particularly relevant for the simulation of heavy ion collisions (see the more detailed argument in the paragraph below).
In the opposite case Re (T ) 1, Re (y) 1, and for κ 1, all terms containing derivatives with respect to y can be neglected and Eq. (5.2) can be seen as a set of equations describing fluid motion in 2 + 1 dimensions with rapidity entering only as a parameter. With respect to the transverse coordinates one would expect a turbulent flow pattern. We remark that κ as defined in Eq. (5.6) decreases with time τ and the ratio of the two Reynolds numbers (which is independent of κ) increases with time Therefore, in contrast to the case Re (T ) 1, Re (y) 1 discussed above, the case Re (T ) 1, Re (y) 1 can persist at late times, and it can be reached dynamically. For completeness, we mention finally the region κ 1 when derivatives with respect to y enter in the damping term. Considered as a function of the transverse coordinates x 1 and x 2 , the velocity field will be damped by locally varying rates due to irregularities with respect to the rapidity argument. However, this damping is dominated by the first non-linear term in Eq. (5.2) so that the fluid will again behave turbulent with respect to the transverse coordinates.
The above discussion indicates that the case Re (T ) 1, Re (y) 1, κ 1 may be of particular relevance for the discussion of the onset of non-linear turbulent behavior in heavy ion collisions, since it can be reached dynamically and since it persists at late times. For a large set of initial conditions, we therefore expect that the fluid dynamics of heavy ion collisions evolves towards a system with effectively two-dimensional turbulent behavior. To address the issue to what extent the evolution towards turbulence could be completed within the finite duration of a heavy ion collision, let us finally put some numbers to this parametric discussion. Because of the late time limit of eq. (5.10), we focus on the transverse Reynolds number We consider typical values for a heavy ion collision at LHC energies, e.g. T ≈ 0.3 GeV and a length scale l ≈ 5 fm. Taking for the transverse velocity a fraction of the velocity of light, say v T = 0.1 c, one finds v T l T ≈ 1 and For small values of the normalized viscosity η/s < 1, one therefore expects a transverse Reynolds number Re (T ) that is larger than unity but not many orders of magnitude larger than unity (Re (T ) < 100). Such values are not sufficiently large to expect fully developed turbulence. A value Re (T ) > 1 indicates, however, that the system can be driven outside the region of validity of a laminar evolution and that it may display features indicative of the onset of turbulent behavior.

Qualitative features of turbulence
In the previous section, we have discussed the general conditions under which the time-evolution of fluctuations in a Bjorken expansion scenario can be expected to lead to the onset of turbulent behavior. Despite the relativistic nature of the system under study, we found that upon coordinate transformation, the time evolution of fluctuations is governed by an equation that takes the form of a two-dimensional Navier-Stokes equation at late times. Although the Navier-Stokes equation presents still many deep problems, much is known about fully developed turbulence in this system. Classical achievements include in particular the scaling theory by Kolmogorov [59] for three-dimensional turbulence and its extension to the two-dimensional case mainly by Kraichnan [60] and Batchelor [61]. For reviews of this field, see [62] To the best of our knowledge, the Bjorken scenario considered here is the first example of a relativistically expanding three-dimensional scenario that under suitable initial conditions evolves dynamically into a system with effectively two-dimensional turbulent dynamics described by a nonrelativistic Navier-Stokes equation. Turbulence in the two-dimensional case is known to display some characteristic qualitative differences in comparison to the three-dimensional case. Although the present section contains, strictly speaking, no novel results, our finding of an effectively twodimensional fluid dynamic propagation of fluctuations at late times prompts us to discuss here the pertinent features of fully developed turbulence that apply to the Bjorken scenario for sufficiently large values of Re (T ) . In particular, we point to the phenomenon of an inverse cascade that exists only in the case of two-dimensional turbulence and that provides a unique mechanism for enhancing fluctuations on large spatial scales during the evolution towards turbulence. To set the stage for this discussion, we introduce first shortly a statistical description of fluctuations in the fluid velocities and energy densities.

Fluctuation spectra
In general, fluctuations in the fluid velocities and the energy densities can be described statistically in terms of a τ -dependent probability distribution Eq. (6.1) describes an ensemble of events with equal "macroscopic" properties such as nucleon number, center of mass energy and impact parameter. Bjorkens model of a unique fluid velocity and energy density is then recovered by taking the probability distribution in (6.1) to be infinitely narrow.
We consider a generalization of Bjorkens model where the assumed symmetries (boost invariance in the longitudinal direction and translation and rotation invariance in the transverse plane) are broken by the fluctuations for a particular event but hold in a statistical sense. This means that the probability distribution (6.1) is invariant under these symmetries. This implies in particular that the expectation value of velocity is given by for all times τ > τ 0 and for all values of x 1 , x 2 , y. Similarly the expectation values of thermodynamic scalars such as , p, s, T etc. will be a function of τ , only. In general, however, the τ -dependence of these expectation values will differ from the ones obtained by Bjorken due to non-linear effects of fluctuations. Beyond the expectation values, the probability distribution (6.1) can be characterized in terms of correlation functions. In particular, the two-point correlation function of velocities at equal time τ is defined by (i, j = 1, 2, y) Due to translational symmetries this depends only on the differences in the spatial coordinates. Similarly, we define and the cross-correlation The generalization to other scalar quantities such as energy density or pressure p or to un-equal time arguments is obvious. It is useful to introduce also the Fourier decomposition This generalizes trivially to the other functions. We note that for symmetry reasons, one has G ij u (τ ) = 0 for i = j and G 11 u (τ ) = G 22 (τ ). For a solenoidal fluid with ∂ j u j = 0 one has In this framework, characterizing the fluid evolution of a heavy ion collision amounts to make statements about the form of correlation functions such as G ij u (τ, k) etc. at some given time τ . It is clear that the form of these correlations depends strongly on the initial conditions. For example, if the initial fluctuations at time τ 0 are small enough (or, equivalently, the viscosities large enough to have small Reynolds numbers) so that the linearized equations (4.3) -(4.6) can be applied, one can derive from them linear evolution equations for the set of correlation functions G ij u , G d and G j ud . The form of these functions at time τ is then directly linked to the corresponding functions at time τ 0 .
The situation is much more complicated in the presence of non-linear contributions to time evolution, even if the system is far from the conditions of fully developed turbulence. Indeed, if one tries to use the non-linear equations (3.14) and (3.15) to derive evolution equations for the two-point correlation functions, one finds that three-point correlations get involved, as well. The evolution equation for these involve even higher correlations and so on. This is an instant of the well-known closure problem in the statistical description of fluids. The same problem appears in nonperturbative formulations of quantum and statistical field theories. No exact analytical solutions are known and advanced techniques from statistical mechanics and field theory are needed to find approximate ones. The scaling theory of Kolmogorov provides some insight into these problems.

Scaling theory of turbulence at large Reynolds number
Any fluctuation present in the initial conditions of a relativistic heavy ion collision can be regarded as a source of non-thermal, say 'mechanical' energy. Most generally, one would like to understand to what extent this mechanical energy dissipates to thermal energy, and to what extent it does not but leaves characteristic, dynamically evolved fluctuations visible in final state observables. In section 4, we have shown examples of the dissipation (or amplification) of some fluctuation modes in a linearized description that applies to very small Reynolds numbers. Here, we want to comment on the opposite case of very large Reynolds number or very small viscosity. For large Reynolds numbers, we have found in section 5 that fluctuations on a Bjorken background field are governed by an evolution equation that takes for sufficiently late times the form of a nonrelativistic Navier-Stokes equation. Also, for situations of sufficiently small Mach number, there is a rationale for setting θ → 0. Therefore, our discussion of the case of large Reynolds number will reduce to recalling pertinent features of Kolomogorov's theory of homogeneous turbulence for a non-relativistic, incompressible fluid at very large Reynolds numbers.
For a non-relativistic fluid, it is common to denote the fluid kinetic energy per unit mass by , where the square brackets stand for spatial averaging. The rate of dissipation to thermal energy for a three-dimensional evolution is given by This shows that energy dissipation is mainly due to fine structures of the velocity field for which the gradients are large; dissipation is mainly taking place at large wave-vectors k. If one now decreases the viscosity, one finds finer and finer structures emerge so that the energy dissipation rate ε diss remains positive. In fact, the mean-square vorticity 1 2 [ ω 2 ] (also called enstrophy) grows ∼ 1/ν for ν → 0. This is possible since 1 2 [ ω 2 ] is not only given by the vorticity present at some initial time or generated from an external driving force. Rather, it can be generated also by non-linear terms in the three-dimensional Navier-Stokes equation. The mechanical energy is cascaded from the large length structures to the smaller ones by virtue of the non-linear terms in the Navier-Stokes equation. This is the famous cascade picture of Richardson [63]. In his words: "Big whorls have little whorls, Which feed on their velocity; And little whorls have lesser whorls, And so on to viscosity." Let us now come to the situation in two spatial dimensions. The most important difference to the three-dimensional case concerns the evolution of mean-square vorticity in the absence of an external driving force (vorticity ω = ∂ 1 v 2 − ∂ 2 v 1 has now only a single component) d dt (2-dim. case) (6.11) This shows that 1 2 [ω 2 ] never increases as a function of time. This in turn implies that energy per unit mass is conserved for vanishing vorticity, d dt A cascade of mechanical energy from large structures into smaller ones where it is finally dissipated is therefore not possible 7 . Without going further into the detailed mechanism of two-dimensional turbulence let us now attempt to transfer some of the insights that have been gained in this field to heavy ion physics, in particular fluctuations around Bjorken flow. We consider the case of large fluctuations and small kinematic viscosity ν so that the Reynolds number Re is large. Also, we concentrate on the limit of large time τ where (5.2) becomes two-dimensional. Similar to the non-relativistic case one has now d dt [ For a fixed time t and rapidity y one can characterize the correlations of the velocities in the transverse plane by (m, n = 1, 2) (6.14) Adapting to the standard notation used in the literature about turbulence, we write the Fourier transform of this as The tensor structure of (6.15) follows from rotational invariance and from the solenoidal constraint (5.3). The function E(t, k) depends on k 1 and k 2 only in the combination k = k 2 1 + k 2 2 . The normalization in (6.15) is chosen such that For a non-relativistic fluid, the function E(t, k) describes how the fluid kinetic energy per unit mass is distributed over the different wave vectors. We emphasize that in the relativistic setup considered here, E(t, k) is not directly representing kinetic energy. Instead it simply parameterizes the contribution to the fluctuating transverse velocity field from different wave vectors. The contribution of these fluctuations to kinetic energy, for example in the laboratory frame, can be determined but the resulting relation is more-complicated than in the non-relativistic case.
For the case of a relativistic heavy ion collision, the initial distribution E(t 0 , k) would characterize the relative strength with which different length scales 1/k are represented in the fluctuating initial conditions. The question of how this distribution of fluctuations evolves amounts then to studying the time-dependence of E(t, k). Here, we point only to one remarkable feature of the time-dependence of a two-dimensional fluid at large Reynolds number, that can be understood in 7 It has been argued that a cascade can take place for enstrophy 1 2 [ω 2 ] instead of energy 1 2 [ v 2 ], however. This is possible if [( ∇ω) 2 ] grows ∼ 1/ν for ν → 0 due to non-linear terms. the scaling theory of freely decaying turbulence in two dimensions developed by Batchelor [61]. This theory is based on the assumption that at sufficiently late times, the function E(t, k) remembers only a single number from its initialization, namely its average fluid velocity λ defined by λ 2 = 1 2 v 2 1 + v 2 2 . From dimensional reasoning it follows then that It follows from (6.16) that the dimensionless function h(x)is normalized to unity, ∞ 0 dx h(x) = 1. Interestingly, if one assumes that E(t, k) is dominated by the region around some characteristic wave vector k M , then this scale will change with time according to This implies that kinetic energy is shifted from small length scales to larger ones, in contrast to the Richardson cascade in three dimensions [63]. We would like to close this section on a cautious but speculative note: We recall first that Batchelor's theory of freely decaying turbulence was developed for very large Reynolds numbers that may not be realized in heavy ion collisions. However, the above considerations may make it conceivable that non-linear effects in the fluid dynamic evolution related to the onset of turbulence can move fluctuations in the initial kinetic energy to larger spatial scales. This phenomenon would be a distinct characteristics of an effectively two-dimensional turbulent evolution, and it may be identified experimentally by finding fluctuations related to length scales that are inconceivable to be present in fluctuating initial conditions. 7 The sensitivity of particle spectra on velocity correlation functions In the previous section 6, we have seen that the fluid dynamic evolution of fluctuations can be characterized efficiently in terms of correlation functions of fluid velocities. Here we discuss how information about such velocity correlations enters the particle spectra that are experimentally accessible in heavy ion collisions. The starting point of our discussion is the 'freeze-out' phase space distribution f (x, p) that parametrizes the matter distribution at the time τ fo when the particles decouple from the fluid dynamic evolution. In the following, we shall view f (x, p) as describing an event-averaged smooth fluid system supplemented by event-specific fluctuations. We shall then ask how these fluctuations are reflected in observables. This logic should apply to arbitrary choices of f (x, p). For the purpose of illustration, however, we shall restrict our discussion to a simple ansatz 8 that describes a locally approximately thermal distribution consistent with the Bjorken background field of section 3, and that allows for the implementation of local fluctuations on top of this background field. A simple choice with these properties is the Boltzmann distribution where the normalization d is fixed by the spin and flavor degeneracy of the degrees of freedom that decouple from the system. Assuming that the freeze-out takes place at some proper time τ fo when the average temperature drops below some freeze-out temperature T fo , the hadronic spectra can be calculated using the Cooper-Frye freeze-out prescription Here, the freeze-out volume is determined by p µ dΣ µ = m T τ cosh(η − y)dx 1 dx 2 dy for the case of Bjorken expansion. In practice, the spectra (7.2) measured in heavy ion collisions include averaging over many events. On the level of the freeze-out distribution f (x, p), this event average corresponds to an ensemble average with respect to the fluid velocity, energy density or temperature fields, respectively. Denoting the corresponding averaging by triangular brackets, we replace therefore in the calculation of (7.2) the function f by Denoting the particle four-momentum by (p 0 , p 1 , p 2 , p 3 ) = (m T cosh y, p 1 , p 2 , m T sinh y) with transverse mass squared m 2 T = p 2 T + m 2 , p 2 T = (p 1 ) 2 + (p 2 ) 2 , we expand f (x, p) for small fluctuations around the velocity profile of Bjorken (u 0 , u 1 , u 2 , u y ) = (cosh η, 0, 0, sinh η) and the constant freezeout temperature T fo The terms linear in T − T fo or u j vanish by definition or due to symmetry reasons and similar the cross-terms ∼ u j u i with i = j. Also, due to translational symmetry the variances at τ = τ fo are actually independent of the coordinates x 1 , x 2 and y, such that (u 1 ) 2 = G 11 u (τ fo ), (u y ) 2 = G yy u (τ fo ), (T − T fo ) 2 = G T (τ fo ). (7.5) This allows one to write the one-particle spectrum for small fluctuations around Bjorken flow as In the simple model studied here, information about velocity correlations enters the spectrum only via the coordinate-independent three numbers G 11 u (τ fo ), G yy u (τ fo ) and G T (τ fo ). For another model choice, the information may be slightly different. For instance, if one would replace the sharp freeze-out at τ fo by a decoupling at times τ around τ fo , then the three numbers G 11 u (τ fo ), G yy u (τ fo ) and G T (τ fo ) would be replaced by averages over τ . Irrespective of such model-dependent nuances, however, it is a generally known feature that the single particle spectrum (7.2) is only sensitive to space-time averages over the distribution f (x, p) and therefore does not contain information about correlations between different space time points.

Generalization to identical two-particle correlations
As discussed in section 6, the dependence of velocity correlations on the wave-numbers k 1 , k 2 and k y allows for a detailed characterization of fluid dynamic behavior, including information about the dissipation of fluctuations and the manifestations of turbulence. The one-particle spectra discussed so far contain only the information (7.5) about the correlations of fluid fields G ij u (τ, x 1 , x 2 , y), G j uT (τ, x 1 , x 2 , y), G T (τ, x 1 , x 2 , y) (7.8) at equal positions (x 1 = x 2 = y = 0). Here, we point out that identical (Bose-Einstein) two-particle correlation functions are linear functionals of (7.8) and may thus provide information about the wave number dependence of velocity correlations. Two-particle spectra for pairs of identical bosons (s B/F = 1) or fermions (s B/F = −1) of 4-momenta p A , p B respectively, can be written as [64] ). (7.9) Here, the first term is what one would expect from kinetic theory for classical particles while the second term results from the quantum statistics of identical particles. Data about two-particle spectra are typically normalized by a mixed-event technique that corresponds to forming a normalized correlation function (7.10) As it stands, equation (7.9) is valid for an event-specific realization of the velocity and temperature fields. Experimental data are for event samples that correspond to averages over the hydrodynamic fields. This amounts to replacing in equation (7.9) the products of phase-space distributions f by the corresponding event averages f (x, ) . Paralleling the arguments employed for the calculation of the one-particle spectra via (7.4), one should then expand the arguments of f (x, p A )f (x , p B ) for small fluctuations around the event-averaged background fields. In general, the arguments of these averages depend on space-time difference x − x , and they contain information about the relative position dependence of the correlations (7.8) in the fluid fields.
In practice, the way in which information about (7.8) enters the two-particle correlation functions may depend significantly on model-specific choices for f (x, p). For illustrative purposes, we explore here the particularly simple Bjorken-like model without spatial constraints in the transverse direction. We ignore all correlations involving the temperature field on the ground that these are for a compressionless situation formally of higher order in the fluctuating velocities, see Eq. (5.4). For the discussion in the following, we also neglect the rapidity-dependence of the correlation functions thus eliminating many terms proportional to G ij u with i = y or j = y (or both), as well as terms involving G y uT (The assumption of a vanishing rapidity dependence is relaxed in Appendix A). With these approximations, we concentrate therefore on the correlation functions of velocities G mn u with m, n = 1, 2. Due to rotational invariance in the transverse plane, the velocity correlation in Fourier space can be written in the form G mn u (τ, k 1 , k 2 , k y ) = 2π δ(k y ) [δ mn g 1 (τ, k) + k m k n g 2 (τ, k)] , (7.11) with k = k 2 1 + k 2 2 . In terms of the pair momentum P = 1 2 (p A + p B ) and the relative momentum q = p A − p B , the two-particle correlation function at mid-rapidity η A = η B = 0 takes then the form (see appendix A for details of the derivation) with P T = (P 1 , P 2 ), q T = (q 1 , q 2 ) and q T = q 2 T .
Realistic phase space distributions f (x, p) have support in a finite transverse area A T only, and this would lead to a fall-off of the correlation function to unity in the relative transverse momentum q T on a scale of order 1/ √ A T . Often, this fall-off is parametrized by a Gaussian ansatz in terms of HBT radius parameters, so that the first two terms of (7.12) would take the form C(P, q) . For the simplified model of infinite transverse extension discussed here, this contribution is singular ∝ (2π) 2 δ (2) ( q T ) and we have written it in a formal way normalized to unity at q T = 0.
For the purpose of the following discussion, the transverse translational invariance of the present toy model presents the technical advantage that the q T -dependence of the correlation function (7.12) is solely dependent on the wave number dependence of velocity correlations. Effects that could confound the interpretation of the q T -dependence in practice, such as effects from a finite transverse geometry and from transverse (event-average) velocity gradients, are not included in the present model. Therefore, the following discussion allows us to illustrate how information about velocity correlations enters two-particle correlation functions, but it limits our discussion of how such information could be disentangled from other effects in a phenomenologically relevant scenario. Keeping this caveat in mind, we observe that the first term in (7.9) can be viewed as an incoherent superposition of single-particle spectra and therefore does not furnish information that is not yet contained in single-particle spectra. In contrast, the quantum-statistical second term ∝ s B/F furnishes novel information about the wave number dependence of g 1 (τ fo , q T ) and g 2 (τ fo , q T ).
We note as a curious aside that for a situation of fully developed two-dimensional turbulence, one can predict the form of g 1 (τ, k) and g 2 (τ, k) from the scaling theory of Kraichnan and Batchelor. In particular, comparing (6.14), (6.15) and (7.11) one can express g 1 and g 2 as a function of E(t, k). Moreover, from Batchelors scaling theory of freely decaying two-dimensional turbulence one finds then the following scaling in the inertial range g 1 (τ, k) = c τ 10/3 k 4 , g 2 (τ, k) = − c τ 10/3 k 6 , (7.13) with a common but unpredicted constant c that reflects the absolute scale of velocity fluctuations. For the correlations function (7.12), this results in an additive term with a very slow power-law fall-off of the form ∝ Interestingly, if the velocity correlations occurs on scales that are significantly smaller than the transverse extension √ A T of the system, then the slow power-law q Tdependence persists at relative momentum scales that are significantly larger than the typical scales 1/R set by HBT radius parameters. It is an exciting possibility that a measurement of a power-law 1/q n T -dependence in two-particle correlation functions may provide a characteristic signature for turbulent distributions in the fluid dynamically evolved velocity fields of a heavy ion collision. We caution that the model discussed here is a simplified one; also, for intermediate Reynolds numbers one expects corrections to the case of fully developed turbulence, see e.g. [62]. What may persist in a phenomenologically realistic scenario, however, is the general idea that velocity correlations on small scales l T induce two-particle correlations on large scales q T ∼ 1/l T , and that these correlations are expected to be governed by a power-law fall-off.

Discussion and Conclusion
Recent data analyses from RHIC and LHC have given support to arguments that soft hadron spectra may result from the fluid dynamic response to initial conditions with significant event-by-event fluctuations. Motivated by this suggestion, we have studied here how event-by-event fluctuations propagate on top of an event-averaged fluid dynamical background of Bjorken type. The choice of a Bjorken background field is a simplification that retains essential elements of the expected fluid dynamical evolution of relativistic heavy ion collisions. As shown in the present paper, it allows for a particularly explicit, partly analytical discussion of the propagation of fluctuations in an expanding fluid dynamic system. In particular, we have found for the case of a laminar evolution explicit expressions for the attenuation or amplification of all fluid dynamic modes over the time scale relevant in heavy ion collisions. And we have specified the general conditions for non-linear effects in the dynamics of fluctuations, finding in particular that the late time dynamics evolves towards an essentially two-dimensional system with turbulent behavior. Here we discuss our main findings in more detail: To discuss the propagation of fluctuations, one needs to specify first the nature of the fluctuations that are propagated. Recent studies have focussed mainly on fluctuations in the energy density (or, equivalently, entropy density), where the Glauber model provides a phenomenologically supported basis for assuming local fluctuations on a particular transverse scale. However, it had been pointed out already that fluctuations in the velocity field may be present in the initial conditions for fluid dynamic evolution. Velocity fluctuations could arise from pre-equilibrium evolution, or they could be a natural consequence of fluctuations in the primary interactions of elementary constituents. To the best of our knowledge, there is no a priori argument that fluctuations in energy density dominate over fluctuations in other fluid dynamic fields. Also, it requires studies allowing for all possible fluctuations to address the question whether and how fluctuations in velocity and energy density can be disentangled. In section 2, we argued that a discussion of the fluid dynamic response to fluctuating initial conditions should be based on a formulation that allows for fluctuations in all fluid dynamic fields. In particular, we supported with a model study the idea that fluctuations in the velocity fields may carry significant vorticity. For a fluid with conserved charges such as baryon number and electric charge one should take fluctuations in these quantities into account, as well.
In general, a separation of fluid dynamical fields into background and fluctuations is not necessary. For instance, recent studies of event-by-event fluid dynamics propagate event samples of initial conditions numerically without separating fluctuations from background fields. In principle, such full fluid dynamical simulations allow to explore under the most versatile model assumptions the fluid dynamic response to fluctuating initial conditions. But an explicit mode-by-mode formulation of the dynamics of fluctuations around a fluid background field seems well-suited to study which fluctuating modes in the initial conditions can survive the strong dynamical evolution in a heavy ion collision unattenuated, whether there are mechanisms that may amplify some modes, and which modes are 'filtered out' by the medium due to dissipative effects in the fluid dynamic evolution. To address such questions in a simplified framework that accounts for the main features of fluid ex-pansion, we have formulated in section 3 the fluid dynamical evolution of local fluctuations around average fluid fields of Bjorken type. For this system, we have found a peculiar rescaling of the time variable, t ∝ τ 4/3 , that allowed us to write the relativistic fluid dynamic evolution of fluctuations in a form resembling a non-relativistic Navier-Stokes equation.
If the Reynolds number of a fluid system is not too large, then important elements of the dynamics may be understood in a linearized treatment. In section 4, we observed that in this laminar case, fluctuations in vorticity decouple from sound modes and fluctuations in energy density. For a parameter set that is characteristic for heavy ion collisions, we have then studied how different modes are amplified or attenuated over the time scale of order 10 fm/c of a heavy ion collision. Remarkably, the longitudinal vorticity modes have an algebraic enhancement factor that can overcompensate on this time scale the typical exponential decay due to dissipative effects. Therefore, vorticity, if not present in the initial conditions, will not be generated as long as the dynamical evolution is laminar, and it is therefore unlikely to be generated in a sizable amount for small Reynolds numbers. However, if present in the initial conditions, some vorticity modes will be amplified significantly during the dynamical evolution. To the best of our knowledge, this phenomenon has not been studied yet in numerical simulations, and it would be interesting to see how it manifests itself in the presence of other background fields. In addition to the vorticity modes, we have also explored the propagation of sound modes that result from initial fluctuations. In general, we find that fluctuations of sufficiently long wave-length pass unattenuated over time scales relevant for heavy ion collisions, while short wave-lengths that reflect finer structures in the fluctuating initial conditions, are dissipated on shorter length scales. There is also a characteristic difference between transverse components, and the components that propagate in the longitudinal direction in which the system expands according to Bjorken's model.
Outside the regime of validity of a linearized treatment, the discussion of solutions of fluid dynamics is very complicated and typically requires numerical techniques. For our problem of fluctuations around Bjorken flow, we are in the special and fortunate case that we can relate the full relativistic dynamics of fluctuations in rescaled coordinates to a non-relativistic Navier-Stokes equation. This allows us to discuss the possibility of a turbulent evolution in terms of those concepts and parametric estimates that have proven useful in characterizing turbulent phenomena of non-relativistic systems. In section 5, we introduce both a longitudinal and a transverse Reynolds number to characterize the non-linear dynamics of fluctuations on top of a Bjorken background field that shows strong dynamical expansion only in the longitudinal direction. We find in particular that the late time dynamics will evolve a large set of initial conditions into a regime where the dynamical evolution is effectively two-dimensional, and where the transverse Reynolds number can be sizable. This indicates a window for a two-dimensional, non-linear evolution towards turbulent behavior. Motivated by this observation, we have summarized in section 6 characteristic features of turbulent behavior, detailing the differences between the cases of two-dimensional and threedimensional evolution. We recall from this discussion in particular that in three dimensions, fluid kinetic energy thermalizes typically by dissipating into vorticity modes of increasing wave number, i.e. decreasing length scales. As first observed by Kraichnan, this mechanism is not possible for a two-dimensional fluid system where one finds an inverse cascade: kinetic energy gets propagated to larger length scales.
The transverse Reynolds numbers estimated in section 5 do not support the assumption that heavy ion collisions create systems with fully developed turbulence. However, our discussion in section 5 shows that realistic Reynolds numbers are not small enough to neglect non-linear effects in the dynamical evolution. Therefore, while we have no rationale to expect that the correlation functions of fluid dynamic fields generated in heavy ion collisions satisfy the scaling laws of Kolmogorov's theory of fully developed turbulence, we do expect that a non-linear dynamics that can be regarded as the onset of turbulent behavior may result in interesting (power-law) wave-number dependences of correlations of fluid fields. By supplementing a standard blast-wave model with event-by-event fluctuations in fluid fields, we have established in Sect. 7 how such fluid dynamic correlation functions manifest themselves in one-and identical two-particle spectra. For one-particle spectra, we find that fluctuations are a confounding factor in interpreting the m T -dependence of spectra in a fluid dynamic scenario. For two-particle correlations, we observe a dependence that may provide an experimentally accessible signature for the onset of turbulence in heavy ion collisions. The observation is that identical two-particle correlations at large relative momentum are sensitive to spatial scales that are much smaller than the transverse size of the particle producing source. If two-point velocity correlations show turbulent behavior on these small scales then this will translate into a characteristic power-law tail of identical two-particle correlations at large relative momentum.
Let us close with a short outlook of how observations made in this study could be pursued further. Our study was largely motivated by the question of how initial conditions with significant velocity fluctuations (in the solenoidal an in the irrotational part) propagate fluid dynamically. Here, full fluid dynamic simulations including initial velocity fluctuations could provide further insight, for instance by evolving fluctuations around average fluid fields with more realistic transverse dependencies, and by quantifying to what extent initial velocity fluctuations could contribute to the observed azimuthal asymmetries in momentum space. Full fluid dynamic simulations including initial velocity fluctuations could also allow for a detailed characterization of how the scale dependence of fluid correlation functions of the type (7.8) builds up during the fluid dynamic evolution. This would provide in particular insight into the question on which time scales and over which range of wave vectors correlation functions of fluid fields may develop power-law dependences that can be regarded as precursors of turbulent phenomena 9 . Also, the semi-analytical approach used in the present work may be pursued further. In the present work, we have seen how single vorticity modes, and sound modes are amplified or filtered out by the dynamical evolution, depending on their wavenumber. We plan to investigate, whether a similar understanding can be gained for a more realistic background field by expanding fluctuating fluid fields in terms of appropriate sets of functions so that one can calculate explicitly how the fluid dynamic evolution mixes different components in the evolution. We expect that such studies could provide an intuitive understanding e.g. of the time scales on which density fluctuations feed sound waves, or on which vorticity modes cascade to other scales. This may help significantly in the interpretation of the complex fluid dynamic phenomena that we expect to find realized in heavy ion collisions.
For k y = 0, this expression reduces to the first two terms in the first line of (7.12).
Let us now consider the second term ∼ s B/F in Eq. (7.9). This term is of the form of a single integral over the freeze-out volume times its complex conjugate. Expanding it up to terms that are quadratic in the fluctuations, one obtains two distinct contributions, In one case, the bilinear terms in the fluctuating fields are written at different points x, x , in the other case they are both taken at the same point. We consider both cases separately. We work in the following with average pair momentum P µ = 1 2 (p µ A + p µ B ) = (m AB T cosh η AB , P 1 , P 2 , m AB T sinh η AB ) and we define a transverse mass defined m AB T = (P 0 ) 2 − (P 3 ) 2 and a rapidity η AB = arctanh(P 3 /P 0 ). This allows us to write q µ x µ = −[m A T cosh(η A − y) − m B T cosh(η B − y)]τ + q T x T with q T = (q 1 , q 2 ) and x T = (x 1 , x 2 ). The contribution to the second term in Eq. (7.9) that is quadratic in fluctuations at the same space-time point, can then be written as Here, the integral over the transverse coordinates leads to a factor (2π) 2 δ (2) ( q T ) (where (2π) 2 δ (2) (0) = A T is understood). This is a consequence of the fact that in the present model the transverse extension of the particle emitting source is not limited. At mid-rapidity, η A = η B = η AB = 0, the term (A.6) simplifies to d τ fo m AB T (2π) 3 (2π) 2 δ (2) ( q T ) 1 + With the approximations used in section 7, this expression reduces to the last terms in the first line of (7.12). We now turn to the contribution ∼ s B/F where the fluctuating fields have different space-time argument. This term is of similar structure as (A.2) and reads