# NON-GAUSSIAN FLUCTUATION DYNAMICS††thanks: Presented at the XXIXth
International Conference on Ultra-relativistic Nucleus-Nucleus Collisions
(Quark Matter 2022), April 4-10, 2022, Kraków, Poland.

Xin An  Department of Physics and Astronomy, University of North Carolina,
Chapel Hill, NC 27599, USA

###### Abstract

Recent progress of a general deterministic approach to the non-Gaussian
fluctuation dynamics is reviewed, with an emphasis on the derivation of the
fluctuation evolution equations and their phenomenological implication in
heavy-ion collision experiments.

##  1 Introduction

Fluctuations are ubiquitous phenomena emerging on all length scales, from the
cosmological length all the way to the size of the quark-gluon plasma created
in heavy-ion collisions, where they could potentially play essential roles. In
many situations, fluctuations are resulted from addition of independent and
random sources on smaller scales and are approximately described by the
Gaussian distribution, per the central limit theorem. In addition, the time
scale of system evolution is usually much larger than the time scale of the
fluctuation equilibration process, so that fluctuations can also be well
approximated as in equilibrium. Nonetheless, there are certain situations
where fluctuations can significantly deviate from their Gaussian and/or
equilibrium distribution, and these deviations arouse more and more interests
on various physical subjects in recent years. A typical example, which will be
focused in this work, is the event-by-event fluctuations of particle
multiplicities in heavy-ion collision experiments. The shape of the
multiplicity distribution for the observed particles (e.g., net protons),
characterized by its cumulants, carries important information of the collision
process. For example, it was argued that in equilibrium the higher-order
cumulants are more sensitive to the critical point [1, 2], which may give rise
to the non-monotonic energy dependence of the particle fluctuations in heavy-
ion collisions – an intriguing hint of the QCD critical point [3]. However,
since the fireballs, created by the collisions of two ultrarelativistic
nuclei, are short lived and rapidly expanding, it would be natural to expect
that the fluctuations may fall out of equilibrium. Thus the argument in Refs.
[1, 2] based on equilibrium theory needs to be amended to incorporate the non-
equilibrium effects of non-Gaussian fluctuations – the major motivation of
Refs. [4, 5]. In this proceeding, the progress achieved in Ref. [4] will be
reviewed in Sec. 2. Additional numerical results will be presented and
discussed in Sec. 3. These works altogether constitute an integral part of the
BEST framework [6] for interpreting the forthcoming experimental results from
the second phase of the RHIC Beam Energy Scan program.

##  2 Fluctuation evolution equations

The Langevin equations for a set of stochastic fields $\psi_{i}$’s, with the
index $i=(i_{1},i_{2},\dots)$ labeling both the variables and their associated
space coordinates (if applicable), read

$\partial_{t}\psi_{i}=F_{i}+\xi_{i}$ (1)

where $F_{i}$ and $\xi_{i}$ are the drift force and random force (noise)
respectively. The noise is assumed to be Gaussian, i.e.,

$\langle\xi_{i_{1}}(t_{1})\xi_{i_{2}}(t_{2})\rangle=2Q_{i_{1}i_{2}}\delta(t_{1}-t_{2})\,,$
(2)

and its amplitude $Q_{ij}$, known as the Onsager matrix, is set by the
fluctuation-dissipation relation. One can convert the Langevin equations (1)
into the Fokker-Planck equation for the probability distribution function of
the stochastic fields, $P(\psi_{i};t)$ [7]:

$\partial_{t}P=\left(-F_{i}P+\left(Q_{ij}P\right)_{,\,j}\right)_{,\,i}$ (3)

where $(\dots)_{,i}=\delta(\dots)/\delta\psi_{i}$.111This equation is
invariant under $Q_{ij}\to Q_{ij}+\Omega_{ij}$ where
$\Omega_{ij}=-\Omega_{ji}$ is the antisymmetric matrix obtained from the
Poisson bracket of stochastic fields [7]. The Fokker-Planck equation solves
two problems the Langevin equation suffered from: the problem of infinite
noise due to the divergence of $\delta$-function in Eq. (2) when $t_{1}\to
t_{2}$; and the problem of multiplicative noise due to the ambiguity of the
noise amplitude $Q_{ij}$ in the continuous limit $\Delta t\to 0$ where $\Delta
t$ is the discretized time increment of the evolution. The divergence in the
first problem can be regularized analytically using the renormalization
techniques [8, 9, 10], while the ambiguity in the second problem can be
eliminated by specifying the discretization prescription, e.g., Ito’s
prescription as in Eq. (3).

The cumulant generating functional, $\mathcal{W}[\mu_{i}]$, is defined through

$e^{\mathcal{W}[\mu_{i}]}=\langle
e^{\mu_{i}\psi_{i}}\rangle\equiv\int\mathcal{D}\psi Pe^{\mu_{i}\psi_{i}}\,,$
(4)

which can be expanded in a series of the external sources $\mu_{i}$’s:

$\mathcal{W}=\sum_{n=1}^{\infty}\frac{1}{n!}G^{\rm c}_{i_{1}\dots
i_{n}}\mu_{i_{1}}\dots\mu_{i_{n}}$ (5)

where the expansion coefficient $G^{\rm c}_{i_{1}\dots i_{n}}\equiv G^{\rm
c}_{n}\equiv\delta\mathcal{W}/\delta\mu_{i_{1}}\dots\delta\mu_{i_{n}}|_{\mu=0}$
is the $n$-th cumulant. Using Eqs. (3) and (4) one arrives at the evolution
equation for $\mathcal{W}[\mu_{i}]$:

$\partial_{t}\mathcal{W}=e^{-\mathcal{W}}\left(\mu_{i}F_{i}+\mu_{i}\mu_{j}Q_{ij}\right)e^{\mathcal{W}}$
(6)

where $F_{i}=F_{i}(\delta/\delta\mu_{i})$,
$Q_{ij}=Q_{ij}(\delta/\delta\mu_{i})$. Substituting Eq. (5) into (6) one can
readily obtain the evolution equation for cumulant $G^{\rm c}_{n}$. However,
these equations are not in closed form, i.e., the equations for lower-order
cumulants also depend on higher-order cumulants. Fortunately, in certain
regimes where some parameters are legitimately small, these equations can be
systematically truncated. We introduced two such parameters. The parameter
$\varepsilon$, which is inversely proportional to the uncorrelated degrees of
freedom, controls the loop expansion. Another parameter $\delta$, identified
as the Knudsen number in hydrodynamics, controls the gradient expansion.
Although each of them plays its own role, they are not independent and are
indeed related by $\varepsilon\sim\delta^{3}$ in hydrodynamics. Assigning the
following power counting:

$F_{i}\sim\delta^{2}\,,\qquad Q_{ij}\sim\delta^{2}\varepsilon\,,\qquad G^{\rm
c}_{n}\sim\varepsilon^{n-1}\,,$ (7)

and keeping the leading terms ($\sim\delta^{2}\varepsilon^{n-1}$) that turn
out to be connected tree-level diagrams, we arrive at the truncated and
iteratively solvable equations, schematically formulated in the form

$\partial_{t}G^{\rm c}_{n}=\mathscr{F}_{\rm tree}[\\{G^{\rm c}_{n},G^{\rm
c}_{n-1},\dots,G^{\rm c}_{2},\langle\psi\rangle\\}]$ (8)

and diagrammatically represented in Fig. 1. Their explicit forms can be found
in Refs. [4, 5].

Figure 1: Diagrammatic representation of Eq. (8) (or Eqs. (3)) for $n=2,3,4$
in terms of the diagrammatic ingredients introduced in this figure. The dot on
the left hand side denotes the time derivative $\partial_{t}$.

If the variables $\psi_{i_{1}},\dots,\psi_{i_{n}}$ are fields defined locally
in continuous real space at $\bm{x}_{1},\dots,\bm{x}_{n}$, it may be more
convenient to analyze Eq. (8) in the wave-vector space. To this end we
introduced a novel multi-point Wigner function [4]:

$W^{\rm
c}_{n}(\bm{x},\bm{q}_{1},\dots,\bm{q}_{n})=\mathop{\text{\Large$\int$}}\nolimits\left[\prod_{i=1}^{n}d^{3}\bm{y}_{i}e^{-i\bm{q}_{i}\cdot\bm{y}_{i}}\right]\delta^{(3)}\left(\frac{1}{n}\sum_{i=1}^{n}\bm{y}_{i}\right)G^{\rm
c}_{n}(\bm{x}_{1},\dots,\bm{x}_{n})\,,$ (9)

where $\bm{q}_{i}$ is the wave-vector conjugate to
$\bm{y}_{i}\equiv\bm{x}_{i}-\bm{x}$ and $\bm{x}=\sum_{i=1}^{n}\bm{x}_{i}/n$.
Under this Wigner transform we imposed the constraint
$\sum_{i=1}^{n}\bm{y}_{i}=0$, as a consequence $\bm{q}_{i}$’s are not
independent, and must also sum to zero after an appropriate shift. Obviously,
Eq. (9) reproduces the traditional definition of Winger function when $n=2$.
Using Eq. (9) one can represent Eq. (8) in the wave-vector space. We defer its
presentation in a specific problem discussed in Sec. 3.

For confluent formulation of non-Gaussian fluctuations in relativistic fluid,
see Ref. [5].

##  3 Non-Gaussian fluctuation dynamics of diffusive charge

In this section we apply our general formalism to a specific problem – the
evolution of diffusive charge, using the translation given by Table. 1. In
this problem the stochastic variable, charge density $n(\bm{x})$, is defined
in continuous space, the Onsager matrix $Q_{ij}$ is related to conductivity
$\lambda$ through the fluctuation-dissipation relation, and the drift force
$F_{i}$ is given by the divergence of a diffusion current whose constitutive
relation, $\lambda\bm{\nabla}\alpha$ where $\alpha$ is the chemical potential
per temperature, can be determined by the second law of thermodynamics.

quantities |  general |  diffusive charge
---|---|---
variable |  $\psi_{i}$ |  $n(\bm{x})$
variable index |  $i,j,k,$ etc. |  $\bm{x},\bm{y},\bm{z},$ etc.
Onsager matrix |  $Q_{ij}$ |  $\bm{\nabla}_{\bm{x}}\lambda\bm{\nabla}_{\bm{y}}\,\delta^{(3)}_{\bm{xy}}$
drift force |  $F_{i}$ |  $\bm{\nabla}_{\bm{x}}\lambda\bm{\nabla}_{\bm{x}}\alpha$
Table 1: Translation of general formalism to the problem of diffusive charge.

Applying the translation in Table. 1 to Eq. (8) and using Eq. (9), one
immediately obtains the evolution equations for the diffusive charge cumulants
in the wave-vector space. The first few equations for $n=2,3,4$ read (cf. Fig.
1):

$\displaystyle\partial_{t}W_{2}(\bm{q}_{1})=$
$\displaystyle-2\left[\gamma\bm{q}_{1}^{2}W_{2}(\bm{q}_{1})+\lambda\bm{q}_{1}\cdot\bm{q}_{2}\right]_{\overline{12}}\,,$
$\displaystyle\partial_{t}W_{3}(\bm{q}_{1},\bm{q}_{2})=$
$\displaystyle-3\left[\gamma\bm{q}_{1}^{2}W_{3}(\bm{q}_{1},\bm{q}_{2})+\gamma^{\prime}\bm{q}_{1}^{2}W_{2}(\bm{q}_{2})W_{2}(\bm{q}_{3})+2\lambda^{\prime}\bm{q}_{1}\cdot\bm{q}_{2}W_{2}(\bm{q}_{3})\right]_{\overline{123}}\,,$
$\displaystyle\partial_{t}W^{\rm c}_{4}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3})=$
$\displaystyle-4\left[\gamma\bm{q}_{1}^{2}W^{\rm
c}_{4}(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3})+3\gamma^{\prime}\bm{q}_{1}^{2}W_{2}(\bm{q}_{2})W_{3}(q_{3},\bm{q}_{4})+\gamma^{\prime\prime}\bm{q}_{1}^{2}W_{2}(\bm{q}_{2})W_{2}(\bm{q}_{3})W_{2}(\bm{q}_{4})\right.$
$\displaystyle\left.+3\lambda^{\prime}\bm{q}_{1}\cdot\bm{q}_{2}W_{3}(\bm{q}_{3},\bm{q}_{4})+3\lambda^{\prime\prime}\bm{q}_{1}\cdot\bm{q}_{2}W_{2}(\bm{q}_{3})W_{2}(\bm{q}_{4})\right]_{\overline{1234}}\,,$
(10)

where $\gamma=\lambda\alpha^{\prime}$,
$\alpha^{\prime}\equiv\partial\alpha/\partial n$ and $\overline{1\dots n}$
denotes the sum over all $n!$ permutations of $\bm{q}_{1},\dots,\bm{q}_{n}$
divided by $n!$. Eqs. (3) are solved by $W_{2}^{\rm eq}=1/\alpha^{\prime}$,
$W_{3}^{\rm eq}=-\alpha^{\prime\prime}/\alpha^{\prime 3}$, $W_{4}^{\rm
c,eq}=(3\alpha^{\prime\prime
2}-\alpha^{\prime}\alpha^{\prime\prime\prime})/\alpha^{\prime 5}$ in
equilibrium, as expected from thermodynamic calculations.

In the critical regime where the correlation length $\xi$ is still much less
than the fluctuation scale $q^{-1}\equiv|\bm{q}|^{-1}$ but becomes much larger
than all other microscopical lengths such as the inverse of temperature, i.e.,
$T^{-1}\ll\xi\ll q^{-1}$, all thermodynamic and transport quantities will
acquire their dependence on the correlation length (which serves as a UV
cutoff of fluctuations). In the dynamical universal class of Model H, we have
approximately [5]

$\alpha^{(n)}\sim\xi^{\frac{n-5}{2}}\,,\qquad\lambda^{(n)}\sim\xi^{\frac{n+2}{2}}\,,\qquad\gamma^{(n)}\sim\xi^{\frac{n-2}{2}}\,,\qquad
W^{\rm c}_{n}\sim\xi^{\frac{5n-6}{2}}\,,\qquad\partial_{t}W^{\rm c}_{n}\sim
q^{2}\xi^{\frac{5n-8}{2}}\,.$ (11)

The last relation in Eq. (11) says that, all terms presented in each equation
of (3) are in the same power of correlation length, demonstrating their equal
importance near the critical point.

To see how Eqs. (3) work in practice, let’s consider the problem of baryon
charge cumulant evolution in the crossover region near the QCD critical point
(see Fig. 2). Our purpose would be illustrating the robust phenomenological
consequence of Eqs. (3). To this end, we postulate the following assumptions
which simplify our simulations without losing the main features of the
results. First, when establishing the dependence of the correlation length
$\xi$ on baryon chemical potential $\mu$ and temperature $T$, we use a mapping
$\xi(\mu,T)=\xi_{\rm Ising}(r(\mu,T),h(\mu,T))$ from Ising model where $r$ is
the Ising temperature and $h$ is the Ising magnetic field, and assume the
mapping is linear and orthogonal around the critical point $(\mu_{\rm
c},T_{\rm c})$, i.e., $\mu-\mu_{\rm c}\sim r$ and $T-T_{\rm c}\sim h$. Second,
the evolution trajectories, rather complicated in reality, are assumed to be
perpendicular to the crossover line (see the right panel of Fig. 2). Third,
the system is assumed to be in equilibrium at the initial time $t_{i}$, i.e.,
the moment when the system just enters the critical region, and assumed to
freeze out before the final time $t_{f}$, i.e., the moment when the system
leaves the critical region. Fourth, the system is assumed to be isotropic such
that $W^{\rm c}_{n}$’s depend only on $n(n-1)/2$ independent scalar
invariants, $\bm{q}_{i}\cdot\bm{q}_{j}$, $i<j$, which are all set to $q^{2}$
for simplicity [4].

Figure 2: QCD phase diagram (left) and the zoomed-in critical region near the
QCD critical point (right). The black dashed and solid curve is the crossover
and first-order phase transition line. The white lines represent the evolution
trajectories of each collision event mapped to the QCD phase diagram, with
arrows indicating the time evolution direction. The yellow dashed curve is the
freeze-out line where observables are measured. The critical region studied in
Figs. 3 and 4 is bounded by the red square.

Simulating Eqs. (3) using the above assumptions, we find the contour maps of
cumulants $W^{\rm c}_{n=2,3,4}$ in the critical region, shown in Fig. 3. The
top panel shows the contours of cumulants in equilibrium, which are symmetric
to the crossover line due to aforementioned simplifications. The bottom panel
shows a particular example (i.e., $q=0.1$)222Here we used arbitrary units,
which are not specified thereby. where the cumulants are out of equilibrium.
In this case the contours are significantly distorted as if they are dragged
toward the time evolution direction. In other words, although the dynamic
cumulants keep approaching their instantaneous equilibrium values, they still
retain memories of their past evolution history.333Similar study can be found
in Ref. [11] for the homogeneous mode that is independent of the scale $q$.

Figure 3: Contour maps of cumulants $W^{\rm c}_{n=2,3,4}$ in equilibrium (top
panel, $q\to\infty$) and out of equilibrium (bottom panel, $q=0.1$).

Focusing on the fourth cumulant $W^{\rm c}_{4}$, we plot the contours for
different $q$’s in Fig. 4. We find that cumulants are closer to equilibrium on
smaller scale (e.g., $q=0.3$), while retain a longer memory of past on larger
scale (e.g., $q=0.1$). This is a consequence of causality: on larger spatial
scale, fluctuations need longer time to equilibrate. When $q\to 0$ (infinitely
large scale), fluctuations will be completely suppressed due to the
conservation of charges, and remains as if they were at initial time (which,
of course, are also suppressed).

The fourth cumulant $W^{\rm c}_{4}$ is of particular interests to us, since
its equilibrium value possesses more sensitivity to the critical point compare
to the lower-order cumulants. If one reads off its value along the freeze-out
line (yellow dashed curve) in Fig. 4, one finds its expected non-monotonicity
as function of baryon chemical potential in Fig. 5. Noting that the baryon
chemical potential is inversely related to the collision energy, Fig. 5 is
qualitatively consistent with the experiment data reported in Ref. [3]. More
interestingly, one finds that the curve for smaller $q$ (red colored) might be
shifted toward large baryon chemical potential more significantly compared to
the one for larger $q$ (blue colored). That is to say, the non-equilibrium
effect, largely coming from the fluctuations with wave-vectors around $q\sim
0.1-0.3$, might shift the curve predicted by equilibrium theory toward the
lower collision energy. As a reminder, one should not compare the values here
quantitatively, as our results are only for illustration purpose.

Figure 4: Contour maps of fourth cumulant $W^{\rm c}_{4}$ at $q\to\infty$
(left), $q=0.3$ (middle) and $q=0.1$ (right) in the critical region. The
freeze-out line (yellow dashed curve) is arbitrarily chosen for illustration
purpose.

Figure 5:  The dependence of $W^{\rm c}_{4}$ on baryon chemical potential
$\mu$ along the freeze-out line in Fig. 4, at $q\to\infty$ (equilibrium),
$q=0.3$ and $q=0.1$ respectively.

##  4 Summary

In this proceeding, we reviewed a general deterministic approach to non-
Gaussian fluctuation dynamics. This novel approach is established with
controllable perturbative loop/gradient expansion and multi-point Wigner
function. In this approach we derived the general evolution equations for
cumulants, and demonstrated their robust phenomenological implication via the
problem of diffusive charge near the QCD critical point. Our numerical results
suggest that the non-equilibrium effect to cumulants needs to be taken into
account for a rigorous quantitative comparison with the experimental results.
As for the future, the cumulant evolution equations need to be implemented in
a more realistic setup, and our formalism is yet to incorporate all non-
Gaussian hydrodynamic fluctuations with their freeze-out prescriptions [5, 8,
9, 12].

## References

  * [1] M. Stephanov, “Non-Gaussian fluctuations near the QCD critical point,” Phys. Rev. Lett. 102 (2009) 032301, arXiv:0809.3450 [hep-ph].
  * [2] M. Stephanov, “On the sign of kurtosis near the QCD critical point,” Phys. Rev. Lett. 107 (2011) 052301, arXiv:1104.1627 [hep-ph].
  * [3] STAR Collaboration, J. Adam et al., “Nonmonotonic energy dependence of net-proton number fluctuations,” Phys. Rev. Lett. 126 (2021) 092301, arXiv:2001.02852 [nucl-ex].
  * [4] X. An, G. Başar, M. Stephanov, and H.-U. Yee, “Evolution of non-gaussian hydrodynamic fluctuations,” Phys. Rev. Lett. 127 (2021) 072301, arXiv:2009.10742 [hep-th].
  * [5] X. An, G. Başar, M. Stephanov, and H.-U. Yee, “Non-Gaussian fluctuation dynamics in relativistic fluid.” to appear.
  * [6] X. An et al., “The BEST framework for the search for the QCD critical point and the chiral magnetic effect,” Nucl. Phys. A 1017 (2022) 122343, arXiv:2108.13867 [nucl-th].
  * [7] J. Zinn-Justin, “Quantum field theory and critical phenomena,” Int. Ser. Monogr. Phys. 113 (2002) 1–1054.
  * [8] X. An, G. Başar, M. Stephanov, and H.-U. Yee, “Relativistic Hydrodynamic Fluctuations,” Phys. Rev. C 100 (2019) 024910, arXiv:1902.09517 [hep-th].
  * [9] X. An, G. Başar, M. Stephanov, and H.-U. Yee, “Fluctuation dynamics in a relativistic fluid with a critical point,” Phys. Rev. C 102 (2020) 034901, arXiv:1912.13456 [hep-th].
  * [10] X. An, “Relativistic dynamics of fluctuations and QCD critical point,” Nucl. Phys. A 1005 (2021) 121957, arXiv:2003.02828 [hep-th].
  * [11] S. Mukherjee, R. Venugopalan, and Y. Yin, “Real time evolution of non-Gaussian cumulants in the QCD critical regime,” Phys. Rev. C 92 (2015) 034912, arXiv:1506.00645 [hep-ph].
  * [12] M. Pradeep, K. Rajagopal, M. Stephanov, and Y. Yin, “Freezing out fluctuations in Hydro+ near the QCD critical point,” arXiv:2204.00639 [hep-ph].