*Article* **Relativistic Dissipative Fluid Dynamics from the Non-Equilibrium Statistical Operator**

**Arus Harutyunyan 1,\* , Armen Sedrakian <sup>2</sup> and Dirk H. Rischke 1,3**


Received: 26 April 2018; Accepted: 19 June 2018; Published: 21 June 2018

**Abstract:** We present a new derivation of second-order relativistic dissipative fluid dynamics for quantum systems using Zubarev's formalism for the non-equilibrium statistical operator. In particular, we discuss the shear-stress tensor to second order in gradients and argue that the relaxation terms for the dissipative quantities arise from memory effects contained in the statistical operator. We also identify new transport coefficients which describe the relaxation of dissipative processes to second order and express them in terms of equilibrium correlation functions, thus establishing Kubo-type formulae for the second-order transport coefficients.

**Keywords:** relativistic fluid dynamics; statistical operator; non-equilibrium states; transport coefficients; correlation functions

#### **1. Introduction**

Fluid dynamics is a powerful tool to describe low-frequency and long-wavelength phenomena in statistical systems [1]. It finds numerous applications in astrophysics, cosmology, heavy-ion physics, and other areas. In particular, it has been successfully applied to describe the collective behavior of hot and dense strongly interacting matter created in heavy-ion collision experiments at Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC). In these experiments, a new state of matter, the quark-gluon plasma (QGP), was discovered, which behaves almost like a perfect fluid.

There are two main approaches which can be used to derive the equations of motion of fluid dynamics and the pertaining transport coefficients from the underlying microscopic theory. For weakly interacting systems, one commonly relies on kinetic theory based on the Boltzmann equation for the quasi-particle distribution function [2–5]. For strongly interacting quantum systems, where the quasi-particle picture breaks down and/or the quantum nature of the fields itself is important, kinetic theory is no longer applicable, and a full quantum-statistical approach based on the Liouville equation for the non-equilibrium statistical operator is required.

In this work, we adopt the method of the non-equilibrium statistical operator (NESO) [6,7] to obtain the relativistic fluid-dynamical equations of motion for strongly correlated matter, such as the QGP, in the non-perturbative regime. The method was applied to quantum fields [8] and has been since extended to treat systems in strong magnetic fields [9]. It is based on a generalization of the Gibbs canonical ensemble to non-equilibrium states, i.e., the statistical operator is promoted to a non-local functional of the thermodynamic parameters and their space-time derivatives. Assuming that the thermodynamic parameters are sufficiently smooth over the correlation lengths characterizing the system, the statistical operator is expanded into a series in gradients of these parameters to the desired order. The fluid-dynamical equations for the dissipative fluxes emerge then after statistically averaging the relevant quantum operators. An advantage of the NESO method is that the transport coefficients of the system are automatically obtained in the form of Kubo-type relations, i.e., they are related to certain correlation functions of the underlying field theory in the strong-coupling regime. There exist a number of formulations of relativistic fluid dynamics in terms of near-equilibrium quantities which are related to the NESO method employed by us; for recent work, see References [10–12].

This contribution provides a concise presentation of our recent work on the derivation of second-order dissipative fluid dynamics via the NESO method [13,14]. As is well known, relativistic fluid dynamics describes the state of a fluid in terms of its energy-momentum tensor and currents of conserved charges, which in the relevant low-frequency and long-wavelength limit can be expanded around their equilibrium values. The zeroth-order expansion corresponds to ideal (non-dissipative) fluid dynamics. At first order, dissipative relativistic fluid dynamics emerges from a truncation that keeps the terms of linear order in gradients [1,15]. Second-order relativistic theories have also been constructed [16,17] to avoid the acausality of the first-order theory and the resulting numerical instabilities. In second-order theories, the dissipative fluxes satisfy relaxation equations, which describe the process of their relaxation towards their Navier–Stokes values at asymptotically large times. While the general structure of second-order fluid dynamics is known, different results have been obtained for the coefficients entering these equations (see, e.g., References [18,19]). The various versions of second-order fluid dynamics and the pertaining relaxation equations are reviewed and compared to each other, e.g., in the review articles [2,4,5], to which we refer the reader for more detailed expositions.

This work is structured as follows. Section 2 gives a brief summary of Zubarev's formalism for the NESO [6,7]. Section 3 recapitulates Navier–Stokes theory and the Kubo formulae for the first-order transport coefficients. The second-order transport equations are discussed in Section 4 and a summary is given in Section 5. We work in flat space-time described by the metric tensor *<sup>g</sup>μν* <sup>=</sup> diag(+, <sup>−</sup>, <sup>−</sup>, <sup>−</sup>).

#### **2. Non-Equilibrium Statistical Operator and Correlation Functions**

The fluid-dynamical state of a relativistic quantum system is described by the operators of the energy-momentum tensor *T*ˆ *μν*(*x*) and the conserved particle current *N*ˆ *<sup>μ</sup>*(*x*). For example, in the case of Dirac fermions, the particle current is given by *N*ˆ *<sup>μ</sup>* = *ψγ*ˆ¯ *μψ*ˆ, where *ψ*ˆ is the Dirac field operator, and *γ<sup>μ</sup>* are the Dirac matrices. The equations of relativistic fluid dynamics consist of the covariant conservation laws for these quantities

$$
\partial\_{\mu}\hat{T}^{\mu\nu}(\mathbf{x}) = 0, \qquad \partial\_{\mu}\hat{N}^{\mu}(\mathbf{x}) = 0. \tag{1}
$$

Here, we assume that the fluid consists of only one particle species. The generalization to the case of several conserved species is straightforward and is given elsewhere [13].

In general, the fluid-dynamical description is applicable, if the actual state of a given system does not deviate too much from local thermodynamic equilibrium. This allows one to introduce a fictitious local-equilibrium reference state, characterized by space-time dependent thermodynamic parameters, such as temperature *<sup>T</sup>*(*x*) <sup>≡</sup> *<sup>β</sup>*−1(*x*), chemical potential *<sup>μ</sup>*(*x*), and fluid 4-velocity *<sup>u</sup>ν*(*x*). The deviation of the actual state from this fictitious reference state is then taken to be proportional to gradients of these fields. The assumption that the deviation from local equilibrium is small is equivalent to assuming that these fields are slowly varying functions of the space-time coordinates *x* ≡ (*x*, *t*). Note that, in this context, "slowly" means that the characteristic *macroscopic* scales over which the fluid-dynamical quantities change in space and time should be much larger than the characteristic *microscopic* scales of the system, e.g., for quasi-particles the mean free path between collisions. In terms of the thermodynamic parameters defined above, we define new auxiliary functions

$$
\beta^v(\mathbf{x}) = \beta(\mathbf{x}) u^v(\mathbf{x}), \qquad u(\mathbf{x}) = \beta(\mathbf{x}) \mu(\mathbf{x}).\tag{2}
$$

Now, consider the NESO given by Huang et al. [9]:

$$\rho(t) = Q^{-1}e^{-\vec{A} + \vec{B}}, \qquad Q = \text{Tr}e^{-\vec{A} + \vec{B}}, \tag{3}$$

where

$$\hat{A}(t) \quad = \int d^3 \mathbf{x} \left[ \mathcal{J}^\nu(\mathbf{x}) \hat{\mathcal{I}}\_{0\nu}(\mathbf{x}) - a(\mathbf{x}) \hat{\mathcal{N}}^0(\mathbf{x}) \right],\tag{4}$$

$$\hat{B}(t) = \int d^3x\_1 \int\_{-\infty}^t dt\_1 e^{\varepsilon(t\_1 - t)} \hat{\mathcal{C}}(x\_1),\tag{5}$$

$$
\hat{\mathcal{L}}(\mathbf{x}) \quad = \quad \hat{T}\_{\mu\nu}(\mathbf{x})\partial^{\mu}\mathcal{J}^{\nu}(\mathbf{x}) - \hat{N}^{\mu}(\mathbf{x})\partial\_{\mu}a(\mathbf{x}), \tag{6}
$$

with *ε* → +0 taken after the thermodynamic limit. The NESO satisfies the quantum Liouville equation with an infinitesimal source term *ε*, which for positive values selects the retarded solution [6,7]. The operators *A*ˆ(*t*) and *B*ˆ(*t*) correspond to the equilibrium and non-equilibrium parts of the statistical operator, where the operator *C*ˆ(*x*) stands for the thermodynamic "force" as it involves the gradients of the thermodynamic variables, i.e., temperature, chemical potential, and fluid 4-velocity. We also define the *local-equilibrium statistical operator* as

$$\phi\_l(t) = Q\_l^{-1} e^{-\hat{A}}, \qquad Q\_l = \text{Tr}e^{-\hat{A}},\tag{7}$$

which is the analog of the Gibbs distribution involving local thermodynamic parameters.

Before proceeding, we remark that the thermodynamic variables are well-defined quantities only in an equilibrium state, but not for a non-equilibrium state. The reason they appear at all in our discussion is the introduction of a fictitious local-equilibrium state, from which the actual non-equilibrium state should not deviate too much. The freedom in choosing this fictitious state can be exploited to determine the parameters *α*(*x*), *β*(*x*), and the fluid 4-velocity *uν*(*x*) characterizing this state. For this purpose, we first define the operators of the energy and particle densities via ˆ(*x*) = *uμ*(*x*)*uν*(*x*)*T*ˆ *μν*(*x*) and *n*ˆ(*x*) = *uμ*(*x*)*N*ˆ *<sup>μ</sup>*(*x*). These simply imply that ˆ(*x*) and *n*ˆ(*x*) are the time-like eigenvalues of the energy-momentum tensor and the particle current, respectively, measured by a local observer comoving with a fluid element. The local values of the Lorentz-invariant thermodynamic parameters *β*(*x*) and *α*(*x*) can now be fixed by requiring that the average values of the operators ˆ(*x*) and *n*ˆ(*x*) match the local-equilibrium values of these quantities. These so-called *Landau matching conditions* [6,7] are then written as

$$
\langle \mathfrak{k}(\mathbf{x}) \rangle = \langle \mathfrak{k}(\mathbf{x}) \rangle\_{l\prime} \qquad \langle \mathfrak{n}(\mathbf{x}) \rangle = \langle \mathfrak{n}(\mathbf{x}) \rangle\_{l\prime} \tag{8}
$$

where for an arbitrary operator *X*ˆ(*x*) the non-equilibrium and local-equilibrium statistical averages are defined as

$$
\langle \hat{X}(\mathbf{x}) \rangle = \text{Tr}[\not p(t)\hat{X}(\mathbf{x})], \qquad \langle \hat{X}(\mathbf{x}) \rangle\_l = \text{Tr}\left[\not p\_l(t)\hat{X}(\mathbf{x})\right]. \tag{9}
$$

Finally, the fluid 4-velocity *u<sup>μ</sup>* can be determined by relating it to a particular physical current. For example, in the Landau–Lifshitz frame, the 4-velocity is parallel to the fluid 4-momentum or, equivalently, to the energy flow, i.e., *<sup>u</sup>μT*<sup>ˆ</sup> *μν* <sup>=</sup> <sup>ˆ</sup>*u<sup>ν</sup>* [1]. In the Eckart frame, the fluid velocity is associated with the particle flow via *N*<sup>ˆ</sup> *<sup>μ</sup>* <sup>=</sup> *n*ˆ*u<sup>μ</sup>* [15]. However, in the following, we keep the fluid velocity generic without specifying any particular reference frame.

The next step is to expand the NESO around the local-equilibrium value in Equation (7) treating the non-equilibrium part, which is described by the operator *B*ˆ, as a perturbation

$$
\not\!\!\!\/ = \not\!\!\/ \_l + \not\!\!\/ \_l + \not\!\!\/ \_2 \tag{10}
$$

where the first-order term is given by

$$\left[\rho\_1(t) = \int d^4x\_1 \int\_0^1 d\tau \left[\mathcal{C}\_\tau(x\_1) - \langle \mathcal{C}\_\tau(x\_1) \rangle\_l \right] \rho\_{l\prime} \tag{11}$$

while the second-order term is

$$
\begin{split}
\partial\_{2}(t) &= \frac{1}{2} \int d^{4} \mathbf{x}\_{1} d^{4} \mathbf{x}\_{2} \int\_{0}^{1} d\tau \int\_{0}^{1} d\lambda \left[ \mathcal{T} \{ \mathbf{\hat{C}}\_{\lambda}(\mathbf{x}\_{1}) \mathbf{\hat{C}}\_{\tau}(\mathbf{x}\_{2}) \} - \langle \mathcal{T} \{ \mathbf{\hat{C}}\_{\lambda}(\mathbf{x}\_{1}) \mathbf{\hat{C}}\_{\tau}(\mathbf{x}\_{2}) \} \rangle\_{l} \right. \\ & \left. - \langle \mathbf{\hat{C}}\_{\lambda}(\mathbf{x}\_{1}) \rangle\_{l} \mathbf{\hat{C}}\_{\tau}(\mathbf{x}\_{2}) - \mathbf{\hat{C}}\_{\lambda}(\mathbf{x}\_{1}) \langle \mathbf{\hat{C}}\_{\tau}(\mathbf{x}\_{2}) \rangle\_{l} + 2 \langle \mathbf{\hat{C}}\_{\lambda}(\mathbf{x}\_{1}) \rangle\_{l} \langle \mathbf{\hat{C}}\_{\tau}(\mathbf{x}\_{2}) \rangle\_{l} \right] \delta\_{l} . \end{split} \tag{12}
$$

Here, *T*˜ is the anti-chronological operator acting on the variables *τ* and *λ* and we used the short-hand notations

$$\int d^4x\_1 = \int d^3x\_1 \int\_{-\infty}^t dt\_1 e^{x(t\_1-t)}, \qquad \hat{X}\_a = e^{-a\hat{A}} \hat{X} e^{a\hat{A}}, \qquad a \in \tau, \lambda. \tag{13}$$

The expansion of Equation (10) implies that the statistical average of any operator *X*ˆ(*x*) can be decomposed into three terms

$$
\langle \hat{X}(\mathbf{x}) \rangle = \langle \hat{X}(\mathbf{x}) \rangle\_l + \langle \hat{X}(\mathbf{x}) \rangle\_1 + \langle \hat{X}(\mathbf{x}) \rangle\_{2\prime} \tag{14}
$$

where the first-order term is given by

$$
\langle \hat{\mathcal{X}}(\mathbf{x}) \rangle\_1 = \int d^4 \mathbf{x}\_1 \left( \hat{\mathcal{X}}(\mathbf{x}), \hat{\mathcal{C}}(\mathbf{x}\_1) \right), \tag{15}
$$

with

$$
\delta \left( \hat{X}(\mathbf{x})\_\prime \hat{Y}(\mathbf{x}\_1) \right) = \int\_0^1 d\tau \langle \hat{X}(\mathbf{x}) \left[ \hat{Y}\_\tau(\mathbf{x}\_1) - \langle \hat{Y}\_\tau(\mathbf{x}\_1) \rangle\_I \right] \rangle\_I \tag{16}
$$

being the two-point correlation function between two arbitrary operators [8,9]. The second-order term in Equation (14) can be written as

$$
\langle \hat{X}(\mathbf{x}) \rangle\_2 = \int d^4 \mathbf{x}\_1 d^4 \mathbf{x}\_2 \Big( \hat{X}(\mathbf{x}), \hat{\mathcal{L}}(\mathbf{x}\_1), \hat{\mathcal{L}}(\mathbf{x}\_2) \Big), \tag{17}
$$

where we introduced the three-point correlation function of the operators *X*ˆ , *Y*ˆ, and *Z*ˆ as

$$
\begin{split}
\langle\hat{\mathcal{X}}(\mathbf{x}),\hat{\mathcal{Y}}(\mathbf{x}\_{1}),\hat{\mathcal{Z}}(\mathbf{x}\_{2})\rangle &= \frac{1}{2} \int\_{0}^{1} d\tau \int\_{0}^{1} d\lambda \langle \hat{T}\hat{\mathcal{X}}(\mathbf{x}) \left[\hat{\mathcal{Y}}\_{\lambda}(\mathbf{x}\_{1})\hat{\mathcal{Z}}\_{\tau}(\mathbf{x}\_{2}) - \langle \hat{T}\hat{\mathcal{Y}}\_{\lambda}(\mathbf{x}\_{1})\hat{\mathcal{Z}}\_{\tau}(\mathbf{x}\_{2})\rangle\_{l} \right. \\ &\left. - \langle \hat{\mathcal{Y}}\_{\lambda}(\mathbf{x}\_{1})\rangle\_{l}\hat{\mathcal{Z}}\_{\tau}(\mathbf{x}\_{2}) - \hat{\mathcal{Y}}\_{\lambda}(\mathbf{x}\_{1})\langle \hat{\mathcal{Z}}\_{\tau}(\mathbf{x}\_{2})\rangle\_{l} + 2\langle \hat{\mathcal{Y}}\_{\lambda}(\mathbf{x}\_{1})\rangle\_{l}\langle \hat{\mathcal{Z}}\_{\tau}(\mathbf{x}\_{2})\rangle\_{l} \right] \}\_{l}.
\end{split} \tag{18}
$$

#### **3. Relativistic Fluid Dynamics at First Order in Gradients**

To examine specific dissipative processes, i.e., heat conduction, particle diffusion, and shear and bulk stresses, the energy-momentum tensor and the particle current are decomposed as

$$\hat{T}^{\mu\nu} \quad = \,\,\hat{\epsilon}\imath\iota^{\mu}\iota^{\nu} - \hat{p}\Delta^{\mu\nu} + \hat{q}^{\mu}\iota^{\nu} + \hat{q}^{\nu}\iota^{\mu} + \hat{\pi}^{\mu\nu},\tag{19}$$

$$
\hat{N}^{\mu}{}\_{\quad} \quad = \quad \hat{n}u^{\mu} + \hat{j}^{\mu}{}\_{\quad} \tag{20}
$$

where the fluid velocity *<sup>u</sup><sup>μ</sup>* is normalized as *<sup>u</sup>μu<sup>μ</sup>* <sup>=</sup> 1, and <sup>Δ</sup>*μν* <sup>=</sup> *<sup>g</sup>μν* <sup>−</sup> *<sup>u</sup>μu<sup>ν</sup>* is the projection operator onto the 3-space orthogonal to *uμ*. The energy-momentum tensor in Equation (19) is assumed to be symmetric with respect to its indices. We remind that *u<sup>μ</sup>* and Δ*μν* appearing in these expressions are classical fields, i.e., *c*-numbers, whereas the rest of the quantities are (microscopic) quantum operators, the physical identification of which is as follows: ˆ is the operator of the energy density; *n*ˆ is the operator of particle (number) density; *p*ˆ is the operator of the pressure; and the dissipative terms *π*ˆ *μν*, *q*ˆ *<sup>μ</sup>*, and ˆ*j <sup>μ</sup>* are the shear-stress tensor, the energy-diffusion flux, and the particle-diffusion flux, respectively.

The tensor decompositions in Equations (19) and (20) have the most general form that can be constructed from the fluid velocity and the tensor Δ*μν*. The operators of the physical quantities on the right-hand sides of Equations (19) and (20) can be written as certain projections of the energy-momentum tensor and the particle current,

$$\hat{\epsilon} = \mu\_{\mu}\mu\_{\nu}\hat{\mathcal{T}}^{\mu\nu}, \qquad \hat{n} = \mu\_{\mu}\hat{\mathcal{N}}^{\mu}, \qquad \hat{p} = -\frac{1}{3}\Delta\_{\mu\nu}\hat{\mathcal{T}}^{\mu\nu}, \tag{21}$$

$$
\hbar^{\mu\nu} = \Delta^{\mu\nu}\_{a\beta} \Upsilon^{a\beta}, \qquad \hbar^{\mu} = u\_a \Delta^{\mu}\_{\beta} \Upsilon^{a\beta}, \qquad \dot{\jmath}^{\dot{\nu}} = \Delta^{\nu}\_{\mu} \hat{N}^{\mu}, \tag{22}
$$

where

$$
\Delta\_{\mu\upsilon\rho\sigma} = \frac{1}{2} \left( \Delta\_{\mu\rho} \Delta\_{\upsilon\sigma} + \Delta\_{\mu\sigma} \Delta\_{\upsilon\rho} \right) - \frac{1}{3} \Delta\_{\mu\upsilon} \Delta\_{\rho\sigma} \tag{23}
$$

is a traceless rank–four projector orthogonal to *uμ*. The dissipative quantities satisfy the conditions

$$
u\_{\nu}\hat{q}^{\nu} = 0, \qquad \boldsymbol{u}\_{\nu}\hat{\boldsymbol{\upbeta}}^{\nu} = 0, \qquad \boldsymbol{u}\_{\nu}\hat{\boldsymbol{\uppi}}^{\mu\nu} = 0, \qquad \hat{\boldsymbol{\pi}}^{\mu}\_{\mu} = 0. \tag{24}$$

In local equilibrium, the averages of these operators vanish [20]:

$$
\langle \langle \mathfrak{f}^{\mu} \rangle\_{l} = 0, \qquad \langle \hat{j}^{\mu} \rangle\_{l} = 0, \qquad \langle \mathfrak{H}^{\mu \nu} \rangle\_{l} = 0,\tag{25}
$$

and one recovers the limit of ideal fluid dynamics. The local-equilibrium pressure is given by the equation of state, i.e., *p*ˆ*<sup>l</sup>* ≡ *p* = *p*(, *n*), which closes the set of ideal fluid-dynamical equations of motion.

Consider next fluid dynamics at first order in gradients. Quite generally, the fluid-dynamical quantities *πμν*, *qμ*, and *j <sup>μ</sup>* are obtained as the statistical averages of the corresponding operators over the NESO according to Equations (14)–(18). Keeping only the first-order terms in Equation (14), we obtain the relativistic Navier–Stokes equations

$$
\pi\_{\mu\nu} = 2\eta\sigma\_{\mu\nu} \qquad \Pi = -\zeta\theta , \qquad \mathcal{J}\_{\mu} = \kappa \left(\frac{nT}{\hbar}\right)^2 \nabla\_{\mu} a\_{\nu} \tag{26}
$$

where *πμν* ≡ *π*ˆ *μν*, the bulk viscous pressure Π ≡ *p*ˆ−*p*ˆ*<sup>l</sup>* is the difference between the first-order average of the pressure operator and the local-equilibrium value of pressure, *h* =  + *p* is the enthalpy density, and

$$
\Box \mathcal{J}\_{\mu} = j\_{\mu} - \frac{n}{h} q\_{\mu} \tag{27}
$$

is the irreversible particle flow, i.e., the particle flow with respect to the energy flow [16,17]. On the right-hand sides of Equation (26), *σμν* <sup>=</sup> *<sup>∂</sup>*<*αuβ*> is the shear tensor, where angular brackets denote the projection with the projector in Equation (23), i.e., *<sup>A</sup>*<*μν*> <sup>=</sup> <sup>Δ</sup>*αβ μνAαβ*, *θ* = *∂μu<sup>μ</sup>* is the expansion scalar, and <sup>∇</sup>*<sup>α</sup>* <sup>=</sup> <sup>Δ</sup>*αβ∂<sup>β</sup>* is the covariant spatial derivative. The coefficients *<sup>η</sup>*, *<sup>ζ</sup>*, and *<sup>κ</sup>* are the transport coefficients of the shear and bulk viscosities, and the thermal conductivity, respectively.

These transport coefficients can be expressed through two-point correlation functions via the following Kubo formulae [8–10]

$$\eta\_{\mu\nu} = -\frac{\beta}{10} \int d^4 \mathbf{x}\_1 \left( \mathfrak{A}\_{\mu\nu}(\mathbf{x}), \mathfrak{A}^{\mu\nu}(\mathbf{x}\_1) \right),\tag{28}$$
 
$$\mathfrak{A}\_{\mu\nu} = \frac{\beta}{10} \int d^4 \mathbf{x}\_1 \left( \mathfrak{A}\_{\mu\nu, \zeta\_1, \zeta\_2, \zeta\_3, \zeta\_4} \right)\tag{29}$$

$$\mathcal{J}\_{\epsilon} = -\oint\_{\mathcal{O}} d^4 \mathbf{x}\_1 \left( \hat{p}^\*(\mathbf{x}), \hat{p}^\*(\mathbf{x}\_1) \right), \tag{29}$$

$$\|\mathbf{x}\|\\_{\ast} = \|-\frac{\beta^2}{3} \int d^4 \mathbf{x}\_1 \left(\hat{h}\_{\mu}(\mathbf{x}), \hat{h}^{\mu}(\mathbf{x}\_1)\right),\tag{30}$$

where

$$
\mathfrak{p}^\* = \mathfrak{p} - \gamma \mathfrak{e} - \delta \mathfrak{h}, \qquad \mathfrak{h}^\mu = \mathfrak{q}^\mu - \frac{h}{n} \mathfrak{j}^\mu,\tag{31}
$$

and

$$\gamma = \left(\frac{\partial p}{\partial \epsilon}\right)\_n , \qquad \delta = \left(\frac{\partial p}{\partial n}\right)\_\epsilon . \tag{32}$$

The correlation functions in Equations (28)–(30) are evaluated in a uniform background, i.e., as if the system was in *global thermodynamical equilibrium*. They can be expressed in terms of the two-point retarded *equilibrium Green functions* as [8,9]

$$\eta\_{\parallel} = -\frac{1}{10} \frac{d}{d\omega} \text{Im} G^{R}\_{\hbar\_{\mu\nu}\hbar^{\mu\nu}}(\omega)\Big|\_{\omega=0} \tag{33}$$

$$\mathcal{L} \quad = \left. -\frac{d}{d\omega} \text{Im} G\_{\vec{p}^\* \vec{p}^\*}^{\text{R}}(\omega) \right|\_{\omega = 0} \, \_{\omega = 0} \, \_{\cdot} \tag{34}$$

$$\text{(\text{x } = \frac{1}{3T} \frac{d}{d\omega} \text{Im} G^R\_{\hat{h}\_{\mu} \hat{h}^{\mu}}(\omega)\Big|\_{\omega = 0} . \tag{35}$$

where, for any two operators *X*ˆ and *Y*ˆ,

$$\mathcal{G}\_{\hat{X}\hat{Y}}^{\mathcal{R}}(\omega) \equiv -i \int\_0^\infty dt e^{i\omega t} \int d^3 \mathbf{x} \langle \left[ \mathcal{X}(\mathbf{x}, t), \hat{Y}(\mathbf{0}, 0) \right] \rangle\_{\mathcal{I}}.\tag{36}$$

Equations (33)–(35) represent a particularly suitable form for the Kubo formulae, which lends itself to evaluation using the methods of equilibrium finite-temperature field theory.

Before closing this section, it is useful to clarify the relation between the expansions in powers of the thermodynamic forces and in powers of the Knudsen number *K* = *l*/*L*, where *l* and *L* are typical microscopic and macroscopic length scales, respectively. To obtain the relations in Equation (26) from Equation (15), we used Curie's theorem. It states that, in an isotropic medium, the correlations between operators of different rank vanish [21]. The integrands in Equations (28)–(30) are mainly concentrated within the range |*x*<sup>1</sup> − *x*| *l*, where *l* is the mean correlation length, which in the weak-coupling limit is of the order of the particle mean free path. The fluid-dynamical regime implies *l L*, where *L* is the typical length scale over which the parameters *β<sup>ν</sup>* and *α* vary in space. Therefore, the thermodynamic forces *∂μβν* and *∂μα* involved in Equation (6) can be factored out from the integral in Equation (15) with their average values at *x*, i.e., the *non-locality* of the thermodynamic forces can be neglected in this approximation. Because <sup>|</sup>*σρσ*||*uρ*|/*L*, the relations in Equation (26) obtained from the gradient expansion in Equation (10) of the NESO are consistent with the expansion scheme in powers of the Knudsen number.

#### **4. Relativistic Fluid Dynamics at Second Order in Gradients**

We have computed systematically all second-order corrections to the dissipative quantities *πμν*, Π, and J *<sup>μ</sup>* based on Equations (14)–(18) [13,14]. Before presenting the results, we note that the second-order contributions arise not only from Equation (17), which is quadratic in the thermodynamic force *C*ˆ, but also from Equation (15), where the non-local nature of the thermodynamic forces in space and time should be carefully taken into account. The non-local effects generate finite relaxation terms in the fluid-dynamical equations, which are required for causality. To see that these corrections are of second order in the Knudsen number, note that they involve the differences of the thermodynamic forces, e.g., *∂μβν*, at the points *x*<sup>1</sup> and *x* (see Equations (6) and (15)). Therefore, we can approximate *<sup>∂</sup>μβν*(*x*1) <sup>−</sup> *<sup>∂</sup>μβν*(*x*) *∂λ∂μβν*(*x*)(*x*<sup>1</sup> <sup>−</sup> *<sup>x</sup>*)*<sup>λ</sup>* <sup>∼</sup> *<sup>K</sup>∂μβν*(*x*), because *<sup>x</sup>*<sup>1</sup> <sup>−</sup> *<sup>x</sup>* <sup>∼</sup> *<sup>l</sup>* and *<sup>∂</sup>* <sup>∼</sup> *<sup>L</sup>*−1, as already done in Section 3. Thus, these corrections contain an additional power of the Knudsen number *K* as compared to the first-order expressions in Equation (26) and, therefore, are of second order.

Here, we restrict ourselves to the second-order expression for the shear-stress tensor and compare it with the results of References [18,22].

#### *4.1. Second-Order Corrections to the Shear-Stress Tensor*

As explained above, we now keep the NESO at second order in small perturbations from local equilibrium and, in addition, we retain terms which are of second order in the gradients of thermodynamic forces. In this manner, we find the shear-stress tensor at second order as

$$
\pi\_{\mu\nu} = 2\eta\sigma\_{\mu\nu} - 2\eta\tau\_{\pi}(\dot{\sigma}\_{\mu\nu} + \gamma\theta\sigma\_{\mu\nu}) + \lambda\_{\pi}\sigma\_{a<\mu}\sigma\_{\nu>}^{a} + 2\lambda\_{\pi\Pi}\theta\sigma\_{\mu\nu} + \lambda\_{\pi\_{\Box}\not\mathcal{J}}\nabla\_{<\mu}a\nabla\_{\nu>}a,\tag{37}
$$

where *<sup>σ</sup>*˙ *μν* <sup>≡</sup> <sup>Δ</sup>*μνρσDσρσ*, with *<sup>D</sup>* <sup>=</sup> *<sup>u</sup>μ∂μ* being the comoving derivative, and *τπ*, *λπ*, *λπ*Π, and *λπ*<sup>J</sup> represent four new coefficients associated with the second-order corrections to the shear stress. The first term on the right-hand side of Equation (37) is easily recognized as the first-order (Navier–Stokes) contribution. The second-order terms collected in the parentheses (i.e., those ∝ *τπ*) represent the non-local corrections to Equation (15), whereas the last three terms stand for the nonlinear corrections arising from the three-point correlation functions in Equation (17). The first non-local correction describes *memory effects* due to its non-locality in time. The relevant transport coefficient *τπ*, which has the dimension of time, measures how long the information remains in the "memory" of the shear-stress tensor *πμν*. Therefore, it is natural to associate it with the relaxation time of the shear stresses towards their asymptotic Navier–Stokes values. The second term involves a product of *σμν* with *θ* = *∂μu<sup>μ</sup>* and can be regarded as a (scalar) measure of the spatial "non-locality" in the fluid-velocity field. This term describes how the shear-stress tensor is distorted by uniform expansion or contraction of the fluid.

We find that the relaxation time *τπ* is related to the frequency derivative of the corresponding first-order transport coefficient, i.e., the shear viscosity, by a Kubo formula

$$\eta \tau \tau\_{\pi} = -i \frac{d}{d\omega} \eta(\omega) \Big|\_{\omega=0} = \frac{1}{10} \frac{d^2}{d\omega^2} \text{Re} G^{R}\_{\hbar\_{\mu\nu}\hbar^{\mu\nu}}(\omega) \Big|\_{\omega=0} \tag{38}$$

where *<sup>η</sup>* <sup>≡</sup> *<sup>η</sup>*(0) is given by Equation (33), the retarded Green's function *<sup>G</sup><sup>R</sup> <sup>π</sup>*<sup>ˆ</sup> *μνπ*<sup>ˆ</sup> *μν* is defined in Equation (36), and the frequency-dependent shear viscosity *η*(*ω*) is given by Equation (33) for non-vanishing *ω*. Similar expressions for the relaxation times were obtained previously in References [18,22–24].

The physical meaning of the Equation (38) for *τπ* is easy to understand. As mentioned above, the relaxation terms originate from the non-local (memory) effects encoded in the non-equilibrium statistical operator. In the case where these memory effects are neglected (first-order theory), the proportionality between *πμν* and *σμν* is given by the zero-frequency (static) limit of the shear viscosity, as seen from Equations (26) and (33). The memory effects imply a time delay, which translates

into a frequency dependence of the shear viscosity [25]. At leading order, this is accounted for by the first-order frequency derivative of *η*(*ω*) as Equation (38) demonstrates.

The last three terms in Equation (37) contain all combinations of the thermodynamic forces *σμν*, *<sup>θ</sup>*, and <sup>∇</sup>*μα* which are allowed by the symmetries to quadratic order. These are *θσμν*, *σρ*<*μσ<sup>ρ</sup> ν*>, and <sup>∇</sup><*μα*∇*ν*>*α*. The second-order transport coefficients associated with these terms can be expressed via three-point correlation functions according to

$$
\lambda\_{\pi} = \frac{12}{35} \beta^2 \int\_{\gamma} d^4 \mathbf{x}\_1 d^4 \mathbf{x}\_2 \left( \mathfrak{H}\_{\mu}^{\nu}(\mathbf{x}), \mathfrak{H}\_{\nu}^{\lambda}(\mathbf{x}\_1), \mathfrak{H}\_{\lambda}^{\mu}(\mathbf{x}\_2) \right), \tag{39}
$$

$$
\lambda\_{\pi\Pi} = -\frac{\beta^2}{5} \int d^4x\_1 d^4x\_2 \left( \mathfrak{H}\_{\mu\nu}(\mathbf{x}), \mathfrak{H}^{\mu\nu}(\mathbf{x}\_1), \mathfrak{p}^\*(\mathbf{x}\_2) \right), \tag{40}
$$

$$
\lambda\_{\pi,\mathcal{J}} = -\frac{1}{5} \int d^4 \mathbf{x}\_1 d^4 \mathbf{x}\_2 \Big( \mathfrak{H}\_{\mu\nu}(\mathbf{x}), \mathcal{J}^{\mu}(\mathbf{x}\_1), \mathcal{J}^{\nu}(\mathbf{x}\_2) \Big), \tag{41}
$$

where Jˆ*<sup>μ</sup>* is the operator corresponding to the 4-current (27). In analogy to the leading-order coefficient *η*, which is given by the two-point correlation of the shear-stress tensor, the second-order coefficient *λπ* is given by the three-point correlation of the shear-stress tensor. The coefficient *λπ*<sup>Π</sup> describes the nonlinear coupling between shear- and bulk-viscous processes and is given by a three-point correlation function between two shear-stress tensors and the bulk viscous pressure. Finally, the coefficient *λπ*<sup>J</sup> describes the nonlinear coupling between the shear and the diffusion processes. Similarly, this coefficient is given by a three-point correlation function between two diffusion currents and the shear-stress tensor. Note that, in Equation (37), the term ∝ *λπ*<sup>Π</sup> and the second term in parenthesis have the same gradient structure, but they have different origins and physical interpretation. The term ∝ *τπ* originates from *non-local* effects in the statistical distribution, whereas the term ∝ *λπ*<sup>Π</sup> stands purely for the *nonlinear* coupling between the bulk- and the shear-viscous effects. In this sense, it is natural to regard as nonlinear only the term ∝ *λπ*Π, but not the term ∝ *τπ*. A similar classification of the second-order terms was suggested earlier in Reference [19].

#### *4.2. Comparison with Previous Studies*

For the sake of simplicity we consider here a fluid without conserved charges. In this case, Equation (32) implies *<sup>γ</sup>* <sup>≡</sup> *<sup>c</sup>*<sup>2</sup> *<sup>s</sup>*, where *cs* is the speed of sound, and Equation (37) reduces to

$$
\pi\_{\mu\nu} = 2\eta\sigma\_{\mu\nu} - 2\eta\tau\_{\pi}(\dot{\sigma}\_{\mu\nu} + c\_s^2 \theta \sigma\_{\mu\nu}) + \lambda\_{\pi}\sigma\_{\pi}\iota\_{\gamma}\sigma\_{\nu\geq\ast}^a + 2\lambda\_{\pi\Pi}\theta\sigma\_{\mu\nu}.\tag{42}
$$

Baier et al. [18] found in this case and for conformal fluids

$$
\Delta \pi\_{\mu\nu}^{B} = 2\eta \sigma\_{\mu\nu} - 2\eta \pi\_{\pi} \left( \sigma\_{\mu\nu} + \frac{1}{3} \theta \sigma\_{\mu\nu} \right) + \lambda\_1 \sigma\_{\kappa < \mu} \sigma\_{\nu > \nu}^{a} \tag{43}
$$

where we have neglected terms proportional to the vorticity tensor *wαβ* = (∇*αu<sup>β</sup>* − ∇*βuα*)/2. (Note that Baier et al. [18] and Romatschke [22] used a metric convention opposite to ours, and their definition of the shear viscosity differs from ours by a factor of 2). Because *c*<sup>2</sup> *<sup>s</sup>* = 1/3 for a conformal fluid, we recover from Equation (42) the term proportional to *τπ* in Equation (43). Furthermore, because conformal invariance implies a vanishing bulk viscous pressure, the correlations involving the operator *p*ˆ<sup>∗</sup> (see Equations (29) and (31)) vanish, i.e., *λπ*<sup>Π</sup> = 0 in this case. Finally we see that *λ*<sup>1</sup> ≡ *λπ*.

In the case of non-conformal fluids, the second-order expression for the shear-stress tensor was found, e.g., in Reference [22] in the absence of conserved charges. Again, neglecting the vorticity tensor and assuming flat space-time,

$$
\Delta \boldsymbol{\pi}\_{\mu\nu}^{R} = 2\eta \sigma\_{\mu\nu} - 2\eta \boldsymbol{\tau}\_{\pi} \left( \boldsymbol{\sigma}\_{\mu\nu} + \frac{1}{3} \theta \boldsymbol{\sigma}\_{\mu\nu} \right) + \lambda\_1 \boldsymbol{\sigma}\_{a<\mu} \boldsymbol{\sigma}\_{\nu>}^{a} - \frac{2}{3} \eta \boldsymbol{\tau}\_{\pi}^{\*} \theta \boldsymbol{\sigma}\_{\mu\nu} + \lambda\_4 \boldsymbol{\nabla}\_{<\mu} \ln \boldsymbol{s} \boldsymbol{\nabla}\_{\nu>} \ln \operatorname{s.} \tag{44}$$

The term ∝ *τ*∗ *<sup>π</sup>* has the same gradient structure as the non-local term −2*ητπθσμν*/3. Comparing Equation (44) with our expression in Equation (42), we identify *τ*∗ *<sup>π</sup>* = *τπ*(3*c*<sup>2</sup> *<sup>s</sup>* −1) −3*λπ*Π/*η*, and *λ*<sup>4</sup> = 0.

We also note that Equation (37) does not contain terms proportional to the vorticity. To derive such terms, one needs to include an initial non-zero angular momentum in the local-equilibrium distribution [26].

#### *4.3. Relaxation Equation for the Shear-Stress Tensor*

A relaxation-type equation for *πμν* can now be derived from Equation (37). For this purpose, we replace the first-order expression 2*σρσ* <sup>→</sup> *<sup>η</sup>*−1*πρσ* in the second term on the right-hand-side of Equation (37), as has also been done in References [18,27,28]. This substitution is justified up to second order in space-time gradients. We then obtain

$$-2\eta \tau\_{\pi} \dot{\sigma}\_{\mu \nu} \simeq -\tau\_{\pi} \dot{\pi}\_{\mu \nu} + \tau\_{\pi} \beta \eta^{-1} \left(\gamma \frac{\partial \eta}{\partial \beta} - \delta \frac{\partial \eta}{\partial a}\right) \theta \pi\_{\mu \nu} \simeq -\tau\_{\pi} \dot{\pi}\_{\mu \nu} + 2\tau\_{\pi} \beta \left(\gamma \frac{\partial \eta}{\partial \beta} - \delta \frac{\partial \eta}{\partial a}\right) \theta \sigma\_{\mu \nu}, \tag{45}$$

where *π*˙ *μν* = Δ*μνρσDπρσ*. The terms in brackets contain the corresponding partial derivatives of *η*, which in general are not small and should not be neglected. In Equation (45), we employed the relations *Dβ* = *βθγ* and *Dα* = −*βθδ* [9]. Combining Equations (37) and (45) and introducing the coefficient

$$
\lambda = \lambda\_{\pi \text{II}} - \gamma \eta \tau\_{\pi} + \tau\_{\pi} \beta \left( \gamma \frac{\partial \eta}{\partial \beta} - \delta \frac{\partial \eta}{\partial a} \right), \tag{46}
$$

we finally obtain

$$
\pi\_{\pi}\pi\_{\mu\nu} + \pi\_{\mu\nu} = 2\eta\sigma\_{\mu\nu} + 2\lambda\theta\sigma\_{\mu\nu} + \lambda\_{\pi}\sigma\_{\rho<\mu}\sigma\_{\nu>}^{\rho} + \lambda\_{\pi\_{\Box}\mathcal{J}}\nabla\_{\circ\langle\mu}\mathbf{a}\nabla\_{\nu>\mu}\mathbf{a}.\tag{47}
$$

The time-derivative term on the left-hand side describes the relaxation of the shear-stress tensor towards its Navier–Stokes value on the characteristic time scale *τπ*. Indeed, for vanishing right-hand side the relaxation is exponential, *πμν* ∝ exp(−*t*/*τπ*), with a characteristic relaxation time scale *τπ*. We would like to stress that the exponential relaxation over a time scale *τπ* is a direct consequence of the substitution 2*σρσ* <sup>→</sup> *<sup>η</sup>*−1*πρσ* made above; it is not a manifestation of a specific relaxation process on the microscopic level. However, a direct way to obtain such a relaxation term is via the method of moments applied to the Boltzmann equation [3]. In this case, the time scale *τπ* is an intrinsic property of the collision kernel in the Boltzmann equation.

#### **5. Summary**

This work concisely presents the derivation of second-order relativistic dissipative fluid dynamics within Zubarev's NESO formalism – a method which is well-suited for treatments of strongly correlated systems. The simple case of a one-component fluid without electromagnetic fields or vorticity in flat space-time was considered here.

Our analysis shows that the second-order dissipative terms arise from: (i) the quadratic terms in the Taylor expansion of the statistical operator; and (ii) the linear terms of the same expansion which include memory and non-locality in space. In particular, we find that the type-(ii) terms describe the relaxation in time of the dissipative fluxes, which is essential for the causality of the fluid-dynamical theory.

Using the NESO method and the example of the shear-stress tensor, we demonstrated that the second-order transport coefficients can be expressed in terms of certain two- and three-point equilibrium correlation functions. A discussion of the transport coefficients associated with other thermodynamic fluxes can be found elsewhere [14]. Furthermore, we have shown that Kubo-type formulae for the relaxation times of the dissipative fluxes can be obtained within the NESO formalism (see Equation (38)). These are given by the zero-frequency limit of the derivatives of the corresponding first-order transport coefficients with respect to the frequency. These can be computed from the theory of quantum fields in equilibrium at non-zero temperature as, for example, was done by us for the QGP within the Nambu–Jona–Lasinio model [29,30].

**Author Contributions:** Conceptualization, A.H., A.S., D.H.R.; Methodology, A.H., A.S., D.H.R.; Investigation, A.H., A.S., D.H.R.; Writing—Original Draft Preparation, A.H., A.S., D.H.R.; Writing—Review & Editing, A.H., A.S., D.H.R.; Funding Acquisition, A.S., D.H.R.

**Funding:** This research was funded by HGS-HIRe graduate program at Goethe University, Deutsche Forschungsgemeinschaft (Germany) through grants SE 1836/4-1 and CRC-TR 211, the State Administration of Foreign Experts Affairs (China) through Expert project GDW20167100136.

**Acknowledgments:** We are grateful to Ulrich Heinz and Xu-Guang Huang for discussions. A.H. acknowledges support from the HGS-HIRe graduate program at Goethe University. A.S. was supported by the Deutsche Forschungsgemeinschaft (Grant No. SE 1836/4-1). D.H.R. acknowledges support by the High-End Visiting Expert project GDW20167100136 of the State Administration of Foreign Experts Affairs (SAFEA) of China and by the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC-TR 211 "Strong-interaction matter under extreme conditions".

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
