General equilibrium second-order hydrodynamic coefficients for free quantum fields

TL;DR

This paper calculates second-order thermal vorticity corrections to stress-energy and currents of free quantum fields, revealing their quantum origin and providing simplified Kubo formulas for these coefficients.

Contribution

It systematically derives non-dissipative second-order hydrodynamic coefficients for free boson and Dirac fields using correlators, clarifying their quantum nature.

Findings

Corrections are non-dissipative and survive at thermodynamic equilibrium.

Axial current of free Dirac field receives vorticity-proportional corrections.

Results compare favorably with previous calculations and highlight quantum origins.

Abstract

We present a systematic calculation of the corrections of the stress-energy tensor and currents of the free boson and Dirac fields up to second order in thermal vorticity, which is relevant for relativistic hydrodynamics. These corrections are non-dissipative because they survive at general thermodynamic equilibrium with non vanishing mean values of the conserved generators of the Lorentz group, i.e. angular momenta and boosts. Their equilibrium nature makes it possible to express the relevant coefficients by means of correlators of the angular-momentum and boost operators with stress-energy tensor and current, thus making simpler to determine their so-called "Kubo formulae". We show that, at least for free fields, the corrections are of quantum origin and we study several limiting cases and compare our results with previous calculations. We find that the axial current of the free Dirac field receives corrections proportional to the vorticity independently of the anomalous term.

Tables

Table 1. Table 1: The stress-energy tensor coefficients ( 27 ) for a free Boson field: the first column reports the coefficients in the massless case with μ = 0 𝜇 0 \mu=0 , the second column reports the n 𝑛 n th term of the expansion ( 74 ), the third column reports the asymptotic expansion at low temperature ( 75 ). Our result for W 𝑊 W in the massless case agrees with that obtained in ref.  moore for λ 3 subscript 𝜆 3 \lambda_{3} (see eq. ( IX ))

Table 2. Table 2: The stress-energy tensor coefficients ( 27 ) for a free Dirac field: the first column reports the coefficients in the massless case, the second column reports the coefficient of the n 𝑛 n th term of the expansion ( 74 ), the third column reports the asymptotic expansion at low temperature ( 75 ).

Table 3. Table 3: The vector current coefficients ( 57 ) for an ideal boson field: the first column is the result at m = 0 𝑚 0 m=0 and μ = 0 𝜇 0 \mu=0 , the second column is the generic terms of the series expansion ( 76 ), the third column is the asymptotic expansion for low temperature ( 77 ).

Table 4. Table 4: The vector current coefficients ( 57 ) for an ideal Dirac field: the first column is the result at m = 0 𝑚 0 m=0 , the second column is the generic terms of the series expansion ( 76 ), the third column is the asymptotic expansion for low temperature ( 77 ).

Table 5. Table 5: The coefficients of the energy momentum tensor ( 27 ), the vector current ( 57 ) and axial current ( 66 ) for the Dirac field at T = 0 𝑇 0 T=0 and finite chemical potential μ 𝜇 \mu , where b = μ / m 𝑏 𝜇 𝑚 b=\mu/m .

Equations

$$T^{\mu \nu }=(\rho +p)u^{\mu }{u}^{\nu }-g^{\mu \nu }p,$$

$$\hat{\rho }=\frac{1}{Z}\exp \left[-{\int }_{\Sigma }\mathrm{d}{\Sigma }_{\mu }\left({\hat{T}}^{\mu \nu }(x){\beta }_{\nu }(x)-\zeta (x){\hat{j}}^{\mu }(x)\right)\right],$$

$$T=\frac{1}{\sqrt{{\beta }^2}}$$

$$u_{\mu }(x)=\frac{{\beta }_{\mu }(x)}{\sqrt{{\beta }^2}},$$

$${\nabla }_{\mu }{\beta }_{\nu }+{\nabla }_{\nu }{\beta }_{\mu }=0{\nabla }_{\mu }\zeta =0,$$

$${\beta }_{\mu }(x)=b_{\mu }+{\varpi }_{\mu \nu }{x}^{\nu }\zeta =const.$$

$${\varpi }_{\mu \nu }=-\frac{1}{2}({\partial }_{\mu }{\beta }_{\nu }-{\partial }_{\nu }{\beta }_{\mu }).$$

$$\hat{\rho }=\frac{1}{Z}\exp \left[-b_{\mu }{\hat{P}}^{\mu }+\frac{1}{2}{\varpi }_{\mu \nu }{\hat{J}}^{\mu \nu }+\zeta \hat{Q}\right],$$

$${\hat{J}}^{\mu \nu }={\int }_{\Sigma }\mathrm{d}{\Sigma }_{\lambda }\left(x^{\mu }{\hat{T}}^{\lambda \nu }-x^{\nu }{\hat{T}}^{\lambda \mu }\right).$$

$$\hat{\rho }=\frac{1}{Z}\exp [-b\cdot \hat{P}+\zeta \hat{Q}],$$

$$\hat{\rho }=\frac{1}{Z}\exp \left[-\hat{H}/T_0+{\omega }_0{\hat{J}}_z/T_0\right],$$

$$⟨\hat{O}(x)⟩=\frac{1}{Z}\mathrm{tr}(\hat{\rho }\hat{O}(x){)}_{\mathrm{ren}},$$

$${\hat{J}}_x^{\mu \nu }\equiv {\int }_{\Sigma }\mathrm{d}{\Sigma }_{\lambda }\left[(y^{\mu }-x^{\mu }){\hat{T}}^{\lambda \nu }(y)-(y^{\nu }-x^{\nu }){\hat{T}}^{\lambda \mu }(y)\right]={\hat{J}}^{\mu \nu }-x^{\mu }{\hat{P}}^{\nu }+x^{\nu }{\hat{P}}^{\mu }=\hat{\mathsf{T}}(x){\hat{J}}^{\mu \nu }\hat{\mathsf{T}}(x{)}^{-1},$$

$$\hat{\rho }=\frac{1}{Z}\exp \left[-{\beta }_{\mu }(x){\hat{P}}^{\mu }+\frac{1}{2}{\varpi }_{\mu \nu }{\hat{J}}_x^{\mu \nu }+\zeta \hat{Q}\right],$$

$$⟨\hat{O}(x)⟩=\frac{1}{Z}\mathrm{tr}\left[\exp \left(-\beta (x)\cdot \hat{P}+\frac{1}{2}\varpi :{\hat{J}}_x+\zeta \hat{Q}\right)\hat{O}(x)\right].$$

$$⟨\hat{O}(x){⟩}_{\beta (x)}=\frac{\mathrm{tr}\left[\exp \left(-\beta (x)\cdot \hat{P}+\zeta \hat{Q}\right)\hat{O}(x)\right]}{\mathrm{tr}\left[\exp \left(-\beta (x)\cdot \hat{P}+\zeta \hat{Q}\right)\right]}.$$

$${\mathrm{e}}^{\hat{A}+\hat{B}}={\mathrm{e}}^{\hat{A}}\left[1+\sum _{n=1}^{\infty }{\int }_0^1\mathrm{d}{\lambda }_1{\int }_0^{{\lambda }_1}\mathrm{d}{\lambda }_2\cdots {\int }_0^{{\lambda }_{n-1}}\mathrm{d}{\lambda }_n\hat{B}({\lambda }_1)\hat{B}({\lambda }_2)\cdots \hat{B}({\lambda }_n)\right],$$

$$\hat{B}(\lambda )={\mathrm{e}}^{-\lambda \hat{A}}\hat{B}{\mathrm{e}}^{\lambda \hat{A}}.$$

$${\mathrm{e}}^{-\beta (x)\cdot \hat{P}+\zeta \hat{Q}+\frac{1}{2}\varpi :\hat{J}}={\mathrm{e}}^{-\beta (x)\cdot \hat{P}+\zeta \hat{Q}}\left[1+\sum _{n=1}^{\infty }\frac{{\varpi }^n}{2^n}{\int }_0^1\mathrm{d}{\lambda }_1{\int }_0^{{\lambda }_1}\mathrm{d}{\lambda }_2\cdots {\int }_0^{{\lambda }_{n-1}}\mathrm{d}{\lambda }_n{\hat{J}}_{x-\mathrm{i}\beta {\lambda }_1}{\hat{J}}_{x-\mathrm{i}\beta {\lambda }_2}\cdots {\hat{J}}_{x-\mathrm{i}\beta {\lambda }_n}\right],$$

$${\hat{J}}_{x-\mathrm{i}\lambda \beta }^{\mu \nu }={\mathrm{e}}^{\lambda (\beta \cdot \hat{P}-\zeta \hat{Q})}{\hat{J}}_x^{\mu \nu }{\mathrm{e}}^{-\lambda (\beta \cdot \hat{P}-\zeta \hat{Q})}={\mathrm{e}}^{\lambda (\beta \cdot \hat{P}-\zeta \hat{Q})}{\hat{J}}_x^{\mu \nu }{\mathrm{e}}^{-\lambda (\beta \cdot \hat{P}-\zeta \hat{Q})}=\hat{\mathsf{T}}(x-\mathrm{i}\lambda \beta ){\hat{J}}^{\mu \nu }\hat{\mathsf{T}}(x-\mathrm{i}\lambda \beta {)}^{-1}.$$

$${\rm T}_{\lambda}\big{(}\widehat{O}_{1}(\lambda_{1})\widehat{O}_{2}(\lambda_{2})\cdots\widehat{O}_{N}(\lambda_{N})\big{)}\equiv\widehat{O}_{p_{1}}(\lambda_{p_{1}})\widehat{O}_{p_{2}}(\lambda_{p_{2}})\cdots\widehat{O}_{p_{N}}(\lambda_{p_{N}})$$

$$p:\{1,2,\dots ,N\}\to \{1,2,\dots ,N\}$$

$${\lambda }_{p_1}\le {\lambda }_{p_2}\le \cdots \le {\lambda }_{p_N}.$$

$$\begin{array}{rl}{\mathrm{e}}^{-\beta (x)\cdot \hat{P}+\zeta \hat{Q}+\frac{1}{2}\varpi :\hat{J}(x)}=& {\mathrm{e}}^{-\beta (x)\cdot \hat{P}+\zeta \hat{Q}}\left[1+\sum _{n=1}^{\infty }\frac{{\varpi }^n}{2^{n}n!}{\int }_0^1\mathrm{d}{\lambda }_1\mathrm{d}{\lambda }_2\cdots \mathrm{d}{\lambda }_n{\mathrm{T}}_{\lambda }\left({\hat{J}}_{x-\mathrm{i}{\lambda }_1\beta }{\hat{J}}_{x-\mathrm{i}{\lambda }_2\beta }\cdots {\hat{J}}_{x-\mathrm{i}{\lambda }_n\beta }\right)\right]\\ =& {\mathrm{e}}^{-\beta (x)\cdot \hat{P}+\zeta \hat{Q}}{\mathrm{T}}_{\lambda }\left[\exp \left(\frac{{\varpi }_{\mu \nu }}{2}{\int }_0^1\mathrm{d}\lambda {\hat{J}}_{x-\mathrm{i}\lambda \beta }^{\mu \nu }\right)\right],\end{array}$$

$$⟨\hat{O}(x)⟩=\sum _{N=0}^{\infty }\frac{{\varpi }^N}{2^{N}N!}{\int }_0^1\mathrm{d}{\lambda }_1\mathrm{d}{\lambda }_2\cdots \mathrm{d}{\lambda }_N⟨{\mathrm{T}}_{\lambda }\left({\hat{J}}_{x-\mathrm{i}{\lambda }_1\beta }{\hat{J}}_{x-\mathrm{i}{\lambda }_2\beta }\cdots {\hat{J}}_{x-\mathrm{i}{\lambda }_n\beta }\hat{O}(x)\right){⟩}_{\beta (x),c},$$

$$\begin{array}{rl}⟨\hat{O}{⟩}_c& =⟨\hat{O}⟩\\ ⟨\hat{J}\hat{O}{⟩}_c& =⟨\hat{J}\hat{O}⟩-⟨\hat{J}⟩⟨\hat{O}⟩\\ ⟨{\hat{J}}_1{\hat{J}}_2\hat{O}{⟩}_c& =⟨{\hat{J}}_1{\hat{J}}_2\hat{O}⟩-⟨{\hat{J}}_1⟩⟨{\hat{J}}_2\hat{O}{⟩}_c-⟨{\hat{J}}_2⟩⟨{\hat{J}}_1\hat{O}{⟩}_c-⟨\hat{O}⟩⟨{\hat{J}}_1{\hat{J}}_2{⟩}_c-⟨J_1⟩⟨J_2⟩⟨\hat{O}⟩\\ & ⋮\end{array}$$

$$\begin{array}{rl}⟨\hat{O}(x)⟩=& ⟨\hat{O}(x){⟩}_{\beta (x)}+\frac{{\varpi }_{\mu \nu }}{2}{\int }_0^1\mathrm{d}\lambda ⟨{\mathrm{T}}_{\lambda }\left({\hat{J}}_{x-\mathrm{i}\lambda \beta }^{\mu \nu }\hat{O}(x)\right){⟩}_{\beta (x),c}\\ & +\frac{{\varpi }_{\mu \nu }{\varpi }_{\rho \sigma }}{8}{\int }_0^1\mathrm{d}{\lambda }_1\mathrm{d}{\lambda }_2⟨{\mathrm{T}}_{\lambda }\left({\hat{J}}_{x-\mathrm{i}{\lambda }_1\beta }^{\mu \nu }{\hat{J}}_{x-\mathrm{i}{\lambda }_2\beta }^{\rho \sigma }\hat{O}(x)\right){⟩}_{\beta (x),c}+\mathcal{O}({\varpi }^3).\end{array}$$

$$\begin{array}{rl}⟨\hat{O}(x)⟩=& ⟨\hat{O}(x){⟩}_{\beta (x)}+\frac{{\varpi }_{\mu \nu }}{2|\beta |}{\int }_0^{|\beta |}\mathrm{d}\tau ⟨{\mathrm{T}}_{\tau }\left({\hat{J}}_{x-\mathrm{i}\tau u}^{\mu \nu }\hat{O}(x)\right){⟩}_{\beta (x),c}\\ & +\frac{{\varpi }_{\mu \nu }{\varpi }_{\rho \sigma }}{8|\beta {|}^2}{\int }_0^{|\beta |}\mathrm{d}{\tau }_1\mathrm{d}{\tau }_2⟨{\mathrm{T}}_{\tau }\left({\hat{J}}_{x-\mathrm{i}{\tau }_{1}u}^{\mu \nu }{\hat{J}}_{x-\mathrm{i}{\tau }_{2}u}^{\rho \sigma }\hat{O}(x)\right){⟩}_{\beta (x),c}+\mathcal{O}({\varpi }^3).\end{array}$$

$$⟨{\hat{O}}_1(x_1){\hat{O}}_1(x_2)\cdots {\hat{O}}_n(x_n){⟩}_{\beta (x)}=⟨{\hat{O}}_1(x_1-x_n){\hat{O}}_2(x_2-x_n)\cdots {\hat{O}}_n(0){⟩}_{\beta (x)},$$

$$\begin{array}{rl}⟨\hat{O}(x)⟩=& ⟨\hat{O}(0){⟩}_{\beta (x)}+\frac{{\varpi }_{\mu \nu }}{2|\beta |}{\int }_0^{|\beta |}\mathrm{d}\tau ⟨{\mathrm{T}}_{\tau }\left({\hat{J}}_{-\mathrm{i}\tau u}^{\mu \nu }\hat{O}(0)\right){⟩}_{\beta (x),c}\\ & +\frac{{\varpi }_{\mu \nu }{\varpi }_{\rho \sigma }}{8|\beta {|}^2}{\int }_0^{|\beta |}\mathrm{d}{\tau }_1\mathrm{d}{\tau }_2⟨{\mathrm{T}}_{\tau }\left({\hat{J}}_{-\mathrm{i}{\tau }_{1}u}^{\mu \nu }{\hat{J}}_{-\mathrm{i}{\tau }_{2}u}^{\rho \sigma }\hat{O}(0)\right){⟩}_{\beta (x),c}+\mathcal{O}({\varpi }^3).\end{array}$$

$${\varpi }_{\mu \nu }={ϵ}_{\mu \nu \rho \sigma }{w}^{\rho }{u}^{\sigma }+{\alpha }_{\mu }{u}_{\nu }-{\alpha }_{\nu }{u}_{\mu }$$

Full text

I Introduction

Relativistic hydrodynamics is an effective dynamical theory of systems at local thermodynamic equilibrium. The condition of local thermodynamic equilibrium, makes it possible to describe the dynamics of interacting quantum fields with classical partial differential equations involving few thermodynamic fields (temperature, velocity, chemical potentials) provided that there is a separation between macroscopic and microscopic scales: the thermodynamic fields should vary significantly over distances which are much (or sufficiently) larger than microscopic scales of the quantum theory, such as mean free path, Compton wavelength etc. The branches of physics for which relativistic hydrodynamics is an effective tool include astrophysics, cosmology as well relativistic heavy ion collisions Baier:2006um ; Huovinen:2006jp ; Florkowski ; Gale:2013da ; Heinz:2013th ; kodama . In the past few years, the derivation of hydrodynamic equations has drawn much attention: from kinetic theory York:2008rr ; Muronga:2006zx ; Tsumura:2006hn ; Betz:2008me ; Monnai:2010qp ; Van:2011yn ; denicol ; jaiswal , from fluid/gravity correspondence baier ; Bhattacharyya:2008jc ; Natsuume:2007ty ; Hubeny:2011hd , from the phenomenological extension of the non equilibrium thermodynamics Koide:2006ef ; Fukuma:2011pr , from the projection operator method Koide:2008nw ; Minami:2012hs , from non equilibrium statistical operator method Hayata:2015 and from imposing some symmetries on the system Banerjee:2012iz ; kovtun . The success of relativistic hydrodynamics in heavy ion collisions lately arose fundamental questions about its domain of applicability Romatschke:2016hle .

Relativistic hydrodynamics is based on the expansion of the mean stress-energy tensor as a function of the thermodynamic fields and their gradients. The zero-order term is the well known ideal form of the mean stress-energy tensor:

[TABLE]

which is an exact expression at the global thermodynamic equilibrium with constant parameters $T,\mu$ and $u$. The first-order corrections to the (1) in the gradients of $u$, $T$ ad $\mu$ are dissipative, that is they imply the irreversible increase of entropy. The second-order corrections - including quadratic terms in first-order derivatives as well as second order derivatives - to the ideal form were classified in refs. baier ; Romatschke:2009kr ; Moore:2010bu and discussed in refs. Bhattacharyya:2008jc ; York:2008rr ; moore ; kovtun ; Becattini:2015nva ; Panerai .

Among the second order terms in the gradients, there are quadratic terms in acceleration and vorticity which do not contribute to the entropy increase, i.e. they are non-dissipative. Such terms were dubbed in ref. moore as thermodynamic, and it was shown in ref. Becattini:2015nva that they appear in the most general form of thermodynamic equilibrium in flat spacetime involving non-vanishing acceleration and vorticity, such as the rotating fluid with constant angular velocity landau ; vilenkin . Furthermore, in ref. Becattini:2015nva was shown that, at least for free fields, these terms are of quantum origin, namely they vanish in the $\hbar\to 0$ limit, as it was already argued in ref. Romatschke:2009kr . Their quantum origin is also born out by the fact that they are missing in the stress-energy tensor expression from the classical expansion of the Boltzmann distribution function in the kinetic theory jaiswal . The interest for such terms has been certainly reinforced by the recent observation of a thermal vorticity of the order of some percent by the STAR experiment through the measurement of hyperon polarization star .

The scalar coefficient multiplying the term $\omega^{\mu}\omega^{\nu}$, where $\omega$ is the kinematic vorticity, was calculated in ref. moore with a suitable ``Kubo formula" involving stress-energy tensor three-point functions. In ref. Becattini:2015nva it was shown that, on the other hand, those coefficients can be expressed as correlators of conserved generators of the Poincaré group and the stress-energy tensor at the usual thermodynamic equilibrium with constant $T$ and $\mu$. This property involves a remarkable simplification of the "Kubo formulae", insofar as the constancy of generators allows to remove time integration, unlike in those of dissipative coefficients, like shear viscosity. The ultimate reason thereof is the relation of these terms to global equilibrium configurations. In fact, in this paper we will use a different method, where the time integration survives in imaginary time.

In ref. Becattini:2015nva their explicit expression was found for the real scalar free field. In this work, we extend that calculation to the free charged scalar field and the Dirac field including a finite chemical potential. We will be using an operator formalism throughout instead of the functional approach kovtun ; Bhattacharyya:2012nq ; Banerjee:2012iz . As it will be shown, the operator formalism is very convenient for a simple and systematic derivation of the "Kubo formulae" of the non-dissipative coefficients. On the other hand, the functional approach allows to obtain relations between those coefficients kovtun ; Bhattacharyya:2012nq ; Banerjee:2012iz , which are more difficult to extract in the operator formalims.

Notation

In this paper we use the natural units, with $\hbar=c=k_{B}=1$.

The Minkowskian metric tensor is ${\rm diag}(1,-1,-1,-1)$; for the Levi-Civita symbol we use the convention $\epsilon^{0123}=1$.

We will use the relativistic notation with repeated indices assumed to be summed over, however contractions of indices will be sometimes denoted with dots, e.g. $u\cdot T\cdot u\equiv u_{\mu}T^{\mu\nu}u_{\nu}$. Operators in Hilbert space will be denoted by a large upper hat, e.g. ${\widehat{T}}$ (with the exception of Dirac field operator that is denoted by ${\sf{\Psi}}$) while unit vectors with a small upper hat, e.g. $\hat{v}$. The stress-energy tensor is assumed to be symmetric with an associated vanishing spin tensor.

II Global thermodynamic equilibrium density operator

The general covariant form of the local thermodynamic equilibrium density operator was introduced in zubarev ; weert ; weldon and lately discussed in detail in refs. becacov ; becalocal ; Becattini:2015nva ; Hayata:2015 ; Hongo:2016mqm :

[TABLE]

where $\widehat{T}^{\mu\nu}$ and $\widehat{j}^{\mu}$ are the stress-energy tensor and a conserved current operators, $\beta^{\mu}$ is the four-temperature vector such that

[TABLE]

is the proper comoving temperature and $\zeta$ is the ratio of comoving chemical potential and temperature $\zeta=\mu/T$; $\Sigma$ is a spacelike 3-D hypersurface. This operator is obtained by maximizing the entropy with fixed mean energy, momentum and charge density becalocal . If time-like, the four-temperature vector $\beta$ defines a local four-velocity of the fluid (the $\beta$* frame* becalocal ; Van:2013sma ):

[TABLE]

and the magnitude of $\beta$ is the proper temperature measured by a thermometer moving with four-velocity $u$.

If the four-temperature $\beta$ is a Killing vector becacov :

[TABLE]

and $\zeta$ is constant, the density operator (2) is a proper global thermodynamic equilibrium density operator insofar as it becomes independent of the hypersurface $\Sigma$ with suitable conditions at a timelike boundary. The general solution of the eq. (4) in Minkowsky space-time is known and can be expressed in terms of a constant four-vector $b_{\mu}$ and a constant antisymmetric tensor $\varpi_{\mu\nu}$:

[TABLE]

The antisymmetric tensor $\varpi_{\mu\nu}$ is called thermal vorticity; from the above equation it turns out that it is the antisymmetric part of the gradient of $\beta$:

[TABLE]

Once the eq. (5) is plugged into the eq. (2), the general expression of the equilibrium density operator in flat space-time is obtained Becattini:2015nva .

[TABLE]

where $\widehat{P}$ is the four momentum operator, $\widehat{Q}$ the conserved charge and $\widehat{J}^{\mu\nu}$ are the generators of the Lorentz transformations:

[TABLE]

This form shows that the equilibrium state in Minkowski spacetime is described by the 10 constant parameters comprised by $b$ and $\varpi$, that is as many as the generators of its maximal symmetry group, the Poincaré group.

The familiar form of the equilibrium density operator can be recovered as a special case of eq. (7), with $b$ timelike and $\varpi=0$ in eq. (5):

[TABLE]

which, in the rest frame of $b$, reduces to well known grand-canonical ensamble density operator. As the above density operator is invariant by translation, we denote this familiar kind of equilibrium as homogeneous thermodynamic equilibrium. Another form of global equilibrium, which is a special case of eq. (7) is the pure rotation, with $b=(1/T_{0},{\bf 0})$ and $\varpi\propto\omega_{0}/T_{0}$ such that:

[TABLE]

which has recently raised much attention for fermions ambrus ; chernodub .

In order to represent a physical fluid at equilibrium, $\beta$ must be a timelike vector. However, the equilibrium $\beta$ vector field in eq. (5) is timelike everywhere only if $\varpi=0$, corresponding to the homogeneous equilibrium with constant temperature, as we have seen. Indeed, if $\varpi\neq 0$ there are spacetime regions where $\beta$ is spacelike or lightlike, like in the rotating case (9). Nevertheless, this does not undermine the possibility to calculate the mean value of local operators in the regions where $\beta$ is timelike. The idea is to expand the exponent in the density operator (7) from the same point $x$ of a local operator $\widehat{O}(x)$, whose mean value is to be calculated, in powers of the supposedly small thermal vorticity $\varpi$. Thereby, the zero-order term is the same mean value as at homogeneous global equilibrium with four-temperature equal to the $\beta$ field in $x$ and corrections arise form a power series in $\varpi$. This method was presented and applied in ref. Becattini:2015nva for local operators of the free scalar field and it will be outlined in the next section.

III Mean value of a local operator

The mean value of a local operator $\widehat{O}(x)$ is defined as:

[TABLE]

where the subscript indicates the need of a renormalization procedure; for free fields, renormalization involves the subtraction of the vacuum expectation value, or, tantamount the use of normal ordered operators. At global equilibrium, the density operator is given by eq. (7) and we can take advantage of this transformation rule for the angular momentum operator:

[TABLE]

where $\widehat{\sf T}(x)=\exp[\mathrm{i}x\cdot\widehat{P}]$ is the translation operator, to rewrite it as:

[TABLE]

where $x$ is the point where the operator $\widehat{O}(x)$ is to be reckoned and the eq. (5) has been used. Thereby, the mean value (10) turns into:

[TABLE]

If $\varpi\ll 1$ (it is adimensional) one can expand the mean value (13) in powers of $\varpi$, the leading order being its homogeneous equilibrium value (i.e. the grand-canonical ensemble average) with a constant four-temperature $\beta$ equal to the four-temperature field in the point $x$:

[TABLE]

In general, the coefficients of the expansion of (13) in $\varpi$ can be expressed in terms of correlators at the homogeneous thermodynamic equilibrium. In this work, we will use the well-known expansion formula of the exponential of a sum of two non commuting operators $\widehat{A}+\widehat{B}$:

[TABLE]

where $B(\lambda_{1})$ is defined as:

[TABLE]

For the mean value in eq. (13) $\widehat{A}=-\beta(x)\cdot\widehat{P}+\zeta\widehat{Q}$ and the small term is $\widehat{B}=\frac{1}{2}\varpi:\widehat{J}_{x}$. Therefore, by using the expansion (15), we can write:

[TABLE]

where $\varpi^{n}$ must be understood as a product of $n$ tensors $\varpi^{\mu\nu}$ with indices fully contracted with the coupled angular momentum operators and $\widehat{J}_{x-\mathrm{i}\lambda\beta}$ stands for:

[TABLE]

This expansion can be written in a more compact form introducing the ${\rm T}_{\lambda}$-ordered product, where the time path is defined on the imaginary direction $\mathrm{i}\beta$:

[TABLE]

with $p$ the permutation that orders $\lambda$ by value:

[TABLE]

Using the above definition and changing the integration limits accordingly, the equation (16) becomes:

[TABLE]

where - as usual - the T-ordered exponential is a shorthand notation of the Taylor expansion. The mean value of an operator (13) can now be calculated as a power series of $\varpi$ and a cumulant expansion is obtained, namely:

[TABLE]

where the correlators are defined respect to the density operator in (14), $\langle\cdots\rangle_{\beta(x)}=Z^{-1}{\rm tr}[e^{-\beta(x)\cdot\widehat{P}+\zeta\widehat{Q}}\cdots]$ and the subscript $c$ means that only the connected correlators are involved, namely:

[TABLE]

So, for instance, at the second order in $\varpi$:

[TABLE]

This expression can be further worked out by using the local four velocity $u^{\mu}$ in eq. (3), the inverse proper temperature $|\beta(x)|=1/T=\sqrt{\beta^{2}}$ and changing the integration variable to $\tau=|\beta|\lambda$:

[TABLE]

Since the density operator is the homogeneous one, defined in (14), which is invariant by translation, one can write:

[TABLE]

whence in the (17):

[TABLE]

This expression makes it apparent that the $x$-dependence of the various terms is determined only by the four-temperature vector.

IV Acceleration and vorticity decomposition

To proceed, it is useful to decompose the thermal vorticity $\varpi$ by using the four-velocity $u$ in the same fashion as the electromagnetic field tensor is decomposed into comoving electric and magnetic field. Specifically:

[TABLE]

where $\alpha$ and $w$ are two spacelike vectors defined as:

[TABLE]

The physical meaning of $\alpha$ and $w$, at the global equilibrium, is that of the acceleration field and vorticity vector field divided by the proper temperature. To prove this, one has just to show, by double contracting eq. (4) with $\beta$, and using the eq. (3), that the proper temperature does not change along the flow:

[TABLE]

so that:

[TABLE]

Similarly, for $w$:

[TABLE]

being $\omega_{\mu}=-\frac{1}{2}\sqrt{\beta^{2}}\epsilon_{\mu\nu\rho\sigma}\partial^{\nu}u^{\rho}u^{\sigma}$ the local vorticity vector. It is also useful to define a fourth four-vector:

[TABLE]

where $\Delta_{\mu\nu}=g_{\mu\nu}-u_{\mu}u_{\nu}$ is the transverse projector to the four-velocity. The four-vector $\gamma$ is non-vanishing if $\alpha$ and $w$ are linearly independent and is, by construction, orthogonal to the other three four-vectors. More relations involving the derivatives of these vectors can be found in the Appendix B.

In the local rest frame, the above four-vectors can be expressed as:

[TABLE]

where $\bm{a}$ is the acceleration seen by the local rest frame and $\bm{\omega}$ the angular velocity. Hence, restoring the physical constants:

[TABLE]

which shows the quantum nature of the adimensional parameters $\alpha,w,\gamma$.

The generators of the Lorentz group $\widehat{J}$ can be similarly decomposed into local boosts $\widehat{K}$ and local angular momenta $\widehat{J}$ by using the four velocity $u_{\mu}$:

[TABLE]

where the operators $\widehat{K}$ and $\widehat{J}$ are defined as:

[TABLE]

Therefore, by using the above decomposition and the eq. (19) one can write:

[TABLE]

where $\{\cdots,\cdots\}$ is the anti-commutator. Plugging the eqs. (23) into the (18):

[TABLE]

Now, advantage can be taken of the transformation properties under rotation, reflection and time reversal of the boost and angular momentum operators to classify all corrections to the mean value of any operator at thermodynamic equilibrium.

V The stress-energy tensor

The mean value of the energy momentum tensor $\widehat{T}^{\mu\nu}$ receives contributions only from second order terms in acceleration and rotation in eq. (LABEL:tauexp3) because the first order terms vanish due to reflection and time reversal symmetry. Thus, the expansion (LABEL:tauexp3) becomes:

[TABLE]

where, according to eq. (14), the mean value of the stress-energy tensor at homogeneous thermodynamic equilibrium $\langle\widehat{T}^{\mu\nu}(x)\rangle_{\beta(x)}$ coincides with its ideal form (1):

[TABLE]

Note that in the above equation we have spelled out all $x$ dependencies and that $\rho,p$ are the thermodynamic equilibrium functions of temperature and chemical potential usually obtained in thermodynamics.

One can now decompose the above correlators into irreducible tensors under rotation. By taking advantage of the rotational invariance of the homogeneous equilibrium density operator in eq. (14), the number of actual coefficients in eq. (25) can be reduced, and it can be shown Becattini:2015nva that:

[TABLE]

with $\gamma$ as in eq. (21). Thereby, the general expression of the stress-energy tensor at the second order in thermal vorticity reads Becattini:2015nva :

[TABLE]

The coefficients in eq. (26) are Lorentz scalars depending only on the magnitude of the local inverse four temperature $|\beta|$, that is the proper temperature $T$, and the proper chemical potential $\mu$, so they can be calculated in any frame. The rest frame is the most convenient choice because $\beta^{\mu}=(1/T,{\bm{0}})$ and the homogeneous equilibrium density operator takes on the familiar grand-canonical form. We denote the mean values in rest frame with a subscript $T$, that is:

[TABLE]

Hence, the coefficients in (26) are found as specific combinations Becattini:2015nva of thermal connected correlators:

[TABLE]

The correlation functions in (27) can be expressed as Euclidean three-point functions of the stress-energy tensor and can be calculated with the imaginary time formalism. The shifted boost and angular momentum generators can be written as (see eq. (11)):

[TABLE]

So one can expand the integral expressions of the generators of the Lorentz group and write the basic structure appearing in all coefficients in (27) as:

[TABLE]

where

[TABLE]

is the imaginary time evolved stress-energy tensor operator. By using the eq. (28), the eqs. (27) can be rewritten in terms of these auxiliary quantities:

[TABLE]

The above coefficients are not all independent, in fact there are relations between them stemming from the conservation equation:

[TABLE]

which states that the mean value of $\widehat{T}^{\mu\nu}$ is conserved if the corresponding operator is conserved and the density operator $\widehat{\rho}$ is time-independent. The relation between the above second-order coefficients can be obtained in a functional approach Bhattacharyya:2012nq ; kovtun ; Banerjee:2012iz by taking advntage of the invariance of the generating functional by diffeomorfisms; in the present operatorial approach, one has to enforce the continuity equation of the stress-energy tensor, which leads to these equalities (see Appendix B):

[TABLE]

where all derivatives are to be taken with fixed $\zeta=\mu|\beta|$. Thus, only $D_{\alpha},D_{w},A,W$ are actually independent, while $G,U_{\alpha},U_{w}$ can be obtained directly from the eqs. (31). Nevertheless, in this work, we have calculated all the coefficients by using the eqs. (27) and have used the relations in eq. (31) as a consistency check.

V.1 Free complex scalar field

We now calculate the coefficients (27) for a complex scalar field at finite temperature and chemical potential, using the imaginary time formalism. In this case we have:

[TABLE]

where $\xi=0$ corresponds to the canonical stress-energy tensor and $\xi=1/6$ to the improved stress-energy tensor callan .

For finite chemical potential, it is convenient to switch to the Euclidean time formalism with the modified hamiltonian $\widehat{H}-\mu\widehat{Q}$ and write the field evolution as:

[TABLE]

Indeed, as the stress-energy tensor commutes with the charge operator, the formulae obtained in the previous section are unaffected. Defining:

[TABLE]

and:

[TABLE]

where $\omega_{n}=2\pi n/|\beta|$ are the Matsubara frequencies, the propagator in the imaginary time reads Kapusta:2006pm ; Laine:2016hma :

[TABLE]

where:

[TABLE]

and ${P^{\pm}}^{2}$ is the Euclidean squared magnitude of the four vector $P^{\pm}$, implicitly defined in eqs. (34). The eqs. (34) are the building block to evaluate the three-point function in eq. (28) because:

[TABLE]

where $\widehat{T}^{\mu\nu}(X)$ is given by (29). The (36) can be computed using the point splitting procedure; first, consider:

[TABLE]

where the form of the differential operators $\mathcal{D}^{\mu\nu}(\partial_{X_{1}},\partial_{X_{2}})$ can be inferred from the eq. (32):

[TABLE]

In the eq. (38), the scalar product has the Euclidean signature, that is $\partial_{X_{1}}\cdot\partial_{X_{2}}=\partial_{\tau_{1}}\partial_{\tau_{2}}+\sum_{i}\partial_{x_{i1}}\partial_{x_{i2}}$ and the D'Alembertian as well: $\square_{XE}=\partial_{\tau}^{2}+\sum_{i}\partial_{i}^{2}$. The imaginary unit in front of the differential operator $(-\mathrm{i})^{\delta_{0\alpha}+\delta_{0\beta}}$ is a consequence of the the Wick rotation. The evaluation of the three-point function can be done by employing the standard Wick theorem and since only its connected part appears in (37), only the following two terms survive:

[TABLE]

Inserting the Fourier decomposition of eq. (34), the differential operators in (37) give rise to a polynomial in momentum; thereby, the limits can be readily done and one obtains:

[TABLE]

where $\Delta$ is defined in (35) and the function $F_{1}$ and $F_{2}$ are just polynomials of momenta:

[TABLE]

Now we take the $Z\to 0$ limit in eq. (39), as prescribed by eq. (36) and, separating the integration from the sum over frequencies, we get:

[TABLE]

where:

[TABLE]

The functions $F_{1}$ and $F_{2}$ are polynomials, hence analytic and the sum over the frequencies can be carried out by using the formula Landsman:1986uw :

[TABLE]

where $\theta(\tau)$ is the Heaviside function and $-|\beta|<\tau<|\beta|$; $n_{B}$ is the Bose distribution function:

[TABLE]

The eq. (42) is needed to work out the function in eq. (41):

[TABLE]

where $E_{p}=\sqrt{{\bm{p}}^{2}+m^{2}}$, and we have defined

[TABLE]

and similarly for $\tilde{Q}(s),\,\tilde{K}(s)$. Note that, after the frequency summation, the arguments of the functions $F_{1}$ and $F_{2}$ no longer depend on the chemical potential, and that, thanks to the symmetry properties of the polynomials, the functions $F_{1}$ and $F_{2}$ become indeed the same. Thus, the eq. (41) becomes:

[TABLE]

Now we can take advantage of the formula:

[TABLE]

to integrate over the coordinates $x$ and $y$ in the eq. (36) and, by using the eq. (40) we obtain:

[TABLE]

One can now plug the (45) into the (46), integrate over $\tau_{1}$ and $\tau_{2}$ so as to express the $C^{\alpha\beta|\gamma\rho|\mu\nu|ij}$ as an integral over momentum of combinations of derivatives of the Bose distribution function. After setting the appropriate indices in (46) according to the (30), the second-order coefficients of the stress-energy tensor are finally obtained:

[TABLE]

where $n_{B}^{\prime\prime}$ denotes the second derivative of the Bose distribution function with respect to the energy.

It is worth noticing that the above expressions (47) fulfill the relations (31).

V.2 Free Dirac field

In this section we calculate the coefficients (27) for the free Dirac field at finite temperature and chemical potential. The thermodynamic properties of the Dirac field can be deduced from the Euclidean Lagrangian $\widehat{{\cal L}}_{E}=-\widehat{\mathcal{L}}(t=-\mathrm{i}\tau)$:

[TABLE]

It is convenient to write the (48) with the so-called Euclidean Dirac matrices:

[TABLE]

fulfilling the relations

[TABLE]

so that the eq. (48) can be written:

[TABLE]

The thermal propagator of the Dirac field reads 111Latin characters $a,b,\cdots$ denote spinorial indicesKapusta:2006pm ; Laine:2016hma :

[TABLE]

where $X,Y,P^{\pm}$ are defined in eq. (33) and $\Delta(P^{\pm})$ in eq. (35); the sum runs over the fermionic Matsubara frequencies $\omega_{n}=2\pi(n+\frac{1}{2})/|\beta|$; $\not{P}^{\pm}=\tilde{\gamma}_{\mu}P_{\mu}^{\pm}$ is the standard contraction between the (Euclidean) Dirac matrices $\tilde{\gamma}_{\mu}$ and the (Euclidean) four-momenta $P^{\pm}=(p_{n}\pm\mathrm{i}\mu,{\bm{p}})$. The Euclidean canonical stress-energy tensor (see eq. (29)) reads:

[TABLE]

where the $\mathrm{i}^{\delta_{\nu 0}}$ factor stems from Wick rotation. The Belinfante-symmetrized stress-energy tensor is the symmetric part of the canonical one:

[TABLE]

which can be expressed according to the point-splitting procedure as:

[TABLE]

where:

[TABLE]

The stress-energy tensor three-point correlation function (36) needed to extract the various coefficients is similar to that in (37):

[TABLE]

Like for the boson case, thanks to the Wick theorem and the presence of the connected part, only two contractions survive:

[TABLE]

Plugging the expression of the propagators (49), after some algebra and integration variable manipulation, we get:

[TABLE]

where the functions $G_{1}$ and $G_{2}$ are defined as:

[TABLE]

The trace is to be carried out over spinorial indices by using the Euclidean $\gamma$ matrices properties:

[TABLE]

We finally obtain:

[TABLE]

where ${\cal S}$ stands for a full symmetrization of the subscript indices (without factorials). This expression is very similar to that for the boson field obtained in the previous subsection, with the proviso that now the Matsubara frequencies to be summed involve odd integers. We can then trace the bosonic procedure, first setting $Z$ to zero and defining the auxiliary function $S^{F}$:

[TABLE]

with:

[TABLE]

The sums over fermionic frequencies can be done by using a formula corresponding to (42) for the boson field:

[TABLE]

$n_{F}$ being the Fermi-Dirac distribution function:

[TABLE]

Like for the boson case, $S^{F}$ comprises 8 terms:

[TABLE]

with $\tilde{P}(s_{1}),\,\tilde{Q}(s_{2})$ and $\tilde{K}(s_{3})$ defined in (44), and the polynomial $G$ no longer depends on chemical potential after frequency summation. Furthermore,

[TABLE]

and so we can rewrite $S^{F}$ in eq. (52):

[TABLE]

The spatial integrations are straightforward, and so we can finally write down the general three-point function for the free Dirac field:

[TABLE]

with $S^{F}$ given by the (53). The coefficients of the symmetrized stress-energy tensor can now be expressed by using the relations (30) like in the boson case, as integrals of the second derivative of the Fermi-Dirac distribution function weighted with polynomials of momentum and mass:

[TABLE]

where $n_{F}^{\prime\prime}(E_{p}\pm\mu)={\rm d}^{2}n_{F}(E_{p}\pm\mu)/{\rm d}E_{p}^{2}$.

Also in the fermionic case, it can be shown that the coefficients (54) fulfill the relations (31).

VI The vector current

Also conserved currents can receive corrections in general thermodynamic equilibrium with non-vanishing thermal vorticity (12). For a vector current $\widehat{j}^{\mu}(x)$, the expansion (LABEL:tauexp3) yields:

[TABLE]

The leading order term $\langle\widehat{j}^{\mu}(x)\rangle_{\beta(x)}$ is simply the homogeneous equilibrium current $nu^{\mu}$ where $n$ is the proper charge density. The first order corrections in $\alpha$ and $w$ are zero due to the time-reversal and parity symmetries, just like for the stress-energy tensor; hence, again, the first non-vanishing corrections are quadratic in thermal vorticity. The invariance under rotation selects only three allowed tensor combinations with $\alpha$ and $w$:

[TABLE]

where $n$ is mean value of the charge density at the homogeneous equilibrium. By comparing the eqs. (56) and (55) and taking into account the rotational invariance, we can identify in the local rest frame the following formulae for $N_{\alpha}$, $N_{w}$ and $G_{V}$:

[TABLE]

The right hand sides of the above equalities involve the three-point thermal functions of the stress-energy tensor (twice) and the current operator. Defining:

[TABLE]

the coefficients in eq. (57) can be obtained as linear combinations of (58):

[TABLE]

VI.1 Free complex scalar field

The current of the free scalar field reads:

[TABLE]

and its Euclidean counterpart:

[TABLE]

The general coefficient in eq. (58) can be calculated by using the point splitting procedure:

[TABLE]

where the differential operator $\mathcal{D}^{\alpha\beta}$ is defined in (38). The operator $\mathcal{J}^{\mu}$ is found inserting the current operator into the correlation function:

[TABLE]

Thus:

[TABLE]

where the function $S$ is now:

[TABLE]

and the polynomials in the momenta are:

[TABLE]

Since the polynomial $\mathcal{J}^{\mu}(X,Y)$ is antisymmetric by argument swap, after the frequency summation the two polynomials $H_{1}$ and $H_{2}$ no longer depend on the chemical potential and are opposite, i.e. $H_{1}=-H_{2}$.

Then, summing over Matsubara frequencies and reminding definition (44):

[TABLE]

Taking the derivative of $S$ respect to $\bm{q}$ and $\bm{k}$, integrating over $\tau_{1}$ and $\tau_{2}$, according to eq. (61), the coefficients (59) turn out to be:

[TABLE]

VI.2 Free Dirac field

The vector current of Dirac field $\widehat{j}^{\mu}={\sf{\bar{\Psi}}}\gamma^{\mu}{\sf{\Psi}}$ in its Euclidean version reads:

[TABLE]

where $\tilde{\gamma}$ are the Euclidean gamma matrices. The generic coefficient (58) in this case can be written as:

[TABLE]

where the matrix associated to the stress-eenrgy tensor $\mathcal{D}^{\alpha\beta}(\partial_{X_{1}},\partial_{X_{2}})_{ab}$ is defined in eq. (50), and $\mathcal{J}^{\mu}(\partial_{Z_{1}},\partial_{Z_{2}})_{ef}$ corresponding to the vector current is:

[TABLE]

The function $S$ analogous to that in eq. (62) is:

[TABLE]

where, in this case, the vertex functions corresponding to the stress-energy tensor and the vector current insertion into the correlation function (64) are denoted as $M_{1}$ and $M_{2}$:

[TABLE]

The sum over frequencies yields:

[TABLE]

where, like in the boson case, the two functions $M_{1}$ and $M_{2}$ turn out to be opposite:

[TABLE]

After performing the integrations in $\tau$, taking the derivative with respect to the momenta $\bm{q}$ and $\bm{k}$ and setting the appropriate indices according the eq. (59), the coefficients for the Dirac field are finally obtained:

[TABLE]

VII The axial current

It is also worth studying the lowest-order expression of the mean axial current $\langle\widehat{j}_{A}^{\mu}\rangle$ at thermodynamic equilibrium with vorticity and acceleration. For the free field without interactions, the axial current $\widehat{j}_{A}^{\mu}={\sf{\bar{\Psi}}}\gamma^{\mu}\gamma^{5}{\sf{\Psi}}$ is conserved only in the massless limit, being, as it is well known:

[TABLE]

Unlike the vector current, because of its properties under space reflection and time-reversal, its mean value vanishes at homogeneous thermodynamic equilibrium, but it has a first-order correction proportional to vorticity:

[TABLE]

The coefficient $W^{A}$ can be calculated from the two-point function between the angular momentum operator and the axial current operator; in the rest frame its formula is:

[TABLE]

To calculate the latter two terms, as usual, we write the correlation functions in the Euclidean form. The Euclidean axial current reads:

[TABLE]

where $\tilde{\gamma}^{5}$ is defined as:

[TABLE]

and it is equal to the usual $\gamma^{5}$ matrix. Then, we can write:

[TABLE]

By using the Wick theorem and the momentum representation of the Dirac propagator we can write:

[TABLE]

where the function $A$ results from of the composite operator in the correlation functions, in this case:

[TABLE]

After having defined $S_{\textrm{axial}}$ in the same fashion as in previous sections:

[TABLE]

we get:

[TABLE]

In the above expression, the two sums over frequencies are independent; using (51):

[TABLE]

with the function $A$ defined in (68). Recalling the definition of $W^{A}$ in the eq. (67) the integration over the spatial coordinates can be done as a derivative respect to a loop momenta:

[TABLE]

Choosing the suitable indices as in eq. (67), taking the derivative with respect to the momentum $\bm{q}$ and integrating over both $\tau$ and the angles, we finally obtain the coefficient $W^{A}$ for the free Dirac field:

[TABLE]

This result can be checked by using the conservation of the axial current in the massless case. In general, the divergence of the mean axial current, at first order in $w$:

[TABLE]

where we have used the fact that $W^{A}$ depends only on the magnitude of $\beta$ and the expressions of the gradients of $w^{\mu}$ that can be found in Appendix B. Specifically, the above combination gives rise to:

[TABLE]

which is manifestly vanishing for a massless Dirac field.

VII.1 Discussion: axial current and anomalies

The equation (66) states that in rotating fermion gas there is a non-vanishing axial current and consequently the right and left chiral fermions get separated. This current can be pictorially understood with a simple argument, which is strictly valid only for massless particles. In a rotating system the fermions spin tend to align with the direction of rotation $\bm{\omega}$ independently of their charge Becattini:2007nd . The right handed particles have their momentum aligned with the spin, consequently will move in the direction of the spin, i.e. we get a net right handed particles flow in the direction of $\bm{\omega}$. On the other hand the left handed particles will move in the opposite direction, giving a net left handed particles flow opposite to $\bm{\omega}$. Together these flows give an axial current $\bm{j}_{\textup{A}}=n_{R}\bm{v}_{R}-n_{L}\bm{v}_{L}\propto\bm{\omega}$. The very reason of the non-vanishing axial current is simply that it is allowed by the symmetry of the statistical operator, as $\bm{\widehat{j}}_{A}$ has the same properties of the angular momentum operator $\widehat{\bf J}$ under reflection, time reversal and charge conjugation; so:

[TABLE]

Using the method described in the next section VIII one can obtain the value of the coefficient $W^{\textup{A}}$ (69) for $m=0$:

[TABLE]

This result is precisely the same found in recent calculations of the mean axial current related to the chiral vortical effect (CVE), which is the onset of a vector current along the vorticity due to the axial anomaly (see Kharzeev:2015znc and references therein). It should be pointed out, though, that the onset of an axial current along vorticity is conceptually a distinct phenomenon, so we denote it by axial vortical effect (AVE) following ref. Kalaydzhyan:2014bfa . The equality of the coefficients found with our method and the anomaly-related method is a consequence of the equality of the Kubo formulae Amado:2011zx ; Landsteiner:Holo ; ChiralKinetic ; Landsteiner:2011cp ; Avkhadiev:2017fxj ; Golkar:2012kb obtained with either approach. In view of the explicit absence of anomaly in our calculation - we deal with free fields in flat space-time from the outset and we ignore gauge interactions - one may wonder whether this equality is accidental or if our derivation is somehow equivalent to the anomalous one for some reason.

The AVE was originally addressed for neutrinos emitted by rotating black hole in ref. VilenkinCVE . Lately, as has been mentioned, this effect was addressed in the context of CVE and quantum anomalies PhysRevD.70.074018 ; Son:2009tf ; Sadofyev:2010pr . Several terms were found to contribute to the proportionality coefficient between ${\bf j}_{A}$ and ${\bf\omega}$, in the massless limit: a term proportional to the chiral potential $\mu_{5}^{2}$, also confirmed in holographic models TorabianYee ; Erdmenger:2008rm ; Banerjee:2008th ; Amado:2011zx ; a term proportional to $T^{2}$ Neiman:2010zi ; Landsteiner:Holo ; ChiralKinetic whose existence was attributed to the gravitational anomaly Jensen:2012kj ; Landsteiner:2011cp ; Avkhadiev:2017fxj and to the modular anomaly Golkar:2015oxw . The relation to anomalies is made stronger in ref. Golkar:2012kb where it is proven that the CVE does not receive corrections from a Yukawa type interactions. Certainly, when dynamical degrees of freedom are considered instead of external fields, all kind of anomalous transport coefficients (CVE and AVE) are no more related nor constrained by anomalies, as first shown in Golkar:2012kb and then in Jensen:2013vta ; Braguta:2014gea . It is also worth pointing out that in fact, all calculations of the $T^{2}$ term in eq. (70) in refs. VilenkinCVE ; Landsteiner:Holo ; ChiralKinetic give the same result, i.e. $T^{2}/6$ like in eq. (70). We would like to point out that in the derivation of Vilenkin VilenkinCVE it is clear that the effect is caused by the modified density operator in the presence of rotation, exactly like in our case, and indeed we recover the same result in eq. (70). We also stress that in this framework, the coefficient $W_{A}$ is non-vanishing also for massive fields, as it turns out from the general formula (69), when the chiral symmetry is explicitly broken. Furthermore, it is likely that the addition of a term $\mu_{5}\widehat{Q}_{5}/T$, where $\widehat{Q}_{5}$ is a conserved axial charge, in the exponent of the statistical operator (7), will lead to a term $\mu_{5}^{2}/2\pi^{2}$ as found in PhysRevD.70.074018 ; Son:2009tf ; TorabianYee ; Erdmenger:2008rm ; Banerjee:2008th ; Amado:2011zx ; ChiralKinetic .

Finally, it should also be pointed out that the formula (69) was implicitly obtained in ref. Becattini:2014ova by means of the relativistic single-particle distribution function of particles with spin $1/2$. Therein, the mean spin tensor of the free Dirac field:

[TABLE]

was found to be, with a non-vanishing thermal vorticity $\varpi$:

[TABLE]

where:

[TABLE]

Noticing that the derivative of $\beta^{2}$ inside $F$ can be recast as a derivative over energy of the Fermi distribution $n_{F}$, after integration by parts it can be shown that

[TABLE]

according to eq. (69). Since the axial current is proportional to the dual of the spin tensor $\mathcal{S}$:

[TABLE]

using the eq. (72) and (71) the mean axial current turns out to be:

[TABLE]

where we used the (20). Hence, the eq. (66) is recovered, with the same coefficient $W_{A}$ in (69). 222While we were completing this work, a paper by Flachi and Fukushima appeared fukushima where they calculated the axial current for a rotating system in curved spacetime. It can be shown that the small mass limit of (69) corresponds to their result in flat spacetime.

VIII Limiting cases

In this section we discuss some limiting cases which may be of interest for various physical situations. It is important to stress that all of these corrections, at least for free fields, are of quantum origin, and thus are expected to contribute in limiting cases where temperature is very low and/or the chemical potential stays finite, in presence of vorticity and acceleration. Any coefficient $Y$ among those of eqs. (47),(54) for the stress-energy tensor can be generally expressed as:

[TABLE]

with $A$, $B$, $C$ some numerical constant and:

[TABLE]

where $\eta$ is $+1$ for bosons $-1$ for fermions. Also, the axial current coefficient (69) can be recast in the form (73) after integration by parts:

[TABLE]

$\underline{\bm{m=0}}$

The simplest case is the massless one, where the integral can be calculated analytically. For the free massless Boson field the only physical chemical potential value is $\mu=0$, whereas, for fermions a non-vanishing $\mu$ is possible:

[TABLE]

where $\Gamma$ is the Euler gamma function and $\mathrm{Li}_{k}$ are the Polylogarithm function Abram . The coefficients for massless Boson field at $\mu=0$ are reported in table 1 and for massless Dirac field in table 2.

$\underline{\bm{m>|\mu|>0}}$

For relativistic massive fields, one has two cases: $|\mu|<m$ and $|\mu|>m$ 333The relativistic chemical potential $\mu$ is related to the one used in non-relativistic statistical mechanics $\mu_{NR}$ by $\mu=\mu_{NR}+m$. For $|\mu|<m$ it is possible to expand the distribution function:

[TABLE]

because $|\beta|(E_{p}\pm\mu)>0$. Under this condition the distribution function $n_{F,B}$ can be expressed as a geometric series of the Boltzmann one. Introducing this expansion, changing the integration variable to the rapidity $y$ with $p=m\sinh y$, defining $x=m/T=m|\beta|$, the eq. (73) can be rewritten as:

[TABLE]

The integration can be carried out by using the integral representation of the modified Bessel function of the second kind, or McDonald functions $K_{\nu}(x)$ Abram :

[TABLE]

where the coefficients $a_{n}$ for bosons ($S=0,\,\eta=1$) are shown in table 1 and for fermions ($S=1/2,\,\eta=-1$) in table 2.

The above series is well suited to study the non-relativistic limit of the coefficients, what happens when the mass is much larger than the temperature, that is $x\gg 1$. So, by using the asymptotic expansion of the McDonald functions Abram :

[TABLE]

In this limit $m\gg T$ either the particle or antiparticle contribution can be neglected, and since $|\mu|<m$, the first term in the series (74) is the dominant one and quantum statistics reduces to the Boltzmann limit. Thus, one can write the coefficient (73) in the non-relativistic limit with $|\mu|<m$ for, e.g. particles with $\mu\geq 0$, to a very good approximation as:

[TABLE]

where $f$ is a polynomial of $1/x=T/m$ and its first leading terms are shown in tables 1 and 2.

Note that eq. (75) can be rewritten as:

[TABLE]

where ${\rm d}N/{\rm d}^{3}{\rm x}$ is the classical expression of particle density in the Boltzmann limit. This has an important consequence, that is all coefficients of the stress-energy tensor in eqs. (47) and (54) have a finite classical limit whose leading term is proportional to either mass times particle density or temperature times density. Therefore, the second-order corrections to the ideal form of the stress-energy tensor appearing in eq. (26) at thermodynamic equilibrium must be quantum Becattini:2015nva , as they vanish in the limit $\hbar\to 0$ according to the (22); this explains why genuine quadratic terms in thermal vorticity are not found in the Boltzmann kinetic approach to the gradient expansion of the stress-energy tensor jaiswal .

The same conclusion applies to the vector currents coefficients, see eqs. (63) and (65), because at the same conditions they can be expressed in a fashion similar to (75). Since vector current is odd under charge conjugation, relevant coefficients are an odd function of $\mu$ and vanish at zero chemical potential. The general coefficient $Y_{\text{curr}}$ among those of eqs. (63) and (65) can be written as, for $|\mu|<m$:

[TABLE]

where the hyperbolic sine (odd function of $\mu$) replaces hyperbolic cosine. The non-relativistic limit is very similar to that of (75):

[TABLE]

The coefficients (57) for bosons are reported in table 3, and for fermions in table 4.

$\underline{\bm{|\mu|>m>0}}$

As has been mentioned, the previous expansion is possible only when $|\mu|<m$. For the Boson gas, at fixed charge (or particle, in the non-relativistic limit) density, at some very low temperature $T\sim 0$ the chemical potential attains the limiting value $\mu=m$ ($\mu_{NR}=0$ in the non-relativistic framework), implying the onset of Bose-Einstein condensation.

In the fermion case, at very low temperature, the case $|\mu|>m$ is the so-called degenerate case. Indeed, when $\mu>m$, the Fermi-Dirac distribution function at $T=0$ becomes a step function:

[TABLE]

and the antiparticle contribution vanishes. The coefficients at zero temperature, in the degenerate case, can be expressed in terms of the parameter of $b=\mu/m$ ad they are shown in table 5.

Although all second-order coefficients vanish in the limit $T\to 0$, it is worth pointing out that a corresponding $1/T^{2}$ or $1/T$ factor appears in the quadratic terms in thermal vorticity, recalling that:

[TABLE]

as seen in Sect. IV. Therefore, all quadratic corrections to the stress-energy tensor in the (26) remain finite in the zero temperature limit in acceleration and vorticity. Particularly, all corrections to the stress-energy tensor, from table 5 turn out to be of the form $\mu^{2}a^{2}F(b)$ or $\mu^{2}\omega^{2}F(b)$ where $F(b)$ is a function of $b=\mu/m$.

In principle, these corrections might be phenomenologically relevant for very cold fermion stars, if their magnitude was comparable to the ideal term $\mu^{4}$ of energy density and pressure (see table 5). However, the typical values of the baryon chemical potential (Fermi energy), spinning frequency and gravitational acceleration of a neutron star imply a tiny ratio $a/\mu\approx\omega/\mu\approx 10^{-19}\div 10^{-27}$ and the functions $F(b)$ remain finite even in the $b\to\infty$ limit, that is for a very light fermion. Therefore, these corrections, at least for a free field, are negligible.

In general, one can argue that these corrections might be relevant for very cold massive gases subject to large accelerations and rotations. Particularly, from table 1 and 2 it can be realized that in the non-relativistic limit $m/T\gg 1$ the ratio $D_{w}/p$ is of the order 1, so that at finite $a$ or $\omega$ in the $T\to 0$ limit the contribution of the corrections to the pressure in the stress-energy tensor (26) blows up. Obviously, when the ratios $a/T$ or $\omega/T$ are ${\cal O}(1)$ the whole expansion method breaks down, but this behaviour points to a relevance of the quantum effects in the low $T$ limit for sufficiently large $a$ and $\omega$.

These corrections may also play a role in high energy nuclear collisions and Quark Gluon Plasma (QGP) physics. The very recent measurement of the $\Lambda$ hyperon polarization with respect to the reaction plane indicates a magnitude of the thermal vorticity $\varpi$ at the hadronization stage of the order of $10^{-2}$ at a centre-of-mass energy ${\cal O}(10)$ GeV. As thermal vorticity is presumably much larger in the early stage of the QGP expansion, the second-order non-dissipative corrections may be of some relevance for the hydrodynamic evolution and could compete with the first-order dissipative terms.

IX Comparison with previous determinations

As has been already mentioned, second order coefficients for the stress-energy tensor were classified in baier ; Romatschke:2009kr ; Moore:2010bu in the Landau frame. Following their notation, for a non-conformal fluid in flat spacetime, the relevant coefficients are $\lambda_{3},\,\lambda_{4},\,\xi_{3},$ and $\xi_{4}$, dubbed as thermodynamical in moore because these terms survive at thermodynamic equilibrium with rotation or acceleration. Since we adopted the $\beta$-frame and not the Landau frame, before comparing the coefficients we have to change the hydrodynamic frame. Up to second order in vorticity the relation between the fluid four velocity in the Landau frame $u_{L}$ and in $\beta$-frame $u$ is found diagonalizing the stress-energy tensor expansion in $\beta$-frame (26), see Appendix A:

[TABLE]

Moreover in baier ; Romatschke:2009kr ; Moore:2010bu a different definition of temperature and chemical potential is introduced, such that the energy and particle density are the same functions of temperature and chemical potential as in homogeneous thermodynamic equilibrium

[TABLE]

The relations between $T,\mu$ and $T^{\prime},\mu^{\prime}$ are reported in eqs. (86) and (87).

The stress-energy tensor and current expansions at the second order in thermal vorticity in the Landau frame read (see Appendix A for details):

[TABLE]

where the definitions of $D^{\prime}_{\alpha}$ and $D^{\prime}_{w}$ in terms of the $\beta$-frame coefficients are:

[TABLE]

A comparison of the above decomposition with the one in baier ; Romatschke:2009kr ; Moore:2010bu allows to extract the relations between $\lambda_{3},\,\lambda_{4},\,\xi_{3}$ and $\xi_{4}$ and our set of independent coefficients $D^{\prime}_{\alpha},\,D^{\prime}_{w},\,A$ and $W$ 444Please note that these transformation relations differs from those quoted in ref. Becattini:2015nva as the temperature and chemical potential redefinition was not taken into account therein..

[TABLE]

The above equalities can be inverted to give:

[TABLE]

The coefficients $\lambda_{3}$ and $\lambda_{4}$ were reported in refs. moore and ref. megias for free massless bosons and fermions. Our result of $\lambda_{3}=W/T^{2}$ for bosons reported in the first column of table 1 agrees with that quoted in ref. moore for a stress-energy tensor with $\xi=0$. Furthermore, in ref. moore $\lambda_{4}$ is argued to be zero in conformal case and correspond to our result for the improved stress-energy tensor $\xi=1/6$. Then, our results of $\lambda_{3}=W/T^{2}$ for massless fermions reported in the first column of table 2 -both vanishing - agree with the results quoted in ref. megias whereas they are in disagreement with those quoted in ref. moore at $\mu=0$, equal to $\lambda_{3}=T^{2}/12$.

X Conclusions

In conclusion, we have studied quantum relativistic free fields of spin 0 and $1/2$ at general thermodynamic equilibrium with non-vanishing acceleration and vorticity and we have calculated the thermodynamic coefficients of a second-order expansion of the stress-energy tensor in thermal vorticity tensor, which includes acceleration and vorticity vectors, also with a finite value of the chemical potential. We have also determined the leading order coefficients for the vector and axial currents. Such corrections may be phenomenologically relevant for system with very high acceleration, or vorticity as in the early stage of relativistic heavy ion collisions star .

We have shown, like in ref. Becattini:2015nva , that our method is very convenient to determine the coefficient of these non-dissipative (i.e. persisting in thermodynamic equilibrium) terms involving vorticity and acceleration, envisaged in the general hydrodynamic expansion of the stress-energy tensor. We have extended the results of our previous work and we have compared our results with the definitions used in the so-called Landau frame. We reinforce our previous conclusion Becattini:2015nva that these terms are of quantum nature.

We have studied the relation between axial current and vorticity known as Axial Vortical Effect (66) for the free Dirac field. The coefficient found for the massless field (70), which, in our calculation, is simply an effect of rotation at equilibrium, coincides with those quoted in literature and attributed to the gauge and gravitational anomalies as the pertaining Kubo formulae are identical. We cannot, for the present, demonstrate that the two derivations are equivalent.

Acknowledgements.

Useful discussions with K. Fukushima, A. Jaiswal and P. Kovtun are gratefyully acknowledged. E. Grossi was supported by the DFG Collaborative Research Centre “SFB 1225 (ISOQUANT)” and carried out most of his work as a post-doctoral fellow in the University of Florence supported by the Ente Cassa di Risparmio di Firenze, grant no. 20140754. This work was partly supported by the University of Florence grant Fisica dei plasmi relativistici: teoria e applicazioni moderne.

Appendix A Landau frame

In this work we have used the so-called $\beta$ becalocal or thermodynamic kovtun ; Van:2013sma frame, determined by the eq. (3) where $\beta$ is the four-temperature fulfilling the Killing equation (4). It has been shown in ref. becalocal that this frame does not coincide with Landau's frame, where the four-velocity $u_{L}$ is defined by the timelike eigenvector of the stress-energy tensor, in thermodynamic equilibrium situations where the thermal vorticity is non-vanishing. This is also apparent from the expansion of the stress-energy tensor at the second order in thermal vorticity, eq. (26) that is rewritten here for the sake of clarity:

[TABLE]

The term involving the $\gamma$ four-vector makes the four-velocity $u$, defined in eq. (3) no longer an eigenvector of $T^{\mu\nu}$, hence not the Landau frame velocity $u_{L}$. Of course, it is possible to change the frame definition and rewrite the stress-energy tensor expression (26) in the new frame, e.g. the Landau frame. This entails a transformation rule for the second-order coefficients $U_{\alpha},\,A,\ldots$ as well.

The transformation to the Landau frame requires the diagonalization of the energy momentum tensor and the determination of its unique time-like eigenvector $u_{L}$:

[TABLE]

where $\rho_{L}$ is the eigenvalue, that is the proper energy density in the Landau frame. By looking at the stress-energy tensor expression (26) it can be readily realized that the eigenvector $u_{L}$ ought to be a linear combination of the $u$ (the $\beta$ frame velocity) and $\gamma$:

[TABLE]

where $|\gamma|\equiv\sqrt{-\gamma^{2}}$ and $a$ and $b$ two unknown constants such that $a^{2}-b^{2}=1$. Contracting the eq. (26) with $u_{L}$ and using (80),(26) one obtains:

[TABLE]

which implies:

[TABLE]

where $\rho_{\rm eff}=\rho-\alpha^{2}U_{\alpha}-w^{2}U_{w}$ and $p_{\rm eff}=p-\alpha^{2}D_{\alpha}-w^{2}D_{w}$.

One can algebraically solve the equations (81) to determine $a,b,\rho_{L}$ taking the constraint $a^{2}-b^{2}=1$ into account. However, since we are dealing with an expansion of the stress-energy tensor to second order in $\varpi$, it is sufficient and more convenient to find an approximate solution of the equations at the same order in this parameter. This can be done by observing that the parameter $b$ must be ``small'' as for $\varpi\to 0$ one expects Landau and $\beta$ frame to coincide in eq. (80). Therefore, from the first of (81), one has:

[TABLE]

and, from the second:

[TABLE]

where we have kept only second-order terms in $\varpi$ keeping in mind that $\gamma={\cal O}(\varpi^{2})$. At the same order of approximation, the coefficient $a\simeq 1$.

The eq. (80) relating the Landau and $\beta$ frames is thus:

[TABLE]

and the stress-energy tensor in the Landau frame at the second order in thermal vorticity reads:

[TABLE]

We can also write the conserved vector current (56) in the Landau frame by using (82):

[TABLE]

It has become customary in the literature to include in the specification of the Landau frame a redefinition of the temperature and the chemical potential in order to avoid corrections to the equilibrium energy and charge density kovtun ; minwalla :

[TABLE]

The possibility to redefine the temperature is usually advocated in out-of-equilibrium situations, and yet its application in global equilibrium situation with $\varpi\neq 0$ is conceptually very questionable (see also discussions in refs. Becattini:2015nva ; becalocal ) because it deprives temperature one of its key relativistic features - crucial to define equilibrium - that is being the inverse of the magnitude of a Killing vector field. If one pursues the implementation of (85) anyway,

[TABLE]

where we have used the previous results, that is $\rho_{L}=\rho_{\rm eff}$ and the eq. (84). In order to find the $T^{\prime}$ and $\mu^{\prime}$ we can perform a Taylor expansion of the temperature and chemical potential in powers of $\alpha^{2}$ and $w^{2}$:

[TABLE]

which, once plugged into the (86) yield, after a Taylor expansion at the second order in $\varpi$:

[TABLE]

Equating the coefficients of $\alpha^{2}$ and $w^{2}$ on both sides, we obtain the solution:

[TABLE]

We can replace these derivatives into the (87) to obtain the relation between the proper temperature and chemical potential and the new $T^{\prime}$ and $\mu^{\prime}$. These relations allow to express all the thermodynamic functions with the new arguments. Clearly, we can neglect any term beyond the second order in $\varpi$. Particularly, in the stress-energy tensor expression at the second order (83), the only relevant thermodynamic function which gets modified is the one involving pressure,

[TABLE]

so that (83) and (84) become:

[TABLE]

where $D^{\prime}_{\alpha}$ and $D^{\prime}_{w}$ read:

[TABLE]

Appendix B Relations between coefficients

Herein we derive the relations between coefficients (31) enforcing the continuity equation for the mean value of the stress-energy tensor $T_{\mu\nu}$ at second order in the thermal vorticity $\varpi$

[TABLE]

Firstly, we observe that scalar thermodynamic functions depend on spacetime coordinates only through the magnitude of the four-temperature $F(|\beta|)=F(\sqrt{\beta^{2}(x)})$, thus:

[TABLE]

where we have used the definition of $u$ in eq. (3), the eqs. (6) and (20).

We start by reckoning the gradients of the four-vectors $\{u,\alpha,w,\gamma\}$ which are needed to calculate the stress-energy tensor divergence. First, we observe that, because of the Killing equation (4), and using the eqs. (3), (6) and (20):

[TABLE]

so that:

[TABLE]

and, contracting with $u$ again:

[TABLE]

what we already saw in section (IV). This equation implies that any scalar function, whose argument are $\beta^{2}$ and $\xi$, has a vanishing derivative along the flow, that is:

[TABLE]

Now, we can find the gradient of $u$ as defined by the eq. (3):

[TABLE]

where we have used again the eq. (3) and the eq. (90). The divergence of $u$ then is:

[TABLE]

because of the Killing vector equation (4), which obviously imply $\partial_{\mu}\beta^{\mu}=0$ and the (91). Instead the derivative along its direction is

[TABLE]

Then, let us calculate the derivative of $\alpha$, keeping in mind that $\varpi$ is constant:

[TABLE]

using the thermal vorticity decomposition (19) and the Levi-Civita tensor properties we can express the previous formula in terms of the tetrad vectors

[TABLE]

whence we obtain the divergence:

[TABLE]

as well as the gradient of $\alpha^{2}$:

[TABLE]

taking into account that $\alpha\cdot u=\alpha\cdot\gamma=0$.

Likewise, we can calculate the derivative of $w$ by using its definition (20) and the (93):

[TABLE]

Replacing the eq. (19) in the first term of the right hand side and using the properties of the Levi-Civita tensor, it can be shown that:

[TABLE]

so that its divergence is

[TABLE]

and the gradient of $w^{2}$:

[TABLE]

In order to calculate the derivative of the last relevant vector field $\gamma$ (21), one can first show that, by using $\Delta_{\lambda\mu}=g_{\lambda\mu}-u_{\lambda}u_{\mu}$ and (19), it can be expressed as:

[TABLE]

so that its divergence vanishes

[TABLE]

where we have used the eq. (96). Another useful relation involving $\gamma$ is the contraction:

[TABLE]

where we have used the decomposition (19). Then, because of the Levi-Civita tensor properties we find:

[TABLE]

We also need the derivative of the transverse projector $\Delta$:

[TABLE]

whence:

[TABLE]

Finally, we observe that

[TABLE]

as they are scalar functions.

To calculate the divergence of the stress-energy tensor we need to work out some intermediate relations involving the derivatives of the four-vectors $u,\alpha,w,\gamma$. The first relation can be obtained by using the (103) and (98):

[TABLE]

The second by using the (103) and (101):

[TABLE]

Moreover with (97) and orthogonality properties:

[TABLE]

and, using (99):

[TABLE]

Futhermore, taking advantage of (93) and (102):

[TABLE]

The last needed relation involves the gradients of $\gamma$ in eq. (21). By expanding its definition and using the previous eqs. (93,97,99):

[TABLE]

where the known contractions of Levi-Civita tensor have been employed in the last equality.

We are now in a position to enforce the continuity equation $\partial^{\mu}T_{\mu\nu}=0$, to the expression of the stress energy tensor (88). The first term involving energy density gives rise to:

[TABLE]

using eqs. (89), (94) and (95). Then, keeping in mind the (92) which implies that $DU_{\alpha}=DU_{w}=\ldots=0$ for all coefficients, and using (89), (94), (95), as well as $D\alpha^{2}=Dw^{2}=0$:

[TABLE]

Now, using eqs. (89) and (103):

[TABLE]

and, thanks to eqs. (89) and (104):

[TABLE]

likewise, because of (89) and (105):

[TABLE]

Furthermore, using (89), (94) and (106):

[TABLE]

and, similarly, by means of (89), (100) and (107):

[TABLE]

Finally, using again (89) and the relations (108) and (109) we can determine the term involving the gradients of $\gamma$:

[TABLE]

The divergence of the stress-energy tensor is obtained by summing all right-hand sides of the eqs. (110-118). As it can be readily checked, the resulting expression involves the sum of terms multiplying $\alpha_{\nu},\,\alpha^{2}\alpha_{\nu},\,w^{2}\alpha_{\nu}$ and $(\alpha\cdot w)w_{\nu}$. As the divergence should vanish independently of the vectors $w$ and $\alpha$, which have as many degrees of freedom as the thermal vorticity (that is 6), the conclusion is that each coefficient of the above combinations must be zero. As a result, four equations are obtained. The coefficient of $\alpha_{\nu}$ gives rise to:

[TABLE]

which is but the well known thermodynamic relation between energy density and pressure. The other three coefficients yield three conditions on the second-order coefficients:

[TABLE]

[TABLE]

and:

[TABLE]

Subtracting eq. (120) from eq. (119) we can simplify some terms and obtain:

[TABLE]

which completes the proof of the relations in eq. (31).