Non-Gaussian fluctuation dynamics in relativistic fluids

We consider non-equilibrium evolution of non-Gaussian fluctuations within relativistic hydrodynamics relevant for the QCD critical point search in heavy-ion collision experiments. We rely on the hierarchy of relaxation time scales, which emerges in the hydrodynamic regime near the critical point, to focus on the slowest mode such as the fluctuations of specific entropy, whose equilibrium magnitude, non-Gaussianity and typical relaxation time are increasing as the critical point is approached. We derive evolution equations for the non-Gaussian correlators of this diffusive mode in an arbitrary relativistic hydrodynamic flow. We compare with the simpler case of the stochastic diffusion on a static homogeneous background and identify terms which are specific to the case of the full hydrodynamics with pressure fluctuations and flow.


I. INTRODUCTION
The physics of thermal fluctuations [1] in hydrodynamics [2] has received renewed interest recently in the context of relativistic heavy-ion collision experiments.Hydrodynamics proved to be remarkably successful in describing the data from such collisions [3,4].One of the major goals of heavy-ion collision experiments is the discovery of the QCD critical point -the end point of the conjectured first-order transition separating hadron gas and the quark-gluon plasma phases [5].This point is characterized by a certain singular behavior of fluctuations -a universal feature of the critical points [6][7][8].Therefore, understanding of fluctuations and, in particular, of their non-equilibrium evolution in the environment of the hydrodynamically expanding QCD fireball is important to enable the experimental search for the QCD critical point [9,10].
The parameter controlling the importance of fluctuations at wavelengths of order ℓ is the ratio of the "correlation volume" of the size of the correlation length ξ to the "homogeneity volume" of order ℓ 3 .We shall generically refer to this parameter as ε.The central limit theorem suppresses the magnitude of fluctuations as well as their non-Gaussianity by a power of ε.In a typical condensed matter experiment ε is extremely small and thus fluctuations often play a negligible role.In heavy-ion collisions, where the system size (O(10 fm)) is large compared to the typical correlation length ξ (a fraction of fm), the parameter ε is small but not negligible.Fluctuations are observable in heavy-ion collisions and play important role in understanding the thermodynamic properties of the QCD fireball.
Furthermore, near the critical point, as the correlation length becomes longer, the fluctuations grow and play even more important role.Non-Gaussianity of fluctuations, also controlled by ε, is increasing even faster at the critical point than their magnitude [11], and is expected to show a specific nonmonotonic dependence on the collision energy [12], which in heavy-ion collisions controls the thermodynamic conditions at freezeout.Therefore, non-Gaussian fluctuation measures have emerged as observables of primary experimental interest in the search for the QCD critical point [13].
Another effect of the critical point is to slow down the relaxation towards local thermodynamic equilibrium.The critical slowing down makes it increasingly important to consider non-equilibrium evolution of fluctuations near the critical point, as earlier estimates and model calculations demonstrate [7,14,15].To make quantitatively reliable predictions, however, one needs to consider the non-equilibrium behavior of fluctuations in the "first principle" hydrodynamic approach.
There has been considerable recent effort and progress towards understanding non-equilibrium evolution of fluctuations in the context of the hydrodynamics of the heavy-ion collisions with the aim of mapping the QCD phase diagram (see Ref. [10] and references therein).The focus of this paper is on the approach where fluctuations are described in terms of the correlation functions obeying deterministic evolution equations, i.e., the so-called deterministic (also known as hydrokinetic) approach.This approach was introduced and developed recently in the context of relativistic heavy-ion collisions [16][17][18], more generally in [19][20][21][22][23], and much earlier, in the non-relativistic condensed matter context, in [24,25].
Within the deterministic approach most of the work has been done on Gaussian fluctuation measures: two-point correlators, or their Wigner transforms.In this approach, the first step towards the study of the non-Gaussian hydrodynamic fluctuations out of equilibrium was made in Ref. [22], 1 where the generalized n-point Wigner transform was introduced and the equations for the corresponding n-point Wigner functions were derived.Ref. [22] considered the simplest hydrodynamic system: nonlinear charge diffusion at constant temperature.Such a system is characterized by a single hydrodynamic field: conserved charge density.Full hydrodynamics necessary to describe heavy-ion collisions involves (at least) five hydrodynamic fields: conserved densities of energy, baryon charge and three-momentum of the fluid.The full theory of hydrodynamic fluctuations is an ambitious goal and in this paper we shall present another step towards it by considering fluctuations of both energy and baryon density in the regime relevant to the critical point search, where the correlation length becomes large compared to typical microscopic scales (while still being smaller than ℓ) and fluctuations of a certain critical mode dominate.
In this regime an important hierarchy of scales emerges.The slowest mode is the diffusive (nonpropagating) mode m = s/n, i.e., entropy to baryon number ratio, which we shall refer to as specific entropy.This mode does not mix with propagating sound oscillations, which makes it purely diffusive (unlike fluctuations of energy and baryon densities on their own).Furthermore, unlike the case of other diffusive modes, such as transverse momentum densities, the diffusion constant for m vanishes at the critical point making m the slowest mode in this regime.Such a hierarchy of scales was exploited in the construction of Hydro+ theory in Ref. [29] as well as in the Hydro++ theory in Ref. [20], where the next-to-slowest momentum diffusion modes were also included.In both cases only two-point correlators were considered.The goal of this work is to extend the Hydro+ formalism to the non-Gaussian fluctuation measures, i.e., n-point correlation functions.
As in Ref. [29] we shall focus on fluctuations of m, which are important for two related reasons.First, as explained above, this mode and its fluctuations are slowest to relax, and therefore are most in need of nonequilibrium description.Second, the magnitude and non-Gaussianity measures of the fluctuations of this mode diverge with the correlation length most strongly (with the largest critical exponent), thus providing the most sensitive observable signatures of the critical point.
Furthermore, in Ref. [22] we considered fluctuations in a stationary fluid without any flow.Here we shall treat the most general relativistic flow in a fully Lorentz covariant formalism necessary for applications to heavy-ion collisions.We review the formalism which was introduced in Ref. [19,20] for two-point correlators with a specific focus on generalizing it to non-Gaussian fluctuations.We use the n-point wave number dependent correlation functions which we introduced in Ref. [22] by generalizing the well-known Wigner transform and show how the confluent formalism of Refs.[19,20] naturally extends to these objects in Section II D.
While the simple diffusion problem in Ref. [22] contains only one fluctuating field variable (charge density n), we find that-even in the regime where the fluctuations of m are the slowest and the fluctuations of the faster hydrodynamic modes, such as pressure p, can be considered as equilibrated-we cannot neglect these fluctuations, as they modify equations for the evolution of fluctuations of m via nonlinearities in the equation of state, or mode coupling.We discuss this nontrivial part of the derivation in Section III B.
We should point out again that we set up the formalism in the most general form necessary to tackle the full system of hydrodynamic equations with fluctuations of all hydrodynamic variables.Our focus on the slowest variable m allows us to make the logical first step towards considering the full system.This focus simplifies the calculations and allows us to hone the tools needed to tackle this ambitious goal.
We shall use the approach of Refs.[19,20], i.e., expand the stochastic hydrodynamic equations in powers of fluctuation magnitude.Unlike Refs.[19,20] we shall not consider feedback of fluctuations and assume, based on the analysis of these contributions in Refs.[19,20], that the major effect of these (UV sensitive) contributions have been absorbed into renormalization of hydrodynamic variables, equation of state and transport coefficients which, therefore, take physical (cutoff independent) values.In the diagrammatic representation we introduced in Ref. [22] such feedback contributions are represented by loop diagrams.As in Ref. [22], we shall consider only "tree-level" terms in the evolution equations for non-Gaussian correlators.Again, this is not a limitation of the approach, but a natural simplifying first step in its development.
In Sec.II we introduce the general formalism for the dynamical evolution equations for the n-point correlation functions in the presence of background hydrodynamic flow.In Sec.III we focus on the specific entropy fluctuations for reasons explained above and, by using the results developed in Sec.II, we derive the evolution equations for the two-, three-and four-point functions of the specific entropy fluctuations.Thus we arrive at the main results of this paper which can be found in Sec.III B. We summarize our findings and discuss the outlook for future developments in Sec.IV.

II. GENERAL DETERMINISTIC FORMALISM FOR FIELD FLUCTUATIONS
In this section we extend the formalism for the dynamical evolution equation of relativistic hydrodynamic fluctuations we developed in Refs.[19,20] to n-point correlation functions.In particular, we show that the confluent formalism, which allows us to describe fluctuations in the local rest frame of the fluid in a natural and fully Lorentz covariant way, generalizes naturally to n-point functions.We also use the multipoint Wigner transform, which we introduced in Ref. [22], to express the n-point correlators in terms of one spatial coordinate of the midpoint and n − 1 independent wave vectors.The results of this section are general, i.e., not limited to a particular fluctuating variable of variables.

A. Evolution of correlation functions
Let us first derive the evolution equation for the n-point correlation functions of a set of generic stochastic field variables ψi where the subscript i labels different local hydrodynamic fields (such as entropy per baryon m ≡ s/n, pressure p, and fluid velocity u as in Refs.[19,20]).Each of those variables satisfies its own Langevin-type equation that can be generically written in a covariant form, where the symbol " ˘" denotes a stochastic quantity.Functions or functionals, such as ȗ ≡ u( ψi ) or F ≡ F ( ψ) inherit the stochastic symbol " ˘" from their arguments.The four-velocity obeys ȗ2 = −1, and should be understood as a four-vector function of ψ i 's (or can also be chosen to be among the variables ψ i as in Refs.[19,20]).Fi is the drift "force" and ξ i is the noise (random "force") for the variable ψ i , expressed in terms of the canonically normalized local Gaussian noise 2 : We include multiplicative noise, since H = H( ψ), and define the product Hη in terms of the Itō calculus. 3 The Onsager matrix (operator Denoting fluctuation of a given stochastic quantity X as and introducing the fluctuation of the variables ψ i , we can expand the stochastic equations (2.1) in powers of ϕ i , using, e.g., 2 As discussed in Ref. [22], in hydrodynamics, the contribution of the non-Gaussianity of the noise will not appear in the leading order in the hydrodynamic (i.e., gradient) expansion. 3In practice that means that H( ψ) and η are considered uncorrelated.Under time discretization this corresponds to evaluating H(ψ(t)) and η(t) at the same time point t.Different discretization prescriptions, such as Stratonovich, where H = H((ψ(t) + ψ(t + ∆t))/2), can be used to describe the same physics with a given equilibrium distribution of fluctuating variables, as long as the drift term is chosen accordingly.Below we shall verify a posteriori that our equations reproduce the correct equilibrium values of the fluctuations given by thermodynamics.This provides a check of the consistency of the implementation of the Ito prescription in our approach.
where F ≡ F (ψ), etc. , and obtain the evolution equation for fluctuation field ϕ i : acts on ϕ j1 only, and the Einstein summation rule over repeated indices is implied throughout this paper.The subscript 1 . . .n denotes the "averaging over permutations", i.e., the sum over all n! permutations of all composite index-position labels (j 1 , x 1 ), . . ., (j n , x n ) divided by n!, . . .
where we used the notation of Ref. [22], P 1...n , to denote the sum over permutations.Furthermore, since Eqs.(2.1) that we consider are differential equations, each expansion coefficient such as L i, j1j2...jn must be treated as a multilinear differential operator.These operators are linear in each of their arguments, which are identified by the matching repeated indices.It is helpful, in this regard, to think of the arguments x i of the field ϕ i (x i ) as a part of the composite label (i, x i ) denoting position as well as the name or index of the variable.
Having this at hand it is straightforward to write down the evolution equation for the "raw" n-point correlation functions by taking the time derivative in the rest frame of the fluid at midpoint x (2.12) 4 Eq.(2.11) assumes that the velocity u varies slowly on the characteristic scale of the correlations we study, as is the case in hydrodynamics.This property of hydrodynamics underlines the approach known as "hydrokinetics" [16][17][18][19][20]. Eqs.(2.10)-(2.12),for n = 2 and truncated at m = 1, upon Wigner transform reproduce the equations for Gaussian fluctuation correlators known as "hydrokinetic equations".Our equations are more general and include non-Gaussian fluctuations and multiplicative noise.
By the first equality in Eq. (2.10) we defined the differential operator, which can be understood as the derivative with respect to the midpoint x (see also Eq. (2.33)).The subscripts 1 . . .n in Eq. (2.10), and 2 . . .n in Eq. (2.11), denote the "averaging over permutations" defined in Eq. (2.8).The factors n = n!/(n − 1)! and n(n − 1) = n!/(n − 2)! in Eq. (2.10) conveniently match the number of the nonidentical terms among the n! permutations.The identical terms arise due to the invariance of the correlators in Eq. (2.9) under the permutation P 1...n .The function Q(x 1 , x 2 ) is the integral kernel of the Onsager operator, i.e., Q(x 1 , x 2 ) ≡ Qδ (3) (x 1 − x 2 ), where the delta-function support is a line parallel to u(x), i.e., δ (3) (y) = 0 unless u(x) • y = 0.The first term on the right-hand side of Eq. (2.10) arises because the velocity at midpoint u(x) on the left-hand side is different from the velocity u(x 1 ) involved in the equation of motion for ϕ(x 1 ) by the amount proportional to y 1 .

B. Power counting and perturbation theory
Evolution equations (2.10) form an infinite system of equations for an infinite set of correlation functions.However, in the regime of applicability of hydrodynamics, there is a natural hierarchy in this system which allows us to systematically truncate and solve these equations.In this section we discuss this hierarchy and the small power counting parameter(s) which control it (see also Refs.[20,22]).
Hydrodynamics is an effective theory that emerges when the thermal state of the system is sufficiently homogeneous.More precisely, hydrodynamic variables (i.e., conserved densities) characterizing the local thermal state of the system must vary on a scale ℓ much longer than the microscopic scale ℓ mic , which typically is of order of the correlation length ξ. 5 Indeed, without such a scale separation one could not describe the state of the system entirely in terms of thermodynamic (conserved) quantities and, e.g., the thermodynamic concept of equation of state would lose meaning.The scale separation ℓ ≫ ℓ mic is the essential property of the hydrodynamic regime and it gives rise to a small parameter ℓ mic /ℓ (Knudsen number) which controls the gradient expansion of hydrodynamic equations and ensures their locality.
When we consider fluctuations in hydrodynamics the relevant characteristic inhomogeneity scale is the wavelength of the fluctuations, or the inverse of their wave number q, and the corresponding small parameter is ε q ≡ qℓ mic .
The small parameter which controls the magnitude and importance of fluctuations is closely related to ε q .This parameter is the inverse of the typical number of the uncorrelated microscopic cells (of volume ξ 3 ) in the homogeneity region (of volume ℓ 3 , or 1/q 3 ): ε ≡ (ξ/ℓ) 3 ∼ (qξ) 3 .The central limit theorem guarantees that (1) the magnitude of the fluctuations averaged over the scale of ℓ ∼ 1/q is suppressed by √ ε and (2) the non-Gaussian connected correlators of order n > 2 are suppressed by ε n−1 .
Both ε q and ε are small in the hydrodynamic regime, when qℓ mic ≪ 1, and ε ∼ ε 3 q .However, for the sake of generality, and to emphasize the difference between the gradient expansion corrections and the fluctuation corrections to hydrodynamics, we shall treat these two parameters as independent. 6he corresponding power counting for various quantities involved in our calculation is as follows: ) where [n/2] = n/2 for even n and (n + 1)/2 for odd n, while G c denote connected correlation functions defined below in Eq. (2.17) or (2.20).The factor ( √ ε) n in G is a consequence of the suppression of the magnitude of fluctuations averaged over hydrodynamic scale by a factor √ ε.The factor ε n−1 in G c results from the requirement that n − 1 internal correlations are needed for a fully connected n-point correlator (with each internal two-point correlator being of order ε).The factor ε q in L reflects the fact that the rate of ideal hydrodynamic evolution for fluctuations is first order in spatial derivatives, while the O(ε 2 q ) term in L represents dissipative terms in hydrodynamic equations which are second order (as in diffusion).The ε q term is absent for purely diffusive variables.The factor ε 2 q in the Onsager operator Q controls the magnitude of the noise and is determined by the fluctuation-dissipation theorem, while the factor ε is due to the locality (short range correlation) of the noise and its suppression on the hydrodynamic scale by the central limit theorem.
The first few equations (with n = 2, 3, 4) of Eq. (2.10) are truncated in the double expansion of ε and ε q , i.e., ) ) and equations for n ⩾ 5 can be obtained accordingly.Here we only consider the first-order hydrodynamics.Thus the equations for all multipoint functions are truncated at order ε 2 q .7While we keep only the leading order terms in ε in Eqs.(2.14b) and (2.14c), both leading and next-to-leading order terms are kept in Eq. (2.14a), for the purpose of deriving the evolution equation for the four-point connected function, Eq. (2.16c).
Our purpose is to derive the evolution equations for connected correlation functions, which can be directly related to the correlations of particles in experiments [31].To this end we would need to express the correlation functions G in terms of the connected correlation functions G c (and vice versa).The equations relevant for our calculation of evolution equations up to n = 4 are given by (2.15) Using Eq. (2.14) and (2.15) and keeping the leading order terms in ε, we obtain the evolution equations for the "raw" connected functions for n = 2, 3, 4: ) (2.16c) The diagrammatic representation of Eqs.(2.16) is shown in Fig. 1.We find that the leading terms (of order ε 2 q ε n−1 ) are represented by connected tree diagrams.
( ) ( ) ( ) n either one new leg to filled vertices or one new G2 correlator to open vertices.This recursive rule can be readily verified using Eq.(2.23); see also Fig. 2. Similar diagrammatic representation applies to evolution equations (3.17) we derive below.

C. Arbitrary n
Our derivation can be extended beyond n ⩽ 4 to higher-order n-point functions, which are also of interest from the experimental point of view [32].In this subsection we discuss generalization of Eqs.(2.15) and (2.16) to multipoint connected functions of arbitrary order n.In this work we do not need to use such general-n equations, but these equations could facilitate the extension of our results to higher n in future research.They also shed more light on the structure of the n ⩽ 4 equations we do use.
To start, we need to express correlators G in terms of connected correlators G c , similarly to Eq. (2.15).The correlator G i1...in is a sum of all possible nonequivalent products of G c with indices i 1 . . .i n divided between the G c factors in all possible ways.As in Eq. (2.15) we can group these terms into sets within each of which the difference between the terms is a permutation of the indices; e.g., is such a set.Using the permutation average notation Eq. (2.8) we can represent each such a set of terms by a single term.The resulting combinatorial factors will simply count the number of terms in each such "equivalent by permutation" set.We thus find the generalization of Eq. (2.15) in the form where the inner sum is over all ordered sets of integer numbers {n 1 , . . ., in a way that each term in the sum in Eq. (2.17) is different.The factor {n 1 , . . ., n k }! is the order of the group of permutations of n indices i 1 , . . ., i n which leaves the product of G c 's in Eq. (2.17) unchanged due to the symmetry of each G c with respect to its own indices as well as the commutativity of the product of G c 's themselves.In other words, the factor counts how many identical terms of a given type appear in the sum over all n! permutations of n indices: where k m is the count of times a given integer m appears in the given set {n This ensures that all terms in Eq. (2.17) with the same k are of the same order in ε, i.e., O(ε n−k ), according to Eq. (2.13).
The inverse of Eq. (2.17) gives connected correlator G c in terms of correlators G's: (2.20)This relation can be obtained by noting that the connected correlator generating function, n µ n , and the correlator generating function, where the logarithmic function is expanded in the last equality.Equating the terms of order µ n on both sides, and noting that the contributions from ( m 1 m! G m µ m ) k to the terms of order µ n are simply given by where The diagrammatic representation of Eq. (2.23) is sketched in Fig. 2.

D. Confluent formalism for multipoint equal-time correlators
An important ingredient for our derivation of the fluctuation evolution equations in relativistic hydrodynamics is the confluent formalism, which allows us to covariantly describe fluctuations in the local rest frame of the fluid.We follow the approach introduced in Ref. [19] for two-point correlators, but generalize it here to arbitrary n-point correlators.Most of the notations in this section are similar to the ones introduced in Ref. [19] and are summarized in Appendix D for convenience.
The four-velocities u of the flow in two space-time points x and x + ∆x can be related by a Lorentz boost Λ which we define as8 Λ(∆x)u(x + ∆x) = u(x) . (2.24) For infinitesimal ∆x, the four-velocities are different by infinitesimal amount ∆u = (∆x • ∂)u and the boost can be represented by the following matrix close to unity For n-point function the fluctuations of the n variables are evaluated at different points (x 1 , . . ., x n ) with different local velocities u(x i ).Before comparing or correlating these local variables it is natural to boost all of them to the same local rest frame at the midpoint (2.26) If the variables are components of a four-vector (such as δu µ ), then the boost mixes those components using matrix Λ given explicitly in Eq. (2.25).If the variables are scalars (such as δm or δp) the boost is trivial (identity), or ∆Λ = 0. To treat these two cases simultaneously in our formalism, we introduce (as in Ref. [19]) an object u i whose components equal to u µ when i = 0, . . ., 3, and zero for all other values of i, corresponding to scalar hydrodynamic variables.Then we can express infinitesimal boost from x + ∆x to x, as a matrix where we defined the confluent connection as in Ref. [19], Boosting all variables into the local rest frame at the midpoint using Λ(x i − x), we can then define the confluent n-point correlator Ḡ and express it in terms of the "raw" correlators G: where in the second line we have kept only the leading term in the expansion in powers of x i − x.The idea of the confluent correlator is illustrated in Fig. 3.In order to describe the rate of change of a correlation function such as Ḡn with respect to the midpoint position x, we want to compare the values of Ḡn at two sets of arguments x i and x ′ i .The ordinary derivative would correspond to x ′ i = x i + ∆x with the same ∆x for all i.A more natural measure in the context of hydrodynamics, however, would take into account the change of the four-velocity between the old midpoint x and the new midpoint x ′ = x + ∆x.Specifically, we want to evaluate the function at a new set of arguments (points x ′ i ), which are located relative to the new midpoint x ′ in exactly the same way (in the sense to be precisely defined below) as they were around x in the rest frame at midpoint.This is essential if we want to evaluate the rate of change of an equal-time correlator.
That rest frame defined by the four-velocity u(x + ∆x), in general, is different from u(x).The relative position of the points can be described by four-vectors (2.30) In order to preserve the relative positions in the rest frame at the midpoint we shall define new relative positions using the same boost as in Eq. (2.24), i.e., y i (x + ∆x) = Λ(∆x) −1 y i (x).This would ensure, in particular, that the time components of the relative four-vectors y i (x) in the rest frame are preserved: Which means that if we define equal-time correlator by y i (x) • u(x) = 0, the same relation remain true at the new point: y i (x + ∆x) • u(x + ∆x) = 0.In other words, if the points x i are equal-time in the frame at their midpoint, so are the points x ′ i .Therefore, we define the confluent derivative via the following relation, where ∆x is infinitesimal: Note that the variables which are being correlated are also boosted accordingly to make sure that only their change with respect to the local rest frame is measured, and not the change of the components of these variables due to the change of the rest frame itself.Taking the limit ∆x → 0 we can write Eq. (2.31) in terms of the partial derivatives of G and the boost connection ω: At this point it is convenient to introduce the derivative with respect to the midpoint, x which was already defined in Eq. (2.10), as well as the derivatives with respect to separation vectors y i at fixed midpoint x: Note that y i variables are not independent, since n i=1 y i = 0, so the derivative ∂ (yi) has unusual, but simple, properties, e.g.,  Another derivative which will be even more useful when discussing the Wigner transform is obtained by introducing the local tetrad consisting of vector u(x) and a triad e a (x), a = 1, 2, 3, thus expressing vector y in components: y(x) = e a (x)y a + u(x)y u , where y u = −u(x) • y is the time component in the local rest frame at point x.In this paper we shall only consider equal-time correlators, so we work with y u = 0 (cf.Fig. 5(a)).We can then define a derivative where y a and y u = 0 are fixed, i.e., ∆x • ∂ Ḡi1...in ≡ Ḡi1...in (x ′ 1 , . . ., x ′ n ) − Ḡi1...in (x 1 , . . ., x n ) , where x ′ i = x + ∆x + e a (x + ∆x)y a i and x i = x + e a (x)y a i .(2.35) Taking the limit ∆x → 0 we find: where we also used Eq.(2.28) and u • y i = 0. We also defined the derivative which is essentially the derivative with respect to three-component vector y a , which, similarly to ∂ (2.40) The confluent derivative in Eq. (2.40) incorporating the two connections defined in Eqs.(2.28) and (2.37) is illustrated in Fig. 4.
Considering now a connected correlation function, G c (x 1 , . . ., x n ) (see Eq. (2.15)), and using the generalized Wigner transform introduced in Ref. [22] we define the Wigner function where q a denote components of three-vector q = {q a } ∈ R 3 with a = 1, 2, 3 which is the wave number conjugate to vector y a .To avoid clutter in our notation, we replaced the indices i 1 , . . ., i n , which are the same on W i1...in and Ḡc i1...in , with a single index n.The inverse transformation of Eq. (2.41) is given by Ḡc n (x + e a y a 1 , . . ., x + e a y a n ) = n i=1 q i 2π W n (x; q 1 , . . ., q n ) .
(2.42) Some of the features of this transformation of variables are illustrated in Fig. 5.In particular we note that, due to the constraint i y i = 0 implemented by the delta function in Eq. (2.41), the Wigner function is invariant with respect to the shift of all wave numbers q i by the same vector.That means we can constrain the value of i q i without losing any information about the dependence of W on its arguments (of course, this is related to the fact that W n has one argument more than G n ).The natural choice of the y-space x 1 x 2 x n x 5. An illustration of the vectors relevant to the Wigner transform of equal-time confluent correlators living in coordinate y-space (shown in panel (a)) to Wigner functions living in wave number q-space (shown in panel (b)).The wave numbers drawn in the figure are four-vectors e µ a q a ≡ q µ .The nature of the constraint i qi = 0 is discussed in the text.
constraint is i q i = 0 as in Eq. (2.42).An intuitive way to understand this is to think of W n as an n-point "amplitude" and of q i as the corresponding "momenta" flowing in.Then the constraint is simply a reflection of "momentum" conservation.
It is easy to see that the derivative of the Wigner function with respect to x at fixed q i 's, is the Wigner transform of ∂ Ḡc ; i.e., the derivative ∂ commutes with the Wigner transform.
It is also easy to see that the derivative defined in Eqs.(2.39) and (2.33) upon Wigner transform becomes simply the multiplication by the corresponding iq a : (2.44) We shall define ∇W as Wigner transform of ∇ Ḡc .Using Eq. (2.34) it is then easy to show that ∇µ

.45)
Now let us turn to the derivation of the evolution equation of connected Wigner functions.Applying the confluent derivative in Eq. (2.45) along u(x) to Eq. (2.41), we express the result in terms of the "raw" connected correlation functions G c as follows: The explicit form of the evolution equations for Wigner functions will depend on the explicit form of multilinear operators L and Q in Eqs.(2.16) that need to be substituted into Eq.(2.46).Below we shall obtain this explicit form for an important subset of correlators in hydrodynamics.

III. SPECIFIC ENTROPY FLUCTUATIONS IN RELATIVISTIC HYDRODYNAMICS A. Stochastic equation for specific entropy fluctuations
With the generic established, we now apply it to a more specific hydrodynamic framework governed by the local equations for the conservation of energy, momentum, and charge: supplemented by constitutive relations for the stress tensor T µν and charge current J µ in the Landau frame for the stochastic relativistic hydrodynamics [33], where ϵ and n are the energy density and charge density respectively, measured in the local rest frame; u is the four-velocity already introduced in Eq. (2.1); ∆ µν = g µν + u µ u ν is the transverse projection operator satisfying ∆ µν u ν = 0; the pressure p, appearing as the coefficient of ∆ µν in the ideal part of the stress tensor, is determined by chosen primary hydrodynamic variables (such as ϵ and n) through the equation of state.The explicit forms of the dissipative parts, denoted by Π µν and ν µ for stress tensor and charge current respectively, are determined by applying the second law of thermodynamics [2]: where α ≡ βµ ≡ µ/T is the chemical potential to temperature ratio, and gradients of velocity are decomposed into the traceless and trace parts where with a = 1, 2, 3 since there are three independent noises (random currents) in the local rest frame, and e µ a is an arbitrary spatial basis triad in the local rest frame: e a • u = 0 (as already defined and discussed in Ref. [19] and in the previous section).As a result (H I µ a )H I ν a = λ∆ µν (note that the choice of the triad e a does not matter, since e µ a e ν a = ∆ µν ).We define Sµν accordingly.Since there are six independent noise variables (corresponding to random stress tensor in local rest frame), the components of Sµν are expressed in terms of a symmetric rank-three random matrix η S ab , i.e., Sµν = HS µν ab η S ab , where ⟨η S ab (x and Similarly to Eqs. (3.5) and (3.6), isotropy requires that the random current and stress noises are statistically independent: We now turn our focus to the fluctuations of specific entropy (i.e., ratio of entropy density s to charge density n), m ≡ s/n, which is parametrically the slowest and also the most significant hydrodynamic mode near a liquid-gas critical point, as we already discussed in Section I, such a focus also underlies the approach in Ref. [29].
Specifically, we are going to derive deterministic equations of motion for the correlators of m.To achieve this, we first obtain the stochastic equation for the evolution of m starting from the conservation equation (3.1).We find where are expressed in terms of independent thermodynamic derivatives chosen in the set I given by Table I in Appendix.A, and Eq. (3.8) has a form similar to Eq. (2.1) with ψ = m.Therefore the formalism of the previous section can be applied.To obtain the equation corresponding to Eq. (2.6), we expand Eq. (3.8) in the fluctuation field, δm, as well as in derivatives applying to fluctuation fields, such as ∂δm = O(q), up to O(q 2 ). 9 It is sufficient to truncate the expansion at order δm 3 to obtain equations for the four-point connected functions. 10hile we are interested in correlators of δm, the evolution of δm also depends on fluctuations of other hydrodynamic variables, in particular, on fluctuations of pressure, δp. 11As a result, the evolution of the correlators of δm will depend on correlators of δp as well. 12As we shall see below, to the leading order in hydrodynamic (gradient) expansion and in the regime we consider, we only need to include terms linear in δp.Correspondingly, if we only write the terms which will contribute to the correlator evolution equations in the regime we consider (e.g., Eqs.(3.15) or (3.17)), we find The multilinear operators L m,m...m serving as expansion coefficients are more explicitly written as13 where In Eqs.(3.12) derivative ∂ applies to all factors to the right of it, and, since we only consider fluctuations of m, i.e., i, j 1 , . . ., j n = m and u ,m...m = 0 in Eq. (2.7), we have L m,m...m = F m,m...m .The noise introduced in Eq. (3.10) satisfies where in the absence of velocity fluctuation δu, terms involving stress noises S µν are of a higher order in the gradient expansion, i.e., they are at least of order kq ∼ ε 3 q whereas we truncated Eq. (3.14) at order q 2 ∼ k ∼ ε 2 q .
In order to close Eqs.(3.17) we need to supply the values of W mmp and W mmmp , i.e., crosscorrelators involving fluctuations of pressure.We have not written terms in which the linear (Gaussian) correlator W mp appears, because this correlator vanishes upon averaging over the timescales longer than the period of sound oscillation, as already observed in Ref. [20].The vanishing of the linear correlator W mp significantly simplifies equations (3.17).However, we cannot drop the nonlinear correlators W mmp and W mmmp , since they do not average to zero.This can be easily seen by computing (cf.Eqs.(3.23) and (A8)) the equilibrium values of W mmp and W mmmp , which are not zero, unlike the equilibrium value of W mp .
Outside of the regime we consider (i.e., at faster timescales) these nonlinear cross-correlators obey dynamic equations which will be a part of the full system of equations for fluctuation correlators.In this paper we do not intend to derive this full system.Instead, we use the fact that different correlators relax on parametrically different timescales, with correlators of the specific entropy fluctuations being the slowest.We use this hierarchy of scales and observe that pressure fluctuations (upon averaging over the sound oscillations) relax parametrically faster (i.e., on the timescale of sound attenuation) than the timescale of the evolution of the specific entropy fluctuations on which we focus.This parametric separation of scales enables the regime we consider, i.e., Hydro+ regime with only one parametrically slow nonhydrodynamic mode.Effectively, we can consider the fluctuations δm frozen when we determine the fluctuations of δp.This means that, while in complete equilibrium, by definition, δp = δm = 0, we can also consider partial equilibrium where δm is not zero, i.e., not in equilibrium, and determine the "equilibrium" value of δp under these (slowly varying) conditions.Since fluctuations δp relax faster than fluctuations δm, the value of δp becomes a function of δm, and not an independent variable, in the regime we consider.That function is the partial equilibrium value (δp) peq , which we derive in Appendix B. Consequently, the correlators W mmp and W mmmp can be expressed in terms of the specific entropy correlators as follows: where S ,... are derivatives, given by Eqs.(A4), of the entropy functional S that is maximized in equilibrium, as discussed in Appendix A. These equations are represented diagrammatically in Fig. 7. Substituting partial equilibrium values for the crosscorrelators W pmm and W pmmm from Eq. (3.22) into Eqs.(3.17) we now obtain a closed system of equations for correlators of the specific entropy. 14n order to check the validity of Eqs.(3.17 by the following space-time and wave number independent values: which, of course, can be obtained from an independent calculation based on thermodynamics (see Appendix A).This is a highly nontrivial check of the evolution equations, since it involves cancellations of many terms, including the contribution of the pressure fluctuations discussed above.
Eqs. (3.23) are expressed in terms of independent thermodynamic derivatives from set I in Table I.This is useful for verifying that equilibrium correlators in Eqs.(3.23) satisfy evolution equations.An alternative representation of Eqs.(3.23) using the independent thermodynamic derivatives from set II is given by Eqs.(C2).

IV. CONCLUSIONS AND OUTLOOK
We have generalized the fully Lorentz covariant deterministic approach to relativistic fluctuating hydrodynamics to non-Gaussian fluctuations.While the full system of equations involving correlators of all hydrodynamic variables is still a work in progress, here we demonstrate how this approach allows us to derive the relativistically covariant equations for the fluctuations of the slowest hydrodynamic mode, which is specific entropy, or m = s/n, in a fluid with arbitrary relativistic flow.Such fluctuations are of special significance near a critical point, in particular, the QCD critical point, for two reasons.First, the equilibrium fluctuations of m are the largest near the critical point: ⟨δm 2 ⟩ ∼ c p ∼ ξ 2−η .In that sense, this is the "soft" direction and the source of the most prominent critical point signatures.Second, this mode is also the slowest (not only it is diffusive, but its diffusion coefficient γ λ = κ/c p vanishes at the critical point as c p diverges) and, thus, it is the furthest from equilibrium.Therefore, the non-equilibrium dynamics of this mode of fluctuations are the most consequential from the point of view of predicting the dynamical effects on the signatures of the QCD critical point.
Similar evolution equations for non-Gaussian correlators in a static (nonflowing) fluid at given temperature were derived in Ref. [22].In this paper we consider diffusion of the slowest diffusive mode in a flowing fluid, more relevant for the description of the dynamics of critical fluctuations in heavy-ion collisions.
Comparing evolution equations (3.17) to the results of Ref. [22] we observe many similarities, which were already anticipated in Ref. [22].In fact, equations (3.17) can be obtained from those in Ref. [22] by a substitution, Eq. (3.26), as conjectured in Ref. [22].However, the terms containing derivatives of ln m n appearing in the coefficients given by Eqs.(3.19) cannot be obtained that way.Their appearance is related to the qualitative difference between the density n, whose space integral is a conserved quantity, and the specific entropy m, whose space integral (i.e., q = 0 mode) is generally not conserved.We have also found nontrivial contributions of pressure fluctuations to the evolution of the specific entropy correlators due to non-linearities of the equation of state (i.e., mode coupling).
As a very nontrivial check of the evolution equations (3.17) and (3.19), we verified that the equilibrium (thermodynamic) correlators in Eq. (3.23) solve the evolution equations.This requires multiple cancellations which cannot be achieved without the terms with derivatives of ln m n as well as the contributions of pressure fluctuations.
As we pointed out in the Introduction, our equations include only tree-level contributions (as did the equations derived in Ref. [22]).It could potentially be interesting to extend this analysis to one-loop order and consider "long-time tail" effects on the evolution of fluctuations.The resulting theory would represent the generalization of Hydro+ [29] to non-Gaussian fluctuations.
Although we limited our analysis to the slowest diffusive mode, we presented our derivation in a sufficiently general form to facilitate the extension to fluctuations of faster hydrodynamic modes.The next-to-slowest modes correspond to fluctuations of the transverse velocity of the fluid.The theory describing fluctuations of all diffusive hydrodynamic modes (specific entropy and transverse velocity) was referred to as Hydro++ in Ref. [20].Finally, the full set of hydrodynamic modes includes pressure/longitudinal velocity fluctuations, i.e., they include sound -the fastest hydrodynamic mode.We leave the development of such a full hydrodynamic theory of fluctuations to further work.
Despite the focus on the slowest hydrodynamic mode, we believe the equations presented in this paper are valuable for more realistic simulations of the dynamical evolution of fluctuations in heavy-ion collisions (as a first step one can consider extending exploratory Hydro+ calculations, such as in Refs.[34][35][36], to non-Gaussian fluctuations).The ingredients required are the same as in the usual hydrodynamic simulation: equation of state, e.g., as in Refs.[37][38][39], and kinetic coefficients.
For a system with two independent thermodynamic variables (such as m and p), there are n+1 independent nth-order thermodynamic derivatives of entropy.It is convenient for calculations to choose a basis set of independent derivatives for each order to make sure that cancellations are easier to carry out.Table I provides two such basis sets (set I and II) for derivatives up to fourth order.In set I and II we have primarily chosen the derivatives of α and ln m n , respectively.Of course, one can relate the independent thermodynamic derivatives from set I to those from set II. Taking the second-order derivative as an example, we have , (ln mn),mmm, (ln mn),mmp, (ln mn),mpp, (ln mn),ppp As a consequence, all thermodynamic derivatives can be expressed solely in terms of the independent ones such as those listed in Table I.In terms of the independent derivatives from set I given by Table I, useful expressions for some second-, third-, and fourth-order thermodynamic derivatives are presented below: FIG. 1. Diagrammatic representation of the evolution equations (2.16) for multipoint connected correlation functions.The dot on the left-hand side of each diagrammatic equation represents material derivative u • ∂.The open half-circle with two legs represents Li 1 , j 1 − δi 1 j 1 (y1 • ∂u) • (∂/∂x1), i.e., it includes the first term on the right-hand side of each equation in (2.16).It is straightforward to generate the evolution equation for G c n+1 recursively by adding in each diagram representing equation for G cn either one new leg to filled vertices or one new G2 correlator to open vertices.This recursive rule can be readily verified using Eq.(2.23); see also Fig.2.Similar diagrammatic representation applies to evolution equations (3.17) we derive below.
1 , . . ., n k }; e.g., for n = 7 set {2, 2, 3} the counts are k 1 = 0, k 2 = 2 and k 3 = 1.The denominator in Eq. (2.19) can be viewed as a symmetry factor corresponding to a diagram with k vertices (factors of G c ) with n 1 , . . ., n k equivalent legs (indices) each.In this picture, k m is the number of equivalent vertices with m legs.By definition, k m = 0 for m > n and m mk m = n.For a given k in the outer sum in Eq. (2.17) all terms must obey k 1

2 .
Diagrammatic representation of the evolution equations for generic multipoint connected correlation functions.Similarly to Fig. 1, the first term on the right-hand side of Eq. (2.23) is included into the diagram with an open semicircle with two legs.The last two terms are drawn in this figure, with all possible combinatorial arrangements at leading order.See Fig. 1 for the meaning of the diagram elements, and for the explicit representation of these diagrams when n = 2, 3, 4.

FIG. 3 .
FIG.3.An illustration, in space-time, of the process involved in constructing the confluent n-point function Ḡ from the "raw" n-point function G.All n fluctuation fields are boosted from their respective local rest frames at xi to the same frame at the midpoint x.

FIG. 4 .
FIG.4.The origin of two connections ω and ω in the confluent derivative applying to confluent n-point functions.Panel (a) explains the origin of boost connection ω, brought out by boosting the fluctuation variable ϕ(x + ∆x) from the local rest frame of u(x + ∆x) to the local rest frame of u(x), i.e., Λ(∆x)ϕ(x + ∆x), using the same boost that takes velocity u(x + ∆x) to u(x).The action of the boost on ϕ is determined by the tensor rank of ϕ and is illustrated in the figure for the case of a four-vector (e.g., velocity fluctuation).Panel (b) explains the origin of spin connection due to the confluent motion of equal-time surface and the arbitrariness of the local triad ea choice.Under the boost from the local rest frame of u(x) to u(x + ∆x) every coordinate vector y in the local equal-time three-hyperplane (represented by the blue rectangle), subject to the constraint u(x) • y = 0, must also be boosted simultaneously by Λ(∆x) −1 .This includes the basis ea(x) in that hyperplane (red vectors, only two dimensions a = 1, 2 are shown).Upon the boost the basis vectors map onto Λ(∆x) −1 ea(x) (dashed red vectors) at point x + ∆x .The additional rotation (indicated by the red circular arrow) in the equal-time plane needed to align the boosted basis (dashed red) with the local basis ea(x + ∆x) (solid red) gives rise to the ω connection.

FIG. 6 .
FIG. 6. Diagrammatic representation of the evolution equations (3.17).Similarly to Fig. 1, the solid circles with k legs represent k-point Wigner functions W k , while the open half-circles represent Lm,... (the open half-circles with two solid legs in the equation for W k represent ∼ kLm,m + (k − 1)θ), and open triangles represent Qmm,....The dot on the left-hand side of the diagrammatic equation represents the Liouville operator in Eq. (3.18).A dashed line leg corresponds to index p replacing index m in each of the above objects.In the regime we consider, the Wigner functions involving pressure fluctuations (p index, or dashed line leg) can be replaced by their partial equilibrium values given in Eq. (3.22) and represented diagrammatically in Fig. 7.Note that the diagrams with two-point correlator Wmp are not shown, since this correlator vanishes in the regime we consider.

FIG. 7 .
FIG.7.Diagrammatic representation for the partial equilibrium values of the correlators involving pressure fluctuation, as given by Eqs.(3.22).An open circle with k legs represents the derivative of entropy respect to k variables, e.g., S,mmp for k = 3, with solid and dashed lines representing derivatives with respect to m and p, respectively.The solid gray circle with two dashed line legs represents −S −1 ,pp = W eq pp , the equilibrium value of Wpp.

. 26 )
While the substitution (3.26) reproduces the terms in Eqs.(3.17)  or(3.25), the expressions for the coefficients in these equations given by Eqs.(3.19)  contain terms with derivatives of ln m n which are not reproduced by the simple map(3.26

TABLE I .
Two sets of independent thermodynamic derivatives chosen in this work.The derivative of quantity X with respect to m (or p) while keeping p (or m) fixed is denoted by Xm (or Xp).