#  Relativistic stochastic hydrodynamics from quantum nonlinear projection
operator

Jin Hu<EMAIL_ADDRESS>Department of Physics, Fuzhou University, Fujian
350116, China

###### Abstract

We systematically derive relativistic stochastic hydrodynamics based on the
method of quantum nonlinear projection operator. Morozov’s nonlinear
projection operator is a generalization of the well-known linear Mori-Zwanzig
projection operator method, from which one can account for the nonlinear
interaction between macroscopic modes. The quantum generalized Fokker-Planck
and Langevin equations are also obtained using this formalism, which are
fundamentally important in non-equilibrium statistical physics. As an
application, the relativistic stochastic hydrodynamic equations with Gaussian
noises are derived, which are applicable in studying anomalous transport
phenomena near critical points. The possible extension to include
multiplicative noises is also discussed.

##  I Introduction

In the phase diagram of quantum chromodynamics (QCD), the critical point (CP)
is the end point of the first-order phase transition line in the temperature-
baryon chemical potential plane. Its property is one of the most fascinating
features of the QCD phase diagram and draws the extensive attention of the
community of heavy-ion collisions all the time. Since the concept of the QCD
critical point was first put forward, constant theoretical efforts have been
put into the related studies. Nowadays its existence is predicted by various
QCD effective models and suggested by lattice QCD simulations [1, 2, 3, 4, 5,
6, 7, 8, 9, 10]. Based on these existing researches, the Beam Energy Scan
(BES) program was launched by Relativistic Heavy Ion Collider (RHIC) at
Brookhaven National Laboratory with one physical goal of searching for this
conjectured QCD critical point. Recent relativistic heavy-ion collision
experiments at RHIC [11], Facility for Antiproton and Ion Research (FAIR) [12,
13], and Nuclotron-base Ion Collider fAcility (NICA) [14, 15] are devoted to
finding the QCD critical point and have seen much experimental progress [16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. At
this point very large fluctuation will occur, which could induce a critical
phenomenon of divergent bulk viscosity [35] as a potential signal of the
existence of the critical point in experiments. However, as pointed out in
[36], the universal mechanism causing the critical divergence of the transport
coefficients has not been elucidated clearly. In fact, it is the nonlinear
fluctuations between the fluctuations of macroscopic variables that lead to
the divergence: the transport coefficients thus can be separated into two
parts, the rapidly relaxing part of microscopic time scale decided by a
microscopic theory, such as the Boltzmann equation [37] and the long-time tail
of macroscopic time scale. Only by taking into account nonlinear interactions
between macroscopic modes can the divergent behavior of the transport
coefficients be fully understood [38, 39], especially when the thermal
fluctuations of long-wavelength modes are necessarily large near the critical
point, which motivates the present study on the relativistic fluctuating
hydrodynamics.

Another important motivation is the introduction of the method of projection
operator for studying the dynamic evolution of some relativistic systems. The
projection operator, initially proposed by Mori and Zwanzig [40, 41], has long
been a powerful tool for the extraction of slow dynamics. In many cases of our
interest, especially when the behavior of a system over long distances and
time scales is concerned, there is a clear separation in the time scales
between the relevant slow variables and the vast number of fast microscopic
degrees of freedom. The distinction of what is slow or fast is not that
arbitrary but dictated by the fundamental physical principles. For instance,
the conservation laws give rise to slow variables. According to the Noether
theorem, every continuous symmetry corresponds to a conserved observable that
generates the symmetry. The dynamics of the local density associated with the
conserved generator exhibits a slow behavior, as is manifestly shown in
hydrodynamics [42]. Another main source of slow variables concerns the phase
transition. As one passes into the ordered state in a phase transition, the
slow variables associated with the breaking of a continuous symmetry appear
known as Nambu–Goldstone modes [43, 44]. For instance, the structure of the
fluctuation function of the transverse staggered fields in isotropic
antiferromagnet is changed qualitatively by the broken symmetry and looks much
like the hydrodynamic correlation functions in a fluid, indicating that the
Nambu–Goldstone modes associated with the spontaneous symmetry breaking have
the dynamical effect of simulating a conservation law [45]. A third example of
a slow variable is the order parameter near a second-order phase transition,
during which the system is experiencing a phenomenon called critical slowing
down [46]. Besides, the time scale separation, as observed in Brownian motion
[47, 48], typically results from the significant difference in masses between
the pollen particle and the background fluid particles, which is also a very
classical example. Note that the formalism of the projection operator has been
widely applied to study slow dynamics listed above in the context of condensed
matter physics and modern statistical physics, see [41, 49, 45, 50, 51] and
the reference therein for a comprehensive overview. However, there are fewer
related works in the context of relativistic statistics [52, 53, 54].
Therefore, it is worthwhile to give discussions on the application of the
projection operator techniques to relativistic situations such as heavy-ion
collisions, cosmology, and astrophysics, as we will devote our efforts to
below.

In this paper, we utilize the nonlinear projection operator method developed
by Morozov [55, 56] to derive the nonlinear evolution equations of
relativistic fluctuating hydrodynamics in a fully quantum manner. Before
passing into the next stage, it is necessary to pause and give comments on the
adopted formalism. Mori’s original projection operator is inherently linear in
the relevant variables and when used to describe a fluid system, it is only
applicable to systems in the linear regime [53]. To circumvent this problem,
new types of projection operators were proposed with notable contributions in
this regard [38, 57, 58, 59]. The Dirac delta functions with macrovariables as
arguments are key elements for technically resolving these issues. As will be
demonstrated in the main text, this elegant approach not only allows for the
incorporation of macrovariable dynamics into a single equation but also
facilitates the derivation of both the generalized Fokker-Planck equation and
the generalized Langevin equation simultaneously, while maintaining exactness.
The nonlinear projection operator proposed by Morozov is a quantum treatment
including the non-commutativity of operators. In this case, the Dirac delta
function is enhanced to the Dirac delta operator according to the Weyl
correspondence principle [60, 61] in accordance with its classical form. In
the meanwhile, the c-number symbols of operators is used to characterize the
dynamical variables. We leave the details about the Weyl symbol description to
the following sections.

This paper is organized as follows. In Sec. II, we present a brief review of
Morozov’s quantum nonlinear projection operator. In Sec. III and subsequent
IV, the generalized Fokker-Planck equation and the generalized Langevin
equation are derived based on the formalism of nonlinear projection operator.
After finishing these, Sec. V serves as a quantum derivation of fluctuating
hydrodynamics with Gaussian noise. The possible extension to include
multiplicative noises in the present framework is discussed in Sec. VI.
Finally, we give a summary and outlook in Sec. VII.

Natural units $\hbar=k_{B}=c=1$ are used. The metric tensor here is given by
$g^{\mu\nu}=\operatorname{diag}(1,-1,-1,-1)$, while $\Delta^{\mu\nu}\equiv
g^{\mu\nu}-u^{\mu}u^{\nu}$ is the projection tensor orthogonal to the four-
vector fluid velocity $u^{\mu}$. In addition, we employ the symmetric
shorthand notations:

$\displaystyle X^{(\mu\nu)}$ $\displaystyle\equiv$
$\displaystyle(X^{\mu\nu}+X^{\nu\mu})/2,$ (1) $\displaystyle
X^{\langle\mu\nu\rangle}$ $\displaystyle\equiv$
$\displaystyle\bigg{(}\frac{\Delta^{\mu}_{\alpha}\Delta^{\nu}_{\beta}+\Delta^{\nu}_{\alpha}\Delta^{\mu}_{\beta}}{2}-\frac{\Delta^{\mu\nu}\Delta_{\alpha\beta}}{3}\bigg{)}X^{\alpha\beta}.$
(2)

##  II Nonlinear quantum projection operator method

In this section, we give a short review on Morozov’s nonlinear projection
operator method [55], which can be used elegantly to extract the slow dynamics
from microscopic Hamiltonian dynamics. According to quantum mechanics, the
evolution of an operator at time $t$,
$\hat{B}(t)=e^{i\hat{H}t}\hat{B}(0)e^{-i\hat{H}t}$, is ruled by the Heisenberg
equation,

$\partial_{t}\hat{B}(t)=i[\hat{H},\hat{B}(t)]\equiv iL\hat{B}(t),$ (3)

where the Hamiltonian $\hat{H}$ encodes all the information of microscopic
dynamics. Formally, the Liouville operator $L$ is defined and
$\hat{B}(t)=e^{iLt}\hat{B}(0)$. For the need of discussing a many-body
statistical system, a density operator $\hat{\rho}(\Gamma,t)$ describing the
phase space $\Gamma$ distribution of a given ensemble system is indispensable.
The Liouville equation governing the evolution of the density operator is
written as

$\displaystyle\partial_{t}\hat{\rho}(t)=-iL\hat{\rho}(t)$ (4)

where the dependence on $\Gamma$ has been omitted. Note as an aside,
$\hat{\rho}(t)$ is sometimes called the nonequilibrium statistical operator.
Formally, $\hat{\rho}(t)$ takes the form
$\hat{\rho}(t)=e^{-iLt}\hat{\rho}(0)$. The average of $\hat{B}$ over the
density operator $\rho(t)$ is defined as

$\displaystyle\langle\hat{B}\rangle=\mathop{\mathrm{Tr}}(\hat{\rho}(t)\hat{B})=\mathop{\mathrm{Tr}}(\hat{\rho}\hat{B}(t)\,)$
(5)

and the second equality indicates that the same expectation value is given in
both Schrödinger picture and Heisenberg picture. Eq.(5) also means that $iL$
is an anti-self-adjoint operator in the sense of

$\displaystyle\mathop{\mathrm{Tr}}\big{(}\hat{B}iL\hat{C}\big{)}=-\mathop{\mathrm{Tr}}\big{(}\hat{C}iL\hat{B}\big{)}.$
(6)

In the remainder of this text, we choose to work in Heisenberg picture,
contrary to the one used in [55].

The physical contemplation must be initiated right from the outset, which lies
in a sensible consideration about how to choose the basis variables that are
projected onto. A first key point is the identification of the slow dynamic
variables. Considering that the conservation laws can be naturally made a
guiding principle, a set of coarse-grained collective variables for describing
slowly-varying macroscopic processes in the system are chosen as the basis
vectors
$\\{\hat{A}(t,\bm{x})\\}=\\{\hat{A}_{1},\hat{A}_{2},...,\hat{A}_{N}\\}$,
namely, all the local densities of correspondent conserved quantities are
included in the list, which are closely related to the slow modes as a result
of their conserved property. As a supplement, these properly chosen collective
variables are all Hermitian.

Mori’s linear projection operator formalism possess inherent shortcomings and
has no access to the nonlinear interactions between the slow modes. To
overcome it, a Dirac delta operator is introduced

$\displaystyle\hat{f}(a)=\delta(\hat{A}-a)=\frac{1}{(2\pi)^{N}}\int
dx\exp(ix\cdot(\hat{A}-a))$ (7)

which is obtained according to Weyl correspondence rule [55, 60, 61]
consistent with its corresponding classical form

$\displaystyle f(a)=\delta(A-a)=\Pi_{n=1}^{N}\delta(A_{n}-a_{n})$ (8)

where $N$ is the number of basis variables. With the Dirac delta operator at
hand, an arbitrary operator can be represented in the form of the following
identity with the help of the Weyl symbol,

$\displaystyle g(\hat{A})=\int dag(a)\delta(\hat{A}-a)$ (9)

where $g(a)$ is the Weyl symbol of $g(\hat{A})$. A key observation should be
notice that $\delta(\hat{A}-a)$ can generate all nonlinear couplings among the
variables $\hat{A}_{i}$ effortlessly as long as $g(a)$ is nonlinear function
of $a_{i}$.

In addition, we introduce

$\displaystyle
f(a,t)\equiv\mathop{\mathrm{Tr}}\big{(}\hat{\rho}(t)\hat{f}(a)\big{)},$ (10)

which is the macroscopic probability density function in $a$ space. Once it is
integrated with $\hat{f}(a)$ over $a$, the distribution function of $\hat{A}$
will be obtained. For pedagogical reason, it is more convenient to rewrite
$f(a,t)$ in a classical form

$\displaystyle f(a,t)\equiv\int d\Gamma\delta(A(\Gamma)-a)\rho(\Gamma,t)$ (11)

where $\Gamma$ denotes the phase space variables and the integration covers
the phase space points $\Gamma$ lying on the compatible hypersurface with a
constraint $A(\Gamma)=a$. One can also transfer to Heisenberg picture without
doubts

$\displaystyle
f(a,t)=\mathop{\mathrm{Tr}}\big{(}\hat{\rho}\hat{f}(a,t)\big{)}$ (12)

where $\hat{f}(a,t)=e^{iLt}\hat{f}(a,0)=e^{iHt}\hat{f}(a,0)e^{-iHt}$ is an
operator in Heisenberg picture.

It is useful to recall the retarded Liouville equation taking the form

$\displaystyle\big{(}\frac{\partial}{\partial
t}+iL\big{)}\hat{\rho}(t)=-\epsilon(\hat{\rho}(t)-\hat{\rho}_{q}(t))$ (13)

where the quasiequilibrium statistical operator $\hat{\rho}_{q}$ is introduced
and $\epsilon\rightarrow 0$ after taking thermodynamic limit [62, 63].
Compared to Eq.(4), the presence of the infinitesimal source $\epsilon$ in the
right-hand side (rhs) of this equation breaks the time reversal symmetry of
Liouville equation and selects causal and retarded solutions. The form of the
quasi-equilibrium statistical operator $\hat{\rho}_{q}(t)$ is determined by
the condition of maximum informational entropy and will be specified soon. In
this script, we assume that $\hat{\rho}_{q}(t)$ and $\rho(t)$ are expressed
through the distribution function of basic variables $f(a,t)$. Furthermore,
the statistical operator $\hat{\rho}_{q}$ also needs to satisfy the self-
consistency condition

$\displaystyle
f(a,t)=\mathop{\mathrm{Tr}}\big{(}\hat{\rho}\hat{f}(a,t)\big{)}=\mathop{\mathrm{Tr}}\big{(}\hat{\rho}_{q}\hat{f}(a,t)\big{)},$
(14)

which is equivalent to the Landau matching conditions widely used in
relativistic kinetic theory for defining the quasiequilibrium distribution
function [42, 37]. Note we have shifted to the Heisenberg picture in Eq.(14).

Using the Weyl correspondence rule, $\hat{\rho}_{q}$ can be put in the
following form

$\displaystyle\hat{\rho}_{q}=\int
da^{\prime}G(a^{\prime})\hat{f}(a^{\prime}).$ (15)

The task in this stage turns into the determination of the Weyl symbol of
$\hat{\rho}_{q}$ based on the self-consistency condition Eq.(14). To make it,
we’ll multiply Eq.(15) by $\hat{f}(a)$ and use Eq.(14). As a result, we obtain

$\displaystyle f(a,0)=\int da^{\prime}W(a,a^{\prime})G(a^{\prime})$ (16)

with the definition

$\displaystyle
W(a,a^{\prime})=\mathop{\mathrm{Tr}}(\hat{f}(a)\hat{f}(a^{\prime})\,)$ (17)

possessing the following properties

$\displaystyle W(a)=\int
da^{\prime}W(a,a^{\prime})=\mathop{\mathrm{Tr}}(\hat{f}(a)\,),$ $\displaystyle
W(a^{\prime})=\int
daW(a,a^{\prime})=\mathop{\mathrm{Tr}}(\hat{f}(a^{\prime})\,)$ (18)

where $W(a)$ denotes a quantum-mechanical generalization of the classical
density of states with given values of the coarse-grained variables $a$.

From Eq.(16), we get

$\displaystyle G(a)=\int da^{\prime}W_{-1}(a,a^{\prime})f(a^{\prime},0)$ (19)

where the following identity is satisfied

$\displaystyle\int
da^{\prime\prime}W(a,a^{\prime\prime})W_{-1}(a^{\prime\prime},a^{\prime})=\delta(a-a^{\prime}).$
(20)

After substituting Eq.(19) into (15), the quasiequilibrium statistical
operator is brought to the form

$\displaystyle\hat{\rho}_{q}=\int
dada^{\prime}\hat{f}(a)W_{-1}(a,a^{\prime})f(a^{\prime},0).$ (21)

Combining Eqs.(II) and (20), it is not hard to show that $W(a,a^{\prime})$ and
$W_{-1}(a,a^{\prime})$ have singular parts that can be conveniently isolated.
To that end, we assume

$\displaystyle W(a,a^{\prime})$
$\displaystyle=W(a)\big{(}\delta(a-a^{\prime})-R(a,a^{\prime})\big{)}$ (22)
$\displaystyle W_{-1}(a,a^{\prime})$
$\displaystyle=W^{-1}(a^{\prime})\big{(}\delta(a-a^{\prime})+r(a,a^{\prime})\big{)}$
(23)

then the following relations should be satisfied by analyzing Eq.(II)

$\displaystyle\int daW(a)R(a,a^{\prime})$ $\displaystyle=\int
da^{\prime}R(a,a^{\prime})=0,$ (24)

and the substitution of Eqs.(22) and (23) into Eq.(20) leads to

$\displaystyle r(a,a^{\prime})=R(a,a^{\prime})+\int
da^{\prime\prime}R(a,a^{\prime\prime})r(a^{\prime\prime},a^{\prime}).$ (25)

Given $R(a,a^{\prime})$, $r(a,a^{\prime})$ is obtained by iteration. The use
of Eq.(24) gives

$\displaystyle\int daW(a)r(a,a^{\prime})$ $\displaystyle=\int
da^{\prime}r(a,a^{\prime})=0,$ (26)

where the last equality sign follows from the iteration of the last equality
sign of Eq.(24). The regular terms appearing in Eqs.(22) and (23) originate
from the non-commutativity of quantum operators.

If focusing on the derivation in [55], one will find that Eq.(21) will play a
central role in getting the definition of the projection operator. The author
utilizes this to calculate the evolution function of $\hat{\rho}_{q}(t)$ (in
Schrödinger picture) and the projection operator is naturally introduced to
rewrite the resulting equation in a compact form. Following [55], the
projection operator $P$ is defined as

$\displaystyle P\hat{B}=\int
dada^{\prime}\hat{f}(a)W_{-1}(a,a^{\prime})\mathop{\mathrm{Tr}}(\hat{B}\hat{f}(a^{\prime})\,),$
(27)

as can be seen clearly, this definition is independent of which picture one
works in. If $W_{-1}(a,a^{\prime})\sim\delta(a-a^{\prime})$ is local, then $P$
degenerates into a similar form to classical delta projection operator [64].
This matching is not coincidental: if one assumes the initial state is the
microcanonical distribution $f(a^{\prime},0)=\delta(a^{\prime}-a_{0})$, then
in the local approximation

$\displaystyle\hat{\rho}_{q}=\delta(\hat{A}-a_{0})/W(a_{0})$ (28)

which is nothing but Eq.(12) of [64].

It is easy to verify $P$ is idempotent, i.e., $P^{2}=1$. Moreover, the
projection operator satisfies

$\displaystyle P\hat{f}(a)=\hat{f}(a),$ (29) $\displaystyle
P\hat{\rho}=\hat{\rho}_{q},$ (30)
$\displaystyle\mathop{\mathrm{Tr}}(\hat{B}P\hat{C})=\mathop{\mathrm{Tr}}(\hat{C}P\hat{B}),$
(31) $\displaystyle PG(\hat{A})=G(\hat{A}),$ (32)

where the first three equalities are almost self-evident, $G(\hat{A})$ is a
nonlinear function of $\hat{A}$ and a proof of the last one is given in
Eq.(106). The fourth equation Eq.(32) reveals the nonlinearity of $P$, which
generalizes $P\hat{A}=\hat{A}$ to $PG(\hat{A})=G(\hat{A})$ such that the
neglected nonlinear components in linear projection operator is retrieved.

##  III Generalized Fokker-Planck equation

This section is devoted to the derivation of generalized Fokker-Planck
equation. We start from

$\displaystyle\frac{\partial\hat{f}(a,t)}{\partial t}=e^{iLt}iL\hat{f}(a).$
(33)

Recall the Dyson-Duhamel identity [65]

$\displaystyle
e^{iLt}=e^{iLt}P+\int_{0}^{t}due^{iLu}PiL(1-P)e^{iL(1-P)(t-u)}+(1-P)e^{iL(1-P)t},$
(34)

this identity can be derived via Laplace transformation just as we will show
below. First, the following decomposition is employed:

$\displaystyle\partial_{t}e^{iLt}$
$\displaystyle=e^{iLt}iL=e^{iLt}PiL+e^{iLt}(1-P)iL.$ (35)

Next, the Laplace transform of $\exp({iLt})$ lead us to,

$\int_{0}^{\infty}dte^{-zt}e^{iLt}=\frac{1}{z-iL}.$ (36)

Then, a decomposition of Eq. (36) into

$\displaystyle\frac{1}{z-iL}$
$\displaystyle=\frac{1}{z-iL}(z-(1-P)iL)\frac{1}{z-(1-P)iL}$
$\displaystyle=\frac{1}{z-iL}(z-iL+PiL)\frac{1}{z-(1-P)iL}$
$\displaystyle=\frac{1}{z-(1-P)iL}+\frac{1}{z-iL}PiL\frac{1}{z-(1-P)iL}$ (37)

is utilized. After performing the inverse Laplace transform, a crucial
identity is reached,

$\displaystyle e^{iLt}=e^{(1-P)iLt}+\int_{0}^{t}dse^{iL(t-s)}PiLe^{(1-P)iLs},$
(38)

which follows from the derivation in [53]. Noticing that

$\displaystyle(1-P)e^{iL(1-P)t}=e^{(1-P)iLt}(1-P),$ (39)

the multiplication of eq.(38) by $(1-P)$ from the right will give us the
desired Dyson-Duhamel identity. The proof has been completed.

It is time to apply Dyson-Duhamel identity to the rhs of Eq.(33)

$\displaystyle\frac{\partial\hat{f}(a,t)}{\partial
t}=e^{iLt}PiL\hat{f}(a)+\int_{0}^{t}due^{iLu}PiL(1-P)e^{iL(1-P)(t-u)}iL\hat{f}(a)+(1-P)e^{iL(1-P)t}iL\hat{f}(a).$
(40)

Then we introduce a useful definition

$\displaystyle iL\hat{f}(a)\equiv-\frac{\partial\hat{J}_{i}(a)}{\partial
a_{i}}$ (41)

where Einstein summation conventions are respected hereafter and

$\displaystyle\hat{J}(a)\equiv\frac{1}{(2\pi)^{N}}\int
dxe^{ix(\hat{A}-a)}\int_{0}^{1}d\tau e^{-i\tau x\hat{A}}iL\hat{A}e^{i\tau
x\hat{A}}.$ (42)

In deriving this equation, an important identity, i.e., Kubo identity [66]

$\displaystyle[e^{\hat{B}},\hat{H}]=e^{\hat{B}}\int_{0}^{1}d\tau
e^{-\tau\hat{B}}[\hat{B},\hat{H}]e^{\tau\hat{B}}$ (43)

has been utilized.

The first term in Eq.(40) is

$\displaystyle e^{iLt}PiL\hat{f}(a)$ $\displaystyle=-\int
da^{\prime\prime}da^{\prime}W_{-1}(a^{\prime\prime},a^{\prime})\mathop{\mathrm{Tr}}(\frac{\partial\hat{J}_{i}(a)}{\partial
a_{i}}\hat{f}(a^{\prime})\,)\hat{f}(a^{\prime\prime},t)$
$\displaystyle=-\frac{\partial}{\partial a_{i}}\int
da^{\prime}v_{i}(a,a^{\prime})\hat{f}(a^{\prime},t)$ (44)

with the nonlocal streaming velocity

$\displaystyle v_{i}(a,a^{\prime})$ $\displaystyle\equiv\int
da^{\prime\prime}W_{-1}(a^{\prime},a^{\prime\prime})\mathop{\mathrm{Tr}}(\hat{J}_{i}(a)\hat{f}(a^{\prime\prime})\,).$
(45)

Define

$\displaystyle\hat{X}_{i}(a)=(1-P)\hat{J}_{i}(a),$ (46)

then the noise term is

$\displaystyle(1-P)e^{iL(1-P)t}iL\hat{f}(a)=-\frac{\partial}{\partial
a_{i}}\hat{X}_{i}(a,t).$ (47)

The remaining diffusion term is

$\displaystyle
PiL(1-P)e^{iL(1-P)(t-u)}iL\hat{f}(a)=PiLe^{i(1-P)L(t-u)}(1-P)iL\hat{f}(a)$
$\displaystyle=\frac{\partial}{\partial a_{i}}\int
da^{\prime}K_{ij}(a,a^{\prime},t-u)\frac{\partial}{\partial
a^{\prime}_{j}}\int
da^{\prime\prime}\hat{f}(a^{\prime\prime})W_{-1}(a^{\prime\prime},a^{\prime}).$
(48)

where
$K_{ij}(a,a^{\prime},t)\equiv\mathop{\mathrm{Tr}}(\hat{X}_{i}(a,t)\hat{X}_{j}(a^{\prime})\,)$
and $\hat{X}(a,t)\equiv e^{(1-P)iLt}\hat{X}(a)$. The relation between the
diffusion kernel $K$ and noise function $X$ is known as a kind of generalized
fluctuation-dissipation theorem, which is of great theoretical significance in
statistical physics. In the process of deriving the diffusion kernel $K_{ij}$,
the properties of $P$ listed in Eqs.(29) to (31) are frequently used.

Putting them together, we are reaching the quantum Fokker-Planck equation in
operator form

$\displaystyle\frac{\partial\hat{f}(a,t)}{\partial t}=$
$\displaystyle-\frac{\partial}{\partial a_{i}}\int
da^{\prime}v_{i}(a,a^{\prime})\hat{f}(a^{\prime},t)+\int_{0}^{t}du\frac{\partial}{\partial
a_{i}}\int da^{\prime}K_{ij}(a,a^{\prime},t-u)\frac{\partial}{\partial
a^{\prime}_{j}}\int
da^{\prime\prime}\hat{f}(a^{\prime\prime},u)W_{-1}(a^{\prime\prime},a^{\prime})$
$\displaystyle-\frac{\partial}{\partial a_{i}}\hat{X}_{i}(a,t),$ (49)

which is consistent with the Fokker-Planck equation obtained in a fully
classical way [38]. The fluctuating force $\frac{\partial}{\partial
a_{i}}\hat{X}_{i}(a,t)$ is uncorrelated with an arbitrary function
$G(\hat{A})$ of $\hat{A}$,

$\displaystyle\mathop{\mathrm{Tr}}(G(\hat{A})\hat{X}(a,t))=0.$ (50)

A comment on these three terms in the rhs of Eq.(III) will be left where the
nonlinear Langevin equations are derived.

Then we calculate the average of Eq.(III) with respect to $\hat{\rho}$,

$\displaystyle\frac{\partial f(a,t)}{\partial t}=$
$\displaystyle-\frac{\partial}{\partial a_{i}}\int
da^{\prime}v_{i}(a,a^{\prime})f(a^{\prime},t)+\int_{0}^{t}du\frac{\partial}{\partial
a_{i}}\int da^{\prime}K_{ij}(a,a^{\prime},t-u)\frac{\partial}{\partial
a^{\prime}_{j}}\int
da^{\prime\prime}f(a^{\prime\prime},u)W_{-1}(a^{\prime\prime},a^{\prime})$
(51)

which is in agreement with the Fokker-Planck equation obtained in [55] but in
a rather different way (directly solve $\hat{\rho}(t)$ from retarded Liouville
equation and next construct $f(a,t)$ based on it). Identical results are given
in two distinct quantum pictures, as it should be. The similar Fokker-Planck
equation has been extensively studied with different methods [64]. There are
also a large number of other works devoting to this, but due to space
limitations, the reference here may not cover them. In this script, we name
the equation as the Fokker-Planck equation in function form. A quick
observation can be made that the noise term does not survive the expectation
because

$\displaystyle\mathop{\mathrm{Tr}}[\hat{\rho}(1-P)e^{iL(1-P)t}iL\hat{f}(a)]=$
$\displaystyle\mathop{\mathrm{Tr}}[\hat{\rho}_{q}(1-P)e^{iL(1-P)t}iL\hat{f}(a)]=\mathop{\mathrm{Tr}}[((1-P)\hat{\rho}_{q})e^{iL(1-P)t}iL\hat{f}(a)]=0.$
(52)

Recall that in Heisenberg picture $\hat{\rho}=\hat{\rho}_{q}$ if the system is
initially in equilibrium, as is implicitly prescribed consistent with the
retarded Liouville equation, see Eq.(4.3) of ref [55]. Due to this identical
setting, Eq.(51) matches the resulting equation derived in [55].

Let us pause for a while to find a very useful relation. The trace over
$\hat{f}(a,t)$ should be independent of time

$\displaystyle\mathop{\mathrm{Tr}}(\hat{f}(a,t)\,)=\mathop{\mathrm{Tr}}(\hat{f}(a)\,)=W(a)$
(53)

owing to the cyclical symmetry of trace. One should not confuse
$f(a,t)=\mathop{\mathrm{Tr}}(\hat{\rho}\hat{f}(a,t))$ with
$W(a)=\mathop{\mathrm{Tr}}(\hat{f}(a,t)\,)$. Then the trace over the whole
Fokker-Planck equation (III) is

$\displaystyle\frac{\partial W(a)}{\partial t}=$
$\displaystyle\,0=-\frac{\partial}{\partial a_{i}}\int
da^{\prime}v_{i}(a,a^{\prime})W(a^{\prime})+\int_{0}^{t}du\frac{\partial}{\partial
a_{i}}\int da^{\prime}K_{ij}(a,a^{\prime},t-u)\frac{\partial}{\partial
a^{\prime}_{j}}\int
da^{\prime\prime}W(a^{\prime\prime})W_{-1}(a^{\prime\prime},a^{\prime})$
$\displaystyle-\frac{\partial}{\partial
a_{i}}\mathop{\mathrm{Tr}}(\hat{X}_{i}(a,t)\,).$ (54)

Now the remaining task awaiting us is to prove the second equality sign to see
whether it holds naturally or offers extra constraints.

The first term in the rhs is nothing but the divergence condition [67],

$\displaystyle\frac{\partial}{\partial a_{i}}\int
da^{\prime}v_{i}(a,a^{\prime})W(a^{\prime})=\frac{\partial}{\partial
a_{i}}\int da^{\prime}\int
da^{\prime\prime}\mathop{\mathrm{Tr}}(\hat{J}_{i}(a)\hat{f}(a^{\prime\prime})\,)W_{-1}(a^{\prime},a^{\prime\prime})W(a^{\prime})$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial a_{i}}\int
da^{\prime}\mathop{\mathrm{Tr}}(\hat{J}_{i}(a)\hat{f}(a^{\prime})\,)=\frac{\partial}{\partial
a_{i}}\mathop{\mathrm{Tr}}(\hat{J}_{i}(a)\,)=-\mathop{\mathrm{Tr}}(iL\hat{f}(a)\,)=i\mathop{\mathrm{Tr}}([\hat{f}(a),\hat{H}])=0.$
(55)

Because this one holds unconditionally, we should promote it to the divergence
theorem instead of treating it as a condition. The same conclusion is also
drawn in the context of classical statistical [45]. The second one vanishes
for

$\displaystyle\int
da^{\prime\prime}W(a^{\prime\prime})W_{-1}(a^{\prime\prime},a^{\prime})=1+W^{-1}(a^{\prime})\int
da^{\prime\prime}W(a^{\prime\prime})r(a^{\prime\prime},a^{\prime})=1$ (56)

is a constant, where the condition Eq.(26) is used.

As for the last noise term, it also vanishes as exhibited in appendix.B. To
conclude, $\frac{\partial W(a)}{\partial t}=0$ does not produce extra
constraints. These extra constraints are very key in the construction of
fluctuating nonlinear hydrodynamics for normal fluids in a classical context
[67]. We shall comment on this when talking about multiplicative noises.

##  IV Generalized Langevin equation

In this section, we derive the generalized Langevin equation, which can be
derived from the generalized Fokker-Planck equation. Multiply Eq.(III) with
$a$ and then integrate $a$, and we will find

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t)}{\partial t}=$
$\displaystyle-\int daa_{\alpha}\frac{\partial}{\partial a_{i}}\int
da^{\prime}v_{i}(a,a^{\prime})\hat{f}(a^{\prime},t)+\int
daa_{\alpha}\int_{0}^{t}du\frac{\partial}{\partial a_{i}}\int
da^{\prime}K_{ij}(a,a^{\prime},t-u)$
$\displaystyle\times\frac{\partial}{\partial a^{\prime}_{j}}\int
da^{\prime\prime}\hat{f}(a^{\prime\prime},u)W_{-1}(a^{\prime\prime},a^{\prime})-\int
daa_{\alpha}\frac{\partial}{\partial a_{i}}\hat{X}_{i}(a,t)$ $\displaystyle=$
$\displaystyle\int da\int
da^{\prime}v_{\alpha}(a,a^{\prime})\hat{f}(a^{\prime},t)+\int_{0}^{t}du\int
da^{\prime}(\frac{\partial}{\partial a^{\prime}_{j}}K_{\alpha
j}(a^{\prime},t-u)\,)\int
da^{\prime\prime}\hat{f}(a^{\prime\prime},u)W_{-1}(a^{\prime\prime},a^{\prime})+\hat{R}_{\alpha}(t)$
(57)

where $\hat{R}_{k}(t)\equiv e^{(1-P)iLt}(1-P)iL\hat{A}_{k}$, and recalling
that

$\displaystyle\hat{X}_{i}(a)=(1-P)\hat{J}_{i}(a)=\frac{1}{(2\pi)^{N}}(1-P)\int
dx\int_{0}^{1}d\tau e^{ix(\hat{A}-a)}e^{-i\tau
x\hat{A}}(iL\hat{A}_{i})e^{i\tau x\hat{A}},$ (58)

then we have

$\displaystyle K_{ij}(a^{\prime},t)=$ $\displaystyle\int
daK_{ij}(a,a^{\prime},t)=\int
da\mathop{\mathrm{Tr}}(\hat{X}_{i}(a,t)\hat{X}_{j}(a^{\prime})\,)$
$\displaystyle=\frac{1}{(2\pi)^{2N}}\int da\int dxdx^{\prime}\int d\tau
d\tau^{\prime}\mathop{\mathrm{Tr}}(e^{i(1-P)Lt}(1-P)e^{ix(\hat{A}-a)}e^{-i\tau
x\hat{A}}(iL\hat{A}_{i})e^{i\tau x\hat{A}}\,\,$
$\displaystyle\times(1-P)e^{ix^{\prime}(\hat{A}-a^{\prime})}e^{-i\tau^{\prime}x^{\prime}\hat{A}}(iL\hat{A}_{j})e^{i\tau^{\prime}x^{\prime}\hat{A}}\,)$
$\displaystyle=\frac{1}{(2\pi)^{N}}\int
dx^{\prime}\int_{0}^{1}d\tau^{\prime}\mathop{\mathrm{Tr}}\big{(}(e^{i(1-P)Lt}(1-P)iL\hat{A}_{i})\,\,(1-P)e^{ix^{\prime}(\hat{A}-a^{\prime})}e^{-i\tau^{\prime}x^{\prime}\hat{A}}(iL\hat{A}_{j})e^{i\tau^{\prime}x^{\prime}\hat{A}}\,\big{)}$
$\displaystyle=\mathop{\mathrm{Tr}}(\hat{R}_{i}(t)\hat{X}_{j}(a^{\prime})).$
(59)

Because of Eq.(50), the Langevin fluctuating force is also statistically
independent of the collective variables, namely,

$\displaystyle\mathop{\mathrm{Tr}}(G(\hat{A})\hat{R}_{i}(a,t))=0.$ (60)

Eq.(IV) is the full generalized quantum Langevin equation without any
approximations. The first term in the rhs is the drift term characterized by a
nonlocal streaming velocity $v(a,a^{\prime})$, whose nonlocality stems from
the non-commutativity of quantum operators and thus is a quantum effect absent
in the classical Langevin equation. The second one is a diffusion term with
non-Markov memory effects. If looking at the specific form of the diffusion or
memory kernel $K_{ij}(a,t)$, one can also find that the quantum non-
commutativity is involved hindering the direct comparison with its classical
correspondence. The third one represents fast motion of noise, whose relation
with the diffusion kernel is so-called fluctuation-dissipation theorem but not
that evident at a simple glance. In order to correctly account for the
physical effects brought by the nonlinear fluctuations, the stochastic noise
terms, trigger by the hydrodynamic excitations generated by thermal
fluctuations within the fluid, are indispensable [68].

In order to compare with the existing classical results, we have to turn to
the local approximation, that is,

$\displaystyle v(a,a^{\prime})=v(a)\delta(a-a^{\prime}),\quad
W_{-1}(a^{\prime\prime},a^{\prime})=W^{-1}(a^{\prime})\delta(a^{\prime}-a^{\prime\prime})$
(61)

with a local averaged velocity

$\displaystyle v(a)=\mathop{\mathrm{Tr}}(\hat{f}(a)iL\hat{A})W^{-1}(a).$ (62)

In fact, a regular part can also be separated from $v(a,a^{\prime})$ in the
same fashion to $W(a,a^{\prime})$ and $W_{-1}(a,a^{\prime})$ of which the
details can be found in [55]. Then the full Langevin equation becomes

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t)}{\partial t}=$ $\displaystyle
v_{\alpha}(\hat{A},t)+\int_{0}^{t}du\int da^{\prime}(\frac{\partial}{\partial
a^{\prime}_{j}}K_{\alpha
j}(a^{\prime},t-u)\,)W^{-1}(a^{\prime})\hat{f}(a^{\prime},u)+\hat{R}_{\alpha}(t),$
(63)

which matches the result of [38]. Furthermore, if we introduce the definition
$\tilde{K}_{ij}(a,t)\equiv K_{ij}(a,t)/W(a)$, then the equation turns into

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t)}{\partial t}=$ $\displaystyle
v_{\alpha}(\hat{A},t)+\int_{0}^{t}du\int da\frac{\partial}{\partial
a_{j}}(\tilde{K}_{\alpha
j}(a,t-u)W(a)\,)W^{-1}(a)\hat{f}(a,u)+\hat{R}_{\alpha}(t)$ $\displaystyle=$
$\displaystyle v_{\alpha}(\hat{A},t)+\int_{0}^{t}du\int
da(\frac{\partial}{\partial a_{j}}\tilde{K}_{\alpha
j}(a,t-u)\,)\hat{f}(a,u)+\int_{0}^{t}du\int da\tilde{K}_{\alpha
j}(a,t-u)F_{j}(a)\hat{f}(a,u)+\hat{R}_{\alpha}(t)$ (64)

in good agreement with the form present in [64]. Here the conjugate variable
of $a_{k}$ is defined, in a conventional way, as the derivative of $\ln W(a)$
with respect to $a_{k}$

$\displaystyle F_{k}(a)=\frac{\partial}{\partial a_{k}}\ln W(a).$ (65)

Since $W(a)$ is the structure function providing a complete thermodynamic
description of a system in thermal equilibrium subject to constraints giving
rise to the specified values ”$\hat{A}=a$”, it can be seen as a partition
function of the microcanonical ensemble. More intuitively, Eq.(62) can be
rewritten as

$\displaystyle
v(a)=\mathop{\mathrm{Tr}}(\hat{f}(a)iL\hat{A})W^{-1}(a)=\mathop{\mathrm{Tr}}(\hat{\rho}(a)iL\hat{A})\,\,.$
(66)

Note $\mathop{\mathrm{Tr}}\hat{\rho}(a)=1$ and
$\mathop{\mathrm{Tr}}(\hat{\rho}(a)\hat{A}_{i}\,)=a_{i}$, then $\hat{\rho}(a)$
acts like a density operator in microcanonical ensemble with fixed values of
the coarse-grained variables. In the case where the coarse-grained variables
are the densities of thermodynamic quantities, one can take advantage of the
equivalence of statistical ensembles and replace the averaging in Eq.(62) with
a canonical or grand-canonical ensemble. Notably, the third term in the last
line of Eq.(IV) resembles the famous Onsager relation and can be viewed as a
generalized Onsager relation in a convolution form with memory effects.

The second term of Eq.(IV) is key and its presence indicates that the
diffusion kernel $\tilde{K}$ has in general nontrivial dependence on the
coarse-grained variables, which will be commented and stressed afterwards. In
this place, we want to show a simplified but very helpful form of the
generalized Langevin equation. By replacing the diffusion kernel with its
average

$\displaystyle\tilde{K}_{\alpha j}(a,t-u)\rightarrow\int
daW(a)\tilde{K}_{\alpha j}(a,t-u)=\int daK_{\alpha
j}(a,t-u)=\mathop{\mathrm{Tr}}(\hat{R}_{\alpha}(t-u)\hat{R}_{j}),$ (67)

where $W(a)$ is taken as a equilibrium distribution function of $a$ up to a
normalized factor

$\displaystyle\int daW(a)=\int
da\mathop{\mathrm{Tr}}(\hat{f}(a)\,)=\mathop{\mathrm{Tr}}(1)=\text{const}.$
(68)

This constant can be completely absorbed into the definition of other
coefficients so there is no need to worry about it. Such a procedure is also
implemented in [38] to seek a simpler equation. After finishing the
replacement, we will get

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t)}{\partial
t}=v_{\alpha}(\hat{A},t)+\int_{0}^{t}du\mathop{\mathrm{Tr}}(\hat{R}_{\alpha}(t-u)\hat{R}_{j})F_{j}(\hat{A}(u))+\hat{R}_{\alpha}(t)$
(69)

which recovers our familiar Langevin equation irrespective of the memory
structure. Moreover, the diffusion kernel is expressed as a two-point time
correlation function of the noise term.

As is seen from above, the local approximation is very crucial to the
simplification operations. At least, this is valid for classical systems.
Furthermore, the validity may also hold for some quantum systems. For example,
the contribution of the integral terms in the Fokker-Planck and Langevin
equations can be minor as long as the commutators of the basis operators are
parametrically small. In this case, the nonlocality represented by
$r(a,a^{\prime})$ are expressed in terms of a series expansion in the small
parameter, which is also the same case for the nonlocality of
$v(a,a^{\prime})$. On the other hand, if the non-equilibrium state of the
system is such that the distribution function $f(a,t)$ and the corresponding
density of states $W(a)$ change slowly when varying the distance
$|a-a^{\prime}|$ compared to the variation of $r(a,a^{\prime})$, then in the
first approximation, they can be considered constant when we calculate the
integral term in the Fokker-Planck and Langevin equations. In this case, the
contribution of regular parts becomes minor due to Eq.(26). Especially, in
some physical system such as the relativistic heavy ion collisions of our
interest, the created physical environment is a super hot medium. In this
case, the quantum effects are suppressed by thee scale of high temperature. In
this case, the quantum effects are often neglected and the local approximation
is acceptably useful.

##  V A quantum derivation of fluctuating hydrodynamics with Gaussian noise

###  V.1 Markov approximation

When casting the generalized Langevin equation into fluctuating hydrodynamics,
the Markov approximation is necessary to remove the delay effects appearing in
the integral term. Considering that the operator of projection $1-P$ in the
diffusion kernel $K(a,t)$ or $K(a,a^{\prime},t)$ eliminates the slow evolution
associated with the coarse-grained variables, the effect of time delays can be
neglected in the evolution equations when a clear separation of time scales
exists. This means that the time scales for an appreciable change of
$\hat{A}(t)$ are distinctly larger than the typical time scales involved in
the Langevin fluctuating force $\hat{R}_{k}(t)$. The simplified Langevin
equation is taken as the starting point

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t)}{\partial t}=$ $\displaystyle
v_{\alpha}(\hat{A}(t))+\int_{0}^{\infty}du\mathop{\mathrm{Tr}}(\hat{R}_{\alpha}(u)\hat{R}_{j})F_{j}(\hat{A}(t))+\hat{R}_{\alpha}(t)$
(70)

where the time hierarchy $t\gg u$ is assumed. A further transformation leads
to

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t)}{\partial t}=$ $\displaystyle
v_{\alpha}(\hat{A}(t))+\gamma_{\alpha j}F_{j}(\hat{A}(t))+\hat{R}_{\alpha}(t)$
(71)

where $\gamma$ is the bare kinetic coefficient defined as

$\displaystyle\gamma_{ij}\equiv\int_{0}^{\infty}du\mathop{\mathrm{Tr}}(\hat{R}_{i}(u)\hat{R}_{j}).$
(72)

The first two terms in the rhs of Eq.(70) generally contain nonlinear coupling
of $\hat{A}$ due to the presence of $v(\hat{A})$ and $F(\hat{A})$ while
$\hat{R}_{\alpha}(t)$ is a fluctuating force satisfying

$\displaystyle\mathop{\mathrm{Tr}}(\hat{R}_{i}(t)\hat{R}_{j}(t^{\prime})\,)=2\gamma_{ij}\delta(t-t^{\prime}).$
(73)

It can be shown in [38] that the nonlinear couplings contained in the drift
term will lead to the renormalization of the bare kinetic coefficient. Last
but not least, if choosing to write the $\bm{x}$ dependence explicitly, we get

$\displaystyle\frac{\partial\hat{A}_{\alpha}(t,\bm{x})}{\partial t}=$
$\displaystyle v_{\alpha}(\hat{A}(t,\bm{x}))+\frac{1}{V}\int
d\bm{x}^{\prime}\gamma_{\alpha
j}(\bm{x-x^{\prime}})F_{j}(\hat{A}(t,\bm{x}^{\prime}))+\hat{R}_{\alpha}(t,\bm{x})$
(74)

with

$\displaystyle\gamma_{ij}(\bm{x-x^{\prime}})\equiv\int_{0}^{\infty}du\mathop{\mathrm{Tr}}(\hat{R}_{i}(u,\bm{x})\hat{R}_{j}(0,\bm{x}^{\prime})\,).$
(75)

###  V.2 Relativistic hydrodynamics

The evolution of a relativistic fluid is governed by the continuity equations

$\displaystyle\partial_{\mu}T^{\mu\nu}=0\,,$ (76)
$\displaystyle\partial_{\mu}N^{\mu}=0\,,$ (77)

where $T^{\mu\nu},N^{\mu}$ are the energy momentum tensor and conserved
current respectively. In this subsection, the hat symbol of the operators is
temporarily neglected.

By utilizing symmetry analysis and the second law of thermodynamics,
$T^{\mu\nu}$ and $N^{\mu}$ can be conveniently written as in Landau frame

$\displaystyle T^{\mu\nu}=eu^{\mu}u^{\nu}-p\Delta^{\mu\nu}+\tau^{\mu\nu}\,,$
(78) $\displaystyle N^{\mu}=nu^{\mu}+j^{\mu}\,,$ (79)

with $e,n$ and $p$ being the energy density, the charge density and the
pressure. The metric tensor here is given by
$g^{\mu\nu}=\operatorname{diag}(1,-1,-1,-1)$, while $\Delta^{\mu\nu}\equiv
g^{\mu\nu}-u^{\mu}u^{\nu}$ is the projection tensor orthogonal to the four-
vector fluid velocity $u^{\mu}$. Also, the dissipative terms, $\pi$ and $j$,
are given by

$\displaystyle\tau^{\mu\nu}=-2\eta\nabla^{\langle\mu}u^{\nu\rangle}-\zeta\theta\Delta^{\mu\nu},$
(80) $\displaystyle j^{\mu}=-\frac{\sigma
T(e+p)}{n}\nabla^{\mu}\frac{\mu}{T},$ (81)

where $\mu$ is the chemical potential conjugate the conserved charge density
$n$ and $\nabla_{\mu}\equiv\Delta_{\mu\nu}\partial^{\nu}$ is the spatial
derivative. A shorthand for the thermodynamic force $\theta=\nabla\cdot u$
denotes the expansion rate. Three transport coefficients $\eta$, $\zeta$ and
$\kappa$ are the bare shear, bulk viscosities and the bare thermal
conductivity, respectively.

In the linear regime, hydrodynamic modes can be sought via the linearized
equations. A linear mode analysis is carried out on top of the background of
thermal equilibrium. In a relativistic fluid, the quiescent equilibrium system
is perturbed according to

$\displaystyle e(t,\bm{x})$ $\displaystyle=e_{0}+\delta e(t,\bm{x}),\quad
p(t,\bm{x})=p_{0}+\delta p(t,\bm{x}),$ $\displaystyle n(t,\bm{x})$
$\displaystyle=n_{0}+\delta n(t,\bm{x}),\quad
u^{\mu}(t,\bm{x})=(1,0)+(0,\delta v^{i}(t,\bm{x})),$ $\displaystyle
T(t,\bm{x})$ $\displaystyle=T_{0}+\delta
T(t,\bm{x}),\quad\frac{\mu}{T}=\frac{\mu_{0}}{T_{0}}+\delta(\frac{\mu}{T}).$
(82)

Due to $u\cdot u=1$, $(\delta v)^{2}\ll 1$, the terms of $O(\delta v^{2})$ is
neglected afterward. Then we have

$\displaystyle\frac{\partial\delta n}{\partial t}=$ $\displaystyle-
n_{0}\nabla\cdot\delta v-\frac{\sigma
T_{0}(e_{0}+p_{0})}{n_{0}}\nabla^{2}\delta\biggl{(}\frac{\mu}{T}\biggr{)},$
(83) $\displaystyle\frac{\partial\delta e}{\partial t}=$
$\displaystyle-(e_{0}+p_{0})\nabla\cdot\delta v,$ (84)
$\displaystyle\frac{\partial\pi_{i}}{\partial t}=$
$\displaystyle-\nabla_{i}\delta
p-(\zeta+\frac{1}{3}\eta)\nabla_{i}(\nabla\cdot\delta v)-\eta\nabla^{2}\delta
v_{i}.$ (85)

where $\pi^{i}\equiv\delta T^{0i}$.

###  V.3 Relativistic fluctuating hydrodynamics

The nonlinear interactions among these hydrodynamic modes is absent in
conventional hydrodynamics in the linear regime [53], more importantly,
without stochastic noises, conventional hydrodynamics ignores the inherent
hydrodynamic excitations triggered by thermal fluctuations. When the
thermodynamic fluctuations are large, the effects of nonlinear couplings
become significantly enhanced, which is the case in the critical regime. The
impacts of nonlinear fluctuations of the coarse-grained variables can be
properly accounted for in the framework of the Langevin equation. In this
subsection, we derive the relativistic fluctuating hydrodynamics from the
Langevin equation step by step.

Inspire by the relativistic hydrodynamic equations, the coarse-grained slow
variables set is taken to be $\\{\hat{A}(t,\bm{x})\\}=\\{\delta
e(t,\bm{x}),\delta n(t,\bm{x}),\pi_{i}(t,\bm{x})\\}$. First, we need to figure
out the streaming term

$\displaystyle
v(a)=\mathop{\mathrm{Tr}}(\hat{\rho}(a)iL\hat{A}(t,\bm{x}))=\mathop{\mathrm{Tr}}(\hat{\rho}(a)\partial_{t}\hat{A}(t,\bm{x})).$
(86)

If looking at the continuity equation,

$\displaystyle\partial_{t}\hat{A}(t,\bm{x})=\nabla\cdot\hat{J}_{A},$ (87)

the derivative of $\hat{A}$ with respect to $t$ gives exactly the related
reversible currents and dissipative currents vanishes when averaged over a
thermodynamic ensemble shown in the rhs of Eq.(86). Therefore, we can directly
read the streaming velocities,

$\displaystyle v_{n}=-\nabla\cdot(\hat{n}\delta\hat{v}),\quad
v_{e}=-\nabla\cdot\hat{\pi},\quad v_{\pi}=-\nabla_{i}\delta\hat{p},$ (88)

where the background is set to be rest and the nonlinear couplings between
fluctuations is manifestly shown in $v_{n}$ and $v_{\pi}$ (noticing that
$\hat{p}$ may contain nonlinear couplings in terms of these basis
fluctuations, see also [45]).

Next, recalling the given definition of the thermodynamic force,

$\displaystyle F_{k}(\hat{A})=\frac{\delta}{\delta\hat{A}_{k}}\ln W(\hat{A})$
(89)

we can translate $\delta(\mu/T)$ and $\delta v$ into
$\frac{1}{V}\frac{\delta\ln W(\hat{A})}{\delta\hat{n}}$ and
$\frac{T_{0}}{V}\frac{\delta\ln W(\hat{A})}{\delta\hat{\pi}}$ according to
Boltzmann relation $S\equiv\ln W(\hat{A})$ with $S,V$ being the total entropy
and volume. The factor $1/V$ comes without surprise because it counteracts the
same factor in Eq.(74).

Comparing the relativistic hydrodynamic equation Eqs.(76) to (81) with the
Langevin equation, the kinetic coefficients can also be read without efforts

$\displaystyle\gamma_{nn}$ $\displaystyle=-\frac{\sigma
T_{0}(e_{0}+P_{0})}{n_{0}}\nabla^{2},$ (90)
$\displaystyle\gamma_{\pi_{i}\pi_{j}}$
$\displaystyle=-\big{(}\,(\zeta+\frac{1}{3}\eta)\nabla_{i}\nabla_{j}+\eta\delta_{ij}\nabla^{2}\big{)},$
(91)

which demonstrates that

$\displaystyle\mathop{\mathrm{Tr}}(\hat{R}_{n}(t,\bm{x})\hat{R}_{n}(t^{\prime},\bm{x}^{\prime})\,)=-\frac{2\sigma
T_{0}(e_{0}+P_{0})}{n_{0}}\nabla^{2}\delta(t-t^{\prime})\delta(\bm{x}-\bm{x}^{\prime}),$
(92)
$\displaystyle\mathop{\mathrm{Tr}}(\hat{R}_{\pi_{i}}(t,\bm{x})\hat{R}_{\pi_{j}}(t^{\prime},\bm{x}^{\prime})\,)=-2\big{(}(\zeta+\frac{1}{3}\eta)\nabla_{i}\nabla_{j}+\eta\delta_{ij}\nabla^{2}\big{)}\delta(t-t^{\prime})\delta(\bm{x}-\bm{x}^{\prime}),$
(93)

where the spatial coordination $\bm{x}$ dependence is explicitly recovered.
Note through out the script, $\bm{x}$ dependence is frequently suppressed for
compactness sometimes.

The parameterization given in Eq.(90) and (91) is expected to successfully
reproduce the Ornstein–Zernike result [45, 36]. Attention is supposed to be
paid to the vanishing of the noise term $R_{e}$ and its corresponding
diffusion kernel. This is a nontrivial consistency cross check. The reason for
its vanishing is that its reversible current is also a local density of the
conserved quantity, i.e., $\delta\hat{T}^{0i}$, which is also in the list of
the chosen coarse-grained slow variables,

$\displaystyle\hat{R}_{e}\equiv(1-P)iL\hat{e}=(1-P)\partial_{i}\hat{\pi}_{i}=\partial_{i}(1-P)\hat{\pi}_{i}=0,$
(94)

where $P\hat{A}=\hat{A}$ is proven in the appendix.A.

The conclusion that the energy density does not dissipate depends on the
choice of the definition of the fluid velocity, or the choice of the set of
coarse-grained slow variables in the language of projection operator . If
working in Eckart frame, the charge density, rather than the energy density,
does not dissipate .

After all the dust has set down, we eventually obtain

$\displaystyle\frac{\partial\delta\hat{n}}{\partial t}=$
$\displaystyle-\nabla\cdot(\hat{n}\delta\hat{v})-\frac{\sigma
T_{0}(e_{0}+P_{0})}{n_{0}}\nabla^{2}\frac{\delta\ln
W(\hat{A})}{\delta\hat{n}}+\hat{R}_{n}(t,\bm{x}),$ (95)
$\displaystyle{\tiny}\frac{\partial\delta\hat{e}}{\partial t}=$
$\displaystyle-\nabla\cdot\hat{\pi},$ (96)
$\displaystyle\frac{\partial\hat{\pi}_{i}}{\partial t}=$
$\displaystyle-\nabla_{i}\delta\hat{p}-(\zeta+\frac{1}{3}\eta)\nabla_{i}(\nabla\cdot\frac{\delta\ln
W(\hat{A})}{\delta\hat{\pi}_{i}})-\eta\nabla^{2}\frac{\delta\ln
W(\hat{A})}{\delta\hat{\pi}_{i}}+\hat{R}_{\pi_{i}}(t,\bm{x}).$ (97)

Note our results here are almost identical to the ones derived in [36] in a
classical way except the streaming term of Eq.(97), where the author neglects
$\delta p$ and uses a potential condition to rewrite that streaming term.
Besides, $\nabla\cdot\pi$ is replaced by
$\nabla\cdot((\hat{e}+\hat{p})\delta\hat{v})$ therein. But we prefer
$\nabla\cdot\pi$ because the linearity in the basis fluctuations is
unambiguously revealed without causing extra confusion.

##  VI Possible extension to include multiplicative noises

In this section, we comment on the possible extension to include
multiplicative noises. The kinetic coefficients in Eqs.(95) to (97) is
independent of the fluctuating fields, therefore, the noise $\hat{R}$ is so-
called Guassian noise, which is not always the case. In view of the fact that
the kinetic coefficients characterize the intrinsic properties of the system
with given thermodynamic states $(e_{0},p_{0}\cdots)$. If the perturbation is
not too large to change the intrinsic properties, the description with
fluctuation-independent kinetic coefficients suffices, in line with the logic
hidden in the linear response theory. However, if the perturbation is strong
enough to fundamentally alter the intrinsic properties and the thermodynamic
state of the system, for instance, when the critical behavior in vicinity of
the QCD critical point is studied, the critical fluctuation is so large that
the Langevin equation with multiplicative noises and field-dependent kinetic
coefficients is essentially needed.

Such a fluctuating nonlinear hydrodynamics for normal fluids is given in the
context of classical nonrelativistic statistics [67]. The authors construct a
phenomenological model to achieve this. In their derivation, a very alike
formula $\frac{\partial W(a)}{\partial t}=0$ appears and provides three
nontrivial constraint conditions, among which there is a detailed balance
condition. The detailed balance condition plays an essential role in the
construction of the fluctuating nonlinear hydrodynamics, from which the
multiplicative noises and the diffusion kernel can be worked out just by
”reading” as is done in the last subsection.

Just as is commented around Eq.(III), there are no constraint conditions
coming from Eq.(III) according to our derivation step by step, which is also
confirmed by one of the authors of [67] in a later textbook [45] using a
classical treatment. Therefore, there is no shortcut to the possible extension
with multiplicative noises, we have to literally work out Eq.(IV), which is
left as a future work in progress.

##  VII Summary and Outlook

In this work, we present a systematic derivation of relativistic fluctuating
hydrodynamics employing the quantum nonlinear projection operator method.
Morozov’s nonlinear projection operator, as an extension of the well-
established linear Mori-Zwanzig projection operator, enables the consideration
of nonlinear interactions among macroscopic modes. This approach yields the
quantum Fokker-Planck and Langevin equations, which are key in investigating
anomalous transport phenomena in the vicinity of critical points. We then
apply these equations to the coarse-grained hydrodynamic variables, resulting
in the formulation of the relativistic fluctuating hydrodynamics with Gaussian
noise. Additionally, the potential extension of this framework to incorporate
the effects of multiplicative noises is also discussed.

There are some possible extensions which deserve to be focused on in future.
The first one, as noted above, is to bring the multiplicative noises in and
re-derive the relativistic fluctuating hydrodynamics. The multiplicative
noises shall play an important role in studying the critical behavior near the
QCD critical point and may shed light on the research for QCD critical point
in the experimental aspect. A comprehensive analysis by the dynamic
renormalization group based on the relativistic fluctuating hydrodynamics can
also be sought, which is a mature technique in studying the critical behavior
in condensed matter physics. In addition, the obtained generalized Fokker-
Planck and Langevin equations with nonlocal quantum effects can be helpful in
condense matter physics and low-temperature physics, where quantum effects are
non-negligible and of much research interest.

##  VIII Acknowledgments

J.H. is grateful to Yuki Minami for helpful correspondence and Shuzhe Shi for
helpful discussions during the early time of preparing this work.

##  Appendix A The proof of useful identities

In this appendix, we prove some useful identities associated with the
nonlinear projection operator. The first one is $P\hat{A}=\hat{A}$. We have

$\displaystyle P\hat{A}$ $\displaystyle=\int
dada^{\prime}\hat{f}(a)W_{-1}(a,a^{\prime})\mathop{\mathrm{Tr}}(\hat{A}\hat{f}(a^{\prime})\,)$
$\displaystyle=\int
dada^{\prime}a^{\prime}\hat{f}(a)W_{-1}(a,a^{\prime})W(a^{\prime})$
$\displaystyle=\int
dada^{\prime}a^{\prime}\hat{f}(a)(\delta(a-a^{\prime})+r(a,a^{\prime})\,)$
$\displaystyle=\hat{A}+\int da\hat{f}(a)\int
da^{\prime}a^{\prime}r(a,a^{\prime}),$ (98)

and

$\displaystyle\int da^{\prime}a^{\prime}W(a,a^{\prime})=\int
da^{\prime}a^{\prime}\mathop{\mathrm{Tr}}(\hat{f}(a)\hat{f}(a^{\prime})\,)=\mathop{\mathrm{Tr}}(\hat{f}(a)\hat{A}\,)=aW(a).$
(99)

On the other hand

$\displaystyle\int da^{\prime}a^{\prime}W(a,a^{\prime})=\int
da^{\prime}a^{\prime}W(a)(\delta(a-a^{\prime})+R(a,a^{\prime})\,)=aW(a)\rightarrow\int
da^{\prime}a^{\prime}R(a,a^{\prime})=0.$ (100)

Now we get

$\displaystyle a^{\prime}r(a,a^{\prime})=a^{\prime}R(a,a^{\prime})+\int
da^{\prime\prime}R(a,a^{\prime\prime})a^{\prime}r(a^{\prime\prime},a^{\prime}).$
(101)

By iterating Eq.(101), we reach

$\displaystyle\int da^{\prime}a^{\prime}r(a,a^{\prime})=0,$ (102)

which follows from the same reason for getting the second identity in Eq.(26).
Finally, we obtain $P\hat{A}=\hat{A}$.

This can be also derived in a simpler way, noticing that

$\displaystyle P\hat{f}(a)=\hat{f}(a),$ (103)

then a integral over $a$ constructed as

$\displaystyle\int daG(a)P\hat{f}(a)=\int daG(a)\hat{f}(a).$ (104)

Because $P$ only act upon the operator, it can be safely factorized out

$\displaystyle\int daG(a)P\hat{f}(a)=P\int daG(a)\hat{f}(a),$ (105)

this eventually gives us

$\displaystyle PG(\hat{A})=G(\hat{A}).$ (106)

Here $G(\hat{A})$ is a nonlinear function of $\hat{A}$.

##  Appendix B The vanishing of noise contribution

The trace of the noise term is

$\displaystyle\mathop{\mathrm{Tr}}[(1-P)e^{iL(1-P)t}iL\hat{f}(a)]=$
$\displaystyle\mathop{\mathrm{Tr}}[e^{iL(1-P)t}iL\hat{f}(a)]-\mathop{\mathrm{Tr}}[Pe^{iL(1-P)t}iL\hat{f}(a)]$
$\displaystyle=\mathop{\mathrm{Tr}}[\hat{Y}(a,t)]-\mathop{\mathrm{Tr}}[\int
da^{\prime\prime}da^{\prime}\hat{f}(a^{\prime\prime})W_{-1}(a^{\prime\prime},a^{\prime})\mathop{\mathrm{Tr}}(\hat{Y}(a,t)\hat{f}(a^{\prime})\,)]$
$\displaystyle=\mathop{\mathrm{Tr}}[\hat{Y}(a,t)]-\int
da^{\prime\prime}da^{\prime}W_{-1}(a^{\prime\prime},a^{\prime})\mathop{\mathrm{Tr}}\big{(}\hat{f}(a^{\prime\prime})\big{)}\mathop{\mathrm{Tr}}\big{(}\hat{Y}(a,t)\hat{f}(a^{\prime})\big{)}$
$\displaystyle=\mathop{\mathrm{Tr}}[\hat{Y}(a,t)]-\int
da^{\prime\prime}da^{\prime}W_{-1}(a^{\prime\prime},a^{\prime})W(a^{\prime\prime})\mathop{\mathrm{Tr}}\big{(}\hat{Y}(a,t)\hat{f}(a^{\prime})\big{)}$
$\displaystyle=\mathop{\mathrm{Tr}}[\hat{Y}(a,t)]-\mathop{\mathrm{Tr}}[\hat{Y}(a,t)]=0$
(107)

where the shorthand notation is $\hat{Y}(a,t)\equiv e^{iL(1-P)t}iL\hat{f}(a)$
and Eq.(26) has been used.

## References

  * Stephanov et al. [1998] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998), eprint hep-ph/9806219.
  * Stephanov et al. [1999] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999), eprint hep-ph/9903292.
  * Fodor and Katz [2004] Z. Fodor and S. D. Katz, JHEP 04, 050 (2004), eprint hep-lat/0402006.
  * Allton et al. [2005] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, and K. Redlich, Phys. Rev. D 71, 054508 (2005), eprint hep-lat/0501030.
  * Gavai and Gupta [2008] R. V. Gavai and S. Gupta, Phys. Rev. D 78, 114503 (2008), eprint 0806.2233.
  * de Forcrand and Philipsen [2008] P. de Forcrand and O. Philipsen, PoS LATTICE2008, 208 (2008), eprint 0811.3858.
  * Hatsuda et al. [2006] T. Hatsuda, M. Tachibana, N. Yamamoto, and G. Baym, Phys. Rev. Lett. 97, 122001 (2006), eprint hep-ph/0605018.
  * Braun-Munzinger et al. [2016] P. Braun-Munzinger, V. Koch, T. Schäfer, and J. Stachel, Phys. Rept. 621, 76 (2016), eprint 1510.00442.
  * Braun-Munzinger and Wambach [2009] P. Braun-Munzinger and J. Wambach, Rev. Mod. Phys. 81, 1031 (2009), eprint 0801.4256.
  * Fukushima and Sasaki [2013] K. Fukushima and C. Sasaki, Prog. Part. Nucl. Phys. 72, 99 (2013), eprint 1301.6377.
  * Odyniec [2013] G. Odyniec, PoS CPOD2013, 043 (2013).
  * [12] _https://www.gsi.de/work/forschung/cbmnqm/cbm.htm_.
  * Senger [2022] P. Senger, Phys. Scripta 97, 064003 (2022), eprint 2204.01056.
  * Kekelidze et al. [2016] V. D. Kekelidze, R. Lednicky, V. A. Matveev, I. N. Meshkov, A. S. Sorin, and G. V. Trubnikov, Eur. Phys. J. A 52, 211 (2016).
  * Burmasov et al. [2022] N. Burmasov, E. Kryshen, D. Ivanishchev, D. Kotov, M. Malaev, V. Riabov, and Y. Ryabov, Phys. Scripta 97, 054013 (2022).
  * Alt et al. [2008] C. Alt et al. (NA49), Phys. Rev. C 78, 034914 (2008), eprint 0712.3216.
  * Adamczyk et al. [2014] L. Adamczyk et al. (STAR), Phys. Rev. Lett. 112, 032302 (2014), eprint 1309.5681.
  * Maćkowiak-Pawłowska and Wilczek [2014] M. Maćkowiak-Pawłowska and A. Wilczek (NA61/SHINE), J. Phys. Conf. Ser. 509, 012044 (2014), eprint 1402.0707.
  * Aduszkiewicz et al. [2016] A. Aduszkiewicz et al. (NA61/SHINE), Eur. Phys. J. C 76, 635 (2016), eprint 1510.00163.
  * Adare et al. [2016] A. Adare et al. (PHENIX), Phys. Rev. C 93, 011901 (2016), eprint 1506.07834.
  * Luo [2015] X. Luo (STAR), PoS CPOD2014, 019 (2015), eprint 1503.02558.
  * Luo and Xu [2017] X. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017), eprint 1701.02105.
  * Esha [2017] R. Esha (STAR), Nucl. Phys. A 967, 457 (2017).
  * Adamczyk et al. [2018] L. Adamczyk et al. (STAR), Phys. Lett. B 785, 551 (2018), eprint 1709.00773.
  * Prokhorova [2018] D. Prokhorova (NA61/SHINE), KnE Energ. Phys. 3, 217 (2018), eprint 1801.06690.
  * Xu [2018] N. Xu, Springer Proc. Phys. 203, 1 (2018).
  * Andronov [2018] E. Andronov (NA61/SHINE), KnE Energ. Phys. 3, 226 (2018), eprint 1801.03711.
  * Behera [2019] N. K. Behera (ALICE), Nucl. Phys. A 982, 851 (2019), eprint 1807.06780.
  * Adam et al. [2021] J. Adam et al. (STAR), Phys. Rev. Lett. 126, 092301 (2021), eprint 2001.02852.
  * Adam et al. [2020] J. Adam et al. (STAR), Phys. Rev. C 102, 024903 (2020), eprint 2001.06419.
  * Abdallah et al. [2021a] M. Abdallah et al. (STAR), Phys. Rev. C 104, 024902 (2021a), eprint 2101.12413.
  * Abdallah et al. [2022] M. S. Abdallah et al. (STAR), Phys. Rev. Lett. 128, 202303 (2022), eprint 2112.00240.
  * Abdallah et al. [2021b] M. Abdallah et al. (STAR), Phys. Rev. Lett. 127, 262301 (2021b), eprint 2105.14698.
  * Aparin [2022] A. Aparin (STAR), Phys. Part. Nucl. 53, 127 (2022).
  * Karsch et al. [2008] F. Karsch, D. Kharzeev, and K. Tuchin, Phys. Lett. B 663, 217 (2008), eprint 0711.0914.
  * Minami [2011] Y. Minami, Phys. Rev. D 83, 094019 (2011), eprint 1102.5485.
  * De Groot [1980] S. R. De Groot, _Relativistic Kinetic Theory. Principles and Applications_ (1980).
  * Mori and Fujisaka [1973] H. Mori and H. Fujisaka, Progress of Theoretical Physics 49, 764 (1973), ISSN 0033-068X, eprint https://academic.oup.com/ptp/article-pdf/49/3/764/5262040/49-3-764.pdf, URL https://doi.org/10.1143/PTP.49.764.
  * Kawasaki [1970] K. Kawasaki, Annals of Physics 61, 1 (1970), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/0003491670903751.
  * Mori [1965] H. Mori, Progress of Theoretical Physics 33, 423 (1965), ISSN 0033-068X, eprint https://academic.oup.com/ptp/article-pdf/33/3/423/5428510/33-3-423.pdf, URL https://doi.org/10.1143/PTP.33.423.
  * Zwanzig [2001] R. Zwanzig, _Nonequilibium statistical mechanics_ (Oxford University Press, 2001).
  * Landau and Lifshitz [1987] L. Landau and E. Lifshitz, _Fluid Mechanics_ (Butterworth Heinemann, Oxford, UK, 1987), second edition ed.
  * Goldstone [1961] J. Goldstone, Nuovo Cim. 19, 154 (1961).
  * Nambu and Jona-Lasinio [1961] Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961).
  * Mazenko. [2006] G. F. Mazenko., _Nonequilibrium Statistical Mechanics_ (2006).
  * Hohenberg and Halperin [1977] P. C. Hohenberg and B. I. Halperin, Review of Modern Physics 49, 435 (1977).
  * Chandrasekhar [1943] S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943).
  * Wang and Uhlenbeck [1945] M. C. Wang and G. E. Uhlenbeck, Rev. Mod. Phys. 17, 323 (1945), URL https://link.aps.org/doi/10.1103/RevModPhys.17.323.
  * Onuki [2007] A. Onuki, _Phase Transition Dynamics_ (Cambridge University Press, 2007).
  * Balucani et al. [2003] U. Balucani, M. H. Lee, and V. Tognetti, Physics reports 373, 409 (2003).
  * Grabert [1982] H. Grabert, _Projection Operator Techniques in Nonequilibrium Statistical Mechanics_ (Springer Tracts in Modern Physics, 95, 1982).
  * Huang et al. [2011] X.-G. Huang, T. Kodama, T. Koide, and D. H. Rischke, Phys. Rev. C 83, 024906 (2011), eprint 1010.4359.
  * Minami and Hidaka [2013] Y. Minami and Y. Hidaka, Phys. Rev. E 87, 023007 (2013), eprint 1210.1313.
  * Hayata and Hidaka [2015] T. Hayata and Y. Hidaka, Phys. Rev. D 91, 056006 (2015), eprint 1406.6271.
  * Morozov [1981] V. G. Morozov, TMF 48, 373 (1981).
  * Morozov and Mukhai [1992] V. G. Morozov and A. N. Mukhai, Theoretical and Mathematical Physics 90, 189–204 (1992).
  * Nordholm and Zwanzig. [1975] S. Nordholm and R. Zwanzig., Journal of Statistical Physics 13, 347–371 (1975).
  * H. Grabert and Talkner. [1980] P. H. i. H. Grabert and P. Talkner., Journal of Statistical Physics 22, 537–552 (1980).
  * Zubarev and Morozov [1983] D. Zubarev and V. Morozov, Physica A: Statistical Mechanics and its Applications 120, 411 (1983), ISSN 0378-4371, URL https://www.sciencedirect.com/science/article/pii/0378437183900626.
  * Balescu [1975] R. Balescu, _Equilibrium and Nonequilibrium Statistical Mechanics_ (A Wiley-Interscience, 1975).
  * Weyl [1928] H. Weyl, _Gruppentheorie und Quantenmechanik_ (Leipzig: Hirzel, 1928).
  * Zubarev [1974] D. N. Zubarev, _Nonequilibrium Statistical Thermodynamics_ (Plenum, New York, 1974).
  * Hu [2023] J. Hu, Phys. Rev. C 107, 024915 (2023), eprint 2209.10979.
  * Zwanzig [1961] R. Zwanzig, Phys. Rev. 124, 983 (1961), URL https://link.aps.org/doi/10.1103/PhysRev.124.983.
  * Schilling [2022] T. Schilling, Physics Reports 972, 1 (2022), ISSN 0370-1573, coarse-Grained Modelling Out of Equilibrium, URL https://www.sciencedirect.com/science/article/pii/S0370157322001314.
  * Kubo. [1959] R. Kubo., _Lectures in theoretical physics, Volumn I_ (Interscience, 1959).
  * Kim and Mazenko. [1991] B. Kim and G. F. Mazenko., Journal of Statistical Physics 64, 631–652 (1991).
  * Pomeau and Résibois [1975] Y. Pomeau and P. Résibois, Physics Reports 19, 63 (1975), ISSN 0370-1573, URL https://www.sciencedirect.com/science/article/pii/0370157375900198.


# Geodesic congruences in acoustic spacetimes and the role of Raychaudhuri
equation

Akshat<EMAIL_ADDRESS>
Department of Physics, Shiv Nadar Institution of Eminence
Greater Noida, Uttar Pradesh-201314, India

###### Abstract

It has been known that the propagation of sound in fluids can be used to model
acoustic spacetimes. These acoustic spacetimes offer analogue models for
gravity. We use the Raychaudhuri equation to study the propagation of sound in
these fluids, which, via the Eikonal approximation, correspond to null
geodesic congruences in the acoustic spacetimes. We explore this within the
acoustic analogues of black holes and cosmological spacetimes. The robustness
of the Raychaudhuri equation and the limits of the acoustic analogue are
emphasised.

##  1 Introduction

There has been a continual interest in the analogue gravity literature to
describe fluids in terms of a Lorentzian metric in which an acoustic coupling
is explicit [1, 2, 3]. Within this framework, the propagation of sound is
described by by a scalar field $\phi$ obeying the massless minimally coupled
scalar field equation

$\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\phi\right)=0$ (1)

Here $g_{\mu\nu}$ is the Lorentzian acoustic metric describing the acoustic
spacetime. The line element corresponding to the acoustic metric is

$\mathrm{d}s^{2}=\frac{\rho_{0}(\vec{x},t)}{c_{s}(\vec{x},t)}\left[-c_{s}^{2}\mathrm{~{}d}t^{2}+\delta_{ij}\left(\mathrm{d}x^{i}-v_{0}^{i}\mathrm{~{}d}t\right)\left(\mathrm{~{}d}x^{j}-v_{0}^{j}\mathrm{~{}d}t\right)\right]$
(2)

Here, $\rho_{0}$ is the background fluid density, $c_{s}$ is the speed of
sound and $v^{i}$ is the fluid velocity field. Imposing additional conditions
on these variables will correspond to specific analogue spacetimes.

Further, the Eikonal approximation will be used in describing the scalar
field. Within this approximation, we expect the sound waves to travel along
null geodesics while the congruences obey the Raychaudhuri equation [4, 5, 6]
for null case. Written in full, in (3+1) dimensions, the latter can be cast in
the form

$\frac{d\Theta}{d\lambda}=-\frac{1}{2}\Theta^{2}-\sigma^{\mu\nu}\sigma_{\mu\nu}+\omega^{\mu\nu}\omega_{\mu\nu}-R_{\mu\nu}k^{\mu}k^{\nu}$
(3)

where $\Theta$ is the expansion scalar, $\sigma_{\mu\nu}$ is the shear tensor,
$\omega_{\mu\nu}$ is the rotation tensor, and the $R_{\mu\nu}$ is the usual
Ricci tensor. Further, $\lambda$ is the affine parameter, and $k^{\mu}$ is the
vector field which is tangent to the geodesics. Note that $\Theta$ can be
represented by the divergence form $\Theta=\nabla_{\mu}k^{\mu}$. The acoustic
coupling, which is constructed within the confines of non-relativistic fluid
mechanics, gives way to analogue models as counterparts of general relativity
(GR) spacetimes that embrace the kinematic effects of GR [2]. The purpose of
this work is to look at the behaviour of the congruence of geodesics in such
analogue spacetimes. Another interesting query that deserves a close study is
the role of Raychaudhuri equation which is basically concerned with the
background geometry of the manifold rather than addressing any specific
feature of GR which dictates, in particular, the kinematics of the model
involved. [8, 9, 10]

The present paper is organized as follows. In section 2, we briefly sketch how
the Eikonal approximation. This lets us describe the sound waves using null
geodesics. In section 3, we address acoustic black holes, in particular the
canonical acoustic black hole. We look at null geodesic congruences, derive
the Raychaudhuri equation and show that it is in accordance with what is
obtained from the Schwarzschild metric. We briefly discuss other acoustic
black hole models. In section 4, we discuss the acoustic analogue of the FLRW
cosmology, wherein the curvature of the spacetime is obtained by varying the
speed of sound with time. In this context, we explore null geodesic
congruences and the Raychaudhuri equation and try to interpret their physical
meaning. Finally, a summary of our work is presented.

##  2 The Eikonal Approximation

Known also as the geometrical optics approximation, the Eikonal approximation
is an analytical tool to approximate waves with rays, and hence with geodesics
when the space(time) is curved. For a scalar field the approximation gives:

$\phi(r,t)=\mathcal{A}(r,t)\exp[\mp i\varphi(r,t)]$ (4)

The spatial and temporal dependence of the phase $\varphi(r,t)$ can be further
separated, given the geometry is slowly evolving
($\omega\gg\max\\{|\dot{c}/c|,|\dot{v}/v|\\}$)

$\varphi(r,t)=\omega t-\int^{r}K\left(r^{\prime}\right)\mathrm{d}r^{\prime}$
(5)

The condition of restricting the analysis to high values of $\omega$ is at the
core of the approximation.

Plugging for $\phi$ into the minimally coupled scalar field equation,

$\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\phi\right)=0$ (6)

The equation corresponding to $\varphi$ turns out to be

$g^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi=0$ (7)

Upon making the identification $\partial_{\mu}\varphi=k_{\mu}$, we get
$k_{\mu}$ to be a null vector such that

$k^{\mu}k_{\mu}=0$ (8)

Acting with the covariant derivative on equation(8) and exploiting the fact
that it is symmetric $\nabla_{\mu}k_{\nu}=\nabla_{\nu}k_{\mu}$, we end up with

$k^{\nu}\nabla_{\nu}k_{\mu}=0$ (9)

These in turn are the geodesic equations. We see that within the validity of
the Eikonal approximation, the waves travel along null geodesics. Choosing an
entire non-intersecting family of such geodesics, we can construct null
geodesic congruences. These congruences, representing high frequency sound
wavefronts, would thus obey the Raychaudhuri equation for the null case, as we
shall explore in further sections.

##  3 Acoustic Black holes

We will begin by looking at acoustic analogues to black hole spacetimes. In
this case, the analogy is theoretically well established with the fluid
velocity being responsible for the curvature. Here we look at the two most
well known models – the Draining bathtub and the Canonical acoustic
spacetimes. The (2+1) dimensional draining bath-tub, which is a rotating
solution, is especially relevant for its experimental feasibility. The
canonical acoustic solution is stationary and spherically symmetric, and thus
offers a close analogue for Schwarszchild-like spacetimes.

###  3.1 The (2+1)-dimensional draining bathtub spacetime

Starting with equation (2) a draining bathtub model in (2+1) dimensions can be
constructed that has a sink at the origin. The conditions imposed on the fluid
are consistent with the velocity potential

$\vec{v}=\frac{(A\hat{r}+B\hat{\theta})}{r}$ (10)

It leads to the following form of the acoustic metric

$ds^{2}=-dt^{2}+\left(dr-\frac{A}{r}dt\right)^{2}+\left(rd\theta-\frac{B}{r}dt\right)^{2}$
(11)

which can be reset as

$ds^{2}=-\left(1-\frac{A^{2}+B^{2}}{r^{2}}\right)dt^{2}-2\frac{A}{r}drdt-2Bd\theta
dt+dr^{2}+r^{2}d\theta^{2}$ (12)

We must note, that we have implicitly assumed cylindrical symmetry and chosen
a constant $z$ slice to reduce the number of effective dimensions. It is
worthwhile to explore the behaviour of the null geodesic congruences in such a
scenario essentially because it admits both the features of rotation, and a
vortex, and is the closest known analogue to Kerr black holes.

Geodesic congruences in the context of Draining bathtub spacetimes of a
similar kind, have been extensively discussed in the literature by Dolan et al
[11, 12, 13].

However, for completeness, we will briefly sketch the form of the null
Raychaudhuri equation. Firstly we note that the rotation tensor
$\omega_{\mu\nu}$ vanishes. This is because, in the case of outward (or
inward) pointing geodesics, the congruences are hypersurface-orthoghonal. This
is, in fact, closely related to the definition of the rotation tensor to be
symmetric. Further, the shear tensor vanishes identically. This is related to
the fact that the shear tensor $\sigma_{\mu\nu}$ is traceless and in (2+1)
dimensions the null congruence is codimension one. Therefore, there is no
trace-free contribution.

The Raychaudhuri equation in (2+1) dimensions thus reduces to

$\frac{d\Theta}{d\lambda}=-\Theta^{2}-R_{\mu\nu}k^{\mu}k^{\nu}$ (13)

Further, the Ricci tensor vanishes. This is to be expected, as we are looking
at analogues to spacetimes which are vaccumm solutions to the Einstein
equations. We are thus left with

$\frac{d\Theta}{d\lambda}=-\Theta^{2}$ (14)

$\Theta$ upon solving comes out to be $\frac{1}{r}$. The shortcomings of the
expansional scalar in this context and the resolution in terms of the Van
Vleck determinant are discussed in [14].

###  3.2 The Canonical Acoustic Black hole spacetime

We start with the incompressible acoustic metric (2) which reads

$\mathrm{d}s^{2}=\frac{\rho_{0}}{c_{s}}\left[-c_{s}^{2}\mathrm{~{}d}t^{2}+\delta_{ij}\left(\mathrm{d}x^{i}-v_{0}^{i}\mathrm{~{}d}t\right)\left(\mathrm{~{}d}x^{j}-v_{0}^{j}\mathrm{~{}d}t\right)\right]$
(15)

For this case the speed of sound $c_{s}$ is taken to a constant independent of
time. This leads to the time independence of $\rho_{0}$, as $c_{s}$ is merely
a relation between fluid density and the corresponding pressure. Further,
imposing spherical symmetry it implies that $\rho_{0}$ has to be position
independent. This, from the assumption of the barotropic nature of the fluid,
points to the position-independence of pressure as well as that of the speed
of the sound. Since from the continuity equation it transpires that
$v_{0}\propto 1/r^{2}$,

We can project, up to a normalisation constant $r_{0}$, the fluid velocity in
the manner

$v_{0}=c_{s}\frac{r_{0}^{2}}{r^{2}}$ (16)

Ignoring an irrelevant position dependent factor, we can express the acoustic
metric in the form

$ds^{2}=-c_{s}^{2}dt^{2}+\left(dr\pm
c_{s}\frac{r_{0}^{2}}{r^{2}}dt\right)^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right).$ (17)

At this point, we observe that on making the following coordinate
transformation

$d\tau=dt\mp\frac{r_{0}^{2}/r^{2}}{c_{s}\left[1-\left(r_{0}^{4}/r^{4}\right)\right]}dr$
(18)

We can map the line element (17) to the canonical structure [2]

$ds^{2}=-c_{s}^{2}\left[1-\left(r_{0}^{4}/r^{4}\right)\right]d\tau^{2}+\frac{dr^{2}}{1-\left(r_{0}^{4}/r^{4}\right)}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)$ (19)

Remark: The $r$ dependence of the above line element is different from that of
the Schwarszchild line element. However, in the near horizon approximation
$(r\approx r_{0})$, the $r$ dependence, upon factorisation can be put in a
form proportional to $1/r$.

For the rest of the section we assume, without loss of generality, the
convention $c_{s}=1$ and $r_{0}=1$

We now define a quantity $\beta$ as given by

$\beta(r)=1-\frac{1}{r^{4}}$ (20)

This transforms (19) to

$ds^{2}=-\beta(r)d\tau^{2}+(\beta(r))^{-1}dr^{2}+r^{2}d\Omega^{2}$ (21)

For radial null geodesics, $d\theta=d\phi=0$ which means that the rotation
($\omega_{\mu\nu}$) and stress ($\sigma_{\mu\nu}$) tensors do not contribute
and the line element is simplified to

$ds^{2}=-\beta(r)d\tau^{2}+(\beta(r))^{-1}dr^{2}=0$ (22)

Introducing the so-called tortoise coordinate $r_{*}$ satisfying the
differential equation

$\frac{dr_{*}}{dr}=\sqrt{-\frac{g_{rr}}{g_{00}}}=\sqrt{\frac{1}{\beta^{2}(r)}}=\frac{1}{\beta(r)}$
(23)

and which integrates to

$r_{*}=\frac{-\ln(|r+1|)-2\tan^{-1}(r)+\ln(|r-1|)}{4}+r$ (24)

the line element in terms of $r_{*}$ becomes

$ds^{2}=\beta(r)\left(-d\tau^{2}+dr_{*}^{2}\right)=0$ (25)

In terms of Eddington-Finkelstein coordinates which defines a pair of
coordinates $(u,v)$

$u=\tau-r_{*},\quad v=\tau+r_{*}$ (26)

we easily notice that $u$ is constant for outgoing null geodesics, while $v$
is constant for incoming null geodesics. Focusing on the outgoing geodesics,
we define the corresponding tangent vector field

$k_{\mu}=-\partial_{\mu}u=(-1,(\beta(r))^{-1},0,0)$ (27)

such that

$k^{\mu}=g^{\mu\nu}k_{\nu}=\left((\beta(r))^{-1},1,0,0\right)$ (28)

The calculations are relevant for the consequences we draw at the end. Hence
we will put in several of the steps. The expansion scalar $\Theta$ can now be
calculated by using its representation

$\Theta=\nabla_{\mu}k^{\mu}=\partial_{\mu}k^{\mu}+\Gamma_{\mu_{\alpha}}^{\mu}k^{\alpha}$
(29)

where the non-vanishing terms are

$\Theta=\Gamma_{01}^{0}k^{1}+\Gamma_{11}^{1}k^{1}+\Gamma_{21}^{2}k^{1}+\Gamma_{31}^{3}k^{1}$
(30)

Noting that

$\Gamma_{01}^{0}k^{1}=\frac{1}{2}\left((\beta(r))^{-1}\right)\frac{d\beta(r)}{dr}$

$\Gamma_{11}^{1}k^{1}=-\frac{1}{2}\left((\beta(r))^{-1}\right)\frac{d\beta(r)}{dr}$

$\Gamma_{21}^{2}k^{1}=\frac{1}{r}$ $\Gamma_{31}^{3}k^{1}=\frac{1}{r}$

$\Theta$ is easily estimated to be

$\Theta=\frac{2}{r}$ (31)

A similar calculation can be performed for incoming geodesics with $k_{\mu}$
defined as

$k_{\mu}=-\partial_{\mu}v$ (32)

Then $\Theta$ turns out to be

$\Theta=-\frac{2}{r}$ (33)

Notice this is twice of what we mentioned for the (2+1) dimensional spacetime.
This is because in this case, there is an extra contribution in the sum from
the third spatial dimension.

As, in the previous case, $\sigma_{\mu\nu}$, $\omega_{\mu\nu}$ and
$R_{\mu\nu}$ all vanish.

Raychaudhuri equation for either case yields

$\frac{d\Theta}{d\lambda}=-\frac{2}{r^{2}}$ (34)

Equations (31) and (33) and equation (34) match exactly with the result
following from the Schwarzschild metric [7]. This is because, in the present
calculation the explicit form of $\beta(r)$ was not required to be
implemented. Note that $\beta(r)=1-\frac{2M}{r}$ would correspond to the
Schwarzschild case. Thus our result for null congruences would hold more
generally for any time independent, spherically symmetric metric which
satisfies the form prescribed in equation (21).

As a trivial consequence, we see that the above result holds for $\beta(r)=1$.
Within the acoustic framework, this can be obtained by taking the fluid
velocity 3-field $v_{0}=0$ everywhere in equation (2) and further taking
$c_{s}$ to be a constant in time. This is analogous to the Minkowski
spacetime. Therefore, in homogenous fluids, the dynamics of high frequency
sound waves can be studied using the Raychaudhuri equation. It would be
interesting to study non-radial geodesics for this case, which we shall try to
explore in a future work.

##  4 Acoustic FLRW cosmologies

We now turn towards an analogue cosmological model obtained from the acoustic
metric. As Barcelo et al [15] pointed out, the appropriate acoustic analogue
of a flat FLRW metric, is obtained by keeping the background fluid at rest
$(v_{0}=0)$ and instead varying the speed of the excitations ($c_{s}(t)$)
within them.

Starting with equation (2), in spherical coordinates two more conditions are
imposed.

  1. 1.

The background fluid is considered to be at rest, $v_{z}=0$. The continuity
equation gives that $\rho_{0}$ is a constant.

  2. 2.

Spatially homogenity is assumed. This implies that $c_{s}$ is independent of
$r$.

These lead to

$\mathrm{d}s^{2}=\frac{\rho_{0}}{c_{s}(t)}\left[-c_{s}^{2}\mathrm{~{}d}t^{2}+\mathrm{d}\vec{x}^{2}\right].$
(35)

Further, the metric can be rescaled by a constant factor
$\frac{c_{0}}{\rho_{0}}$. Here $c_{0}$ is some convenient reference speed. The
line element thus becomes

$\mathrm{d}s^{2}=\frac{c_{0}}{c_{s}(t)}\left[-c_{s}^{2}(t)\mathrm{~{}d}t^{2}+\mathrm{d}z^{2}\right]=-c_{0}c_{s}(t)\mathrm{~{}d}t^{2}+.\frac{c_{0}}{c_{s}(t)}\mathrm{d}\vec{x}^{2}$
(36)

In order to convert this line element into the traditional FLRW form, a pseudo
time is introduced. This pseudo time $\tau$ is related to the laboratory time
$t$ such that

$d\tau=dt\sqrt{c_{s}(t)/c_{0}}$ (37)

Equation (36) becomes

$\mathrm{d}s^{2}=-c_{0}^{2}\mathrm{~{}d}\tau^{2}+\frac{c_{0}}{c_{s}(t)}\mathrm{d}\vec{x}^{2}=-c_{0}^{2}\mathrm{~{}d}\tau^{2}+\frac{c_{0}^{2}}{\bar{c}_{s}^{2}(\tau)}\mathrm{d}\vec{x}^{2}$
(38)

Here, $\bar{c}_{s}(\tau)$ is the speed of sound in terms of the pseudo time.
This can be obtained via the chain rule

$\bar{c}_{s}(\tau)=\frac{d\vec{x}}{d\tau}=\frac{dt}{d\tau}\frac{d\vec{x}}{dt}=\sqrt{c_{0}c_{s}(t)}$
(39)

Equation (38) is completely equivalent to the flat FLRW line element upon
making the identification

$c_{0}d\tau\sim cdt_{F}$

and

$\frac{c_{0}^{2}}{\bar{c}_{s}^{2}(\tau)}\sim a^{2}(t_{F})$

Here $c$ is the speed of light and $t_{F}$ is the comoving FLRW time. It is
worth re-emphasising that within this model the curvature is obtained by
adding time dependence to the speed of sound. This is unlike the case for
acoustic black holes where the speed of sound was taken to be constant in
time, and the curvature was a consequence of inhomogenous fluid flow. This
particular analogue model is also interesting because here we have an analogy
with a non-vaccumm solution of GR. For works exploring analogue cosmologies
see, [16, 18, 17, 19].

For convenience of notation we will choose $c_{0}=1$. Having obtained the
metric, we would like to construct the various geometric quantities out of it.
We shall, again, focus on radial null geodesics. For outgoing geodesics, we
define the tangent vector field

$k_{\mu}=-\partial_{\mu}\left(\bar{c}_{s}d\tau-
dr\right)=\left(-\bar{c}_{s}(\tau),1,0,0\right)$ (40)

Such that

$k^{\mu}=\left(\bar{c}_{s}(\tau),\bar{c}^{2}_{s}(\tau),0,0\right)$ (41)

This ensures, $k_{\mu}k^{\mu}=0$

Upon solving for the expansion factor, we end up with

$\Theta=2\bar{c}_{s}^{2}(\tau)\left(-\frac{\dot{\bar{c}}_{s}(\tau)}{\bar{c}_{s}(\tau)}+\frac{1}{r}\right)$
(42)

The dot represents a derivative with respect to the Pseudo time $\tau$. Note
that equation (39) tells us that this pseudo time is linearly related to the
laboratory time $t$. This is important as in this case, the qualitative
features of the dynamics will get carried over, even when using $t$.

There are a couple of things to notice about the expression for $\Theta$.
Firstly, for $\dot{\bar{c}}_{s}(\tau)=0$, the expression reduces to

$\Theta=\bar{c}_{s}^{2}(\tau)\frac{2}{r}$

This is upto a constant, the same as equation (26). This was expected as
$\dot{\bar{c}}_{s}(\tau)=0$ would correspond to a static (and flat) acoustic
spacetime. Further, the minus sign in the $\dot{\bar{c}}_{s}(\tau)$ term in
(44) ensures that in order to get an expanding (analogue) spacetime, a
decreasing speed of sound is required.

Now, we move on to calculate the Ricci tensor. The relevant non-vanishing
terms are

$R_{00}=3\left(\frac{\ddot{\bar{c}}_{s}}{c_{s}}-\frac{1}{2}\frac{\dot{\bar{c}}^{2}_{s}}{\bar{c}^{2}_{s}}\right)$
(43)

and

$R_{11}=\frac{5}{4}\frac{\dot{\bar{c}}^{2}_{s}}{\bar{c}^{4}_{s}}-\frac{\ddot{\bar{c}}_{s}}{c^{3}_{s}}$
(44)

Here, $\bar{c}_{s}(\tau)$ is written as $\bar{c}_{s}$ to avoid cluttering the
equations.

Spatial homogeneity and isotropy together with the fact that we are looking at
hypersurface orthogonal geodesics ensure that $\sigma_{\mu\nu}$ and
$\omega_{\mu\nu}$ will vanish.

Writing out the Raychaudhuri equation in full gives us

$\displaystyle\frac{d\Theta}{d\lambda}$
$\displaystyle=-\frac{1}{2}\Theta^{2}-R_{\mu\nu}k^{\mu}k^{\nu}$ (45)
$\displaystyle=\ddot{\bar{c}}_{s}\left(\frac{1}{\bar{c}_{s}}-3\bar{c}_{s}\right)+\dot{\bar{c}}_{s}\left(\frac{7}{4}\dot{\bar{c}}_{s}-2c^{2}_{s}\dot{\bar{c}}_{s}+\frac{4}{r\bar{c}^{3}_{s}}\right)+\frac{2\bar{c}^{4}_{s}}{r^{2}}$

This equation puts forward the various competing terms generating the dynamics
of the expansion scalar. As expected there is a contribution from the second
derivative of the scale factor (inverse of $c_{s}$). It could be interesting
to think about this equation in a fluid mechanical context. It might give
alternate geometry based explanations to known results.

However, this equation does not quite resemble the form for the null
Raychaudhuri equation as is usually seen in cosmology textbooks. This is so
because, in writing down that particular form of the Raychaudhuri equation,
the Freidmann equation is implicitly used. Our inability to do that comes from
the fact that the acoustic analogue of gravity can mimic only the kinematic
effects of GR. That is, within the acoustic framework the dynamics of GR,
namely, Einstein equations and consequently the (first) Freidmann equation
cannot be reproduced. Hence we do not expect here to obtain the form involving
terms from the energy-momentum tensor.

##  5 Summary

In this work, we emphasised the robustness of the Raychaudhuri equation by
pointing out its relevance in an analogue spacetime situation where we did not
have to make any use of Einstein equations or the null energy condition [20,
21]. We noted the similarities and differences in the results we obtained here
with those in GR. Therefore, we pointed out the limits of the acoustic
analogue for gravity, which is well known to mimic only the kinematic effects
of GR[22]. Finally, we provided an illustrative example of how the
mathematical methods constructed for GR could be used to understand the
behaviour of sound in fluids.

We would also like to mention a query which we believe deserves some
investigation. It is the following — Upon quantising the fluid excitations, we
end up with phonons. Analysing the validity and consequences of the Quantum
Raychaudhuri Equation (QRE) [23] in this framework could be valuable for
pracitioners of both analogue gravity and phonon-based physics. It could
possibly be of interest to physicists studying the foundations of quantum
mechanics as the QRE is set within the Bohmian [24, 25] formulation of Quantum
Mechanics.

## Acknowledgements

I wish to thank Prof. Bijan Bagchi for the insightful discussions and for his
comments on the early draft of this manuscript. I also wish to thank Rahul
Ghosh and Sauvik Sen for their valuable feedback.

## References

  * [1] Unruh W. G. 1981 Phys. Rev. Lett. 46, 1351.
  * [2] Visser M. 1998 Class. Quant. Grav. 15, 1767.
  * [3] Barceló C., Liberati S., and Visser M. 2011 Living Rev. Relativ. 14 1-159.
  * [4] Raychaudhuri A. 1955  Phys. Rev. 98, 1123.
  * [5] Lifshitz E.M. and Khalatnikov I. M. 1963 Adv. Phys. 12, 185.
  * [6] Ellis G. F. 2007  Pramana, 69.
  * [7] Poisson E. 2004 A Relativist’s Toolkit, 250.
  * [8] Kar S. 2001 Phys. Rev. D64, 105017.
  * [9] Kar S. and SenGupta S. 2007 Pramana 69, 49.
  * [10] Kar S. 1996 Phys. Rev. D54, 6408.
  * [11] Torres, T., Lloyd M., Dolan S.R. and Weinfurtner S. 2022 Phys. Rev. Research 4(3).
  * [12] Dempsey D. and Dolan S.R. 2016 Int. J. Mod. Phys. D 25 (09).
  * [13] Dolan S.R., Oliveira L.A. and Crispino L.C. 2012 Phys. Rev. D 85 (4).
  * [14] Dempsey D. 2017 Wave propagation on black hole spacetimes.
  * [15] Barceló C., Liberati, S. and Visser M. 2003 Int. J. Mod. Phys. D12, 09.
  * [16] Jain P, Weinfurtner S., Visser M. and Gardiner C. W. 2007 Phys. Rev. A 76, 033616.
  * [17] Pandey A. 2023 A note on Analogue semi-classical gravity in (1+1) dimensions. (unpublished).
  * [18] Lang S. and Schützhold R. 2019 Phys. Rev. D 100, 065003.
  * [19] Barceló C., Liberati S. and Visser M. 2003 Phys. Rev. A 68, 053613.
  * [20] Visser M. 1996 Lorentzian Wormholes. From Einstein to Hawking.
  * [21] Wald R. 1984 General Relativity.
  * [22] Visser M. 2003 Int. J. Mod. Phys., D12.
  * [23] Das S. 2014 Physical Review D89 8\.
  * [24] Bohm D. 1952 Phys. Rev. 85 166\.
  * [25] Bohm D., Hiley B. J. and Kaloyerou P. N. 1987 Phys. Rep. 144 321\.