*Article* **Reworking Zubarev's Approach to Nonequilibrium Quantum Statistical Mechanics**

**Francesco Becattini 1,\*, Matteo Buzzegoli <sup>1</sup> and Eduardo Grossi <sup>2</sup>**


Received: 3 February 2019; Accepted: 19 March 2019; Published: 8 April 2019

**Abstract:** In this work, the nonequilibrium density operator approach introduced by Zubarev more than 50 years ago to describe quantum systems at a local thermodynamic equilibrium is revisited. This method, which was used to obtain the first "Kubo" formula of shear viscosity, is especially suitable to describe quantum effects in fluids. This feature makes it a viable tool to describe the physics of Quark–Gluon Plasma in relativistic nuclear collisions.

**Keywords:** non-equilibrium statistical operator; quantum statistical mechanics; relativistic hydrodynamics; Kubo formulae

#### **1. Introduction**

One of the authors (F.B.) would like to start this paper with a personal recollection. I first ran across Zubarev's papers when I was studying the derivation by A. Hosoya et al. [1] of the shear viscosity in quantum-field theory, a result widely known as "Kubo formula", like many of the same sort. This derivation was overtly based on Zubarev's method of nonequilibrium density (or statistical) operator, and I surmised that this method must have been a very important and renowned tool in quantum statistical mechanics. In fact, surprisingly, it could hardly be found in textbooks and in the recent literature, and I did not quite understand why the founding method of such an important formula was that overlooked. After some more self-education, I realized that, perhaps part of the problem was that Zubarev himself did not put the right emphasis on the crucial feature that his proposed operator should possess: to be stationary, hence well-suited to be used in relativistic quantum-field theory as a density operator in the Heisenberg representation. A nonequilibrium stationary density operator sounds somewhat contradictory, but this is not the case if we deal with a system that, at some time, is known to be in local thermodynamic equilibrium, as we see in more detail in Section 3.

In this work, we would like to not just summarize Zubarev's method [2–4], but also to make a critical appraisal and to provide a reformulation thereof that highlights the important features of this approach in a hopefully clear fashion. I also hope that this work could do justice to Zubarev and his remarkable achievement.

#### *Notation*

In this paper we use natural units, with ¯*h* = *c* = *K* = 1.

The Minkowskian metric tensor is diag(1, −1, −1, −1); for the Levi-Civita symbol, we use convention <sup>0123</sup> = 1.

Operators in the Hilbert space are denoted by a large upper hat, e.g., *T*", while unit vectors with a small upper hat, e.g., *v*ˆ. Scalar products and contractions are sometimes denoted with a dot, e.g., *<sup>A</sup>μB<sup>μ</sup>* <sup>=</sup> *<sup>A</sup>* · *<sup>B</sup>*.

#### **2. Local Thermodynamic Equilibrium**

Zubarev formalism can be used in nonrelativistic as well as in relativistic quantum statistical mechanics. We can then start from the latter, more general case, which is applicable to relativistic fluids out of equilibrium [5]. The relativistic version of the nonequilibrium density operator was first put forward by Zubarev himself and his collaborators in 1979 [6], and later reworked by Van Weert in Reference [7].

The starting point is the definition of the local equilibrium density operator. In relativity, this notion needs the specification of a one-parameter family of 3D spacelike hypersurfaces Σ(*τ*) (see Figure 1), also known as foliation of spacetime [6–9]. "Time" *τ* does not necessarily coincide with the proper time marked by comoving clocks. Local equilibrium density operator *<sup>ρ</sup>*"LE is obtained by maximizing the total entropy:

$$S = -\text{tr}(\hat{\rho}\log(\hat{\rho})) \tag{1}$$

with constrained values of energy–momentum and charge density, which should be equal to the actual values. In a covariant formulation, these densities are obtained by projecting the mean values of the stress–energy tensor and current onto the normalized vector perpendicular to Σ:

$$n\_{\mu}\text{tr}\left(\hat{\rho}^{\hat{T}^{\mu\nu}}\right) = n\_{\mu}T^{\mu\nu}, \qquad n\_{\mu}\text{tr}\left(\hat{\rho}^{\hat{\jmath}^{\mu}}\right) = n\_{\mu}j^{\mu}.\tag{2}$$

where *Tμν* and *j <sup>μ</sup>* are the true values of the stress–energy and current fields. The operators in Equation (2) are in the Heisenberg representation. In addition to the energy, momentum, and charge densities, one should include the angular momentum density, but if the stress–energy tensor is the Belinfante [10], this further constraint is redundant and can be disregarded.

The resulting operator is the Local Equilibrium Density Operator (LEDO):

$$\hat{\rho}\_{\rm LE} = \frac{1}{Z\_{\rm LE}} \exp\left[-\int\_{\Sigma(\mathbf{r})} \mathbf{d}\Sigma \, n\_{\mu} \left(\hat{T}^{\mu\nu}(\mathbf{x}) \beta\_{\nu}(\mathbf{x}) - \zeta(\mathbf{x}) \hat{j}^{\mu}(\mathbf{x})\right)\right] \tag{3}$$

where *β* and *ζ* are the relevant Lagrange multiplier functions for this problem, whose meaning is the four-temperature vector and the ratio between local chemical potential and temperature, respectively [8], dΣ is the measure of the hypersurface induced by the Minkowskian metric, and fields *<sup>β</sup>* and *<sup>ζ</sup>* are the solution of Constraints (2) with *<sup>ρ</sup>*" <sup>=</sup> *<sup>ρ</sup>*"LE, namely:

$$n\_{\boldsymbol{\mu}}\text{tr}\left(\widehat{\rho}\_{\text{LE}}\,\widehat{T}^{\mu\nu}\right) = n\_{\boldsymbol{\mu}}T^{\mu\nu}\_{\text{LE}}[\beta\_{\prime}\zeta\_{\prime}n] = n\_{\boldsymbol{\mu}}T^{\mu\nu}, \qquad n\_{\boldsymbol{\mu}}\text{tr}\left(\widehat{\rho}\_{\text{LE}}\,\widehat{j}^{\mu}\right) = n\_{\boldsymbol{\mu}}j^{\mu}\_{\text{LE}}[\beta\_{\prime}\zeta\_{\prime}n] = n\_{\boldsymbol{\mu}}j^{\mu}.\tag{4}$$

These equations indeed define a vector field *β*, which, in turn, can be used as a hydrodynamic frame, *β* [8] or thermodynamic frame [11], by identifying the four-velocity with

$$
\mu = \frac{\beta}{\sqrt{\beta^2}} = T\beta,\tag{5}
$$

which somehow inverts the usual definition.

It is important to stress that the LEDO in Equation (3) is not stationary because operators are generally time-dependent. The sufficient condition for stationarity is that *β* is a killing vector field, and *ζ* a constant; in this case, the LEDO becomes the general global thermodynamic equilibrium operator [12].

#### **3. Nonequilibrium Density Operator Revisited**

The true density operator in the Heisenberg representation must be stationary by definition, whereas the LEDO is not. The solution of how to work it out (which is an amendment of Zubarev's original idea) is overly simple: if, at some initial time *τ*0, the system is known to be in local

thermodynamic equilibrium, the actual, stationary, nonequilibrium density operator (NEDO) is *<sup>ρ</sup>*"LE(*τ*0). Therefore, the true mean values of quantum operators should be calculated as:

$$\langle \hat{O} \rangle \equiv \text{tr}(\hat{\rho}\hat{O}) = \text{tr}(\hat{\rho}\_{\text{LE}}(\pi\_0)\hat{O})$$

One can rewrite *<sup>ρ</sup>*"LE(*τ*0) in terms of the operators at present "time" *<sup>τ</sup>* by means of Gauss' theorem, taking into account that *T*" and "*j* are conserved. Defining

$$\mathrm{d}\Sigma\_{\mu} = \mathrm{d}\Sigma \, n\_{\mu\nu}$$

and dΩ being the measure of a 4D region in spacetime, we have

$$-\int\_{\Sigma(\tau\_0)} \mathrm{d}\Sigma\_{\mu} \left( \widehat{T}^{\mu\nu} \, {\beta\_{\nu}} - \widehat{j}^{\mu} \zeta \right) = -\int\_{\Sigma(\tau')} \mathrm{d}\Sigma\_{\mu} \left( \widehat{T}^{\mu\nu} \, {\beta\_{\nu}} - \widehat{j}^{\mu} \zeta \right) + \int\_{\Omega} \mathrm{d}\Omega \, \left( \widehat{T}^{\mu\nu} \nabla\_{\mu} \beta\_{\nu} - \widehat{j}^{\mu} \nabla\_{\mu} \zeta \right), \tag{6}$$

where ∇ is the covariant derivative. Region Ω is the portion of spacetime enclosed by two hypersurfaces Σ(*τ*0) and Σ(*τ*) and the timelike hypersurface at their boundaries, where the flux of (*T*"*μνβν*(*x*) <sup>−</sup> "*<sup>j</sup> μζ*(*x*)) is supposed to vanish (see Figure 1). Consequently, the stationary NEDO reads:

$$\hat{\rho} = \frac{1}{Z} \exp\left[ - \int\_{\Sigma(\tau\_0)} d\Sigma\_{\mu} \left( \hat{T}^{\mu\nu} \beta\_{\nu} - \hat{\beta}^{\mu} \zeta \right) \right] = \frac{1}{Z} \exp\left[ - \int\_{\Sigma(\tau)} d\Sigma\_{\mu} \left( \hat{T}^{\mu\nu} \beta\_{\nu} - \hat{\beta}^{\mu} \zeta \right) + \int\_{\Omega} d\Omega \, \left( \hat{T}^{\mu\nu} \nabla\_{\mu} \beta\_{\nu} - \hat{\beta}^{\mu} \nabla\_{\mu} \zeta \right) \right] \tag{7}$$

This expression is the generally covariant form of the one used in Reference [1] (Equation (2.9) therein), with the only difference that factor exp[*ε*(*t* − *τ*)] does not appear in the second term. In Section 5, we see that such a factor is not necessary to obtain the correct "Kubo" formulae.

**Figure 1.** Spacelike hypersurfaces Σ(*τ*), Σ(*τ*0) and their normal unit vector *n* defining local thermodynamical equilibrium for a relativistic fluid in Minkwoski spacetime. At timelike boundary Σ*l*, the flux is supposed to vanish.

The NEDO can be worked out perturbatively by identifying the two terms in the exponent of Equation (7):

$$\hat{A} = -\int\_{\Sigma(\tau)} \mathbf{d}\Sigma\_{\mu} \left( \hat{T}^{\mu\nu} \beta\_{\nu} - \hat{j}^{\mu} \mathcal{J}\_{\nu} \right) \tag{8}$$

and

$$
\widehat{B} = \int\_{\Omega} \mathrm{d}\Omega \,\, \left( \widehat{T}^{\mu\nu} \nabla\_{\mu} \beta\_{\nu} - \widehat{j}^{\mu} \nabla\_{\mu} \zeta \right) , \tag{9}
$$

and assuming that *B*" is small compared to *A*"; this happens if the system has small correlation length, and if the gradients in Equation (9) are small, which is the hydrodynamic limit. We can then use identity

$$\exp[\hat{A} + \hat{B}] = \exp[\hat{A}] + \int\_0^1 \mathrm{d}z \, \exp[z(\hat{A} + \hat{B})] \hat{B} \exp[-z\hat{A}] \exp[\hat{A}]$$

The expansion of exp[*A*" + *B*"] can be iterated in the integrand, and one obtains an operator expansion in *B*". Taking into account that

$$Z = \text{tr}(\exp[\bar{A} + \hat{B}])$$

at the lowest order in *B*" (linear response):

$$
\hat{\rho} \simeq \hat{\rho}\_{\text{LE}} + \int\_0^1 \text{d}z \, \exp[z\hat{A}] \hat{B} \exp[-z\hat{A}] \hat{\rho}\_{\text{LE}} - \langle \hat{B} \rangle\_{\text{LE}} \hat{\rho}\_{\text{LE}} \tag{10}
$$

which is the starting point to obtain the "Kubo" formulae.

It should be pointed out that the original Zubarev formulae were somewhat different [6]. We work it out by using Cartesian coordinates and hyperplanes as hypersurfaces. Zubarev modified the equation for the NEDO in Heisenberg representation

$$\frac{d\vec{\rho}}{dt} = -\epsilon(\hat{\rho} - \hat{\rho}\_{\text{LE}}) \tag{11}$$

*<sup>ε</sup>* <sup>&</sup>gt; 0 being a real parameter whose limit *<sup>ε</sup>* <sup>→</sup> 0 is to be taken *after* the thermodynamic limit. The solution of the above equation at the present time, which can be chosen to be *t* = 0, reads:

$$
\hat{\rho}(0) = \hat{\rho}\_{\rm LE} - \int\_{-\infty}^{0} \mathbf{d}t \,\mathrm{e}^{\rm t} \frac{\mathbf{d}\hat{\rho}\_{\rm LE}}{\mathbf{d}t} \tag{12}
$$

One can now use the general expression for the derivative of an exponential to calculate:

$$\frac{\mathbf{de}^{\hat{A}}}{\mathbf{d}t} = \int\_0^1 \mathbf{d}z \,\, \mathbf{e}^{z\hat{A}} \frac{\mathbf{d}\hat{A}}{\mathbf{d}t} \mathbf{e}^{(1-z)\hat{A}}$$

with *A*" given by Equation (8). This implies

$$\frac{dZ\_{\rm LE}}{dt} = \frac{d}{dt} \text{tr}(\mathbf{e}^{\hat{A}}) = \text{tr}\left(\frac{\mathbf{d}\hat{A}}{\mathbf{d}t} \mathbf{e}^{\hat{A}}\right) = Z\_{\rm LE} \langle \frac{\mathbf{d}\hat{A}}{\mathbf{d}t} \rangle\_{\rm LE}$$

so that

$$\frac{d\hat{\rho}\_{\rm LE}}{dt} = \int\_0^1 \mathrm{d}z \,\mathrm{e}^{z\tilde{A}} \frac{\mathrm{d}\ddot{A}}{\mathrm{d}t} \mathrm{e}^{-z\tilde{A}} \hat{\rho}\_{\rm LE} - \langle \frac{\mathrm{d}\ddot{A}}{\mathrm{d}t} \rangle\_{\rm LE} \hat{\rho}\_{\rm LE} \tag{13}$$

If the surface boundary terms vanish, we have

$$\frac{d\hat{A}}{dt} = -\int \mathrm{d}^3 \mathbf{x} \, \frac{\partial}{\partial t} (\hat{T}^{0\nu} \beta\_{\nu}) = -\int \mathrm{d}^3 \mathbf{x} \, \partial\_{\mu} (\hat{T}^{\mu\nu} \beta\_{\nu}) = -\int \mathrm{d}^3 \mathbf{x} \, \hat{T}^{\mu\nu} \partial\_{\mu} \beta\_{\nu} \tag{14}$$

By plugging Equations (13) and (14) into Equation (12), we have:

$$
\hat{\rho}(0) - \hat{\rho}\_{\rm LE} = \int\_0^1 \text{d}z \,\text{e}^{z\overline{A}} \int\_{-\infty}^0 \text{d}^4 x \,\text{e}^{tt} \hat{\mathcal{I}}^{\mu\nu} \partial\_{\mu} \beta\_{\nu} \text{e}^{-z\overline{A}} \,\hat{\rho}\_{\rm LE} - \int\_{-\infty}^0 \text{d}^4 x \,\text{e}^{tt} \langle \hat{\mathcal{I}}^{\mu\nu} \rangle\_{\rm LE} \partial\_{\mu} \beta\_{\nu} \,\hat{\rho}\_{\rm LE} \,\text{e}^{zz}
$$

Taking Equation (9) into account, the above equation is basically linear Approximation (10) with an extra factor exp(*εt*) in the integrand. In a sense, Zubarev Assumption (11) of a small source term in the density operator evolution equation in the Heisenberg representation leads to the linear approximation of the fully stationary density operator operator (7). However, it should be emphasized that such an extra factor is not necessary. The Heisenberg equation for the true density operator is d*ρ*"/d*<sup>t</sup>* <sup>=</sup> 0 does not need any modification for the derivation of the Kubo formulae or any other result depending on local thermodynamic equilibrium, as we will show in Section 5. A fully relativistic viewpoint with the application of the Gauss theorem makes the derivation of the NEDO expression (7) straightforward, transparent and simple.

#### **4. Entropy Production**

A remarkable consequence of this approach is the derivation of a general equation for the entropy production rate, which was reported in References [6,7]. Let us start with the assumption that *S* is an integral of entropy current *sμ*:

$$S(\tau) = -\text{tr}(\hat{\rho}\_{\text{LE}}(\tau) \log \hat{\rho}\_{\text{LE}}(\tau)) = \int\_{\Sigma(\tau)} \text{d}\Sigma\_{\mu} \text{ s}^{\mu}$$

On the other hand, entropy can be expanded by using Equation (3):

$$S(\tau) = -\text{tr}(\hat{\rho}\_{\text{LE}}(\tau) \log \hat{\rho}\_{\text{LE}}(\tau)) = \log Z\_{\text{LE}} + \int\_{\Sigma(\tau)} \text{d}\Sigma\_{\mu} \, \langle \hat{T}^{\mu\nu} \rangle\_{\text{LE}} \beta\_{\nu} - \zeta \langle \hat{j}^{\mu} \rangle\_{\text{LE}}$$

$$= \log Z\_{\text{LE}} + \int\_{\Sigma(\tau)} \text{d}\Sigma\_{\mu} \, \left( T^{\mu\nu} \beta\_{\nu} - \zeta j^{\mu} \right) \tag{15}$$

where we have used Constraints (2), taking into account that dΣ*<sup>μ</sup>* = dΣ *nμ*.

The derivative with respect to *τ* can be computed by taking advantage of a general expression for the variation of an integral between two infinitesimally closed hypersurfaces:

$$\frac{d\mathbf{S}}{d\tau} = \int\_{\Sigma(\tau)} \mathbf{d}\Sigma(\mathbf{n} \cdot \mathcal{U}) \nabla \cdot \mathbf{s} + \frac{1}{2} \int\_{\partial\Sigma(\tau)} \mathbf{d}\mathcal{S}\_{\mu\nu} (\mathbf{s}^{\mu} \mathcal{U}^{\nu} - \mathbf{s}^{\nu} \mathcal{U}^{\mu}) \tag{16}$$

where *∂*Σ is the 2D boundary of Σ and *U<sup>μ</sup>* = *∂xμ*/*∂τ*; *S*˜ is the dual of the surface element. Now, assume that the boundary term does not contribute, and calculate the same derivative by using Expression (15):

$$\begin{split} \frac{d\mathbf{S}}{d\tau} &= \frac{d\log Z\_{\rm LE}}{d\tau} + \int\_{\Sigma(\tau)} \mathbf{d}\Sigma(n \cdot \mathcal{U}) \nabla\_{\mu} \left( T^{\mu\nu} \beta\_{\nu} - \zeta j^{\mu} \right) \\ &= \frac{d\log Z\_{\rm LE}}{d\tau} + \int\_{\Sigma(\tau)} \mathbf{d}\Sigma(n \cdot \mathcal{U}) T^{\mu\nu} \nabla\_{\mu} \beta\_{\nu} - j^{\mu} \nabla\_{\mu} \zeta \end{split} \tag{17}$$

where we have taken advantage of the conservation of the *exact* values *Tμν* and *j <sup>μ</sup>*. The remaining task is to calculate the derivative of log *Z*LE, which can be done by using its definition

$$\frac{\operatorname{d}\log Z\_{\rm LE}}{\operatorname{d}\tau} = \frac{1}{Z\_{\rm LE}} \frac{\operatorname{d}}{\operatorname{d}\tau} \operatorname{tr}(\exp[\widehat{A}]),$$

with *A*" in Equation (8). By using the same formula as the derivative of a *τ*-dependent integral in Equation (17), and assuming that the boundary term vanishes:

$$\hat{\mathcal{V}}\_{\text{L\text{\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\text\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\}\/\}\}\/\}\}\/\}\}\/\}\}\/\}\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\}\/\/\}\/\/\}\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/)\tag\text{!}\text{ }\text{\textquotedbl0\text{ }\text{\textquend\text\text$$

By plugging Equation (18) into Equation (17) and comparing with Equation (16), taking into account that the equation should hold for any *τ*, we have

$$\nabla \cdot \mathbf{s} = \left( T^{\mu \nu} - T\_{\text{LE}}^{\mu \nu} \right) \nabla\_{\mu} \beta\_{\nu} - \left( j^{\mu} - j\_{\text{LE}}^{\mu} \right) \nabla\_{\mu} \zeta\_{\nu} \tag{19}$$

which was found in Reference [6], and tells us that the deviations of the conserved currents actual values from those at the local thermodynamic equilibrium are responsible for entropy production.

#### **5. Kubo Formulae**

Let us now apply the expansion of the NEDO (10) to calculate the mean value of a local operator *O*" at present time *t*:

$$
\langle \hat{O}(\mathbf{x}) \rangle \simeq \langle \hat{O}(\mathbf{x}) \rangle\_{\rm LE} - \langle \hat{O}(\mathbf{x}) \rangle\_{\rm LE} \langle \hat{B} \rangle\_{\rm LE} + \int\_0^1 \mathrm{d}z \,\langle \hat{O}(\mathbf{x}) \mathbf{e}^{z\hat{A}} \hat{B} \mathbf{e}^{-z\hat{A}} \rangle\_{\rm LE} \tag{20}
$$

where *A*" and *B*" are in Equations (8) and (9), respectively. To work out Formula (20), it is customary to approximate the *A*" in the *z* integral on the right-hand side with the global equilibrium expression. In a covariant fashion, this means making a zero-order approximation of the Taylor expansion of the thermodynamic fields from point *x*, where operator *O*" is to be calculated:

$$\hat{A} = -\int\_{\Sigma(\tau)} d\Sigma\_{\mu} \left( \hat{T}^{\mu\nu} \beta\_{\nu} - \hat{\beta}^{\mu} \zeta \right) \simeq -\beta\_{\nu}(\tau, \sigma) \int\_{\Sigma(\tau)} d\Sigma\_{\mu} \, \hat{T}^{\mu\nu} + \zeta(\tau, \sigma) \int\_{\Sigma(t)} d\Sigma\_{\mu} \, \hat{\beta}^{\mu} = -\beta\_{\nu}(\mathbf{x}) \hat{P}^{\nu} + \zeta(\mathbf{x}) \hat{Q} \tag{21}$$

where *P*" is the total four-momentum and *Q*" the total charge. Hence,

$$
\hat{\rho}\_{\rm LE} \simeq \frac{1}{Z\_{\rm LE}} \exp[\hat{A}] \simeq \frac{1}{Z} \exp[-\beta(\mathbf{x}) \cdot \hat{P} + \zeta(\mathbf{x}) \hat{Q}] \equiv \hat{\rho}\_{\rm eq(\mathbf{x})} \tag{22}
$$

that is, *<sup>ρ</sup>*"eq(x) is the global equilibrium density operator having the same vector at point *<sup>x</sup>* and similarly for *ζ* as constant inverse temperature four-vector.

Furthermore, we replace the integration region enclosed by the two LTE hypersurfaces at *t* and *t*<sup>0</sup> with the spacelike tangent hyperplanes at points *x* = (*τ*, *σ*) and *x*<sup>0</sup> = (*τ*0, *σ*), respectively, whose normal versor is *n*. This allows to carry out the integration over Minkowski spacetime by using Cartesian coordinates, that is, time *t* marked by an observer moving with velocity *n*, and a vector of coordinates **x** for the hyperplanes. These approximations make it possible to replace covariant derivatives with usual partial derivatives in Cartesian coordinates:

$$\int\_{\Omega} \mathrm{d}\Omega \,\left(\hat{T}^{\mu\nu}\nabla\_{\mu}\beta\_{\nu} - \hat{\beta}^{\mu}\nabla\_{\mu}\zeta\right) \to \int\_{T\Omega} \mathrm{d}^{4}\mathbf{x} \,\left(\hat{T}^{\mu\nu}\partial\_{\mu}\beta\_{\nu} - \hat{\beta}^{\mu}\partial\_{\mu}\zeta\right) \tag{23}$$

where *T*Ω is the region encompassed by the two hyperplanes. Thereby, and provided that *n*(*x*) = *β*ˆ(*x*), that is, that the local equilibrium hypersurface is locally normal to the flow velocity defined by the four-temperature vector [8], Formula (20) can be turned into a more manageable one (see Appendix A for a summary of the derivation) involving the commutators of operator *O*", with the stress–energy tensor and current operators

$$
\langle \langle \hat{O}(\mathbf{x}) \rangle - \langle \hat{O}(\mathbf{x}) \rangle\_{\rm LE} \simeq iT \int\_{I\_0^\mathbf{d}} \mathbf{d}^4 \mathbf{x}' \int\_{I\_0}^{\mathbf{r}'} \mathbf{d} \theta \left( \langle [\hat{O}(\mathbf{x}), \hat{T}^{\mu\nu}(\theta, \mathbf{x}')] \rangle\_{\hat{\mathbb{P}}(\mathbf{r})} \mathfrak{d}\_{\mu} \mathfrak{E}\_{\mathbf{v}}(\mathbf{x}') - \langle [\hat{O}(\mathbf{x}), \hat{f}^{\mu}(\theta, \mathbf{x}')] \rangle\_{\hat{\mathbb{P}}(\mathbf{r})} \mathfrak{d}\_{\mu} \mathfrak{J}(\mathbf{x}') \right) \tag{24}
$$

where *T* = 1/ *β*2, and subscript *β*(*x*) stands for averaging with the density operator in Equation (22). It is important to stress the different time arguments for the operators and the thermodynamic fields in Equation (24).

From Equation (24), it turns out that the deviation from LTE of the mean value of *O*" at any time depends on the whole history of thermodynamic fields *β* and *ζ*. However, the correlation length between *O*"(*x*) and both *T*"(*x* ), "*j*(*x* ) is typically much smaller than the distance over which the gradients

of *β* and *ζ* have significant variations. This statement amounts to assuming a separation between the typical microscopic interaction scale and the macroscopic hydrodynamical scale. One would then be tempted to take the gradients out of the integral in Equation (24). However, much care should be taken in this because the derivation of Formula (24), more precisely nonequilibrium density operator (7), required the vanishing of the flux of *<sup>T</sup>*"*μνβν* <sup>−</sup> "*<sup>j</sup> μζ* at the boundary timelike hypersurface. If one expands the perturbation of the thermodynamic fields with respect to their equilibrium value, by definition those at point *x*, that is,

$$
\delta\beta \equiv \beta - \beta\_{\text{eq}} = \beta - \beta(\mathbf{x}) \qquad \delta\zeta \equiv \zeta - \zeta\_{\text{eq}} = \zeta - \zeta(\mathbf{x})
$$

in a Fourier series, the only relevant components in the hydrodynamical limit for Integral (24) are those with very small frequency *ω* and wave vector **k**. At the same time, the vanishing of the flux can be achieved by enforcing periodicity of the perturbations in **x** − **x** . Taking these requirements into account, perturbations only include smallest wave four-vector *K*:

$$
\delta\beta\_V(\mathbf{x'}) \simeq A\_V \frac{1}{2i} (\mathbf{e}^{iK \cdot (\mathbf{x'} - \mathbf{x})} - \mathbf{e}^{-iK \cdot (\mathbf{x'} - \mathbf{x})}) \tag{25}
$$

with *Aν* being a real constant, the amplitude of the smallest wave four-vector Fourier component. The above form fulfils *δβ*(*x* ) = 0 as well as the request of vanishing flux, provided that *K<sup>i</sup>* = *π*/*Li*, with *Li* being the size of the compact domain in direction *i*. Hence, after the use of Equation (25), limit *K* → 0 is to be taken, which is equivalent to the limit of infinite volume. The gradient of Equation (25) (keep in mind that, in Equation (24), *∂μ* = *∂*/*∂xμ*) can then be written as:

$$
\partial\_{\mu}\beta\_{\nu} \simeq K\_{\mu}A\_{\nu}\frac{1}{2}(\mathbf{e}^{i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} + \mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})}) = \partial\_{\mu}\beta\_{\nu}(\mathbf{x})\text{Re }\mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} = \text{Re }\partial\_{\mu}\beta\_{\nu}(\mathbf{x})\mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} \tag{26}
$$

Plugging Equation (26) in Equation (24), in limit *K* → 0, one obtains:

$$
\begin{split}
\langle\langle\hat{O}(\mathbf{x})\rangle\rangle - \langle\hat{O}(\mathbf{x})\rangle\_{\rm LE} &\simeq \\
\langle\langle\mathbf{P}(\mathbf{x})\rangle\_{\rm LE} &\simeq \\
& -\partial\_{\boldsymbol{\mu}}\tilde{\zeta}(\mathbf{x}) \lim\_{K\to 0} \operatorname{Im} T \int\_{t\_{0}}^{t} \mathbf{d}^{4} \mathbf{x}' \int\_{t\_{0}}^{t'} \mathbf{d}\theta \,\langle\langle\hat{T}^{\mu\nu}(\theta,\mathbf{x}'),\hat{O}(\mathbf{x})\rangle\rangle\_{\beta(\mathbf{x})} \mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} \\
& -\partial\_{\boldsymbol{\mu}}\tilde{\zeta}(\mathbf{x}) \lim\_{K\to 0} \operatorname{Im} T \int\_{t\_{0}}^{t} \mathbf{d}^{4} \mathbf{x}' \int\_{t\_{0}}^{t'} \mathbf{d}\theta \,\langle\langle\hat{T}^{\mu}(\theta,\mathbf{x}'),\hat{O}(\mathbf{x})\rangle\rangle\_{\beta(\mathbf{x})} \mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} \end{split} \tag{27}
$$

As the macroscopic time scale *t* − *t*<sup>0</sup> and the microscopic time scale inherent in the correlators are so different, one can take limit *t*<sup>0</sup> → −∞. If functions

$$\int \mathbf{d}^3 \mathbf{x}' \, \langle [\hat{X}(\theta, \mathbf{x}'), \hat{O}(x)] \rangle\_{\mathcal{F}(x)}$$

with *X*" = *T*", "*j* remain finite for *θ* → −∞, then Equation (27), after integration by parts in *t* , can be turned into:

$$
\begin{split}
\langle\langle\hat{O}(\mathbf{x})\rangle\rangle - \langle\hat{O}(\mathbf{x})\rangle\_{\rm LE} &\simeq \qquad \left.\partial\_{\mu}\beta\_{\nu}(\mathbf{x})n^{a}\frac{\partial}{\partial\mathbf{K}^{a}}\Big|\_{\mathbf{x}\cdot\mathbf{K}=0} \lim\_{\mathbf{K}\rightarrow 0} \operatorname{Im} iT\int\_{-\infty}^{t} \mathbf{d}^{4}\mathbf{x}'\,\langle\langle\hat{O}(\mathbf{x}),\hat{T}^{\mu\nu}(\mathbf{x}')\rangle\rangle\_{\beta(\mathbf{x})}\mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} \\ & \qquad \quad - \partial\_{\mu}\tilde{\gamma}(\mathbf{x})n^{a}\frac{\partial}{\partial\mathbf{K}^{a}}\Big|\_{\mathbf{x}\cdot\mathbf{K}=0} \operatorname{Im} iT\int\_{-\infty}^{t} \mathbf{d}^{4}\mathbf{x}'\,\langle\langle\hat{O}(\mathbf{x}),\hat{T}^{\mu}(\mathbf{x}')\rangle\rangle\_{\beta(\mathbf{x})}\mathbf{e}^{-i\mathbf{K}\cdot(\mathbf{x}'-\mathbf{x})} \end{split} \tag{28}
$$

where *KT* is the projection of *K* orthogonal to *n*. This, as it becomes clear later, is the covariant form of the same formula obtained in Reference [1], with the (important) addition of the current term. In other words, it is the well-known formula expressing the transport coefficients as derivatives with

respect to the frequency of the retarded correlators of stress–energy components, the so-called Kubo formula. Defining:

$$\begin{split} \langle (\hat{X}, \hat{Y}) \equiv & \quad n^{a} \frac{\partial}{\partial K^{a}} \Big|\_{n:k=0} \lim\_{k\_{\Gamma} \to 0} \operatorname{Im} iT \int\_{-\infty}^{t} \mathrm{d}^{4} \mathbf{x}' \, \langle [\hat{X}(\mathbf{x}), \hat{Y}(\mathbf{x}')] \rangle\_{\hat{\mathbb{P}}(\mathbf{x})} \mathbf{e}^{-i\mathbf{K} \cdot (\mathbf{x}' - \mathbf{x})} \\ & = & \quad n^{a} \frac{\partial}{\partial K^{a}} \Big|\_{n:k=0} \lim\_{k\_{\Gamma} \to 0} \operatorname{Im} iT \int\_{-\infty}^{0} \mathrm{d}^{4} \mathbf{x}' \, \langle [\hat{X}(0), \hat{Y}(\mathbf{x}')] \rangle\_{\hat{\mathbb{P}}(\mathbf{x})} \mathbf{e}^{-i\mathbf{K} \cdot \mathbf{x}'}, \end{split} \tag{29}$$

which is bilinear in *X*" and *Y*", one can write the deviations of the stress–energy tensor from its LTE value as:

$$
\langle \hat{T}^{\mu\nu}(\mathbf{x}) \rangle - \langle \hat{T}^{\mu\nu}(\mathbf{x}) \rangle\_{\rm LE} \equiv \delta T^{\mu\nu}(\mathbf{x}) \simeq \left( \hat{T}^{\mu\nu}, \hat{T}^{\rho\sigma} \right) \partial\_{\rho} \beta\_{\sigma}(\mathbf{x}) - \left( \hat{T}^{\mu\nu}, \hat{f}^{\rho} \right) \partial\_{\rho} \zeta(\mathbf{x}).\tag{30}
$$

Similarly, the deviation of the current with respect to its value at LTE reads:

$$
\langle \hat{j}^{\mu}(\mathbf{x}) \rangle - \langle \hat{j}^{\mu}(\mathbf{x}) \rangle\_{\rm LE} \equiv \delta \hat{j}^{\mu}(\mathbf{x}) = \left( \hat{j}^{\mu}, \hat{T}^{\rho \sigma} \right) \partial\_{\rho} \beta\_{\sigma}(\mathbf{x}) - \left( \hat{j}^{\mu}, \hat{j}^{\rho} \right) \partial\_{\rho} \mathbb{Z}(\mathbf{x}).\tag{31}
$$

The next step is to decompose the correlators and the gradients of the relativistic fields into irreducible components under rotations, a procedure leading to the identification of familiar transport coefficients: shear and bulk viscosities, thermal conductivities, etc. We are not going to show how this is accomplished, but we would just like to point out, for the purpose of identifying the transport coefficients, that the gradients of *β* can be turned into gradients of velocity field *u* by using Equation (5). Having defined

and

$$
\Delta\_{\mu\nu} = \mathcal{g}\_{\mu\nu} - \mathcal{u}\_{\mu}\mathcal{u}\_{\nu}
$$

$$D = \mu \cdot \partial \qquad \qquad \nabla\_T^{\mu} = \partial^{\mu} - \mu^{\mu} D\_{\nu}$$

the transverse gradients of velocity field <sup>∇</sup>*<sup>μ</sup> <sup>T</sup>u<sup>ν</sup>* can be written as follows:

$$\begin{split} \nabla\_{T\mu}\mu^{\nu} &= \nabla\_{T\mu}\frac{\beta^{\nu}}{\sqrt{\beta^{2}}} = -\frac{1}{2}\beta^{\nu}(\beta^{2})^{-3/2}\nabla\_{T\mu}\beta^{2} + \frac{1}{\sqrt{\beta^{2}}}\nabla\_{T\mu}\beta^{\nu} \\ &= \frac{1}{\sqrt{\beta^{2}}} \left( -\frac{\beta^{\nu}\beta^{\rho}}{\beta^{2}}\nabla\_{T\mu}\beta\_{\rho} + \nabla\_{T\mu}\beta^{\nu} \right) = \frac{1}{\sqrt{\beta^{2}}}\Delta^{\rho\nu}\nabla\_{T\mu}\beta\_{\rho\nu} \end{split} \tag{32}$$

where we have used Relation (5). Thereby, the Navier–Stokes shear term can be fully expressed in terms of inverse temperature four-vector *β* and its gradients. The same transformation can be proven for the other terms [8].

#### **6. Outlook**

The nonequilibrium statistical operator method introduced by D. Zubarev more than 50 years ago has been a very important achievement in statistical physics and has not received deserved attention. It can be used in all physical problems where a local thermodynamic equilibrium is reached, and it can be quite straightforwardly extended to relativistic statistical mechanics. In this work, we presented an amendment of its original formulation that reproduces known results and makes its application easier to relativistic hydrodynamics problems. Since it is a fully fledged quantum framework, this approach is especially suitable for the calculation of quantum effects. Among various applications, recent evidence of nonvanishing polarization in quark–gluon plasma [13] makes it the ideal tool to deal with this newly found phenomenon.

**Author Contributions:** Conceptualization, F.B.; Methodology, F.B., M.B., E.G.; Investigation, F.B., M.B. and E.G.; Writing–original draft preparation, F.B.; Writing–review & editing, F.B. and M.B.

**Funding:** This research was funded by INFN under the grant SIM.

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

#### **Appendix A. Supplementary Notes on the Derivation of the Kubo Formula**

Working out Equation (20) requires Equations (8) and (9). By also using Approximations (21) and (23), Equation (20) turns into:

$$
\begin{split}
\langle\langle\hat{\mathcal{O}}(\mathbf{x})\rangle-\langle\hat{\mathcal{O}}(\mathbf{x})\rangle\_{\mathrm{LE}} &\rightharpoonup -\int\_{l\_{0}}^{l} \mathrm{d}\mathbf{f}' \int d\mathbf{\hat{x}}' \,\langle\hat{\mathcal{O}}(\mathbf{x})\rangle\_{\mathrm{LE}} \left(\langle\hat{\mathcal{T}}^{\mu\nu}(t',\mathbf{x}')\rangle\_{\mathrm{LE}}\mathfrak{a}\_{\mu}\mathfrak{b}\_{\nu} - \langle\hat{\mathcal{P}}^{\mu}(t',\mathbf{x}')\rangle\_{\mathrm{LE}}\mathfrak{a}\_{\mu}\mathfrak{b}\_{\nu}\right) \\ &+\int\_{0}^{1} \mathrm{d}\mathbf{z} \int\_{l\_{0}}^{l} \mathrm{d}\mathbf{f}' \int d\mathbf{\hat{x}}' \int \langle\langle\hat{\mathcal{O}}(\mathbf{x})\mathbf{e}^{-\mathbf{z}\langle\boldsymbol{\beta}(\mathbf{x})\cdot\hat{\mathcal{P}}-\zeta(\mathbf{x})\langle\hat{\mathcal{Q}}\rangle}\hat{\mathcal{T}}^{\mu\nu}(t',\mathbf{x}')\mathfrak{e}^{\mathbf{z}\langle\boldsymbol{\beta}(\mathbf{x})\cdot\hat{\mathcal{P}}-\zeta(\mathbf{x})\langle\hat{\mathcal{Q}}\rangle}\rangle\_{\mathrm{LE}}\mathfrak{a}\_{\mu}\mathfrak{b}\_{\nu}
\\ &-\langle\hat{\mathcal{O}}(\mathbf{x})\mathbf{e}^{-\mathbf{z}\langle\boldsymbol{\beta}(\mathbf{x})\cdot\hat{\mathcal{P}}-\zeta(\mathbf{x})\langle\hat{\mathcal{Q}}\rangle}\hat{\mathcal{P}}^{\mu}(t',\mathbf{x}')\hat{\mathcal{P}}^{\nu}(t',\mathbf{x}')\hat{\mathcal{P}}^{\nu}(t',\mathbf{x}')\hat{\mathcal{Q}}^{\nu}\rangle\_{\mathrm{LE}}\mathfrak{b}\_{\mu}\mathfrak{b}\_{\nu}
\end{split}
$$

With [*Q*", *T*"(*x*)] = 0 and [*Q*", "*j*(*x*)] = 0, one can also write

$$\mathbf{e}^{-z(\boldsymbol{\beta}(\mathbf{x})\cdot\boldsymbol{P}-\boldsymbol{\zeta}(\mathbf{x})\bar{\boldsymbol{Q}})} \hat{\mathbf{X}}(t',\mathbf{x'}) \,\mathbf{e}^{z(\boldsymbol{\beta}(\mathbf{x})\cdot\boldsymbol{P}-\boldsymbol{\zeta}(\mathbf{x})\bar{\boldsymbol{Q}})} = \hat{\mathbf{X}}(t'+iz\sqrt{\boldsymbol{\beta}^2},\mathbf{x'})$$

with *<sup>X</sup>*" <sup>=</sup> *<sup>T</sup>*", "*j*, where, in the last expression, we tacitly assumed that *<sup>n</sup>* <sup>=</sup> *<sup>β</sup>*ˆ, i.e., that the local equilibrium hypersurface coincides, locally around *x*, with the hypersurface normal to *β* [8]. Hence, the last term in the right-hand side of Equation (A1) can be rewritten as:

$$\int\_0^1 \mathrm{d}z \int\_{t\_0}^t \mathrm{d}t' \int \mathrm{d}^3 \mathbf{x}' \left( \langle \hat{O}(\mathbf{x}) \hat{T}^{\mu \nu} (t' + iz \sqrt{\beta^2}, \mathbf{x}') \rangle\_{\mathrm{LE}} \partial\_{\mu} \beta\_{\nu} - \langle \hat{O}(\mathbf{x}) \hat{f}^{\mu} (t' + iz \sqrt{\beta^2}, \mathbf{x}') \rangle\_{\mathrm{LE}} \partial\_{\mu} \zeta \right) d\zeta$$

provided that *n* = *β*ˆ, that is, if the local equilibrium hypersurface locally coincides with the hypersurface normal to *β* [8]. Operator *X*" = *T*", "*j*, can be rewritten as

$$\begin{split} \hat{X}(t' + iz\sqrt{\beta^2}, \mathbf{x'}) &= \quad \hat{X}(t\_0 + iz\sqrt{\beta^2}, \mathbf{x'}) + \int\_{t\_0}^{t'} \mathbf{d}\theta \, \frac{\partial}{\partial \theta} \hat{X}(\theta + iz\sqrt{\beta^2}, \mathbf{x'}) \\ &= \quad \hat{X}(t\_0 + iz\sqrt{\beta^2}, \mathbf{x'}) + \int\_{t\_0}^{t'} \mathbf{d}\theta \, \frac{1}{i\sqrt{\beta^2}} \frac{\partial}{\partial z} \hat{X}(\theta + iz\sqrt{\beta^2}, \mathbf{x'}) \end{split} \tag{A2}$$

and integrating in *z*:

$$\begin{split} \int\_{0}^{1} \mathrm{d}z \langle \hat{\mathcal{O}}(\mathbf{x}) \hat{\mathcal{X}}(t' + iz \sqrt{\beta^2}, \mathbf{x'}) \rangle\_{\mathrm{LE}} &= \int\_{0}^{1} \mathrm{d}z \, \langle \hat{\mathcal{O}}(\mathbf{x}) \hat{\mathcal{X}}(t\_{0} + iz \sqrt{\beta^2}, \mathbf{x'}) \rangle\_{\mathrm{LE}} \\ + \int\_{t\_{0}}^{t'} \mathrm{d}\theta \, \frac{1}{i\sqrt{\beta^2}} \left( \langle \hat{\mathcal{O}}(\mathbf{x}) \hat{\mathcal{X}}(\theta + i\sqrt{\beta^2}, \mathbf{x'}) \rangle\_{\mathrm{LE}} - \langle \hat{\mathcal{O}}(\mathbf{x}) \hat{\mathcal{X}}(\theta, \mathbf{x'}) \rangle\_{\mathrm{LE}} \right). \end{split} \tag{A3}$$

Now, we use the same approximation of *A*" as in Equation (21), and LTE mean values are calculated at equilibrium, with density Operator (22). Thus,

$$\begin{split} & \langle \langle \hat{O}(\mathbf{x}) \hat{\mathcal{X}} (\theta + i \sqrt{\beta^2}, \mathbf{x'}) \rangle\_{\mathrm{LE}} - \langle \hat{O}(\mathbf{x}) \hat{\mathcal{X}} (\theta, \mathbf{x'}) \rangle\_{\mathrm{LE}} \simeq \langle \hat{O}(\mathbf{x}) \hat{\mathcal{X}} (\theta + i \sqrt{\beta^2}, \mathbf{x'}) \rangle\_{\beta(\mathbf{x})} - \langle \hat{O}(\mathbf{x}) \hat{\mathcal{X}} (\theta, \mathbf{x'}) \rangle\_{\beta(\mathbf{x})} \\ &= \langle \langle \hat{O}(\mathbf{x}) \mathbf{e}^{-\beta(\mathbf{x}) \cdot \hat{\mathcal{P}} + \zeta(\mathbf{x}) \hat{\mathcal{Q}}} \hat{\mathcal{R}} (\theta, \mathbf{x'}) \mathbf{e}^{\beta(\mathbf{x}) \cdot \hat{\mathcal{P}} - \zeta(\mathbf{x}) \hat{\mathcal{Q}}} \rangle\_{\beta(\mathbf{x})} - \langle \hat{O}(\mathbf{x}) \hat{\mathcal{X}} (\theta, \mathbf{x'}) \rangle\_{\beta(\mathbf{x})} \\ &= \langle \hat{\mathcal{X}} (\theta, \mathbf{x'}) \hat{\mathcal{O}}(\mathbf{x}) \rangle\_{\beta(\mathbf{x})} - \langle \hat{O}(\mathbf{x}) \hat{\mathcal{X}} (\theta, \mathbf{x'}) \rangle\_{\beta(\mathbf{x})} = \langle [\hat{\mathcal{X}} (\theta, \mathbf{x'}), \hat{O}(\mathbf{x})] \rangle\_{\beta(\mathbf{x})} \end{split}$$

Substitution into Equation (A3) yields:

$$\int\_{0}^{1} \mathrm{d}z \langle \hat{O}(\mathbf{x}) \hat{X}(t + iz\sqrt{\beta^{2}}, \mathbf{x}') \rangle\_{\mathrm{LE}} \simeq \int\_{0}^{1} \mathrm{d}z \, \langle \hat{O}(\mathbf{x}) \hat{X}(t\_{0} + iz\sqrt{\beta^{2}}, \mathbf{x}') \rangle\_{\hat{\mathbb{P}}(\mathbf{x})} + \frac{1}{i\sqrt{\beta^{2}}} \int\_{t\_{0}}^{t'} \mathrm{d}\theta \, \langle [\hat{X}(\theta, \mathbf{x}'), \hat{O}(\mathbf{x})] \rangle\_{\hat{\mathbb{P}}(\mathbf{x})}.\tag{A4}$$

Using this result for *X*" = *T*", "*j* allows to turn Equation (A1) into:

*O*"(*x*)−*O*"(*x*)LE *<sup>t</sup> <sup>t</sup>*<sup>0</sup> d*t* d3x <sup>1</sup> <sup>0</sup> d*z O*"(*x*)*T*"*μν*(*t*<sup>0</sup> <sup>+</sup> *izβ*2, **<sup>x</sup>** )*β*(*x*) − *O*"(*x*)*β*(*x*)*T*"*μν*(*<sup>t</sup>* , **x** )*β*(*x*) *∂μβν O*"(*x*)"*j <sup>μ</sup>*(*t*<sup>0</sup> + *izβ*2, **<sup>x</sup>** )*β*(*x*) − *O*"(*x*)*β*(*x*)"*j <sup>μ</sup>*(*t* , **x** )*β*(*x*) *∂μζ* + *iT <sup>t</sup> <sup>t</sup>*<sup>0</sup> d*t <sup>t</sup> <sup>t</sup>*<sup>0</sup> d*θ* d3x [*O*"(*x*), *<sup>T</sup>*"*μν*(*θ*, **<sup>x</sup>** )]*β*(*x*)*∂μβν*(*x* ) − [*O*"(*x*), "*j <sup>μ</sup>*(*θ*, **x** )]*β*(*x*)*∂μζ*(*x* ) . (A5)

In the paper by A. Hosoya et al. [1], in limit *t*<sup>0</sup> → −∞, the first of the two integral terms is shown to be vanishing, based on the idea that lim*t*0→−<sup>∞</sup> *<sup>X</sup>*"(*t*<sup>0</sup> <sup>+</sup> *izβ*2, **<sup>x</sup>** ) lim*t*0→−<sup>∞</sup> *X*"(*t*0, **x** ) and that correlation between operator *O*" at time *t* and *X*" at an infinitely remote past is 0, that is,

$$\lim\_{t\_0 \to -\infty} \langle \hat{O}(\mathbf{x}) \hat{T}^{\mu\nu}(t\_0, \mathbf{x}') \rangle\_{\beta(\mathbf{x})} \simeq \lim\_{t\_0 \to -\infty} \langle \hat{O}(\mathbf{x}) \rangle\_{\beta(\mathbf{x})} \langle \hat{T}^{\mu\nu}(t\_0, \mathbf{x}') \rangle\_{\beta(\mathbf{x})} = \langle \hat{O}(\mathbf{x}) \rangle\_{\beta(\mathbf{x})} \langle \hat{T}^{\mu\nu}(t', \mathbf{x}') \rangle\_{\beta(\mathbf{x})}$$

where, in the last equality, we took advantage of the fact that the mean value of any operator is constant at equilibrium. Therefore, the first integral on the right-hand side of Equation (A5) vanishes, and we are left only with the second integration, that is, Equation (24).

Let us now define the *λ*-dependent partition function:

$$Z\_{\rm LE}(\lambda) = \text{tr}\left(\exp\left[-\lambda \text{d}\Sigma \ n\_{\mu} \left(\widehat{T}^{\mu\nu}\beta\_{\nu} - \widehat{j}^{\mu}\zeta\right)\right]\right)$$

If we take the derivative with regard to *λ*, we have

$$\frac{\partial}{\partial \lambda} \log Z\_{\rm LE}(\lambda) = \int \mathrm{d}\Sigma \, n\_{\mu} \left( \langle \hat{T}^{\mu \nu} \rangle\_{\rm LE}(\lambda) \beta\_{\nu} - \zeta \langle \hat{j}^{\mu} \rangle\_{\rm LE}(\lambda) \right),$$

whence

$$\log Z\_{\rm LE}(1) - \log Z\_{\rm LE}(\lambda\_0) = \int\_{\lambda\_0}^1 \mathrm{d}\lambda \int \mathrm{d}\Sigma \ n\_{\mu} \left( \langle \hat{T}^{\mu\nu} \rangle\_{\rm LE}(\lambda) \beta\_{\nu} - \zeta \langle \hat{j}^{\mu} \rangle\_{\rm LE}(\lambda) \right)$$

If we find *λ*0, such that log *Z*LE(*λ*0) = 0, which happens under some reasonable assumptions *λ*<sup>0</sup> = +∞, and inverting the order of integrations:

$$\log Z\_{\rm LE} = \int \mathrm{d}\Sigma \, n\_{\mu} \int\_{\lambda\_{0}}^{1} \mathrm{d}\lambda \left( \langle \hat{\mathcal{I}}^{\mu\nu} \rangle\_{\rm LE}(\lambda) \beta\_{\nu} - \zeta \langle \hat{\mathcal{J}}^{\mu} \rangle\_{\rm LE}(\lambda) \right) \equiv \int \mathrm{d}\Sigma\_{\mu} \phi^{\mu} \,, \tag{A6}$$

which shows that log *Z*LE is extensive, i.e., it can be written as an integral over the 3D hypersurface of a current, the thermodynamic potential current. At the same time, the above equation provides a formula to calculate *φ<sup>μ</sup>* as an integral in *λ*.

#### **References**


c 2019 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/).
