Relativistic Hydrodynamic Fluctuations

We present a general systematic formalism for describing dynamics of fluctuations in an arbitrary relativistic hydrodynamic flow, including their feedback (known as long-time hydrodynamic tails). The fluctuations are described by two-point equal-time correlation functions. We introduce a definition of equal time in a situation where the local rest frame is determined by the local flow velocity, and a method of taking derivatives and Wigner transforms of such equal-time correlation functions, which we call confluent. We find that the equations for confluent Wigner functions not only resemble kinetic equations, but that the kinetic equation for phonons propagating on an arbitrary background nontrivially matches the equations for Wigner functions, including relativistic inertial and Coriolis forces due to acceleration and vorticity of the flow. We also describe the procedure of renormalization of short-distance singularities which eliminates cutoff dependence, allowing efficient numerical implementation of these equations.


A. Motivation and overview
Hydrodynamics -the universal theory describing macroscopic motion of fluids -hardly needs an introduction. The range of its applications is extraordinarily wide -from molecular biology to astrophysics. The well-established conceptual text-book framework of hydrodynamics [1] has received considerable renewed attention and further development from various points of view recently. One of the drivers of the recent interest is the necessity to develop tools for quantitative analysis of heavy-ion collisions [2,3]. A major ingredient which is needed is relativistic hydrodynamics with fluctuations. In many common contexts fluctuations in hydrodynamics could be considered negligible, such as in truly macroscopic systems with O(10 24 ) particle degrees of freedom. Heavy-ion collisions, however, occupy a "sweet spot" in terms of the size: with O(10 2−4 ) particle degrees of freedom the relevant systems are large enough to be treated hydrodynamically but small enough for fluctuations to be important and directly observable via event-by-event measurements. In particular, fluctuations are expected to be further enhanced if the matter created in the collisions is in a state close to a critical point. In this case fluctuations can serve as signatures of the critical point [4][5][6][7] in the beam energy scan experiments [8].
In the classic Landau-Lifshitz [9] approach to hydrodynamic fluctuations the local noise due to microscopic degrees of freedom is introduced into constitutive equations. Generalizing this formalism to relativistic hydrodynamics and applying it to relativistically expanding solutions is one of the approaches pursued in recent literature [10][11][12]. The main drawback of this approach is that practical implementation (e.g., for realistic heavy-ion collision simulations) requires introducing local noise whose amplitude needs to be taken to infinity as the coarse-graining distance scale (hydrodynamic cell size) is sent to zero. Nonlinearities lead to divergent noise-induced corrections to equation of state as well as transport coefficients and make numerical simulations difficult if not outright infeasible.
An alternative way of describing dynamical effects of fluctuations was introduced by Andreev in the 1970s who considered evolution of two-point equal-time correlation functions [13]. This approach has the advantage of being formulated in terms of deterministic equations, avoiding the "infinite noise" problem in the implementation of the stochastic approach. More precisely, the effects of the "infinite noise" can be isolated and absorbed into "renormalization" of the equation of state and transport coefficients in a close analogy with the renormalization in quantum field theories. Of course, the stochastic and the correlation function (deterministic) approaches are equivalent and complementary, in ways very similar to Langevin and Fokker-Plank description of stochastic processes, and the ultimate choice is to be made based on practicality, in particular, for numerical simulations. The deterministic approach (also referred to as hydro-kinetic approach due to the similarity of some of the additional equations to kinetic equations for phonons 1 ) has been recently discussed in a relativistic context [14][15][16] for a special case of Bjorken boost invariant solution where symmetries allow to reduce the effective dimensionality of the problem and simplify the analysis.
A more general approach is needed in order to lift the limitations of the static or boost invariant solution and to enable practical simulations of relativistic hydrodynamics with fluctuations in a general inhomogeneous three-dimensional background characteristic of heavy-ion collisions. Such an approach should, for example, capture the effects of vorticity, absent in the Bjorken solution, but important in heavy-ion collisions [17]. The aim of this paper is to develop such a universal approach.

B. Variables and scales
Hydrodynamic variables are macroscopically averaged values of densities of conserved quantities, such as energy and momentum. The macroscopic averaging is done at fixed time t over a region of linear size b (hydrodynamic cell) around a point with spatial coordinates x. In order to be macroscopic, the length b must greatly exceed microscopic scales, mic , such as mean-free path in a weakly-coupled system, or thermal length 1/T in a relativistic strongly coupled system b mic . (1.1) The resulting coarse-grained variables can be used to describe evolution of inhomogeneities at larger length scales To facilitate the discussion let us refer to the hydrodynamic variables defined via coarse-graining discussed above asψ A (t, x), where index A labels a variable. Since we are describing a thermal system, the variablesψ A are stochastic -fluctuating between members of the statistical ensemble describing our system (in heavyion collisions -between collision events). To be more precise, variablesψ A are operators. However, due to macroscopic averaging involved in their construction they behave as classical (commuting) stochastic variables. Their quantum fluctuations are negligible compared to (classical) thermal fluctuations. 2 Due to coarse graining, fluctuations at scales shorter than b are averaged out, i.e., suppressed. In this sense, Λ = 1/b plays the role of the ultraviolet (wave-vector) cutoff.
In order to describe fluctuations we introduce the ensemble averages of the variables ψ A ≡ ψ A . The ensemble averages ψ A obey deterministic hydrodynamic equations. In addition to these usual hydrodynamic variables (one-point functions) we must introduce two-point functions which are ensemble averaged equaltime products at two space-time points: φ A (t, x 1 )φ B (t, x 2 ) , where φ =ψ − ψ is the fluctuating part of the variable, as usual. 3 In equilibrium, the correlators φ A (t, x 1 )φ B (t, x 2 ) are translationally invariant, i.e., independent of the midpoint x ≡ (x 1 + x 2 )/2 at fixed separation y ≡ x 1 − x 2 . The dependence of correlation functions of operators in equilibrium on separation y is characterized by exponential fall-off at distances larger than correlation length ξ: e −|x1−x2|/ξ . Correlation length ξ is a microscopic scale, typically: ξ ∼ mic . 4 From the point of view of the coarse-grained variables φ A , therefore, the distance ξ b is negligible and the equilibrium correlator φ A (t, x 1 )φ B (t, x 2 ) is essentially a multiple of the delta-function δ 3 (x 1 − x 2 ).
Hydrodynamics, however, describes systems which are not in complete equilibrium: variations, or gradients, of the variables over macroscopic scales L lead to evolution (flow) characterized by time scale τ ev ∼ L/c s , 1 Despite this similarity to kinetic theory, the correlation function approach does not rely on validity of any underlying microscopic kinetic description. As hydrodynamics itself, the approach is applicable for either weakly or strongly coupled quantum field theories. The quasiparticles described by "hydro-kinetic" equations are macroscopic hydrodynamic excitations, such as phonons. 2 The precise condition for that is that the quantum uncertainty of the energy due to finite characteristic time of the evolution of these variables is much smaller than their typical thermal energy, T . The fastest evolving degrees of freedom after coarse graining are sound modes with wave-length b. Their frequency cs/b must therefore be much smaller than T , i.e., b cs/T . 3 More generally, the variable φ A could be (and will be) a linear combination ofψ B − ψ B . 4 Near a critical point correlation length is large, e.g., ξ 1/T , and additional hierarchy of scales emerges. In this case our analysis applies if the macroscopic size b is taken to be much greater than ξ: b ξ. Dynamics of fluctuations near a critical point in the opposite regime, b ξ, is characterized by dynamical scaling and is discussed in [18].
FIG. 1. Schematic illustration of various scales described in the text. The scale L of the variation of the background ψ(x) is the longest in the problem. In equilibrium, the fluctuation correlator G(x, y) = φ(x + y/2)φ(x − y/2) becomes a function (illustrated by a sharp peak on the figure) whose width in y is narrower than the shortest hydrodynamic scale -the coarse graining scale b. If the system evolves, the correlations at scales eq ∼ √ τev ∼ √ L are not yet completely vanishing, giving G(x, y) a finite width of order eq L. A negative contribution from additional correlations is necessary to satisfy conservation laws x φ(x) = 0, i.e., y G(x, y) = 0 (this integration does not commute with τev → ∞, i.e., equilibrium limit, in which the integral is equal to a susceptibility).
where c s is the sound speed. The (re)equilibration, as the system evolves, requires transport of conserved quantities, which is a diffusive process. Therefore equilibrium can be established only over scales which can be reached by diffusion over time of order τ ev : where γ is an appropriate diffusion constant (typically, γ ∼ 1/T ). 5 This means for distances b |y| eq or, more precisely, for wave-vectors q conjugate to y such that 1/ eq |q| Λ the equilibration is complete. However, at scales around q ∼ 1/ eq ∼ c s k/γ, where k ∼ 1/L, the equilibration is ongoing, as it is trying to catch up with the evolution of the system. It is this competition between the equilibration and evolution that we will be describing.
Since the relevant values of y ∼ eq are parametrically shorter than L, we can consider the equilibration process as local, occurring on a slowly varying background set by local values of ψ A (t, x). For that reason it is also naturally convenient to use the mixed Fourier (i.e., Wigner) transform of the correlation function φ A (t, x + y/2)φ B (t, x − y/2) ≡ G AB (x, y) with respect to separation vector y, which we shall denote W AB (x, q). The relevant values of q will satisfy k, γq 2 /c s q Λ T , (1.5) where, for simplicity, we took −1 mic ∼ T . The condition k q allows us to treat background as smooth when describing the relaxation of correlations W AB (x, q) to equilibrium. However, we must retain non-zero gradients, ∂ µ ψ A , (proportional to k) of the background variables in the equations for Wigner functions W AB , since those gradients drive the deviations of W AB from equilibrium.
The fluctuations described by W AB , in turn, feed back into constitutive equations which determine the evolution of the background flow. One must, therefore, solve the equations for the background flow together with the equations for W AB that we are going to derive in this work.
Fluctuations with all wave vectors q up to the cutoff Λ contribute to this feedback. The integral of the contributions over q is divergent, i.e., it depends polynomially on the cutoff Λ. This would cause difficulties in numerical implementation of the hydrodynamic equations and is the manifestation of the "infinite noise" problem. In a remarkable similarity to the renormalizaton of wave-functions and couplings in quantum field theories, the fluctuations in hydrodynamics renormalize variables (energy density and flow velocity) and parameters (equation of state and transport coefficients). The hydrodynamic renormalization absorbs the leading large-q terms in W (x, q) responsible for divergences and thus removes the polynomial dependence on the cutoff Λ, allowing efficient numerical implementation.
Once these cutoff-dependent contributions are absorbed into the "renormalized" hydrodynamics, the true (observable) feedback of the out-of-equilibirum fluctuations comes predominantly from modes with wave vectors q ∼ 1/ eq . Most importantly, it is finite and cutoff independent. The magnitude of these nonequilibrium effects can be estimated as the phase-space volume d 3 q ∼ −3 eq ∼ (c s k/γ) 3/2 . The power 3/2 indicates that these effects are non-local. They are known as "long time tails" of hydrodynamic response [13,[19][20][21]. Their contribution is typically more important than that of the second order O(k 2 ) terms in the hydrodynamic derivative expansion (unless suppressed by a microscopic parameter, such as e.g., number of colors in a gauge theory) 6 .
The discussion of the correlation function above has glossed over an important issue: "equal time" in the definition of G AB implies a certain choice of the frame with respect to which equality of time, i.e., simultaneity, is to be determined. This problem does not arise in non-relativistic hydrodynamics, but in the case of heavy-ion collisions it is essential, since the relative velocities at different points in the fireball are comparable to the speed of light. If the fluid moves as a whole, with the same velocity (in the lab frame), the rest frame of such a fluid, not the lab frame, is the natural choice. In the cases of interest, such as relativistically expanding fluid, the local rest frame of the fluid is a function of space and time. We can describe it, as usual, by the 4-velocity u(x) (macroscopically averaged as described above). Therefore, to define the equal-time correlation function we consider correlator and evaluate it at points where 4-vector y lies in the hyperplane orthogonal to u(x): u(x) · y = 0. 7 The corresponding wave vector q in W AB (x, q) also resides in a hyperplane orthogonal to u(x), which is x-dependent. We find it useful to introduce a type of space-time derivatives which account for this x-dependence due to inhomogeneous flow and which we call "confluent" derivatives. The paper is organized as follows: In Section II we introduce stochastic hydrodynamics and expand its constitutive equations up to quadratic order in fluctuations around an arbitrary background. We use linearized hydrodynamic equations for fluctuations to derive equations obeyed by two-point correlators. In Section III we introduce confluent derivative and Wigner function which allow us to write the equations obeyed by "equal-time" correlators. These equations are presented and studied in Section IV.
In Section IV B we observe that some components of W AB oscillate at frequencies of order c s q, which are faster than the evolution of the background and thus, for most practical purposes, can be averaged out by introducing additional temporal coarse-graining scale b t 1/(c s q). The equations for remaining, slower components simplify.
In Section V we consider in detail the fluctuation contributions due to nonlinearities, and review a general procedure of renormalization of first order hydrodynamics. We study the asymptotic behavior of W AB at large q and identify the parts of W AB that lead to renormalization of the equation of state and the transport coefficients.
In Section VI we obtain equations of motion for a phonon in a non-trivial flow using variational principle, find the corresponding kinetic Liouville operator, and show that it exactly matches, in several nontrivial ways, the kinetic equation derived in Section IV B.
Several Appendices contain useful supplementary information. In particular, we assemble a list of our notation choices used throughout the paper in Appendix D.

A. Stochastic hydrodynamics
Hydrodynamic equations express conservation and transport of energy and momentum densities: 8 To simplify notations later in the paper we label fluctuating hydrodynamic quantities with an accent, as in T µν , to distingish them, where necessary, from quantities which are not fluctuating. The four conservation equations (2.1) are solved for the same number of hydrodynamic variables. A convenient covariant choice for them is the fluid velocityȗ µ (normalized asȗ ·ȗ = −1) and the energy density˘ in the rest frame of the fluid, that are defined by the Landau's matching condition To form a closed system, we need six additional (constitutive) equations to express all components of T µν in terms of and u µ . For macroscopically large scale dynamics of hydrodynamic variables, T µν can be expanded in gradients of and u µ . The first-order (Landau-Lifshitz) hydrodynamics corresponds to truncating this expansion at first order in gradients: where p( ) is pressure as a function of -also known as the equation of state and w( ) = + p( ) is the enthalpy. The viscous tensor is linear in gradients of u: where shear and bulk viscosities are denoted as η and ζ, respectively, and is the projection operator to the spatial hypersurface orthogonal to u, in terms of which we define: However, the constitutive equations (2.3) relating T µν and the hydrodynamic variables are valid only on average, and there exist random local thermal noiseS µν which makes Eq. (2.1) a stochastic differential equation withT The functions of hydrodynamic variables such as w, p, Π µν , etc. in Eq. (2.8) are the same as in Eq. (2.3) and Eq. (2.4) but they are evaluated for fluctuating variablesȗ and˘ . The hydrodynamic variables in (2.4) fluctuate as they are driven by the random noiseS µν , and we need to consider statistical ensemble average over these fluctuations for any observables on macroscopic scales. We write our stochastic hydrodynamic variablesȗ µ and˘ as a sum of their averages, u ≡ ȗ , ≡ ˘ , and linear fluctuations around them as:ȗ = u + δu,˘ = + δ . These fluctuations are driven by the noise termS µν with S µν (x) = 0, the strength of which is set by the fluctuation-dissipation theorem 9 , In principle, it is possible to numerically solve the stochastic equation ∂ µT µν = 0 with some coarse-graining, or wave vector cutoff Λ, which regularizes the infinite amplitude of the noise arising from the δ (4) (x − x ) term. However, as we already mentioned in the introduction, the results would depend sensitively on the cutoff Λ due to non-linearity of hydrodynamic equations.
We follow an alternative approach, that is, we include fluctuation contributions to T µν by expandingT µν to second order in fluctuations. The fluctuation contributions to T µν are given by two-point correlators of the fluctuations, and to describe their evolution we derive a separate set of equations. After proper renormalization that absorbs cutoff dependence into physical parameters, the equation of motion ∂ µ T µν = 0 along with the equations for the two-point functions defines a deterministic coupled time evolution of the averaged variables and of the correlation functions that can be solved numerically.
Due to non-linearities in the relation between the variables ( , u) and T µν in Eq. (2.3), including nonlinearities in the equation of state, such as where c 2 s = dp( )/d is the square of sound speed,T µν in Eq. (2.8) contains terms that are nonlinear in fluctuations. ExpandingT µν up to second order in fluctuations and taking the average, we have: Note that we neglected the fluctuations of the viscous part Π µν , which are parametrically smaller than the terms kept in the above expansion 10 . In the last line we introduced the collective notation for the fluctuating modes, δ and δu µ : where the scalar δeδe , vector δeδg µ , and tensor δg µ δg ν components of the two-point correlation function are expressed compactly as where A ∈ (e, 0, 1, 2, 3) 11 . In terms of our definition of the correlator G AB (x, y) in Eq. (1.6), We can express the fluid velocity in the collective notation as well: Due to the presence of gradients, the system is slightly out of equilibrium and the fluctuation-dissipation relation given in Eq. (2.11) contains corrections proportional to the gradients. However the effects of these corrections are higher order (in k/q) in the fluctuation expansion as well as the kinetic equation that we discuss in this paper. Therefore we can safely use relation Eq. (2.11) with T and w being functions of x in the remainder of the paper. 10 We rely on γq ∼ q/T 1, according to Eq. (1.5), where q is the typical wave vector of the fluctuations. 11 The mixed index A ∈ (e, 0, 1, 2, 3) is raised and lowered by the "metric", diag(1, −1, 1, 1, 1). However the object u A is not a vector, rather an array that conveniently combines scalar and vector modes.
It should be noted that not all five variables φ A are independent since, due to normalizationȗ ·ȗ = −1, we have a constraint u A φ A = 0. Correspondingly, where s = w/T is the average entropy density. The functions G AB (x) in Eq. (2.13) are, in general, non-local functionals of the background fields and u. In the next section, we will derive the evolution equation for them by using the linearized hydrodynamics equation of motion for fluctuations.

B. Linearized stochastic equations for fluctuations
In this section we derive the stochastic equation that governs the dynamics of the linearized fluctuations of the hydrodynamic modes, δe and δg µ . This equation is the building block for the evolution equation for the two-point function G AB (x, y) and its Wigner transform, which we call "kinetic equation". The energy-momentum tensor expanded to linear order in fluctuations is given by where γ η = η/w and γ ζ = ζ/w. (2.21) In this expansion the first two terms are zero'th order in gradients and the third term Π µν is of first order or, equivalently, of order k. These three terms constitute the average background value without fluctuation contributions, i.e., T µν in Eq. (2.3). The remaining terms are linear in fluctuations. We consistently neglected several terms (e.g., fluctuations of viscosities) that are suppressed by either a factor of k/q 1 or γk ∼ k/T 1 compared to the terms being kept, according to our hierarchy of scales in Eq. (1.5) (recall that k is the scale of background gradients, and q is the wave-vector of fluctuations).
The stochastic equation for the linearized modes follows from the energy momentum conservation ∂ µT µν = 0, where we also neglect several terms based on similar considerations discussed above. By averaging both sides we obtain ∂ µ T µν = 0 to leading order in fluctuation expansion, which we insert back into Eq. (2.22) to arrive at a stochastic differential equation for the linearized fluctuations. In terms of the notation φ A introduced in Eq. (2.14), the equation for the linearized fluctuations reads: where L, Q, and K are 5 × 5 matrix operators. The operators L and Q are the ideal and dissipative terms, respectively, K contains the corrections due to the first-order gradients of background flow, and ξ denotes the random noise. Explicitly 12 where a µ = u · ∂u µ is the fluid acceleration. Note thatċ s terms arise due to space-time variation of c s via its dependence on (x).

C. Equations of motion for the fluctuation correlators
Having equipped ourselves with the equations of motion for the linearized fluctuations, we now derive the evolution equation for G AB (x, y) with respect to the midpoint variable x at fixed y. This equation will be used in Section IV to eventually obtain the evolution equation for its Wigner transform W AB (x, q). Using the definition of G AB (x, y) in Eqs. (1.6) and (1.7) and noting that we apply Eq. (2.23) and keep only the leading order in derivative expansion (i.e., retaining terms of order γq 2 or k, but not qk, consistently with Eq. (1.5)), to obtain where we converted the independent space-time variables from (x + , x − ) to (x, y) by using Eq. (1.7) and used superscripts '(y)' on the operators to specify that the derivatives involved are to be taken with respect to y at fixed x. In particular, the operator comes from the conversion of x ± derivatives into x, y derivatives, and collects terms proportional to y, which resulted from the y-dependence of u(x ± ) and c s (x ± ). The last term in Eq. (2.26) follows from the usual procedure of stochastic calculus, by keeping the random noise two-point function in double integrals over the time interval δt and using the correlation given in Eq. (2.11).
In order to convert this equation into an equation for the Wigner transform W AB of the correlator G AB we need to define the Wigner transform more carefully than was necessary until now. We also find it necessary to introduce a concept of derivative adjusted for the boost by flow, which one can call "flow-adjusted derivative" or "confluent derivative". The concept of frame transformation (boost) involved in its definition bears some resemblance to the parallel transport in differential geometry and the derivative itself is similar to covariant derivative.

III. CONFLUENT DERIVATIVE, CONNECTION AND WIGNER FUNCTION
In this section we discuss several ingredients which we need to translate equation (2.26) into an equation for the Wigner function. We begin by discussing how to take derivative of an equal-time correlator in a situation where the concept of equal time is different in different space-time points.
In Eq. (1.6) we defined the equal-time correlator of hydrodynamic variables as a function of the mid-point x and the separation vector y as where the domain of y is the 3-dimensional plane orthogonal to u(x), i.e., y is purely spatial in the local rest frame at x. We want to define a partial x derivative of such a function at "fixed" y. This is not straightforward, as the following expression illustrates: In G(x + ∆x, y) the orthogonality condition u(x + ∆x) · y = 0 is, in general, false, given u(x) · y = 0 is true: vector y spatial in the frame u(x) is not spatial in u(x + ∆x) (see Fig. 2). To preserve the relationship between u and y we need to transform vector y by the same boost that takes u(x) to u(x + ∆x). Defining this boost as Λ −1 (∆x) (inverse for later convenience), i.e.: 13 we can then define a derivative at "fixed" y as We have so far suppressed indices A and B in G AB which label hydrodynamic variables being correlated. These variables transform covariantly under Lorentz boosts (five components of φ A contain a scalar and a 4-vector according to Eq. (2.14)). It is natural to define a derivative which measures the changes of the hydrodynamic variables with respect to the local rest frame defined by flow velocity u. I.e., we are not interested in the changes between φ A (x + ∆x) and φ A (x) which are simply due to the difference in the local velocity u, i.e., induced by boost transformation from frame u(x) to u(x + ∆x). In other words, we are interested in the "internal" state of the variables, not affected by frame choice. The corresponding derivative could be constructed by boosting the variable φ(x + ∆x) in the same way as u in Eq. (3.3) before comparing to φ(x), i.e., 14 With respect to such a derivative, by construction, the flow vector field u(x) is "constant": according to Eq. (3.3). We shall refer to such a derivative as "confluent" to distinguish it from a common covariant derivative. Using explicit form of the infinitesimal boost defined by Eq. (3.3): where ∆u = u(x + ∆x) − u(x), we obtain the explicit expression for the derivative: where the connection associated with the boost created by flow gradients is given bȳ Note that this connection is antisymmetric with respect to µν, reminiscent of a spin connection. In a sense, it is a spin connection for a tangent space spanned by hydrodynamic variables φ A at point x. In that sense confluent derivative is a covariant derivative for the connection given by flow gradients in Eq. (3.9). To unify equations we can extend the range of indices to accommodate the full 5-dimensional space of variables and write∇ including the case when A or B is e. The corresponding connection is, of course, zero, since φ e = c s δ is a scalar.
Following the same logic that led us to Eq. (3.5), we would also like to eliminate the dependence of the correlator on the difference of the flow velocities between points x + = x + y/2 and x − = x − y/2. Therefore, we define a confluent correlation function by boosting both variables φ A (x + y/2) and φ B (x − y/2) into the rest frame at the midpoint, x, i.e, As a result, the confluent correlator, in contrast to Eq. (2.18), satisfies a simpler orthogonality condition: Now combining the three ingredients given by Eqs. (3.4), (3.8) and (3.11) we define the confluent derivative in the following way: This expression may be more useful for numerical integration of equations we derive, where derivatives need to be discretized. The expression which is used in analytical manipulations is obtained by Taylor expanding in ∆x:∇ (3.14) Another connection,ω b µa (a, b = 1, 2, 3), appears because we need to define a tangent space at each point x and introduce coordinates, such as y a , in this space to describe vector y and to keep them fixed, as we take x derivative (and to take derivatives with respect to y a at fixed x). To do this we choose an arbitrary local basis triad, e µ a (x), at each point x, such that u · e a = 0. When we keep vector y "fixed" (in the sense described above), its local coordinates y a = e a µ (x)y µ may still change due to the rotation of the triad, i.e., not only e a (x + ∆x) = e a (x), but, in general, also (Λ(∆x)e) a (x + ∆x) = e a (x), in contrast to Eq. (3.3) (see Appendix A). The last term in Eq. (3.14) makes sure this change of the basis is corrected for. In other words, we need additional connection to make e µ a confluently constant (like u already is without additional connection), i.e.,∇ λ e µ a ≡ ∂ λ e µ a +ω µ λν e ν a −ω c λa e µ c = 0 , (3.15) so that∇ λ y =∇ λ (e µ a (x)y a ) = 0. The second term in Eq. (3.15) accounts for the boost of e a as the one in Eq. (3.8) and illustrated in Fig. 2. In other words, and this is important for applications, the x-deriviative in ∂ µḠAB in Eq. (3.14) is taken at fixed y a and the boost needed to keep y = e a y a orthogonal to u is taken care of by e a . The last term in Eq. (3.15) describes the possible additional rotation of the triad basis vector in the tangent hyperplane. This rotation depends on our (arbitrary) choice of the local triad e a (x).
Equation (3.15) can be solved for the connectionω a λb by multiplying by dual basis vector e b µ such that where we used u · e b = u · e a = 0 and the definition ofω connection in Eq. (3.9). In Appendix A we provide a simple explicit example of the local triad e a with the corresponding connection.
We can now define the Wigner transform of the equal-time correlatorḠ AB (x, y) by integrating over the 3d hyperplane normal to u(x) at each point x. The integral can be expressed explicitly as the integral over coordinates y a , in which form it can be practically evaluated in numerical applications, or, more formally, as an integral over y constrained to a plane by u · y = 0 condition, i.e., ( 3.17) Thus we arrive at the definition of the confluent Wigner function: Note that due to the delta-function constraint the Wigner function W AB (x, q) does not depend on the component of q along u (energy/frequency in local rest frame). Therefore, we only need 3 independent components for vector q. We shall use the triad basis we already introduced above for vector y (see also Appendix A) to express 4-vector q µ in terms of its 3 independent components q a as q µ = e a µ q a . It is now straighforward to write the expresion for the confluent x derivative of the Wigner function at q fixed. We need to use the rules of the Fourier transform to replace y a → i∂/∂q a and ∂/∂y b → iq b : where we also took into account antisymmetry ofω a µb . The partial derivative ∂ µ W (x, q) is to be taken at fixed q a (not fixed q = e a (x)q a , since that vector has to get boost adjusted to maintain q · u = 0).
To simplify notations below we shall also use the following expression involving derivatives with respect to components of q: (3.21)

IV. FLUCTUATION KINETIC EQUATIONS
The two-point functions W AB (x, q) can be viewed as degrees of freedom additional to the hydrodynamic fields ψ A (i.e., and u) in ways similar to phase-space distribution functions in kinetic theory. This is not just a vague similarity. A certain linear combination of W AB (x, q) can be quantitatively interpreted as phonon distribution function satisfying Boltzmann equation for a particle with momentum q and energy E = c s |q| as will be shown in Section VI. Regardless of this interpretation, these additional degrees of freedom satisfy a coupled differential (matrix) equation which we call somewhat loosely the "fluctuation kinetic equation" or simply "kinetic equation". The kinetic equations have to be supplemented by the usual hydrodynamic field equations of motion (with fluctuation feedback), ∂ µ T µν = 0, to obtain a closed set of equations (somewhat similar to Vlasov equations) to be solved simultaneously. In this section we derive these fluctuation kinetic equations, i.e., equations for W AB .

A. Matrix equation for the Wigner function
We return to Eq. (2.26) for G AB and use it to derive the evolution of the Wigner function defined in the previous section, expressing all derivatives in terms of the confluent derivative. Both definitions of the Wigner functions and of the confluent derivative bring additional terms, but they also lead to many nontrivial cancellations. After a rather lengthy and tedious algebra we find: and its symmetric partner θ µν defined in Eq. (2.6). Note that, within the order of approximation we are working, we can further use the ideal hydrodynamic equation wa µ = −∂ ⊥µ p to eliminate the time-like derivatives of velocity, i.e., a µ , on the right-hand side of Eq. (4.1). This may be useful for numerical solution of the equations which would require solving for time evolution of u(x) simultaneously.

B. Diagonalization and averaging out fast modes
The matrix L (q) in the right hand side of the kinetic equation Eq. (4.1) gives the dominant contribution since it is of order of q whereas the remaining terms are either order k or γq 2 both of which are assumed to be much smaller than q according to our hierarchy of scales in Eq. (1.5). Therefore it is useful to express the kinetic equation in the basis where L (q) is diagonal. L (q) has five eigenvalues: corresponding to 5 eigenvectors ψ A with A = +, −, T 1 , T 2 , . The eigenvectors form a 5 × 5 matrix whereq = q/|q| is the unit vector along q and t (1) and t (2) are two transverse unit vectors that satisfy In other words t (1) , t (2) andq span the spatial hyperplane orthogonal to u, The basis vectors in Eq. (4.5) correspond to the eigenmodes of ideal hydrodynamic equations. Their eigenvalues in Eq. (4.4) correspond to positive and negative frequency sound waves and two degenerate transverse momentum modes. The last zero mode is a consequence of the orthogonality condition Eq. (3.12). The transverse modes are degenerate and the basis in this two-dimensional subspace can be chosen arbitrarily. A convenient explicit choice for t (i) is given in Appendix B.
We can now transform the kinetic equation (4.1) into the diagonal basis of L (q) by the orthogonal transformation M → ψ T M ψ 15 and express the equation in terms of new variables: The modes W A , W B , and W are constrained to vanish by Eq. (3.12). We can therefore view W AB as effectively a 4 × 4 matrix. Furthermore, since in the diagonal basis, ten of the modes, namely W ±∓ , W ±Ti and W Ti± , oscillate with the frequency of order c s q, which is much faster than the background evolution frequency of order c s k. We can use this separation of time scales to introduce (in addition to spatial coarse graining at scale b described in the Introduction) averaging over time intervals of order b t such that After such averaging only six components of the matrix W survive and equations simplify considerably (as noted in [14]). As a result we are left with six modes which can be classified into two sound modes W ±± , two transverse modes W T1T1 and W T2T2 , and two shear modes W T1T2 and W T2T1 . The sound modes are completely decoupled and satisfy 16 where The confluent derivative of W ± is defined as follows: The transverse and shear modes mix and satisfy 2 × 2 matrix equation 17 where and we introduced a covariant q-derivative taking into account rotation of the basis t (i) (x, q) of the transverse modes due to change of q: (4.16) 15 Note that since there are derivatives with respect to x and q in Eq. (4.1), one needs to use ψ T dM ψ = d(ψ T M ψ) + [ψ T dψ, ψ T M ψ]. 16 For notational simplicity we denote W ++ and W −− simply as W + and W − respectively. 17 Here W represents the 2 × 2 matrix W T i T j . Similarly, the ij indices of the 2 × 2 matrices K ij , Ω ij , ω ij µ , ω ij µ and 1 ij = δ ij are suppressed.
The confluent derivative in Eq. (4.14) also includes additional terms associated with ij indices of W (i.e., of W TiTj ), which are due to the x-dependence of the basis vectors t (i) : In Appendix B we propose a simple and intuitive choice for the t (i) basis suitable for applications, and compute corresponding connections ω ij µ and ω ij µ . These fluctuation kinetic equations are the central result of this work. By considering these equations together with the conservation equation for the energy-momentum tensor, including contribution from the fluctuations, such as, e.g., G µν (x)/w in Eq. (2.13), we obtain a closed system of equations that determines the dynamics of both the background flow and the fluctuation correlators self-consistently. In order for this program to work in practice, we need to deal with the singularity of G µν (x) which is manifested in the ultraviolet divergence of the the wave-vector integral relating W AB to G AB . To eliminate the resulting unphysical cutoff dependence we shall absorb ultra-violet divergent contributions of fluctuations into renormalization of a finite number of physical parameters that define first order viscous hydrodynamics, i.e. the equation of state and transport coefficients (viscosities). The remaining part of fluctuation contributions is physical, well-defined and insensitive to the cutoff. In the next section, we describe in detail how this renormalization procedure works.

A. Short-range singularities and renormalization
The "infinite noise" that one has to introduce via the δ-function in Eq. (2.11), which causes cutoff dependence in solutions of the stochastic hydrodynamic equations, does have its counterpart in our deterministic approach. The main advantage of our approach is that it allows us to analytically separate the effects of the cutoff and to absorb them into renormalization of hydrodynamic variables and u, as well as the equation of state p( ) and first-order transport coefficients. This procedure has been discussed by Andreev [13] in non-relativistic context, and also recently in boost invariant Bjorken background in Ref. [15]. In this section we describe how this can be done in relativistic hydrodynamics in arbitrary background flow.
First of all, due to the non-linearity of the energy-momentum tensor in fluctuations, the rest frame defined by the averaged energy-momentum tensor T µν in Eq. (2.13) via Landau's matching is not given by u, i.e., u R = u due to fluctuation contributions. Similarly, the energy density R in the rest frame u R is different from . We shall refer to the hydrodynamic variables (u R , R ) as "renormalized" variables, since they take into accounts the effects of fluctuations. The u R can be found by first observing that u can be shifted to eliminate the terms proportional to G eµ (x) in Eq. (2.13). We also need to keep in mind that due to non-linearities in the constraintȗ ·ȗ = −1, u is not properly normalized; namely Therefore we find the renormalized fluid velocity as 18 For notational simplicity we will denote u R simply as u in the following. Expressing T µν in Eq. (2.13) in terms of u R instead of u using Eq. (5.3), substituting into Eq. (5.1) and multiplying both sides by u Rµ we obtain In terms of these renormalized quantities, we have Note that p = p( ) here is still expressed in terms of "bare" energy density. As usual, due to contribution of short-wavelength fluctuations the coincident point correlators such as G AB (x) ≡ G AB (x, 0) are divergent, i.e., dependent on the wave-vector cutoff Λ. These divergences fall into two classes which are easier to disentangle using the Wigner transform of G AB (x, y), i.e., Fourier transform with respect to y: W AB (x, q) that we define in Section III. In order to study the short-distance cutoff dependence of G AB (x) we need to look at the large-q behavior of W (x, q), since The leading singularity is apparent even in static homogeneous equilibrium, since within our coarsegrained resolution G AB (x, y) ∼ δ 3 (y) (i.e., correlation length is negligible) and thus G The solution of our kinetic equation (4.1) in equilibrium is simply given by where ∆ AB = diag(1, ∆ µν ). Because W (0) is q-independent, the integral in Eq. (5.6) is divergent, i.e., cutoff-dependent: The corresponding contributions to the energy-momentum tensor in Eq. (5.5) can be absorbed into a renormalization of the pressure, i.e., the equation of state. The renormalized pressure is then given by Written in terms of the renormalized energy density given in Eq. (5.4), we obtain the renormalized equation of state as It is worth emphasizing that, even though Λ is an ultraviolet cutoff for the wave-vector of the fluctuating modes, it is still considered small compared to the microscopic scale, i.e. Λ T (see Eq. (1.5)). Therefore the "divergent" contributions of the fluctuations are still small corrections to the averaged background variables that are of order T 4 . However, in numerical simulations, where this separation of scales in not ideal, these corrections will introduce noticeable cutoff dependence. Therefore, we would like to remove these divergent terms not only as a matter of principle, but also as a practical matter. These considerations are not dissimilar in quantum field theories.
In the presence of background gradients, W AB (x, q) deviates from the equilibrium q-independent value. This Wigner function is a solution to an equation we derive in this paper (Eq. (4.1)) -a linear differential equation with coefficients linear in the gradients of velocity. As such, W AB is a non-local functional of those gradients. The fact that allows us to remove divergences by redefining physical parameters (as in quantum field theories) is that the divergent contributions are simply local functions of the velocity gradients.
Since we are interested in the behavior at large q, responsible for divergences, we shall expand in inverse powers of q: where the first and leading term is the equilibrium value (5.7). Using Eq. (5.6) we find correspondingly, The contributions from the first term to the energy-momentum tensor in Eq. (5.5) can be absorbed into the renormalization of the equation of state as shown in Eq. (5.10). The second term in the large q expansion in Eq. (5.11) is a local linear function of velocity gradients. We shall determine W (1) AB in detail in the next section. Here we only need to know that it is of order k/q 2 or, schematically, where γ represents relaxation constants proportional to the viscosities and ∂u represents velocity gradients. Since only scalar, G ee , and tensor, G µν , components appear in the expansion (5.5), we focus on those. The phase space integration in Eq. (5.6) leads to terms linear in Λ which can be decomposed into shear (i.e., traceless) and bulk viscous terms: where the coefficients C shear , C bulk , and C ee are given explicitly in the next section (see Eq. (5.36)). Note that only terms satisfying the orthogonality condition u A (x)G AB (x) = 0 (according to Eq. (2.18)) can appear in Eq. (5.14).
The O(Λ) terms in Eq. (5.5) due to G (1) AB in Eq. (5.14) can be absorbed by the renormalized transport coefficients (namely shear and bulk viscosities) and pressure. The shear (traceless) term in Eq. (5.14) can be absorbed by a renormalization of shear viscosity (it has the same form as the shear part of the viscous term Π µν given in Eq. (2.4)): The remaining terms are related to the renormalization of bulk viscosity. To see how this works, let us look at the trace of the energy momentum tensor given in Eq. (5.5). Separating the Λ 3 and Λ terms explicitly we have It might be tempting to think that the G (0) terms that are independent of flow gradients renormalize the pressure and the G (1) terms that are proportional to θ renormalize the bulk viscosity. However, this is not entirely correct, because we also have to take into account that the relation between and R given by Eq. (5.4) contains G (1) terms, which via the renormalized equation of state p R ( R ) in Eq. (5.10), contribute to what we mean by bulk viscous term separated from the renormalized pressure. Explicitly, we first need to express the bare pressure, p( ), as a function of the renormalized R defined by Eq. (5.4), After inserting this into Eq. (5.16) and using the renormalized equation of state in Eq. (5.10), we obtain The leading term is the correct ideal part in terms of the renormalized equation of state, and the G (1) terms in the second parentheses in Eq. (5.18) can now be absorbed by a renormalization of bulk viscosity 19 It is satisfying to see (and is a non-trivial check) that a conformal symmetry would preserve both the vanishing bulk viscosity ζ = 0 and the conformal equation of state p = /3, according to Eq. (5.10), under fluctuation corrections, by the virtue of c 2 s = 1/3. Finally, the last term in Eq. (5.11), W AB (x, q), is asymptotically of order k 2 /q 4 and does not lead to any large-q divergence in G AB (x) via integration in Eq. (5.6). Unlike G (1) AB (x), G AB (x) is a non-local functional of velocity gradients and is responsible for the physical effects known as long-time tails in hydrodynamic response [13,20]. These terms are finite and constitute the leading corrections to the hydrodynamic derivative expansion. In terms of k these corrections are of noninteger order k 3/2 which are formally in between the viscous first-order O(k) terms, and the second order O(k 2 ) terms in constitutive equations.
After the above renormalization of first order viscous hydrodynamics, we finally obtain the cutoff independent expression for the energy momentum tensor (constitutive equation): where we dropped subscript "R" on the right hand side. It should be understood that all quantities in Eq. (5.20) and in the kinetic equations such as Eq. (4.1) are renormalized. For example, the pressure, p( ), in Eq. (5.20) is given by the physical equation of state, which could be, e.g., obtained from a lattice calculation. The conservation equations for this tensor together with the fluctuation kinetic equations (4.11) and (4.14) form a closed set of cutoff-independent hydrodynamic equations. 20 The necessary components of G AB (x) can be obtained by solving the kinetic equations (4.11) and (4.14) for W and subtracting the leading large-q contributions W (0) and W (1) determined as local functions of hydrodynamic variables and gradients of velocity by expressions we shall derive in the next section (Eqs. (5.23), (5.29), (5.30) and (5.32)). Alternatively, with the explicit expressions for W (0) and W (1) given below, one can substitute (5.11) into Eqs. (4.11) and (4.14) and solve the resulting equations for W directly. Using Eq. (5.6) we can then determine

B. Large-q behavior of Wigner functions
As we discussed in the previous section, in order to obtain finite, cutoff independent equations we need to separate the leading and subleading large-q terms from the Wigner functions, since these terms should be absorbed by renormalization of equation of state and kinetic coefficients.
The leading term is easy to find, since for large q the x-dependence of the background can be neglected and only the relaxation term γq 2 (W − T w) in Eqs. (4.11) and (4.14) should be kept, with the solution given simply by local equilibrium values of fluctuations: Ti,Tj = T wδ ij .
(5.23) 19 In principle, the renormalization of equation of state p( ) leads to a corresponding renormalization of temperature T ( ). However, since T itself appears only in the noise amplitude and therefore only in fluctuation-induced correction, this is a higher order effect. 20 The usual concerns about the acausal response and associated instabilities in this equation can be addressed by the standard Israel-Stewart treatment introducing relaxation dynamics for the viscous tensor. This modification affects the regime beyond the domain of applicability (k T ) of hydrodynamics [24], and we shall leave it outside the scope of this paper, as part of the set of established procedures (such as, e.g., discretization) needed for numerical implementation. This is the consequence of the thermal noise in Eq. (2.11) satisfying fluctuation-dissipation theorem.
In the presence of gradients, the solution deviates from local equilibrium at x: We can substitute this ansatz into Eq. (4.11) and Eq. (4.14) and use the ideal hydrodynamic equations and thermodynamic relations s = w/T , dp = sdT to evaluate derivatives of W (0) = T w to leading order in flow gradients: As a result we obtain equations for W (neq) AB (x, q). However, as discussed in the previous section, W AB integrated over q, still produces an ultraviolet divergence (albeit linear in Λ and not Λ 3 ). To isolate this divergence we write where we define W AB (x, q) as the leading term in W (neq) AB (x, q) in the large q limit. To find this term we note that the terms in the equation for W (neq) AB obtained from Eqs. (4.11) and (4.14) by substituting Eq. (5.24) fall into two classes: the terms proportional to W (neq) (or its derivatives) and the terms independent of W (neq) . Within each class we identify the leading terms in the limit of q → ∞ and require that these terms cancel when we replace W (neq) with its leading term, W (1) . This gives us the following equations: which are easily solved as: As expected (see Eq. (5.13)) these terms are of order ∂u/(γq 2 ) and lead to order Λ terms after the q integration in Eq. (5.6). The remaining terms in W (neq) in Eq. (5.26), i.e., For finite q (not satisfying γq 2 k) the dependence of W on q can be represented to linear order in ∂u, schematically, as where v = ±c sq or 0 depending on which mode we are considering. Integration over d 3 q leads to, also schematically, G(x) ∼ k 1/2 ∂u/γ 3/2 . The non-integer power of k represents the fact that G is a non-local functional of the gradients of ∂u. These non-local terms give rise to power-law (in space and/or time) tails in hydrodynamic response [13]. 21

C. Renormalization of transport coefficients
Let us now calculate the renormalization of shear and bulk viscosities using Eq. (5.29). First of all, we convert back into the original e, µ basis in order to plug them into the energy momentum tensor, Eq. (5.5). The conversion is given by: In particular, (5.29) is converted into With the help of the integrals Finally inserting these expressions into Eq. (5.15) and (5.19) we obtain the renormalized shear and bulk viscosities These expressions agree with the results which were computed via different methods in the earlier literature (e.g., Eq. (51), Eq. (A14) and footnote 7 in Ref. [15]). The positivity of correction to ζ R is remarkably non-trivial. It follows from appropriately renormalizing the energy density as well as the equation of state as described by Eqs. (5.10) and (5.17). It is satisfying to see that the corrections to viscosities are positivedefinite in agreement with the Second Law of Thermodynamics.

VI. PHONON INTERPRETATION OF THE FLUCTUATION KINETIC EQUATION
A. Phonon kinetic equation Consider a classical particle whose motion is described in terms of the space-time vector x µ and 4momentum p µ with dispersion relation given by some condition F (p) = 0. For example, for a massive particle in vacuum F = p 2 − m 2 . A phonon dispersion relation is given by p 0 = E(p) ≡ c s |p| in the rest frame of the fluid. This can be represented by F + (p) = p · u + E(p ⊥ ), (6.1) where u is the the 4-velocity of the fluid rest frame E = c s |p ⊥ | and The classical action can be then written as where λ is a Lagrange multiplier. Variation of the action is given by: Classical trajectory is then given by equations of motioṅ where dot denotes d/(λdτ ) (or one can use reparametrization invariance to set λτ to equal coordinate time x 0 in frame u) and (where we used ∂p ⊥ν /∂p µ = ∆ µ ν ) as well aṡ together with the condition F + = 0. We consider local properties of the fluid to be varying (sufficiently slowly) in space and time. I.e., u µ = u µ (x), as well as E = E(x, p ⊥ ), which in the case of a phonon means c s = c s (x). The corresponding Liouville operator acting on a function N (x, p) is given by Note that L[F + ] = 0. This property is important because it allows us to restrict the 8-dimensional phase space to the 7-dimensional subspace defined by F + = 0, i.e., to consider functions of the form 9) where N is the usual phase-space distribution function (of 7 variables only). In other words In order to write the kinetic equation in terms of the distribution function N (x, p ⊥ ) we need to express x derivatives in L[N ] at fixed p (∂/∂x µ in Eq. (6.8)) in terms of x derivatives at fixed p ⊥ . These derivatives are not the same because the relationship between p and p ⊥ depends on x (via u(x) in Eq. (6.2)). One finds where we denoted by∇ µ the x derivative at p ⊥ fixed 22 . Corresponingly, the last term in Eq. (6.7) should be written as Similarly, the p derivatives at fixed x should be expressed as p ⊥ derivatives Substituting Eqs. (6.5), (6.7), (6.11), (6.10) and (6.12) into Eq. (6.8) we find Finally, using ∂ µ p ⊥ν = −E∂ µ u ν + u ν ∂ µ (p · u), we can write the Liouville operator as 14) The expression in the square brakets is (the negative of) the force acting on the phonon. 23 The two terms in parenthses multiplied by E are easily recognized as the inertial force due to acceleration a and the Coriolis force due to rotation ω µν , respectively. The force −p ⊥ν ∂ ⊥µ u ν is easier to understand by considering isotropic Hubble-like expansion, i.e., such that ∂ ⊥µ u ν = H∆ ν µ , where H is the rate of expansion (Hubble constant). This term then describes the rescaling of the momentum p ⊥ (stretching of the sound wave) due to expansion of the background medium, leading to the "red shift" of the phonon spectrum, similar to the photon red shift in the expanding universe. The last term is the force due to the dependence of energy on the location of the phonon via the coefficient c s in its dispersion relation: Remarkably, upon changing the notation for the phonon momentum the Liouville operator in Eq. (6.14) with E = c s |p ⊥ | is identical to the one in Eq. (4.11) obtained using completely different (but apparently complementary) considerations. The two signs in front of c s in Eq. (4.11) correspond to positive and negative frequency sound waves, or positive/negative energy solutions of the condition where F ± = (p · u) ± E and the positive energy solution is given by F + = 0 in Eq. (6.1). Curiously, for linear dispersion, E = c s |p ⊥ |, the condition in Eq. (6.17) can be written as g µν p µ p ν = 0 using flow induced effective "metric tensor" g µν = −u µ u ν + c 2 s ∆ µν . Since d(F + F − ) = F − dF + + F + dF − and δ(F ) = δ(F + )/F − + δ(F − )/F + , we see that the equations of motion localized on the F + = 0 surface are given by Eqs. (6.5) and (6.7) up to rescaling of proper time. On the other hand, the equations of motion with the constraint F + F − = 0 are given bẏ from which one can derive the "geodesic" equation of motion by taking additional time derivative to the first equation and using these equations once more. From this point of view the forces in Eq. (6.14) can be viewed as "gravitational" forces. Perhaps even more remarkably than the matching of the Liouville operators in Eqs. (6.14) and (4.11), the identification W ± (x, q) = c s |q|wN ± (x, q) (6.19) leads to nontrivial cancellation of the whole second line in Eq. (4.11) (i.e., of all terms proportional to the background gradients θ µν and a µ times W ± ) leaving simply the relaxation term in Eq. (4.11): (6.20) where L ± are different by the sign in front of c s in Eq. (4.11). Note that the equilibrium value of N ± , N ± /(c s |p ⊥ |w), equals T /E as expected for the low-energy limit of the phonon Bose-Einstein distribution function.
In contrast to Eqs. (4.11) for longitudinal modes which reduces to a simple form Eq. (6.20) upon rescaling given by Eq. (6.19), Eq. (4.14) for transverse modes cannot be simplified in this way. This may be related to the fact that there is no quasiparticle interpretation for these non-propagating, diffusive modes.

B. Phonon contributions to stress-energy tensor
It is also remarkable that certain contributions of the fluctuations to stress-energy tensor can be related directly to the stress-energy tensor of the phonon gas. This provides a justification to the two-fluid picture (hydrodynamic fluid plus gas of phonons) which guided the original approach by Andreev [19].
Let us start with the expression for the stress tensor for one particle moving along a trajectory specified by x(τ ): This means for a gas of such particles and holes with distribution functions N ± (x, p) we have ('s' for 'sound'): (Minus here is because a hole is the absence of a particle). Using equation of motionẋ = u ± v (Eq. (6.5)) with ± for particles/holes we obtain: Using now E = c s |p ⊥ | for the phonon, we find: Comparing to Eq. (5.32) (neglecting transverse modes) and identifying N ± = W ± /(c s |p ⊥ |w) (as in Eq. (6.19)) we can write: We can recognize the last two terms as the last two terms in Eq. (2.13). The first term contributes to the shift/renormalization of the energy density (see Eq. (5.4)). Phonons contribute to pressure, i.e., renormalize equation of state. The equilibrium pressure of the phonon gas equals Together with the contribution from the shift of the energy due to phonons, −c 2 s G (0) µ µ (x)/w, this reproduces the last term before the second equal sign in Eq. (5.10).

VII. CONCLUSIONS AND OUTLOOK
We have derived a set of equations which describe coupled evolution of hydrodynamic variables and u and the correlation functions of their fluctuations, collectively denoted by φ A . The correlation functions are expressed as Wigner transform W AB (x, q) of the equal-time correlator φ A φ B .
An essential feature of this approach, which distinguishes it from the stochastic approach also pursued in the literature, is the possibility to cleanly separate the short-range singularities due to "infinite noise" and absorb them into renormalization of the variables, equation of state and transport coefficients. The success of this procedure relies on the locality of the short-range singularities in the same way as in quantum field theoretical renormalization due to ultraviolet singularities. The resulting constituitive equations for stress-energy tensor Eq. (5.20) contain only finite contributions from fluctuations.
One of the crucial new issues we tackled is the treatment of the "equal-time" in the definition of the correlators. Frame-dependence of simultaneity is a quintessential relativistic effect and has not been an issue in the earlier work [13] where similar equations have been considered in non-relativistic context. Our analysis led us to a definition of correlators, Wigner functions and x-derivatives adjusted for the changing local rest frame of the fluid. We refer to the objects which account for the flow in this way as confluent.
The success of our approach relies significantly on the separation of scales inherent in the hydrodynamic regime. The scales of homogeneity, L, must be longer (equivalently, the corresponding wave vector k = 1/L, must be softer) than the range of the correlations. In hydrodynamics, this range, eq , is of order γL/c s L. The corresponding "hard" momentum q ∼ 1/ eq is much larger than the "soft" momentum k.
Given that this separation of scales is similar to the separation of scales which leads to kinetic regime in weakly-coupled quantum field theories, it is not a coincidence that the equations for the Wigner functions we obtain are similar to kinetic equations. It is, nevertheless, remarkable that the equations we obtain by focusing on the longitudinal mode fluctuations completely coincide with kinetic equations describing phonon gas on a background with arbitrary non-uniform flow. This includes nontrivial inertial, Coriolis, and "Hubble" forces -Eq. (4.11) vs Eq. (6.14).
The originial approach by Andreev postulated Hamiltonian kinetic equations for a phonon distribution function without derivation [13]. We derive these equations directly from stochastic hydrodynamics and confirm that the collision/relaxation term has the simplest form assumed in Ref. [13] and does not depend on the gradients of the flow. This result emerges after nontrivial cancellations (Eq. (6.20) vs Eq. (4.11)), which appear even more nontrivial given that, in contrast, the "kinetic" equation for the transverse modes does contain flow gradients among the relaxation terms (Eq. (4.14)). These gradient terms were assumed to be absent in Ref. [13]. We plan to explore possible implications of these gradient-dependent relaxation terms in the future.
Although it would be interesting to study possible analytical solutions of our equations and their consequences, we also hope that these equations will find application in numerical simulations of the evolution of heavy-ion collisions or other relativistic many-body systems where fluctuations are important. Our equations can be directly simulated for any flow, not necessarily limited by Bjorken boost-invariance assumption (as in, e.g., Ref. [14]). In particular, the effects of vorticity, absent in the Bjorken flow, but important in heavy-ion collisions [17], can be studied.
The approach based on the evolution of correlation functions has been also introduced recently to describe the dynamics near a critical point [25]. The extension of hydrodynamics, or Hydro+, by a slow mode describing evolution of critical fluctuations towards equilibrium is a particular example of the correlation function (hydro-kinetic) approach to hydrodynamic fluctuations. It would be interesting to re-derive Hydro+ formalism in the framework used in this paper. In order to do this we must generalize the formalism to include conserved current (baryon current in QCD), which we defer to future work 24 . Such a generalization in essential for the hydrodynamic modeling of fluctuations and its effects in the the beam-energy scan program [8] conducted at RHIC.
Finally, recent advances in formulating hydrodynamics as an effective field theory on Schwinger-Keldysh path-integration contour, in principle, allows using powerful field-theoretical methods and insights, including a diagrammatic machinery for calculating real-time correlation functions (see, e.g., Refs. [26,27] for reviews). However, the practical usage is so far mostly limited to correlation functions in simple backgrounds, such 24 Special cases of static and boost invariant backgrounds for a conformal fluid with conserved charge have been discussed in Ref. [16].
as, e.g., static equilibrium, and it is not yet clear how to apply this approach to a realistic heavy-ion collision simulation. It would be interesting to establish an explicit connection between the formalism we present here and the Schwinger-Keldysh effective field theory. This could provide better understanding of some conceptual issues, such as renormalisation at higher (or even all) orders in hydrodynamic expansion, and help generalize the approach to tackle higher-order correlation functions (beyond two-point functions discussed in this paper.).
For certain flow configurations u(x) it may be possible to find a choice of triad fields e a (x) which makes the spin connectionω vanish. This requires integrability of Eq. (3.15) withω = 0, which means that the change of vector e a obtained by integrating Eq. (3.15) withω = 0 between two points should not depend on the path, i.e., dx λωµ λν e ν a = 0.
We can then define equal-time correlator in such a way that points x and x ± = x ± y/2 lie on the same curved surface τ = const. Using 3-dimensional vector y in the tangent plane to u(x) to parameterize points on such a surface, we can write explicitly where y · u(x) = 0 and the last term describes the curvature of the surface. Using this definition of x ± instead of Eq. (1.7) will change the definition of the "equal-time" correlator and of the Wigner function.
In what follows in this section we shall use that modified definition, but retain the same notation for simplicity. Due to the modification described above, we must replace L (y) defined in Eq. (2.27) by where Θ λ ≡ c s 2 0 θ νλ θ µλ 0 .
As a consequence, Eq. (2.26) and (4.1) are modified and take the form Note, that the only change compared to Eq. (4.1) is in the double commutator term. 26 For the Bjorken flow, a µ = ω µν = 0, and all the terms on the second line in Eq. (C9) vanish.
To complete the comparison, for the boost-invariant flow, we perform the coordinate transformation from (t, x, y, z) to (τ, x, y, Y ) given by t = τ cosh Y, x = x, y = y, z = τ sinh Y , where τ is the proper time and Y is the space-time rapidity. One can easily check that for the Bjorken flow u ·∇ = ∂ τ , θ = 1/τ , a µ = ω µν = 0. Thus Eq. (C7) is reduced to ∂ τ W (x; q) = − iL (q) + K (a) , W − 1 2L + Q (q) + K (s) , W + 2T wD (q) + 1 τ W + q z τ ∂W ∂q z . (C9) FIG. 3. Left: Illustration of the surface orhogonal to the conservative flow u at each point. Boost is represented by ordinary rotation, preserving angles, for clarity. Right: The same is not possible for non-conservative flow, i.e., for nonzero vorticity. However, it is possible to make the normal vector to the surface (not shown) and the flow vector u (shown) at the same point be different by a purely vortical vector: v µ = ∂ µ τ − u µ , such that ∂ · v = 0.
Since q Y = τ q z where q Y is the wave vector conjugate to Y , we define W B (x; q Y ) = W (x; q z )/τ to take into account the change in the measure of the momentum integration. Using where the last two terms in Eq. (C9) were eliminated by the momentum rescaling. Similarly, one can check that our Eqs. (4.11) and (4.14), rewritten in terms of W B , will reduce to Eq. (A7) in Ref. [15] exactly. k -wave vector Fourier conjugate to midpoint vector x; N (x, p) -phonon phase-space distribution function -Eq. (6.9); N ± (x, p) -phase-space distribution function for positive/negative frequency phonons; p -phonon momentum (not to be confused with pressure p( )) -Section VI; p ≡ p( ) -pressure at average energy density ; q -wave vector Fourier conjugate to separation vector y, also q · u = 0 -Section III; T -temperature (local value at (x)); t (i) µ -(i = 1, 2) momentum space diad vector orthogonal to q µ -Eq. W ij -2 × 2 matrix of Wigner functions for transverse modes -Eq. (4.14); x -in a two-point correlator -the midpoint space-time vector; x ± -arguments of the 2-point correlator; y -in an equal-time two-point correlator -the separation vector, constrained by u(x) · y = 0.