*Article* **Calculation of Acceleration Effects Using the Zubarev Density Operator**

**Georgy Prokhorov 1,\* , Oleg Teryaev 1,2,3 and Valentin Zakharov 2,4,5**


Received: 16 November 2019; Accepted: 23 December 2019; Published: 3 January 2020

**Abstract:** The relativistic form of the Zubarev density operator can be used to study quantum effects associated with acceleration of the medium. In particular, it was recently shown that the calculation of perturbative corrections in acceleration based on the Zubarev density operator makes it possible to show the existence of the Unruh effect. In this paper, we present the details of the calculation of quantum correlators arising in the fourth order of the perturbation theory needed to demonstrate the Unruh effect. Expressions for the quantum corrections for massive fermions are also obtained.

**Keywords:** Zubarev operator; Unruh effect; acceleration

#### **1. Introduction**

There are wonderful quantum-field effects associated with non-uniform motion of the medium. A well-known example of such an effect is the Unruh effect, according to which an accelerated observer perceives the Minkowski vacuum as a medium filled with particles with a temperature depending on the acceleration [1]. This temperature is called the Unruh temperature, and it is equal to

$$T\_{ll} = \frac{\hbar|a|}{2\pi c\mathbf{k}}\,. \tag{1}$$

The Unruh effect is similar to the Hawking effect, since it is also associated with the appearance of the event horizon in the accelerated system. This effect continues to be the focus of theorists [2–5]. The possibility of experimental observation of the Unruh effect needs the generation of ultrahigh acceleration in a system, which is relevant, in particular, for particle collisions [6,7] and systems with two-level atoms in quantum optics [8–10].

There is a universal fundamental statistical approach to describing the equilibrium thermodynamics of quantized fields. This approach is based on the relativistic form of the Zubarev density operator [11,12]. It has recently been shown that this approach allows to study in a regular way the effects of rotation and acceleration in a medium of relativistic particles [13–15].

Using the Zubarev operator method, various effects associated with the motion of the medium are shown. In particular, the well-known chiral vortical effect is shown and corrections to this effect are calculated [13,14,16]. Since the chiral vortical effect is associated with the axial electromagnetic anomaly [17–19], as well as with the gravitational anomaly [20,21], it turns out that the approach with the Zubarev operator carries information about the most fundamental properties of matter.

A remarkable observation made recently is that the Unruh effect can also be obtained from the Zubarev density operator [22,23]. Relativistic quantum statistical mechanics considers a continuous medium filled with particles described by quantized fields. This medium in equilibrium is characterized by a number of thermodynamic parameters, such as temperature, energy density, pressure, and others. Non-trivial aspect is connected with the need for normalization of the thermodynamic quantities of the system to a specific vacuum, as a rule, the Minkowski vacuum. With such a statement of the problem, a direct consequence of the Unruh effect from the point of view of quantum statistical mechanics is the vanishing of the observables, in particular, the energy-momentum tensor, at a proper temperature equal to the Unruh temperature [24,25]. This is exactly what was found in [23].

This means that, in the Zubarev approach, nontrivial gravitational effects, associated with the occurrence of an event horizon, and the changes in vacuum properties depending on the reference system, are reproduced. This observation seems even more surprising because the corresponding calculation was carried out in ordinary flat Minkowski space-time, that is, by observing an accelerated medium from an inertial frame. Nevertheless, nontrivial physics associated with Unruh effect is reproduced.

Moreover, as discussed in [26], the Zubarev density operator exactly reproduces quantum corrections that were derived in space of a cosmic string, characterized by a conical singularity [25,27]. The existence of such exact duality means that the Zubarev operator of the accelerated medium leads to emergent conical geometry.

To justify the Unruh effect in [23], it was necessary to calculate a five-point correlator with boost operators and energy-momentum tensor. This calculation in [23] was made for the massless Dirac field. The method we used was developed in a series of works [13–15], where the perturbation theory with the boost operator was developed and corrections up to the second order were calculated. It is well known [13] that the boost operator does not commute with the Hamiltonian of the system. Because of this, with each subsequent order of the perturbation theory, the complexity of calculation of the corresponding quantum correlators increases. The fourth order found in [23] is currently a record one. In the present paper, we describe a never before given scheme for calculating higher orders of the perturbation theory with the boost operator and also derive expressions for the fourth-order corrections to the energy-momentum tensor at nonzero mass.

To date, the Unruh effect has been considered from various points of view. In particular, in the framework of quantum optics [8–10], the Unruh effect manifests itself in the thermal distribution with the Unruh temperature, in the probability of absorption and emission of gamma quanta by accelerated two-level atoms. It is necessary to consider the interaction of atoms with an electromagnetic field in the framework of perturbation theory with respect to the coupling constant, while acceleration effects can be taken into account in a nonperturbative way through Rindler coordinates.

Despite the difference in approaches, a parallel can be established between our consideration and the usual approach to the Unruh effect, as well as quantum optics. In particular, in the statistical approach we also obtained a term in the energy density (which is the last term in Equation (3.1) in [23]), which corresponds to the Bose distribution of gamma quanta at the Unruh temperature.

The paper has the following structure. Section 2 introduces the basic concepts of the method of Zubarev density operator. An algorithm of constructing a perturbation theory in acceleration is also discussed. In the Section 3 we describe in details the calculation of the corrections of the fourth-order in acceleration to the energy-momentum tensor. The interpretation associated with the Unruh effect is given in Section 4. In Section 5 the conclusions are given. In the Appendix A the formulas for the coefficients at finite mass are presented.

The system of units ¯*h* = *c* = *k* = 1 is used.

#### **2. Perturbation Theory in Acceleration Based on the Zubarev Density Operator**

In this section, we introduce the basic concepts related to the density operator and describe how the acceleration perturbation theory can be constructed. In general, in this section we follow the paper [13]. In [11,12], a relativistic form of the Zubarev density operator was obtained for a medium in a state of local thermodynamic equilibrium

$$\hat{\rho} = \frac{1}{Z} \exp\left\{-\int\_{\Sigma} d\Sigma\_{\mu} [\hat{T}^{\mu\nu}(\mathbf{x})\beta\_{\nu}(\mathbf{x}) - \xi(\mathbf{x})\hat{\jmath}^{\mu}(\mathbf{x})] \right\},\tag{2}$$

where integration over the three-dimensional hypersurface <sup>Σ</sup> is performed. Here, *βμ* <sup>=</sup> *<sup>u</sup><sup>μ</sup> <sup>T</sup>* is the inverse temperature 4-vector, *T* is the proper temperature, *ξ* = *<sup>u</sup> <sup>T</sup>* is the ratio of the chemical potential in the co-moving frame to temperature, *T*ˆ *μν* and ˆ*j <sup>μ</sup>* are the energy-momentum tensor and current operators. The conditions of global thermodynamic equilibrium for a medium with rotation and acceleration, that is, conditions under which the density operator (2) becomes independent on the choice of the hypersurface Σ, over which the integration occurs, thus acquiring the properties of a density operator in a state of global thermodynamic equilibrium, have the form [13,15,28,29]

$$\begin{aligned} \beta\_{\mu} &= b\_{\mu} + \mathcal{O}\_{\mu\nu} \mathbf{x}\_{\nu}, \quad b\_{\mu} = \text{const}, \quad \mathcal{O}\_{\mu\nu} = \text{const}, \\ \mathcal{O}\_{\mu\nu} &= -\frac{1}{2} (\partial\_{\mu} \beta\_{\nu} - \partial\_{\nu} \beta\_{\mu}), \quad \xi^{\tag{3}} = \text{const}, \end{aligned} \tag{3}$$

where *μν* is the thermal vorticity tensor. In the general case, integration over the hypersurface is to be done and the quantum statistical theory should be projected to this hypersurface [30–32]. So under the condition (3), the density operator (2) becomes the global equilibrium density operator [13,15,22]

$$\boldsymbol{\hat{\rho}} = \frac{1}{Z} \exp\left\{-\boldsymbol{\beta}\_{\mu}(\mathbf{x})\boldsymbol{\hat{P}}^{\mu} + \frac{1}{2}\boldsymbol{\sigma}\_{\mu\nu}\boldsymbol{\hat{I}}\_{\mu}^{\mu\nu} + \boldsymbol{\xi}^{\mu}\boldsymbol{\hat{Q}}\right\},\tag{4}$$

where *P*ˆ is the 4-momentum operator, *Q*ˆ is the charge operator, and ˆ*Jx* are the generators of the Lorentz transformations shifted to the point *x*

$$\mathcal{J}\_{\mathbf{x}}^{\mu\nu} = \int d\Sigma\_{\lambda} \left[ (y^{\mu} - \mathbf{x}^{\mu}) \mathcal{I}^{\lambda\nu}(y) - (y^{\nu} - \mathbf{x}^{\nu}) \mathcal{I}^{\lambda\mu}(y) \right]. \tag{5}$$

The technique of calculating the mean values of physical quantities based on (4) was developed in [13,15], in which second-order corrections in the thermal vorticity tensor were calculated to various thermodynamic quantities for scalar and Dirac fields.

Note that the condition (3) also lead to a system of kinematic equations of motion, solving which, we can construct trajectories of motion. Particular cases of this solution are the rotation of the medium as a solid, as well as uniformly accelerated motion.

Following [13], we introduce the thermal acceleration vector *αμ* and the thermal vorticity pseudo-vector *wμ*

$$
\alpha\_{\mu} = \mathcal{O}\_{\mu\nu} u^{\nu}, \quad w\_{\mu} = -\frac{1}{2} \epsilon\_{\mu\nu\alpha\beta} u^{\nu} \mathcal{O}^{a\beta}. \tag{6}
$$

Drawing a parallel with the electrodynamics, *αμ* and *w<sup>μ</sup>* can be called the "electrical" and "magnetic" components of the tensor . The tensor *μν* can be decomposed into these components as follows

$$
\omega\_{\mu\nu} = \epsilon\_{\mu\nu\alpha\beta} \varpi^{\mu} \iota^{\beta} + \omega\_{\mu} \iota\_{\nu} - \alpha\_{\nu} \iota\_{\mu} \,. \tag{7}
$$

*Particles* **2020**, *3*

The meaning of the vectors *αμ* and *w<sup>μ</sup>* becomes clear when considering the case of global thermodynamic equilibrium, in which they are proportional to the usual kinematic 4-acceleration *aμ* and vorticity *ωμ*

$$
\mu^{\mu} = \mathcal{\alpha}^{\mu}{}\_{\nu} u^{\nu} = \mu^{\nu} \partial\_{\nu} \beta^{\mu} = \frac{1}{T} \mu^{\nu} \partial\_{\nu} u^{\mu} = \frac{a^{\mu}}{T} \tag{8}
$$

and for thermal vorticity, we get

$$w\_{\mu} = -\frac{1}{2}\epsilon\_{\mu\nu\alpha\beta}u^{\nu}\omega^{a\beta} = -\frac{1}{2}\epsilon\_{\mu\nu\alpha\beta}u^{\nu}\partial^{\beta}\beta^{a} = \frac{1}{2T}\epsilon\_{\mu\nu\alpha\beta}u^{\nu}\partial^{a}u^{\beta} = \frac{\omega\_{\mu}}{T}.\tag{9}$$

In the rest frame, *a<sup>μ</sup>* and *ω<sup>μ</sup>* are expressed in terms of three-dimensional vectors

$$a^{\mu} = \begin{pmatrix} 0, \mathbf{a} \end{pmatrix}, \quad \omega^{\mu} = \begin{pmatrix} 0, \mathbf{w} \end{pmatrix}, \tag{10}$$

where **a** and **w** are three-dimensional acceleration and angular velocity.

The density operator (4) allows one to find corrections related to thermal vorticity in the framework of perturbation theory. To do this, it is necessary to expand (4) in a series taking into account the fact that we are constructing a perturbation theory with non-commuting operators. According to [13] we have

$$
\langle \mathcal{O}(\mathbf{x}) \rangle = \langle \mathcal{O}(0) \rangle\_{\beta(\mathbf{x})} + \sum\_{N=1}^{\infty} \frac{\mathcal{O}^N}{2^N N! |\beta|^N} \int\_0^{|\beta|} d\tau\_1 d\tau\_2 \dots d\tau\_N \langle T\_\tau f\_{-i\tau\_1 \mathbf{u}} \dots f\_{-i\tau\_N \mathbf{u}} \mathcal{O}(0) \rangle\_{\beta(\mathbf{x}), \mathbf{c}'} \tag{11}
$$

where it is assumed that each of the thermal vorticity tensors is contracted with the tensor ˆ*J* so that *μν* ˆ*Jμν*. Equation (11) includes only connected correlators, all disconnected correlators are reduced due to the contribution of the denominator 1/*Z* to (4). This fact is reflected in the subscript *c*; the subscript *β*(*x*) means that the mean values are taken at = 0, that is, averaging is performed over a grand canonical distribution. The *<sup>T</sup><sup>τ</sup>* operator orders in imaginary time *<sup>τ</sup>*, and <sup>|</sup>*β*<sup>|</sup> <sup>=</sup> *βμβμ* <sup>=</sup> <sup>1</sup> *T* .

It is convenient to introduce the boost operator *K*ˆ and the angular momentum operator ˆ*J*

$$
\hat{J}^{\mu\nu} = u^{\mu}\hat{\mathbb{K}}^{\nu} - u^{\nu}\hat{\mathbb{K}}^{\mu} - \epsilon^{\mu\nu\rho\sigma}u\_{\rho}\hat{\mathbb{J}}\_{\sigma}.\tag{12}
$$

From (7) and (12), it follows that scalar products with vorticity tensor in (4) and (11) decompose into terms with boost and angular momentum

$$
\omega\_{\mu\nu}\hat{f}\_{\mathbf{x}}^{\mu\nu} = -2\alpha\_{\mu}\hat{\mathcal{K}}\_{\mathbf{x}}^{\mu} - 2w\_{\mu}\hat{f}\_{\mathbf{x}}^{\mu} \,. \tag{13}
$$

Further, we will consider uniformly accelerated media without vorticity and chemical potential; therefore (4), transforms to the density operator of the form

$$\beta = \frac{1}{Z} \exp\left\{-\beta\_{\mu} \mathcal{P}^{\mu} - a\_{\mu} \mathcal{K}\_{x}^{\mu}\right\},\tag{14}$$

and the perturbation theory in (11) takes the form of the series in acceleration

$$
\langle \hat{O}(\mathbf{x}) \rangle = \langle \hat{O}(0) \rangle\_{\hat{\mathcal{P}}(\mathbf{x})} + \sum\_{N=1}^{\infty} \frac{(-1)^N a^N}{N!} \int\_0^{|\mathcal{S}|} d\tau\_1 d\tau\_2 ... d\tau\_N \langle T\_\tau \hat{\mathbb{K}}\_{-i\tau\_1 \mathcal{U}} ... \hat{\mathbb{K}}\_{-i\tau\_N \mathcal{U}} \hat{O}(0) \rangle\_{\hat{\mathcal{P}}(\mathbf{x}), \mathbf{c}},\tag{15}
$$

#### **3. Calculation of Fourth-Order Coefficients in Acceleration**

The second-order coefficients in acceleration in the energy-momentum tensor of the Dirac field were calculated in [13,14]. In this section, we present the details of calculation of the fourth-order coefficient.

The operator form of the energy-momentum tensor of the mass-less Dirac fields is well known. We will use the symmetrized Belinfante energy-momentum tensor

$$\hat{T}^{\mu\nu} = \frac{i}{4} \left( \overline{\Psi}\gamma^{\mu}\partial^{\nu}\Psi - \partial^{\nu}\overline{\Psi}\gamma^{\mu}\Psi + \overline{\Psi}\gamma^{\nu}\partial^{\mu}\Psi - \partial^{\mu}\overline{\Psi}\gamma^{\nu}\Psi \right) \,. \tag{16}$$

As follows from (15), the calculation of the necessary correlators is performed in imaginary time—a time shift is made along the imaginary axis. Thus, it is necessary to pass to the Euclidean formalism in imaginary time. The Euclidean version of the energy-momentum tensor (16) has the form

$$\mathcal{T}^{\mu\nu} = \frac{\mathbf{r}^{\delta\_{0\mu} + \delta\_{0\nu}}}{4} \left( \overline{\Psi} \gamma^{\mu} \partial^{\nu} \Psi - \partial^{\nu} \overline{\Psi} \gamma^{\mu} \Psi + \overline{\Psi} \gamma^{\nu} \partial^{\mu} \Psi - \partial^{\mu} \overline{\Psi} \gamma^{\nu} \Psi \right) \,, \tag{17}$$

where *γ*˜ are the Euclidean Dirac matrices

$$
\hat{\gamma}\_{\mu} = i^{1-\delta\_{0\mu}} \gamma\_{\mu}, \quad \hat{\gamma}^{\mu} = i^{1-\delta\_{0\mu}} \gamma^{\mu}, \quad \{\hat{\gamma}\_{\mu}\hat{\gamma}\_{\nu}\} = 2\delta\_{\mu\nu}.\tag{18}
$$

and derivatives are also taken in Euclidean space-time, so that

$$
\tilde{\partial}\_{\mu} = (-i)^{\delta\_{0\mu}} \partial\_{\mu} \,. \tag{19}
$$

However, we will omit the tilde sign for derivatives. Consider the mean value of the energy-momentum tensor in the fourth order of the perturbation theory in acceleration using (15)

$$
\begin{split}
\langle\langle\mathcal{T}^{\mu\nu}(\mathbf{x})\rangle\rangle &=& \langle\langle\mathcal{T}^{\mu\nu}(0)\rangle\_{\beta(\mathbf{x})} + \frac{a\_{\rho}a\_{\sigma}}{2} \int\_{0}^{|\beta|} d\tau\_{1}d\tau\_{2} \langle\mathcal{T}\_{\mathbf{r}}\mathcal{K}^{\rho}\_{-i\tau\_{1}u}\mathcal{K}^{\sigma}\_{-i\tau\_{2}u}\mathcal{T}^{\mu\nu}(0)\rangle\_{\beta(\mathbf{x}),\varepsilon} \\ &+ \frac{8a\_{\rho}a\_{\sigma}a\_{\gamma}a\_{\eta}}{4!} \int\_{0}^{|\beta|} d\tau\_{1}d\tau\_{2}d\tau\_{3}d\tau\_{4} \langle\mathcal{T}\_{\mathbf{r}}\dot{\mathcal{K}}^{\rho}\_{-i\tau\_{1}u}\dot{\mathcal{K}}^{\sigma}\_{-i\tau\_{2}u}\dot{\mathcal{K}}^{\gamma}\_{-i\tau\_{3}u}\dot{\mathcal{K}}^{\mu\nu}(0)\rangle\_{\beta(\mathbf{x}),\varepsilon} + \mathcal{O}(a^{6}).
\end{split}
\tag{20}
$$

Symmetry and parity considerations fix the form of the energy-momentum tensor in the fourth order of perturbation theory

$$\begin{aligned} \langle \hat{T}^{\mu\nu} \rangle &= \langle \rho\_0 - A\_1 T^2 a^2 + A\_2 a^4 \rangle u^\mu u^\nu - (p\_0 - A\_3 T^2 a^2 + A\_4 a^4) \Delta^{\mu\nu} \\ &+ (A\_5 T^2 - A\_6 a^2) a^\mu a^\nu + \mathcal{O}(a^6) \qquad \Delta^{\mu\nu} = g^{\mu\nu} - u^\mu u^\nu, \end{aligned} \tag{21}$$

where *a*<sup>2</sup> = *aμaμ*. As already mentioned, 2-order coefficients were calculated earlier in [13,14]. Our goal is to calculate coefficients of the 4th order *A*2, *A*4, *A*6. Comparing (20) with (21), we obtain

$$\begin{split} A\_2 a^4 u^\mu u^\nu - A\_4 a^4 \Delta^{\mu \nu} - A\_6 a^2 a^\mu a^\nu &= \frac{a\_\rho a\_\sigma a\_\tau a\_\eta}{4!} \int\_0^{|\mathcal{S}|} d\tau\_1 d\tau\_2 d\tau\_3 d\tau\_4 \\ &\times \langle T\_\tau \mathcal{K}^\rho\_{-i\tau\_1 \mu} \mathcal{K}^\sigma\_{-i\tau\_2 \mu} \mathcal{K}^\gamma\_{-i\tau\_3 \mu} \mathcal{K}^\eta\_{-i\tau\_4 \mu} \mathcal{I}^{\mu \nu}(0) \rangle\_{\mathcal{S}(x),\varepsilon}. \end{split} \tag{22}$$

*Particles* **2020**, *3*

The coefficients *A*2, *A*4, *A*<sup>6</sup> are Lorentz invariants, and the relation (22) is valid for any choice of the vectors *uμ*, *aμ*. Therefore, to determine the coefficient, we can choose the vectors *uμ*, *a<sup>μ</sup>* in any form convenient for us. To determine *<sup>A</sup>*2, we choose *<sup>a</sup><sup>μ</sup>* = (0, 0, 0, <sup>|</sup>*a*|) and *<sup>u</sup><sup>μ</sup>* = (1, 0, 0, 0) and consider the components *<sup>μ</sup>* <sup>=</sup> 0, *<sup>ν</sup>* <sup>=</sup> 0, to determine *<sup>A</sup>*<sup>4</sup> we choose *<sup>a</sup><sup>μ</sup>* = (0, 0, <sup>|</sup>*a*|, 0) and *<sup>u</sup><sup>μ</sup>* = (1, 0, 0, 0) and consider the components *<sup>μ</sup>* <sup>=</sup> 3, *<sup>ν</sup>* <sup>=</sup> 3, and to determine *<sup>A</sup>*<sup>6</sup> we choose *<sup>a</sup><sup>μ</sup>* = (0, 0, 0, <sup>|</sup>*a*|) and *u<sup>μ</sup>* = (1, 0, 0, 0) and consider the components *μ* = 3, *ν* = 3. As a result, we obtain

$$\begin{split} A\_{2} &=& -\frac{1}{4!} \int\_{0}^{|\beta|} d\tau\_{1} d\tau\_{2} d\tau\_{3} d\tau\_{4} \langle T\_{\tau} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{1}u} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{2}u} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{3}u} \mathring{\mathcal{K}}^{40}\_{-i\tau\_{4}u} \mathring{\mathcal{I}}^{40}(0) \rangle\_{\beta(\mathbf{x})\boldsymbol{\varepsilon}} \, \, \/ \\ A\_{4} &=& \frac{1}{4!} \int\_{0}^{|\beta|} d\tau\_{1} d\tau\_{2} d\tau\_{3} d\tau\_{4} \langle T\_{\tau} \mathring{\mathcal{K}}^{2}\_{-i\tau\_{1}u} \mathring{\mathcal{K}}^{2}\_{-i\tau\_{2}u} \mathring{\mathcal{K}}^{2}\_{-i\tau\_{3}u} \mathring{\mathcal{K}}^{33}\_{-i\tau\_{4}u} \mathring{\mathcal{I}}^{33}(0) \rangle\_{\beta(\mathbf{x})\boldsymbol{\varepsilon}} \, \, \/ \\ A\_{6} &=& -A\_{4} + \frac{1}{4!} \int\_{0}^{|\beta|} d\tau\_{1} d\tau\_{2} d\tau\_{3} d\tau\_{4} \langle T\_{\tau} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{1}u} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{2}u} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{3}u} \mathring{\mathcal{K}}^{3}\_{-i\tau\_{4}u} \mathring{\mathcal{I}}^{33}(0) \rangle\_{\beta(\mathbf{x})\boldsymbol{\varepsilon}} \, \, \, \/ \end{split} \tag{23}$$

We now use the representation of the boost operator through the energy-momentum tensor. According to (5) and (12), we have

$$\begin{aligned} \mathcal{K}\_{-i\tau u}^{3} &= \mathcal{J}\_{-i\tau u}^{03} = \int d^3 \mathbf{x} (-1) \mathbf{x}^3 \mathcal{I}^{00} (\tau, \mathbf{x}) \,, \\ \mathcal{K}\_{-i\tau u}^{2} &= \mathcal{J}\_{-i\tau u}^{02} = \int d^3 \mathbf{x} (-1) \mathbf{x}^2 \hat{\mathcal{I}}^{00} (\tau, \mathbf{x}) \,, \end{aligned} \tag{24}$$

Substituting (24) into (23), we come to the need of calculating quantities of the form

$$\begin{split} \mathbb{C}^{a\_{1}a\_{2}|a\_{3}a\_{4}|a\_{5}a\_{6}|a\_{7}a\_{8}|a\_{9}a\_{10}|ijkl} &= \int\_{0}^{|\mathcal{S}|} d\tau\_{x}d\tau\_{y}d\tau\_{z}d\tau\_{f}d^{3}x d^{3}y d^{3}z d^{3}f \\ \geq \mathbf{x}^{i}y^{j}z^{k}f^{l}\langle T\_{\tau}\mathcal{T}^{a\_{1}a\_{2}}(\tau\_{x},\mathbf{x})\mathcal{T}^{a\_{3}a\_{4}}(\tau\_{y},\mathbf{y})\mathcal{T}^{a\_{5}a\_{6}}(\tau\_{z},\mathbf{z})\mathcal{T}^{a\_{7}a\_{8}}(\tau\_{f},\mathbf{f})\mathcal{T}^{a\_{9}a\_{10}}(0)\rangle\_{\mathcal{S}(\mathbf{x}),\mathcal{L}} .\end{split} \tag{25}$$

In particular, from (23), we have

$$A\_2 = \frac{1}{4!} \mathcal{C}^{0|00\rangle|00\rangle|00\rangle|00\rangle|333}, \ A\_4 = \frac{1}{4!} \mathcal{C}^{0|00\rangle|00\rangle|00\rangle|33|2222}, \ A\_6 = -A\_4 + \frac{1}{4!} \mathcal{C}^{0|00\rangle|00\rangle|00\rangle|33|333} \,. \tag{26}$$

Next, we will focus on calculating the coefficient in energy *A*2; the remaining coefficients can be calculated by analogy.

We represent the energy-momentum tensor (17) in a split form

$$\begin{array}{rcl}\mathcal{T}^{a\notin}(X) &=& \lim\_{\chi\_{1},\chi\_{2}\to X} \mathcal{D}\_{ab}^{a\notin}(\partial\_{X\_{1}},\partial\_{X\_{2}})\Psi\_{a}(X\_{1})\Psi\_{b}(X\_{2}),\\\mathcal{D}\_{ab}^{a\notin}(\partial\_{X\_{1}},\partial\_{X\_{2}}) &=& \frac{i^{\delta\_{0x}+\delta\_{0\beta}}}{4}[\tilde{\gamma}\_{ab}^{a}(\partial\_{X\_{2}}-\partial\_{X\_{1}})^{\beta}+\tilde{\gamma}\_{ab}^{\beta}(\partial\_{X\_{2}}-\partial\_{X\_{1}})^{a}],\end{array} \tag{27}$$

and substitute it in (26). As a result, we get for the corresponding correlator

$$
\begin{split}
\langle\!\langle\!T\_{\mathsf{T}}\hat{\mathcal{I}}^{00}(X)\hat{\mathcal{I}}^{00}(Y)\hat{\mathcal{I}}^{00}(Z)\hat{\mathcal{I}}^{00}(F)\hat{\mathcal{I}}^{00}(0)\rangle\_{\beta(\mathbf{x}),\mathcal{L}} &= \lim\_{\begin{subarray}{c}\mathcal{X}\_{1},\mathcal{X}\_{2}\rightarrow\mathcal{X} \\ \mathcal{X}\_{1},\mathcal{X}\_{2}\rightarrow\mathcal{X} \\ \mathcal{X}\_{1},\mathcal{X}\_{2}\rightarrow\mathcal{X} \\ \mathcal{Y}\_{1},\mathcal{Y}\_{2}\rightarrow\mathcal{Y} \end{subarray}} \mathcal{D}^{00}\_{a\_{1}a\_{2}}(\mathcal{O}\_{\mathbf{X}\_{1}})\_{\mathcal{E}} \\ \times \langle\!\mathcal{D}^{00}\_{a\_{3}a\_{4}}(\mathcal{O}\_{\mathbf{Y}\_{1}},\mathcal{O}\_{\mathbf{Y}\_{2}})\mathcal{D}^{00}\_{a\_{3}a\_{4}}(\mathcal{O}\_{\mathbf{Z}\_{1}},\mathcal{O}\_{\mathbf{Z}\_{2}})\mathcal{D}^{00}\_{a\_{4}a\_{5}}(\mathcal{O}\_{\mathbf{H}\_{1}},\mathcal{O}\_{\mathbf{Z}\_{2}})\mathcal{D}^{00}\_{a\_{6}a\_{10}}(\mathcal{O}\_{\mathbf{H}\_{1}},\mathcal{O}\_{\mathbf{H}\_{2}})\langle\!\!\langle\!\mathcal{T}\_{\mathsf{T}}\mathcal{R}\_{a\_{1}}(\mathcal{X}\_{1})\!\!\!\!\langle\!\mathcal{R}\_{a\_{2}}(\mathcal{X}\_{2})\!\!\!\rangle\_{\beta(\mathbf{x}),\mathcal{X}}\rangle\rangle\_{\beta(\mathbf{x}),\mathcal{X}}.
\end{split}
$$

*Particles* **2020**, *3*

When calculating the correlator with 10 Dirac fields of the form (28), it is necessary to use an analogue of Wick theorem for field theory at finite temperatures. Then, the five-point correlator in (28) leads to the product of mean values of quadratic combinations of Dirac fields, that is, thermal propagators. For short, we denote Ψ*an* → *n*, and Ψ*an* → *n*¯ and omit *T<sup>τ</sup>* and the index *β*(*x*). Then, after extraction on the connected part in (28) according to Wick theorem, we obtain 24 terms

$$\begin{split} & (\mathbf{7}\_1 \mathbf{\bar{7}}\_{a1} (\mathbf{1}\_1) \mathbf{\bar{9}}\_{a2} (\mathbf{2}\_2) \mathbf{\bar{7}}\_{a3} (\mathbf{1}\_1) \mathbf{\bar{9}}\_{a4} (\mathbf{1}\_2) \mathbf{\bar{9}}\_{a5} (\mathbf{2}\_1) \mathbf{\bar{9}}\_{a6} (\mathbf{2}\_2) \mathbf{\bar{9}}\_{a7} (\mathbf{1}\_1) \mathbf{\bar{9}}\_{a8} (\mathbf{2}\_2) \\ & \mathbf{\bar{9}}\_{a9} (\mathbf{1}\_1) \mathbf{\bar{9}}\_{a0} (\mathbf{2}\_{10}) \mathbf{\bar{9}}\_{\bar{9}} (\mathbf{2}\_2) &= \langle 12 \mathbf{\bar{4}} \langle 2 \rangle \langle 3 \rangle \langle \mathbf{3} \rangle \langle \mathbf{8} \rangle \langle \mathbf{7} \rangle \langle \mathbf{1} \rangle \\ & + \langle 14 \rangle \langle 27 \rangle \langle 3 \rangle \langle \mathbf{5} \rangle (\mathbf{1} \mathbf{0}) \mathbf{\bar{8}} \rangle + \langle \overline{1} \langle 14 \rangle \langle 2 \rangle \langle \langle \mathbf{8} \rangle \langle \mathbf{1} \mathbf{0} \rangle \langle \mathbf{6} \rangle \rangle + \langle \overline{1} \langle \mathbf{4} \rangle \langle \mathbf{2} \rangle \langle \langle \mathbf{3} \rangle \langle \mathbf{4} \rangle \langle \mathbf{7} \rangle \langle \mathbf{1} \rangle \\ & + \langle 14 \rangle \langle 27 \rangle \langle \langle \mathbf{1} \mathbf{0} \rangle \langle \mathbf{8} \rangle \langle \mathbf{6} \rangle - \langle \overline{1} \langle \mathbf{4} \rangle \langle \mathbf{2} \rangle \langle \langle \mathbf{3} \rangle \langle \mathbf{6} \rangle \langle \$$

where signs correspond to the number of permutations of anti-commuting fields. Thermal propagators have a standard form [13,33]

$$\begin{split} \mathcal{G}\_{\mathfrak{a}\_{1}\mathfrak{a}\_{2}}(X\_{1},X\_{2}) &= \langle \!\langle \!\!/\tau\_{\mathtt{T}}\Psi\_{\mathfrak{a}\_{1}}(X\_{1})\!\!\!/\mathfrak{V}\_{\mathfrak{a}\_{2}}(X\_{2})\!\rangle\_{\beta(\mathtt{x})} = \sum\_{\mu}^{\mathfrak{f}} e^{i\varGamma^{+}(X\_{1}-X\_{2})} (-i\varGamma^{+}\_{\mu}\tilde{\gamma}\_{\varmu} + m)\_{\mathfrak{a}\_{1}\mathfrak{a}\_{2}} \Delta(P^{+}) \rangle, \\ \mathcal{G}\_{\mathfrak{a}\_{1}\mathfrak{a}\_{2}}(X\_{1},X\_{2}) &= \langle T\_{\mathsf{T}}\Psi\_{\mathfrak{a}\_{1}}(X\_{1})\!\!\!/\mathfrak{V}\_{\mathfrak{a}\_{2}}(X\_{2})\rangle\_{\beta(\mathtt{x})} = -\langle T\_{\mathsf{T}}\Psi\_{\mathfrak{a}\_{2}}(X\_{2})\!\!\!/\mathfrak{V}\_{\mathfrak{a}\_{1}}(X\_{1})\rangle\_{\beta(\mathtt{x})} \\ &= -\sum\_{\mathsf{P}}^{\mathfrak{f}} e^{i\varGamma^{-}(X\_{1}-X\_{2})} (i\varGamma^{-}\_{\mu}\tilde{\gamma}\_{\varmu} + m)\_{\mathfrak{a}\_{2}\mathfrak{a}\_{1}} \Delta(P^{-}) \rangle, \end{split} \tag{30}$$

where integration over the three-dimensional components of the momentum and summation over the Matsubara frequencies of fermion field appear. The following notations are used in (30): *P*± = (*p*± *<sup>n</sup>* , **p**), *p*± *<sup>n</sup>* = *π*(2*n* + 1)/|*β*| ± *iμ*, *n* = 0, ±1, ±2, ··· , *X* = (*τ*, **x**), <sup>∑</sup> *<sup>P</sup>* <sup>=</sup> <sup>1</sup> <sup>|</sup>*β*<sup>|</sup> <sup>∑</sup><sup>∞</sup> *n*=−∞ *<sup>d</sup>*<sup>3</sup> *<sup>p</sup>* (2*π*)<sup>3</sup> , and <sup>Δ</sup>(*P*) = <sup>1</sup> *<sup>P</sup>*2+*m*<sup>2</sup> , where the square is taken with the Euclidean metric, as also in *P*± *<sup>μ</sup> <sup>γ</sup>*˜*<sup>μ</sup>* <sup>=</sup> /*P*<sup>±</sup> (unlike *<sup>P</sup>*+(*X*<sup>1</sup> <sup>−</sup> *<sup>X</sup>*2) according to [33]). Since we consider mass-less field at zero chemical potential, the mass and chemical potential must be set equal to zero *m* = 0, *μ* = 0. Nevertheless, we retain the notation *P*±, bearing in mind the possibility of generalization to the case with nonzero chemical potential in the future.

Next, substitute (30) in (29). We will describe the calculations for the first term in (29), while all other terms can be calculated by analogy

<sup>−</sup> lim *<sup>X</sup>*1,*X*2→*<sup>X</sup> Y*1,*Y*2→*Y Z*1,*Z*2→*Z F*1,*F*2→*F H*1,*H*2→*H*=0 <sup>D</sup><sup>00</sup> *<sup>a</sup>*1*a*<sup>2</sup> (*∂X*<sup>1</sup> , *<sup>∂</sup>X*<sup>2</sup> )D<sup>00</sup> *<sup>a</sup>*3*a*<sup>4</sup> (*∂Y*<sup>1</sup> , *<sup>∂</sup>Y*<sup>2</sup> )D<sup>00</sup> *<sup>a</sup>*5*a*<sup>6</sup> (*∂Z*<sup>1</sup> , *<sup>∂</sup>Z*<sup>2</sup> )D<sup>00</sup> *<sup>a</sup>*7*a*<sup>8</sup> (*∂F*<sup>1</sup> , *<sup>∂</sup>F*<sup>2</sup> )D<sup>00</sup> *<sup>a</sup>*9*a*<sup>10</sup> (*∂H*<sup>1</sup> , *∂H*<sup>2</sup> ) <sup>×</sup>*G*¯ *<sup>a</sup>*1*a*<sup>4</sup> (*X*1,*Y*2)*Ga*2*a*<sup>9</sup> (*X*2, *<sup>H</sup>*1)*G*¯ *<sup>a</sup>*3*a*<sup>6</sup> (*Y*1, *<sup>Z</sup>*2)*G*¯ *<sup>a</sup>*5*a*<sup>8</sup> (*Z*1, *<sup>F</sup>*2)*G*¯ *<sup>a</sup>*7*a*<sup>10</sup> (*F*1, *<sup>H</sup>*2) = − ∑ {*P*,*Q*,*K*,*R*,*L*} *e* <sup>−</sup>*i***p**(**x**−**y**)−*i***qx**−*i***k**(**y**−**z**)−*i***r**(**z**−**f**)−*i***lf***e ip*− *<sup>n</sup>* (*τx*−*τy*)+*iq*<sup>+</sup> *<sup>n</sup> τx*+*ik*<sup>−</sup> *<sup>n</sup>* (*τy*−*τ<sup>z</sup>* )+*ir*<sup>−</sup> *<sup>n</sup>* (*τz*−*τf*)+*il*<sup>−</sup> *<sup>n</sup> τ<sup>f</sup>* <sup>×</sup>Δ(*P*−)Δ(*Q*+)Δ(*K*−)Δ(*R*−)Δ(*L*−) ×tr (−*i*/*L*−)D00(*iL*−, <sup>−</sup>*iR*−)(−*iR*/−)D00(*iR*−, <sup>−</sup>*iK*−)(−*iK*/−) ×D00(*iK*−, <sup>−</sup>*iP*−)(−*i*/*P*−)D00(*iP*−, *iQ*+)(−*iQ*/+)D00(−*iQ*+, <sup>−</sup>*iL*−) , (31)

where it was necessary to arrange all the matrices under the trace in accordance with the order of the spinor indices. To calculate (31), it is necessary to find a trace of the form

$$\times \left[ \mathcal{J}\_{\mathbb{T}} \mathcal{D}^{00} (P\_2, P\_3) \mathcal{J}\_{\mathbb{4}} \mathcal{D}^{00} (P\_5, P\_6) \mathcal{J}\_{\mathbb{T}} \mathcal{D}^{00} (P\_8, P\_9) P\_{\mathbb{f}0} \mathcal{D}^{00} (P\_{11}, P\_{12}) P\_{\mathbb{f}3} \mathcal{D}^{00} (P\_{14}, P\_{15}) \right] . \tag{32}$$

The subsequent calculations are more convenient to carry out using special software applications. Calculation (32) requires finding the trace of 10 Euclidean Dirac matrices

$$\text{tr}\left[\tilde{\gamma}^{\mathbb{A}\_1}\tilde{\gamma}^{\mathbb{A}2}\tilde{\gamma}^{\mathbb{A}2}\tilde{\gamma}^{\mathbb{A}3}\tilde{\gamma}^{\mathbb{A}4}\tilde{\gamma}^{\mathbb{A}5}\tilde{\gamma}^{\mathbb{A}6}\tilde{\gamma}^{\mathbb{A}7}\tilde{\gamma}^{\mathbb{A}8}\tilde{\gamma}^{\mathbb{A}9}\tilde{\gamma}^{\mathbb{A}9}\tilde{\gamma}^{\mathbb{A}10}\right].\tag{33}$$

Using the definition (18), this trace can be easily transformed to the trace of ordinary Dirac matrices, which can be found using standard methods. We denote the trace in (32) as *A*(*P*1, *P*2, *P*3, *P*4, *P*5, *P*6, *P*7, *P*8, *P*9, *P*10, *P*11, *P*12, *P*13, *P*14, *P*15). Then (31) will be presented in the form

$$\begin{split} & - \int \frac{d^3 p d^3 q d^3 k d^3 r d^3 l}{(2\pi)^{15}} e^{-i\mathbf{p}(\mathbf{x}-\mathbf{y}) - i\mathbf{q}\mathbf{x} - i\mathbf{k}(\mathbf{y}-\mathbf{z}) - i\mathbf{r}(\mathbf{z}-\mathbf{f}) - i\mathbf{l}\mathbf{f}} \\ & \times \sum\_{p,q,\mu,k,\mu,l,\mu} \frac{1}{|\beta|^5} e^{i p^+\_{\text{a}}(\mathbf{\tau}\_{\mathbf{x}}-\mathbf{\tau}\_{\mathbf{y}}) + i q^+\_{\text{a}}\tau\_{\mathbf{z}} + i\mathbf{k}^-\_{\text{a}}(\tau\_{\mathbf{y}}-\tau\_{\mathbf{z}}) + i\mathbf{r}^-\_{\text{n}}(\tau\_{\mathbf{z}}-\tau\_{\mathbf{y}}) + i\mathbf{l}^-\_{\text{n}}\tau\_{\mathbf{y}} \\ & \times \Delta(P^-)\Delta(Q^+)\Delta(K^-)\Delta(R^-)\Delta(L^-) \\ & \times A(-i\mathbf{L}^-, i\mathbf{L}^-, -i\mathbf{R}^-, -i\mathbf{R}^-, i\mathbf{R}^-, -i\mathbf{K}^-, -i\mathbf{K}^-, -i\mathbf{K}^-, i\mathbf{K}^-, -i\mathbf{P}^-, -i\mathbf{P}^-, -i\mathbf{P}^-, i\mathbf{P}^-, i\mathbf{Q}^+) \\ & -i\mathbf{Q}^+, -i\mathbf{Q}^+, -i\mathbf{L}^-). \end{split} \tag{34}$$

Next, one needs to sum over the Matsubara frequencies in (34) using the relation

$$\frac{1}{|\beta|}\sum\_{\omega\_n} \frac{(\omega\_n \pm i\mu)^k e^{i(\omega\_n \pm i\mu)\tau}}{(\omega\_n \pm i\mu)^2 + E^2} = \frac{1}{2E} \sum\_{s=\pm 1} (-isE)^k e^{\tau s E} [\theta(-s\tau) - n\_F(E \pm s\mu)],\tag{35}$$

where *<sup>E</sup>* = **p**<sup>2</sup> + *<sup>m</sup>*2, *nF*(*E*) = 1/(<sup>1</sup> + *<sup>e</sup>E*/*T*) is the Fermi distribution, and *<sup>θ</sup>* is the Heaviside theta function. Again, we can take *m* = 0, *μ* = 0. As a result, we obtain

*Particles* **2020**, *3*

− *d*<sup>3</sup> *pd*3*qd*3*kd*3*rd*3*l* (2*π*)<sup>15</sup> *<sup>e</sup>* −*i***p**(**x**−**y**)−*i***qx**−*i***k**(**y**−**z**)−*i***r**(**z**−**f**)−*i***lf** <sup>×</sup> <sup>1</sup> <sup>32</sup>*EpEqEkErEl* <sup>∑</sup> *<sup>s</sup>*1,*s*2,*s*3,*s*4,*s*<sup>5</sup> *e* (*τx*−*τy*)*s*1*Ep*+*τxs*2*Eq*+(*τy*−*τ<sup>z</sup>* )*s*3*Ek*+(*τz*−*τf*)*s*4*Er*+*τ<sup>f</sup> s*5*El* ×*A*(−*i L*, *i L*, −*iR*, −*iR*, *iR*, −*iK*, −*iK*, *iK*, −*iP*, −*iP*, *iP*, *iQ*, −*iQ*, −*iQ*, −*i L*) × % *θ* −*s*1(*τ<sup>x</sup>* − *τy*) − *nF* % *Ep* && %*<sup>θ</sup>* [−*s*2] <sup>−</sup> *nF* % *Eq* && × % *θ* −*s*3(*τ<sup>y</sup>* − *τz*) − *nF* (*Ek*) & *θ* −*s*4(*τ<sup>z</sup>* − *τf*) − *nF* (*Er*) (*θ* [−*s*5] − *nF* (*El*)) . (36)

Here, following [13], the notations *P* = *P*(*s*1)=(−*is*1*Ep*, **p**), *Q* = *Q*(*s*2), ··· are introduced. We return now to the formula for *A*<sup>2</sup> (26) with spatial integrals and calculate the contribution of the term (36). This contribution has the form

$$A\_{2} = \int \frac{d\tau\_{\mathbf{x}}d\tau\_{\mathbf{y}}d\tau\_{\mathbf{z}}d\tau\_{f}d^{3}zd^{3}zd^{3}d^{3}f d^{3}p d^{3}d^{3}d^{3}d^{3}l}{4!(2\pi)^{15}}e^{-i\mathbf{p}(\mathbf{x}-\mathbf{y})-i\mathbf{q}\mathbf{x}-i\mathbf{k}(\mathbf{y}-\mathbf{z})-i\mathbf{r}(\mathbf{z}-\mathbf{f})-i\mathbf{f}} \tag{37}$$

where the ellipsis indicates the contribution of the remaining 23 terms from (29), and *D* equals to

*<sup>D</sup>* <sup>=</sup> <sup>−</sup> <sup>1</sup> <sup>32</sup>*EpEqEkErEl* <sup>∑</sup> *<sup>s</sup>*1,*s*2,*s*3,*s*4,*s*<sup>5</sup> *e* (*τx*−*τy*)*s*1*Ep*+*τxs*2*Eq*+(*τy*−*τ<sup>z</sup>* )*s*3*Ek*+(*τz*−*τf*)*s*4*Er*+*τ<sup>f</sup> s*5*El* ×*A*(−*i L*, *i L*, −*iR*, −*iR*, *iR*, −*iK*, −*iK*, *iK*, −*iP*, −*iP*, *iP*, *iQ*, −*iQ*, −*iQ*, −*i L*) × % *θ* −*s*1(*τ<sup>x</sup>* − *τy*) − *nF* % *Ep* && %*<sup>θ</sup>* [−*s*2] <sup>−</sup> *nF* % *Eq* && × % *θ* −*s*3(*τ<sup>y</sup>* − *τz*) − *nF* (*Ek*) & *θ* −*s*4(*τ<sup>z</sup>* − *τf*) − *nF* (*Er*) (*θ* [−*s*5] − *nF* (*El*)) . (38)

Next, one needs to rewrite the product of spatial coordinates in the integral through derivatives using the formula

 *d*<sup>3</sup> *pd*3*qd*3*kd*3*rd*3*xd*3*yd*3*zd*<sup>3</sup> *f F*(**p**, **q**, **k**,**r**, **l**)*e* <sup>−</sup>*i***p**(**x**−**y**)−*i***qx**−*i***k**(**y**−**z**)−*i***r**(**z**−**f**)−*i***lf***x*3*y*3*z*<sup>3</sup> *f* <sup>3</sup> = (2*π*)<sup>12</sup> *d*<sup>3</sup> *p* <sup>−</sup> *<sup>∂</sup>*<sup>3</sup> *<sup>∂</sup>q*3*∂l*3*∂p*3*∂r*<sup>3</sup> <sup>−</sup> *<sup>∂</sup>*<sup>3</sup> *∂q*3*∂l*3*∂p*3*∂l*<sup>3</sup> + *∂*3 *<sup>∂</sup>q*3*∂l*3*∂r*3*∂q*<sup>3</sup> <sup>+</sup> *∂*3 *∂q*3*∂q*3*∂l*3*∂l*<sup>3</sup> *F*(**p**, **q**, **k**,**r**, **l**) # # # # # **l**=**p r**=**p k**=**p q**=−**p** , (39)

resulting from integration by parts and properties of the delta function. After that, (37) is converted to the form

$$\begin{split} A\_{2} &= \quad \frac{1}{4!(2\pi)^{3}} \int d\tau\_{3}d\tau\_{3}d\tau\_{5}d\tau\_{f}d^{3}p \left( -\frac{\partial^{3}}{\partial q^{3}\partial l^{3}\partial p^{3}\partial r^{3}} - \frac{\partial^{3}}{\partial q^{3}\partial l^{3}\partial p^{3}\partial l^{3}} \right. \\ &\left. + \frac{\partial^{3}}{\partial q^{3}\partial l^{3}\partial r^{3}\partial q^{3}} + \frac{\partial^{3}}{\partial q^{3}\partial q^{3}\partial l^{3}\partial l^{3}} \right) D(\mathbf{p}, \mathbf{q}, \mathbf{k}, \mathbf{r}, \mathbf{l}) \bigg| \begin{array}{ll} \frac{1-\mathbf{p}}{\mathbf{r}-\mathbf{p}} & + \cdots \\ \frac{\mathbf{k}-\mathbf{p}}{\mathbf{q}-\mathbf{p}} & + \cdots \end{array} . \tag{40} \end{split} \tag{40}$$

*Particles* **2020**, *3*

Now, it remains to integrate over the imaginary time and also over the last momentum, which can be done directly in spherical coordinates *<sup>d</sup>*<sup>3</sup> *<sup>p</sup>* <sup>=</sup> <sup>|</sup>**p**<sup>|</sup> <sup>2</sup>*d*|**p**<sup>|</sup> sin(*θ*)*dθdφ*. The sequence of actions in this case, from the point of view of calculation speed, will be most convenient as follows: first one needs to make differentiations with respect to the four momentum variables in (40), then make the corresponding changes of the variables following from the delta functions, then sum over the indices *sn* from (38), then integrate over the angles in *d*<sup>3</sup> *p*, and then integrate over imaginary time variables, which requires careful handling of theta functions. The transformations with each of the 24 terms in (29) can be performed independently and using parallel computing tools. We do not give the described intermediate steps, since they are most conveniently performed using the program, and the intermediate formulas themselves are extremely long, while the calculations themselves are not difficult from a mathematical point of view and are done directly. As a result, we obtain the following integral

$$\begin{split} A\_{2} &= \quad \int\_{0}^{\infty} d\hat{p} e^{\frac{6\tilde{p}}{2}} \hat{p}^{3} \Big( 5600 \vec{\tilde{p}} \left( 49 \hat{p}^{2} - 95 \right) \cosh\left( \frac{\tilde{p}}{2} \right) + 2016 \vec{p} \left( 25 - 119 \hat{p}^{2} \right) \cosh\left( \frac{3\tilde{p}}{2} \right) \\ &+ 53200 \left( \sinh\left( \frac{3\tilde{p}}{2} \right) - 11 \sinh\left( \frac{\tilde{p}}{2} \right) \right) \cosh^{4}\left( \frac{\tilde{p}}{2} \right) + \hat{p} \left( -224 \left( \vec{p}^{2} + 25 \right) \cosh\left( \frac{7\tilde{p}}{2} \right) \right) \\ &+ 224 \left( 119 \hat{p}^{2} + 575 \right) \cosh\left( \frac{5\tilde{p}}{2} \right) + 18 \vec{p} \sinh\left( \frac{\tilde{p}}{2} \right) \left( -5786 \hat{p}^{2} + \left( \vec{p}^{2} + 210 \right) \cosh\left( 3\tilde{p} \right) \right) \\ &- 6 \left( 41 \hat{p}^{2} + 1890 \right) \cosh\left( 2\vec{p} \right) + 3 \left( 1349 \hat{p}^{2} + 9450 \right) \cosh\left( \vec{p} \right) \\ &+ 39900 \left( \right) \left( (50400 \pi^{2} \left( \vec{e}^{\theta} + 1 \right)^{9} )^{-1}, \tag{41} \tag{42} \end{split}$$

where the dimensionless variable *p*˜ = |**p**|/*T* was introduced. This integral converges and can be found analytically:

$$A\_2 \quad = \quad -\frac{17}{960\pi^2} \,\text{.}\tag{42}$$

Repeating the entire calculation algorithm for the coefficients *A*4, *A*<sup>6</sup> in (26), we obtain at *m* = 0

$$A\_4 = -\frac{17}{2880\pi^2}, \quad A\_6 = 0.\tag{43}$$

Saving the mass in all formulas, in particular, in the propagators (30), we get more complicated expressions for the coefficients at finite mass given in the Appendix A.

#### **4. Discussion**

In the previous section, we described the details of the calculation of the corrections of the fourth order in acceleration to the energy-momentum tensor of the Dirac field, first obtained in [23]. Taking into account (42) and (43), we obtain the next formula for the energy-momentum tensor at *m* = 0

$$\left< \langle \mathcal{T}^{\mu\nu} \rangle \right> = \left( \frac{7\pi^2 T^4}{60} + \frac{T^2 |a|^2}{24} - \frac{17|a|^4}{960\pi^2} \right) \mu^{\mu} u^{\nu} - \left( \frac{7\pi^2 T^4}{180} + \frac{T^2 |a|^2}{72} - \frac{17|a|^4}{2880\pi^2} \right) \Delta^{\mu\nu} + \mathcal{O}(a^6), \tag{44}$$

where the notation <sup>|</sup>*a*<sup>|</sup> <sup>=</sup> −*aμa<sup>μ</sup>* is used.

As discussed in [22,24,25,27], the mean value of the energy-momentum tensor calculated in this way should vanish at the proper temperature equal to Unruh temperature. Since the energy-momentum tensor is normalized with respect to the Minkowski vacuum, such a vanishing is a direct consequence of the Unruh effect—an accelerated medium with Unruh temperature corresponds to the Minkowski vacuum. It is easy to see that energy-momentum tensor (44) satisfies this condition following from the Unruh effect

$$
\langle \hat{T}^{\mu\nu} \rangle (T = T\_{ll}) = 0 \,. \tag{45}
$$

Moreover, as discussed in [26], from the presentation of the result (44) in the form of Sommerfeld integrals, as well as comparison with the field theory in a space with a conical singularity, it follows that the calculated fourth order of perturbation theory is maximal; that is, <sup>O</sup>(*a*6) = 0 at least at *<sup>T</sup>* <sup>&</sup>gt; *TU* [26]. Thus, Equation (44) is an exact non-perturbative formula in this region.

We also note that expression (44) can be obtained from the point of view of another approach, where field theory in a space with a conical singularity is considered [25,27]. As discussed in [26], this indicates the duality of the statistical and geometrical approaches to the description of accelerated media.

#### **5. Conclusions**

The Zubarev density operator provides a powerful fundamental theoretical method for studying quantum-field effects in the accelerated medium. This makes it possible to obtain information about such a medium from the point of view of an inertial observer and there is no need to go to the curvilinear coordinates of the accelerated frame and consider the features of nontrivial space with a boundary. All effects can be calculated in ordinary flat space described by the Minkowski metric using standard Green functions at finite temperature. In this case, the effects of acceleration are calculated in a regular way in the framework of perturbation theory with the boost operator. However, it is possible to obtain exact non-perturbative expressions in the chiral limit, since the first few orders of the perturbation theory are to give a complete perturbative series.

In particular, earlier in [23], the Unruh effect for fermions was demonstrated by calculating fourth-order quantum corrections. In the language of the statistical approach with the Zubarev operator, the Unruh effect should lead to the vanishing of the energy-momentum tensor at the proper temperature equal to the Unruh temperature. Thus, the Zubarev density operator allows one to obtain information about the effects associated with the occurrence of an event horizon in an accelerated system and the radiation associated with it.

In more usual formulation or from the point of view of modern developments in the quantum optics [8–10], the Unruh effect should be manifested in the thermal distribution of photons with Unruh temperature. However, it can be shown that the formula we obtained (44) also contains such a distribution with the Unruh temperature [23].

In this paper, we described the details of the calculations of the coefficients with acceleration in the energy-momentum tensor given in [23], focusing on the calculation of the quantum correction to the energy density. The calculation of this correction consists in finding the mean value of the product of the boost operators and operator of the energy-momentum tensor. Applying Wick theorem, one can transform the average of the product of operators to the product of five thermal propagators. Each of the propagators adds one summation over the Matsubara frequencies and a three-dimensional integral over the momentum, and also each boost operator adds three-dimensional integral over the coordinate and one integral over the imaginary time. The procedure for calculating these sums and integrals is described. In addition, expressions for the coefficients at a finite mass are given.

*Particles* **2020**, *3*

The effects of acceleration we are discussing are of interest from the experimental point of view, in particular, in heavy-ion collisions, where large acceleration can occur. A systematic study of the effects of acceleration requires calculating the acceleration resulting from particle collisions, similar to calculating the vorticity [34–36]. Since the vorticity turns out to be significant in the collision of particles, acceleration, being another combination of derivatives, is also expected to affect the observables. We predict that the effects of acceleration should be significant at early stages of the collision, when the system is not yet fully thermalized and the terms with acceleration are not suppressed with respect to temperature. In this case, non-equilibrium processes can arise that are associated with instability at the Unruh temperature, which were discussed in [26]. One can also make a prediction that the discussed electron-ion collider (EIC) can become a good laboratory for studying effects of acceleration [37]. An elementary particle like an electron, colliding with an ion, behaves like a wave, which allows us to separate the effects of acceleration from the effects of vorticity.

**Author Contributions:** Conceptualization, G.P., O.T. and V.Z.; investigation, G.P., O.T. and V.Z.; draft preparation, G.P.; supervision, O.T. and V.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported by Russian Science Foundation Grant No 16-12-10059.

**Acknowledgments:** Useful discussions with F. Becattini are gratefully acknowledged.

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

#### **Appendix A. The Coefficients** *a*<sup>4</sup> **at Finite Mass**

The coefficient *A*<sup>2</sup> at a finite mass is described by the expression

*A*<sup>2</sup> = ∞ 0 *dp*˜*p*˜ 2*e* 9*E*˜ *p* 2 *E*˜ *p* 9 <sup>51450</sup> <sup>−</sup> <sup>15619</sup>*m*˜ <sup>2</sup> *p*˜ <sup>4</sup> <sup>+</sup> <sup>175</sup> <sup>2450</sup>*m*˜ <sup>2</sup> <sup>−</sup> <sup>361</sup> *p*˜ 2 +175*m*˜ <sup>2</sup> <sup>392</sup>*m*˜ <sup>2</sup> <sup>−</sup> <sup>285</sup> − 140571*p*˜ 6 sinh  *E*˜ *p* 2 + 27*E*˜ *<sup>p</sup>* <sup>27</sup> <sup>53</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>490</sup> *p*˜ 4 <sup>+</sup><sup>175</sup> <sup>70</sup>*m*˜ <sup>2</sup> <sup>−</sup> <sup>19</sup> *p*˜ <sup>2</sup> + 35*m*˜ <sup>2</sup> <sup>56</sup>*m*˜ <sup>2</sup> <sup>−</sup> <sup>75</sup> + 1431*p*˜ 6 sinh  3*E*˜ *<sup>p</sup>* 2 <sup>−</sup>*E*˜ *<sup>p</sup>* 9 <sup>247</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>11550</sup> *p*˜ <sup>4</sup> <sup>+</sup> <sup>175</sup> <sup>550</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>133</sup> *p*˜ <sup>2</sup> + 175*m*˜ <sup>2</sup> <sup>88</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>105</sup> +2223*p*˜ 6 sinh  5*E*˜ *<sup>p</sup>* 2 + *E*˜ *<sup>p</sup>* 9 *<sup>m</sup>*˜ <sup>2</sup> <sup>+</sup> <sup>210</sup> *p*˜ <sup>4</sup> <sup>+</sup> <sup>175</sup> <sup>10</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>19</sup> *p*˜ 2 +35*m*˜ <sup>2</sup> <sup>8</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>75</sup> + 9*p*˜ 6 sinh  7*E*˜ *<sup>p</sup>* 2 − 8 *m*˜ <sup>2</sup> + *p*˜ 2 5*m*˜ <sup>2</sup> 2*p*˜ <sup>2</sup> <sup>+</sup> <sup>63</sup> +28*p*˜ 2 *p*˜ <sup>2</sup> <sup>+</sup> <sup>25</sup> cosh  7*E*˜ *<sup>p</sup>* 2 <sup>−</sup> <sup>504</sup> *m*˜ <sup>2</sup> + *p*˜ 2 5*m*˜ <sup>2</sup> 34*p*˜ <sup>2</sup> <sup>−</sup> <sup>9</sup> +476*p*˜ <sup>4</sup> <sup>−</sup> <sup>100</sup>*p*˜ 2 cosh  3*E*˜ *<sup>p</sup>* 2 <sup>+</sup> <sup>56</sup> *m*˜ <sup>2</sup> + *p*˜ 2 5*m*˜ <sup>2</sup> 34*p*˜ <sup>2</sup> <sup>+</sup> <sup>207</sup> +476*p*˜ <sup>4</sup> + 2300*p*˜ 2 cosh  5*E*˜ *<sup>p</sup>* 2 <sup>+</sup> <sup>1400</sup> *m*˜ <sup>2</sup> + *p*˜ 2 *m*˜ <sup>2</sup> 70*p*˜ <sup>2</sup> <sup>−</sup> <sup>171</sup> +4*p*˜ 2 49*p*˜ <sup>2</sup> <sup>−</sup> <sup>95</sup> cosh  *E*˜ *p* 2 (50400*π*<sup>2</sup> *eE*˜ *<sup>p</sup>* + 1 <sup>9</sup>*E*˜2 *<sup>p</sup>*)−<sup>1</sup> , (A1)

*Particles* **2020**, *3*

where the dimensionless quantities *<sup>m</sup>*˜ = *<sup>m</sup>*/*T*, *<sup>E</sup>*˜ *<sup>p</sup>* = *p*<sup>2</sup> + *<sup>m</sup>*2/*<sup>T</sup>* are introduced. The coefficient *<sup>A</sup>*<sup>4</sup> at a finite mass has the form

*A*<sup>4</sup> = ∞ 0 *dp*˜ *p*˜ 4*e* 9E˜ *p* 2 1960 sinh  E˜ *p* 2 cosh2 E˜ *p* 2 8*m*˜ <sup>2</sup> <sup>+</sup> <sup>15</sup> cosh % 2E˜ *<sup>p</sup>* & −8 <sup>56</sup>*m*˜ <sup>2</sup> <sup>+</sup> <sup>15</sup> cosh % E˜ *p* & <sup>+</sup> <sup>984</sup>*m*˜ <sup>2</sup> <sup>−</sup> <sup>135</sup> + *p*˜ 2 408170 − 421713*p*˜ 2 sinh  E˜ *p* 2 +27*p*˜ 2 4293*p*˜ <sup>2</sup> <sup>+</sup> <sup>11662</sup> sinh  3E˜ *<sup>p</sup>* 2 − *p*˜ 2 6669*p*˜ <sup>2</sup> <sup>+</sup> <sup>91630</sup> sinh  5E˜ *<sup>p</sup>* 2 +*p*˜ 2 27*p*˜ <sup>2</sup> <sup>+</sup> <sup>1666</sup> sinh  7E˜ *<sup>p</sup>* 2 <sup>+</sup> <sup>29400</sup> 14*p*˜ <sup>2</sup> <sup>−</sup> <sup>19</sup> E˜ *<sup>p</sup>* cosh  E˜ *p* 2 <sup>−</sup><sup>168</sup> 2*p*˜ <sup>2</sup> <sup>+</sup> <sup>35</sup> E˜ *<sup>p</sup>* cosh  7E˜ *<sup>p</sup>* 2 <sup>−</sup> <sup>10584</sup> 34*p*˜ <sup>2</sup> <sup>−</sup> <sup>5</sup> E˜ *<sup>p</sup>* cosh  3E˜ *<sup>p</sup>* 2 <sup>+</sup><sup>1176</sup> 34*p*˜ <sup>2</sup> <sup>+</sup> <sup>115</sup> E˜ *<sup>p</sup>* cosh  5E˜ *<sup>p</sup>* 2 1058400*π*<sup>2</sup> *e* E˜ *<sup>p</sup>* + 1 9E˜ *<sup>p</sup>* −<sup>1</sup> . (A2)

The coefficient *A*<sup>6</sup> is zero both for massless and massive Dirac fields

*A*<sup>6</sup> = 0 . (A3)

#### **References**


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