compat=1.1.0

# On the correlation functions in stable first-order relativistic
hydrodynamics

Navid<EMAIL_ADDRESS>Sara<EMAIL_ADDRESS>
__
aSchool of Nuclear Science and Technology, Lanzhou University,
222 South Tianshui Road, Lanzhou 730000, China
__ b Institute of Modern Physics, Chinese Academy of Sciences,
Lanzhou 730000, China

First-order relativistic conformal hydrodynamics in a general frame is
characterized by three transport coefficients. Within a certain range of these
coefficients, the equilibrium is stable and propagation is causal. By
explicitly computing relevant hydrodynamic correlation functions in the linear
response regime, we derive associated Kubo formulas for the above mentioned
coefficients. Then we go beyond the linear response to explore the analytic
structure of response functions in stable frames. A method for numerically
finding the branch cut structure from the well-known Landau equations is
developed. We then apply our method to the shear channel and find the branch
cuts of a certain response function, without computing the response function
itself. We also solve the Landau equations analytically and find the threshold
singularities of the same response function. The results turn out to be in
complete agreement with the loop calculations in effective field theory.

###### Contents

  1. 1 Introduction
  2. 2 First-order stable conformal hydrodynamics
  3. 3 Linear response
    1. 3.1 Correlation functions in longitudinal channel
    2. 3.2 Correlation functions in transverse channel
    3. 3.3 Kubo Formulas
  4. 4 Correlation functions beyond linear response
    1. 4.1 Correlation functions of the shear stress
    2. 4.2 Solving the Landau equations
  5. 5 Discussion and outlook
  6. A Equilibrium fluctuations
  7. B Hydrodynamic fluctuations
  8. C Landau equations

##  1 Introduction

It is well known that relativistic hydrodynamics beyond zeroth order in
derivatives can be constructed in various ways, depending on how the fluid
velocity and temperature are defined when the system is not in equilibrium
[1]. Any particular choice of definition for these out of equilibrium
variables is called a choice of “frame”333It does not have to be confused with
the concept of Lorentz frames.. In the famous Eckart [2] and Landau-Lifshitz
[3] frames, hydrodynamics exhibits some unphyiscal features: the thermal
equilibrium is unstable [4] and the theory supports the propagation of
superluminal modes [5].

Problems of instability and acausality are associated with modes outside the
range of validity of hydrodynamics, i.e., beyond the long wavelength low
frequency limit. In other words, the problem is due to the UV modes. To remedy
this, one way is to introduce extra dynamical degrees of freedom into
hydrodynamics, thereby modifying the behavior of the theory in the UV. This is
actually the idea of the Israel-Stewart theory [6, 7].

A more systematic way of dealing with this problem is to find certain out of
equilibrium definitions for hydrodynamic variables; something other than those
used by Eckart or Landau-Lifshitz. This program was initiated by [8, 9], and
built upon the work of Bemfica, Disconzi and Noronha, in the case of conformal
hydrodynamics in [10]. Finally, by constructing the theory of relativistic
hydrodynamics within a general frame, Kovtun introduced a ‘stable first-order
theory of hydrodynamics’ [11]. These theories ([10, 11]) are sometimes
referred to as BDNK in the literature 444After [10] and [11], many works have
investigated relativistic hydrodynamics in general frames as well as
stability/causality problems. See [27, 28, 12, 19, 31, 29, 30, 13, 17, 18, 25,
21, 16, 26, 24, 23, 22, 20, 14, 15] for an incomplete list of related works. .

In practice, causal relativistic hydrodynamics 555As mentioned earlier,
causality can be realized in the Muller-Israel-Stewart theory, or through the
second order relativistic hydrodynamics [32]. can describe the evolution of
hot QCD matter formed in heavy ion collisions (see [33, 34] for a review). One
may naturally expect the BDNK theory would also be able to describe the
evolution of quark-gluon plasma (QGP). Taking QGP as the conformal system, it
turns out that BDNK theory introduces two new transport coefficients in
addition to the shear viscosity coefficient. The latter two are responsible
for the relaxation of energy density and momentum density. Any numerical
modeling of QGP based on BDNK equations then requires specification of these
transport coefficients. The equations of BDNK and the transport coefficients
will be reviewed briefly in § 2.

The first goal of this work is to derive the Kubo formulas related to the
transport coefficients of the BDNK theory. To this end, in § 3, we first
compute the correlation function of the energy and momentum densities in the
linear response regime, for both sound and shear channels. Then from the
correlation functions of the sound channel, we derive a set of Kubo formulas
for the above transport coefficients. These Kubo formulas can be tested
experimentally. Let us assume that in a particular system the correlation
function of the energy and momentum densities can be extracted from some
experimental data, in the long wavelength limit. The value of the BDNK
transport coefficients associated with this system can then be found using the
Kubo formulas. If it is determined that the two relaxation coefficients are
non-zero and satisfy both the stability and causality conditions, it can be
concluded that the BDNK theory is capable of truly describing the stable and
causal dynamics of the system in the long wavelength limit.

Apart form some numerical works (including [22, 28]), little is known about
BDNK theory in the ‘nonlinear’ regime. The second part of this work is devoted
to the analysis of the correlation functions beyond the linear response
regime. The main question is how the interactions between hydrodynamic modes
affects the correlation functions. In a systematic approach, this can be
explored within the framework of the effective field theory of hydrodynamics.
However, we choose a different route here. We show that the correlation
function of certain operators in the nonlinear regime factorize if the
distribution of fluctuations is taken to be Gaussian. In quantum field theory
(QFT), the factorized correlation function is in the form of a Feynman
integral, representing a simple loop diagram with two lines, i.e., a bubble
diagram. We then follow the famous Landau method to find the singularities of
the Feynman integral. Instead of computing the integral explicitly, Landau
proposed a set of equations, the so-called “Landau equations”, which are
solved to specify the threshold singularities of the Feynman integral.

Keeping in mind the ideas discussed in the previous paragraph, we explore the
analytical structure of the shear stress response function in § 4. We derive
the related Landau equations and solve them analytically. In addition, we
develop a method to numerically find the location of branch cuts. We show that
our results are in complete agreement with recent effective field theory (EFT)
results obtained from explicit loop calculations.

Finally, in § 5 we end with a review of the results, mentioning possible
applications and discussing some followup directions.

Note: Throughout this paper we will refer to “Landau” in two contexts:

  1. 1.

“Landau-Lifshitz frame” (or Landau-Lifshitz hydrodynamics), which refers to
some specific way of defining the fluid four-velocity beyond zeroth order in
the derivative expansion.

  2. 2.

“Landau equations/conditions” refer to those equations that allow finding the
singularity of the Feynman integral without explicitly performing the
integration.

It should be emphasized that the above two items have nothing to do with each
other.

##  2 First-order stable conformal hydrodynamics

In an uncharged conformal fluid, the constitutive relations in the most
general frame are determined by four dimensionless numbers $\bar{p}$,
$\bar{\pi}$, $\bar{\theta}$ and $\bar{\eta}$ [11]. Specializing to $D=4$
dimensions, one writes

$\begin{split}T^{\mu\nu}=&\,T^{3}\left(\bar{p}\,T+\,\bar{\pi}\frac{u^{\lambda}\partial_{\lambda}T}{T}+\,\bar{\pi}\frac{\partial_{\lambda}u^{\lambda}}{3}\right)(\eta^{\mu\nu}+\,4\,u^{\mu}\,u^{\nu})\\\
&+\bar{\theta}\,T^{3}\left[\left(u^{\lambda}\partial_{\lambda}u^{\mu}+\,\frac{P^{\mu\lambda}\partial_{\lambda}T}{T}\right)u^{\nu}+\,(\mu\leftrightarrow\nu)\right]-\,\bar{\eta}\,T^{3}\,\sigma^{\mu\nu}\,.\end{split}$
(2.1)

Note that $p=\bar{p}T^{4}$ is the equilibrium pressure and
$\eta=\bar{\eta}T^{3}$ is the shear viscosity, both depending of the
microscopic of the fluid. Two new transport coefficients, $\pi=\bar{\pi}T^{3}$
and $\theta=\bar{\theta}T^{3}$ are responsible for regulating the theory in
the UV limit.

The stability and causality of the first order hydrodynamics has been
discussed extensively in the literature, over the past few years [10, 11, 27,
28, 12, 19, 31, 29, 30, 13, 17, 18, 25, 21, 16, 26, 24, 23, 22, 20, 14, 15].
In the case of conformal fluid in $D=4$ dimensions, these conditions are given
by [10]

$1-\frac{3\bar{\eta}}{\bar{\theta}}-\frac{\bar{\eta}}{\bar{\pi}}>0\,,\,\,\,\,\,\bar{\pi}>4\bar{\eta}\,.$
(2.2)

To satisfy (2.2), it is sufficient to take $\bar{\theta}>4\bar{\eta}$ and
$\bar{\pi}>4\bar{\eta}$ [11]. These conditions ensure that in any equilibrium
state specified by $u^{\mu}=\gamma(1,\textbf{v}_{0})$ and $T=T_{0}$,
linearised evolution of (2.1) is always dissipative and in the light-cone. Let
us denote that the relativistic hydrodynamics in a general frame, i.e., (2.1),
together with the conditions (2.2) is usually referred to as the “stable
first-order relativistic hydrodynamics” theory.

##  3 Linear response

As mentioned in the Introduction, the focus of this work is on correlation
functions. In this section, we explore the correlation functions in the linear
regime. To this end, we take $T^{\mu\nu}$ from (2.1) and linearize the
hydrodynamic equations, namely $\partial_{\mu}T^{\mu\nu}$=0, around the
equilibrium state specified by $u^{\mu}=\gamma(1,\textbf{v}_{0})$ and
$T=\text{constant}$. The linearised equations then govern the coupled linear
evolution of $\delta\textbf{v}$ and $\delta T$, as follows.

In conventional hydrodynamic treatments, taking $\textbf{v}_{0}\neq 0$ is
especially important to make the instabilities of the equilibrium manifest
[4]. However, in this work, we restrict our study to the linearised equations
when $\textbf{v}_{0}=0$.

The conservation of energy, $\partial_{\mu}T^{\mu t}=0$ gives

$\underline{\partial_{t}\frac{\delta
T}{T}+\frac{1}{3}\boldsymbol{\nabla}\cdot\delta\textbf{v}}+\frac{1}{12\bar{p}T}\left[\bar{\theta}\,\boldsymbol{\nabla}\cdot\left(\underline{\underline{\partial_{t}\delta\textbf{v}+\boldsymbol{\nabla}\frac{\delta
T}{T}}}\right)+3\bar{\pi}\,\partial_{t}\left(\underline{\partial_{t}\frac{\delta
T}{T}+\frac{1}{3}\boldsymbol{\nabla}\cdot\delta\textbf{v}}\right)\right]=\,0$
(3.1)

The equations of momentum conservation, $\partial_{\mu}T^{\mu i}=0$, give

$\underline{\underline{\partial_{t}\delta\textbf{v}+\boldsymbol{\nabla}\frac{\delta
T}{T}}}+\frac{1}{4\bar{p}T}\left[\bar{\theta}\partial_{t}\left(\underline{\underline{\partial_{t}\delta\textbf{v}+\boldsymbol{\nabla}\frac{\delta
T}{T}}}\right)+\bar{\pi}\boldsymbol{\nabla}\left(\underline{\partial_{t}\frac{\delta
T}{T}+\frac{1}{3}\boldsymbol{\nabla}\cdot\delta\textbf{v}}\right)-\frac{4}{3}\bar{\eta}\boldsymbol{\nabla}^{2}\delta\textbf{v}\right]=\,0$
(3.2)

The first two terms in each equation come from $T^{\mu\nu}$ of order zero,
while the terms in brackets are associated with $T^{\mu\nu}$ of order one. It
is well known, and evident in the above equations, that the stable theory of
hydrodynamics is constructed in such a way that the zeroth order equation
emerges “with new transport coefficients $\bar{\theta}$ and $\bar{\pi}$” at
the first order. If one proceeds to solve the equations perturbatively, that
is, order by order in the derivative expansion, the above two coefficients
will be washed out at the first order simply because the zeroth order equation
is used. However, this differs from the concept of “stable relativistic
hydrodynamics”; instead, the equations must be solved non-perturbatively.

###  3.1 Correlation functions in longitudinal channel

An important aspect of the “stable relativistic hydrodynamics” is that it
allows us to treat fluid velocity and temperature as the hydrodynamic
variables, employing the same concepts we know in thermodynamics; and similar
to what is done in non-relativistic fluid dynamics. Then this suggests that in
addition to the correlation functions of conserved charges (or currents), the
correlation function of hydrodynamic variables, such as $\langle
v_{i}(t_{1},\textbf{x}_{1})v_{j}(t_{2},\textbf{x}_{2})\rangle$, may also have
sensible description. It should be emphasized that our main physical arguments
in this paper will be made based on the correlation function of conserved
densities. However, as we will explain below, each of the latter correlators
is a linear combination of former ones.

keeping in mind the above discussion, we first compute the correlation
function of the hydrodynamic variables. Considering
$\delta\phi_{a}=(\frac{\delta T}{T},\delta v_{i}^{\perp},\delta
v^{\parallel})$, let us define

$\langle\delta\phi_{a}^{\omega\textbf{k}}\delta\phi_{b}^{\omega^{\prime}\textbf{k}^{\prime}}\rangle=\,(2\pi)^{4}\delta(\omega+\omega^{\prime})\delta(\textbf{k}+\textbf{k}^{\prime})\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}_{\omega\textbf{k}}\,.$
(3.3)

Here
$\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}_{\omega\textbf{k}}$ is
the correlation function of $\phi_{a}$ and $\phi_{b}$:

$\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}_{\omega\textbf{k}}=\int_{-\infty}^{+\infty}dt\int
d^{3}\text{x}\,e^{i\omega
t}e^{-i\textbf{k}\cdot\textbf{x}}\langle\delta\phi_{a}(t,\textbf{x})\delta\phi_{b}(0,\textbf{0})\rangle\,.$
(3.4)

We also find it useful to define dimensionless quantities as

${\mathfrak{w}}=\frac{\bar{\eta}}{\bar{p}}\,\frac{\omega}{T}\,,\,\,\,\,\,\,{\mathfrak{q}}=\frac{\bar{\eta}}{\bar{p}}\,\frac{\textbf{k}}{T}\,,\,\,\,\,\,\,\pi_{\eta}=\frac{\bar{\pi}}{\bar{\eta}}\,,\,\,\,\,\,\,\theta_{\eta}=\frac{\bar{\theta}}{\bar{\eta}}\,.$
(3.5)

Then we follow the standard textbook method to calculate various correlation
functions
$\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}$,
for the longitudinal and transverse channels; see Appendix B for details.
Interestingly we find that 666Similar to (3.6), one can easily compute
$\big{\langle}\frac{\delta T}{T}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}$ and $\big{\langle}\delta
v_{i}\delta v_{j}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}$. More details
can be found in Appendix B.

$\begin{split}\big{\langle}\delta v_{i}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}=&\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}({\mathfrak{w}},{\mathfrak{q}})}\,iq_{i}\,(i\,\boxed{\pi_{\eta}}\,{\mathfrak{w}}-4)\,\big{(}4-i(\theta_{\eta}+\pi_{\eta}){\mathfrak{w}}\big{)}\\\
&+\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}^{*}({\mathfrak{w}},{\mathfrak{q}})}\,(-iq_{i})\,(-i\,\boxed{\theta_{\eta}}\,{\mathfrak{w}}-4)\,\big{(}4+i(\theta_{\eta}+\pi_{\eta}){\mathfrak{w}}\big{)}\,,\\\
\big{\langle}\frac{\delta T}{T}\delta
v_{i}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}=&\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}({\mathfrak{w}},{\mathfrak{q}})}\,iq_{i}\,(i\,\boxed{\theta_{\eta}}\,{\mathfrak{w}}-4)\,\big{(}4-i(\theta_{\eta}+\pi_{\eta}){\mathfrak{w}}\big{)}\\\
&+\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}^{*}({\mathfrak{w}},{\mathfrak{q}})}\,(-iq_{i})\,(-i\,\boxed{\pi_{\eta}}\,{\mathfrak{w}}-4)\,\big{(}4+i(\theta_{\eta}+\pi_{\eta}){\mathfrak{w}}\big{)}]\,.\end{split}$
(3.6)

Regardless of whether the stability and causality conditions are satisfied, it
is clear from (3.6) that for general values of $\theta_{\eta}$ and
$\pi_{\eta}$:

$\big{\langle}\delta v_{i}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\neq\big{\langle}\frac{\delta
T}{T}\delta v_{i}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\,.$ (3.7)

If these correlators are treated as microscopic correlation functions, then
“time-reversal invariance” requires that they must be equal. The latter leads
to

$\pi_{\eta}=\,\theta_{\eta}\,.$ (3.8)

This constraint simply tells us that the time-reversal-invariant stable frames
occupy a zero-measure subspace in the entire space of stable frames identified
by (2.2), i.e., like a curve lying on the plane. However, this is not the
case. If we implement derivative expansion, the above two correlation
functions are equal to first order in the derivative, independent of the
values of $\pi_{\eta}$ and $\theta_{\eta}$. In other words, the inequality
(3.7) occurs outside the domain where only first-order dissipative corrections
are sufficient to describe the dynamics. Below, for convenience, we restrict
the computation to frames satisfying (3.8), whether stable or not.

Now let us move on to calculating the energy density correlation function,
namely $G_{T_{tt}T_{tt}}({\mathfrak{w}},{\mathfrak{q}})=\langle
T^{tt}T^{tt}\rangle_{{\mathfrak{w}}{\mathfrak{q}}}$. Using (2.1), one writes

$\begin{split}G_{T_{tt}T_{tt}}({\mathfrak{w}},{\mathfrak{q}})=\bar{p}^{2}T^{8}&\bigg{[}9(\pi_{\eta}^{2}{\mathfrak{w}}^{2}+16)\big{\langle}\frac{\delta
T}{T}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}+\pi_{\eta}^{2}\,q_{i}q_{j}\big{\langle}\delta
v_{i}\delta v_{j}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\\\
\,\,\,\,\,\,\,\,\,\,\,&\,\,\,\,\,-3\pi_{\eta}q_{j}\bigg{(}\big{(}\pi_{\eta}{\mathfrak{w}}+4i\big{)}\big{\langle}\frac{\delta
T}{T}\delta
v_{j}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}+\big{(}\pi_{\eta}{\mathfrak{w}}-4i\big{)}\big{\langle}\delta
v_{j}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\bigg{)}\bigg{]}\,.\end{split}$
(3.9)

Then by using the corresponding correlation functions of the hydrodynamic
variables (Appendix B), we find:

$\begin{split}G_{T_{tt},T_{tt}}=\frac{-6\bar{\eta}T^{4}{\mathfrak{q}}^{2}}{|\mathcal{D}_{L}({\mathfrak{w}},{\mathfrak{q}})|^{2}}&\bigg{[}+\left(\pi_{\eta}-4\right)\pi_{\eta}^{2}\theta_{\eta}^{2}\,{\mathfrak{q}}^{6}+96\left(\pi_{\eta}^{2}-4\pi_{\eta}+8\right)\theta_{\eta}\,{\mathfrak{q}}^{4}+2304\left(\pi_{\eta}-4\right){\mathfrak{q}}^{2}\\\
&\,\,\,\,\,\,+3{\mathfrak{w}}^{2}\bigg{(}-2304\pi_{\eta}+\pi_{\eta}\theta_{\eta}\left(4\theta_{\eta}^{2}-3\left(\pi_{\eta}-4\right)\pi_{\eta}\theta_{\eta}-8\left(\pi_{\eta}-2\right)\pi_{\eta}\right){\mathfrak{q}}^{4}\\\
&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+48\left(-2\left(\pi_{\eta}-2\right)\pi_{\eta}\theta_{\eta}+\left(\pi_{\eta}-4\right)\theta_{\eta}^{2}+\left(\pi_{\eta}-8\right)\pi_{\eta}^{2}\right){\mathfrak{q}}^{2}\bigg{)}\\\
\,+9{\mathfrak{w}}^{4}\big{(}4\pi_{\eta}^{3}&\theta_{\eta}{\mathfrak{q}}^{2}+\pi_{\eta}\theta_{\eta}^{2}\left(\pi_{\eta}\left(3\pi_{\eta}-4\right){\mathfrak{q}}^{2}-48\right)-4\pi_{\eta}^{3}\left(\pi_{\eta}{\mathfrak{q}}^{2}+12\right)\big{)}-27{\mathfrak{w}}^{6}\,\pi_{\eta}^{3}\theta_{\eta}^{2}\bigg{]}\,.\end{split}$
(3.10)

In figure 1, this correlation function has been plotted for several different
cases. The red curve represents the familiar result in the Landau-Lifshitz
frame. The two peaks correspond to two sound modes, the width of which is
proportional to the shear viscosity. The green and blue curves correspond to
hydrodynamics in the general frame; the first one in the unstable frame while
the second is in the stable frame (by means of (2.2)).

Figure 1: Energy density correlation function in the Landau-Lifshitz frame
(red), in an unstable frame (green), and in a stable frame (blue).

As it is seen, the correlation function develops a range of negative values in
both the unstable and stable frames. The situation shown in the figure is for
${\mathfrak{q}}=1.5$. The smaller the value of ${\mathfrak{q}}^{2}$, the
narrower the range of negative values around ${\mathfrak{w}}=0$. Close to
${\mathfrak{w}}=0$ we have

$\lim_{{\mathfrak{w}}\rightarrow
0}G_{T_{tt}T_{tt}}({\mathfrak{w}},{\mathfrak{q}})=-\frac{6\,\bar{\eta}\,T^{4}\,\big{(}\pi_{\eta}^{2}\theta_{\eta}\,{\mathfrak{q}}^{2}+48(\pi_{\eta}-4)\big{)}}{\theta_{\eta}(\pi_{\eta}-4)\,{\mathfrak{q}}^{2}+48}\,.$
(3.11)

As it is clearly seen, in a stable frame ($\theta_{\eta}\,,\pi_{\eta}>4$)
777The sufficient conditions for stability and causality, discussed earlier,
take a simpler form when expressed in terms of dimensionless coefficients:
$\pi_{\eta}>4$ and $\theta_{\eta}>4$., (3.11) is negative for any value of
${\mathfrak{q}}^{2}$ 888However, for unstable frames
($\theta_{\eta}=\pi_{\eta}<4$), this negative range disappears at sufficiently
small values of ${\mathfrak{q}}^{2}$: $\lim_{{\mathfrak{q}}\rightarrow
0}\lim_{{\mathfrak{w}}\rightarrow
0}G_{T_{tt}T_{tt}}({\mathfrak{w}},{\mathfrak{q}})=6\bar{\eta}T^{4}(4-\pi_{\eta})$..
Of course, $G_{T_{tt}T_{tt}}$ has such a negative value interval, it does not
mean that the correlation function of a hermitian operator can become
negative, but that the relativistic hydrodynamics in the general frame has
pathological characteristics that leads to this result999In [22], another kind
of pathology in BDNK theory was observed in the evolution of $\langle
T^{tt}(t,\textbf{x})\rangle$ in plane-symmetric configurations of a conformal
fluid in Minkowski spacetime..

Implementing the derivative expansion cures this pathology. To see this, let
us consider equations (3.1) and (3.2) at the zeroth order of the hydrodynamic
derivative expansion

$\partial_{t}\frac{\delta
T}{T}+\frac{1}{3}\boldsymbol{\nabla}\cdot\delta\textbf{v}=\mathcal{O}(\partial^{2})\,,\,\,\,\,\,\,\,\,\,\partial_{t}\delta\textbf{v}+\boldsymbol{\nabla}\frac{\delta
T}{T}=\mathcal{O}(\partial^{2})\,.$ (3.12)

By substituting these equations into the first order equations (3.1) and
(3.2), the overall effect of general frames, namely the $\bar{\theta}$ and
$\bar{\pi}$ terms, is washed out. Since this is equivalent to take
$\theta_{\eta}=\pi_{\eta}=0$, it is then clear that the results associated
with unstable and stable frames will reduce to the Landau-Lifshitz frame
result. It should be noted that this feature is actually specific to conformal
fluids: implementing derivative expansion together with using zeroth order
equations make the ‘general-frame’ hydrodynamic equations in a conformal fluid
fall into the Landau-Lifshitz hydrodynamics 101010In a non-conformal fluid,
however, this is not the case (see Appendix A of [11])..

###  3.2 Correlation functions in transverse channel

Among various correlation functions associated with this channel, we choose to
explicitly express
$G_{T_{ti_{\perp}}T_{ti_{\perp}}}({\mathfrak{w}},{\mathfrak{q}})=\langle
T^{ti_{\perp}}T^{ti_{\perp}}\rangle_{{\mathfrak{w}}{\mathfrak{q}}}$, $i=x,y$.
To this end, we first use the energy momentum tensor at linear order and write

$G_{T_{ti_{\perp}}T_{tj_{\perp}}}({\mathfrak{w}},{\mathfrak{q}})=\delta_{ij}\,\bar{p}^{2}T^{8}\gamma^{2}\bigg{[}\bigg{(}\gamma(\textbf{v}_{0}^{2}-\theta_{\eta}){\mathfrak{w}}+\gamma(\theta_{\eta}-1){\mathfrak{q}}\cdot\textbf{v}_{0}\bigg{)}^{2}+16\bigg{]}\langle\delta
v_{i_{\perp}}\delta v_{j_{\perp}}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}$
(3.13)

Note that we have set the equilibrium fluid velocity to be
$u^{\mu}=(1,\textbf{v}_{0})$. Then by utilizing the correlation function of
the transverse components of the velocity, found in Appendix B, we arrive at

$\begin{split}G_{T_{ti_{\perp}},T_{tj_{\perp}}}=\,\delta_{ij}\frac{-2\bar{\eta}T^{4}}{|\mathcal{D}_{T}({\mathfrak{w}},{\mathfrak{q}})|^{2}}&\bigg{[}\bigg{(}\big{(}\gamma^{2}(\theta_{\eta}-1)+1\big{)}{\mathfrak{w}}-\gamma(\theta_{\eta}-1){\mathfrak{q}}\cdot\textbf{v}_{0}\bigg{)}^{2}+16\bigg{]}\\\
&\,\,\,\,\,\,\,\times\bigg{[}\big{(}\gamma^{2}(\theta_{\eta}-1)-\theta_{\eta}\big{)}{\mathfrak{q}}^{2}+\big{(}1-\gamma^{2}(\theta_{\eta}-1)\big{)}({\mathfrak{q}}\cdot\textbf{v}_{0})\,{\mathfrak{w}}\bigg{]}\end{split}$
(3.14)

where $\mathcal{D}_{T}$ is the spectral function of the transverse (shear)
channel (given by (B.8)). We find that

  1. 1.

In Landau-Lifshitz frame, $G_{T_{ti_{\perp}},T_{tj_{\perp}}}$ is positive
definite only when $\textbf{v}_{0}=0$. The non-positivity of this correlation
function at $\textbf{v}_{0}\neq 0$ may correspond to the existence of a gapped
mode in the upper half of the complex frequency plane, and thus to an
instability of the equilibrium. Approval of this statement requires more
information, especially from sound channel. We leave more investigation of the
issue to future work.

  2. 2.

In a general frame, we find that the above correspondence fails to exist. For
both stable and unstable frames, the correlation function develops a range of
negativity. This is similar to our previous result in the sound channel.

###  3.3 Kubo Formulas

Having calculated the correlation functions, we can now derive the
corresponding Kubo formulas for the three transport coefficients $\bar{\eta}$,
$\bar{\pi}$, and $\bar{\theta}$. As mentioned in [11], in conformal fluids
$\bar{\pi}T^{3}$ and $\bar{\theta}T^{3}$ can be considered as the relaxation
time for energy density and momentum density, respectively. This then suggests
that, as an example, $\bar{\pi}$ can be derived from (3.10). We find

$\bar{\pi}=\,\frac{1}{2T^{4}}\,\lim_{\textbf{k}\rightarrow
0}\,\frac{\omega^{2}}{\textbf{k}^{2}}G_{T_{tt},T_{tt}}(\omega,\textbf{k})\,.$
(3.15)

To find $\bar{\theta}$, we should compute the correlation function of the
energy density-momentum density. In terms of dimensionless quantities defined
in (3.5), we find

$\begin{split}G_{T_{tt}T_{ti}}({\mathfrak{w}},{\mathfrak{q}})=\bar{p}^{2}T^{8}&\bigg{[}3\theta_{\eta}q_{i}\big{(}\pi_{\eta}{\mathfrak{w}}+4i\big{)}\big{\langle}\frac{\delta
T}{T}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}+\big{(}\pi_{\eta}{\mathfrak{w}}+4i\big{)}\,q_{j}\big{\langle}\delta
v_{i}\delta v_{j}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\\\
\,\,\,\,\,\,\,\,\,\,\,&\,\,\,\,\,-\bigg{(}3\big{(}\pi_{\eta}^{2}{\mathfrak{w}}^{2}+16\big{)}\big{\langle}\frac{\delta
T}{T}\delta
v_{i}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}+\pi_{\eta}\theta_{\eta}{\mathfrak{q}}^{2}\big{\langle}\delta
v_{i}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\bigg{)}\bigg{]}\,.\end{split}$
(3.16)

We can then use the correlation function for hydrodynamic variables (see
Appendix B) to find an explicit expression for the above correlation function.
Doing so, we find

$\bar{\theta}=\,-\bar{\pi}-\frac{1}{T^{4}}\,\lim_{{\mathfrak{w}}\rightarrow
0}\lim_{\textbf{k}\rightarrow
0}\frac{\omega\,k_{i}}{\textbf{k}^{2}}G_{T_{tt},T_{ti}}(\omega,\textbf{k})\,.$
(3.17)

Finally, $\bar{\eta}$ can be found from the momentum density correlation
function. In terms of the dimensionless quantities defined in (3.5), we have

$\begin{split}\frac{q_{i}q_{j}}{{\mathfrak{q}}^{2}}G_{T_{ti}T_{tj}}({\mathfrak{w}},{\mathfrak{q}})=\bar{p}^{2}T^{8}&\bigg{[}\theta_{\eta}^{2}{\mathfrak{q}}^{2}\big{\langle}\frac{\delta
T}{T}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}+\big{(}\pi_{\eta}^{2}{\mathfrak{w}}^{2}+16\big{)}\frac{q_{i}q_{j}}{{\mathfrak{q}}^{2}}\big{\langle}\delta
v_{i}\delta v_{j}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\\\
\,\,\,\,\,\,\,\,\,\,\,&\,\,\,\,\,+\bigg{(}\big{(}4i-\theta_{\eta}{\mathfrak{w}}\big{)}q_{i}\big{\langle}\frac{\delta
T}{T}\delta
v_{i}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}-\big{(}4i+\theta_{\eta}{\mathfrak{w}}\big{)}q_{i}\big{\langle}\delta
v_{i}\frac{\delta
T}{T}\big{\rangle}_{{\mathfrak{w}}{\mathfrak{q}}}\bigg{)}\bigg{]}\,.\end{split}$
(3.18)

Again, we use correlation functions of hydrodynamic variables to find explicit
expressions for the above correlation function. From the latter we find

$\bar{\eta}=\,-\frac{1}{4}\,\bar{\theta}+\frac{3}{8T^{4}}\,\lim_{\textbf{k}\rightarrow
0}\frac{\omega^{2}k_{i}k_{j}}{\textbf{k}^{4}}G_{T_{ti},T_{tj}}(\omega,\textbf{k})\,.$
(3.19)

The above Kubo formulas were found from the correlation functions of the
conserved charges in the sound channel. One can also find Kubo formulas from
the “current” correlators, in both sound and shear channels.

Equations (3.15), (3.17), and (3.19) are the first results of the Kubo formula
for the transport coefficients in BDNK theory. In a given system, the
correlation functions involved in the above Kubo formulas may be extracted
from some experimental data. These formulas then open a window to examine
whether ‘stable hydrodynamics’ can truly describe dynamics of the system in
the long wavelength limit. This is evident if $\bar{\pi}$ and $\bar{\theta}$
are found to be non-zero and satisfy (2.2).

##  4 Correlation functions beyond linear response

Hydrodynamics is essentially a nonlinear theory. Nonlinearities manifest as
the interactions between hydrodynamic modes. A systematic approach for
accounting these interactions is to construct a hydrodynamic effective field
theory [35, 36, 37, 38]. Then from this effective field theory, the corrected
correlation function can be found [39], resulting in a large renormalization
of the transport coefficients [40] as well as the long time tails [41].

These effects can also be explained by some general hydrodynamic statements
[42]. Let us recall that hydrodynamic variables are macroscopically averaged
values of densities over regions of size $b$, where $b\gg\ell_{mic}$.
Typically, the microscopic length scale $\ell_{mic}$ is of order of the
correlation length: $\ell_{mic}\sim\xi$. In other words, a hydrodynamic
fluctuation, such as $\delta v$, can be regarded as averaging over
$\big{(}b/\xi\big{)}^{3}$ independent correlation volumes. According to the
central limit theorem, the distribution of such average will approach a
Gaussian as $\big{(}b/\xi\big{)}^{3}\gg 1$. The distribution is Gaussian, the
non-Gaussian contribution is suppressed by an additional factor of
$\big{(}b/\xi\big{)}^{-3}$, and knowledge of the two-point function at the
level of the linear response is sufficient to find correlation functions
beyond the linear response [41].

Following the above logic, we would like to investigate in this section the
structure of the correlation functions beyond the linear response. Instead of
constructing the effective field theory, we assume that the distribution of
fluctuations is Gaussian.

###  4.1 Correlation functions of the shear stress

To reduce technical complexity, we restrict this section to the analysis of
the correlation function of the shear stress,
$G_{T_{xy}T_{xy}}=\big{\langle}T_{xy}T_{xy}\big{\rangle}$. The reason for this
selection is as follows. At linear order, the shear stress $T_{xy}$ vanishes;
it just starts to contribute in quadratic order:

$T_{xy}=\,4\bar{p}T^{4}\,\delta v_{x}\delta v_{y}\,;$ (4.1)

so the correlation function
$G_{T_{xy}T_{xy}}=\big{\langle}T_{xy}T_{xy}\big{\rangle}$ is given by (n.l.
stands for nonlinear)

$G^{(\text{n.l.})}_{T_{xy}T_{xy}}(t,\textbf{x})=\,w^{2}\,\big{\langle}\delta
v_{x}(t,\textbf{x})\delta v_{y}(t,\textbf{x})\delta v_{x}(0,\textbf{0})\delta
v_{y}(0,\textbf{0})\big{\rangle}\,,$ (4.2)

with $w=4\bar{p}T^{4}$. Since the distribution of fluctuations is assumed to
be Gaussian, the above correlation can be factorized as

$G^{(\text{n.l.})}_{T_{xy}T_{xy}}(\omega,\textbf{k})=w^{2}\int\frac{d\omega^{\prime}}{2\pi}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}\,G_{\delta
v_{x}\delta v_{x}}(\omega^{\prime},\textbf{k}^{\prime})G_{\delta v_{y}\delta
v_{y}}(\omega-\omega^{\prime},\textbf{k}-\textbf{k}^{\prime})\,.$ (4.3)

The two correlation functions in the integrand are those already calculated in
the linear response regime. It is now clear why we chose to use shear stress:
this is actually a special case, and to study it the only information about
the linear response we need is the correlation function of velocity in the
‘transverse channel’ (see Appendix B).

Taking $p=(\omega,\textbf{k})$, the integrand can be written as

$\bigg{(}\underbrace{G^{(+)}_{\delta v_{x}\delta
v_{x}}(p^{\prime})}_{①}+\underbrace{G^{(+)}_{\delta v_{x}\delta
v_{x}}(-p^{\prime})}_{②}\bigg{)}\bigg{(}\underbrace{G^{(+)}_{\delta
v_{y}\delta v_{y}}(p-p^{\prime})}_{③}+\underbrace{G^{(+)}_{\delta v_{y}\delta
v_{y}}(-p+p^{\prime})}_{④}\bigg{)}$ (4.4)

where $(+)$ represents the one-sided Fourier transformation (see Appendix B
for details). Of the four possible multiplicative terms, $①\times④$ and
$②\times③$ do not contribute to the integral. The reason is simply that, for
each of them, all pole singularities lie in only half of the $\omega^{\prime}$
complex plane; either in the upper half or in the lower half. Regarding the
remaining two terms, $②\times④$ only contributes to the upper half plane
structure of $G^{(\text{n.l.})}_{T_{xy}T_{xy}}(p)$, the advanced Green’s
function. And finally this is $①\times③$ term which contributes to the lower
half plane structure of $G^{(\text{n.l.})}_{T_{xy}T_{xy}}(p)$. Therefore, to
find the analytic structure of the “response function” of $T_{xy}$, i.e.,
$G^{R}_{T_{xy}T_{xy}}$, it is sufficient to consider the following integral

$\begin{split}\bar{w}^{2}&\int\frac{d\omega^{\prime}}{2\pi}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}\,G^{(+)}_{\delta
v_{x}\delta v_{x}}(\omega^{\prime},\textbf{k}^{\prime})G^{(+)}_{\delta
v_{y}\delta v_{y}}(\omega-\omega^{\prime},\textbf{k}-\textbf{k}^{\prime})\\\
&\sim\int\frac{d{\mathfrak{w}}^{\prime}}{2\pi}\int\frac{d^{3}q^{\prime}}{(2\pi)^{3}}\,\frac{\mathcal{N({\mathfrak{w}}^{\prime})}}{\mathcal{D}_{T}({\mathfrak{w}}^{\prime},{\mathfrak{q}}^{\prime})}\frac{\mathcal{N({\mathfrak{w}}-{\mathfrak{w}}^{\prime})}}{\mathcal{D}_{T}({\mathfrak{w}}-{\mathfrak{w}}^{\prime},{\mathfrak{q}}-{\mathfrak{q}}^{\prime})}\,.\end{split}$
(4.5)

The expression $\mathcal{D}_{T}$ in the denominator is the spectral function
of the transverse channel:

$\mathcal{D}_{T}\big{(}{\mathfrak{w}},{\mathfrak{q}},\textbf{v}_{0}=\textbf{0}\big{)}=\theta_{\eta}{\mathfrak{w}}^{2}+4i{\mathfrak{w}}-{\mathfrak{q}}^{2}=0$
(4.6)

and encodes the dispersion relations of the two gapped modes associated with
this channel

${\mathfrak{w}}_{1,2}=\,-\frac{i}{\theta_{\eta}}\big{(}2\mp\sqrt{4-\theta_{\eta}{\mathfrak{q}}^{2}}\big{)}\,.$
(4.7)

Note that we are considering fluctuations on top of an equilibrium state with
$u^{\mu}=(1,\textbf{0})$.

Needless to say, the integral (4.5) is a kind of one-loop integral that we
usually deal with in quantum field theory:

$\big{(}{\mathfrak{w}},{\mathfrak{q}}\big{)}$$\big{(}{\mathfrak{w}}^{\prime},{\mathfrak{q}}^{\prime}\big{)}$

It is convenient to use Feynman parameters to rewrite (4.5) as

$\mathcal{I}({\mathfrak{w}},{\mathfrak{q}})=\int\frac{d{\mathfrak{w}}^{\prime}}{2\pi}\int\frac{d^{3}q^{\prime}}{(2\pi)^{3}}\,\int_{0}^{1}d\alpha_{1}\int_{0}^{1}d\alpha_{2}\,\delta(\alpha_{1}+\alpha_{2}-1)\,\frac{\mathcal{N({\mathfrak{w}}^{\prime})}N({\mathfrak{w}}-{\mathfrak{w}}^{\prime})}{\big{[}\alpha_{1}\mathcal{D}_{T}({\mathfrak{w}}^{\prime},{\mathfrak{q}}^{\prime})+\alpha_{2}\mathcal{D}_{T}({\mathfrak{w}}-{\mathfrak{w}}^{\prime},{\mathfrak{q}}-{\mathfrak{q}}^{\prime})\big{]}^{2}}$
(4.8)

Then based on general complex analysis arguments, it is possible to discuss
the analytic properties of $\mathcal{I}({\mathfrak{w}},{\mathfrak{q}})$,
without performing the frequency/momentum integral. This is basically the idea
proposed by Landau in the famous paper [44]. Here, our goal is to use the
Landau equations/conditions derived in [44] to find poles and branch cuts of
$\mathcal{I}({\mathfrak{w}},{\mathfrak{q}})$ in (4.8), which are actually the
singularities of $G^{R}_{T_{xy}T_{xy}}$ as well.

###  4.2 Solving the Landau equations

Applying Landau conditions (see Appendix C for a brief review) to
$\mathcal{I}({\mathfrak{w}},{\mathfrak{q}})$, one can find the (threshold)
singularities of $G^{R}_{T_{xy}T_{xy}}$. This is what will be done in the end
of current subsection. Before this, we want to specify the branch cut
structure of the response function. Branch cuts are specified by singular
points of the integrand on the integral contour in (4.5) [45]. The latter is
equivalent to requiring

$\mathcal{D}_{T}({\mathfrak{w}}^{\prime},q^{\prime})=0\,,\,\,\,\,\,\,\,\mathcal{D}_{T}({\mathfrak{w}}-{\mathfrak{w}}^{\prime},q\pm
q^{\prime})=0\,.$ (4.9)

These are actually the on-shell conditions coming from the denominator of
(4.5). Note that we have used
$\int\frac{d^{3}q^{\prime}}{(2\pi)^{3}}=\frac{1}{(2\pi)^{2}}\int q^{\prime
2}dq^{\prime}\int_{-1}^{1}d\cos\theta^{\prime}$ and performed the integration
over $\theta^{\prime}$, which is the angle between ${\mathfrak{q}}$ and
${\mathfrak{q}}^{\prime}$.

Our goal is to solve the conditions (4.9) to find the branch cuts. Let us
assume that the entire solution is identified by the following complex
equation

$s({\mathfrak{w}},{\mathfrak{q}})=0\,.$ (4.10)

For a given ${\mathfrak{q}}$, this equation represents a set of points in the
complex ${\mathfrak{w}}$ plane: the set $\mathcal{S}_{{\mathfrak{q}}}$. In the
following, we introduce a method to solve (4.9), not for the function $s$, but
directly for the set $\mathcal{S}_{{\mathfrak{q}}}$.

Clearly, the two equations in (4.9) are the on-shell conditions for the two
lines of the loop when the internal momenta are parallel. This simply means
that for a given external $({\mathfrak{w}},{\mathfrak{q}})$, each line of the
loop carries one of the two modes (4.7) (in
${\mathfrak{q}}^{\prime}||{\mathfrak{q}}$ conditions). We then conclude that
the set $\mathcal{S}_{{\mathfrak{q}}}$ must be given by

$\mathcal{S}_{{\mathfrak{q}}}=\bigcup_{{\mathfrak{q}}^{\prime}||{\mathfrak{q}}}\,\,\bigcup_{i,j\in\\{1,2\\}}\bigg{\\{}{\mathfrak{w}}_{i}({\mathfrak{q}}^{\prime})+{\mathfrak{w}}_{j}({\mathfrak{q}}-{\mathfrak{q}}^{\prime})\bigg{\\}}\,.$
(4.11)

To make things simpler, we take the external momentum directed in one specific
direction: ${\mathfrak{q}}=(0,0,q)$. Then by reparameterizing the momenta in
the loop as

$\big{(}{\mathfrak{w}},{\mathfrak{q}}\big{)}$$\big{(}\frac{{\mathfrak{w}}}{2}+{\mathfrak{w}}^{\prime},\frac{{\mathfrak{q}}}{2}+{\mathfrak{q}}^{\prime}\big{)}$

equation (4.11) takes the following more applicable form

$\boxed{\mathcal{S}_{{\mathfrak{q}}}=\bigcup_{q^{\prime}\,\geq\,0}\,\,\bigcup_{i,j\in\\{1,2\\}}\bigg{\\{}{\mathfrak{w}}_{i}\big{(}\frac{q}{2}+q^{\prime}\big{)}+{\mathfrak{w}}_{j}\big{(}\frac{q}{2}-q^{\prime}\big{)}\bigg{\\}}}$
(4.12)

Regarding the evaluation of $\mathcal{S}_{{\mathfrak{q}}}$ above, some
comments are in order:

  * •

We perform the evaluation of $\mathcal{S}_{{\mathfrak{q}}}$ in four separate
parts, $\mathcal{S}^{ij}_{{\mathfrak{q}}};i,j\in\\{1,2\\}$, corresponding to
the various values the pair $(i,j)$ may take.

  * •

At a given $q$, for each part, e.g., $(i,j)=(1,1)$, we calculate
${\mathfrak{w}}_{i}\big{(}\frac{q}{2}+q^{\prime}\big{)}+{\mathfrak{w}}_{j}\big{(}\frac{q}{2}-q^{\prime}\big{)}$
in a large number of non-negative values of $q^{\prime}$.

  * •

It turns out that the points
${\mathfrak{w}}_{i}\big{(}\frac{q}{2}+q^{\prime}\big{)}+{\mathfrak{w}}_{j}\big{(}\frac{q}{2}-q^{\prime}\big{)}$
form a curve with a few cusp points in the complex ${\mathfrak{w}}$ plane.

  * •

Depending on $(i,j)$, the points corresponding to
${\mathfrak{w}}_{i}\big{(}\frac{q}{2}+q^{\prime}\big{)}+{\mathfrak{w}}_{j}\big{(}\frac{q}{2}-q^{\prime}\big{)}$
are accumulated around those associated with either $q^{\prime}=0$ or
$q^{\prime}\rightarrow\infty$.

As we will see below, the last two comments above point to the branch cut and
branch point singularities of $\mathcal{I}$ (or equivalently
$G^{R}_{T_{xy}T_{xy}}$), respectively.

Figure 2: Illustration of the $\mathcal{S}^{ij}_{{\mathfrak{q}}}$’s at
${\mathfrak{q}}=(0,0,q=0.5)$ in the lower half of the complex frequency plane.
In each panel, we have shown values of
${\mathfrak{w}}_{i}\big{(}\frac{q}{2}+q^{\prime}\big{)}+{\mathfrak{w}}_{j}\big{(}\frac{q}{2}-q^{\prime}\big{)}$
for 150 different $q^{\prime}$ values between $0\leq q^{\prime}<1.5$.

Figure 3: Left panel: Branch cut structure of $G^{R}_{T_{xy}T_{xy}}$, found
from solving the Landau equations. Four branch point singularities
${\mathfrak{w}}_{11}$, ${\mathfrak{w}}_{22}$, ${\mathfrak{w}}_{12}$, and
${\mathfrak{w}}_{12}$ corresponds to the four parts of
$\mathcal{S}_{{\mathfrak{q}}}$, show in figure 2. Right panel: Branch cut
structure of $G^{R}_{nn}$, found from explicit loop calculations. The right
panel plots are taken from [43].

Specializing to $\theta_{\eta}=16$ in (4.7), we have shown the four parts of
$S_{{\mathfrak{q}}}$, for ${\mathfrak{q}}=0.5$, in figure 2. In each panel,
the four comments above can be clearly seen. Each colored trajectory starts in
purple at $q^{\prime}=0$ and ends in dark red at $q^{\prime}=1.5$. We have
checked that by taking large values of $q^{\prime}$, no significant features
are added to these plots.

The set of the four parts of $\mathcal{S}_{{\mathfrak{q}}}$ depicted in figure
2 have been shown in the top left panel of figure 3,. This figure actually
shows the entire branch cut structure of $G^{R}_{T_{xy}T_{xy}}$, as the
complete solution to (4.9). As it is seen, four branch points have been
induced: ${\mathfrak{w}}_{ij}:i,j=1,2$. Each of these points actually
corresponds to producing a pair of on-shell excitations
${\mathfrak{w}}={\mathfrak{w}}_{i}({\mathfrak{q}})$ and
${\mathfrak{w}}={\mathfrak{w}}_{j}({\mathfrak{q}})$, in the loop. We show the
latter at the end of this section by solving the Landau conditions
analytically.

In the same regard, more recently, ref. [43] announced some analytical results
on the structure of branch cuts in the theory of “UV-regulated nonlinear
diffusion”. The on-shell condition in this theory is given by

$\tau^{2}\omega^{2}+i\tau\omega-{\mathfrak{q}}^{2}=0$ (4.13)

where $\tau$ is the relaxation time. A comparison between this equation and
the on-shell condition in our present case, namely (4.6), shows that the two
cases will be the same if one takes

$\boxed{\theta_{\eta}=16\,\,\,\text{in this
paper}}\,\,\,\,\,\,\,\text{vs.}\,\,\,\,\,\,\boxed{\tau=4\,\,\,\text{in ref.
\cite[cite]{[\@@bibref{Number}{Abbasi:2022aao}{}{}]}}}$ (4.14)

The (top) left panel of figure 3 has been already produced for the value of
$\theta_{\eta}$ mentioned in (4.14). In the (top) right panel of figure 3, we
have shown the branch cut structure of the response function in ref. [43], at
$\tau=4$, and for the same value of ${\mathfrak{q}}$ taken in the left panel.
Note that this reference shows the result in dimensionless frequency
${\mathfrak{w}}=\tau\omega$. Therefore, for comparison with the left panel,
all numbers on both axes must be divided by 4. Doing so, we find that the two
panels are in perfect agreement with each other.

By increasing the value of ${\mathfrak{q}}$, ${\mathfrak{w}}_{11}$ and
${\mathfrak{w}}_{22}$ move closer to each other and collide at
${\mathfrak{q}}=1$. As the second case to illustrate, we have considered
${\mathfrak{q}}=1.25$. Repeating the procedure used at ${\mathfrak{q}}=0.5$,
we solve the Landau equations numerically by (4.12). The corresponding branch
cut singularities are given in the bottom left panel of figure 3. Again, we
see that the results are exactly in line with the most recent results of ref.
[43] in the theory of UV-regulated nonlinear diffusion, shown in the bottom
right panel.

Let us emphasize that the left and right panels of figure 3 have been found
through two different methods; the left figures were obtained by developing a
new numerical method to solve the Landau equations. However, the right panel
results are the result of explicit loop calculations in ref. [43]. The branch-
cut structure shown in the left panel of Figure 3 is the central result of
this section.

Let us now proceed to solve the Landau equation in order to analytically find
the threshold singularity of the response function, i.e.
${\mathfrak{w}}_{ij}:i,j\in\\{1,2\\}$ (see figure 3). Following the discussion
in Appendix C, we rewrite the denominator of the integrand in (4.8) as

$f\equiv\big{[}\alpha\,\mathcal{D}({\mathfrak{w}}^{\prime}\,,{\mathfrak{q}}^{\prime})+(1-\alpha)\,\mathcal{D}({\mathfrak{w}}-{\mathfrak{w}}^{\prime}\,,{\mathfrak{q}}-{\mathfrak{q}}^{\prime})\big{]}^{2}\,.$
(4.15)

In order for the singularity of the integrand to be a pinch singularity on the
contour integral, one needs (see (C.4) in Appendix C),

$\begin{split}\frac{\partial
f}{\partial{\mathfrak{w}}^{\prime}}=0&\,\,\,\,\,\rightarrow\,\,\,\,\,\,{\mathfrak{w}}^{\prime}=(1-\alpha){\mathfrak{w}}+\frac{2i(1-2\alpha)}{\theta}\,,\\\
\frac{\partial
f}{\partial{\mathfrak{q}}^{\prime}}=0&\,\,\,\,\,\rightarrow\,\,\,\,\,\,{\mathfrak{q}}^{\prime}=(1-\alpha){\mathfrak{q}}\,.\end{split}$
(4.16)

Substituting into the on-shell conditions, we arrive at

$\displaystyle f\big{|}_{\alpha=1}=0$
$\displaystyle\,\,\,\,\,\rightarrow\,\,\,\,\,\,0=4-(1-\alpha)^{2}\big{(}\theta{\mathfrak{q}}^{2}+(4-i\theta{\mathfrak{w}})^{2}\big{)}\,,$
(4.17a) $\displaystyle f\big{|}_{\alpha=0}=0$
$\displaystyle\,\,\,\,\,\rightarrow\,\,\,\,\,\,0=-12-\theta{\mathfrak{q}}^{2}+8i\theta{\mathfrak{w}}+\theta^{2}{\mathfrak{w}}^{2}+\alpha(2-\alpha)\big{(}\theta{\mathfrak{q}}^{2}+(4-i\theta{\mathfrak{w}})^{2}\big{)}$
(4.17b)

From (4.17a), we find

$\alpha^{*}=1-\frac{2}{\sqrt{\theta{\mathfrak{q}}^{2}+(4-i\theta{\mathfrak{w}})^{2}}}\,.$
(4.18)

According to ref. [44], the above $\alpha^{*}$ must be associated with the
threshold singularity of the response function. Substituting this value into
the second on-shell condition (4.17b), we arrive at

$-16-\theta{\mathfrak{q}}^{2}+8i\theta{\mathfrak{w}}+\theta^{2}{\mathfrak{w}}^{2}+4\sqrt{\theta{\mathfrak{q}}^{2}+(4-i\theta{\mathfrak{w}})^{2}}=\,0\,.$
(4.19)

This equation has four roots

$\displaystyle{\mathfrak{w}}=$
$\displaystyle-\frac{i}{\theta}\left(4\pm\sqrt{16-\theta{\mathfrak{q}}^{2}}\right)\,,$
(4.20a) $\displaystyle{\mathfrak{w}}=$
$\displaystyle-\frac{4i}{\theta}\pm|{\mathfrak{q}}|\sqrt{\frac{1}{\theta}}\,.$
(4.20b)

It is easy to check that the two thresholds identified by (4.20a) correspond
to ${\mathfrak{w}}_{11}$ and ${\mathfrak{w}}_{22}$ in figure 3, while the
thresholds given by (4.20b) correspond to ${\mathfrak{w}}_{12}$ and
${\mathfrak{w}}_{21}$ in the same figure. Note that threshold singularities
${\mathfrak{w}}_{ij}:i,j\in\\{1,2\\}$ shown in figure 3 were all found
numerically; now, here we see that they can be expressed by two explicit
analytical formulas (4.20a) and (4.20b).

Let us now evaluate (4.18) at (4.20a) and (4.20b). One finds

$\displaystyle\alpha_{1}\,\,\equiv\,\,\alpha^{*}\big{|}_{\text{at}\,\eqref{threshold_1}}\,\,\,\,\,=$
$\displaystyle\,\,\,\,\frac{1}{2}$ (4.21a)
$\displaystyle\alpha_{2}\,\,\equiv\,\,\alpha^{*}\big{|}_{\text{at}\,\eqref{threshold_2}}\,\,\,\,\rightarrow$
$\displaystyle-\infty$ (4.21b)

This simply tells us that ${\mathfrak{w}}_{11}$ and ${\mathfrak{w}}_{22}$ are
normal thresholds, i.e., $0<\alpha_{1}<1$, while ${\mathfrak{w}}_{12}$ and
${\mathfrak{w}}_{21}$ are pseudo-thresholds, because $\alpha_{2}\notin(0,1)$
[45].

Let us summarize. We first numerically specified the branch cut structure of
the shear stress response function, based on the formula (4.12). The results
are given in the left panel of the figure 3. Then to find the threshold
singularities shown in the figure, i.e.,
${\mathfrak{w}}_{i,j}:i,j\in\\{1,2\\}$, we analytically solved the well-known
Landau equations. The four threshold singularities found in (4.20a) and
(4.20b), and the branch cut structure of the figure 3 are exactly consistent
with the results of the loop calculations in ref. [43].

##  5 Discussion and outlook

This work focuses on correlation functions in BDNK theory, especially in
conformal systems. The time-reversal invariance apparently reduces the space
of stable and causal frames from the two-dimensional space identified by (2.2)
to a one-dimensional subspace represented by (3.8). By explicitly computing
$G_{T_{tt}T_{ti}}$ and $G_{T_{ti}T_{tt}}$, it is easy to see that the
situation is worse; only if $\theta_{\eta}=\pi_{\eta}=0$, the latter two
correlators are equal! However, as discussed in the text, these apparent
inconsistencies occur precisely outside the domain of validity of first-order
hydrodynamics.

Our central result in linear response regime is a set of Kubo formulas. One
possible place to apply these formulas is in holography, where turning on a
field perturbation in the bulk of AdS is equivalent to perturbing the boundary
theory by a specific external source. This external source is coupled to some
operator in boundary field theory. In order to use our Kubo formulas, we need
to find various correlation functions between energy density and momentum
density in boundary theory. The corresponding transport coefficients can then
be read using the associated Kubo formulas found in § 3.3. It would be very
interesting to read the BDNK transport coefficients in the far-equilibrium
holographic setup [47, 48, 49, 50, 51] and compare their late time behavior
with the near-equilibrium Kubo results.

In this work, the correlation functions, except for (3.14), were calculated in
a static fluid. It would be useful to repeat the calculations on top of a
moving fluid, i.e., $\textbf{v}_{0}\neq 0$. This might clarify the pathology
observed in the sign of correlation functions (figure 1). More generally, it
would be interesting to derive evolution equations of the correlation
functions on top of a dynamical flow [42]. This can be important when
simulating stable and causal hydrodynamics near the QCD critical point. It is
well-known that hydrodynamic fluctuations fall out of equilibrium near the
critical point [52]. Therefore, the evolution of fluctuations must be taken
into account together with the hydrodynamic evolution.

In the second part of the paper, we borrowed some methods from field theory to
find the structure of the correlation function outside the linear response
regime. By developing a numerical method for solving the on-shell conditions,
we discovered the branch-cut structure of the shear stress response function.
We also solved the Landau conditions analytically to find the threshold
singularity of the same response function. To the best of our knowledge, this
is the first time these ideas have been introduced and applied in a fluid
dynamics framework. It would be interesting to extend this analysis to the
sound channel, where one would hope to see a richer structure of correlation
functions.

Regarding the correlation functions in the nonlinear regime, it would be very
interesting to construct an effective field theory corresponding to the stable
first-order hydrodynamics. Some initial steps in this direction have recently
been taken in [43]. The theory of UV-regulated nonlinear diffusion constructed
in [43] can be regarded as the effective field theory of transverse channel
velocity fluctuations. A more general treatment requires including
fluctuations of the temperature as well as the longitudinal component of the
fluid velocity.

## Acknowledgment

We thank Ali Davody, Matthias Kaminski, and Omid Tavakol for valuable
discussions and comments. We especially thank Holmfridur S. Hannesdottir for
discussing ref. [45]. NA was supported by grant number 561119208 “Double First
Class” start-up funding of Lanzhou University, China. ST was supported by
Younger Scientist Scholarship QN2021043004 with the funding number E11I631KR0.

##  Appendix A Equilibrium fluctuations

The problem of finding the mean value of quantities in equilibrium can be
studied via finding the probability distribution of their deviations from the
associated mean values [53]. In the simplest non-relativistic setting,
quantities of interest are energy density and momentum density. What we need
then is to exploit a fundamental property of entropy: probability of finding a
subsystem at energy density $\mathcal{E}$ and momentum density $\mathcal{P}$
is given in terms of the total entropy of the system together with the medium,
$S_{t}$, and also with its change in fluctuations:

$\text{w}(\mathcal{E},\mathcal{P}_{i})\sim e^{S_{t}}\sim e^{\Delta S_{t}}\sim
e^{-\int_{\text{x}}\mathcal{R}/T}$ (A.1)

where
$\mathcal{R}=\Delta\mathcal{E}-\textbf{v}\cdot\Delta\boldsymbol{\mathcal{P}}-T\Delta\mathcal{S}$,
with v and $T$ being the equilibrium (mean) values of the velocity and the
temperature of the system. $\Delta\mathcal{E}$,
$\Delta\boldsymbol{\mathcal{P}}$, and $\Delta\mathcal{S}$ are the fluctuations
of energy density, momentum density and the entropy density. Expanding the
latter to second order in terms of the former two ones, and using the fact
that first derivatives vanish in equilibrium, we find

$\frac{\mathcal{R}}{T}\sim\frac{1}{2}\bigg{[}-\delta\left(\frac{1}{T}\right)\Delta\mathcal{E}+\delta\left(\frac{\textbf{v}}{T}\right)\cdot\Delta\boldsymbol{\mathcal{P}}\bigg{]}$
(A.2)

where $\delta$ and $\Delta$ indicate the fluctuation in source and the
fluctuation in energy/momentum density, respectively. When
$\textbf{v}=\textbf{0}$ in equilibrium, $\Delta\mathcal{E}=c_{V}\delta T$ and
$\Delta\boldsymbol{\mathcal{P}}=w\delta\textbf{v}$. Therefore (A.2) becomes

$\frac{\mathcal{R}}{T}\sim\frac{1}{2}\bigg{[}\frac{c_{V}}{T^{2}}(\delta
T)^{2}+\frac{w}{T}(\delta\textbf{v})^{2}\bigg{]}$ (A.3)

Substituting this into (A.1), the correlation functions in complete
equilibrium can be read immediately as follows [54]:

$\displaystyle\langle\delta v_{i}(\textbf{x}_{1})\delta
v_{j}(\textbf{x}_{2})\rangle=$
$\displaystyle\,\delta_{ij}\frac{T}{w}\,\delta^{3}(\textbf{x}_{1}-\textbf{x}_{2})\,,$
(A.4a) $\displaystyle\langle\delta T(\textbf{x}_{1})\delta
T(\textbf{x}_{2})\rangle=$
$\displaystyle\,\frac{T^{2}}{c_{V}}\,\delta^{3}(\textbf{x}_{1}-\textbf{x}_{2})\,.$
(A.4b)

##  Appendix B Hydrodynamic fluctuations

It is easy to show that

$\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}_{\omega\textbf{k}}=\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}^{(+)}_{\omega\textbf{k}}+\big{\langle}\delta\phi_{b}\delta\phi_{a}\big{\rangle}^{(+)}_{-\omega\,-\textbf{k}}$
(B.1)

where
$\big{\langle}\delta\phi_{a}\delta\phi_{b}\big{\rangle}^{(+)}_{\omega\textbf{k}}$
is one-sided Fourier transformation of
$\langle\delta\phi_{a}(t,\textbf{x})\delta\phi_{b}(0,\textbf{0})\rangle$,
defined as (3.4), but with the time integral taken from $0$ to $+\infty$ [53,
54]. The strategy is to find equations between one-sided Fourier transforms
and then to calculate correlation functions by use of (B.1).

For convenience, let us first define dimensionless quantities as

${\mathfrak{w}}=\frac{\bar{\eta}}{\bar{p}}\,\frac{\omega}{T}\,,\,\,\,\,\,\,{\mathfrak{q}}=\frac{\bar{\eta}}{\bar{p}}\,\frac{\textbf{k}}{T}\,,\,\,\,\,\,\,\pi_{\eta}=\frac{\bar{\pi}}{\bar{\eta}}\,,\,\,\,\,\,\,\theta_{\eta}=\frac{\bar{\theta}}{\bar{\eta}}$
(B.2)

We also take the momentum to be directed in the $z-$direction
${\mathfrak{q}}=(0,0,q)$. This then allows us to calculate the correlation
functions in two channels. In the longitudinal (sound) channel, the
correlation functions of $\delta T/t$ and $\delta v^{\parallel}$ are involved.
In the transverse (shear) channel, only $\langle\delta
v^{\perp}v^{\perp}\rangle$ needs to be computed.

### Longitudinal channel

Multiplying the two equations (3.1) and the $z-$component of (3.2) by $\delta
T(0,\textbf{0})$ from the right side, and averaging over the local
equilibrium, we find the equations of one-sided Fourier transformation as the
following:

$\begin{split}\left[-i{\mathfrak{w}}-\frac{1}{12}(\theta_{\eta}{\mathfrak{q}}^{2}+3\pi_{\eta}{\mathfrak{w}}^{2})\right]\langle\frac{\delta
T}{T}\frac{\delta
T}{T}\rangle^{(+)}+\bigg{[}\frac{iq}{3}+\frac{1}{12}(\theta_{\eta}+\pi_{\eta})q{\mathfrak{w}}\bigg{]}&\langle\delta
v_{z}\frac{\delta
T}{T}\rangle^{(+)}=\frac{\bar{\eta}(4-i\pi_{\eta}{\mathfrak{w}})}{48\bar{p}^{2}\,T^{4}}\\\
\bigg{[}-i{\mathfrak{w}}+\frac{1}{3}{\mathfrak{q}}^{2}-\frac{1}{12}(\pi_{\eta}{\mathfrak{q}}^{2}+3\theta_{\eta}{\mathfrak{w}}^{2})\bigg{]}\langle\delta
v_{z}\frac{\delta
T}{T}\rangle^{(+)}+\bigg{[}iq+\frac{1}{4}(\theta_{\eta}+\pi_{\eta})q{\mathfrak{w}}\bigg{]}&\langle\frac{\delta
T}{T}\frac{\delta T}{T}\rangle^{(+)}=\,0\end{split}$ (B.3)

Note that the right side of the first equation above is the consequence of
using (A.4a). When multiplying (3.1) and the $z-$component of (3.2) by $\delta
v_{z}(0,\textbf{0})$, and repeating the above process, we arrive at

$\begin{split}\left[-i{\mathfrak{w}}-\frac{1}{12}(\theta_{\eta}{\mathfrak{q}}^{2}+3\pi_{\eta}{\mathfrak{w}}^{2})\right]\langle\frac{\delta
T}{T}\delta
v_{z}\rangle^{(+)}+\left[\frac{iq}{3}+\frac{1}{12}(\theta_{\eta}+\pi_{\eta})q{\mathfrak{w}}\right]&\langle\delta
v_{z}\delta v_{z}\rangle^{(+)}=0\\\
\bigg{[}-i{\mathfrak{w}}+\frac{1}{3}{\mathfrak{q}}^{2}-\frac{1}{12}(\pi_{\eta}{\mathfrak{q}}^{2}+3\theta_{\eta}{\mathfrak{w}}^{2})\bigg{]}\langle\delta
v_{z}\delta
v_{z}\rangle^{(+)}+\bigg{[}iq+\frac{1}{4}(\theta_{\eta}+\pi_{\eta})q{\mathfrak{w}}\bigg{]}&\langle\frac{\delta
T}{T}\delta
v_{z}\rangle^{(+)}=\,\frac{\bar{\eta}(4-i\theta_{\eta}{\mathfrak{w}})}{16\bar{p}^{2}\,T^{4}}\end{split}$
(B.4)

From these equations we find

$\displaystyle\langle\frac{\delta T}{T}\frac{\delta T}{T}\rangle^{(+)}=$
$\displaystyle\,\frac{\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}_{L}({\mathfrak{w}},{\mathfrak{q}})}\,(i\pi_{\eta}{\mathfrak{w}}-4)\,\big{(}(\pi_{\eta}-4){\mathfrak{q}}^{2}+12i{\mathfrak{w}}+3\theta_{\eta}{\mathfrak{w}}^{2}\big{)}$
(B.5a) $\displaystyle\langle\delta v_{z}\frac{\delta T}{T}\rangle^{(+)}=$
$\displaystyle\,\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}_{L}({\mathfrak{w}},{\mathfrak{q}})}\,iq\,(i\pi_{\eta}{\mathfrak{w}}-4)\,\big{(}4-i(\theta_{\eta}+\pi_{\eta}){\mathfrak{w}}\big{)}$
(B.5b) $\displaystyle\langle\frac{\delta T}{T}\delta v_{z}\rangle^{(+)}=$
$\displaystyle\,\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}_{L}({\mathfrak{w}},{\mathfrak{q}})}\,iq\,(i\theta_{\eta}{\mathfrak{w}}-4)\,\big{(}4-i(\theta_{\eta}+\pi_{\eta}){\mathfrak{w}}\big{)}$
(B.5c) $\displaystyle\langle\delta v_{z}\delta v_{z}\rangle^{(+)}=$
$\displaystyle\,\frac{3\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{1}{\mathcal{D}_{L}({\mathfrak{w}},{\mathfrak{q}})}\,(i\theta_{\eta}{\mathfrak{w}}-4)\,\big{(}\theta_{\eta}{\mathfrak{q}}^{2}+12i{\mathfrak{w}}+3\pi_{\eta}{\mathfrak{w}}^{2}\big{)}$
(B.5d)

with

$\begin{split}\mathcal{D}_{L}({\mathfrak{w}},{\mathfrak{q}})=&\,9\pi_{\eta}\,\theta_{\eta}\,{\mathfrak{w}}^{4}+36i\left(\theta_{\eta}+\pi_{\eta}\right)\,{\mathfrak{w}}^{3}-6\left(\pi_{\eta}\left(\theta_{\eta}+2\right){\mathfrak{q}}^{2}+24\right){\mathfrak{w}}^{2}\\\
&-12i\left(\theta_{\eta}+\pi_{\eta}+4\right){\mathfrak{q}}^{2}{\mathfrak{w}}+\left(\left(\pi_{\eta}-4\right)\theta_{\eta}{\mathfrak{q}}^{2}+48\right){\mathfrak{q}}^{2}\end{split}$
(B.6)

Then it is easy to calculate the correlation functions by using (B.1). We do
not represent those expressions here explicitly.

### Transverse channel

The correlation functions in the transverse channel are not complicated, even
at $\textbf{v}_{0}\neq 0$. Therefore, we continue to express the result in
general $\textbf{v}_{0}$.

Due to the rotational symmetry in the transverse plane, we only need to
consider the $x-$component of (3.2). Multiplying by $\delta
v_{x}(0,\textbf{0})$, the one-sided Fourier transformation of $\langle\delta
v_{x}(t,\textbf{x})\delta v_{x}(0,\textbf{0})\rangle$ is found to be

$\langle\delta v_{i}\delta
v_{i}\rangle^{(+)}=\,\frac{i\bar{\eta}}{4\bar{p}^{2}T^{4}}\frac{4i\gamma+{\mathfrak{w}}+\gamma^{2}(\theta_{\eta}-1){\mathfrak{w}}}{\mathcal{D}_{T}({\mathfrak{w}},{\mathfrak{q}})}\,\,\,\,\,\,i=x,y$
(B.7)

where

$\begin{split}\mathcal{D}_{T}({\mathfrak{w}},{\mathfrak{q}})=&\big{(}1+\gamma^{2}(\theta_{\eta}-1)\big{)}{\mathfrak{w}}^{2}+2\gamma\big{(}2i-\gamma(\theta_{\eta}-1){\mathfrak{q}}\cdot\textbf{v}_{0}\big{)}{\mathfrak{w}}\\\
&-4i\gamma{\mathfrak{q}}\cdot\textbf{v}_{0}+\big{(}\gamma^{2}(\theta_{\eta}-1)-\theta_{\eta}\big{)}{\mathfrak{q}}^{2}\end{split}$
(B.8)

Again, it is easy to calculate $\langle\delta v_{i}\delta v_{i}\rangle$ by
using (B.1).

##  Appendix C Landau equations

Let us suppose an arbitrary Feynman diagram (in the scalar theory) is
represented by the following integral

$\int\frac{d^{4}k}{(2\pi)^{4}}\frac{d^{4}l}{(2\pi)^{4}}\cdots\frac{B}{A_{1}\,A_{2}\,\cdots}$
(C.1)

where $A_{i}=q_{i}^{2}+m^{2}$ and $q_{i}$ is the four-momentum corresponding
to the given line in the diagram. In addition, $B$ is a polynomial of the
four-vectors $q_{i}$. Using the well known “Feynman’ parameters” method, one
may write

$\frac{1}{A_{1}A_{2}\cdots
A_{n}}=(n-1)!\int_{0}^{1}\int_{0}^{1}\cdots\int_{0}^{1}\frac{d\alpha_{1}d\alpha_{2}\cdots
d\alpha_{n}{\delta(\alpha_{1}+\alpha_{2}+\cdots+\alpha_{n}-1)}}{\big{(}\alpha_{1}A_{1}+\alpha_{2}A_{2}+\cdots\alpha_{n}A_{n}\big{)}^{n}}$
(C.2)

Landau argues that the singularities of the integral (C.1)111111These
singularities include both normal (discussed in this paper) and anomalous
thresholds. For anomalous thresholds see [55, 56]. are actually a solution to
the following equations/conditions [44] :

  1. 1.

For each propagator $i=1,\cdots,n$:

$q_{i}^{2}+m^{2}=0\,\,\,\,\,\,\,\text{or}\,\,\,\,\,\,\,\,\alpha_{i}=0\,.$
(C.3)

  2. 2.

For each loop momentum integration variable $k_{j}$:

$\frac{\partial}{\partial k_{j}}\sum_{i}\alpha_{i}(q_{i}^{2}+m^{2})=0.$ (C.4)

These conditions are referred to as the “Landau equations” or “Landau
conditions”. The first condition, (C.3), is simply the on-shell condition for
the internal momenta. the second one, given by (C.4), is called the “Landau
loop equation”.

To understand (C.4), let us denote that the singularity of the integrand will
be the singularity of the integral if it cannot be avoided by doing any
contour deformation. One way to achieve this is that the singularity of the
integrand reaches the endpoint of the integration; the latter is called “end-
point singularity” which is actually quite easy to identify. Another
possibility is that at some point, the contour of integration is pinched
between singularities of the integrand. Remember when two curves are tangent;
their normal vectors at the point of intersection are parallel. Similarly,
there will be a contour pinch when the gradients of the on-shell conditions
(with respect to the independent loop momenta) are linearly dependent. The
latter singularity is called “pinch singularity” and determined by equation
(C.4) [45, 46] 121212See also [57, 58]..

## References

  * [1] P. Kovtun, “Lectures on hydrodynamic fluctuations in relativistic theories,” J. Phys. A 45 (2012), 473001 [arXiv:1205.5040 [hep-th]].
  * [2] C. Eckart, “The Thermodynamics of irreversible processes. 3.. Relativistic theory of the simple fluid,” Phys. Rev. 58 (1940), 919-924
  * [3] L. D. Landau and E. M. Lifshitz, “Fluid Mechanics,” Course of Theoretical Physics, Vol 6 (Elsevier Science, 2013).
  * [4] W. A. Hiscock and L. Lindblom, “Generic instabilities in first-order dissipative relativistic fluid theories,” Phys. Rev. D 31 (1985), 725-733.
  * [5] W. A. Hiscock and L. Lindblom, “Linear plane waves in dissipative relativistic fluids,” Phys. Rev. D 35 (1987), 3723-3732.
  * [6] W. Israel, “Nonstationary irreversible thermodynamics: A Causal relativistic theory,” Annals Phys. 100 (1976) 310.
  * [7] W. Israel and J. M. Stewart, “Thermodynamics of nonstationary and transient effects in a relativistic gas,” Phys. Lett. A 58 (1976) 213.
  * [8] P. Van and T. S. Biro, “First order and stable relativistic dissipative hydrodynamics,” Phys. Lett. B 709 (2012), 106-110 [arXiv:1109.0985 [nucl-th]].
  * [9] H. Freistuhler and B. Temple, “Causal dissipation for the relativistic dynamics of ideal gases,” Proceedings of the Royal Society A 473 (2017) 0729.
  * [10] F. S. Bemfica, M. M. Disconzi and J. Noronha, “Causality and existence of solutions of relativistic viscous fluid dynamics with gravity,” Phys. Rev. D 98 (2018) no.10, 104064 [arXiv:1708.06255 [gr-qc]].
  * [11] P. Kovtun, “First-order relativistic hydrodynamics is stable,” JHEP 10 (2019), 034 [arXiv:1907.08191 [hep-th]].
  * [12] F. S. Bemfica, F. S. Bemfica, M. M. Disconzi, M. M. Disconzi, J. Noronha and J. Noronha, “Nonlinear Causality of General First-Order Relativistic Viscous Hydrodynamics,” Phys. Rev. D 100 (2019) no.10, 104020 [erratum: Phys. Rev. D 105 (2022) no.6, 069902] [arXiv:1907.12695 [gr-qc]].
  * [13] A. Das, W. Florkowski, J. Noronha and R. Ryblewski, “Equivalence between first-order causal and stable hydrodynamics and Israel-Stewart theory for boost-invariant systems with a constant relaxation time,” Phys. Lett. B 806 (2020), 135525 [arXiv:2001.07983 [nucl-th]].
  * [14] M. Shokri and F. Taghinavaz, “Conformal Bjorken flow in the general frame and its attractor: Similarities and discrepancies with the Müller-Israel-Stewart formalism,” Phys. Rev. D 102 (2020) no.3, 036022 [arXiv:2002.04719 [hep-th]].
  * [15] F. Taghinavaz, “Causality and Stability Conditions of a Conformal Charged Fluid,” JHEP 08 (2020), 119 [arXiv:2004.01897 [hep-th]].
  * [16] R. E. Hoult and P. Kovtun, “Stable and causal relativistic Navier-Stokes equations,” JHEP 06 (2020), 067 [arXiv:2004.04102 [hep-th]].
  * [17] A. Das, W. Florkowski and R. Ryblewski, “Correspondence between Israel-Stewart and first-order casual and stable hydrodynamics for the boost-invariant massive case with zero baryon density,” Phys. Rev. D 102 (2020) no.3, 031501 [arXiv:2006.00536 [nucl-th]].
  * [18] L. Gavassino, M. Antonelli and B. Haskell, “When the entropy has no maximum: a new perspective on the instability of the first-order theories of dissipation,” Phys. Rev. D 102 (2020) no.4, 043018 [arXiv:2006.09843 [gr-qc]].
  * [19] F. S. Bemfica, M. M. Disconzi and J. Noronha, “First-Order General-Relativistic Viscous Fluid Dynamics,” Phys. Rev. X 12 (2022) no.2, 021044 [arXiv:2009.11388 [gr-qc]].
  * [20] C. Plumberg, D. Almaalol, T. Dore, J. Noronha and J. Noronha-Hostler, “Causality violations in realistic simulations of heavy-ion collisions,” Phys. Rev. C 105 (2022) no.6, L061901 [arXiv:2103.15889 [nucl-th]].
  * [21] H. Freistuhler, “Nonexistence and existence of shock profiles in the Bemfica-Disconzi-Noronha model,” Phys. Rev. D 103 (2021) no.12, 124045 [arXiv:2103.16661 [math.AP]].
  * [22] A. Pandya and F. Pretorius, “Numerical exploration of first-order relativistic hydrodynamics,” Phys. Rev. D 104 (2021) no.2, 023015 [arXiv:2104.00804 [gr-qc]].
  * [23] J. Noronha, M. Spaliński and E. Speranza, “Transient Relativistic Fluid Dynamics in a General Hydrodynamic Frame,” Phys. Rev. Lett. 128 (2022) no.25, 252302 [arXiv:2105.01034 [nucl-th]].
  * [24] S. Mitra, “Causality and stability analysis of first-order field redefinition in relativistic hydrodynamics from kinetic theory,” Phys. Rev. C 105 (2022) no.5, 054910 [arXiv:2106.08510 [nucl-th]].
  * [25] L. Gavassino, “Can We Make Sense of Dissipation without Causality?,” Phys. Rev. X 12 (2022) no.4, 041001 [arXiv:2111.05254 [gr-qc]].
  * [26] R. E. Hoult and P. Kovtun, “Causal first-order hydrodynamics from kinetic theory and holography,” Phys. Rev. D 106 (2022) no.6, 066023 [arXiv:2112.14042 [hep-th]].
  * [27] J. Armas and F. Camilloni, “A stable and causal model of magnetohydrodynamics,” JCAP 10 (2022), 039 [arXiv:2201.06847 [hep-th]].
  * [28] H. Bantilan, Y. Bea and P. Figueras, “Evolutions in first-order viscous hydrodynamics,” JHEP 08 (2022), 298 [arXiv:2201.13359 [hep-th]].
  * [29] R. Biswas, S. Mitra and V. Roy, “Is first-order relativistic hydrodynamics in a general frame stable and causal for arbitrary interactions?,” Phys. Rev. D 106 (2022) no.1, L011501 [arXiv:2202.08685 [nucl-th]].
  * [30] A. Daher, A. Das and R. Ryblewski, “Stability studies of first order spin-hydrodynamic frameworks,” [arXiv:2209.10460 [nucl-th]].
  * [31] F. S. Bemfica, M. M. Disconzi, J. Noronha and R. J. Scherrer, “Cosmological consequences of first-order general-relativistic viscous fluid dynamics,” [arXiv:2210.13372 [gr-qc]].
  * [32] M. Luzum and P. Romatschke, “Conformal Relativistic Viscous Hydrodynamics: Applications to RHIC results at s(NN)**(1/2) = 200-GeV,” Phys. Rev. C 78 (2008), 034915 [erratum: Phys. Rev. C 79 (2009), 039903] [arXiv:0804.4015 [nucl-th]].
  * [33] W. Florkowski, M. P. Heller and M. Spalinski, “New theories of relativistic hydrodynamics in the LHC era,” Rept. Prog. Phys. 81 (2018) no.4, 046001 [arXiv:1707.02282 [hep-ph]].
  * [34] C. Shen and L. Yan, “Recent development of hydrodynamic modeling in heavy-ion collisions,” Nucl. Sci. Tech. 31 no.12, 122 [arXiv:2010.12377 [nucl-th]].
  * [35] M. Crossley, P. Glorioso and H. Liu, “Effective field theory of dissipative fluids,” JHEP 09 (2017), 095 [arXiv:1511.03646 [hep-th]].
  * [36] K. Jensen, N. Pinzani-Fokeeva and A. Yarom, “Dissipative hydrodynamics in superspace,” JHEP 09 (2018), 127 [arXiv:1701.07436 [hep-th]].
  * [37] F. M. Haehl, R. Loganayagam and M. Rangamani, “Topological sigma models & dissipative hydrodynamics,” JHEP 04 (2016), 039 [arXiv:1511.07809 [hep-th]].
  * [38] H. Liu and P. Glorioso, “Lectures on non-equilibrium effective field theories and fluctuating hydrodynamics,” PoS TASI2017 (2018), 008 [arXiv:1805.09331 [hep-th]].
  * [39] X. Chen-Lin, L. V. Delacrétaz and S. A. Hartnoll, “Theory of diffusive fluctuations,” Phys. Rev. Lett. 122 (2019) no.9, 091602 [arXiv:1811.12540 [hep-th]].
  * [40] P. Kovtun, G. D. Moore and P. Romatschke, “The stickiness of sound: An absolute lower limit on viscosity and the breakdown of second order relativistic hydrodynamics,” Phys. Rev. D 84 (2011), 025006 [arXiv:1104.1586 [hep-ph]].
  * [41] P. Kovtun and L. G. Yaffe, “Hydrodynamic fluctuations, long time tails, and supersymmetry,” Phys. Rev. D 68 (2003), 025007 [arXiv:hep-th/0303010 [hep-th]].
  * [42] X. An, G. Basar, M. Stephanov and H. U. Yee, “Relativistic Hydrodynamic Fluctuations,” Phys. Rev. C 100 (2019) no.2, 024910 [arXiv:1902.09517 [hep-th]].
  * [43] N. Abbasi, M. Kaminski and O. Tavakol, “Ultraviolet-regulated theory of non-linear diffusion,” [arXiv:2212.11499 [hep-th]].
  * [44] L. D. Landau, “On analytic properties of vertex parts in quantum field theory,” Nuclear Physics 13, (1959), 181-192.
  * [45] H. S. Hannesdottir, A. J. McLeod, M. D. Schwartz and C. Vergu, “Constraints on Sequential Discontinuities from the Geometry of On-shell Spaces,” [arXiv:2211.07633 [hep-th]].
  * [46] R. J. Eden, P. V. Landshoff, D. I. Olive and J. C. Polkinghorne, “The Analytic S-Matrix.” Cambridge University Press, 1966.
  * [47] Y. Bu and M. Lublinsky, “All order linearized hydrodynamics from fluid-gravity correspondence,” Phys. Rev. D 90 (2014) no.8, 086003 [arXiv:1406.7222 [hep-th]].
  * [48] Y. Bu and M. Lublinsky, “Linearized fluid/gravity correspondence: from shear viscosity to all order hydrodynamics,” JHEP 11 (2014), 064 [arXiv:1409.3095 [hep-th]].
  * [49] P. Romatschke, “Relativistic Fluid Dynamics Far From Local Equilibrium,” Phys. Rev. Lett. 120 (2018) no.1, 012301 [arXiv:1704.08699 [hep-th]].
  * [50] M. Baggioli, L. Li and H. T. Sun, “Shear Flows in Far-from-Equilibrium Strongly Coupled Fluids,” Phys. Rev. Lett. 129 (2022) no.1, 011602 [arXiv:2112.14855 [hep-th]].
  * [51] C. Cartwright, M. Kaminski and M. Knipfer, “Hydrodynamic attractors for the speed of sound in holographic Bjorken flow,” [arXiv:2207.02875 [hep-th]].
  * [52] X. An, M. Bluhm, L. Du, G. V. Dunne, H. Elfner, C. Gale, J. Grefa, U. Heinz, A. Huang and J. M. Karthein, 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]].
  * [53] L. D. Landau and E. M. Lifshitz, “Statistical Physics Part 1,” Course of Theoretical Physics, Vol 5 (Elsevier Science, 2013).
  * [54] L. D. Landau and E. M. Lifshitz, “Statistical Physics Part 2,” Course of Theoretical Physics, Vol 9 (Elsevier Science, 2013).
  * [55] Y. Nambu, “Dispersion relations for form factors.,” Nuovo Cim 9, 610–623 (1958).
  * [56] R. Karplus, C. M. Sommerfield and E. H. Wichmann, “Spectral Representations in Perturbation Theory. 1. Vertex Function,” Phys. Rev. 111 (1958), 1187-1190.
  * [57] A. Windisch, M. Q. Huber and R. Alkofer, “How to determine the branch points of correlation functions in Euclidean space,” Acta Phys. Polon. Supp. 6 (2013) no.3, 887-892 [arXiv:1304.3642 [hep-ph]].
  * [58] M. Q. Huber, W. Kern and R. Alkofer, “On the analytic structure of three-point functions from contour deformations,” [arXiv:2212.02515 [hep-ph]].