Evolution of Non-Gaussian Hydrodynamic Fluctuations

In the context of the search for the QCD critical point using non-Gaussian fluctuations, we obtain the evolution equations for non-Gaussian cumulants to the leading order of the systematic expansion in the magnitude of thermal fluctuations. We develop a diagrammatic technique in which the leading order contributions are given by tree diagrams. We introduce a Wigner transform for multipoint correlators and derive the evolution equations for three- and four-point Wigner functions for the problem of nonlinear stochastic diffusion with multiplicative noise.


INTRODUCTION
The recent resurgence of interest in the classic subject of hydrodynamics in general [1] and hydrodynamic fluctuations [2] in particular has been largely driven by the progress in heavy-ion collision experiments that create and study droplets of hot and dense matter governed by the physics of strong interactions described by quantum chromodynamics (QCD). The importance of fluctuations in heavy-ion collisions is due to the fact that the QCD fireballs created in such experiments, with typical particle multiplicities O(10 2−4 ), while being large enough for hydrodynamics to apply, are not too large for fluctuations to be negligible.
Such fluctuations are observable in heavy-ion collisions via event-by-event measurements. Furthermore, fluctuations are enhanced if the matter created in the collisions is in a state close to a critical point and can serve as signatures of the criticality [3][4][5][6] in the beam energy scan experiments [7,8]. The magnitude of the signatures is determined by the competition between the critical slowing down and the finiteness of the expansion time [4,[9][10][11].
This necessitates a quantitative description of the fluctuation evolution within a hydrodynamic framework, and there have been significant advances in that area recently [12][13][14][15][16][17][18][19][20][21][22][23][24][25]. Most relevant for this work is the formalism describing the evolution of correlation functions coupled to the hydrodynamic background. While the approach was considered long ago in a nonrelativistic context [26], the relativistic formalism has been introduced recently in the boost-invariant Bjorken flow characteristic of heavy-ion collisions in Refs. [16,17], and in general background in Refs. [20,21].
However, so far the formalism has been limited to two-point correlation functions. The description of the higher-point correlators quantifying the non-Gaussianity of the fluctuations has been elusive until now. On the other hand, the experimental search for the QCD critical point relies heavily on such measures of non-Gaussianity [5][6][7][8] (which, similar to fluctuation magnitudes, depend on time evolution [10]). We tackle this crucial gap between the ability of the theory w.r.t. non-Gaussian fluctuations and the demand of the experiment.
We obtain evolution equations for appropriate measures of non-Gaussianity in the hydrodynamic regime, that is, the regime where the ratio of correlation length to typical fluctuation wavelength is small. Such a regime exists even near the critical point, provided the correlation length remains much shorter than the size of the system, as is the case in most physical systems, including heavy-ion collisions.

GENERAL MULTIVARIABLE FORMALISM
To understand better the issues of nonlinearity and multiplicative noise, which is essential for non-Gaussian fluctuations, we begin with a more general formalism for a discrete set of stochastic variablesv i , labeled by index i. The set of stochastic Langevin equations reads where drift F and noise magnitude H are functions ofv i , summation over repeated indices is implied, and η i is the Gaussian white noise, i.e., Eq. (1) suffers from a well-known ambiguity often referred to as the problem of multiplicative noise: it needs additional information to define the product of the stochastic function H[v] and the noise η. From our point of view, this is not really a problem but rather a shortcoming of the notation used in Eq. (1), which does not reflect the necessary information. The ambiguity is removed by discretizing time in Eq. (1) and taking the limit ∆t → 0. The choice of the definition of the product of H and η is a matter of convenience. We shall use the choice known as Ito calculus, where H and η are evaluated at the same moment. It is also well-known that the change of the definition of stochastic calculus (e.g., Stratonovich instead of Ito) is simply equivalent to a shift of the drift term F i [27].
Using Ito calculus, one can derive a Fokker-Plank equation for the probability distribution ofv: where Q ≡ HH T , (. Let us consider the equilibrium solution P eq to the Fokker-Plank equation, i.e., ∂ t P eq = 0. While the divergence of the probability flux on the rhs of Eq. (3) vanishes, the flux itself does not have to and could be equal to the divergence of an antisymmetric 2-form, which we write as where Ω ij = −Ω ji . Thus, F i can be expressed in terms of the equilibrium distribution F i = P −1 eq ((Q + Ω) ij P eq ) ,j = M ij S ,j + M ij,j , (5) where we introduced the Onsager matrix M ≡ Q + Ω and also S ≡ log P eq .
Rather than defining the stochastic process in terms of the function F i , which additionally needs specification of the stochastic calculus rule, it makes more sense to use S and M to define the process [28]. While S describes the equilibrium distribution, and in a certain sense is analogous to entropy, the matrix M describes the dynamics of the stochastic process: the symmetric semipositive definite part Q is responsible for relaxation and the antisymmetric matrix Ω (symplectic form) for the Hamiltonian-like nondissipative motion. These physical properties of the process are independent of the stochastic calculus prescription, while their relationship to drift F depends on such a prescription and is written in Eq. (5) for Ito calculus.
For example, if S is a bilinear function of v i and M ij are constants, (i.e., when Eq. (1) describes a linear Ornstein-Uhlenbeck process) the equations for G n are linear and involve only G n and G n−2 , and therefore can be solved iteratively. The equations for the cumulants, are even simpler -they decouple from each other. Equilibrium is achieved when all n > 2 cumulants vanish, as expected for the Gaussian distribution P eq = e S , while G eq 2 = −(S ) −1 where (S ) ij ≡ S ,ij . In general, S is not a bilinear and M ij are not constants, and we have an infinite system of coupled equations for cumulants. To organize this system into a hierarchy, we develop a perturbation theory.
In many physical systems of interest, particularly in hydrodynamics, the fluctuations around the equilibrium are controllably small. In other words, the probability distribution P is sharply peaked and can be treated as Gaussian in the lowest order of an approximation.
To be systematic, we introduce an expansion parameter ε to control the magnitude of the deviations from equilibrium. We assume S ∼ ε −1 (8) and that this remains true for all higher-order derivatives of S, ensuring that for small ε the probability distribution e S approaches a narrow Gaussian with the characteristic magnitude of fluctuations The odd-order moments G 2n−1 are of the same order as G 2n because δv = 0. On the other hand, the non-Gaussian cumulants are smaller for the same order: This power counting is easily established by a diagrammatic expansion of Eq. (7) using probability e S , with (S ) −1 ∼ ε playing the role of a propagator and each vertex being of order ε −1 .
Because, according to Eq. (9), cumulants of higher orders are progressively suppressed, the hierarchy is now parametrically controllable by ε. Truncating each equation at the leading order, we find, for n = 2, 3, ..in denotes the sum over permutations of indices. The corresponding equation for ∂ t G c i1i2i3i4 contains nine terms (before index permutations), and it is easier to express using diagrammatic representation in Figs. 1 and 2. The symmetry factors, such as 1/2 in Eq. (10b), count the number of permutations that do not produce different terms.
It is notable that the only diagrams appearing at the leading order in ε are "trees." The first term in each equation is the same as in the linear problem, as it involves only S and no derivatives of M . Other terms are due to nonlinearities and/or multiplicative noise, which contribute at the same order in ε. Higher-order terms contain loop diagrams, which we shall not discuss in this Letter.
In most applications, the noise η is almost Gaussian, being a cumulative result of many uncorrelated factors. In our approach, one can treat noise non-Gaussianity systematically when H is a small parameter and the noise cumulants obey the same hierarchy, as in Eq. (9). These assumptions hold in hydrodynamics. The n-th (n > 2) noise cumulant contribution will enter at order H n , while the leading terms we keep in Eq. (10) are of order H 2 . Eq. (10) and its diagrammatic representation are among the main results of this work. We shall apply this approach to hydrodynamics, i.e., a system of stochastic continuum fields, rather than discrete

STOCHASTIC NONLINEAR DIFFUSION
As the simplest example of a stochastic hydrodynamic problem, we shall consider the diffusion of conserved density (of, e.g., particle number or charge) in a slowly evolving (e.g., expanding and/or cooling) medium. This problem carries the most important features we want to address in this work, including multiplicative noise, without the complications of additional degrees of freedom and the Hamiltonian (nondissipative) dynamics of ideal fluid motion. We shall postpone the extension of these results to full stochastic hydrodynamics to future work but will discuss the implications of our findings below.
The diffusion equation is essentially a conservation equation for fluctuating densityn: where the constitutive relation is with the stochastic noise η given by λ = λ(n) is conductivity andα = α(n) is chemical potential (in units of temperature).
To translate the generic multivariable stochastic system into a stochastic hydrodynamic diffusion problem, we use the following dictionary: where α = −∂s/∂n andᾱ is a constant, understood as the chemical potential (in units of temperature) of the particle reservoir controlling the average particle density (similar to a Lagrange multiplier). Since the indices i, j, etc. become continuous coordinates x, y, etc., matrices M and H become operators, which in Eq. (13) are expressed in terms of their integral kernels. They are local, i.e., could be also written as differential operators, i.e., H = ∇ √ λ and M = −∇λ∇. Note that M = HH T , i.e., Ω = 0 for the diffusion problem.
As in Eq. (1), because the noise is multiplied by a function of the stochastic variablen, the stochastic equation [Eq. (11)] is not well-defined as written. As before, we can define the problem by specifying the equilibrium distribution or, equivalently, the thermodynamic entropy S, as well as the Onsager matrix/operator M , given in Eq. (13) [29].
The role of expansion parameter allowing us to organize and systematically truncate the infinite system of coupled equations for correlation functions is played by a certain ratio of scales. Hydrodynamics describes long wavelength modes of the field n x , characterized by wave numbers q ξ −1 , where ξ is the microscopic correlation length. The noise is local, i.e., correlated on the scale ξ, and thus its effect on the long wavelength modes involves averaging over a large number O(qξ) −3 of uncorrelated cells, suppressing fluctuations (see, e.g., Refs. [20,21]). The small parameter (qξ) 3 plays a role similar to parameter ε in Eq. (9). The smallness of H, suppressing the contribution of noise non-Gaussianity, is due to the gradient in conservation equation [Eq. (11a)] and is also controlled by qξ. Near the critical point, ξ becomes large, but as long as ξ is much smaller than the system size, hydrodynamic fluctuation theory applies.
Applying the dictionary in Eq. (13), we can derive evolution equations for the n-point correlators The equations for G n can then be converted into corresponding equations for multipoint Wigner functions after we introduce and define these objects.

MULTIPOINT WIGNER TRANSFORM
In hydrodynamics we consider fluctuations on a smoothly varying background (see, e.g., Refs. [16,20,21]. As a result, two separate scales characterize fluctuation correlators G n . A shorter scale, corresponding to wave number q, characterizes the dependence on the separation between points, while the dependence on the midpoint position occurs at a much longer scale. The well-known method to take advantage of such a scale separation in a two-point function is to work with the Wigner transform (as in the derivation of kinetic theory). In order to do this for n-point functions, we need to generalize the Wigner transform, thus far only known for two-point functions.
We propose to define the symmetric generalization of the Wigner transform and its inverse as follows: Note that, because of the δ function in Eq. (15a), the Wigner function W n is invariant under the shift of all momenta q i by the same vector. Effectively, there are only n − 1 nonredundant wave vector arguments in W n (i.e., the total number of nonredundant arguments is the same as for G n itself). Therefore, it is not surprising that, in Eq. (15b), to obtain G n we only need to evaluate W n for a set of q i 's that sum to zero. The same is true for all expressions that follow.
As an example, consider the case n = 2. Then This is the usual Wigner transform. Because one of the q i 's is always redundant, we shall adopt a simplified notation for W n by dropping the redundant argument, as in Eq. (16). Since, in this work, we consider Wigner functions symmetric w.r.t. their arguments, it does not matter which argument is dropped. Also note that the Wigner transform of a generalized n-point δ function is unity for all n.
The Wigner transform of the partial derivative of G n is given by In hydrodynamics, the gradient term is subleading to the term proportional to q i . As a result, partial derivatives in the equations for G n turn into factors of q and partial differential equations become ordinary differential equations. To simplify notations, below we shall omit the spatial x and the time t arguments for the Wigner functions.
Applying the generalized Wigner transform to the evolution equations for G n we arrive at the following evolution equations for n = 2, 3, 4: where γ = λα . Eq. (18) is the main result of this work. It is easy to map Eq. (18) to diagrams in Figs. 1 and 2 by using the dictionary in Eq. (13) and noting that S ,ij → −α , S ,ijk → −α , M ij → −λq 1 · q 2 , M ij,k → −λ q 1 · q 2 , etc. Eq. (18) describes the relaxation of the Wigner functions to their equilibrium values which, as one can check, agree with thermodynamics. Consistent with the underlying conservation laws, the longer wavelength modes relax more slowly [30] at a rate proportional to the square of their wave number.
Eq. (18) bears a certain similarity to equations in Ref. [10] for cumulants in a uniform finite size system (with 1/q in Eq. (18) and the system size in Ref. [10] playing similar roles). However, the terms with derivatives of λ, related to multiplicative noise, are absent in Ref. [10], where λ is a constant. These terms are not negligible in general and, in particular, near the liquid-gas-type critical point where λ is divergent [31].
We considered the problem of the diffusion of a conserved quantity and left the generalization to full stochastic hydrodynamics, including pressure and flow, to future work. However, we can anticipate some features of this full system based on what we have already learned from the diffusion problem. If we focus on the fluctuations of the slowest hydrodynamic mode, that is, entropy per charge, m ≡ s/n, at constant pressure, we should find a similar diffusion equation with the substitution where κ is thermal conductivity and c p is heat capacity at constant pressure. The equation for a two-point function δmδm was derived by two different methods in Refs. [15,21]. As expected, in the regime we consider, qξ 1, it coincides with our Eq. (18a) upon substitution Eq. (19).
One can argue that this should hold true for higherpoint correlators as well. It would be interesting and important for applications to heavy-ion collisions to combine the approach to non-Gaussianity presented here with the relativistic treatment of fluctuations in Ref. [21]. One can anticipate that, in addition to Eq. (19), the time derivatives will be replaced by corresponding Liouville operators. Additional correlators will appear but, near the critical point, the slowest and thus most out-of-equilibrium fluctuations will be described by cumulants of m. The validity of this argument should be checked by direct derivation, which we defer to future work.

CONCLUSIONS
We found a systematically controllable hierarchy of equations describing the evolution of higher-order, non-Gaussian cumulants of fluctuations in a general multivariable stochastic system and introduced a convenient diagrammatic representation.
We used this approach to tackle the problem of stochastic nonlinear diffusion with densitydependent conductivity, which involves multiplicative noise. We introduced a generalization of a Wigner transform to multipoint correlation function [Eq. (15)] that allows us to take advantage of the separation of scales in hydrodynamics and obtain evolution equations [Eq. (18)]. The equations for the full system of hydrodynamic variables can be derived along the same lines, and we defer this to future work.
It would also be interesting to extend the evolution equations [Eqs. (10) and (18)] beyond the leading order to loop diagrams and study the effects of ultraviolet renormalization and long-time tails, as in, e.g., Ref. [21], now involving multipoint correlations and multiplicative noise.
The formalism and the results we present are very general and would pertain to problems where nonlinearity and non-Gaussian fluctuations are of interest -from cosmology and astrophysics [32][33][34] to the physics of electronic devices [35][36][37], glassy, granular, or colloidal materials [38], and even finance [39], to name only a few examples. Among the most immediate and practical applications is the description of the evolution of the non-Gaussian measures of fluctuations in heavy-ion collisions, which is crucial for the ongoing QCD critical point search.
We thank M. Hongo, D. Pilipovic and Y. Yin for helpful discussions and comments. This work is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration and grant No. DE-FG0201ER41195.

Numerical solution
To illustrate the main features of the evolution of Wigner functions we solve Eq. (18) numerically in a simplified setup relevant to critical point search in heavy-ion collisions. We consider the system that expands along a trajectory which passes through the crossover region near a critical point. The main features of the equation of state are due to the singular behavior of the correlation length which we model approximately by an ansatz which satisfies wellknown scaling properties, i.e., ξ = r −ν f (M r −β ). A simple model of the scaling function, f = (1 + x 2 ) −ν/2β , will be more than sufficient for our purposes. The Ising variables r (reduced temperature) and M (magnetization) are given by linear functions of T − T c and n − n c respectively. We model expansion by assuming that n is a linear function of time t, choosing n(t = 0) = n c . We also use the scaling form of the inverse susceptibility α ∼ ξ −γ/ν and conductivity λ ∼ ξ x λ with well-known critical exponents. We use arbitrary units, since our purpose is to illustrate the major qualitative features.
To further simplify the calculation without losing the main features of the results, we consider an isotropic system such that W n depend only on scalar invariants. We note that the number of such independent invariants of n vectors q i which obey a constraint i q i = 0 is n(n − 1)/2. We can choose them to be q i · q j , i < j, since q 2 i = − j =i q i · q j . We shall demonstrate the dependence of the Wigner functions on time and on the scale of q. For that purpose we shall plot solutions for a set of momenta for which all independent invariants q i ·q j are equal to the same value q 2 . The numerical solutions of Eq. (18) for two different representative values of q are shown in Fig 3. One readily observes two characteristic features. First, due to the relaxation nature of the equations, the evolution of W n lags behind the change of their equilibrium values (dashed curves) driven by the expansion of the system through the critical region. Similar "memory" effects have been observed in various models for cumulants in Refs. [9,10]. Second, unlike previous approaches, Eq. (18) also describes the dependence of the memory effects on the wave number q of the fluctuations. As expected from the conservation laws underlying hydrodynamics, the longer wavelength modes relax slower, i.e., retain "memory" longer. This behavior has important consequences for observation of critical fluctuations in heavy-ion collisions. The strength of the observed signal depends on the fluctuation scale probed. A more quantitative analysis of such effects will require application of Eq. (18) in a more realistic setting of a heavy-ion collision and we leave it to future work.