Relativistic Hydrodynamic Fluctuations

TL;DR

This paper develops a comprehensive formalism for analyzing fluctuations in relativistic hydrodynamic flows, incorporating feedback effects and providing a framework that resembles kinetic equations, with applications to numerical simulations.

Contribution

It introduces a systematic approach to describe relativistic hydrodynamic fluctuations using equal-time correlation functions and confluent Wigner transforms, including renormalization techniques.

Findings

Derived equations for confluent Wigner functions resembling kinetic equations.

Established the equivalence between phonon propagation and Wigner function dynamics.

Presented a renormalization procedure for short-distance singularities.

Abstract

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.

Equations

$$b\gg {\ell }_{\mathrm{mic}}.$$

$$L\gg b.$$

$${\ell }_{\mathrm{eq}}\sim \sqrt{\gamma {\tau }_{\mathrm{ev}}}\sim \sqrt{\gamma L/c_s},$$

$${\ell }_{\mathrm{eq}}\sim \sqrt{L/T}\ll L,$$

$$k,\gamma q^2/c_s\ll q\ll \Lambda \ll T,$$

$$G_{AB}(x,y)\equiv ⟨{\varphi }_A(x^{+}){\varphi }_B(x^{-})⟩$$

$$x^{\pm }=x\pm y/2$$

$${\partial }_{\mu }{\breve{T}}^{\mu \nu }=0.$$

$$-{\breve{T}}_{\nu }^{\mu }{\breve{u}}^{\nu }=\breve{ϵ}{\breve{u}}^{\mu }.$$

$$T^{\mu \nu }(ϵ,u)$$

$${\Pi }^{\mu \nu }$$

$${\Delta }^{\mu \nu }=g^{\mu \nu }+u^{\mu }{u}^{\nu },$$

$${\theta }^{\mu \nu }=\frac{1}{2}\left({\partial }_{\perp }^{\mu }{u}^{\nu }+{\partial }_{\perp }^{\nu }{u}^{\mu }\right),\theta ={\theta }_{\mu }^{\mu },$$

$${\partial }_{\perp \mu }={\Delta }_{\mu }^{\nu }{\partial }_{\nu }.$$

$${\breve{T}}^{\mu \nu }$$

$$\breve{u}=u+\delta u,\breve{ϵ}=ϵ+\delta ϵ.$$

$$⟨\delta u⟩=⟨\delta ϵ⟩=0.$$

$$⟨{\breve{S}}^{\mu \nu }(x){\breve{S}}^{\lambda \kappa }(x^{\prime })⟩$$

$$p(\breve{ϵ})=p(ϵ)+c_s^2\delta ϵ+\frac{1}{2}\frac{d{c}_s^2}{dϵ}\delta {ϵ}^2+\mathcal{O}(\delta {ϵ}^3),$$

$$⟨{\breve{T}}^{\mu \nu }(x)⟩$$

$${\varphi }_A\equiv (\delta e,\delta g_{\mu })\equiv (c_s\delta ϵ,w\delta u_{\mu })$$

$$G_{AB}(x)\equiv ⟨{\varphi }_A(x){\varphi }_B(x)⟩,$$

$$G_{AB}(x)=G_{AB}(x,0).$$

$$u_A\equiv (0,u_{\mu }).$$

$$u^A(x^{+})G_{AB}(x,y)=G_{AB}(x,y)u^B(x^{-})=0.$$

$${\dot{c}}_s=\frac{w}{c_s}\frac{d{c}_s}{dϵ}=\frac{d\log c_s}{d\log s}=\frac{1}{2}T\frac{d{c}_s^2}{dT},$$

$${\breve{T}}^{\mu \nu }$$

$${\gamma }_{\eta }=\eta /w\text{mbox}and{\gamma }_{\zeta }=\zeta /w.$$

$${\partial }_{\mu }{\breve{T}}^{\mu \nu }$$

$$u\cdot\partial\phi_{A}=-\big{(}{\mathbb{L}}+{\mathbb{Q}}+{\mathbb{K}}\big{)}_{AB}\phi^{B}-{\xi}_{A}\,,$$

Full text

I Introduction

I.1 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 Landau:2013fluid 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 Jeon:2015dfa ; Romatschke:2017ejr . 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 $\mathcal{O}(10^{24})$ particle degrees of freedom. Heavy-ion collisions, however, occupy a “sweet spot” in terms of the size: with $\mathcal{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 Stephanov:1998dy ; Stephanov:1999zu ; Stephanov:2008qz ; Stephanov:2011pb in the beam energy scan experiments Aggarwal:2010cw .

In the classic Landau-Lifshitz Landau:2013stat2 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 Kapusta:2011gt ; Kapusta:2012zb ; Young:2014pka . 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 Andreev:1978 . 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 111Despite 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. ) has been recently discussed in a relativistic context Akamatsu:2017 ; Akamatsu:2018 ; Martinez:2018 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 STAR:2017ckg . The aim of this paper is to develop such a universal approach.

I.2 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 $\bm{x}$. In order to be macroscopic, the length ${b}$ must greatly exceed microscopic scales, $\ell_{\rm mic}$, such as mean-free path in a weakly-coupled system, or thermal length $1/T$ in a relativistic strongly coupled system

[TABLE]

The resulting coarse-grained variables can be used to describe evolution of inhomogeneities at larger length scales

[TABLE]

To facilitate the discussion let us refer to the hydrodynamic variables defined via coarse-graining discussed above as ${\breve{\psi}}_{A}(t,\bm{x})$, where index $A$ labels a variable. Since we are describing a thermal system, the variables ${\breve{\psi}}_{A}$ are stochastic – fluctuating between members of the statistical ensemble describing our system (in heavy-ion collisions – between collision events). To be more precise, variables ${\breve{\psi}}_{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.222The 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 $c_{s}/{b}$ must therefore be much smaller than $T$, i.e., ${b}\gg c_{s}/T$. Due to coarse graining, fluctuations at scales shorter than $b$ are averaged out, i.e., suppressed. In this sense, $\Lambda=1/{b}$ plays the role of the ultraviolet (wave-vector) cutoff.

In order to describe fluctuations we introduce the ensemble averages of the variables $\psi_{A}\equiv\langle{\breve{\psi}}_{A}\rangle$. The ensemble averages $\psi_{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 equal-time products at two space-time points: $\langle\phi_{A}(t,\bm{x}_{1})\phi_{B}(t,\bm{x}_{2})\rangle$, where $\phi={\breve{\psi}}-\langle\psi\rangle$ is the fluctuating part of the variable, as usual.333More generally, the variable $\phi_{A}$ could be (and will be) a linear combination of ${\breve{\psi}}_{B}-\langle\psi_{B}\rangle$.

In equilibrium, the correlators $\langle\phi_{A}(t,\bm{x}_{1})\phi_{B}(t,\bm{x}_{2})\rangle$ are translationally invariant, i.e., independent of the midpoint $\bm{x}\equiv(\bm{x}_{1}+\bm{x}_{2})/2$ at fixed separation $\bm{y}\equiv\bm{x}_{1}-\bm{x}_{2}$. The dependence of correlation functions of operators in equilibrium on separation $\bm{y}$ is characterized by exponential fall-off at distances larger than correlation length $\xi$: $e^{-|\bm{x}_{1}-\bm{x}_{2}|/\xi}$. Correlation length $\xi$ is a microscopic scale, typically: $\xi\sim\ell_{\rm mic}$. 444Near a critical point correlation length is large, e.g., $\xi\gg 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 $\xi$: $b\gg\xi$. Dynamics of fluctuations near a critical point in the opposite regime, ${b}\ll\xi$, is characterized by dynamical scaling and is discussed in Hohenberg:1977 . From the point of view of the coarse-grained variables $\phi_{A}$, therefore, the distance $\xi\ll{b}$ is negligible and the equilibrium correlator $\langle\phi_{A}(t,\bm{x}_{1})\phi_{B}(t,\bm{x}_{2})\rangle$ is essentially a multiple of the delta-function $\delta^{3}(\bm{x}_{1}-\bm{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 $\tau_{\rm ev}\sim L/c_{s}$, 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 $\tau_{\rm ev}$:

[TABLE]

where $\gamma$ is an appropriate diffusion constant (typically, $\gamma\sim 1/T$). 555In terms of the notation $k^{*}$ from Ref. Akamatsu:2017 : $\ell_{\rm eq}=1/k^{*}$. This means for distances ${b}\ll|\bm{y}|\ll\ell_{\rm eq}$ or, more precisely, for wave-vectors $\bm{{q}}$ conjugate to $\bm{y}$ such that $1/\ell_{\rm eq}\ll|\bm{{q}}|\ll\Lambda$ the equilibration is complete. However, at scales around ${q}\sim 1/\ell_{\rm eq}\sim\sqrt{c_{s}k/\gamma}$, where $k\sim 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\sim\ell_{\rm eq}$ are parametrically shorter than $L$,

[TABLE]

we can consider the equilibration process as local, occurring on a slowly varying background set by local values of $\psi_{A}(t,\bm{x})$. For that reason it is also naturally convenient to use the mixed Fourier (i.e., Wigner) transform of the correlation function $\langle\phi_{A}(t,\bm{x}+\bm{y}/2)\phi_{B}(t,\bm{x}-\bm{y}/2)\rangle\equiv G_{AB}(x,\bm{y})$ with respect to separation vector $\bm{y}$, which we shall denote ${W}_{AB}(x,\bm{{q}})$. The relevant values of ${q}$ will satisfy

[TABLE]

where, for simplicity, we took $\ell_{\rm mic}^{-1}\sim T$. The condition $k\ll{q}$ allows us to treat background as smooth when describing the relaxation of correlations ${W}_{AB}(x,\bm{{q}})$ to equilibrium. However, we must retain non-zero gradients, $\partial_{\mu}\psi_{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 $\Lambda$ contribute to this feedback. The integral of the contributions over $q$ is divergent, i.e., it depends polynomially on the cutoff $\Lambda$. 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 $\Lambda$, 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}\sim 1/\ell_{\rm eq}$. Most importantly, it is finite and cutoff independent. The magnitude of these non-equilibrium effects can be estimated as the phase-space volume $\int d^{3}{q}\sim\ell_{\rm eq}^{-3}\sim(c_{s}k/\gamma)^{3/2}$. The power $3/2$ indicates that these effects are non-local. They are known as “long time tails” of hydrodynamic response andreev1970twoliquid ; Andreev:1978 ; Kovtun:2003 ; Kovtun:2011np . Their contribution is typically more important than that of the second order $\mathcal{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) 666The second order terms are also used in numerical applications of hydrodynamics to ensure causality and stability muller ; is ..

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

[TABLE]

and evaluate it at points

[TABLE]

where 4-vector $y$ lies in the hyperplane orthogonal to $u(x)$: $u(x)\cdot y=0$. 777In Section C we shall discuss the choice of the equial-time hypersurface which is not a plane, so as to take into account the variation of $u$ between points $x$ and $x_{\pm}$ and see what, if any, modifications of the results this entails. One can anticipate that these modifications will be insignificant because the typical range of the correlation function, $\sqrt{\gamma L/c_{s}}$, is short compared to the scale $L$ over which the background $u$ varies significantly. 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.2 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}\gg 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.2.

Several Appendices contain useful supplementary information. In particular, we assemble a list of our notation choices used throughout the paper in Appendix D.

II Stochastic hydrodynamics and fluctuations

II.1 Stochastic hydrodynamics

Hydrodynamic equations express conservation and transport of energy and momentum densities:888For simplicity we do not consider any additional conserved charge in this paper. This generalization will be addressed in future work.

[TABLE]

To simplify notations later in the paper we label fluctuating hydrodynamic quantities with an accent, as in $\breve{T}^{\mu\nu}$, to distingish them, where necessary, from quantities which are not fluctuating. The four conservation equations (8) are solved for the same number of hydrodynamic variables. A convenient covariant choice for them is the fluid velocity $\breve{u}^{\mu}$ (normalized as $\breve{u}\cdot\breve{u}=-1$) and the energy density $\breve{\epsilon}$ in the rest frame of the fluid, that are defined by the Landau’s matching condition

[TABLE]

To form a closed system, we need six additional (constitutive) equations to express all components of $T^{\mu\nu}$ in terms of $\epsilon$ and $u^{\mu}$. For macroscopically large scale dynamics of hydrodynamic variables, $T^{\mu\nu}$ can be expanded in gradients of $\epsilon$ and $u^{\mu}$. The first-order (Landau-Lifshitz) hydrodynamics corresponds to truncating this expansion at first order in gradients:

[TABLE]

where $p(\epsilon)$ is pressure as a function of $\epsilon$ – also known as the equation of state and $w(\epsilon)=\epsilon+p(\epsilon)\,$ is the enthalpy. The viscous tensor is linear in gradients of $u$:

[TABLE]

where shear and bulk viscosities are denoted as $\eta$ and $\zeta$, respectively, and

[TABLE]

is the projection operator to the spatial hypersurface orthogonal to $u$, in terms of which we define:

[TABLE]

where

[TABLE]

However, the constitutive equations (10) relating $T^{\mu\nu}$ and the hydrodynamic variables are valid only on average, and there exist random local thermal noise $\breve{S}^{\mu\nu}$ which makes Eq. (8) a stochastic differential equation with

[TABLE]

The functions of hydrodynamic variables such as $w$, $p$, $\Pi^{\mu\nu}$, etc. in Eq. (15) are the same as in Eq. (10) and Eq. (11) but they are evaluated for fluctuating variables $\breve{u}$ and $\breve{\epsilon}$.

The hydrodynamic variables in $\eqref{eq:Pimunu}$ fluctuate as they are driven by the random noise $\breve{S}^{\mu\nu}$, and we need to consider statistical ensemble average over these fluctuations for any observables on macroscopic scales. We write our stochastic hydrodynamic variables $\breve{u}^{\mu}$ and $\breve{\epsilon}$ as a sum of their averages, $u\equiv\langle{\breve{u}}\rangle$, $\epsilon\equiv\langle{\breve{\epsilon}}\rangle$, and linear fluctuations around them as:

[TABLE]

By definition, the linear fluctuations vanish upon averaging

[TABLE]

These fluctuations are driven by the noise term $\breve{S}^{\mu\nu}$ with $\langle{\breve{S}^{\mu\nu}(x)}\rangle=0$, the strength of which is set by the fluctuation-dissipation theorem999Due to the presence of gradients, the system is slightly out of equilibrium and the fluctuation-dissipation relation given in Eq. (18) 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. (18) with $T$ and $w$ being functions of $x$ in the remainder of the paper.,

[TABLE]

In principle, it is possible to numerically solve the stochastic equation $\partial_{\mu}\breve{T}^{\mu\nu}=0$ with some coarse-graining, or wave vector cutoff $\Lambda$, which regularizes the infinite amplitude of the noise arising from the $\delta^{(4)}(x-x^{\prime})$ term. However, as we already mentioned in the introduction, the results would depend sensitively on the cutoff $\Lambda$ due to non-linearity of hydrodynamic equations.

We follow an alternative approach, that is, we include fluctuation contributions to $\langle{\breve{T}^{\mu\nu}}\rangle$ by expanding $\breve{T}^{\mu\nu}$ to second order in fluctuations. The fluctuation contributions to $\langle{\breve{T}^{\mu\nu}}\rangle$ 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 $\partial_{\mu}\langle{\breve{T}^{\mu\nu}}\rangle=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 $(\epsilon,u)$ and $T^{\mu\nu}$ in Eq. (10), including non-linearities in the equation of state, such as

[TABLE]

where $c_{s}^{2}=dp(\epsilon)/d\epsilon$ is the square of sound speed, $\breve{T}^{\mu\nu}$ in Eq. (15) contains terms that are nonlinear in fluctuations. Expanding $\breve{T}^{\mu\nu}$ up to second order in fluctuations and taking the average, we have:

[TABLE]

Note that we neglected the fluctuations of the viscous part $\Pi^{\mu\nu}$, which are parametrically smaller than the terms kept in the above expansion101010 We rely on $\gamma{q}\sim{q}/T\ll 1$, according to Eq. (5), where $q$ is the typical wave vector of the fluctuations.. In the last line we introduced the collective notation for the fluctuating modes, $\delta\epsilon$ and $\delta u^{\mu}$:

[TABLE]

where the scalar $\langle{\delta e\delta e}\rangle$, vector $\langle{\delta e\delta g^{\mu}}\rangle$, and tensor $\langle{\delta g^{\mu}\delta g^{\nu}}\rangle$ components of the two-point correlation function are expressed compactly as

[TABLE]

where $A\in(e,0,1,2,3)$111111The mixed index $A\in(e,0,1,2,3)$ is raised and lowered by the ”metric”, $\text{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.. In terms of our definition of the correlator $G_{AB}(x,y)$ in Eq. (6),

[TABLE]

We can express the fluid velocity in the collective notation as well:

[TABLE]

It should be noted that not all five variables $\phi_{A}$ are independent since, due to normalization $\breve{u}\cdot\breve{u}=-1$, we have a constraint $u^{A}\phi_{A}=0$. Correspondingly,

[TABLE]

Finally, we define

[TABLE]

where $s=w/T$ is the average entropy density.

The functions ${G}_{AB}(x)$ in Eq. (20) are, in general, non-local functionals of the background fields $\epsilon$ and $u$. In the next section, we will derive the evolution equation for them by using the linearized hydrodynamics equation of motion for fluctuations.

II.2 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, $\delta e$ and $\delta g^{\mu}$. 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

[TABLE]

where

[TABLE]

In this expansion the first two terms are zero’th order in gradients and the third term $\Pi^{\mu\nu}$ is of first order or, equivalently, of order $k$. These three terms constitute the average background value without fluctuation contributions, i.e., $T^{\mu\nu}$ in Eq. (10). 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\ll 1$ or $\gamma k\sim k/T\ll 1$ compared to the terms being kept, according to our hierarchy of scales in Eq. (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 $\partial_{\mu}\breve{T}^{\mu\nu}=0$,

[TABLE]

where we also neglect several terms based on similar considerations discussed above. By averaging both sides we obtain $\partial_{\mu}T^{\mu\nu}=0$ to leading order in fluctuation expansion, which we insert back into Eq. (29) to arrive at a stochastic differential equation for the linearized fluctuations.

In terms of the notation $\phi_{A}$ introduced in Eq. (21), the equation for the linearized fluctuations reads:

[TABLE]

where ${\mathbb{L}}$, ${\mathbb{Q}}$, and ${\mathbb{K}}$ are $5\times 5$ matrix operators. The operators ${\mathbb{L}}$ and ${\mathbb{Q}}$ are the ideal and dissipative terms, respectively, ${\mathbb{K}}$ contains the corrections due to the first-order gradients of background flow, and ${\xi}$ denotes the random noise. Explicitly121212It is useful to keep in mind the power counting according to our hierarchy of scales in Eq. (5): ${\mathbb{L}}\sim{q}$, ${\mathbb{Q}}\sim\gamma{q}^{2}$ and ${\mathbb{K}}\sim k$.

[TABLE]

where $a_{\mu}=u\cdot\partial u_{\mu}$ is the fluid acceleration. Note that $\dot{c}_{s}$ terms arise due to space-time variation of $c_{s}$ via its dependence on $\epsilon(x)$.

II.3 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. (6) and (7) and noting that

[TABLE]

we apply Eq. (30) and keep only the leading order in derivative expansion (i.e., retaining terms of order $\gamma{q}^{2}$ or $k$, but not ${q}k$, consistently with Eq. (5)), to obtain

[TABLE]

where we converted the independent space-time variables from $(x^{+},x^{-})$ to $(x,y)$ by using Eq. (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

[TABLE]

comes from the conversion of $x^{\pm}$ derivatives into $x$, $y$ derivatives, and

[TABLE]

collects terms proportional to $y$, which resulted from the $y$-dependence of $u(x^{\pm})$ and $c_{s}(x^{\pm})$. The last term in Eq. (II.3) follows from the usual procedure of stochastic calculus, by keeping the random noise two-point function in double integrals over the time interval $\delta t$ and using the correlation given in Eq. (18).

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 (II.3) 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. (6) we defined the equal-time correlator of hydrodynamic variables as a function of the mid-point $x$ and the separation vector $y$ as

[TABLE]

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:

[TABLE]

In $G(x+\Delta x,y)$ the orthogonality condition $u(x+\Delta x)\cdot y=0$ is, in general, false, given $u(x)\cdot y=0$ is true: vector $y$ spatial in the frame $u(x)$ is not spatial in $u(x+\Delta 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+\Delta x)$. Defining this boost as $\Lambda^{-1}(\Delta x)$ (inverse for later convenience), i.e.:131313Strictly speaking $\Lambda$ is also a function of $x$ and should be denoted by $\Lambda(\Delta x,x)$. For notational simplicity we drop the $x$ argument.

[TABLE]

we can then define a derivative at “fixed” $y$ as

[TABLE]

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 $\phi_{A}$ contain a scalar and a 4-vector according to Eq. (21)). 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 $\phi_{A}(x+\Delta x)$ and $\phi_{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+\Delta 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 $\phi(x+\Delta x)$ in the same way as $u$ in Eq. (38) before comparing to $\phi(x)$, i.e.,141414Fermi-Walker transport along a world-line is constructed in a similar way, in which case $\Delta x$ is displacement along the particle’s trajectory. In our case $\Delta x$ can point in any direction, not necessarily along $u$.

[TABLE]

With respect to such a derivative, by construction, the flow vector field $u(x)$ is “constant”:

[TABLE]

according to Eq. (38). 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. (38):

[TABLE]

where $\Delta u=u(x+\Delta x)-u(x)$, we obtain the explicit expression for the derivative:

[TABLE]

where the connection associated with the boost created by flow gradients is given by

[TABLE]

Note that this connection is antisymmetric with respect to $\mu\nu$, reminiscent of a spin connection. In a sense, it is a spin connection for a tangent space spanned by hydrodynamic variables $\phi_{A}$ at point $x$. In that sense confluent derivative is a covariant derivative for the connection given by flow gradients in Eq. (44). To unify equations we can extend the range of indices to accommodate the full 5-dimensional space of variables and write

[TABLE]

including the case when $A$ or $B$ is $e$. The corresponding connection is, of course, zero, since $\phi_{e}=c_{s}\delta\epsilon$ is a scalar.

Following the same logic that led us to Eq. (40), 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 $\phi_{A}(x+y/2)$ and $\phi_{B}(x-y/2)$ into the rest frame at the midpoint, $x$, i.e,

[TABLE]

As a result, the confluent correlator, in contrast to Eq. (25), satisfies a simpler orthogonality condition:

[TABLE]

Now combining the three ingredients given by Eqs. (39), (43) and (46) we define the confluent derivative in the following way:

[TABLE]

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 $\Delta x$:

[TABLE]

Another connection, $\mathring{\omega}^{b}_{\mu 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}^{\mu}(x)$, at each point $x$, such that $u\cdot e_{a}=0$. When we keep vector $y$ “fixed” (in the sense described above), its local coordinates $y^{a}=e^{a}_{\mu}(x)y^{\mu}$ may still change due to the rotation of the triad, i.e., not only $e^{a}(x+\Delta x)\neq e^{a}(x)$, but, in general, also $(\Lambda(\Delta x)e)^{a}(x+\Delta x)\neq e^{a}(x)$, in contrast to Eq. (38) (see Appendix A). The last term in Eq. (49) makes sure this change of the basis is corrected for. In other words, we need additional connection to make $e_{a}^{\mu}$ confluently constant (like $u$ already is without additional connection), i.e.,

[TABLE]

so that $\bar{\nabla}_{\lambda}y=\bar{\nabla}_{\lambda}(e_{a}^{\mu}(x)y^{a})=0$. The second term in Eq. (50) accounts for the boost of $e_{a}$ as the one in Eq. (43) and illustrated in Fig. 2. In other words, and this is important for applications, the $x$-deriviative in $\partial_{\mu}{\bar{G}}_{AB}$ in Eq. (49) 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. (50) 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 (50) can be solved for the connection $\mathring{\omega}^{a}_{\lambda b}$ by multiplying by dual basis vector $e^{b}_{\mu}$ such that $e_{c}\cdot e^{b}=\delta_{c}^{b}$:

[TABLE]

where we used $u\cdot e^{b}=u\cdot e_{a}=0$ and the definition of $\bar{\omega}$ connection in Eq. (44). 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 ${\bar{G}}_{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\cdot y=0$ condition, i.e.,

[TABLE]

Thus we arrive at the definition of the confluent Wigner function:

[TABLE]

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}_{\mu}$ in terms of its 3 independent components ${q}_{a}$ as

[TABLE]

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}\to i\partial/\partial{q}_{a}$ and $\partial/\partial y^{b}\to i{q}_{b}$:

[TABLE]

where we also took into account antisymmetry of $\mathring{\omega}^{a}_{\mu b}$. The partial derivative $\partial_{\mu}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\cdot u=0$).

To simplify notations below we shall also use the following expression involving derivatives with respect to components of $q$:

[TABLE]

IV Fluctuation kinetic equations

The two-point functions ${W}_{AB}(x,{q})$ can be viewed as degrees of freedom additional to the hydrodynamic fields $\psi_{A}$ (i.e., $\epsilon$ 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 $\bm{{q}}$ and energy $E=c_{s}|\bm{{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), $\partial_{\mu}\langle{T^{\mu\nu}}\rangle=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}$.

IV.1 Matrix equation for the Wigner function

We return to Eq. (II.3) 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:

[TABLE]

where

[TABLE]

with $[A,B]=AB-BA$ and $\{A,B\}=AB+BA$. In the expression for the anti-commutator, the usual matrix multiplication rules are assumed and the derivates are assumed to act on ${W}$. The matrices $\mathbb{\Omega}^{(s,a)}$ and ${\mathbb{K}}^{(s,a)}$ encode the terms proportional to the gradients of the flow velocity, including the vorticity $\omega_{\mu\nu}$

[TABLE]

and its symmetric partner $\theta_{\mu\nu}$ defined in Eq. (13).

Note that, within the order of approximation we are working, we can further use the ideal hydrodynamic equation $wa_{\mu}=-\partial_{\perp\mu}p$ to eliminate the time-like derivatives of velocity, i.e., $a_{\mu}$, on the right-hand side of Eq. (57). This may be useful for numerical solution of the equations which would require solving for time evolution of $u(x)$ simultaneously.

IV.2 Diagonalization and averaging out fast modes

The matrix ${\mathbb{L}}^{({q})}$ in the right hand side of the kinetic equation Eq. (57) gives the dominant contribution since it is of order of ${q}$ whereas the remaining terms are either order $k$ or $\gamma{q}^{2}$ both of which are assumed to be much smaller than ${q}$ according to our hierarchy of scales in Eq. (5). Therefore it is useful to express the kinetic equation in the basis where ${\mathbb{L}}^{({q})}$ is diagonal. ${\mathbb{L}}^{({q})}$ has five eigenvalues:

[TABLE]

corresponding to 5 eigenvectors $\psi_{\bf A}$ with ${{\bf A}}=+,-,T_{1},T_{2},\parallel$. The eigenvectors form a $5\times 5$ matrix

[TABLE]

where ${\hat{q}}={q}/|{q}|$ is the unit vector along ${q}$ and ${t}^{(1)}$ and ${t}^{(2)}$ are two transverse unit vectors that satisfy

[TABLE]

In other words ${t}^{(1)}$, ${t}^{(2)}$ and ${\hat{q}}$ span the spatial hyperplane orthogonal to $u$,

[TABLE]

The basis vectors in Eq. (61) correspond to the eigenmodes of ideal hydrodynamic equations. Their eigenvalues in Eq. (60) 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. (47). 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 (57) into the diagonal basis of ${\mathbb{L}}^{({q})}$ by the orthogonal transformation $M\rightarrow\psi^{T}M\psi$ 151515Note that since there are derivatives with respect to $x$ and $q$ in Eq. (57), one needs to use $\psi^{T}dM\psi=d(\psi^{T}M\psi)+[\psi^{T}d\psi,\psi^{T}M\psi]$. and express the equation in terms of new variables:

[TABLE]

The modes ${W}_{{\bf A}\parallel}$, ${W}_{\parallel{\bf B}}$, and ${W}_{\parallel\parallel}$ are constrained to vanish by Eq. (47). We can therefore view ${W}_{{\bf A}{\bf B}}$ as effectively a $4\times 4$ matrix. Furthermore, since in the diagonal basis,

[TABLE]

ten of the modes, namely ${W}_{\pm\mp}$, ${W}_{\pm T_{i}}$ and ${W}_{T_{i}\pm}$, 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

[TABLE]

After such averaging only six components of the matrix ${W}$ survive and equations simplify considerably (as noted in Akamatsu:2017 ). As a result we are left with six modes which can be classified into two sound modes ${W}_{\pm\pm}$, two transverse modes ${W}_{T_{1}T_{1}}$ and ${W}_{T_{2}T_{2}}$, and two shear modes ${W}_{T_{1}T_{2}}$ and ${W}_{T_{2}T_{1}}$. The sound modes are completely decoupled and satisfy161616For notational simplicity we denote ${W}_{++}$ and ${W}_{--}$ simply as ${W}_{+}$ and ${W}_{-}$ respectively.

[TABLE]

where

[TABLE]

The confluent derivative of $W_{\pm}$ is defined as follows:

[TABLE]

The transverse and shear modes mix and satisfy $2\times 2$ matrix equation171717Here $\widehat{W}$ represents the $2\times 2$ matrix ${W}_{T_{i}T_{j}}$. Similarly, the $ij$ indices of the $2\times 2$ matrices $\widehat{K}^{ij}$, $\widehat{\Omega}^{ij}$, $\widehat{\omega}_{\mu}^{ij}$, $\widehat{\mathring{\omega}}_{\mu}^{ij}$ and $\widehat{\mathbb{1}}^{ij}=\delta^{ij}$ are suppressed.

[TABLE]

where

[TABLE]

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$:

[TABLE]

The confluent derivative in Eq. (70) also includes additional terms associated with $ij$ indices of $\widehat{W}$ (i.e., of $W_{T_{i}T_{j}}$), which are due to the $x$-dependence of the basis vectors $t^{(i)}$:

[TABLE]

In Appendix B we propose a simple and intuitive choice for the $t^{(i)}$ basis suitable for applications, and compute corresponding connections $\widehat{\omega}_{\mu}^{ij}$ and $\widehat{\mathring{\omega}}_{\mu}^{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}^{\mu\nu}(x)/w$ in Eq. (20), 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}^{\mu\nu}(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.

V Renormalization of first order hydrodynamics

V.1 Short-range singularities and renormalization

The “infinite noise” that one has to introduce via the $\delta$-function in Eq. (18), 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 $\epsilon$ and $u$, as well as the equation of state $p(\epsilon)$ and first-order transport coefficients. This procedure has been discussed by Andreev Andreev:1978 in non-relativistic context, and also recently in boost invariant Bjorken background in Ref.Akamatsu:2018 . 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 $\langle\breve{T}^{\mu\nu}\rangle$ in Eq. (20) via Landau’s matching

[TABLE]

is not given by $u$, i.e., $u_{R}\neq u$ due to fluctuation contributions. Similarly, the energy density $\epsilon_{R}$ in the rest frame $u_{R}$ is different from $\epsilon$. We shall refer to the hydrodynamic variables $(u_{R},\epsilon_{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\mu}(x)$ in Eq. (20). We also need to keep in mind that due to non-linearities in the constraint $\breve{u}\cdot\breve{u}=-1$, $u$ is not properly normalized; namely

[TABLE]

Therefore we find the renormalized fluid velocity as181818 Due to our hierarchy of scales in Eq. (5) the fluctuation contribution is parametrically small, as we also explain below Eq. (83), which justifies the expansion in Eq. (76). It also means that the difference between $u_{R}$ and $u$ is negligible in subleading (viscous) terms and in kinetic equations such as Eq. (57).

[TABLE]

For notational simplicity we will denote $u_{R}$ simply as $u$ in the following.

Expressing $\langle{\breve{T}^{\mu\nu}}\rangle$ in Eq. (20) in terms of $u_{R}$ instead of $u$ using Eq. (76), substituting into Eq. (74) and multiplying both sides by $u_{R\mu}$ we obtain

[TABLE]

In terms of these renormalized quantities, we have

[TABLE]

Note that $p=p(\epsilon)$ 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)\equiv{G}_{AB}(x,0)$ are divergent, i.e., dependent on the wave-vector cutoff $\Lambda$. 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

[TABLE]

The leading singularity is apparent even in static homogeneous equilibrium, since within our coarse-grained resolution ${G}_{AB}^{(0)}(x,y)\sim\delta^{3}(y)$ (i.e., correlation length is negligible) and thus ${G}_{AB}^{(0)}(x,0)$ is undefined. Of course, this is an artifact of neglecting the finiteness of coarse-graining scale ${b}=1/\Lambda$.

The solution of our kinetic equation (57) in equilibrium is simply given by

[TABLE]

where $\Delta_{AB}=\text{diag}(1,\Delta_{\mu\nu})$. Because ${W}^{(0)}$ is $q$-independent, the integral in Eq. (79) is divergent, i.e., cutoff-dependent:

[TABLE]

The corresponding contributions to the energy-momentum tensor in Eq. (78) can be absorbed into a renormalization of the pressure, i.e., the equation of state. The renormalized pressure is then given by

[TABLE]

Written in terms of the renormalized energy density given in Eq. (77), we obtain the renormalized equation of state as

[TABLE]

It is worth emphasizing that, even though $\Lambda$ is an ultraviolet cutoff for the wave-vector of the fluctuating modes, it is still considered small compared to the microscopic scale, i.e. $\Lambda\ll T$ (see Eq. (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. (57)) – 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$:

[TABLE]

where the first and leading term is the equilibrium value (80). Using Eq. (79) we find correspondingly,

[TABLE]

The contributions from the first term to the energy-momentum tensor in Eq. (78) can be absorbed into the renormalization of the equation of state as shown in Eq. (83). The second term in the large ${q}$ expansion in Eq. (84) 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,

[TABLE]

where $\gamma$ represents relaxation constants proportional to the viscosities and $\partial u$ represents velocity gradients. Since only scalar, ${G}_{ee}$, and tensor, ${G}_{\mu\nu}$, components appear in the expansion (78), we focus on those. The phase space integration in Eq. (79) leads to terms linear in $\Lambda$ which can be decomposed into shear (i.e., traceless) and bulk viscous terms:

[TABLE]

where the coefficients $C_{\rm shear}$, $C_{\rm bulk}$, and $C_{ee}$ are given explicitly in the next section (see Eq. (109)). Note that only terms satisfying the orthogonality condition $u^{A}(x)G_{AB}(x)=0$ (according to Eq. (25)) can appear in Eq. (87).

The $\mathcal{O}(\Lambda)$ terms in Eq. (78) due to ${G}^{(1)}_{AB}$ in Eq. (87) can be absorbed by the renormalized transport coefficients (namely shear and bulk viscosities) and pressure. The shear (traceless) term in Eq. (87) can be absorbed by a renormalization of shear viscosity (it has the same form as the shear part of the viscous term $\Pi^{\mu\nu}$ given in Eq. (11)):

[TABLE]

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. (78). Separating the $\Lambda^{3}$ and $\Lambda$ terms explicitly we have

[TABLE]

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 $\theta$ renormalize the bulk viscosity. However, this is not entirely correct, because we also have to take into account that the relation between $\epsilon$ and $\epsilon_{R}$ given by Eq. (77) contains ${G}^{(1)}$ terms, which via the renormalized equation of state $p_{R}(\epsilon_{R})$ in Eq. (83), contribute to what we mean by bulk viscous term separated from the renormalized pressure. Explicitly, we first need to express the bare pressure, $p(\epsilon)$, as a function of the renormalized $\epsilon_{R}$ defined by Eq. (77),

[TABLE]

After inserting this into Eq. (89) and using the renormalized equation of state in Eq. (83), we obtain

[TABLE]

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. (91) can now be absorbed by a renormalization of bulk viscosity191919In principle, the renormalization of equation of state $p(\epsilon)$ leads to a corresponding renormalization of temperature $T(\epsilon)$. However, since $T$ itself appears only in the noise amplitude and therefore only in fluctuation-induced correction, this is a higher order effect.

[TABLE]

It is satisfying to see (and is a non-trivial check) that a conformal symmetry would preserve both the vanishing bulk viscosity $\zeta=0$ and the conformal equation of state $p=\epsilon/3$, according to Eq. (83), under fluctuation corrections, by the virtue of $c_{s}^{2}=1/3$.

Finally, the last term in Eq. (84), $\widetilde{W}_{AB}(x,{q})$, is asymptotically of order $k^{2}/q^{4}$ and does not lead to any large-$q$ divergence in $\widetilde{G}_{AB}(x)$ via integration in Eq. (79). Unlike ${G}^{(1)}_{AB}(x)$, $\widetilde{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 Andreev:1978 ; Kovtun:2003 . 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 $\mathcal{O}(k)$ terms, and the second order $\mathcal{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):

[TABLE]

where we dropped subscript “R” on the right hand side. It should be understood that all quantities in Eq. (93) and in the kinetic equations such as Eq. (57) are renormalized. For example, the pressure, $p(\epsilon)$, in Eq. (93) is given by the physical equation of state, which could be, e.g., obtained from a lattice calculation. The conservation equations for this tensor

[TABLE]

together with the fluctuation kinetic equations (67) and (70) form a closed set of cutoff-independent hydrodynamic equations.202020 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\ll T$) of hydrodynamics Geroch:1995bx , 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.

The necessary components of $\widetilde{G}_{AB}(x)$ can be obtained by solving the kinetic equations (67) and (70) 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. (96), (102), (103) and (105)). Alternatively, with the explicit expressions for ${W}^{(0)}$ and ${W}^{(1)}$ given below, one can substitute (84) into Eqs. (67) and (70) and solve the resulting equations for $\widetilde{W}$ directly. Using Eq. (79) we can then determine

[TABLE]

V.2 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 $\gamma{q}^{2}(W-Tw)$ in Eqs. (67) and (70) should be kept, with the solution given simply by local equilibrium values of fluctuations:

[TABLE]

This is the consequence of the thermal noise in Eq. (18) satisfying fluctuation-dissipation theorem.

In the presence of gradients, the solution deviates from local equilibrium at $x$:

[TABLE]

We can substitute this ansatz into Eq. (67) and Eq. (70) and use the ideal hydrodynamic equations and thermodynamic relations $s=w/T$, $dp=sdT$ to evaluate derivatives of ${W}^{\rm(0)}=Tw$ to leading order in flow gradients:

[TABLE]

As a result we obtain equations for ${{W}^{\rm(neq)}}_{{\bf A}{\bf B}}(x,{q})$. However, as discussed in the previous section, ${W}_{{\bf A}{\bf B}}$ integrated over ${q}$, still produces an ultraviolet divergence (albeit linear in $\Lambda$ and not $\Lambda^{3}$). To isolate this divergence we write

[TABLE]

where we define ${W}_{{\bf A}{\bf B}}^{(1)}(x,{q})$ as the leading term in ${{W}^{\rm(neq)}}_{{\bf A}{\bf B}}(x,{q})$ in the large ${q}$ limit. To find this term we note that the terms in the equation for ${{W}^{\rm(neq)}}_{{\bf A}{\bf B}}$ obtained from Eqs. (67) and (70) by substituting Eq. (97) fall into two classes: the terms proportional to ${{W}^{\rm(neq)}}$ (or its derivatives) and the terms independent of ${{W}^{\rm(neq)}}$. Within each class we identify the leading terms in the limit of ${q}\to\infty$ and require that these terms cancel when we replace ${{W}^{\rm(neq)}}$ with its leading term, $W^{(1)}$. This gives us the following equations:

[TABLE]

which are easily solved as:

[TABLE]

As expected (see Eq. (86)) these terms are of order $\partial u/(\gamma{q}^{2})$ and lead to order $\Lambda$ terms after the ${q}$ integration in Eq. (79).

The remaining terms in ${{W}^{\rm(neq)}}$ in Eq. (99), i.e.,

[TABLE]

are of order $k\partial u/(\gamma{q}^{2})^{2}$ at large ${q}$ and lead to finite $\widetilde{G}(x)$ in Eq. (93) upon ${q}$ integration in Eq. (79). For finite ${q}$ (not satisfying $\gamma{q}^{2}\gg k$) the dependence of $\widetilde{W}$ on ${q}$ can be represented to linear order in $\partial u$, schematically, as

[TABLE]

where $v=\pm c_{s}\hat{q}$ or [math] depending on which mode we are considering. Integration over $d^{3}{q}$ leads to, also schematically, $\widetilde{G}(x)\sim k^{1/2}\partial u/\gamma^{3/2}$. The non-integer power of $k$ represents the fact that $\widetilde{G}$ is a non-local functional of the gradients of $\partial u$. These non-local terms give rise to power-law (in space and/or time) tails in hydrodynamic response Andreev:1978 .212121 Non-linearity in $\partial u$ leads to natural cutoff of the power-law tails, i.e., of non-analiticity at small $k$. Schematically, $\widetilde{G}(x)\sim(k+\partial u)^{1/2}\partial u$.

V.3 Renormalization of transport coefficients

Let us now calculate the renormalization of shear and bulk viscosities using Eq. (102). First of all, we convert back into the original $e,\mu$ basis in order to plug them into the energy momentum tensor, Eq. (78). The conversion is given by:

[TABLE]

In particular, (102) is converted into

[TABLE]

where $\widehat{\Delta}_{\mu\nu}=\sum_{i=1}^{2}{t}^{(i)}_{\mu}{t}^{(i)}_{\nu}=\Delta_{\mu\nu}-{\hat{q}}_{\mu}{\hat{q}}_{\nu}$. With the help of the integrals

[TABLE]

we obtain

[TABLE]

therefore

[TABLE]

Finally inserting these expressions into Eq. (88) and (92) we obtain the renormalized shear and bulk viscosities

[TABLE]

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. Akamatsu:2018 ). The positivity of correction to $\zeta_{R}$ is remarkably non-trivial. It follows from appropriately renormalizing the energy density as well as the equation of state as described by Eqs. (83) and (90). It is satisfying to see that the corrections to viscosities are positive-definite in agreement with the Second Law of Thermodynamics.

VI Phonon interpretation of the fluctuation kinetic equation

VI.1 Phonon kinetic equation

Consider a classical particle whose motion is described in terms of the space-time vector $x^{\mu}$ and 4-momentum ${p}^{\mu}$ 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(\bm{{p}})\equiv c_{s}|\bm{{p}}|$ in the rest frame of the fluid. This can be represented by

[TABLE]

where $u$ is the the 4-velocity of the fluid rest frame $E=c_{s}|{p}_{\perp}|$ and

[TABLE]

The classical action can be then written as

[TABLE]

where $\lambda$ is a Lagrange multiplier. Variation of the action is given by:

[TABLE]

Classical trajectory is then given by equations of motion

[TABLE]

where dot denotes $d/(\lambda d\tau)$ (or one can use reparametrization invariance to set $\lambda\tau$ to equal coordinate time $x^{0}$ in frame $u$) and

[TABLE]

(where we used $\partial{p}_{\perp\nu}/\partial{p}_{\mu}=\Delta^{\mu}_{\nu}$) as well as

[TABLE]

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^{\mu}=u^{\mu}(x)$, as well as $E=E(x,{p}_{\perp})$, which in the case of a phonon means $c_{s}=c_{s}(x)$.

The corresponding Liouville operator acting on a function $\mathcal{{N}}(x,{p})$ is given by

[TABLE]

Note that $\mathcal{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

[TABLE]

where ${N}$ is the usual phase-space distribution function (of 7 variables only). In other words $\mathcal{L}[\delta(F_{+}){N}]=\delta(F_{+})\mathcal{L}[{N}]$.

In order to write the kinetic equation in terms of the distribution function ${N}(x,{p}_{\perp})$ we need to express $x$ derivatives in $\mathcal{L}[N]$ at fixed ${p}$ ($\partial/\partial x^{\mu}$ in Eq. (119)) in terms of $x$ derivatives at fixed ${p}_{\perp}$. These derivatives are not the same because the relationship between ${p}$ and ${p}_{\perp}$ depends on $x$ (via $u(x)$ in Eq. (113)). One finds

[TABLE]

where we denoted by $\bar{\nabla}_{\mu}$ the $x$ derivative at ${p}_{\perp}$ fixed 222222A more explicit definition involves projections ${p}_{a}$ of ${p}_{\perp}$ on the local triad ${p}_{\perp\mu}=e_{\mu}^{a}{p}_{a}$, in terms of which $\bar{\nabla}_{\mu}{N}=\partial_{\mu}{N}+\mathring{\omega}^{a}_{\mu b}{p}_{a}\partial{N}/\partial{p}_{b}$ (cf. Eqs. (50), (51) and (55)). The projections ${p}^{a}$ are kept fixed while taking $x$ derivative, and connection term accounts for the rotation of the basis triad $e^{a}(x)$ which changes ${p}_{\perp}$ while $p_{a}$ is fixed. Similarly, ${p}_{\perp}$ derivatives at fixed $x$ are more explicitly written as $\partial/\partial{p}_{\perp\mu}=e^{\mu}_{a}\,\partial/\partial{p}_{a}$ (cf. Eq. (56)).. Corresponingly, the last term in Eq. (118) should be written as

[TABLE]

Similarly, the ${p}$ derivatives at fixed $x$ should be expressed as ${p}_{\perp}$ derivatives

[TABLE]

Substituting Eqs. (116), (118), (122), (121) and (123) into Eq. (119) we find

[TABLE]

Finally, using $\partial_{\mu}{p}_{\perp\nu}=-E\partial_{\mu}u_{\nu}+u_{\nu}\partial_{\mu}({p}\cdot u)$, we can write the Liouville operator as

[TABLE]

The expression in the square brakets is (the negative of) the force acting on the phonon. 232323One can also obtain this expression by taking the spatial projection of the rate of change of ${p}_{\perp}$, i.e., the force is $\Delta_{\mu}^{\nu}\dot{p}_{\perp\nu}$, and using equations of motion (116) and (118) together with Eq. (113). 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 $\omega_{\mu\nu}$, respectively. The force $-{p}_{\perp\nu}\partial_{\perp\mu}u^{\nu}$ is easier to understand by considering isotropic Hubble-like expansion, i.e., such that $\partial_{\perp\mu}u^{\nu}=H\Delta_{\mu}^{\nu}$, where $H$ is the rate of expansion (Hubble constant). This term then describes the rescaling of the momentum ${p}_{\perp}$ (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:

[TABLE]

Remarkably, upon changing the notation for the phonon momentum

[TABLE]

the Liouville operator in Eq. (125) with $E=c_{s}|{p}_{\perp}|$ is identical to the one in Eq. (67) obtained using completely different (but apparently complementary) considerations. The two signs in front of $c_{s}$ in Eq. (67) correspond to positive and negative frequency sound waves, or positive/negative energy solutions of the condition

[TABLE]

where $F_{\pm}=(p\cdot u)\pm E$ and the positive energy solution is given by $F_{+}=0$ in Eq. (112).

Curiously, for linear dispersion, $E=c_{s}|p_{\perp}|$, the condition in Eq. (128) can be written as $g^{\mu\nu}{p}_{\mu}{p}_{\nu}=0$ using flow induced effective “metric tensor” $g^{\mu\nu}=-u^{\mu}u^{\nu}+c_{s}^{2}\Delta^{\mu\nu}$. Since $d(F_{+}F_{-})=F_{-}dF_{+}+F_{+}dF_{-}$ and $\delta(F)=\delta(F_{+})/F_{-}+\delta(F_{-})/F_{+}$, we see that the equations of motion localized on the $F_{+}=0$ surface are given by Eqs. (116) and (118) up to rescaling of proper time. On the other hand, the equations of motion with the constraint $F_{+}F_{-}=0$ are given by

[TABLE]

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. (125) can be viewed as “gravitational” forces.

Perhaps even more remarkably than the matching of the Liouville operators in Eqs. (125) and (67), the identification

[TABLE]

leads to nontrivial cancellation of the whole second line in Eq. (67) (i.e., of all terms proportional to the background gradients $\theta_{\mu\nu}$ and $a^{\mu}$ times ${W}_{\pm}$) leaving simply the relaxation term in Eq. (67):

[TABLE]

where $\mathcal{L}_{\pm}$ are different by the sign in front of $c_{s}$ in Eq. (67). Note that the equilibrium value of $N_{\pm}$, ${N}_{\pm}^{(0)}={W}_{\pm}^{(0)}/(c_{s}|{p}_{\perp}|w)$, equals $T/E$ as expected for the low-energy limit of the phonon Bose-Einstein distribution function.

In contrast to Eqs. (67) for longitudinal modes which reduces to a simple form Eq. (131) upon rescaling given by Eq. (130), Eq. (70) 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.

VI.2 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 andreev1970twoliquid .

Let us start with the expression for the stress tensor for one particle moving along a trajectory specified by $x(\tau)$:

[TABLE]

This means for a gas of such particles and holes with distribution functions $N_{\pm}(x,p)$ we have (’s’ for ’sound’):

[TABLE]

(Minus here is because a hole is the absence of a particle). Using equation of motion $\dot{x}=u\pm v$ (Eq. (116)) with $\pm$ for particles/holes we obtain:

[TABLE]

Using now $E=c_{s}|p_{\perp}|$ for the phonon, we find:

[TABLE]

Comparing to Eq. (105) (neglecting transverse modes) and identifying $N_{\pm}={W}_{\pm}/(c_{s}|p_{\perp}|w)$ (as in Eq. (130)) we can write:

[TABLE]

We can recognize the last two terms as the last two terms in Eq. (20). The first term contributes to the shift/renormalization of the energy density (see Eq. (77)).

Phonons contribute to pressure, i.e., renormalize equation of state. The equilibrium pressure of the phonon gas equals

[TABLE]

Together with the contribution from the shift of the energy due to phonons, $-c_{s}^{2}{{G}^{(0)}}^{\mu}_{\mu}(x)/w$, this reproduces the last term before the second equal sign in Eq. (83).

VII Conclusions and outlook

We have derived a set of equations which describe coupled evolution of hydrodynamic variables $\epsilon$ and $u$ and the correlation functions of their fluctuations, collectively denoted by $\phi_{A}$. The correlation functions are expressed as Wigner transform $W_{AB}(x,{q})$ of the equal-time correlator $\langle\phi_{A}\phi_{B}\rangle$.

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. (93) 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 Andreev:1978 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, $\ell_{\rm eq}$, is of order $\sqrt{\gamma L/c_{s}}\ll L$. The corresponding “hard” momentum ${q}\sim 1/\ell_{\rm 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. (67) vs Eq. (125).

The originial approach by Andreev postulated Hamiltonian kinetic equations for a phonon distribution function without derivation Andreev:1978 . We derive these equations directly from stochastic hydrodynamics and confirm that the collision/relaxation term has the simplest form assumed in Ref. Andreev:1978 and does not depend on the gradients of the flow. This result emerges after nontrivial cancellations (Eq. (131) vs Eq. (67)), 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. (70)). These gradient terms were assumed to be absent in Ref. Andreev:1978 . 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. Akamatsu:2017 ). In particular, the effects of vorticity, absent in the Bjorken flow, but important in heavy-ion collisions STAR:2017ckg , 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 Stephanov:2018hydro+ . 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 242424Special cases of static and boost invariant backgrounds for a conformal fluid with conserved charge have been discussed in Ref. Martinez:2018 .. Such a generalization in essential for the hydrodynamic modeling of fluctuations and its effects in the the beam-energy scan program Aggarwal:2010cw 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. Crossley:2015evo ; Haehl:2018lcu for reviews). However, the practical usage is so far mostly limited to correlation functions in simple backgrounds, such 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.).

Acknowledgements.

We thank Boris Spivak for drawing our attention to Ref. andreev1970twoliquid . We thank Mauricio Martinez and Derek Teaney for helpful discussions. This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration and grant Nos. DE-FG0201ER41195 and DE-SC0018209.

Appendix A A local confluent triad

In order to describe the separation vector $y$ (for example, to enable numerical solution of the equations for fluctuation correlators) we need to introduce a basis triad $e^{a}_{\mu}(x)$ for the tangent plane orthoginal to $u(x)$ at each point $x$. The basis is arbitrary and here we shall propose a simple and intuitive choice of $e_{a}(x)$. We choose a (lab) frame $\mathring{u}$ and a fixed triad ($a=1,2,3$) satisfying $\mathring{e}_{a}\cdot\mathring{e}^{b}=\delta^{b}_{a}$ and $\mathring{e}_{a}\cdot\mathring{u}=\mathring{e}^{b}\cdot\mathring{u}=0$. For simplicity we shall consider an orthogonal triad, equivalent to its dual, $e^{a}=e_{a}$.

We can then define $e^{a}(x)$ by a finite boost from $\mathring{u}$ to $u(x)$. The resulting triad vectors at point $x$ are given by explicit algebraic formulas:

[TABLE]

One can check that $e_{a}\cdot u=0$ and $e^{a}\cdot e_{b}=\delta^{a}_{b}$.

Corresponding spin connection is given by Eq. (51)

[TABLE]

In terms of the confluent connection defined in Eq. (44) one can express spin connection as

[TABLE]

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 $\mathring{\omega}$ vanish. This requires integrability of Eq. (50) with $\mathring{\omega}=0$, which means that the change of vector $e_{a}$ obtained by integrating Eq. (50) with $\mathring{\omega}=0$ between two points should not depend on the path, i.e.,

[TABLE]

Using Stokes theorem we see that this is possible if curvature associated with connection $\bar{\omega}^{\mu}_{\lambda\nu}$ vanishes. Using Eq. (44), we find:

[TABLE]

One might say that $\bar{R}_{\alpha\beta}{}^{\mu}{}_{\nu}=0$ means $\bar{\omega}^{\mu}_{\lambda\nu}$ is a “pure gauge” connection.

A nontrivial example of flow with $\bar{\omega}^{\mu}_{\lambda\nu}\neq 0$ but $\bar{R}_{\alpha\beta}{}^{\mu}{}_{\nu}=0$ is the Bjorken flow. In this case our proposed choice of $e_{a}(x)$ in Eq. (138) provides a rotationless (i.e., $\mathring{\omega}=0$) triad field.

Appendix B A basis in the space orthogonal to $\hat{q}$ and $u$ and

monopole connection

A basis in the space orthogonal to $\hat{q}$ and $u$ (cf. Eq. (63)) can be obtained easily by rotating the local confluent basis $e^{a}$ in such a way that one of the vectors, say $e_{3}$, lines up with $\hat{q}$. The result is given by

[TABLE]

satisfying ${t}^{(i)}\cdot{t}^{(j)}=\delta^{ij}$ and ${t}^{(i)}\cdot\hat{q}_{\perp}={t}^{(i)}\cdot u=0$.

Since $t^{(i)}$ depends on $x$ (to maintain $u(x)\cdot t^{(i)}=0$) as well as on ${q}$ (to keep ${\hat{q}}\cdot{t}^{(i)}=0$), there are two types of connections in Eq. (70) defined by Eqs. (73) and Eqs. (72). Applying these definitions to our choice of ${t}^{(i)}$ in Eq. (143) we find for the $x$-derivative connection in Eq. (73)

[TABLE]

where the connection $\mathring{\omega}^{a}_{\lambda b}$ is defined by Eq. (51).

For the $q$-derivative connection, using definition in Eq. (72), one obtains

[TABLE]

The last expression can be easily recognized as the connection describing a monopole at ${q}=0$ and Dirac string along $-e_{3}$. The corresponding curvature 252525 Because the space spanned by $t^{(i)}$ is two-dimensional the connection is abelian, i.e., $[\widehat{\omega}_{\mu},\widehat{\omega}_{\nu}]=0$.

[TABLE]

is the field of a monopole with charge 1 (twice the amount of Berry curvature monopole charge for spin-$1/2$ fermion). The singularity at ${q}=0$ is associated with the ambiguity of ${\hat{q}}$ at ${q}=0$.

Appendix C Comparison to Bjorken flow results

The purpose of this appendix is to compare our equations with the ones for a particular case of Bjorken flow derived in Ref. Akamatsu:2018 .

The first observation we need to make is that the definition of the equal-time correlator in Ref. Akamatsu:2018 is subtly different. The Bjorken flow allows us to define a hypersurface globally which is orthogonal to the flow 4-vector $u(x)$ at each point: the constant proper-time surface $\tau={\rm const}$. It is then natural to define “equal time” correlator in such a way that points $x^{\pm}$ lie on the same proper-time hypersurface as $x$. The difference with our definition is subtle because our equal-time hyperplane is tangential to the equal-$\tau$ hypersurface at point $x$ and the difference is of order $y^{2}$, due to the curvature of the surface. This difference does lead to a subtle change in the last term in Eq. (57), which is necessary to make this equation agree with Ref. Akamatsu:2018 .

To describe this in more detail, let us consider a definition of the equal-time correlator which is slightly different from ours, but will coinside with $\tau={\rm const}$ for Bjorken flow. It is possible to define a hypersurface orthogonal to flow if the flow is conservative, i.e., $u^{\mu}=\partial^{\mu}\tau$ (as is the case for the Bjorken flow, for example). In general it is not possible, however, one can perform a Helmholz decomposition into conservative (potential) and purely vortical flow: $u^{\mu}=\partial^{\mu}\tau+v^{\mu}$, where $\partial\cdot v=0$ (see Fig. 3 for illustration). We will not be interested in doing this globally since we only need to describe the surface near a given point $x$ to quadratic order in $y$. Thus we Taylor expand $u$ to linear order in $\Delta x$:

[TABLE]

The last term is purely vortical, while the first two terms are potential, i.e.,

[TABLE]

We can then define equal-time correlator in such a way that points $x$ and $x^{\pm}=x\pm y/2$ lie on the same curved surface $\tau={\rm const}$. Using 3-dimensional vector $y$ in the tangent plane to $u(x)$ to parameterize points on such a surface, we can write explicitly

[TABLE]

where $y\cdot u(x)=0$ and the last term describes the curvature of the surface. Using this definition of $x^{\pm}$ instead of Eq. (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 ${\mathbb{L}}^{(y)}$ defined in Eq. (34) by

[TABLE]

where

[TABLE]

As a consequence, Eq. (II.3) and (57) are modified and take the form

[TABLE]

and

[TABLE]

respectively, where

[TABLE]

Note, that the only change compared to Eq. (57) is in the double commutator term.262626This is consistent with the fact that upon diagonalization and time-averaging over faster modes this term drops completely. Indeed, the scale of time-averaging, ${b}_{t}$ is much longer than the typical time-like separation between the plane tangent to $u$ and the $\tau={\rm const}$ defined by Eq. (149), which is of order $(\partial u)y^{2}\sim k/q^{2}\ll 1/q$, compared to ${b}_{t}\gg 1/q$ according to Eq. (66).

For the Bjorken flow, $a_{\mu}=\omega_{\mu\nu}=0$, and all the terms on the second line in Eq. (155) vanish.

To complete the comparison, for the boost-invariant flow, we perform the coordinate transformation from $(t,x,y,z)$ to $(\tau,x,y,Y)$ given by $t=\tau\cosh Y,~{}x=x,~{}y=y,~{}z=\tau\sinh Y$, where $\tau$ is the proper time and $Y$ is the space-time rapidity. One can easily check that for the Bjorken flow $u\cdot\bar{\nabla}=\partial_{\tau}$, $\theta=1/\tau$, $a_{\mu}=\omega_{\mu\nu}=0$. Thus Eq. (LABEL:eq:kinetic_eq_conserv) is reduced to

[TABLE]

Since ${q}_{Y}=\tau{q}_{z}$ where ${q}_{Y}$ is the wave vector conjugate to $Y$, we define ${W}_{B}(x;{q}_{Y})={W}(x;{q}_{z})/\tau$ to take into account the change in the measure of the momentum integration. Using

[TABLE]

we obtain

[TABLE]

where the last two terms in Eq. (155) were eliminated by the momentum rescaling. Similarly, one can check that our Eqs. (67) and (70), rewritten in terms of $W_{B}$, will reduce to Eq. (A7) in Ref. Akamatsu:2018 exactly.

Appendix D Notations

• $\bar{\nabla}_{\mu}$ – confluent derivative – Eqs. (45), (55), (73);
• $\gamma$ – generic relaxation constant used in order-of-magnitude expressions;
• $\gamma_{\eta}$, $\gamma_{\zeta}$, $\gamma_{L}$ – shear, bulk and longitudinal (sound) relaxation constants Eqs. (28) and (68);
• $\Delta^{\mu\nu}$ – projector on hyperplane orthogonal to $u$ – Eq. (12);
• $\delta e$, $\delta g^{\mu}$ – scaled fluctuations of energy density (times $c_{s}$) and velocity (times $w$) – Eq. (21);
• $\epsilon$, $u$ – local averaged energy density and fluid velocity – Eq. (16);
• $\breve{\epsilon}$, $\breve{u}$ – fluctuating energy density and fluid velocity – Eq. (16);
• $\theta^{\mu\nu}$ – symmetrized velocity gradients – Eq. (13);
• $\Lambda$ – with no argument – wave-vector cutoff, $\Lambda=1/b$;
• $\Lambda(\Delta x)$ – Lorentz boost bringing fluid at point $x+\Delta x$ to rest with respect to fluid at point $x$ – Eq. (38);
• $\phi_{A}$ – ($A=e,0,1,2,3$) the set of fluctuations of hydrodynamic variables – Eq. (21);
• $\psi_{A}$ – the set of hydrodynamic variables, e.g., $(\epsilon,u^{\mu})$;
• $\omega^{\mu\nu}$ – vorticity – Eq. (59);
• $\bar{\omega}^{\nu}_{\lambda\mu}$ – confluent connection – Eq. (44);
• $\mathring{\omega}^{a}_{\mu b}$ – spin connection for local triad $e^{a}$ – Eq. (51);
• $\widehat{\omega}^{ij}_{\mu}$ – momentum space spin connection for diad $t^{(i)}$ – Eq. (72);
• $\widehat{\mathring{\omega}}^{ij}_{\mu}$ – coordinate space spin connection for diad $t^{(i)}$ – Eq. (73);
• $a^{\mu}$ – local acceleration $a^{\mu}=(u\cdot\partial)u^{\mu}$;
• ${b}$ – coarse-graining scale – Eq. (1);
• ${b}_{t}$ – temporal coarse-graining scale – Eq. (66);
• $\dot{c}_{s}$ – logarithmic rate of dependence of sound speed on entropy – Eq. (26);
• $\partial_{\perp}$ – partial derivative projected on hyperplane orthogonal to 4-velocity – Eq. (14);
• $\partial_{\mu}{W}_{AB}(x,q)$ – partial $x$-derivative at fixed $q^{a}=e^{a}(x)\cdot q$ – Eq. (55);
• $E$ or $E({p}_{\perp})$ – phonon energy – Eq. (112);
• $e^{a}_{\mu}$ – ($a=1,2,3$) local triad basis vector orthogonal to $u(x)$ – Eq. (50), Appendix A;
• $F_{\pm}$ – Eq. (128);
• ${G}_{AB}(x)\equiv{G}_{AB}(x,0)$ – correlator at coincident points – Eq. (22);
• $\widetilde{G}_{AB}(x)$ – finite part of ${G}_{AB}(x)$ – Eq. (95);
• ${G}_{AB}(x,y)$ – equal-time correlator – Eq. (6);
• ${\bar{G}}_{AB}(x,y)$ – confluent equal-time correlator – Eq. (46);
• $k$ – wave vector Fourier conjugate to midpoint vector $x$;
• ${N}(x,p)$ – phonon phase-space distribution function – Eq. (120);
• ${N}_{\pm}(x,p)$ – phase-space distribution function for positive/negative frequency phonons;
• ${p}$ – phonon momentum (not to be confused with pressure $p(\epsilon)$) – Section VI;
• $p\equiv p(\epsilon)$ – pressure at average energy density $\epsilon$;
• ${q}$ – wave vector Fourier conjugate to separation vector $y$, also $q\cdot u=0$ – Section III;
• $T$ – temperature (local value at $\epsilon(x)$);
• $t^{(i)}_{\mu}$ – ($i=1,2$) momentum space diad vector orthogonal to ${q}^{\mu}$ – Eq. (62), Appendix B;
• $u^{\mu}(x)$ – local averaged 4-velocity at point $x$;
• $w\equiv\epsilon+p(\epsilon)$ – enthalpy density at average energy density;
• $W_{AB}(x,{q})$ – Wigner function, Wigner/Fourier transform of ${\bar{G}}_{AB}(x,y)$ – Eq. (53);
• $\widetilde{W}_{AB}(x,{q})$ – $W_{AB}(x,{q})$ after subtraction of leading and subleading large-$q$ terms – Eq. (103).
• $W_{{\bf A}{\bf B}}(x,{q})$ – Wigner function $W_{AB}(x,{q})$ in the basis of ideal hydrodynamics modes – Eqs. (61), (64);
• $W_{\pm}\equiv W_{\pm\pm}$ – Wigner functions for positive/negative frequency sound modes – Eqs. (61), (64);
• $\widehat{W}^{ij}$ – $2\times 2$ matrix of Wigner functions for transverse modes – Eq. (70);
• $x$ – in a two-point correlator – the midpoint space-time vector;
• $x^{\pm}$ – arguments of the 2-point correlator;
• $y$ – in an equal-time two-point correlator – the separation vector, constrained by $u(x)\cdot y=0$.