Fluctuation dynamics in a relativistic fluid with a critical point

To describe dynamics of bulk and fluctuations near the QCD critical point we develop general relativistic fluctuation formalism for a fluid carrying baryon charge. Feedback of fluctuations modifies hydrodynamic coefficients including bulk viscosity and conductivity and introduces nonlocal and non-instantaneous terms in constitutive equations. We perform necessary ultraviolet (short-distance) renormalization to obtain cutoff independent deterministic equations suitable for numerical implementation. We use the equations to calculate the universal non-analytic small-frequency dependence of transport coefficients due to fluctuations (long-time tails). Focusing on the critical mode we show how this general formalism matches existing Hydro+ description of fluctuations near the QCD critical point and nontrivially extends it inside and outside of the critical region.


I. INTRODUCTION
The classic subject of relativistic hydrodynamics [1] has experienced a renaissance of interest in recent years [2,3]. The interest is driven in large part by the progress in heavy-ion collision experiments which allow us to create and study droplets of hot and dense matter governed by physics of strong interaction described by Quantum Chromodynamics. The increasing body of experimental evidence that relativistic hydrodynamics is describing the evolution of the expanding fireball created in these collisions motivates technical developments as well as a closer look at many fundamental theoretical concepts in hydrodynamics.
The subject of hydrodynamic fluctuations is particularly relevant to heavy-ion collisions. The system size L is not astronomically large compared to the typical microscopic scale, mic , (factor 10 at most is a typical scale separation). 1 As a result, fluctuations are large enough to be easily observable in experiments. In addition, since the leading corrections to hydrodynamics are due to the nonlinear feedback of fluctuations, we cannot afford to neglect them -a luxury one is used to in ordinary fluid dynamics. Fluctuations are even more important when enhanced by critical phenomena.
From the modern point of view, hydrodynamics is a systematic expansion in spatial gradients. More precisely, it is the expansion of constitutive equations for stress tensor (and conserved current). The expansion parameter is the ratio of a typical hydrodynamic wave-number k = 1/L to a microscopic scale, say temperature T , or inverse scattering length, or, generically, 1/ mic . In this view, the ideal, non-dissipative (i.e., reversible) hydrodynamics is the truncation of this expansion at lowest (zeroth) order. At first order in gradients (i.e., at order k 1 or, more precisely, (k mic ) 1 ) one recovers standard Landau-Lifshitz or Navier-Stokes hydrodynamics. It is the following order in this expansion that concerns us here. That order is not k 2 , but rather is k 3/2 (or k d/2 in d-dimensions). Such non-analytic in k and, therefore, nonlocal contributions come from fluctuations in hydrodynamics. Thus it is essential to understand the physics of hydrodynamic fluctuations to faithfully describe physics of heavy-ion collisions. 2 Furthermore, in addition to modifying hydrodynamic equations by effectively nonlocal contributions, the fluctuations themselves are measured in heavy-ion collision experiments. In particular, one of the most fundamental questions these experiments aim to answer is the existence and location of the critical point on the QCD phase diagram [4,5]. The signature of this phenomenon is a certain non-monotonous behavior of event-by-event fluctuation measures when the parameters of the collision (such as center of mass energy) is varied in order to "scan" QCD phase diagram [6,7]. This non-monotonous behavior is driven by critical phenomena and thus predictable without being able to determine QCD equation of state at finite density (still an unsolved theoretical problem).
The existing predictions for critical behavior rely significantly on the assumption of local thermal equilibrium. However, near the critical point, the equilibrium is increasingly difficult to achieve due to critical slowing down and a finiteness of expansion time. Essentially, this limitation determines the magnitude of the observable signatures of the critical point [7,8]. Therefore the ability to describe the dynamical evolution of fluctuations during the fireball evolution, in particular, in the proximity the critical point is crucial. The goal of this paper is to provide such a description.
One of the recent advances towards this goal has been the introduction of Hydro+ in Ref. [9], with a recent numerical implementation in a simplified setup reported in Ref. [10]. Focusing on the mode responsible for the critical slowing down, identifying it with the fluctuation correlator of the slowest hydrodynamic mode, the authors of Ref. [9] proposed the evolution equation which describes the relaxation of this non-hydrodynamic mode to equilibrium. Extending hydrodynamics by addition of such a mode one is then able to broaden the range of applicability of hydrodynamics near the critical point and describe the dominant mode of critical fluctuations at the same time. The crucial ingredient of this formalism is a nonequilibrium entropy of fluctuations derived in Ref. [9].
We approach this problem from a different direction. We start with the general formalism of relativistic hydrodynamic fluctuations introduced earlier in Ref. [11] for neutral (chargeless) fluids and extend it to include a crucial ingredient -baryon charge density. QCD critical point, if it exists, is located at finite baryon density. The approach we pursue, in which the two-point correlators of hydrodynamic variables play the role of additional non-hydrodynamic variables, has been introduced and developed recently in the context of heavy-ion collisions, but limited to special types of flow such as longitudinal boost-invariant expansion in Refs. [12][13][14]. In a more general but nonrelativistic case this approach was pioneered by Andreev in the 1970's [15]. The approach is often referred to as 'hydro-kinetic' to acknowldge the similarity between the two-point correlators and the distribution functions in kinetic theory. In particular, the dynamics of the correlators of the pressure fluctuations is essentially equivalent to the kinetics of the phonon gas. This physically intuitive picture was the original source of this formalism [12,16] and was rigorously derived in general relativistic context in Ref. [11].
The hydro-kinetic approach should be also contrasted with the traditional stochastic hydrodynamics where the noise is introduced into hydrodynamic equations [17]. From this point of view, the 'hydro-kinetic' approach could also be called 'deterministic', as it replaces stochastic equations with deterministic equations for the evolution of correlation functions. Of course, the two approaches solve the same system of stochastic equations, but in complementary ways. The advantage of the deterministic approach is that it allows one to deal with the problem of the "infinite noise": the noise amplitude needs to become infinitely large as the hydrodynamic cell size is sent to zero, even though the physical effect of the noise is finite due to its averaging out in a medium whose properties vary slowly in space and time. The effect of the infinite (or more precisely cutoff dependent) noise can be absorbed into renormalization of hydrodynamic equationsa procedure which can be performed analytically in the deterministic approach. This avoids having to deal with numerical cancellations which would otherwise be necessary in a direct implementation of stochastic equations.
Near the critical point the deterministic approach we develop here, although different from Hydro+ in Ref. [9], nevertheless leads to the description of fluctuations in terms of two-point correlators as in Hydro+. In this paper we verify that in the limit of large correlation length the two approaches exactly match. This is a nontrivial check of the validity of both approaches. Furthermore, since the deterministic approach is more general it allows us to extend the Hydro+ approach both closer to the critical point and further away from the critical point to describe also ordinary, noncritical fluctuations.
The paper is organized as follows. In Section II we start from stochastic hydrodynamics and derive linearized stochastic equations for fluctuations of hydrodynamic variables. In Section III we use this result to derive deterministic evolution equations for two-point correlators of hydrodynamic variables. These equations bear resemblance to evolution equations in kinetic theory. In Section IV we expand the stochastic hydrodynamic equations again, now further, to second order in the fluctuations. Upon averaging over noise, we obtain the equations for averages of hydrodynamic variables (one-point functions). These equations, due to nonlinearities, now contain the contributions of two-point functions. These contributions lead to renormalization of "bare" hydrodynamics equations, i.e., they change the "bare" equation of state and "bare" transport coefficients into physical quantities. All cutoff dependence is absorbed at this stage. The remaining contributions are nonlocal and are known as long-time tails. As an example, we work out explicitly the non-analytic small-frequency dependence of transport coefficients. In Section V we study the behavior of the equations we derived near the critical point, perform comparison with Hydro+ and propose an extension to shorter time/length scales which we refer to as Hydro++.

II. STOCHASTIC HYDRODYNAMICS AND FLUCTUATIONS
A. Hydrodynamics with noise The starting point of our analysis is the relativistic hydrodynamics with stochastic noise which drives the thermal fluctuations of hydrodynamic variables. The amplitude of the noise is proportional to the dissipative transport coefficients as a consequence of the fluctuation-dissipation theorem. For macroscopically large systems we are considering, fluctuations are small because they originate from the noise at microscopically short scales and average out over macroscopic distances. This allows us to treat the effects of these fluctuations systematically in a power expansion. We shall discuss the corresponding power-counting in Section VI. In particular, in this paper we will focus on the quadratic order in fluctuations and neglect cubic and higher order. In this section, we lay out this basic set-up, up to the quadratic order of fluctuations that we are working with. This is a generalization of our previous work for fluids without conserved charges [11] to the case of relativistic plasma with a conserved charge.
Stochastic hydrodynamics with a conserved U(1) charge is defined by the conservation equations Here, we follow the conventions in Ref. [11], and distinguish stochastic quantities by the breve accent˘. Defining the stochastic hydrodynamic variables (˘ ,n,ȗ µ ) by the Landau's conditions we can writeT µν andJ µ asT in terms of the "bare" constitutive relations where ∆ µν ≡ g µν + u µ u ν is the standard space-like projection operator. The function p( , n) is the pressure, given by the equation of state. As usual, p, in terms of entropy s: where temperature T and chemical potential µ are defined via derivatives of s (the first law of thermodynamics): Throughout the paper we will also use the enthalpy density, and, most importantly, entropy per charge ratio m ≡ s n . (2.8) The constitutive equations in Eqs. (2.4) are organized as an expansion in powers of spatial gradients. The terms first-order in gradients are given by and (2.11) The shear and bulk viscosities are denoted by η and ζ, and the charge conductivity is denoted by λ 3 . The fluctuations are sourced by random noise terms (S µν ,Ȋ µ ) that are sampled over a Gaussian distribution with an amplitude determined by the fluctuation-dissipation theorem, where λ, T, η and ζ here assumed be functions of averaged thermodynamic variables, such as energy density and number density. Eqs. (2.1) together with constitutive equations (2.4) determine evolution of stochastic variables˘ ,n andȗ. Since we have the freedom to choose an independent pair of scalar variables arbitrarily, we use this freedom to keep our calculations and resulting equations relatively simple. We find the following set of variables particularly convenient:m ≡ m(˘ ,n) andp ≡ p(˘ ,n), (2.13) where m( , n) and p( , n) are the entropy per charge and the pressure, expressed as functions of the energy and charge densities, defined in Eqs. (2.8) and (2.5). This choice simplifies our calculations because the fluctuations of m and p are statistically independent in equilibrium and correspond to two eigenmodes of linearized ideal hydrodynamic equations. We shall denote the ensemble averages of these variables by simply removing the accent, i.e., Having defined variables m and p as average values (one-point functions) of primary variables in Eqs. (2.14) we shall now define other deterministic variables, which appear in our equations, such as and n as functions of m and p obtained via equation of state: ≡ (m, p), n ≡ n(m, p).
(2.15) Note that, due to nonlinearities in these relationships, = ˘ and n = n . In order to describe the evolution of these deterministic quantities we shall perform the ensemble average on the stochastic equations. Although this eliminates the noise terms, because of the nonlinearities in the constitutive equations the averaged equations cannot be simply obtained by substituting stochastic variables by their averages. We shall describe the effect of these nonlinearities on the evolution of average values (i.e., one-point functions in Eq. (2.14)) in Section IV. These effects, to lowest order in the magnitude of the fluctuations, are given in terms of the two-point functions. Our goal in Section III will be to derive evolution equation for these correlators. We should also keep in mind that these two-point functions are of interest in their own right, since they describe the magnitude of the fluctuations and correlations which, in heavy-ion collisions, are measurable.

B. Linearized equations
To obtain equations for the two-point functions we first begin by writing stochastic hydrodynamic equations linearized in deviations of the stochastic variables from their average values in Eq. (2.14): m = m + δm,p = p + δp,ȗ µ = ȗ µ + δu µ . (2.16) To linear order, the fluctuations ofm andp are simply related to fluctuations of˘ = (m,p) andn = n(m,p) by a linear transformation with coefficients given by thermodynamic derivatives. We shall use the following intuitive short-hand notations for these derivatives: d = m dm + p dp, dn = n m dm + n p dp (2.17) whose exact definitions are given in Appendix A. Similarly, we find it useful to express the fluctuations of the thermodynamic functionα = α(˘ ,n) defined in Eq. (2.6) in terms of δm and δp and define corresponding coefficients: Of course, due to nonlinearities in the equation of state the relationship between fluctuations of (˘ ,n,α) and (m,p) is nonlinear, and we shall deal with this in Section IV where we consider the second order terms in the fluctuation expansion. Now we are ready to expand the constitutive equations to linear order in fluctuations: The equations of motion for both the background and the fluctuations are obtained by substituting Eq. (2.19) into Eq. (2.1). By definition, Eq. (2.16), one-point averages of fluctuations vanish, δm = δp = δu = 0. Therefore, upon averaging the equations of motion, ∂ µT µν = ∂ µJ µ = 0, we obtain At leading order in gradients, this gives us equations of ideal hydrodynamics, We also use the thermodynamic relation, and express our equations in terms of gradients of p and α. With the help of above expressions, we find, after some amount of algebra, the following equations of motion for our fluctuating variables: (2.24) We also introduced a useful notation "dot" for the operation defined as: for a given thermodynamic quantity X. Note that since this operation is a logarithmic derivative it satisfies This operator appears in our equations because to leading order (ideal hydrodynamics), (u · ∂) log X = −Ẋθ and u·∂ m = 0. The quantitiesċ s and˙ m involve third order thermodynamic derivatives (i.e., third derivatives of a s( , n)), while the quantityṪ is a second order thermodynamic derivative (i.e., second derivatives of a s( , n)). Indeed, among all the second order thermodynamic quantities defined in Eqs. (2.17), (2.18) and appearing thereafter, only three are independent. In other words, all second order derivatives can be expressed by at most three independent quantities. We find that a convenient choice at this stage of the calculation is α m , α p and c 2 s ≡ (∂p/∂ ) m . However, at a later stage, we shall instead choose another set, c s , c p andṪ (together with independent third order thermodynamic derivatives such asċ s ,ċ p ), to express our final results. The expressions relating those quantities belonging to different sets are (2.27) Using the above relations derived in Appendix A we can express all coefficients in terms of either chosen set. Although we choose α m , α p and c 2 s at this stage, in the intermediate steps of the following calculations, we shall sometimes use terms from both sets, if necessary, in order to keep our expressions relatively simple. Choosing c s , c p andṪ as the independent quantities will also benefit us when considering particular situations. As discussed in Appendix B, the meaning of these quantities are more straightforward if one considers a neutral fluid whereṪ = c 2 s , or a fluid in the presence of conformal equation of state whereṪ = c 2 s = 1/3 andċ s vanishes.
Another commonly known quantity we shall find useful in what follows is heat conductivity: in terms of which the diffusion coefficient is simply We introduce a collective notation for the fluctuating modes, where normalization of the modes is chosen to make resulting matrix equations simpler and more symmetric. We can then write the above equations for the linearized fluctuations (Eq. (2.24)) in a compact matrix form, where L, D, and K are 6 × 6 matrix operators. The operators L and D are the ideal and dissipative terms, respectively, K contains the corrections due to the first-order gradients of background flow, and six-vector ξ A denotes the random noise. Explicitly (2.32) Equation (2.31) for linearized fluctuations provides the foundation for the fluctuation evolution equations for the two-point correlation functions, derived in the next section.

III. FLUCTUATION KINETIC EQUATIONS
The physical effects of fluctuations on hydrodynamic flow manifest themselves through two-point functions. This is because, by definition, the first order fluctuations average to zero (i.e. φ A (x) = 0 via Eq. (2.16)) and the leading order corrections to T µν and J µ come from the second order terms in the fluctuation expansion (i.e. the two-point functions) whose time evolution equation we derive in this section. How these two-point functions modify the hydrodynamic flow, in other words the feedback of fluctuations on background flow, will be discussed in Section IV.
Our strategy is to use equations of motion for linearized fluctuations, Eq. (2.31), to derive an evolution equation for the "equal-time" two-point correlation function of fluctuations, φ A (x + )φ B (x − ) , obtained by averaging over the statistical ensemble generated by the stochastic noises. Before we do so, we discuss some general features of the two-point correlator of hydrodynamic variables In a static homogeneous equilibrium state of the fluid, the correlator is translationally invariant, i.e., depends only on the separation y = x + −x − and not on the midpoint position x = (x + +x − )/2. Furthermore, because equilibrium correlation length is shorter than the coarse grained resolution of hydrodynamics, the equilibrium equal-time correlation function is essentially a delta function of the separation vector y with the magnitude determined by the the well-known functions of average thermodynamic variables and n. However, in a generic relativistic hydrodynamic flow several new observations need to be made. First of all, the concept of "equal time" is no longer obvious, as it depends on the frame of reference. The most natural choice, the rest frame of the fluid, is now locally different in different points. We will discuss how to implement it in the next subsection, following the formalism introduced in Ref. [11].
Furthermore, not only the local thermodynamic conditions, and thus equilibrium magnitude of fluctuations, slowly vary in space-time, but also the fluctuations themselves are driven out of equilibrium. Therefore, not only the fluctuation correlator depends slowly on x, but it also acquires nontrivial y dependence, beyond the equilibrium delta function. It is crucial that the scale of that y dependence is short compared to the scale of the dependence on x.
The estimate of the y-dependence scale can be made by observing that the equilibration of fluctuations of hydrodynamic variables is a diffusive process (since the variables obey conservation equations). This means that the scale of equilibration 4 * is the diffusion length during time interval characteristic of the evolution. For the reciprocal quantities such as fluctuation wave-number q * ≡ 1/ * and the frequency c s k of the sound, one obtains γq 2 * ∼ c s k and thus q * = c s k/γ k. In other words, * L ≡ 1/k. This separation of scales of y and x dependence of the correlation function, or between characteristic wave-numbers q of the fluctuations and k of the background will be used to systematically organize our calculations and results in the form of an expansion in k/q 1 as well as k mic 1. Note that, for the characteristic wave-numbers of the fluctuations and the background, the ratio k/q ∼ (k mic ) 1/2 . In other words, this expansion is controlled by a power of the same small parameter as the hydrodynamic gradient expansion itself.
With this separation of scales in mind, it is convenient to work with the Wigner transform of G AB (x, y), that is essentially the Fourier transform with respect to (spatial components of) y, which we shall label as W AB (x, q). Since q corresponds to the wave vector of fluctuating modes that contribute to G AB , it is similar in concept to the momentum of a particle in quantum mechanics. In this quantum mechanical analogy, the Wigner transform would be the (matrix valued) phase space distribution of the fluctuation modes or a density matrix in phase space, (x, q), in an effective kinetic theory of fluctuation quanta. The evolution equation of W AB (x, q), which is derived in this section, closely resembles a Boltzmann-type kinetic equation for the fluctuation degrees of freedom that are, in the case of hydrodynamics, phonons. In our previous work, Ref. [11], we derived and studied such an equation for relativistic hydrodynamics without a conserved charge. In this paper we present its generalization to the case with a conserved U(1) charge, which contains additional nontrivial features we will discuss in the subsequent sections IV and V. Some of these features, such as the existence of the slow scalar mode, play an important role in the critical dynamics near the QCD critical point which we discuss in Section V. The derivation in this section closely parallels the analysis presented in Ref. [11]. For completeness of this paper we will briefly summarize the key concepts and steps from Ref. [11] here before presenting our final result at the end.

A. Confluent correlator and confluent derivative
The concepts of "equal-time" and "spatial" y coordinates we invoke when defining G AB (x, y) and its Wigner transformation in the above discussion become nontrivial in a general background of relativistic flow. Both concepts require choosing a frame of reference. The most natural choice -the local rest frame of the fluid, characterized by the (average) fluid velocity u µ (x) -varies point to point with x. The change of the frame from point to point is responsible for changing the values of various vector components of hydrodynamic fluctuations φ A , such as δu µ , entering in the definition of G AB in Eq. (3.1). This variation is purely kinematic (Lorentz boost) and has nothing to do with the local dynamics of fluctuations that we are interested in. Our goal it to define a measure of fluctuations and a measure of its changes with space and time to be independent of such mundane kinematic effects. We achieve this by introducing the notions of "confluent correlator" and "confluent derivative" which we describe below, summarizing Ref. [11].
The key to defining these new concepts is a parallel transport or, equivalently a connection, that takes care of the change of u(x) between two points, say, x and x + ∆x. We introduce a boost Λ(∆x) which maps In principle, this boost is not unique. In our previous work, Ref. [11], we propose to use the most natural choice, that is a pure boost without a spatial rotation in the local rest frame of u(x). Note that the boost in Eq. (3.2) is defined for arbitrary ∆x. In practice, however, we only need its infinitesimal form: Next, we introduce the notion of "confluent correlator". This notion arises because a certain property of the "raw" definition of the two-point correlator G AB prevents us from cleanly separating x and y dependence and performing Wigner transform. Specifically where u A ≡ (0, 0, u µ ). These constraints follow from the orthogonality u µ (x ± )δu µ (x ± ) = 0 and relate different vector components of G AB (x, y) in a y-dependent way. We can deal with this problem by using Eq. (3.2) to boost G AB (x, y) in such a way that instead of being orthogonal to u A (x + ) and u B (x − ), it is orthogonal to u A (x) and u B (x). Thus we define the "confluent correlator" as where Λ C A (∆x) = Λ ν µ (∆x) when AC = µν, and an identity transformation otherwise. It is straightforward to check that the confluent correlator indeed satisfies i.e., the constraints are now independent of y. This allows us to meaningfully perform the Wigner transformation of this object with respect to y coordinates without affecting the constraint. Correspondingly, the confluent Wigner function, defined as a Fourier transform on the locally spatial hyper-surface u(x) · y = 0, is given by and obeys Note that, although the wave vector q is a four-vector, W AB depends only on its projection on the hyperplane defined by u(x) · q = 0. In order to eliminate the redundant component along u(x) we can impose the constraint u(x) · q = 0. Because this constraint depends on u(x), a meaningful derivative of W AB (x, q) with respect to x should be then defined with a parallel transport of q by Λ(∆x) from Eq. (3.2) to maintain the constraint. We shall also use the same transport to eliminate the purely kinematic effect of the boost on the vector components of variables φ A . This leads to the notion of "confluent derivative",∇ µ W AB (x, q), which we define as It is straightforward to see that∇ µ W AB (x, q) is equal to the Wigner transformation of the confluent derivativē ∇ µḠAB (x, y) similarly defined by Although the above confluent derivative is well defined conceptually, its practical evaluation requires us to introduce a local basis in the spatial hyper-surface of u(x) · q = 0 at each point x, the triad e µ a (x) (a = 1, 2, 3) satisfying u(x) · e a (x) = 0 and e a · e b = δ ab . The choice of the triad field e a (x) is arbitrary and different choices are related by local SO(3) rotations. Using this basis, we can write q = e a (x)q a with an internal three-vector q = {q a } ∈ R 3 , and consider W (x, q) as a function of q: (3.11) Working out (3.9) explicitly, we obtain where the connection,ω arises from the Lorentz boost acting on indices A and B, whereas the internal SO(3) connection, is due to the fact that Λ(∆x) −1 e a (x) is in general not equal to e a (x + ∆x), but could be related by an additional SO(3) rotation. 5 Note that the partial derivative ∂ µ in (3.12) is taken at fixed q, not fixed q = e a (x)q a . From here on we will use the notation W AB (x, q), with understanding that it is a function of x and q given by Eq. (3.11).

B. Kinetic equation for the Wigner function WAB
Having introduced the necessary mathematical tools, we now derive the evolution equation of W AB (x, q) by using, as the starting point, the stochastic equation of motion for linearized fluctuations, Eq. (2.31). The crux of this derivation is expressing the equation for W AB in terms of the confluent derivatives, which have a clear physical meaning as the derivatives in the co-moving frame. This leads to many nontrivial cancellations.
We start with the evolution equation for the two-point function , with x ± = x ± y/2, and choose y to be spatial in the frame u(x). The time evolution of G AB (x, y) is obtained by where the derivative operator, ∂, always acts on the first argument of the function, such as x in G(x, y), or x ± in φ(x ± ). Derivative with respect to the second argument, if there is any, will be labeled explicitly. Next, we convert the time derivatives in the right hand side of (3.15) into spatial derivatives. In order to do so we have to expand u(x) = u(x ± ) ∓ 1 2 y · ∂u(x) and use the evolution equation for the one point function, Eq. (2.31). To perform the resulting averaging in the r.h.s. of Eq. (3.15) we need to know the average of the two-point function of the noise, which can be calculated using the definition in Eq. (2.32) and Eq. (2.12): Proceeding from Eq. (3.15) along these steps we arrive at where the superscript (y) on an operator indicates that the derivatives within that operator act on y, the second argument of G AB (x, y). For example, The matrix Y, results from the y-dependence in u(x ± ) and c s (x ± ). Note that in deriving Eq. (3.18), we neglected higher order terms in y, based on the scale separation between background wave-number k and fluctuation wavenumber q: (∂u)y ∼ (∂c s )y ∼ k/q 1. Eq. (3.18) for G AB can then be used to derive the evolution equation for the Wigner function W AB , while expressing all derivatives in terms of the confluent derivatives. After some algebraic manipulations, we find the following result, as expressed in matrix form, where [A, B] = AB − BA and {A, B} = AB + BA are the usual matrix (anti) commutators, while a new notation is introduced for the anti-commutator which appears naturally in this context: The matrices that appear in Eq. (3.21) read where K is defined in Eq. (2.32) and ω µν is the fluid vorticity, is the source for random noise, and the matrix D (q) characterizes dissipation. In a static uniform background, the balance between the two gives the equation for the equilibrium value for the Wigner function: which is a fluctuation-dissipation relation. This equation is solved by where we used Eq. (2.27): α −1 m = −c p T /w. Eq. (3.26), taken together with Eq. (2.30), is in agreement with the well-known thermodynamic expectation values: V (δm) 2 = c p /n 2 , V (δp) 2 = c 2 s T w, V (δu) 2 = T /w and δmδp = δmδu = δpδu = 0, where V is the volume of the system. The matrices K , Ω λ , and H λ encode the effects of background gradients, that drive the system out of equilibrium.

C. Averaging out fast modes
Some of the components of W AB (x, q) oscillate fast with a characteristic frequency ω ∼ c s q, due to the L (q) ∼ c s q term in the matrix kinetic equation (3.21). According to our hierarchy of scales, the other terms in (3.21), are of order either k or γq 2 which are smaller than this oscillation frequency, c s q. This separation of time scales leads to a new effective description of the system where the fast components of W AB are eliminated by time averaging and only slow modes remain. The corresponding coarse-graining time scale b t satisfies The slow components of W AB that survive time averaging correspond to effective distribution functions in a Boltzmann-like kinetic theory of fluctuations. Note that this is also similar to how we diagonalize a quantum density matrix to identify the particle distribution functions starting from quantum field theory.
To identify the fast components, we express the kinetic equation in the basis where L (q) is diagonal. L (q) has six eigenvalues: corresponding to six eigenvectors ψ A where A = m, +, −, T 1 , T 2 , . We arrange the eigenvectors to form an orthogonal transformation matrix whereq = q/|q| is the unit vector along q and t (1) and t (2) are two transverse unit vectors that satisfy Note that the last eigenvector is a consequence of orthogonality constraint and is not a physical fluctuation mode. The choice of the dyad t (i) (x, q) is not unique, and is subject to SO(2) rotations that are local in both x and q spaces. This local freedom will bring about additional connections in the confluent derivatives, after we project W AB onto the slow components. We go to the basis where L (q) is diagonal by the orthogonal transformation M → ψ T M ψ, upon which W AB transforms to The spurious components W A , W B , and W vanish automatically due to the constraint, Eq. (3.8), and we are left with 5 × 5 matrix W AB . In the basis, we have which means that the modes with λ A = λ B are the fast modes. They average out on the coarse grained time scale b t and thus can be neglected. The remaining modes are not all independent. In particular, is the longitudinal mode associated with sound fluctuations. The remaining diffusive modes form a 3 × 3 matrix and obey W AB (x, q) = W BA (x, −q), i.e., only six of these modes are independent. These seven independent components, W L and W AB , (A, B = m, T 1 , T 2 ), constitute the degrees of freedom in the new effective kinetic description of fluctuations. Note that the 3 × 3 block of W AB ≡ W (A, B = m, T 1 , T 2 ) still contains off-diagonal components, which reflects the fact that the three modes of A = m, T 1 , T 2 are degenerate and can mix with each other. The kinetic equation for the surviving slow components follows straightforwardly from Eq. (3.21). The sound fluctuation mode completely decouples from other components and satisfies where the sound damping coefficient γ L is given by and γ ζ , γ η and γ p are defined by Eq. (2.22). Here, the confluent derivative of W L is defined as consistent with the fact that W L behaves as a Lorentz scalar. Defining such that its equilibrium value, N (0) L = T /c s |q|, is equal to what one would expect for the distribution function of "phonons" with the dispersion relation ω = c s |q|, Eq. (3.38) can be recast into the form that resembles a Boltzmann kinetic equation for phonons, Remarkably, the advection operator L L [N L ] is precisely equal to the Liouville operator, derived in Ref. [11], in a relativistic kinetic theory of massless particles (phonons), with an effective inverse metric g µν eff (x) = −c 2 s u µ u ν + ∆ µν that gives the local dispersion relation of sound waves ω = c s |q|. It should be emphasized that the Liouville operator emerges after ∂ ⊥ α terms vanish due to rather nontrivial cancellations. The simplicity of the collision (relaxation) term in the right hand side, emerging after cancellation of all the background gradient terms in Eq. (3.34), is equally striking.
The diffusive and transverse shear modes, contained in 3 × 3 matrix W , satisfy the matrix equation Here we introduced a covariant q-derivative that takes into account the rotation of the basis t (i) (x, q) of the transverse modes in q space: The confluent derivative in Eq. (3.39) also includes additional SO(2) connection ω ij µ ≡ t (i) ν ∂ µ t (j)ν , associated with the x-dependence of the basis vectors t (i) : Introducing the rescaled variables and also a Liouville-like operator, we can simplify Eq. (3.39) substantially:  (3.39), are the main results in this section. In the next section, these equations will be used to isolate short-distance singularities in the energy-momentum tensor and the charge current, when we consider two-point correlator contributions to these observables. This procedure can be done analytically, which allows us to absorb the short distance singularities into the renormalization of the equation of state and the first-order transport coefficients. After the renormalization procedure is carried out analytically, the resulting renormalized first-order viscous hydrodynamics with fluctuations will be well-defined and not suffer from short-distance ambiguity (i.e., cutoff dependence). These equations can be then applied to numerical studies of fluctuations in hydrodynamically evolving systems, such as heavy-ion collisions.

IV. FEEDBACK OF FLUCTUATIONS
Having studied the dynamics of the fluctuating modes described by the two-point functions, W AB (x, q), in the previous sections, we now discuss how these fluctuations affect background hydrodynamic flow. Hydrodynamics describes the evolution of the average values, or one-point functions, of hydrodynamic variables, such as , n, u, or more precisely, by our choice, m, p, u. The equations governing this evolution are obtained by averaging conservation equations (2.1). However, the evolution of the one-point functions is affected by the feedback from the higher-point functions. This is because energy momentum-tensor and charge current are nonlinear functions of the fluctuating variables,m,p andȗ, as follows from the constitutive relations, Eq. (2.4), as well as the equation state. Therefore expanding the fluctuating variables in φ A inside T µν and J µ to quadratic order, we get contributions proportional to In this section we discuss two aspects of the fluctuation feedback: (i) the renormalization of the variables, the equation of state and the transport coefficients as well as (ii) the time lagged hydrodynamic response, falling off as a power of time, known as "long-time tails", or, equivalently, non-analytic frequency dependence of the response at low frequencies.
In order to calculate the contribution of the two-point functions, we begin by expanding the energymomentum tensor and the charge current given in Eq. where the coefficients of linear terms were already defined in Eq. (2.17). The coefficients of bilinear terms are third order thermodynamic derivatives and are defined similarly (see Appendix A, Eqs. (A10)). Similarly to expressions for second-order thermodynamic derivatives in terms of three independent ones c s , c p andṪ in Eq. (2.27), the third order thermodynamic derivatives can be also expressed in terms of two independent third order derivativesċ p andċ s as 6 As a result, we obtain the following expansion for T µν and J µ : The mixed term, G mp (x) ∼ δmδp , is dropped because it is a rapidly oscillating component of G whose contribution vanishes after time averaging, as explained in Sec. III C. Furthermore, we neglect the fluctuations of the viscous part Π µν relying on the scale hierarchy γq ∼ q/T 1. The two point functions G AB (x) given by solutions of fluctuation kinetic equations are nontrivial functionals of the background gradients and contain both local and nonlocal terms which are associated with renormalization and long-time tails respectively.

A. Renormalization of variables and hydrodynamic coefficients
The locality of the noise in stochastic hydrodynamics is manifested by the delta functions in Eqs. (2.12). In the coarse-grained picture, this singularity is smeared out and the amplitude of the noise is proportional to b −3/2 where b is the size of the fluid cell. That means taking b → 0 requires infinitely large noise. The fluid cell must be larger than the microscopic correlation length, say T −1 or ξ whichever larger, for hydrodynamic description to be valid, but it is otherwise arbitrary. And because it is arbitrary, the physical results obtained from hydrodynamic equations cannot depend on the cutoff b.
Because of the infinite (delta function) noise, in our deterministic formalism, the singularities appear as infinite contributions to G AB (x), which arise as ultraviolet (UV) divergences in the integrals over the fluctuation wave-vector q in Eq. (4.1). Introducing the UV cutoff Λ = 1/b, we expect that these Λ dependent terms must be absorbed into the renormalized variables, equation of state and transport coefficients in order for the physics to be cutoff independent.
This renormalization procedure has been by now well understood in both nonrelativistic hydrodynamics [15] and relativistic hydrodynamics without conserved charge [11,13] or in some special cases, such as, e.g., conformal fluids, [14]. In this section, we complete this line of developments by performing the renormalization of hydrodynamics of arbitrary fluid with conserved charge in arbitrary backgrounds.
It must be kept in mind that, while Λ is a high wave-number cutoff from the perspective of the scale of fluctuations, q, it is still small compared to the microscopic scales, T or ξ −1 . Therefore even the most dominant UV divergent contribution to G AB /w ∝ Λ 3 T is still a small correction to the average background variables that are of order T 4 . However, in practical numerical simulations, these corrections will introduce a noticeable cutoff dependence, and the elimination of the cutoff dependence via renormalization is not only a matter of principle, but also an issue of practical importance.
Our starting point is to identify the physical, or "renormalized", fluid velocity u R and the physical local energy and charge densities ( R , n R ) which are determined by Landau's matching condition in terms of the "bare" variables u, and n. Although the fluctuating fluid velocity is properly normalized (i.e.ȗ ·ȗ = −1), the average velocity, u ≡ ȗ , is not since We define u R such that it is normalized, u 2 R = −1. Expanding u R to first order in G AB , we obtain 7 From Eqs. (4.5a) and (4.5b) we then find R = + δ R , n R = n + δ R n. (4.8) where the fluctuation corrections to local rest frame energy and charge densities are given by In terms of the R , n R and u R , we have now the following expressions for T µν and J µ : The transformation to physical variables is not yet complete in Eqs. (4.10a) and (4.10b) -the "bare" values and n still appear in, e.g., p( , n), which will need to be expressed in terms of physical R and n R . We shall do this below.
After establishing the expressions for physical energy and charge densities, our next goal is to determine the physical values of pressure and transport coefficients. Their physical values differ from their "bare values" that appear in the constitutive relations Eq. (2.4) and Eq. (2.9) due to fluctuations. The fluctuations contain local terms that are zeroth order (nonvanishing for homogeneous backgrounds) and first order in gradients. We shall denote these as G  AB (x) respectively. The former contributes to the physical value of the pressure and the latter contributes to the physical values of the transport coefficients. The remaining parts of G AB (x), denoted by G(x) AB , are higher order in gradients (in fact, as we shall see, they are nonlocal functionals of hydrodynamic variables): where the superscripts '(0)' and '(1)' denote the terms that are zeroth order and first order in gradient expansion 8 . Similarly since δ R and δ R n in Eq. (4.9) are linear combinations of G AB (x), these quantities can be also expanded: where expressions for δ AB and G AB respectively. By substituting this gradient expansion, Eqs. (4.11) and (4.12) into Eqs. (4.10a) and (4.10b) we can identify the physical values of pressure and transport coefficients by collecting terms zeroth order in gradients into physical (renormalized) pressure p R and terms first order in gradients into physical (renormalized) values of kinetic coefficients: where the zeroth-order terms in gradient expansion are given by , (4.14) 8 Note that G AB (x) still depends on x via terms such as w(x) however it does not contain any gradient terms such as ∂µu or ∂µα and it does not vanish in a homogeneous background. G (1) AB (x) terms are explicitly linear in gradients and do vanish in a homogeneous background. and the first-order terms are given by The remaining, i.e., higher-order in gradients (and nonlocal), contributions to constitutive equations are given by Let us now work out the explicit expressions for the physical pressure and transport coefficients. According to Eq. (4.1), the corresponding decomposition of Eq. (4.11) in phase space is given by (4.17) The zeroth-order contribution G Since W (0) does not depend on q, the integration over q is divergent. We regularize this integral by the wave-number cutoff q < Λ.
Though the tensor part of the two-point function, G (0)µν (x) appearing in Eq. (4.14) is cutoff-dependent, it is proportional to ∆ µν , and thus is absorbed into the definition of the physical (renormalized) pressure. Combining this contribution with the contributions from the terms containing δ R n in Eq. (4.14) we find for the renormalized, i.e., physical, pressure: which can be derived by using Eq. (2.27) and (4.3). This procedure of defining the physical pressure that combines "bare" pressure with the effects of equilibrium fluctuations is similar to the standard renormalization procedure in quantum field theory. Having performed the renormalization of hydrodynamic variables and the equation of state, in what follows, for notational simplicity, we will drop the subscript R on hydrodynamic variables R , n R and u R and thermodynamic functions such as p R .
We now turn to the first-order terms in the gradient expansion given by Eq. (4.15). Since these terms are linear in gradients, they must be combined with the "bare" transport terms into the physical transport terms. It may seem that this procedure, similar to renormalization of pressure, is guaranteed to succeed. It indeed does, but this is not trivial because not all gradient (transport) terms are allowed by second law of thermodynamics. The fact that only those that are allowed arise from fluctuation emerges after delicate cancellations and is a nontrivial test of the conceptual validity of the framework we develop.
Let us begin with calculating W (1) AB (x, q). This calculation essentially follows the same steps given in Ref. [11] but in this case there is an additional mode associated with conserved charge. We begin with inserting the decomposition Eq.  AB (x), we can use the ideal equations of motion to convert the time derivatives into spatial derivatives: In deriving these, we use thermodynamic relations given in Appendix A. Keeping only the terms that are linear in background gradients, the equations for W AB can be solved as (4.23) Note that these expressions are given in the (A, B) basis where L (q) is diagonal and we need to convert them back into the (A, B) basis:  AB (x, q) in components, AB (x) after q integration are given by Finally we substitute the above expressions for G (1) into Eq. (4.15). The resulting contributions are linear in the gradients and have the same form as "bare" viscous terms in Π µν and diffusion term in ν µ . Therefore they can be absorbed into the definitions of viscosities η, ζ and conductivity λ. After straightforward computation, we obtain the renormalized transport coefficients as A couple of comments are in order. First, all the gradients appearing in the expansion of G (1) are matched by the gradients appearing in the first-order terms in the constitutive equations, Π µν and ν µ . For Π µν , this is a simple consequence of the fact that, by construction, Π µν involves all gradients allowed by Lorentz symmetry, so nothing else could have appeared in Eqs. (4.26a) or (4.26c). However, this is less trivial in the case of the corrections to ν µ . This is because there are two linearly independent gradient terms allowed by Lorentz symmetry alone, e.g., ∂ µ α and ∂ µ p, and, naively, any their linear combination could have appeared in the expression for G (1) in Eqs. (4.26b). However, precisely ∂ µ α appears in Eqs. (4.26b), which allows us to absorb the fluctuation contribution into λ R . Any other linear combination would require additional kinetic coefficient to absorb it. However, the second law of thermodynamics only allows the gradient ∂ µ α to appear in ν µ in order to guarantee the semi-positivity of entropy production rate. The way this constraint is respected by fluctuation contributions appears to be highly nontrivial, relying on delicate cancellations that result in rather elegant thermodynamic identities given by Eq. (4.21). Of course, we can view this as one of the many nontrivial checks of the consistency of this approach and the validity of the calculations. Second, in a similarly remarkable deference to the second law of thermodynamics manifested in delicate cancellations, the correction to the bulk viscosity given in Eq. (4.28) is nonnegative. Also, as expected, but similarly achieved through nontrivial cancellations, the fluctuation corrections vanish in the conformal limit, where c 2 s = 1/3 and = 3p, when bulk viscosity must vanish.

B. Long time tails
After all constitutive equations are expressed in terms of the physical, i.e., renormalized, variables, pressure and transport coefficients, the remaining contributions, denoted by T µν are cutoff independent. This is very similar to renormalization in quantum field theory, and it works for a similar reason -the locality of the first-order hydrodynamics (similar to the locality of quantum field theory Lagrangian). On a more technical level, the gradient expansion in W AB is accompanied by the expansion in 1/q 2 . This can be traced back to the power-counting scheme in which k ∼ q 2 . The terms of order k 2 are accompanied by 1/q 4 leading to convergent integrals in G AB .
Thus, expressed in terms of physical quantities, the constitutive equations (4.13) do not contain UV divergences which could lead to cutoff dependence. Together with the conservation equations and the fluctuation evolution equations (3.38) and (3.39), they now form a closed set of cutoff-independent, deterministic equations that describe the evolution of the background flow, including the feedback of the fluctuating modes W . In principle this coupled system of equations can be solved numerically and nonlocal effects of long time tails in an arbitrary background can be studied. We leave such a numerical study for future work. Instead, for the remainder of this section we will describe important analytical properties of the long-time tails in simple backgrounds by solving the fluctuation evolution equations (3.38).
A quick look at the evolution equations, (3.39) and (3.38) leads to the following "impressionistic" expression for the nonequilibrium part of the Wigner function: where v = ±c sq or 0 depending on which mode we are considering and γ and ∂f are schematic notations for the relaxation rate coefficients and terms linear in background gradients respectively. Note that k ∼ ∂ and u · k = ω is the frequency. After subtracting the term linear in the background gradients, which is absorbed into the definitions of renormalized transport coefficients, we obtain a schematic expression for the finite part of the Wigner function: This procedure could be viewed a a subtraction scheme that regulates the phase space integral of the fluctuation modes where the local (and instantaneous) short distance term is subtracted. The integration over q leads to G ∼ k 1/2 ∂f /γ 3/2 ∼ k 3/2 which is a nonlocal functional of the gradients [15]. Notice that k 3/2 in terms of gradient expansion lies in between k (first order, viscous terms) and k 2 (second order terms) After Fourier transformation these terms lead to power-law corrections which correspond to the long-time tails.
To be more quantitative, let us consider a special case and focus on the non-analytic ω dependence, by taking spatial k to zero for simplicity. This means that we only keep the k dependence for the background gradient term ∂f that is in the numerator of Eq. (4.31) which is consistent with the order of gradient expansion that we are working with. In other words we are looking at the frequency dependence of the transport coefficients. From Eq. (4.32) we see that the frequency dependence can be expressed as (4.33) The contribution of the two-point functions to the constitutive relation for the charge current is given in Eq. (4.16b). We can calculate the relevant W (x, q) by using the substitution, Eq. (4.33), in Eq. (4.25). By plugging the resulting expression into Eq. (4.16b), we obtain 34) from which we find the frequency dependent conductivity, λ(ω), to be Here, λ denotes the renormalized value of the zero frequency conductivity. This result is consistent with the already known result for the special case of a conformal, boost invariant plasma with conserved charge given in Eq. (50b) in Ref. [14]. The frequency dependent viscosities can be computed in the same way. The fluctuation contributions to the viscous tensor is: 36) where Π µν stands for Π µν (ω = 0). After substituting the ω dependence in Eq. (4.33) in Eq. (4.36) we obtain from which find the frequency dependent viscosities, (4.38) Here, η and ζ denote the renormalized values of the zero frequency viscosities.

V. FLUCTUATIONS NEAR THE CRITICAL POINT
In the preceding sections we saw that kinetic coefficients, ζ, η and λ receive contributions from fluctuations. These contributions are dominated by the fluctuations at the cutoff scale Λ and therefore depend on the cutoff.
In this section we consider the physics of fluctuations at the critical point. The main feature of the critical point is that the equilibrium correlation length of the fluctuations, ξ, becomes infinite. To maintain the separation between the hydrodynamic scales L ∼ k −1 and microscopic scales, such as ξ, we must limit the domain of applicability of hydrodynamic description to wave-vectors k ξ −1 . However, as emphasized in Ref. [9], this does not mean that hydrodynamics applies until k ∼ ξ −1 . Instead, hydrodynamics breaks down before k reaches that limitation. Hydrodynamics breaks down when the frequency of the fastest hydrodynamic mode (the sound, with ω ∼ c s k) reaches the rate of the relaxation of the slowest nonhydrodynamic mode. Near the critical point this rate vanishes much faster than ξ −1 .
The slowest non-hydrodynamic variable at the critical point is the fluctuation of the slowest hydrodynamic mode (diffusive mode m), given by N mm . The relaxation rate depends on q and equals 2γ λ q 2 for q ξ −1 . Because the contribution of the fluctuations to pressure and kinetic coefficient is UV divergent, it is dominated by the modes near the cutoff, which in the case of the critical point is effectively Λ ∼ ξ −1 . Thus the characteristic rate of non-hydrodynamic relaxation, Γ ξ , is of order γ λ ξ −2 . Together with the fact that γ λ vanishes as a power of ξ, i.e., to a good approximation γ λ ∼ ξ −1 , 9 we find that the hydrodynamic breaks down already when the frequency reaches ω ∼ ξ −3 . For the sound modes this corresponds to k ∼ ξ −3 , much earlier than ξ −1 .
To extend hydrodynamics past k ∼ ξ −3 we need to include the slowest non-hydrodynamic mode, which is the idea behind Hydro+ [9]. In our notations this mode (or modes, labeled by index q) is N mm . In this section we intend to show that in the regime k > ξ −3 our formalism reproduces Hydro+. This is a nontrivial check because Hydro+ formalism was derived in Ref. [9] using a completely different approach by considering a generalized entropy which depends on the non-hydrodynamic variables (2PI entropy).
The formalism of Hydro+, while extending ordinary hydrodynamics beyond the scales k ∼ ξ −3 , in turn, also breaks down well before k reaches k ∼ ξ −1 . The breakdown occurs when the frequency reaches the relaxation rate Γ ξ of the next-to-slowest nonhydrodynamic mode. This mode (or modes) are the fluctuations of velocity transverse to the wave-vector. This relaxation rate is of order γ η q 2 at q ξ −1 . Again, the dominant contribution comes from modes at q ∼ ξ −1 and, since γ η to a good approximation can be treated as finite at the critical point [18], Hydro+ breaks down when frequency reaches ω ∼ Γ ξ ∼ ξ −2 , which for the sound modes corresponds to k ∼ ξ −2 . Near the critical point this scale is still much lower than ξ −1 .
In our formalism the next-to-slowest modes responsible for the breakdown of Hydro+ are N m(i) and N (i)(j) (normalized Wigner functions obeying Eqs. (3.45b) and (3.45c)). Therefore, within our formalism we can extend Hydro+ beyond its limit at k ∼ ξ −2 . In Section V B we shall describe how to do that. Prior to that, in Section V A, we shall verify that in the regime where Hydro+ is applicable, it is in agreement with our more general formalism.

A. Connection to Hydro+
The main ingredient of Hydro+ is the entropy density s (+) of the system in partial equilibrium state where a non-hydrodynamically slow variable ϕ, or more generally, a set of variables ϕ q indexed by a discrete or continuous index q is not equal to the equilibrium value ϕ (0) q ( , n) for given and n. For brevity of notations we shall denote such a set of variables by a bold letter, similar to a vector with components φ q : The equations of motion for ϕ describe relaxation to equilibrium (maximum of s (+) ) accompanied, in general, by dilution due to expansion: Second law of thermodynamics requires (F ϕ ) q = q γ qq π q with semi-positive-definite γ where π q is the thermodynamic "force" defined, as usual, via where π · ϕ = q π q ϕ q . The coefficient A ϕ in Eq. (5.2) describes the response of the variable ϕ to the expansion or compression of the fluid (since θ = ∂ · u is the expansion rate). 10 The hydrodynamic variables and u obey, as usual, equations of the energy-momentum conservation. The equation of state enters into constitutive equations via pressure p (+) which, as a function of , n, and φ, is given by the Legendre transform of s (+) : This relationship between pressure and entropy is dictated by the second law of thermodynamics [9].
Near equilibrium, the deviation of the entropy s (+) ( , n, ϕ) from the equilibrium value s( , n) is quadratic in π, since entropy is maximized in equilibrium. The deviation of pressure p (+) from equilibrium p is linear in π, p (+) = p + p π · π + O(π 2 ) . (5.6) The coefficient p π can be expressed (see Appendix B in Ref. [9]) using Eqs. (5.3) and (5.5), in terms of the equilibrium value of φ at given and n, which we denote by ϕ (0) ( , n), as We wish to show that the constitutive equations in Hydro+ with generalized pressure p (+) are in agreement with the equations we derived by expanding to quadratic order in fluctuations, such as Eq. (4.13).
Application of the Hydro+ approach near the critical point consists of considering the two-point correlation function of the slowest mode (m ≡ s/n): ϕ ∼ δmδm . Essentially, using our notations Due to the reparametrization invariance of Hydro+ (see Appendix C in Ref. [9]), either choice, N mm or W mm , different by a normalization factor in Eq.(3.43), will lead to the same result. The choice of N mm is convenient because in this case the compression coefficient vanishes: A ϕ = 0 (see Eq. (3.45a)).
In order to find non-equilibrium correction to Hydro+ pressure in Eq. (5.6) we need to use the expression for the non-equilibrium contribution to entropy derived in Ref. [9] to determine π: where The equilibrium value N where we used N We should compare this to the nonequilibrium contribution to pressure from N mm (which is dominant near critical point due to being proportional to c p ) in Section IV: 13) which is similar to equilibrium contribution (renormalization of static pressure) found in Eq. (4.20) with index '(0)' replaced by '(neq)'. One can see that Hydro+ reproduces these nonequilibrium contributions exactly. We emphasize that this is a very nontrivial cross-check, involving an elaborate thermodynamic identity for third derivatives of entropy in Eq. (4.21). This is not unexpected since Hydro+ formalism emerged in Ref. [9] via very different route, starting from the derivation of the nonequilibrium entropy functional s (+) in Eq. (5.9).

B. Hydro++
Since, as we already discussed above, the fluctuation contributions are dominated by the modes near the cutoff Λ, and for critical fluctuations the role of this cutoff is played by ξ −1 , the contributions responsible for the breakdown of ordinary hydrodynamics and of Hydro+ are dominated by fluctuations at scale q ∼ ξ −1 . These modes themselves cannot be described by ordinary hydrodynamics. The dynamics of these modes is essentially nonlinear and nonlocal (often referred to as mode-coupling phenomenon). However, this dynamics is universal in the sense of universality of dynamical critical phenomena and is described by model H in the classification of Ref. [18]. We shall, therefore, use the known results from this universality class to describe the dynamics of these fluctuation modes.
Near the critical point, where the correlation length ξ greatly exceeds all other microscopic scales, the description simplifies due to (static and dynamic) scaling. That means the relaxation rates, even though no longer polynomial in q, as in the hydrodynamic regime where gradient expansion applies, depend on the q and ξ via functions of only the dimensionless combination qξ (times a power of ξ). Furthermore, these functions (and the powers of ξ) are universal, i.e., independent of the microscopic composition or properties of the system close to the critical point in a given universality class. The universality class relevant for our discussion is that of model H, defined in Ref. [18] as dynamic universality class of liquid-gas phase transitions.
As we already said, the fluctuation kinetic equations, such as (3.45), do not apply in the regime qξ ∼ 1 as they are. However, a modification of these equations, to match the known results from model H is possible and shall be described below. We must emphasize, that unlike the formalism derived in the preceding sections, which was exact to a certain order in a systematic expansion, here our out goal is to provide the formalism which reproduces the physics of critical point fluctuations correctly, but not necessarily exactly. For once, the exact description would at a minimum require exact solution to model H, which is not available. Our approximation is essentially equivalent to a one-loop approximation introduced by Kawasaki in Ref. [19], which is known to be in good quantitative agreement with experimental data [18]. Similarly to Hydro+ formalism, the purpose of the new extended formalism, which we shall refer to as Hydro++ in this paper, is to provide a practical way of simulating the dynamics near the critical point, e.g., in heavy-ion collisions.
There are two main modifications required. First of all, we need to modify equation for N mm to make sure that the equilibrium correlation function has finite correlation length ξ, i.e., N (0) mm (x, q) must depend on momentum q. We shall express this as where we defined function c p (q) in such a way that c p (0) = c p is the usual thermodynamic quantity (heat capacity at constant pressure). In this work we shall adopt a simple approximation for the momentum dependence: This is known as Ornstein-Zernike form and is consistent with other approximations we are making. 11 A more sophisticated form and a better approximation to the exact correlation function (which is not known exactly as of this writing 12 ) can be used if necessary, see Ref. [9]. The second essential modification is required to correctly describe relaxation rate of the slowest nonhydrodynamic mode, N mm . The critical contribution, ∼ ξ −1 dominates near the critical point. It is given in terms of the Kawasaki function Keeping also noncritical contribution, we can write for the q-dependent rate Note that at small q, i.e., qξ 1, the rate is given by twice the diffusion rate γ λ q 2 , where γ λ = κ/c p with κ being the zero-frequency heat conductivity: It contains a noncritical contribution κ 0 , but is dominated, near the critical point, by the critical contribution due to the fluctuations. The latter increases with ξ as κ ∼ ξ (in Kawasaki approximation). With these two modifications, the equations for Hydro++ we propose read: where again, α p = (1−Ṫ /c 2 s )/T n. The function γ λ (q) is defined in Eq. (5.17). The presence of function c p (q), defined in Eq. (5.15), in Eq. (5.19b) ensures important property of N m(i) in equilibrium -proportionality to ∂α, which follows from the second law of thermodynamics as we already discussed in connection with Eq. (3.45b). 13 Other terms may also contain "formfactors", i.e., functions of qξ, which could be determined from a more detailed calculation of 3-point functions in model-H. We leave such and similar refinements to future work. It is likely that given the general degree of applicability of hydrodynamics in heavy-ion collisions these will be beyond the experimentally relevant precision.
Eq. (5.19a) describes relaxation of the slowest nonhydrodynamic mode, N mm , to equilibrium given by Eq. (5.14). It would be identical to the corresponding Hydro+ equation in Ref. [9], but for the last term describing the coupling to next-to-slowest mode, N m(i) . Because c q ∼ ξ 2 diverges at the critical point, this term is indeed much smaller than the first term sufficiently close to the critical point. However, if we want to interpolate Hydro+ description close to the critical point with dynamics of fluctuations away from the critical point this term has to be kept.

C. Conductivity and its frequency dependence in Hydro++
Let us discuss physics described by Eqs. (5.19) which is pertinent to the breakdown of Hydro+ and its crossover to Hydro++.
We can use Eqs. (5.19) to determine the critical contribution λ ξ to the conductivity λ and verify it diverges as ξ → ∞. Following the procedure of renormalization described in Section IV we now find that W (1) mµ , i.e., the part of W mµ linear in gradients, is given by Eq. (4.25) with a simple substitution c p → c p (q). This, in turn, makes the integral of W mµ (x), finite. The cutoff is now essentially given by 1/ξ, instead of Λ. This means that instead of Eq. (4.26) we find, using c p (q) in Eq. (5.15), Λ → π/(2ξ). Substituting this result into equation (4.29) for the renormalized conductivity we find a contribution to renormalized conductivity which diverges with ξ: -a well-known result [18]. We used the fact that c p ∼ ξ 2 . In particular, since γ λ ∼ λ/c p ∼ ξ −1 we neglected γ λ compared to γ η ∼ ξ 0 . Denoting the noncritical contribution to conductivity by λ 0 we can write the total We can further quantify this description by considering the dependence of W (1) mµ on frequency following the same procedure as in Section IV B. Combining the substitution in Eq. (4.33) with the substitution (5.15) in Eqs. (4.16b) and (4.25) we find for frequency-dependent leading critical contribution to conductivity: We do not need to regularize and subtract a divergence, as we did in Section IV B, because the divergence is tamed by the fall-off of c p (q) at large q. At small ω Γ ξ Eq. (5.24) reproduces the power-law non-analytic dependence characteristic of the longtime tails in Eq. (4.35): λ ξ (ω) − λ ξ (0) ∼ λ ξ (0)ω 1/2 . Not surprisingly, since, compared to Section IV B, we only changed the nature of the cutoff Λ. At large ω we find λ ξ (ω) ∼ ω −1/2 with no ξ dependence as expected from scaling behavior characterizing this regime. 16 The dependence of λ ξ on ω described by Eqs. (5.24) and (5.25) corresponds to the physics we anticipated -the large critical contribution "switches off" when ω Γ ξ .
It may also be helpful to note that while real part of λ ξ (ω) corresponds to (frequency dependent) conductivity, its imaginary part (divided by ω) is the electric permittivity.
One can also understand frequency dependence as a time-delayed medium response to gradient of density, i.e., ∂α. The diffusive current induced by the gradient is given by The delay is given by the Fourier transform of F λ (y): As a function of t − t it has a characteristic width given by 1/Γ ξ ∼ ξ 2 and becomes δ-function in the limit ξ → 0 corresponding to instantaneous response. At large t − t it falls off as (t − t ) −3/2 typical of the long-time hydrodynamic tails. The discussion of the frequency dependence of conductivity here carries many similarities to the discussion of the bulk viscosity in Ref. [9]. For completeness, let us present the calculation of the leading critical contribution to the bulk viscosity in Hydro++, which, of course, gives the same result as Hydro+. Near the critical point the leading contribution of fluctuations to the bulk viscosity comes from G (1) mm in Eq. (5.13). In Hydro++ the corresponding W (1) mm is given in Eq. (4.23), with the substitution of c p with c p (q) as in Eq. (5.15), as well as γ λ with γ λ (q) according to Eq. (5.17). As a result we obtain for the leading critical contribution to bulk viscosity: where we usedċ p = 2ξ (according to scaling c p ∼ ξ 2 ) and Γ ξ = T /(3πηξ 3 ) (according to Eqs. (5.23) and (5.20)). We introduced This is a known result in Kawasaki approximation [9,19]. 17 At ω = 0 Eq. (5.28) gives ζ ξ (0) ∼ ξ 3 (according to the scaling ofξ ∼ ξ 3/2 ). This large critical contribution is "switched off" via function F ζ when ω > Γ ξ . 18 The resulting behavior is illustrated in Fig. 1 together with the behavior of λ(ω). 16 As before (see footnote 9), the exact value of the scaling exponent in λ ξ (ω) ∼ ω −1/2 differs slightly from the rational value −1/2. The exact value in model H following from dynamic scaling −x λ /(2 − xη) = −(4 − η − z)/(d + 2 − z) is approximately −1/2 in the Kawasaki approximation we are using, which corresponds to z ≈ 3 and ηx ≈ 0. 17 As we already discussed, Kawasaki approximation only gives a good approximation to the correct scaling behavior. To match the exact scaling behavior on would need a more elaborate choice of the substitution in Eq. (5.15), see e.g., Refs. [9,19,20].

VI. DISCUSSION, CONCLUSIONS AND OUTLOOK
In this paper we continued the development of the deterministic approach to fluctuation hydrodynamics for an arbitrary relativistic flow. We extended the method developed in Ref. [11] to the case of a fluid carrying conserved charge. In QCD the relevant charge is the baryon number. Our ultimate goal is practical -a formalism which would allow to simulate heavy-ion collisions with dynamical effects of fluctuations, especially relevant for the QCD critical point search.
We would like to emphasize that, despite its practical aim, this formalism is based on a systematic and controllable expansion, similar to the effective field theory formalism in quantum field theory. The expansion parameter in hydrodynamics is the ratio of the wave-number k = 1/L associated with background flow and density gradients to a microscopic scale which sets the scale of hydrodynamic coefficients and which we denote 1/ mic . This allows us to view hydrodynamics as an effective theory.
Instead of directly solving stochastic hydrodynamic equations, we convert them into a hierarchy of equations for equal-time correlation functions, which we truncate at two-point correlators. This truncation is controlled by the same expansion parameter as the gradient expansion in hydrodynamics. One can see how the relevant power counting emerges by considering the effects of fluctuations on the constitutive equations for stress tensor (or conserved current). In stochastic hydrodynamics the noise is local, i.e., it is only correlated inside a hydrodynamic cell, as reflected in the δ-function value of the two-point noise correlator in Eq. (2.12). This locality is the source of short-distance singularities, similar to ultraviolet singularities in quantum field theories. Hydrodynamics is regulated by finiteness of the cell size, which we denote by b mic , equivalent to wave-number cutoff Λ = 1/b. As a function of this regulator, the square variance of the noise in each cell is proportional to Λ 3 -the regulated value of the δ-function. This is, of course, the source of the cutoff dependent contribution to renormalized pressure in Eq. (4.20) and, as such, is not of physical relevance.
The physically consequential contribution comes from the fluctuations whose relaxation time is comparable to the evolution time of the background. Correspondingly, this scale, characterized by wave-number q * , can be estimated by the condition γq 2 * ∼ c s k. The effect of these fluctuations is the delayed or nonlocal response to perturbations of the background (such as long-time tails) and cannot be simply absorbed by renormalization of the local hydrodynamic parameters such as pressure or transport coefficients. Since q * Λ, the noise on these longer distance scales, * = 1/q * , averages out and the magnitude of the fluctuations is effectively reduced by a factor (b 3 / 3 * ) 1/2 = (q * /Λ) 3/2 -the inverse of the square root of the number of uncorrelated cells in a region of linear size * -the familiar random walk factor. Therefore the physically relevant magnitude of the fluctuations, obtained by averaging over scales * is given by Λ 3/2 × (q * /Λ) 3/2 ∼ q 3/2 * ∼ k 3/4 . It is cutoff independent, of course. Therefore, the two-point correlator of these fluctuations contributes at order k 3/2 , suppressed compared to first-order gradients, but more important than second order gradients. Similarly, the contribution of n-point functions, due to higher order nonlinearities in the constitutive equations, would come at order k 3n/4 One can see that the hierarchy of higher-point contributions is controlled by a power of k, or more precisely, a power of dimensionless parameter k mic = mic /L 1. The equations we derive form a closed set of deterministic equations which can be solved numerically. The one-point functions (averaged values of hydrodynamic variables) obey conservation equations (4.30). The constitutive equations (4.13) contain contributions T µν and J µ which are given in terms of the subtracted two-point functions G in Eqs. (4.16). The unsubtracted two-point functions G are evaluated at coinciding points and therefore contain short-range singularities. When unsubtracted G are expressed in terms of the wave-number integrals of the Wigner functions Eq. (4.1), these singularities appear as ultra-violet divergences which need to be subtracted. The unsubtracted Wigner functions are obtained by solving equations (3.38) and (3.45), rescaling according to Eqs. (3.37) and (3.43) and substituting into the matrix in Eq. (4.24). The subtraction of terms of zero and first order in gradients, W (0) and W (1) , given by Eqs. (4.18) and (4.25) respectively, can be done analytically, and either before or after solving equations (3.38) and (3.45), depending on numerical efficiency. The resulting solutions to one-point and two-point equations will describe evolution of the average hydrodynamic variables, their fluctuations, as well as the feedback of the fluctuations on the evolution of average quantities.
With the equations we derived we can now also describe the essential features of the hydrodynamic evolution near the QCD critical point. The critical phenomena are originating from the divergence of the correlation length ξ. The phenomenon of the most consequence for hydrodynamics is the critical slowing down. Since it is caused by the fluctuations of the slowest diffusive mode out of equilibrium, our formalism is ideally suited to accommodate and describe this phenomenon. The formalism of Hydro+ introduced earlier in Ref. [9] is based on the same observation and adds the two-point correlation function of the diffusive mode to hydrodynamics to describe critical slowing down. The approach in the present paper is very different from the derivation in Ref. [9], therefore, the exact agreement between the results is a nontrivial check on the validity of both derivations.
Since, our present approach is more general, we can now connect Hydro+ description of critical fluctuations to description of ordinary fluctuations away from the critical point. Because the validity of Hydro+ is limited by the relaxation rate of the next-to-slowest mode, and this mode, absent in Hydro+, is now a part of our description, we are able to extend the validity of hydrodynamic description closer to the critical point than Hydro+. We propose a set of equations, which we call Hydro++ which could accomplish this. It should be kept in mind that, unlike the systematic approach taken in the rest of the paper, the Hydro++ equations (5.19) are an attempt to interpolate between the description of fluctuations outside of the critical regime and the known properties of the fluctuations in the critical, scaling regime described by model H (in the standard classification of Ref. [18]). While the hydrodynamic description still works for the background gradients for which k mic ∼ kξ 1, it breaks down for critical fluctuations, for which qξ ∼ 1. This means that the coefficients become non-polynomial in q and that the theory becomes fully nonlinear and the truncation to two-point functions is no longer, strictly speaking, controllable. However, it is known from the studies of model H that the results obtained in one-loop (Kawasaki) approximation are in good quantitative agreement with experiment [18]. Therefore we propose a set of equations (5.19) which incorporate the model H physics at the corresponding level of approximation. This approach is similar to the one taken in the derivation of Hydro+ and extends the region of applicability closer to the critical point. More precisely, while Hydro+ breaks down at k ∼ ξ −2 , the validity of Hydro++ extends to k ∼ ξ −1 . The physical phenomenon which leads to breakdown of Hydro+ is the frequency dependence of (i.e., time-lag of) conductivity, which is described by the next-to-slowest mode in Hydro++.
Once the fluctuation hydrodynamics in the deterministic approach is implemented in a fully functional hydrodynamic code 19 , the extension to full Hydro++ approach should be straightforward and will allow eventual comparison with heavy-ion collision experiments not only near, but also away from the critical point. However, additional developments are needed to make this comparison more impactful. First of all, it should be straightforward to generalize this approach to multiple conserved charges. In the case of QCD, of course, fluctuations of isospin are a primary candidate. We have not included these fluctuations in our description because they are not exhibiting singularities near the critical point, unlike the baryon number fluctuations, which lead to signatures of the QCD critical point [21]. Furthermore, the approach must be extended to description of non-Gaussian fluctuations, which are related to most sensitive signatures of the critical point [22,23]. This means going beyond two-point correlators considered in this paper. It would also be interesting and important for comparison with experiment to consider the extension of this approach to the fluctuations near the first-order phase transition, which is, of course, an inseparable part of the physics near a critical point. We defer these and other pertinent developments to future work.
A charged fluid studied in Ref. [14] is (i) conformal and (ii) undergoes a boost-invariant (Bjorken) expansion. Thermodynamic functions of a conformal fluid satisfyṪ = c 2 s = 1/3, m = α p =ċ s = 0,ċ p = 1 as summarized in Table I. The boost-invariant flow implies that a µ = ω µν = 0 and spatial gradients of background scalar fields vanish (e.g., ∂ ⊥µ α = 0). Under these conditions our results are significantly simplified. Since in a boost-invariant Bjorken flow, the charge does not diffuse due to the absence of background gradients forbidden by boost-invariance, in order to generate the dissipative (ohmic) charge current, one needs to apply an external electric field to the system. Adding such a source term is indispensable for obtaining such important results as renormalized or frequency-dependent conductivity in Ref. [14]. We find that except for a few minor typos, our Eqs. (4.27) and (4.29) for renormalized transport coefficients, as well as Eq. (4.35) and (4.38) for frequency-dependent transport coefficients, reduce to Eq. (51) and (50) in Ref. [14] respectively. Notice that ζ = 0 in conformal fluid.
(B2) Then, our equations read The equations for N B (1) (1) and N B (2)(2) match those in Ref. [14]. The remaining equations, although very similar, do not match completely. We believe our results are correct but do not have a definitive explanation for these disagreements.
Unlike the chargeless fluid in Ref. [11], the charged fluid in the present paper can be taken to nonrelativistic limit, where it can be compared with Ref. [15,24]. The most glaring omission in Ref. [15,24] are the G m(i) components of the correlators. It appears they were omitted based on the observation that their equilibrium values vanish. They do, but they are not zero out of equilibrium and are essential, for example, for describing the dominant critical contribution to conductivity as discussed in Section V C.
To compare the equations for remaining components of W AB , we need to rescale our variables as Omitting W m(i) terms as in Ref. [15,24] we find