Helical massive fermions under rotation

The properties of a massive fermion field undergoing rigid rotation at finite temperature and chemical potential are discussed. The polarisation imbalance is taken into account by considering a helicity chemical potential, which is dual to the helicity charge operator. The advantage of the proposed approach is that, as opposed to the axial current, the helicity charge current remains conserved at finite mass. A computation of thermal expectation values of the vector, helicity and axial charge currents, as well as of the fermion condensate and stress-energy tensor is provided. In all cases, analytic constitutive equations are derived for the non-equilibrium transport terms, as well as for the quantum corrections to the equilibrium terms (which we derive using an effective relativistic kinetic theory model for fermions with helicity imbalance) in the limit of small masses. In the context of the parameters which are relevant to relativistic heavy ion collisions, the expressions derived in the massless limit are shown to remain valid for masses up to the thermal energy, except for the axial charge conductivity in the azimuthal direction, which presents strong variations with the particle mass.


Introduction
Quantum systems undergoing rigid rotation have been the object of academic study since the late '70s, when Vilenkin showed that, due to the spin-orbit coupling, more anti-neutrinos would be emitted on the direction which is parallel to the angular momentum vector of a rotating black hole. Conversely, an excess of neutrinos would be emitted on the antiparallel direction [1]. Recently, evidence from relativistic heavy-ion collisions experiments suggests that a global polarisation of the medium formed following the collision can be highlighted by looking at the decay products of the Λ hyperons [2,3]. Such measurements can provide quantitative validation of models describing polarization of strongly interacting matter induced through various mechanisms [4], such as the chiral vortical effect [5][6][7] and the chiral magnetic effect [8,9].
The past decade has seen many attempts to extend the standard thermal field theory, in which the canonical thermodynamic variables are the temperature four-vector β µ = u µ /T (given as the ratio between the local fluid four-velocity, u µ , and its temperature, T ) and chemical potentials related to the electric or baryon charges, to also incorporate polarisationrelated charges. We mention the recent results in Ref. [10], where the excess is taken into account by means of an axial chemical potential, µ A , coupled with the axial charge current; and Ref. [11], where it is suggested that an extra chemical potential, called the spin potential, could explain a polarisation excess through a coupling with the spin tensor.
There are two major drawbacks of the aforementioned approaches. First, the axial current J µ A = ψγ µ γ 5 ψ can be used to consistently define an axial charge, Q A , only for massless fermions. In the case of massive particles, the axial current is no longer conserved. Second, the spin tensor (or its contraction with a constant spin potential) does not commute with the Dirac equation, which implies that its eigenfunctions are not solutions of the Dirac equation.
In this paper, an alternative to the axial charge current and spin tensor extensions of canonical thermal field theory is employed, which is based on the helicity charge current, introduced in Ref. [12]. The advantage of the present approach is that the helicity charge current remains conserved for massive fermions and thus, one can consider the helicity operatorĥ (related to the temporal component of the Pauli-Lubanski vector operator) as a natural alternative to the chirality operator for the regime of non-vanishing mass. It is defined through:ĥ where J and P are the total angular momentum and linear momentum operators, while p = |p| is the momentum magnitude. The eigenvalues ofĥ are 1/2 and −1/2, distinguishing between right-handed and left-handed particles, respectively. At the relativistic quantum mechanics level, the helicity operator commutes with the Dirac equation and the transformation ψ → e 2iαĥ ψ represents a symmetry transformation of the Dirac Lagrangian. These helicity transformations form an abelian group, denoted here using U (1) H . Using Noether's theorem, it is straightforward to use the U (1) H symmetry to define a helicity charge current, J µ H = 2ψγ µ hψ, (1.2) which is conserved (the extra factor of 2 is used to normalise the eigenvalues of h to unity). The associated time-independent charge, Q H , can be assumed to couple with a helicity chemical potential µ H to account for an overall helicity bias. This paper presents an analysis of the quantum thermal expectation values (t.e.v.s) of the vector and helicity charge currents (VCC and HCC), axial charge current (ACC), fermion condensate (FC) and stress-energy tensor (SET) corresponding to the free massive Dirac field under rotation at finite vector and helicity chemical pontentials. Using analytical methods, exact (closed form) results are obtained. Furthermore, a simple kinetic model based on the Fermi-Dirac distribution is proposed, which takes into account the helicity bias by means of µ H . Unsurprisingly, the classical analysis predicts a perfect fluid behaviour, as expected since rigid-rotation is a global thermodynamic equilibrium solution in relativistic kinetic theory. This thermal equilibrium state is fully characterised by means of the vector Q RKT V and helicity Q RKT H charge densities, the energy density E RKT and the pressure P RKT . At the quantum level, the hydrodynamical content of the charge currents and the SET is extracted with the aid of the β (thermometer) frame [13][14][15], in which the four velocity is that corresponding to rigid rotation. In the β frame, the quantum corrections ∆Q V /H , ∆E and ∆P are highlighted with respect to the classical equilibrium quantities Q RKT V /H , E RKT and P RKT , respectively. These quantum corrections are quadratic with respect to the vorticity parameter Ω. We also find terms that deviate from the perfect fluid form, which describe transport phenomena. The transport terms appearing at the level of the charge currents and of the heat flux are characterised using charge and heat conductivities, σ and κ. These conductivities are expressed in terms of the local properties of the fluid (chemical potential and temperature) by means of simple constitutive relations. The validity of these constitutive relations is probed in the case of non-vanishing mass in the parameter regime which is typical of heavy ion collisions and discuss their robustness both with respect to the increase of the mass, chemical potentials, vorticity and temperature.
In anticipation of the results which will be derived later on, the form of the t.e.v.s of the vector and charge currents is reproduced below: as summarised in Eq. (6.31). The conductivity σ τ V /H = µ V /H /6π 2 characterises the charge transport along the vector τ µ , which lies in the t − ϕ plane [see Eq. (3. 12)]. As expected, the vector and helicity current conductivities are proportional to the respective chemical potentials. Contrary to expectation, the conductivities along the vorticity vector ω µ behave quite differently: The above result [also discussed in Eq. (6.32)] shows that a helicity imbalance can induce an electric current and conversely, a net electric charge can induce a helicity current. This kind of coupling can lead to the development and propagation of the helical vortical wave, as indicated in Ref. [12], which are similar to, e.g., the chiral magnetic waves [4]. Since this paper is focussed mainly on the transport properties of free massive Dirac fermions at finite temperature and chemical potentials while undergoing rotation, we defer the analysis of such phenomena for future publications. Before ending the introduction, it is worth pointing out that the downside of our proposed formulation is that the helicity is not a Lorentz invariant property, meaning that the theory becomes frame-dependent. The relevance of such a formulation can be seen in the case of systems that explicitly break Lorentz invariance, such as the confined state encountered in a single nucleus or the QGP undergoing rigid rotation.
The outline of the paper is as follows. The helicity current is introduced in Sec. 2. The kinematic tetrad consisting of the local velocity u µ , acceleration a µ , vorticity ω µ and a fourth vector τ µ is introduced for the case of rigid rotation in Sec. 3. The kinetic theory model taking into account the helicity bias is formulated in Sec. 4. The finite temperature field theory formalism employed in this paper is summarised in Sec. 5 and the t.e.v.s of the vector and helicity charge currents (VCC and HCC), axial charge current (ACC), fermion condensate (FC) and stress-energy tensor (SET) are discussed in Sections 6, 7, 8 and 9, respectively. Section 10 concludes this paper. Throughout this paper, we employ Planck units, such that = k B = c = 1, as well as the (+, −, −, −) signature for the Minkowski metric.

Helicity current of fermions
This section begins with a brief overview of the vector and axial charge currents, derived based on the U (1) V and U (1) A abelian symmetries of the Dirac Lagrangian, which are presented in Subsections 2.1 and 2.2, respectively. Noting that the axial charge current is conserved only in the limit of massless fermions, the helicity charge current is introduced in Subsec. 2.3 by means of another abelian symmetry, denoted U (1) H , which holds for any value of the fermion mass.

Vector charge and current
The free Dirac field is described by the following Lagrangian: where M is the mass of the field quanta, ψ = ψ † γ 0 is the Dirac adjoint, ψ is the Dirac 4-component spinor and γ µ are the 4 × 4 gamma matrices. For definiteness, these matrices are taken in the Dirac representation: where σ i are the Pauli matrices, given by: 3) The Dirac equation for the spinor ψ and its adjoint ψ can be obtained using the Euler-Lagrange formalism, starting from the Lagrangian in Eq. (2.1): where the derivative ← − ∂ µ acts to the left. The Dirac inner product of two 4-spinors ψ and χ can be introduced as: ψ, χ = d 3 x ψγ 0 χ, (2.5) being time independent when ψ and χ satisfy the Dirac equation (2.4). The Dirac Lagrangian is invariant under the following U (1) V abelian transformation: The above symmetry gives rise via Noether's theorem to the vector charge current (VCC) which is conserved when ψ is a solution of the Dirac equation (2.4). Indeed, taking the divergence of J µ V yields ∂ µ J µ V = 0, since γ µ ∂ µ ψ = −iM ψ and ∂ µ ψγ µ = iM ψ by virtue of Eq. (2.4). The total vector charge Q V associated to J µ V , defined by is time independent in the sense that ∂ t Q V = 0. After second quantization, ψ is promoted to a quantum operator Ψ, which can be expressed with respect to a complete set of particle (U j ) and anti-particle (V j ) modes as follows: whereb j andd † are the particle annihilation and antiparticle creation operators, respectively. These operators satisfy canonical anticommutation relations, i.e.
while all other anticommutators vanish.
Demanding that the modes be normalised with respect to the Dirac inner product (2.5), the vector charge operator can be expressed as: where the colons :: indicate normal (Wick) ordering, which, for operators that are quadratic in the one-particle operatorsb j andd j , amounts to subtracting the vacuum expectation value: : Equation (2.12) shows that each (anti-)particle in a given quantum state adds (subtracts) one unit to the expectation value of the vector charge operator, : Q V :.

Axial current for massless fermions
In addition to the vector current (2.7), one may also introduce the axial current: where γ 5 = iγ 0 γ 1 γ 2 γ 3 is the chirality operator, which in the Dirac representation of the gamma matrices (2.3) is given by: The axial current can be derived using Noether's theorem starting from the abelian U (1) A symmetry of the Dirac Lagrangian, which is exact only in the case of massless fermions, when M = 0. In this case, it can be checked that the axial current (2.14) satisfies since γ 5 anti-commutes with the Dirac gamma matrices (γ µ γ 5 = −γ 5 γ µ ), while iγ µ ∂ µ ψ = 0 for massless particles. The total axial charge corresponding to J µ A , is a time independent quantity in the sense that ∂ t Q A = 0. It is worth noting that for non-vanishing mass, M = 0, the axial current (2.14) is no longer conserved at the level of the classical equations of motion: Before promoting the axial charge to a quantum operator, the particle mode solutions U j introduced in Eq. (2.9) are taken to be eigenfunctions of the chirality operator, in the sense that where the chirality number takes the values χ j = ±1. When U j satisfies Eq. (2.20), the corresponding anti-particle modes V j = iγ 2 U * j satisfy In the second-quantization formalism, the Wick-ordered axial charge operator can then be written as follows: : where the expansion (2.9) was used for the field operator.

Helicity current
The helicity of Dirac fermions can be introduced as the eigenvalue of the helicity operator, defined as [16]: where p = |p| represents the momentum magnitude and the spin operator S has the following components: In Eq. (2.23), W 0 is the zeroth component of the Pauli-Lubanski vector, which is defined as: 25) which is expressed via the total angular momentum operator J αβ , the momentum operator P λ = i∂ λ and the Levi-Civita tensor ε µαβλ , which is introduced such that ε 0123 = +1. It can be shown that the eigenvalues of h are λ = ±1/2. The mode solutions of the Dirac equation can be chosen as eigenfunctions of h, satisfying where the antiparticle spinors V j = iγ 2 U * j , obtained via charge conjugation from the particle spinors U j correspond to the same eigenvalue λ j = ±1/2.
Noting that the helicity operator h (2.26) can be written as: and taking into account the Dirac equation (2.4), it can be shown that The presence of the Hamiltonian H = P 0 = P 0 = i∂ t prompts us to consider the mode solutions that are simultaneous eigenvectors of this operator: In the massless limit, when E = p, it can be seen that The above eigenvalue equations show that the set of eigenvectors of the helicity operator coincides with that of the chirality operator (provided they are also eigenvectors of the Hamiltonian). The action of the chirality operator on the helicity eigenvectors in the limit of vanishing mass is shown using explicit mode solutions in Eq. (5.5).
The central element of the present paper is the helicity charge current (HCC): where the factor of 2 was introduced to normalise the eigenvalues of h to ±1. Since [P µ , W ν ] = 0, it can be seen that the divergence of J µ H is given by: Using the expression (2.27) for h, the following commutators can be computed: Thus, it can be seen that The above relation confirms that J µ H is conserved at the level of classical equations of motion, allowing the corresponding time-independent helicity charge to be introduced via: The HCC in Eq. (2.31) can be derived via Noether's theorem by noting that the Dirac Lagrangian (2.1) is invariant under the following helicity phase transformation: The above transformations form an abelian group, denoted U (1) H , which represents a symmetry group of the Dirac Lagrangian for any value of the fermion mass M . Considering now the mode decomposition (2.9) of the field operator Ψ, the helicity charge operatorQ H can be written as: After normal ordering, the following expression is obtained:

Kinematics of rigid rotation
The analysis in this paper is focussed on states undergoing rigid rotation about the z axis.
The properties of such states are most conveniently expressed with respect to cylindrical coordinates. For future convenience, the following tetrad is introduced: ωt =dt, ωρ =dρ, ωφ =ρ dϕ, ωẑ = dz. In what follows, hatted indices are employed to refer to vector (or tensor) components expressed with respect to the above tetrad. The four-velocity of a fluid undergoing rigid rotation can be written as: such that its components with respect to the coordinate indices and with respect to the tetrad are: The Lorentz factor Γ, diverges as ρ → ρ SLS , where ρ SLS is the distance from the rotation axis to the speed of light surface (SLS), on which the rigidly-rotating fluid rotates at the speed of light. It is given by: Aside from the tetrad introduced in Eq. (3.1), it is convenient to define another tetrad, comprised of kinematic quantities, derived from the four-velocity [15,17]. Starting from the expression (3.2) for u = uαeα, the four-acceleration can be defined via: It can be seen that a · u = 0 by construction. Taking into account the following expression for the gradient of the velocity vector, the kinematic vorticity tensor can be computed as follows: The kinematic vorticity vector, ω = ωαeα, can be defined through: It can be seen that ω · u = 0 by construction. Moreover, Eq. (3.9) can be inverted to obtain: Multiplying both sides of Eq. (3.10) by ωβuγuσ gives: where the property aβ = −2uσΩσβ was used, which is valid when ∇αuα = 0. Thus, in the case of rigid rotation (when ∇αuα = 0), ω is orthogonal to both u and a. A fourth vector which is orthogonal to u, a and ω is [15] τα = εαβγσωβaγuσ, τ = −ρΩ 3 Γ 5 (ρΩet + eφ). It can be seen that the vectors u, a, ω and τ comprise an orthogonal tetrad. The norms of these vectors are given below: , (3.13) Figure 1 shows a schematic representation of these four vectors.

Relativistic kinetic theory analysis
Rigid body motion is fundamentally regarded as a solution of the fluid equations for which dissipative processes are absent. From a kinetic theory perspective, this corresponds to a state which in canonical thermodynamics is characterised via the temperature four-vector, βα = T −1 uα. The state corresponds to global thermodynamic equilibrium when βα satisfies the Killing equation. For rigid rotation characterised by the velocity in Eq. (3.2), this can be achieved when T satisfies [18]: where T 0 is the temperature on the rotation axis and Γ is the Lorentz factor given in Eq. (3.4). Recent works proposed the extension of the canonical formulation to account at the kinetic level for the degree of polarisation of the underlying quantum fluid. Starting from the Wigner function formalism, discussed in, e.g., Ref. [19], Becattini et al proposed in Ref. [20] expressions which couple the thermal vorticity, ω µν = − 1 2 (∇ µ β ν − ∇ ν β µ ), to the spin operator. Recently, Florkowski et al applied this formalism to obtain dynamic equations for the macroscopic polarisation in the frame of relativistic fluid dynamics with spin [21,22], however their analysis is performed in the frame of Boltzmann (classical) statistics.
In this section, a simple kinetic model is proposed to account for the equilibrium distribution of Fermi-Dirac particles with two possible helicities (λ = 1/2 and −1/2). In addition to the vector chemical potential, µ V , which distinguishes between particles and anti-particles, a straightforward extension of the Fermi-Dirac distribution is considered that includes the helicity chemical potential, µ H . Requiring that µ H distinguishes between polarisations as indicated by the corresponding quantum helicity charge operator Q H , defined in Eq. (2.39), the following expression is obtained: where µ λ = µ V + 2λµ H is the total chemical potential corresponding to right-handed (µ + = µ V +µ H ) and left-handed (µ − = µ V −µ H ) particles. In the above, pα and λ represent the four-momentum and helicity of the particle, which is assumed to have mass M (p 2 = M 2 ). Only one fermion species is considered, although multiple species or internal degrees of freedom, such as colour, can be accounted for at the kinetic level through an overall degeneracy factor, g s (the spin is already taken into account by considering λ = ±1/2). It is worth pointing out that Eq. (4.2) loses Lorentz covariance since the helicity λ is not a Lorentz scalar. It can be expected that a more fundamental formulation, starting from the theory of Wigner functions, may lead to a manifestly covariant equilibrium distribution, however such an analysis is beyond the scope of the present work. The model proposed herein serves just as a baseline to highlight the effects of taking into account a helicity bias using the standard chemical potential approach at the level of a classical theory.
In the absence of an established theory for this phenomena, we assume that the dynamics of the distribution of particles with spin f q/q;λ ≡ f q/q;λ (x, p, t), is given by the relativistic Boltzmann equation [18]: The collision operator J[f ] leading to the gas thermalisation should be implemented such that maximum entropy is ensured when f q/q;λ = f (eq) q/q;λ , where the equilibrium distributions is given in Eq. (4.2). In global thermodynamic equilibrium, f q/q;λ = f (eq) q/q;λ everywhere and the right hand side of Eq. (4.3) vanishes identically. Assuming that µ H does not depend on the particle momentum p, Eq.  As mentioned at the beginning of this section, the above equations are satisfied when uα/T is proportional to a Killing vector and the ratios µ V /T and µ H /T are constant. Taking the solution corresponding to rigid rotation, given in Eq. (3.2), with the temperature given through Eq. (4.1), it can be seen that the chemical potentials satisfy: where µ V ;0 and µ H;0 represent the values of the vector and helicity chemical potentials on the rotation axis. The distributions (4.2) can be specialised to the case summarised in Eq. (3. 2): where µ λ;0 = µ p;0 + 2λµ H;0 = µ λ /Γ is the total chemical potential on the rotation axis. In the above, p t denotes the co-rotating energy of the particle, while the azimuthal velocity vφ = pφ/pt is written in terms to the azimuthal component of the momentum, pφ = − sin ϕp x + cos ϕp y . Since −1 < vφ < 1, it can be seen that p t satisfies: p t > 0, (4.7) valid for any particle motion, as long as ρΩ < 1. The above inequality will be essential in Sec. 5.2, where the second quantisation of the Dirac field is discussed.
In the zero temperature limit, Eq. (4.6) reduces to: where E F q/q;λ is the Fermi level for the particle/anti-particle distributions. Thus, the particle/anti-particle distributions have non-vanishing values only when the co-rotating energy is below the Fermi level.
Starting from the distributions (4.6), the particle charge current Jα V (VCC), helicity charge current Jα H (HCC) and stress-energy tensor Tαβ (SET) can be computed as follows: where the charge densities Q RKT V /H , energy density E RKT and pressure P RKT are given by: (4.10) -12 - The above results are obtained from Eq. (4.9) in two steps. The first step consists of contracting J RKT;α p/h with uα, as well as Tασ RKT with uαuσ (for E RKT ) and with ηασ (for E RKT −3P RKT ), respectively. Afterwards, a Lorentz transformation is performed, such that u · p → pt. This transformation is permitted by the Lorentz invariant integration measure (d 3 p/pt), however, one must assume at this step that µ λ /T = (µ V + 2λµ H )/T is a Lorentz scalar. The vector part, µ V /T = µ V ;0 /T 0 , clearly satisfies this assumption. However, the second part, 2λµ H /T = 2λµ H;0 /T 0 , is not a Lorentz scalar due to the polarisation prefactor, 2λ. It is known that for massive fermions, the polarisation is frame dependent. This is generally true for the particles with azimuthal velocity pφ/m smaller than uφ = ρΩ. In the vicinity of the rotation axis, these particles populate the infrared sector of the integral in Eq. (4.10), which makes small contributions due to the factor p 2 in the integrand. One can conclude that treating 2λµ H /T as a Lorentz scalar can give a correct order of magnitude assessment of Eq. (4.9), at least in the vicinity of the rotation axis.
In the massless limit, pt = p, the helicity becomes frame-independent, such that the approximation in Eq. (4.10) becomes exact. Using the relations in Eq. (A.11), the integrals in Eq. (4.10) can be performed analytically: where T = T 0 Γ, µ V = µ V ;0 Γ and µ H = µ H;0 Γ, while E RKT = 3P RKT in general for massless constituents. For future reference, the massless limit of the ratio (E − 3P )/M 2 between the trace of the SET and M 2 was also included. This quantity will be compared with the QFT equivalent (also equal to the ratio ψψ/M between the fermion condensate and M ). It can be seen that in Q RKT V and Q RKT H , the roles of the chemical potentials µ V and µ H are reversed, while the expressions for E RKT and (E RKT − 3P RKT )/M 2 are symmetric under µ V ↔ µ H . It can be seen that, in the limit µ H → 0 of vanishing helicity chemical potential, the results for Q RKT V , E RKT and (E RKT − 3P RKT )/M 2 coincide with those presented in Eqs. (26)(27)(28) of Ref. [23], while Q RKT H vanishes. For future reference, the result for the right-handed (Q + ) and left-handed (Q − ) helicities is presented below: For small masses, a perturbative approach can be employed to derive the behaviour of the corrections to Eq. (4.11) due to finite particle mass. First, the integration variable is switched from p to x = pt/M : When the chemical potentials are non-vanishing, the following expansions can be made inside the integrands in Eq. (4.13): (4.14) The Fermi-Dirac factor is now replaced using the following series representation: The derivatives of orders k = 2j + 1 and 2j with respect to x appearing in Eq. (4.14) can be computed automatically, giving factors of (− M/T ) k . Assuming that the sums over j and over commute, the sum over j can be performed by noting that: This allows Eq. (4.14) to be written as:  It is clear that the above series converge only when xM > µ λ , however, the method can be further employed to recover the first correction due to small but non-vanishing M . After noting that the following expressions are obtained: where Q RKT H and µ ± = µ V ± µ H were introduced in Eq. (4.12). The integrals with respect to x in Eq. (4.19) can be performed by employing [24]: where K ν (z) is a modified Bessel function of the third kind and ν is assumed to satisfy Re(ν) > −1/2. Setting ν = n ∈ N yields: where (2n − 1)!! = 1 · 3 · 5 · · · (2n − 1) is the double factorial and we use the convention that (−1)!! = 1. After employing the above steps, Eq. (4.19) can be put in the following form [25]: The massless limit result given in Eq. (4.11) can be recovered from Eq. (4.22) using the limiting behaviour [24] where γ = 0.577216 . . . is the Euler-Mascheroni constant. The above small z expansion can be used to compute the contributions to Eq. (4.22) corresponding to various orders of the mass, M . This procedure can be used only for the first few terms in this series, since the expansion in Eq. (4.23) is not uniformly convergent. Indeed, since z = M/T , it is clear that the infinite sums over of the terms corresponding to positive powers of are divergent. Furthermore, the terms with z 0 in Eq. (4.23) are independent of T , µ V and µ H . With the above discussion in mind, the following O(M 2 ) corrections to the massless results in Eq. (4.11) can be obtained: (4.24) The correction to (E RKT − 3P RKT )/M 2 cannot be obtained using this method. It can be seen that the mass dependence of the pressure and charge densities is induced by terms of the form M/T = M/T 0 Γ, indicating that all higher order mass corrections make subleading contributions in the vicinity of the SLS, where Γ → ∞.

Quantum rigidly-rotating thermal states
This section begins with a review of the mode solutions of the Dirac equation which are helicity eigenvectors, discussed in Subsec. 5.1. The massless limit and the rest frame limit are also discussed therein. The second quantisation procedure and a review of the stationary (Minkowski) and rotating vacua is provided in Subsec. 5.2. The formalism employed for the computation of thermal expectation values (t.e.v.s) is briefly reviewed in Subsec. 5.3. Some analytic techniques which will be employed in the analysis of the t.e.v.s are summarised in Subsec. 5.4.

Mode solutions
At the relativistic quantum mechanics level, the Dirac equation can be solved in terms of a complete set of particle and anti-particle mode solutions, denoted using U j and V j . The particle mode solutions U j can be taken to simultaneously satisfy the following eigenvalue equations: where j = (E j , k j , m j , λ j ) denotes collectively the eigenvalues of the Hamiltonian H = i∂ t , z components of the momentum P z = −i∂ z and total angular momentumM z = −i∂ ϕ + S z , and helicity operatorĥ =M ·P /p. The explicit expressions for U j were obtained in Refs. [26,27] and are reproduced below without derivation: In the above, the sign of the energy E j is left arbitrary. The anti-particle modes V j are obtained via the charge conjugation operation: The modes U j and V j are normalised with respect to the Dirac inner product (2.5) according to: In the massless case, E ± j = 1 and it is easy to see that the particle modes U j and their respective anti-particle modes become eigenstates of the chirality operator γ 5 . For the case when E j > 0, it can be seen that thereby confirming explicitly that the chirality and helicity eigenmodes coincide. It can be seen that the chirality eigenvalues χ j satisfy χ j = 2λ j for the positive energy modes and χ j = −2λ j for the negative energy modes, as anticipated in Eqs. (2.20) and (2.21), respectively. Before concluding this subsection, it is worth discussing the properties of U j and V j in the rest frame, where k j = q j = 0. Since the Bessel functions J n (z) vanish for n = ±1, ±2, . . . , only the modes with m = ±1/2 are non-trivial. In this case, the particle modes U j ≡ U λ j E j ,k j ,m j and their corresponding anti-particle modes are given by: It is clear that the above spinors become eigenmodes of the spin operator, S z : Thus, m = ±1/2 reduces to the spin quantum number. The helicity λ = ±1/2 distinguishes between two linearly independent solutions, as predicted by the orthogonality relation in Eq. (5.4).

Second quantisation
Before promoting the wave function ψ(x) to the quantum operator Ψ, it is instructive to first discuss the definition of the vacuum state. In particular, two vacuum states are considered: the Minkowski (non-rotating) vacuum state, denoted using M , and the co-rotating vacuum state, denoted using Ω. The difference between these vacuum states is explained by Iyer in Ref. [28]. For the stationary state, the natural definition of the vacuum state is to interpret the states with positive Hamiltonian eigenvalues, E j > 0, as particle states. The states with E j < 0 are then anti-particle states. This is achieved by expanding the field operator Ψ with respect to the modes (U j , V j ) while considering only positive values of E j : The above expansion is invertible since the set of solutions where the sum over modes j θ(E j ) defined in Eq.
, it can be seen that the particle creation and annihilation operators,b M ; † j andb M j and their antiparticle counterparts satisfy the following anti-commutation relations: In the case of the rotating vacuum, by analogy to the RKT restriction on p t , given in Eq. (4.7), the particle and anti-particle modes are split with respect to the sign of the co-rotating energy This is achieved at the level of the mode sum by introducing a step function θ( E), as follows [26,28]: where the integral with respect to E runs over both negative and positive values, provided |E| > M . Given the relation (5.3) between the particle and anti-particle modes, it is not difficult to see that the completeness relation (5.9) can be recast for the rotating vacuum as follows: while the particle and anti-particle operators for the rotating vacuum satisfy the canonical anti-commutation relations: Noting that both the Minkowski and the rotating one-particle operators are obtained by taking the inner product between U j and V j with the field operator, Ψ, simple Bogoliubov relations can be established between them, as follows [26]: It is not difficult to see that the expectation values of the products of two one-particle operators corresponding to the rotating vacuum with respect to the Minkowski vacuum is just: Let us now consider an operator F which is quadratic with respect to Ψ. Considering that at the classical level, F can be written as: where F(ψ, χ) is a bilinear form involving ψ and χ, at the quantum level, F can be written as: where "mixed terms" refers to terms involving , respectively. The field operators were expanded with respect to the rotating vacuum one-particle operators, for definiteness.
The concept of normal ordering can be implemented for F with respect to the rotating vacuum state, |0 Ω , as follows: d Ω j : Taking into account that the state |0 Ω is annihilated by b Ω j and d Ω j and using the decomposition in Eq. (5.18), the following expression is obtained: Considering Wick ordering with respect to the Minkowski vacuum, one can derive: where the expectation value of : F : Ω computed in the Minkowsi vacuum state can be obtained using Eq. (5.16): where the flip j →  was considered on the second line. In general, the bilinear forms F(ψ, χ) are scalars from the point of view of the spinor indices, involving the product of ψ and χ (as well as othe spinor and/or differential operators acting on these spinors). By virtue of Eq. (5.3), it can be assumed that In this case, Eq. (5.24) becomes:

Thermal expectation values
In this paper, the finite temperature expectation values are computed by taking the thermal average over the Fock space, as follows [29,30]: where Z = Tr(ρ) is the partition function. The density operatorρ can be obtained in a straightforward fashion, by promoting the microscopic momenta from RKT to quantum operators. Taking into account that the chemical potentials µ V and µ H introduced in Eq. (4.6) are conjugate to the vector and helicity charge operators Q V and Q H introduced in Eqs. (2.12) and (2.39), respectively,ρ can be written as [31]: where H is the Hamiltonian and M z is the z component of the total angular momentum operator. In particular, the t.e.v.s ofb † where µ λ j ;0 = µ p;0 + 2λ j µ H;0 . Considering the limit β 0 → ∞ and µ V ;0 , µ H;0 → 0, it can be seen that Eq. (5.27) recovers the vacuum state only when E j > 0: It is therefore natural to discuss the construction of the thermal expectation values (t.e.v.s) at non-vanishing values of Ω by considering the normal ordering imposed via Eq. (5.12), which corresponds to the rotating vacuum state. Focussing on normal ordering with respect to the rotating vacuum state, when E j > 0 is imposed through the step function θ( E j ) in Eq. (5.12), it can be seen that at vanishing temperatures, the t.e.v.s in Eq. (5.27) vanish when E j is below the Fermi level (4.8), in agreement with the RKT theory.
Taking the t.e.v. of Eq. (5.20) using Eq. (5.27) gives: The sum over j is understood to stand for the summation and integration on the first line of Eq. (5.12). Equation (5.29) can be written such that the domain for the integration with respect to the energy E spans only positive values, i.e. M < E < ∞. To this end, the notation F m,E is introduced through: It can be seen that F m,E ( E) depends explicitly on both m and E and the dependence on E = E − Ωm is taken into account explicitly for future convenience. The integration with respect to E appearing in Eq. (5.30) can be broken into the corresponding positive and negative ranges. For the term containing the negative values of E, the simultaneous transformation (m, E) → (−m, −E) is performed, yielding: : For the particular operators considered in this paper, F(V j , V j ) = ±F(U j , U j ), such that only one of the terms above survives. Equation (5.32) forms the basis for the numerical analysis discussed in the following sections. However, further analysis of Eq. (5.32) is encumbered by the presence of the modulus | E| in the Fermi-Dirac factors, as well as by the sign function sgn( E). This can be achieved by considering the t.e.v. expressed with respect to the Minkowski vacuum, which can be computed starting from Eq. (5.21): Using Eq. (5.24), it can be shown that the t.e.v. of F computed with respect to the Minkowski vacuum is: It is worth emphasising that the difference 0 M | : F : Ω |0 M between Eq. (5.32) and (5.34) is a purely vacuum quantity, depending only on the rotation parameter Ω, which by virtue of Eq. (5.24), vanishes whenever F(V j , V j ) = F(U j , U j ). Previously, the t.e.v.s with respect to the Minkowski vacuum were computed by Vilenkin in Refs. [31,32] and indeed temperatureindependent terms were revealed in the final results. In was argued in a previous publication [26] that such temperature-independent terms are spurious and disappear when the t.e.v.s are computed with respect to the rotating vacuum.

Small mass analysis
In the small mass limit, analytic closed form expressions can be obtained following the methodology introduced in Ref. [26]. The basic idea is illustrated in Subsecs. 6.2, 7.2, 8.2 and 9.2 for the particle and helicity charge currents, the axial current, the fermion condensate and the SET, respectively. The basic idea is to expand the Fermi-Dirac distributions, which depend on the corotating energy E = E − Ωm, in a power series with respect to Ω, following the prescription: When computing the small mass limit, a further expansion with respect to M should be performed. Noting that a function f (E) that depends only on the energy (no extra dependence on p) can be expanded as: it can be seen that Eq. (5.35) can be written as: The sum over m can be performed using the following relations [26]: where the coefficients s + n,j are determined from: Note that in Ref. [26], the parameter m = 0, ±1, ±2, . . . takes integer values, while in this paper, the convention that m = ± 1 2 , ± 3 2 , . . . is an odd half-integer is employed.
It can be seen that s + n,j vanishes when j > n. For small values of n − j ≥ 0, the first few coefficients are given by [26]: For general values of n > j, the following recurrences can be established: Finally, the following formula can be employed for the integration with respect to k:

Vector and helicity charge currents
In this section, the thermal expectation values (t.e.v.s) of the vector charge current (VCC) and helicity charge current (HCC) operators are investigated. The quantum operators are defined starting from the classical expressions (2.7) and (2.31), as follows: The axial charge current (ACC) operator will be considered in Sec. 7. The general expressions that form the basis of the computation of the t.e.v.s are derived in Subsec. 6.1. Analytic expressions are derived in the small mass limit in Subsec. 6.2. The validity of these expressions at finite mass is considered using numerical integration in Subsec. 6.3.

General analysis
In the notation of Sec. 5.4, the bilinear forms introduced in Eq. (5.17) which correspond to the HCC and VCC are: is it easy to show that: where the final equality follows after noting that: Using the explicit expression for the modes, given in Eq. (5.2), the following relations can be obtained: while U j γρU j = 0 for all j. In the above, J * m (qρ) denotes: When substituting the expressions from Eq. (6.6) into Eq. (6.5), the terms proportional to k can be discarded since they are odd with respect to the transformation k → −k. The non-vanishing components of the charge currents can be summarised through: . In the following, it is convenient to refer to the charge currents for the right-and left-handed particles, Jα ± = Jα V ± Jα H , which can be computed as follows: The t.e.v.s of the charge currents can be decomposed with respect to the orthogonal tetrad formed by the vectors uα, aα, ωα and τα, introduced in Eqs. (3.2), (3.6), (3.9) and (3.12), as follows: : where Q ± represents the charge density and the charge current in the rest frame, Jα ± , is by construction orthogonal to uα. There is no term multiplying the acceleration, a = −ρΩ 2 Γ 2 eρ, since : Jρ ± : Ω = 0. The charge densities Q ± , the vortical charge conductivities σ ω ± and the conductivities σ τ ± along the τα vector can be obtained via By comparison to the RKT decomposition in Eq. (4.9), it is clear that, when they are nonvanishing, σ τ ± and σ ω ± lead to non-equilibrium currents which represent quantum corrections to the perfect fluid form.

Small mass limit
We apply the algorithm described at the end of Sec.
where the summations over j and n were interchanged and the transformation n → n + j was subsequently performed in the sum over n, in order to shift the summation range from j ≤ n < ∞ to 0 ≤ n < ∞. The integration variable was changed from E to p. For theφ andẑ components, an integration by parts was performed prior to this change of variable.
The analysis for the small mass regime can be performed starting from Eq. (5.37). Due to the nature of the integrands, it is convenient to discuss first : Jt ± : Ω and : Jφ ± : Ω . Integration by parts can be performed 2j times for the leading order (massless) term and 2j + 1 times for the O(M 2 ) term, allowing Eq. (6.12) to be written as follows: : Focussing on the integration with respect to p and noting that: where the functions I − 0 and I − 1 are given in Eq. (A.11), it can be seen that the series with respect to n appearing in Eq. (6.13) terminates after a finite number of terms. For each value of n, the summation over j can be performed by noting that: This leads to the following exact results: : The results in Eq. (6.16) can be written in terms of Q ± and j τ ± : , (6.17) where the acceleration aα and vorticity ωα vectors are defined in Eqs. (3.6) and (3.9), respectively, while the squares ω 2 = −ω 2 ≥ 0 and a = −a 2 ≥ 0 of their spatial parts are given in Eq. (3.13). In the above, the term Q RKT ± corresponds to the relativistic kinetic theory (RKT) prediction, given up to O(M 4 ) in Eq. (4.24), while ∆Q ± and σ τ ± represent quantum corrections, which are independent of the particle mass up to O(M 4 ).
Turning to : Jẑ ± : Ω , the summation with respect to m and integration with respect to k can be performed as before, yielding: : After performing the integration with respect to p, the following result is obtained: : where the expression for I − 1/2 (a) can be found in Eq. (A.20). As mentioned before, the series does not terminate, indicating that the dependence on Γ is not polynomial. In the following, only the term on the first line together with the term corresponding to n = 0 on the second line of Eq. (6.19) are considered. Performing the sum over j using Eq. (6.15), together with the relation: the following expression can be obtained for σ ω ± : where it is understood that T = T 0 Γ and µ ± = µ ±;0 Γ, while the acceleration a = −ρΩ 2 Γ 3 eρ and kinematic vorticity ω = ΩΓ 2 eẑ are given in Eqs. (3.6) and (3.9), respectively. The high temperature limit of Eq. (6.21) can be computed by noting that: The following result is obtained: (6.23) The result derived in Eq. (6.21) is valid only for small values of the rotation parameter. An expression which is exact on the rotation axis can be derived starting from Eq. (6.8), by noting that J − ±1/2 (0) = ±1, such that, where the following notation was employed: The techniques introduced for the analysis of the case of non-vanishing mass in the RKT formulation, in Sec. 4, can now be employed. First, the integration variable is changed to x = E/M , as indicated in Eq. (4.13). Then, the formula on the second line of Eq. (4.17) can be used to transform the Fermi-Dirac factors, allowing Eq. (6.24) to be written as: : It is not difficult to perform the integration with respect to x, giving: : The sum over can be performed in terms of the polylogarithm function, allowing σ ω ± = : Jẑ ± : Ω /ΩΓ 2 to be expressed on the rotation axis as: where ζ = (ς ± µ ± + ς Ω Ω 2 − M )/T . The above expression is exact for any value of the mass. The large temperature limit can be extracted using Eq. (6.22) together with: The result is: which is consistent with the expression derived in Eq. (6.23). Both Eqs. (6.23) and (6.30) show that at large temperature, σ ω ± scales linearly with the temperature T and the local chemical potential µ ± . It is noteworthy that σ ω H = 1 2 (σ ω + − σ ω − ) 2 π 2 µ V T ln 2 is generated by the vector chemical potential µ V , and vice-versa for σ ω V = 1 2 (σ ω + + σ ω − ) 2 π 2 µ H T ln 2. The result is summarised below: The quantities referring to the vector and helicity charge currents can be obtained by adding or subtracting the above expressions, e.g.
In the case of the vortical conductivity σ ω ± , the ± sign on the right hand side indicates that σ ω V ∼ µ H and σ ω H ∼ µ V , i.e.: At vanishing mass and helicity chemical potential, the results in Eq. (6.31) coincides with those reported in Eq. (4.9) and Table 2 of Ref. [10] when the axial chemical potential is assumed to vanish.

Numerical analysis
In this section, the validity of the constitutive equations derived in the previous subsection for the quantum corrections are analysed at finite mass. The explorations presented in this section are based on the region of the parameter space corresponding to the conditions which are prevalent in ultrarelativistic heavy ion collisions, namely T = 150 MeV [33], µ = 30 MeV [34] and Ω = 6.6 MeV (corresponding to 10 22 s −1 [2,35]). These values of the parameters are labelled HIC throughout this paper. The convention when presenting results for different parameters is to start from these HIC values and change only one parameter per curve. The parameter change is shown in the plot legend in the form parameter × multiplier. E.g., a curve corresponding to a temperature which is twice that corresponding to the HIC parameters (T = 300 MeV) is labelled as T × 2.
The discussion in this section begins by considering the behaviour of the correction ∆Q ± to the charge density and the conductivity σ τ ± as the mass is increased. The analysis in the previous subsection showed that there are no contributions to ∆Q ± and σ τ ± from the mass term up to O(M 4 ). For each curve, only one parameter from the original list is increased. This parameter is indicated in the legend, together with the corresponding multiplier. The second curve (green and empty circles), corresponding to Ω 2 GeV, shows that the validity domain of Eq. (6.31) is enhanced at higher values of Ω, with deviations occurring at M Ω/4. The third curve (blue and filled circles) shows deviations from the massless prediction also when M T = 750 MeV. Finally, the last curve (red line and empty triangles) curresponds to µ ± = 3 GeV T = 150 MeV. In this case, the gas is strongly degenerate, such that a strong suppression can be seen when M exceeds µ. Also, for M c 2 < µ, the constitutive relations for the massless case retain their validity, except in the vicinity M µ. Here, a strong deviation can be seen, indicating an unexpected resonance. While this regime does not seem to be of immediate relevance to the field of relativistic heavy ion collisions, it is worth exploring its origin.
The starting point for the analysis of the vanishing temperature regime is to take the β 0 → ∞ limit of Eq. (6.9). Restricting the analysis to the axis of rotation, the following resutl is obtained for the quantum correction to the charge density: where the last line corresponds to the RKT contribution Q RKT ± , which can be obtained by taking the vanishing temperature limit of Eq. (4.9): When |µ ± | − |Ω| 2 < M < |µ ± | and the term on the second line in Eq. (6.33) cancels, ∆Q ± admits a global maximum with respect to M , which is obtained by solving d∆Q ± /dM = 0: It is understood that this maximum occurs only when |µ| > |Ω|/2. This behaviour is confirmed in Fig. 3(a), where the ratio ∆Q ± (M )/∆Q ± (0) is represented at fixed Ω with respect to the ratio M/Ω, for various values of µ. The dotted black line represents the value of this ratio when |µ| = µ max , where When |Ω/µ| 1 and M max → |µ| and ∆Q ± (M max )/∆Q ± (0) peaks sharply around M = M max : The development of this peak is highlighted in Fig. 3(b). For the case considered in Fig. 2, the ratio between the chemical potential and the rotation parameter is µ ± /Ω 450. In this  case, Eq. (6.37) predicts that ∆Q ± (M max )/∆Q ± (0) 16.4, compared to 1.5 observed in Fig. 3(a). This discrepancy is an indication of the thermal suppression of this effect, which is already significant when T 0 /µ = 0.05.
We now consider the validity of the results derived for the charge current vortical conductivity σ ω ± . An exact result was obtained in Eq. (6.28) for this quantity evaluated on the rotation axis, which is valid at any mass. Figure 4(a) shows that, in the vicinity of the HIC parameters, the relative mass correction 1 − σ ω (M )/σ ω (0) is well represented by the quadratic term 3M 2 /24T 2 ln(2) highlighted in Eqs. (6.23) and (6.30).
Away from the rotation axis, the result for σ ω ± was obtained only up to O(Ω 4 ) in Eq. (6.21). This expression, involving the polylogarithm function, is approximate even on the rotation axis. To test its validity for relativistic heavy ion collisions, various levels of approximation σ ω;n ± are considered based on Eq. (6.21), distinguished using the index n. The n = 1 approximation corresponds to the leading order with respect to the temperature, highlighted in Eq. (6.23) and reproduced below: As can be seen from Fig. 4(b), σ ω;1 ± offers a very good approximation in the conditions that are prevalent during heavy ion collisions. The mass corrections are important only when M T , even for T = 750 MeV. However, when µ or Ω become comparable with T , deviations can be seen. In order to understand the nature of these deviations, the next to leading order is considered below, which can also be read from Eq. (6.23): In the above, the mass contribution was not taken into account. As can be seen from Fig. 4(c), this approximation improves the small M results at moderate values for Ω and µ. However, increasing µ and Ω further brings about new discrepancies. The final approximation, σ ω;3 ± , corresponds to the result in Eq. (6.21) minus the mass term: (6.40) Figure 4(d) shows that the above approximation correctly accounts for the chemical potential, even when µ/T 1. However, no improvement can be seen for the case of large values of Ω. The reason why no improvement can be seen is that Eq. (6.21) was obtained following a truncation at (and including) order O(Ω 2 ).

Axial current
Initially shown to be so by Vilenkin [1], it is known that the local vorticity induces an axial current parallel to the vorticity vector through the chiral vortical effect [4]. This section presents an analysis of the effect of helicity imbalance induced by the helicity chemical potential µ H on the thermal expectation value (t.e.v.) of the axial charge current (ACC) operator Jα A , defined by promoting the wavefunction in Eq. (2.14) to a quantum operator: The analysis begins with a general discussion of the t.e.v. of the ACC, presented in Subsec. 7.1, followed by an analytical analysis of the small mass regime in Subsec. 7.2 and a numerical analysis of the t.e.v.s at finite values of the mass, presented in Subsec. 7.3.

General analysis
The bilinear form introduced in Eq. (5.17) corresponding to the ACC is the t.e.v.s of the components of the ACC taken with respect to the rotating vacuum can be written as: The t.e.v.s expressed with respect to the Minkowski vacuum are: As was the case for the VCC and HCC, the t.e.v. of the ACC can be decomposed as: :

Small mass limit
The t.e.v.s of the components of the ACC taken with respect to the Minkowski vacuum can be put in the form: For the z component of the t.e.v. of the AC, the massless limit can be derived in closed form. The M 2 correction is inaccessible using the method outlined in Eq. (5.37) due to an infrared divergence of the j = 0 term. After changing the integration variable to dp = E p dE, the procedure introduced in Eq. (5.35) can be used to obtain: : where the coefficient of M 2 was obtained after noting that E = p + M 2 2p and 1 E = 1 p − M 2 2p 3 . In the term proportional to M , the shift j → j + 1 must be performed. After performing integration by parts 2j times, the following formula can be used for the integration with respect to p: ∞ 0 dp p d 2n dp 2n 2I + 0 (βµ λ )/β 2 0 , n = 0, 1, n = 1, 0, n > 1, (7.10) where I + 0 (a) is given in Eq. (A.11). It is not difficult to see that the sum over n terminates at n = 1 for the first term and at n = 0. The exact result is: : where a = −ρΩ 2 Γ 2 eρ and ω = ΩΓ 2 eẑ are the acceleration and vorticity vectors, introduced in Eqs. (3.6) and (3.9), respectively, while their squares ω 2 = −ω 2 ≥ 0 and a 2 = −a 2 ≥ 0 are given in Eq. (3.13). The second term inside the square brackets in Eq. (7.11) is a purely vacuum contribution which is not taken into account when expressing the conductivity σ ω A .
It is worth comparing the result in Eq. (7.11) with the one obtained in Ref. [36] using the Wigner function proposed in Ref. [20]: Noting that ω 2 = −ω 2 and a 2 = −a 2 , it can be seen that the parts which depend on T and µ = µ V reproduced in Eq. (7.12) agree with those obtained in Eq. (7.11), when µ H = 0.
There is a discrepancy in the vacuum term, equal to which may be due to a fundamental difference between the formulation based on the Wigner function (employed in Ref. [36]) and the QFT formulation employed in this work.
As opposed to the case of theẑ component, the sum over n in Eq. (7.8) no longer terminates when considering thet andφ components, which is an indication that their dependence on Γ and Ω is non-polynomial (see also the discussion in Subsec. 6.2). Performing the change of integration variable from E to p (dp = EdE/p) in Eq. (7.8) and using Eq. (5.37) gives: : 1 e β 0 (p−µ λ,0 ) + 1 + 1 e β 0 (p+µ λ,0 ) + 1 , : 1 e β 0 (p−µ λ,0 ) + 1 + 1 e β 0 (p+µ λ,0 ) + 1 . (7.14) Taking into account the n = 0 and n = 1 terms gives: : : The charge Q A and the conductivity σ τ A are The sum over λ can be performed by taking into account the expressions for I + ±1/2 given in Eq. (A.20): The high temperature limit of the above expressions is: It can be seen that the mass term makes a significant contribution to σ τ A when Ω is sufficiently small. A more quantitative assessment will be presented in the following subsection.
The results in Eq. (7.17) are valid only at small values of Ω. It is possible to obtain exact expressions for Q A and σ τ A on the rotation axis. These quantities require the evaluation of : Jt A : and 1 ρΩ : Jφ A : in the limit ρ → 0. Noting that otherwise, (7.19) it is not difficult to obtain the following expressions: : : where x = E/M and Eq. (4.17) was used to expand the Fermi-Dirac factors. The integration with respect to x can be readily performed and the sums over yield polylogarithm functions. The result can be summarised as: where ζ = (σ Ω Ω 2 + σ V µ V + σ H µ H − M )/T . The above expressions are exact for any particle mass. Extracting the large temperature limit yields: The above results are in agreement with Eq. (7.18). The high temperature limit of the results obtained in this section can be summarised as follows: The mass term was retained in the coefficient of τα due to its unusual effect. It is worth noting that, in the limit µ H = M = 0, the results in Eq. (7.23) reproduce those reported in Eq. (4.14) and Table 2 of Ref. [10] for the case of a vanishing axial chemical pontential (µ A = 0).

Numerical analysis
In this section, the validity of the results summarised for the small mass limit in Eq. (7.23) is investigated as the mass is increased. As in Sec. 6.3, the discussion is focussed on parameters which are pertinent to heavy ion collisions (HIC). In the case of the axial current, the vector and helicity chemical potentials play equivalent roles. Thus, for simplicity, the convention µ V = µ H = µ is used throughout this section. First, the validity of the constitutive equation for σ ω A derived in the small mass limit in Eq. (7.23), is considered as the mass is increased. At nonvanishing M , σ ω A (M ) is evaluated numerically starting from Eq. (7.4). For simplicity, the analysis is limited to the case of the rotation axis. Figure 5(a) confirms that the relative mass correction, 1 − σ ω A (M )/σ ω A (0), is of fourth order with respect to M . Around the parameters relevant to heavy ion collisions (labelled HIC), the mass term makes a relative contribution of the form αM 4 /Ω 2 T 2 , where α 0.150 is a dimensionless number. Figure 5(b) presents the ratio σ ω A (M )/σ ω A (0). It can be seen that the constitutive relation holds for M 100 MeV, which is below the thermal energy (T = 150 MeV). Increasing the temperature by a factor of 5 (blue line with filled circles) increases the validity up to M 500 MeV, i.e. by a factor of 5. At higher chemical potentials (µ V = µ H = µ = 3 GeV), it can be seen that the constitutive relation remains valid until the mass approaches the Fermi level, which is given by ∼ µ V +µ H = 6 GeV. The ratio σ ω A (M )/σ ω A (0) seems to have a monotonic behaviour, as also observed for the ratio σ ω ± (M )/σ ω ± (0), shown in Fig. 4. An interesting effect can be observed when the rotation parameter is increased by a factor of 300, to Ω 2 GeV. At such a high value, Ω/2 acts like a Fermi level and the maximum observed in Fig. 5 resembles the one highlighted for the ratios ∆Q ± (M )/∆Q ± (0) and σ τ ± (M )/σ τ ± (0) in Fig. 3. However, in the latter case, the maximum was observed at large chemical potential. Now, the axial charge density Q A and the conductivity σ τ A along τ are discussed. As opposed to σ ω A , it is not possible to obtain the radial profiles of these quantities analytically, even in the massless limit. Instead, Eq. (7.17) gives Q A and σ τ A up to order O(Ω 4 , M 4 , Ω 2 M 2 ) anywhere inside the channel, while Eq. (7.21) gives exact expressions for their values on the rotation axis. The latter expressions involve the polylogarithm function. More insightful expressions can be obtained by considering the high temperature expansions (7.18) and (7.22).
In the case of Q A , according to Eq. (7.18), the mass term makes second order contributions of the form O(M 2 /T 2 ). This is confirmed in Fig. 6(a)  are evaluated using Eq. (7.21), which is valid on the rotation axis for any mass. In the vicinity of the HIC parameters, M 2 /8T 2 ln 2 provides a good approximation for the relative mass corrections. Away from the rotation axis, Eq. (7.17) provides an approximation which is obtained using a small Ω expansion. Eq. (7.17) involves the polylogarithm function, and is thus less insightful compared to its high temperature expansion presented in Eq. (7.18). In what follows, three levels of approximation are considered, denoted using Q A;n (1 ≤ n ≤ 3), which are based on Eq. (7.17). Their validity is tested compared with the exact solution in Eq. (7.21). The first two approximations take into account the leading and next-to-leading order terms in the high temperature expansion, given in Eq. (7.18): The third expression is the massless limit of Eq. (7.17): ) . (7.26) Figures 6(b), 6(c) and 6(d) show the ratio Q A (M )/Q A;n (M ) for n = 1, 2 and 3, respectively. It can be seen in Fig. 6(b) that Q A;1 (0) is a good approximation for Q A on the axis in the case of the HIC parameters, up to M 50 MeV. Increasing the temperature by a factor of 5 (red curve and filled circles) pushes this limit by the same factor. At higher values and of the rotation parameter (Ω 660 MeV) or of the chemical potentials (µ V = µ H = 300 MeV), Q A;1 is no longer a good approximation for Q A (0). Panel (c) shows that Q A;2 is much better at these values, though it still presents a relative error of about 10%. The validity of Q A;2 worsens as Ω and µ are further increased. Finally, panel (d) shows that considering the polylogarithms in Eq. (7.26) is sufficient to correctly account for any value of the chemical potential. The agreement between Q A;3 (0) and Q A (0) is very good even at µ V = µ H = µ = 3 GeV. However, the discrepancies at high Ω persist, since Eq. (7.21) is valid only up to second order in Ω. For all parameters studied, it seems that the variations with respect to the mass appear only at M 100 MeV.
Considering now the properties of σ τ A , Eq. (7.18) indicates that it receives corrections due to the mass term which are proportional to M 2 /Ω 2 . This is confirmed in Fig. 7(a), which shows the relative mass correction 1 − σ τ A (M )/σ τ A (0) with respect to the ratio M/Ω. Surprisingly, the analytic prediction 4M 2 /Ω 2 is dominant for parameters in the vicinity of the HIC values at high values of M/Ω, where the mass term dominates by a few orders of magnitude over the massless limit. This behaviour is also confirmed in Fig. 7(b), where the ratio σ τ A (M )/σ τ A (0) is shown with respect to M . It can be seen that this ratio achieves negative values for the HIC parameters. Furthermore, |σ τ A (M )| increases with respect to the massless prediction σ τ A (0) by a few orders of magnitude as M is increased, decreasing to 0 after reaching a maximum. The amplitude of this maximum decreases when Ω is increased, but it increases when either µ V = µ H = µ or T are increased (the two chemical potentials play an equivalent role). Furthermore, the maximum shifts to higher values of M when the chemical potentials or T are increased.

General analysis
The bilinear form introduced in Eq. (5.17) which corresponds to the FC is: After noting that V j V j = −(U j U j ) * , it is easy to see that the term on the second line of Eq. (5.32) makes a vanishing contribution to the t.e.v. of the FC. Furthermore, the expressions in Eq. (5.2) for the particle modes allow the following result to be obtained: allowing the t.e.v. of the FC to be put in the following form: Taken with respect to the Minkowski vacuum, the t.e.v. of the FC becomes: which differs from the t.e.v. in Eq. (8.3) by a vacuum term.

Small mass limit
Employing the same method as in Subsec. 6.2, the t.e.v. of the FC taken with respect to the Minkowski vacuum can be put in the form: The correction due to the mass cannot be obtained using the methodology from Eq. (5.37), due to an infrared divergence of the j = 0 term. Dividing Eq. (8.5) by M and taking the massless limit yields: where T = T 0 Γ, µ V /H = µ V /h,0 Γ, ω = ΩΓ 2 eẑ and a = ρΩ 2 Γ 2 eρ, while ω 2 = −ω 2 and a 2 = −a 2 . The last term is the vacuum contribution corresponding to the difference between the t.e.v.s taken with respect to the rotation and Minkowski vacua, respectively. This contribution must be subtracted in order to obtain the t.e.v. of the FC with respect to the rotating vacuum: It is surprising that Eq. (8.7) coincides with the RKT prediction (E RKT −3P RKT )/M 2 given in Eq. (4.11).

Numerical analysis
At vanishing mass, a comparison between Eqs. (8.7) and (4.11) indicates that the t.e.v. of the FC (divided by M ) is identical with the RKT prediction for the trace of the SET (divided by M 2 ). At finite mass, it can be expected that this equality holds only approximately. Figure 8 presents the relative difference between the quantum and RKT predictions for the FC, evaluated on the rotation axis. It can be seen that there is a rapid increase form 0 in the massless case to some plateau values. The first plateau occurs quite early, for M 5 MeV. The second plateau appears when M T . The values on these plateaus depend on the parameters of the system. Doubling Ω increases the plateaus by a factor of 4, while doubling the temperature decreases the value, again by a factor of 4, as can be seen from Fig. 8(a). Increasing the fermion condensates (µ V = µ H is considered here) seems to lower the plateau in the central region, however, for M 1 GeV, it gives the same result as in the case of the HIC parameters. The increase from the massless limit to the first plateau seems to follow a power law of the form ∼ M 3/2 , as can be seen from Fig. 8(b). HIC Ω x 2 T x 2 µ x 5

Stress-energy tensor
The stress-energy tensor (SET) operator is defined as:

General analysis
The bilinear form defined in Eq. (5.17) which corresponds to the SET operator is Substituting ψ = χ = V j in the above equation and using the relations in Eq. (5.3), it can be shown that Tασ(U j , U j ) are: The terms appearing in Eq. (9.7) which are odd with respect to k j → −k j do not contribute to the t.e.v. of the SET. It is also interesting to note that, starting from the above expressions, the t.e.v.s : Tρρ : and : Tẑẑ : are exactly equal. This can be seen by first noting that The above relations can be derived using the properties of the derivatives of the Bessel functions J m−1/2 (z) and J m+1/2 (z) given in Eq. (9.6). Next, I m (p) is introduced as follows: Since the integrand is even with respect to k, the integration domain can be replaced by 0 ≤ k ≤ p: Replacing now J − m (qρ) using the two expressions in Eq. (9.8) and using integration by parts to remove the derivative term, it can be shown that Rearranging the above expressions, it can be seen that For a given velocity field uα, the SET can be decomposed as follows: : Tασ : Ω = E uαuσ − (P + ω)∆ασ + Πασ + uαWσ + uσWα, (9.14) where E and P are the usual energy density and pressure, Wα represents the heat flux in the local rest frame, ω is the dynamic pressure and Πασ is the pressure deviator, which is traceless by construction. The tensor ∆ασ = ηασ − uαuσ is a projector on the hypersurface orthogonal to uα. The nonequilibrium quantities Πασ and Wα are also orthogonal to uα, by construction. The isotropic pressure P + ω is given as the sum of the hydrostatic pressure P , which is computed using the equation of state of the fluid, and of the dynamic pressure ω, which in general depenends on the divergence of the velocity. With this convention, the pressure deviator is considered to be traceless. It should be noted that for ultrarelativistic fluids, the SET is traceless (this is true in the quantum case as well since the Dirac field is conformally coupled) and E = 3P . Since Tαα = E − 3(P + ω), it can be seen that ω = 0 for ultrarelativistic (massless) fluid constituents. The macroscopic quantities can be extracted from Tασ as follows: Wα =∆ασuλ : Tσλ : Ω , Πασ = : T ασ : Ω , (9.15) where the notation A ασ for a general two-tensor refers to: The nonequilibrium quantities Wα and Πασ can be further decomposed with respect to the orthogonal tetrad formed by u, a, ω and τ , which is depicted in Fig. 1. Noting that W · u = 0 and that : Tρα : Ω = 0 when α = ρ, the most general decomposition for W is Wα =κ τ τα + κ ω ωα, In the case of Πασ, the orthogonality to uα and the tracelessness condition, together with the consideration that : Tαρ : Ω = 0 for all α = ρ, allow Πασ to be written as: Noting that : Tρρ : Ω = P + Πρρ and : Tẑẑ : Ω = P + Πẑẑ are equal by virtue of Eq. (9.13), it can be seen that The tracelessness condition gives Thus, only two degrees of freedom are required to describe Πασ, which are introduced below: where Π V and Π H can be obtained from the components of the SET, e.g., through: 2n + 2j + 1 2j + 3 p 2 p 2 2j + 3 + E 2 + j 4(ρp) 2 p 2 + E 2 (2j + 1)     , (9.23) while : Tẑẑ : M = : Tρρ : M , as shown in Eq. (9.13). The M subscript indicates that the above expressions are computed with respect to the Minkowski vacuum, as discussed in Sec. 5.3. In obtaining the expression for : Ttφ : M , integration by parts and the relation s + n+1,j = s + n,j + (j + 1 2 ) 2 s + n,j−1 were used. After employing integration by parts 2j times in the integral with respect to p, it is not difficult to see that the summation over n terminates at a finite value of n. This is because ∞ 0 dp p 3 d 2n dp 2n 1 e β 0 (p−µ λ,0 ) + 1 + 1 e β 0 (p+µ λ,0 ) + 1 =  The last equality follows from noting that (1 + e −a ) −1 + (1 + e a ) −1 = 1. Taking into account that the n = 2 term makes a purely vacuum contribution, the following exact results are obtained for the t.e.v.s with respect to the rotating vacuum: : : : There is a set of components which is non-vanishing only when µ H,0 = 0. These can be computed using: As in the previous subsections, the sum over n does not terminate at finite n. Instead, the in agreement with Eq. (9.31). In summary, the SET for the Dirac field can be written as: : Tασ : Ω = (E RKT + ∆E) uαuσ − (P RKT + ∆P )∆ασ + Πασ + uαWσ + uσWα, where in the above, only the leading order terms with respect to T were presented. The corrections due to the mass are of order O(M 4 ). It is worth noting that, in the limit when µ H = M = 0, the results in Eq. (9.37) reproduce the vanishing axial potential limit of Eq. (4.5) and Tables 2 and 3 from Ref. [10].

Numerical analysis
In the beginning of this sectrion, an analysis of the constitutive relations for the quantum corrections ∆E and ∆P is presented. The massless limit results are summarised in Eq. In the case of the energy correction ∆E, a spike can be seen in the case when µ = 3 GeV Ω = 6.6 MeV. This spike is akin to the one observed for the charge correction, ∆Q ± , shown in Fig. 2(c).
Turning to the heat conductivity κ τ along the vector τα, whose massless limit is given in Eq. (9.37), Fig. 10(a) confirms that the leading order mass correction, computed as κ τ (M )/κ τ (0) − 1, is of order O(M 4 ). Fig. 10(b) confirms the validity of the constitutive relation for massess up to ∼ 100 MeV for the HIC values of µ, T and Ω. The spike seen when M → µ at large µ is akin to the one observed for ∆Q ± , σ τ ± and ∆E in Figs. 2(c), 2(d) and 9(c). The peak observed in the curve corresponding to Ω 2 GeV is akin to the one observed for the conductivity σ ω A in Fig. 5(b). The heat conductivity κ ω along the vector ωα is considered in Fig. 11. This quantity can be obtained in closed form on the rotation axis, for any mass, and is given in terms of the polylogarithm functions in Eq. (9.35). First, the leading order contribution coming from the mass is considered. Figure 11(a) presents the relative mass correction, 1−κ ω (M )/κ ω (0), with respect to the ratio M/T . The second order contribution predicted in Eqs. (9.31) and (9.36) is confirmed.
The expression for κ ω derived in Eq. (9.35) is valid only on the rotation axis and is thus inapplicable for studying the properties of κ ω , e.g., in the vicinity of the speed of light surface, as well as the interplay between the vorticity ωẑ = ΩΓ 2 and acceleration aρ = −ρΩΓ 2 , since the latter vanishes on the rotation axis. Another drawback concerns the complexity of the polylogarithm functions, which make the physical properties of κ ω difficult to assess. To this end, three levels of approximations are employed, following the approach taken for σ ω ± (M ) and Q A (M ), which are represented in Figs. 4(b-d) and 6(b-d), respectively. The approximations are denoted using κ ω;n , with 1 ≤ n ≤ 3, and are based on Eq. (9.30). The cases n = 1 and n = 2 are obtained from the high temperature expansion given in Eq. (9.31). For the n = 1 term, only the leading order contribution is taken into account: κ ω;1 = 4µ H µ V T π 2 ln 2. (9.38) As can be seen in Fig. 11(b), this approximation works well for the HIC parameters, breaking down when M T . This is confirmed by looking at the curve corresponding to T 750 MeV (blue line and solid circles). However, the results for higher vorticity (Ω = 660 Mev for the green line and empty circles) and higher chemical potential (µ = 300 Mev for the red line and empty triangles) show a discrepancy.
The n = 2 approximation takes into account the T −1 correction appearing in Eq. (9.31), but disregards the mass term: It should be noted that on the rotation axis, Eq. (9.39) coincides with the massless limit of Eq. (9.36), which represents the high temperature expansion of the exact result in Eq. (9.35). Figure 11(c) shows that this second approximation works well for the two problematic cases discussed in Fig. 11(b), i.e., Ω = 660 MeV and µ = 300 MeV. Further increasing the vorticity (Ω 2 GeV for the blue line with filled upper triangles) or the chemical potential . (9.40) Fig. 11(d) shows that this approximation works well at high chemical potential, however, the discrepancy at high Ω = 2 GeV remains similar to that seen in Fig. 11(c). This discrepancy is due to the fact that Eq. (9.30) represents a small Ω expansion which is valid only up to second order in Ω.
Finally, the quantity Π V is considered. The analysis in the previous section indicated that Π V vanishes up to O(M 4 ). Fig. 12 shows that the leading order contribution to Π V is of order (M/Ω) 6 .
In the case of Π H , the same situation as for κ ω is encountered: Eq. (9.35) gives the exact value of Π H on the rotation axis, for any mass, while Eq. (9.30) gives its value up to O(Ω 2 ) for any value of ρ. The mass is predicted through Eqs. (9.31) and (9.36) to contribute a correction of order O(M 2 ). This is confirmed in Fig. 13(a). Next, the high temperature approximation given in Eq. (9.31) is considered, which is labelled according to the convention in this section as Π H;1 : As demonstrated in Fig. 13(b), this approximation is valid for the HIC conditions, but loses applicability when either Ω or µ are increased. The second approximation is the one derived in Eq. (9.30), which is valid up to O(Ω 2 ): . 13(c) shows that this approximation is valid at high values of the chemical potential (µ = 0.9 GeV is shown with black lines and lower, filled triangles), but remains inaccurate at high vorticities. This is because Eq. (9.30) is derived as an O(Ω 4 ) expansion, which can be expected to be inaccurate at high values of Ω.

Conclusions
In this paper, the properties of massive Dirac fermions under rotation at finite temperature and finite chemical potential were considered. Aside from the chemical potential associated to the vector charge, a helicity chemical potential was introduced, with the help of which the polarisation imbalance of massive fermions was accounted for. A simple relativistic kinetic theory model incorporating the helicity imbalance was proposed to serve as a background classical theory. Using quantum field theory at finite temperature, analytic expressions valid for small masses were derived and the constitutive equations that describe the quantum corrections to the classical quantities (energy density, pressure and charge densities) were highlighted. The various transport phenomena driven by the rotation was discussed. The results presented here reduce in the limit of vanishing mass and vanishing helicity chemical potential to the vanishing axial chemical potential limit of the results reported in Ref. [10]. It was shown that, under the conditions which are prevalent in the quark-gluon plasma formed in relativistic heavy ion collisions experiments (the HIC conditions), the constitutive relations are in general robust, in the sense that they remain valid for masses up to the thermal energy (M 150 MeV). The only exception was seen in the case of the axial current conductivity σ τ A along the τ µ vector (which is dual to the velocity four-velocity and lies in the t − ϕ plane), which exhibits a very strong mass dependence even at small values of the mass.
where the explicit dependence of I ± n on a was dropped for brevity. The second terms appearing above can be integrated, yielding: Noting that 1 e x + 1 + 1 e −x + 1 = 1 2 , (A. 6) it can be seen that (e x + 1) −1 admits the following series representation about x = 0: where the exact expression for a is not important. 2 The only even term in the series (A.7) is the leading 1/2, while all other terms are odd with respect to x. Thus, Eq. (A.5) is non-vanishing only when j = 0, such that: Aside from the integrals I ± n (a) considered in Eq. (A.1), the analysis in the main text requires the computation of integrals when the subscript of I is no longer an integer. Retaining the convention that n = 0, 1, 2, . . . is a natural number, the integrals I ± n+1/2 (a) can be defined by analogy to Eq. (A.1) as follows: where Li n (z) = ∞ k=1 z k /n k is the polylogarithm. The following relations were employed: