11institutetext: Matteo Buzzegoli 22institutetext: Universitá di Firenze and
INFN Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino (Firenze),
Italy, 22email<EMAIL_ADDRESS>
# Thermodynamic equilibrium of massless fermions with vorticity, chirality and
electromagnetic field

Matteo Buzzegoli

###### Abstract

We present a study of the thermodynamics of the massless free Dirac field at
equilibrium with axial charge, angular momentum and external electromagnetic
field to assess the interplay between chirality, vorticity and electromagnetic
field in relativistic fluids. After discussing the general features of global
thermodynamic equilibrium in quantum relativistic statistical mechanics, we
calculate the thermal expectation values. Axial imbalance and electromagnetic
field are included non-perturbatively by using the exact solutions of the
Dirac equation, while a perturbative expansion is carried out in thermal
vorticity. It is shown that the chiral vortical effect and the axial vortical
effect are not affected by a constant homogeneous electromagnetic field.

##  1 Introduction

The collective macroscopic behavior of matter in the presence of quantum
anomalies and external fields is an increasingly important subject in several
fields of physics. Specifically, the experiments of relativistic heavy-ion
collisions at RHIC and LHC have posed new and interesting questions about the
theoretical foundations of relativistic collective phenomena. The experimental
data of heavy ion collisions indicates the creation of a deconfined quark
gluon plasma in a strongly coupled regime at extreme conditions of
temperature, density, thermal vorticity STAR:2017ckg  and magnetic fields
Kharzeev:2015znc . Moreover, it was argued Kharzeev:2007jp ; Fukushima:2008xe
that the fluctuations of topological configurations of the QCD vacuum in the
early stages of a heavy ion collision generate a chiral imbalance, which is an
imbalance between the number of right- and left- handed quarks. Despite the
usual relativistic hydrodynamic has been very effective Heinz:2013th  in
reproducing the experimental data for collective flow phenomena, it is now
essential for the interpretation of heavy ion collisions to address
hydrodynamics in the contemporaneous presence of chiral imbalance, thermal
vorticity, and external electromagnetic fields.

The first crucial step towards understanding the hydrodynamics of matter
subject to external fields is to study its thermodynamic properties. It is the
main purpose of this contribution to investigate the effects of an external
electromagnetic field on the thermodynamics of a chiral vorticous fluid. The
effects of electromagnetic fields on (non-chiral non-vorticous) relativistic
quantum fluids were already studied in the past, see for instance HakimBook ;
Israel:1978up  and reference therein and Canuto:1969cs  for the special case
of a constant magnetic field. More recently this topic has been addressed in
Huang:2011dc  using the Zubarev’s non-equilibrium statistical operator, in
Kovtun:2016lfw  using the generating functional method and in
Weickgenannt:2019dks ; Gao:2019znl ; Chen:2012ca  with Wigner function derived
from kinetic theory.

This contribution aims to highlight the modifications caused by chiral
imbalance and by thermal vorticity. The paper is organized in the following
way. In Sec. 2, we introduce the global thermal equilibrium of a chiral system
with the contemporaneous presence of an external electromagnetic field and a
thermal vorticity within Zubarev’s non-equilibrium statistical operator
formalism. In Sec. 2.1 we give a brief overview of the main results for the
case of a chiral Dirac field in the absence of the electromagnetic field. In
Sec. 3 we review the relativistic quantum theory of fermions under the effect
of an external magnetic field. Then, we obtain the exact form of the chiral
fermionic propagator with external constant magnetic field and we obtain the
exact thermal averages of the axial and electric currents. In Sec. 5 we
examine the properties of a system at thermal equilibrium with constant
vorticity and electromagnetic field. The last part of the paper is concerned
with the consequences of an electromagnetic field on the chiral vortical
effect and the axial vortical effect.

## Notation

In this work we use the _natural unit_ system in which $\hslash=c=G=k_{B}=1$.
The _Minkowski metric_ is defined by the tensor $\eta_{\mu\nu}=$
diag$(1,-1,-1,-1)$; for the Levi-Civita symbol we use the convention
$\epsilon^{0123}=+1$.

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
$\Psi$). The stress-energy tensor used to define Poincaré generators is always
assumed to be symmetric with an associated vanishing spin tensor.

##  2 General global equilibrium with electromagnetic field

In this section we introduce the methods to study the thermodynamic
equilibrium of a quantum relativistic system in the presence of a chiral
imbalance and of an external electromagnetic field. For that purpose we review
the Zubarev method of stationary non-equilibrium density operator Zubarev1979
; vanWeert1982  (see also BecaBetaF ; Hayata:2015lga ; Hongo:2016mqm ;
Becattini:2019dxo  for recent developments) and we discuss the inclusion of a
conserved axial current and of an external electromagnetic field.

When we are dealing with a relativistic system we must consider local
quantities in order to address the appropriate covariant properties. To
identify those quantities we use a Arnowitt-Deser-Misner (ADM) decomposition
of space-time Zubarev1979 ; vanWeert1982 ; BecaBetaF ; Hayata:2015lga . Choose
then a foliation of space-time and suppose that the system in consideration
thermalize faster than the evolution of “time” $\tau$ in which we are
interested. At each step of evolution $\mathrm{d}\tau$, the system is at local
thermal equilibrium and the macroscopic behavior of the system is described by
a stress-energy density $T_{\mu\nu}(x)$, an (electric) current density
$j_{\mu}(x)$ and an axial current ${j_{\textup{A}}}_{\mu}(x)$, all lying on a
space-like hyper-surface $\Sigma(\tau)$. We can then describe the thermal
properties of the system with a density operator which lives on
$\Sigma(\tau)$. As in the non-relativistic case, the density operator at local
equilibrium $\widehat{\rho}_{\textup{LE}}$ is obtained as the operator which
maximizes the entropy $S=-{\text{tr}}(\,\widehat{\rho}\log\widehat{\rho}\,)$.
To reproduce the actual thermodynamics on the hyper-surface we maximize the
entropy with the constraints that the mean values of the stress-energy tensor
and of the currents on $\Sigma(\tau)$ corresponds to the values of the
densities $T_{\mu\nu}(x)$, $j_{\mu}(x)$ and ${j_{\textup{A}}}_{\mu}(x)$
vanWeert1982 . To obtain these densities we project the stress–energy tensor
and the current mean values onto $n$, i.e. the normalized four-vector
perpendicular to $\Sigma$:

$\begin{split}n_{\mu}(x)\,{\text{tr}}\left[\,\widehat{\rho}\,\,\widehat{T}^{\mu\nu}(x)\right]=&\,n_{\mu}(x)\,\langle\,\widehat{T}^{\mu\nu}(x)\rangle\equiv
n_{\mu}(x)\,T^{\mu\nu}(x),\\\
n_{\mu}(x)\,{\text{tr}}\left[\,\widehat{\rho}\,\,\widehat{j}^{\,\mu}(x)\right]=&\,n_{\mu}(x)\,\langle\widehat{j}^{\,\mu}(x)\rangle\equiv
n_{\mu}(x)\,j^{\mu}(x),\\\ \end{split}$

and similarly for the axial current. We could also impose a constraint on the
angular momentum density, but since we are choosing the Belinfante operator as
the stress–energy tensor, it turns out that this additional requirement is
automatically taken into account Becattini:2018duy .

The maximum solution $\widehat{\rho}_{\textup{LE}}$ gives the _Local
Equilibrium Density Operator_ (LEDO) BecaBetaF ; Hayata:2015lga :

$\widehat{\rho}_{\textup{LTE}}=\frac{1}{Z}\exp\left[-\int_{\Sigma}\mathrm{d}\Sigma_{\mu}\left(\widehat{T}^{\mu\nu}(x)\beta_{\nu}(x)-\zeta(x)\,\widehat{j}^{\,\mu}(x)-\zeta_{\textup{A}}(x)\,\widehat{j}_{\textup{A}}^{\,\mu}(x)\right)\right],$
(1)

where $\beta^{\mu}$ is the four-temperature vector such that
$T=1/\sqrt{\beta^{2}}$ is the proper comoving temperature, $\zeta$ and
$\zeta_{\textup{A}}$ are the ratio of comoving chemical potentials and the
temperature (e.g. $\zeta=\mu/T$) and $Z$ is the partition function. In the
presence of an external electromagnetic field, we indicate with $A^{\mu}(x)$
the non-dynamical gauge field and with
$F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}$ the electromagnetic
strength tensor. Therefore, the operator relations stemming for conservation
equations are

$\partial_{\mu}\widehat{j}^{\,\mu}=0,\quad\partial_{\mu}\widehat{T}^{\mu\nu}=\widehat{j}_{\lambda}F^{\nu\lambda},\quad\partial_{\mu}\widehat{j}^{\,\mu}_{\textup{A}}=0.$
(2)

Furthermore, the four-momentum operator $\widehat{P}$ and the conserved
charges $\widehat{Q}_{i}$ are obtained by

$\widehat{P}^{\mu}=\int_{\Sigma}\mathrm{d}\Sigma_{\lambda}\,\widehat{T}^{\lambda\mu},\qquad\widehat{Q}_{i}=\int_{\Sigma}\mathrm{d}\Sigma_{\lambda}\,\widehat{j}^{\,\lambda}_{i},$

while the angular momentum is

$\widehat{J}^{\mu\nu}=\int_{\Sigma}\mathrm{d}\Sigma_{\lambda}\left(x^{\mu}\widehat{T}^{\lambda\nu}-x^{\nu}\widehat{T}^{\lambda\mu}\right).$
(3)

Notice that, as we will discuss in details in Sec. 3, the four-momentum
$\widehat{P}$ and the angular momentum $\widehat{J}$ in the presence of
external electromagnetic field are neither conserved nor the generators of
translations and Lorentz transformations.

In the case of Dirac fermions interacting with an external gauge field, the
explicit form of the operators above is

$\begin{split}\widehat{j}^{\,\mu}=&\,q\bar{\Psi}\gamma^{\mu}\Psi,\quad\widehat{j}^{\,\mu}_{\textup{A}}=\,\bar{\Psi}\gamma^{\mu}\gamma^{5}\Psi,\\\
\widehat{T}^{\mu\nu}=&\frac{{\rm
i}}{4}\left[\bar{\Psi}\gamma^{\mu}\overrightarrow{\partial}^{\nu}\Psi-\bar{\Psi}\gamma^{\mu}\overleftarrow{\partial}^{\nu}\Psi+\bar{\Psi}\gamma^{\nu}\overrightarrow{\partial}^{\mu}\Psi-\bar{\Psi}\gamma^{\nu}\overleftarrow{\partial}^{\mu}\Psi\right]-\frac{1}{2}\left(\widehat{j}^{\,\mu}A^{\nu}+\widehat{j}^{\nu}A^{\mu}\right),\end{split}$
(4)

and the stress-energy tensor and the electric current indeed satisfy the
relations (2). Regarding the axial current $\widehat{j}_{\textup{A}}$, we also
have to take into account the chiral anomaly. The chiral anomaly affects the
axial current divergence as follows

$\partial_{\mu}\widehat{j}^{\,\mu}_{\textup{A}}=-\frac{1}{8}\epsilon^{\mu\nu\rho\lambda}\frac{q^{2}}{2\pi^{2}}F_{\mu\nu}F_{\rho\lambda}=-\frac{q^{2}}{2\pi^{2}}(E\cdot
B),$

where $q$ is the electric charge of the fermion, and $E$ and $B$ are comoving
electric and magnetic field, defined by

$F^{\mu\nu}=E^{\mu}u^{\nu}-E^{\nu}u^{\mu}-\epsilon^{\mu\nu\rho\sigma}B_{\rho}u_{\sigma},$

with $u$ the fluid velocity. Even when the product $E\cdot B$ is non-vanishing
and consequently the axial current is not conserved, we can still define a new
conserved “axial” current by means of the Chern-Simons current $K$, whose
divergence gives the chiral anomaly:

$K^{\mu}=\epsilon^{\mu\nu\rho\sigma}A_{\nu}F_{\rho\sigma},\quad\frac{q^{2}}{8\pi^{2}}\partial_{\mu}K^{\mu}=\frac{q^{2}}{2\pi^{2}}(E\cdot
B).$

The new conserved axial current $\widehat{j}_{\textup{CS}}$ is then defined as

$\widehat{j}^{\,\mu}_{\textup{CS}}=\widehat{j}^{\,\mu}_{\textup{A}}+\frac{q^{2}}{8\pi^{2}}K^{\mu},\quad\partial_{\mu}\widehat{j}^{\,\mu}_{\textup{CS}}=0,$

and the axial chemical potential $\mu_{\textup{A}}$ is to be associated to
this current. Since the additional current $K$ depends only on external
fields, it is not a quantum operator and it does not contribute to thermal
averages. Therefore, all the results discussed in the absence of chiral
anomaly will also be valid for the case of equilibrium with conserved Chern-
Simons current. Because there is no difference in the results, we will
continue to denote the current associated to $\mu_{\textup{A}}$ inside the
statistical operator with $\widehat{j}_{\textup{A}}$ even when the chiral
anomaly is non-vanishing.

Let us now move to describe the system at global thermal equilibrium. The
global equilibrium is reached when the statistical operator (1) is time
independent. This occurs when the integrand inside Eq. (1) is divergenceless
BecCov . Then, it is easily proven using relations (2) that global thermal
equilibrium is realized when the following relations are satisfied:

$\partial_{\mu}\beta_{\nu}(x)+\partial_{\nu}\beta_{\mu}(x)=0,\qquad\partial^{\mu}\zeta(x)=F^{\nu\mu}\beta_{\nu}(x),\qquad\partial^{\mu}\zeta_{\textup{A}}(x)=0.$
(5)

The inverse four-temperature and the axial chemical potential solves the
previous conditions if they are given by BecaBetaF :

$\beta_{\mu}(x)=b_{\mu}+\varpi_{\mu\rho}x^{\rho},\quad\zeta_{\textup{A}}=\text{constant},$

where $b$ is a constant time-like four-vector and $\varpi$ is a constant anti-
symmetric tensor. We refer to $\varpi$ as the _thermal vorticity_ because it
is the anti-symmetric derivative of inverse four-temperature:

$\varpi_{\mu\nu}=-\frac{1}{2}\left(\partial_{\mu}\beta_{\nu}-\partial_{\nu}\beta_{\mu}\right)$

and because it contains information about the fluid’s acceleration and
rotation. Indeed, if the $\beta$ four-vector is a time-like vector, then we
can choose the $\beta$-frame as hydrodynamic frame BecCov ; Kovtun:2016lfw .
The unitary four-vector fluid velocity $u$ is therefore identified with the
direction of $\beta$:

$u^{\mu}(x)=\frac{\beta^{\mu}(x)}{\sqrt{\beta_{\rho}(x)\beta^{\rho}(x)}}.$

As long as we are considering physical observables in a region where the
coordinate $x$ are such that $\beta(x)$ is a time-like vector, this definition
provides a proper choice for the fluid velocity. We can decompose the thermal
vorticity into two space-like vector fields, each having three independent
components, by projecting along the time-like fluid velocity $u$:

$\varpi^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}w_{\rho}u_{\sigma}+\alpha^{\mu}u^{\nu}-\alpha^{\nu}u^{\mu}.$

The four-vectors $\alpha$ and $w$ are explicitly written inverting the
previous relation:

$\alpha^{\mu}(x)\equiv\varpi^{\mu\nu}u_{\nu},\quad
w^{\mu}(x)\equiv-\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}\varpi_{\nu\rho}u_{\sigma}.$

The vectors $\alpha$ and $w$ depend on the coordinates, are space-like, and
are orthogonal to $u$. All the quantity $u,\,\varpi,\,\alpha,\,w$ are
dimensionless. From their definitions, we can easily show that $\alpha$ and
$w$ are given by

$\alpha^{\mu}=\sqrt{\beta^{2}}\,a^{\mu},\quad
w^{\mu}=\sqrt{\beta^{2}}\,\omega^{\mu},$

where $a$ and $\omega$ are the local acceleration and rotation of the fluid,
which are given by

$a^{\mu}=u_{\nu}\partial^{\nu}u^{\mu},\quad\omega^{\mu}\equiv\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}u_{\sigma}\partial_{\nu}u_{\rho}.$

Furthermore, it will proves useful to define the projector into the orthogonal
space of fluid velocity:

$\Delta_{\mu\nu}\equiv g_{\mu\nu}-u_{\mu}u_{\nu},$

and the four vector $\gamma$ orthogonal to the other ones: $u,\alpha,w$

$\gamma^{\mu}=\epsilon^{\mu\nu\rho\sigma}w_{\nu}\alpha_{\rho}u_{\sigma}=(\alpha\cdot\varpi)_{\lambda}\Delta^{\lambda\mu}.$

The $\varpi$ decomposition above defines a tetrad $\\{u,\alpha,w,\gamma\\}$
which can be used as a basis for four-vectors. It must be noticed, however,
that the tetrad is neither unitary nor orthonormal, indeed in general we have
$\alpha\cdot w\neq 0$.

Returning to the global equilibrium conditions (5), notice that in the absence
of electromagnetic field also $\zeta$ must be a constant. In that case the
global equilibrium statistical operator takes the following form BecGro ;
Buzzegoli:2018wpy :

$\widehat{\rho}=\frac{1}{Z}\exp\left\\{-b\cdot\widehat{P}+\frac{1}{2}\varpi:\widehat{J}+\zeta\widehat{Q}+\zeta_{\textup{A}}\widehat{Q}_{\textup{A}}\right\\}.$
(6)

The thermodynamics of Dirac fermions which follows from this operator is
quickly reviewed in Section 2.1

For the case of a non-vanishing electromagnetic field instead, we need to
solve the equation:

$\partial^{\mu}\zeta(x)=F^{\sigma\mu}\beta_{\sigma}.$ (7)

To find the solution, we first derive it respect to $\partial^{\nu}$:

$\partial^{\nu}\partial^{\mu}\zeta=\partial^{\nu}(F^{\sigma\mu}\beta_{\sigma}).$
(8)

Since we can exchange the order of the derivatives
$\partial^{\nu}\partial^{\mu}$ on the l.h.s. of (8), it follows that the anti-
symmetrization respect to indices $\mu$ and $\nu$ of (8) must be vanishing:

$\begin{split}\partial^{\nu}\partial^{\mu}\zeta-\partial^{\mu}\partial^{\nu}\zeta&=0=\left[\partial^{\nu}(F^{\sigma\mu}\beta_{\sigma})-\partial^{\mu}(F^{\sigma\nu}\beta_{\sigma})\right]\\\
&=\left[\beta_{\sigma}(\partial^{\nu}F^{\sigma\mu}-\partial^{\mu}F^{\sigma\nu})+(\partial^{\nu}\beta_{\sigma})F^{\sigma\mu}+(\partial^{\mu}\beta_{\sigma})F^{\nu\sigma}\right].\end{split}$

Using the first Bianchi identity
$\partial^{\nu}F^{\sigma\mu}+\partial^{\mu}F^{\nu\sigma}+\partial^{\sigma}F^{\mu\nu}=0$,
we obtain

$\beta_{\sigma}\partial^{\sigma}F^{\mu\nu}+(\partial^{\nu}\beta_{\sigma})F^{\mu\sigma}+(\partial^{\mu}\beta_{\sigma})F^{\sigma\nu}=0.$

We may recognize the Lie derivative of $F$ along $\beta$ in the previous
equation. This constitutes a first condition for global equilibrium, the
system can reach global equilibrium only when

$\mathcal{L}_{\beta}(F)=0,\quad\leftrightarrow\quad\beta_{\sigma}(x)\partial^{\sigma}F^{\mu\nu}(x)=\varpi^{\mu}_{\hphantom{\mu}\sigma}F^{\sigma\nu}(x)-\varpi^{\nu}_{\hphantom{\nu}\sigma}F^{\sigma\mu}(x),$
(9)

that is to say when the electromagnetic field follows the field lines of
inverse-four temperature.

To actually solve Eq. (7) we translate the global equilibrium condition of the
strength tensor (9) to the four-vector potential $A^{\mu}$. We see that the
constraint (9) is satisfied if $A$ solves

$\beta_{\sigma}(x)\partial^{\sigma}A^{\mu}(x)=\varpi^{\mu}_{\hphantom{\mu}\sigma}A^{\sigma}(x)+\partial^{\mu}\Phi(x),$
(10)

where $\Phi$ is a smooth function of $x$. In ref. Hongo:2016mqm  was also
stated that a gauge potential with vanishing Lie derivative along $\beta$
gives a stationary statistical operator, which is condition (10). It is
important to stress that after a gauge transformation the condition (10) still
holds true for the new gauge potential because the function $\Phi$ is also
affected by the gauge transformation. Indeed, let $A^{\mu}$ satisfies Eq.
(10); after the gauge transformation
$A^{\prime\mu}=A^{\mu}+\partial^{\mu}\Lambda$ we find:

$\begin{split}\beta_{\sigma}\partial^{\sigma}A^{\prime\mu}=&\beta_{\sigma}\partial^{\sigma}A^{\mu}+\beta_{\sigma}\partial^{\sigma}\partial^{\mu}\Lambda=\omega^{\mu}_{\,\,\sigma}A^{\sigma}+\partial^{\mu}\Phi+\partial^{\mu}(\beta_{\sigma}\partial^{\sigma}\Lambda)-(\partial^{\mu}\beta_{\sigma})\partial^{\sigma}\Lambda\\\
=&\varpi^{\mu}_{\hphantom{\mu}\sigma}(A^{\sigma}+\partial^{\sigma}\Lambda)+\partial^{\mu}(\Phi+\beta_{\sigma}\partial^{\sigma}\Lambda)=\varpi^{\mu}_{\hphantom{\mu}\sigma}A^{\prime\sigma}+\partial^{\mu}\Phi^{\prime},\end{split}$

which is exactly condition (10) for $A^{\prime\mu}$ and for $\Phi^{\prime}$,
that is $\Phi$ shifted by the transport of $\Lambda$ along $\beta$.

We can now write Eq. (7) by taking advantage of Eq. (10):

$\begin{split}\partial^{\mu}\zeta&=F^{\sigma\mu}\beta_{\sigma}=\beta_{\sigma}(\partial^{\sigma}A^{\mu}-\partial^{\mu}A^{\sigma})=\beta_{\sigma}\partial^{\sigma}A^{\mu}-\partial^{\mu}(\beta_{\sigma}A^{\sigma})+(\partial^{\mu}\beta_{\sigma})A^{\sigma}\\\
&=\varpi^{\mu}_{\,\sigma}A^{\sigma}+\partial^{\mu}\Phi-\varpi^{\mu}_{\,\sigma}A^{\sigma}-\partial^{\mu}(\beta_{\sigma}A^{\sigma}).\end{split}$

We can then collect all the derivatives together into the equation

$\partial^{\mu}\left(\zeta-\Phi+\beta_{\sigma}A^{\sigma}\right)=0,$

from which we immediately get the solution:

$\zeta(x)=\zeta_{0}-\beta_{\sigma}(x)A^{\sigma}(x)+\Phi(x),$ (11)

where $\zeta_{0}$ is a constant. The parameter $\Phi$ is analogous to the
parameter which grants gauge invariance to chemical potential in
Jensen:2013kka . Even though Eq. (11) is given in terms of the gauge potential
it is still gauge invariant. Indeed, we have shown that with a gauge
transformation, $A^{\mu}$ and $\Phi$ transform as

$A^{\prime\mu}=A^{\mu}+\partial^{\mu}\Lambda,\quad\Phi^{\prime}=\Phi+\beta_{\sigma}\partial^{\sigma}\Lambda,$

therefore, the chemical potential $\zeta$ is overall unaffected by gauge
transformations:

$\zeta(x)^{\prime}=\zeta_{0}-\beta_{\sigma}(x)A^{\prime\sigma}+\Phi^{\prime}=\zeta_{0}-\beta_{\sigma}(x)A^{\sigma}+\Phi-\beta_{\sigma}\partial^{\sigma}\Lambda+\beta_{\sigma}\partial^{\sigma}\Lambda=\zeta(x).$

The global equilibrium statistical operator is then obtained from the local
one in Eq. (1) by replacing the global equilibrium form of the thermodynamic
fields $\beta,\zeta,\zeta_{\textup{A}}$:

$\begin{split}\widehat{\rho}=&\frac{1}{Z}\exp\left\\{-\int_{\Sigma}\mathrm{d}\Sigma_{\mu}\left[\left(\widehat{T}^{\mu\nu}(x)+\widehat{j}^{\,\mu}(x)A^{\nu}(x)\right)\beta_{\nu}(x)\right.\right.\\\
&\left.\left.-\left(\zeta_{0}+\Phi(x)\right)\widehat{j}^{\,\mu}(x)-\zeta_{\textup{A}}\,\widehat{j}_{\textup{A}}^{\,\mu}(x)\right]\right\\}.\end{split}$
(12)

This operator, on par with (6), is given in terms of global conserved
quantities. The difference is that for the general form of the external
magnetic field satisfying the constraint (9), the integration over the hyper-
surface $\Sigma$ does not give easily recognizable quantities like the four-
momentum and the angular momenta in Eq. (6). However, the identification of
global conserved operators can be carried out in the special case of constant
homogeneous electromagnetic field and it is discussed in the following
sections.

###  2.1 Vanishing electromagnetic field

Before proceeding with the effects of electromagnetic fields, we briefly
review the thermodynamics properties of a relativistic system in the presence
of thermal vorticity but without an external electromagnetic field. Regarding
thermal equilibrium in the presence of rotation, exact solutions for the free
scalar and Dirac fields are discussed in Ambrus:2014uqa ; Ambrus:2019cvr ;
Becattini:2020qol . Instead the effects of acceleration has been recently
investigated using the Zubarev method in Ref. Becattini:2020qol ;
Becattini:2017ljh ; Becattini:2019poj ; Prokhorov:2019sss ; Prokhorov:2019cik
; Prokhorov:2019yft . Here we want to report the constitutive equations at
second order on thermal vorticity discussed in BecGro ; Buzzegoli:2017cqy  and
in Buzzegoli:2018wpy  including an axial current (see also Hongo:2019rbd  for
first order in thermal vorticity and magnetic field). Using linear response
theory on thermal vorticity, the thermal expectation value of a local operator
$\widehat{O}(x)$ evaluated with statistical operator (9) can be written as
BecGro ; Buzzegoli:2018wpy :

$\begin{split}\langle\widehat{O}(x)\rangle=&\langle\widehat{O}(0)\rangle_{\beta(x)}-\alpha_{\rho}\langle\langle\,\widehat{K}^{\rho}\widehat{O}\,\rangle\rangle-
w_{\rho}\langle\langle\,\widehat{J}^{\rho}\widehat{O}\,\rangle\rangle+\frac{\alpha_{\rho}\alpha_{\sigma}}{2}\langle\langle\,\widehat{K}^{\rho}\widehat{K}^{\sigma}\widehat{O}\,\rangle\rangle\\\
&+\frac{w_{\rho}w_{\sigma}}{2}\langle\langle\,\widehat{J}^{\rho}\widehat{J}^{\sigma}\widehat{O}\,\rangle\rangle+\frac{\alpha_{\rho}w_{\sigma}}{2}\langle\langle\,\\{\widehat{K}^{\rho},\widehat{J}^{\sigma}\\}\widehat{O}\,\rangle\rangle+\mathcal{O}(\varpi^{3}).\end{split}$
(13)

In the previous expression we indicated with double angular bracket the
correlator

$\begin{split}\langle\langle\widehat{K}^{\rho_{1}}\cdots\widehat{K}^{\rho_{n}}\widehat{J}^{\sigma_{1}}\cdots\widehat{J}^{\sigma_{m}}&\widehat{O}\rangle\rangle\equiv\int_{0}^{|\beta|}\frac{\mathrm{d}\tau_{1}\cdots\mathrm{d}\tau_{n+m}}{|\beta|^{n+m}}\times\\\
&\times\langle{\rm T}_{\tau}\left(\widehat{K}^{\rho_{1}}_{-{\rm
i}\tau_{1}u}\cdots\widehat{K}^{\rho_{n}}_{-{\rm
i}\tau_{n}u}\widehat{J}^{\sigma_{1}}_{-{\rm
i}\tau_{n+1}u}\cdots\widehat{J}^{\sigma_{m}}_{-{\rm
i}\tau_{n+m}u}\widehat{O}(0)\right)\rangle_{\beta(x),c},\end{split}$

where $\widehat{J}$ and $\widehat{K}$ are the comoving rotation and boost
generators, identified by

$\widehat{K}^{\mu}=u_{\lambda}\widehat{J}^{\lambda\mu},\quad\widehat{J}^{\mu}=\frac{1}{2}\epsilon^{\alpha\beta\gamma\mu}u_{\alpha}\widehat{J}_{\beta\gamma},$

and the averages $\langle\cdots\rangle_{\beta(x)}$ are evaluated at a fixed
point $x$ with the homogeneous statistical operator

$\widehat{\rho}_{0}=\frac{1}{Z_{0}}\exp\left\\{-\beta(x)\cdot\widehat{P}+\zeta(x)\widehat{Q}+\zeta_{\textup{A}}(x)\widehat{Q}_{\textup{A}}\right\\}.$

The subscript $-{\rm i}\tau u$ indicates an imaginary translation along $u$ as
follows

$J^{\mu}_{-{\rm i}\tau u}\equiv\mathrm{e}^{-{\rm i}\tau
u\cdot\widehat{P}}\widehat{J}^{\mu}\mathrm{e}^{{\rm i}\tau
u\cdot\widehat{P}}.$

Constitutive equations at second order on thermal vorticity of the stress-
energy tensor, the electric current and the axial current can be obtained
using the expansion in Eq. (13). We obtain BecGro ; Buzzegoli:2017cqy ;
Buzzegoli:2018wpy

$\begin{split}\langle\widehat{T}^{\mu\nu}\rangle=&\,\mathbb{A}\,\epsilon^{\mu\nu\kappa\lambda}\alpha_{\kappa}u_{\lambda}+\mathbb{W}_{1}w^{\mu}u^{\nu}+\mathbb{W}_{2}w^{\nu}u^{\mu}\\\
&+(\rho-\alpha^{2}U_{\alpha}-w^{2}U_{w})u^{\mu}u^{\nu}-(p-\alpha^{2}D_{\alpha}-w^{2}D_{w})\Delta^{\mu\nu}\\\
&+A\,\alpha^{\mu}\alpha^{\nu}+Ww^{\mu}w^{\nu}+G_{1}u^{\mu}\gamma^{\nu}+G_{2}u^{\nu}\gamma^{\mu}+\mathcal{O}(\varpi^{3}),\end{split}$
(14)
$\langle\widehat{j}_{\textup{V}}^{\mu}\rangle=n_{\textup{V}}\,u^{\mu}+\left(\alpha^{2}N^{\textup{V}}_{\alpha}+w^{2}N^{\textup{V}}_{\omega}\right)u^{\mu}+W^{\textup{V}}w^{\mu}+G^{\textup{V}}\gamma^{\mu}+\mathcal{O}(\varpi^{3}),$
(15)
$\langle\widehat{j}_{\textup{A}}^{\mu}\rangle=n_{\textup{A}}\,u^{\mu}+\left(\alpha^{2}N_{\alpha}^{\textup{A}}+w^{2}N_{\omega}^{\textup{A}}\right)u^{\mu}+W^{\textup{A}}w^{\mu}+G^{\textup{A}}\gamma^{\mu}+\mathcal{O}(\varpi^{3}).$
(16)

Not all of these coefficients are independent, indeed conservation equations
(2) impose the following relations Buzzegoli:2017cqy  (this is explained in
details in Sec. 5.2):

$\begin{split}U_{\alpha}&=-|\beta|\frac{\partial}{\partial|\beta|}\big{(}D_{\alpha}+A\big{)}-\big{(}D_{\alpha}+A\big{)},\\\
U_{w}&=-|\beta|\frac{\partial}{\partial|\beta|}\big{(}D_{w}+W\big{)}-D_{w}+2A-3W,\\\
G_{1}+G_{2}&=2\big{(}D_{\alpha}+D_{w}\big{)}+A+|\beta|\frac{\partial}{\partial|\beta|}W+3W,\end{split}$

instead, for the first-order coefficients, conservation equations require that

$-2\mathbb{A}=|\beta|\frac{\partial\mathbb{W}_{1}}{\partial|\beta|}+3\mathbb{W}_{1}+\mathbb{W}_{2}.$

For electric and axial current, we find that only the following equations must
be fulfilled:

$|\beta|\frac{\partial
W^{\textup{V}}}{\partial|\beta|}+3W^{\textup{V}}=0,\quad|\beta|\frac{\partial
W^{\textup{A}}}{\partial|\beta|}+3W^{\textup{A}}=0.$ (17)

We can also take advantage of the Lorentz symmetry to show that the thermal
coefficients $\mathbb{S}$ and $\Gamma_{w}$ of the canonical spin tenor
constitutive equation:

$\langle\frac{{\rm
i}}{8}\bar{\Psi}\,\\{\gamma^{\lambda},\,[\gamma^{\mu},\,\gamma^{\nu}]\,\\}\,\Psi\rangle=\mathbb{S}\,\epsilon^{\lambda\mu\nu\rho}u_{\rho}+\Gamma_{w}\,\left(u^{\lambda}\varpi^{\mu\nu}+u^{\nu}\varpi^{\lambda\mu}+u^{\mu}\varpi^{\nu\lambda}\right)+\mathcal{O}(\varpi^{2}),$

satisfy the following relations

$\begin{split}-\left(\frac{\mathbb{S}}{|\beta|}+\frac{\partial\mathbb{S}}{\partial|\beta|}\right)&=2\mathbb{A}^{\textup{Can}},\\\
2\frac{\mathbb{S}}{|\beta|}&=\mathbb{W}_{1}^{\textup{Can}}-\mathbb{W}_{2}^{\textup{Can}},\\\
\frac{\Gamma_{w}}{|\beta|}-\frac{\partial\Gamma_{w}}{\partial|\beta|}&=4\frac{\Gamma_{w}}{|\beta|}=G_{1}^{\textup{Can}}-G_{2}^{\textup{Can}},\end{split}$
(18)

where $\mathbb{A}^{\textup{Can}}$ and $\mathbb{W}_{1,2}^{\textup{Can}}$ are
the thermal coefficients of Eq. (14) related to the mean value of _canonical_
stress-energy tensor. Furthermore, because the axial current is dual to the
spin tensor, we can show that

$\mathbb{S}=\frac{1}{2}n_{\textup{A}},\quad\Gamma_{w}=\frac{1}{2}W^{\textup{A}}.$
(19)

Then, combining Eq.s (19) and (18), the coefficients of canonical stress-
energy tensor and axial current are related by:

$\begin{split}\mathbb{A}^{\textup{Can}}=&-\left(\frac{n_{\textup{A}}}{|\beta|}+\frac{\partial
n_{\textup{A}}}{\partial|\beta|}\right),\\\
\frac{n_{\textup{A}}}{|\beta|}=&\mathbb{W}_{1}^{\textup{Can}}-\mathbb{W}_{2}^{\textup{Can}},\\\
\frac{W^{\textup{A}}}{|\beta|}=&\frac{G_{1}^{\textup{Can}}-G_{2}^{\textup{Can}}}{2},\end{split}$

which expose an interesting connection between the Axial Vortical Effect (AVE)
conductivity $W^{\textup{A}}$ and the second order thermal coefficients of the
canonical stress-energy tensor.

To understand the constraint (17) and the relation between axial vortical
effect and anomalies, we also consider the case of a free massive field. In
that case, axial current is not conserved, but its divergence is given by

$\partial_{\mu}\widehat{j}_{\textup{A}}^{\mu}=2m{\rm
i}\bar{\Psi}\gamma^{5}\Psi.$

It follows that global equilibrium with a conserved axial charge can not be
reached. Then, to still use the previous global equilibrium analysis to
massive field, we simply set the axial chemical potential to zero and we
consider global equilibrium with just thermal vorticity and finite electric
charge. As a consequence, the symmetries impose that all chiral coefficients
(i.e. those which are not parity invariant) must be vanishing. However, the
term in $W^{\textup{A}}$ of axial current decomposition is not chiral and
consequently could be different from zero. Since the conservation equation is
changed, we expect that also the condition (17) will be modified. We then have
to consider the pseudo-scalar operator ${\rm i}\bar{\Psi}\gamma^{5}\Psi$ that
appears on the divergence of axial current. Pseudoscalar thermal expectation
value can be decomposed at second order in thermal vorticity in the same way
as other local operators and we find that it is given by a single term:

$\langle{\rm i}\bar{\Psi}\gamma^{5}\Psi\rangle=(\alpha\cdot w)L^{\alpha\cdot
w}$

where the non chiral thermal coefficient can be obtained by

$L^{\alpha\cdot
w}=\frac{1}{2}\langle\langle\\{\widehat{K}_{3},\widehat{J}_{3}\\}{\rm
i}\bar{\Psi}\gamma^{5}\Psi\rangle\rangle.$ (20)

With this definition we find that the condition on axial vortical effect
conductivity $W^{\textup{A}}$ becomes:

$|\beta|\frac{\partial
W^{\textup{A}}}{\partial|\beta|}+3W^{\textup{A}}=-2mL^{\alpha\cdot w}.$ (21)

Differently from (17) the constraint (21) no longer imposes $W^{\textup{A}}$
to be proportional to the third power of temperature and $W^{\textup{A}}$
acquires terms which depends on the mass of the fields.

As concluding remarks we give some results for these coefficients for the free
massless Dirac field. In that case this method reproduces the well-know
Kharzeev:2015znc  chiral vortical effect and axial vortical effect
conductivities

$W^{\textup{V}}=\frac{\mu\,\mu\,_{\textup{A}}T}{\pi^{2}},\quad
W^{\textup{A}}=\frac{T^{3}}{6}+\frac{(\mu^{2}+\mu_{\textup{A}}^{2})T}{2\pi^{2}}.$
(22)

In the case of massive Dirac fields, global thermal equilibrium with thermal
vorticity and vanishing axial chemical potential is well defined and the axial
currents mean value can be directed along the rotation of the fluid. In that
situation the AVE conductivity for a free massive Dirac field is
Buzzegoli:2017cqy

$W^{A}=\frac{1}{2\pi^{2}|\beta|}\int_{0}^{\infty}\,\mathrm{d}p\left[n_{F}(E_{p}-\mu)+n_{F}(E_{p}+\mu)\right]\frac{2p^{2}+m^{2}}{E_{p}},$
(23)

where $E_{p}=\sqrt{p^{2}+m^{2}}$. This coefficient is related to pseudo-scalar
thermal coefficient $L^{\alpha\cdot w}$ via Eq. (21) and indeed pseudo-scalar
coefficient is given by

$L^{\alpha\cdot
w}=-\frac{m}{4\pi^{2}\beta^{2}}\int_{0}^{\infty}\frac{\mathrm{d}p}{E_{p}}\left[n^{\prime}_{F}(E_{p}-\mu)+n^{\prime}_{F}(E_{p}+\mu)\right],$
(24)

where the prime on distribution functions stands for derivative respect to
$E_{p}$. We can give approximate results for integral in Eq. (23). For high
temperature regime ($T\gg m$), if the gas is non-degenerate ($|\mu|<m$), we
extract the AVE conductivity behavior using the Mellin transformation
technique LandWeer . The result is

$\frac{W^{\textup{A}}}{T}\simeq\frac{T^{2}}{6}+\frac{\mu^{2}}{2\pi^{2}}-\frac{m^{2}}{4\pi^{2}}-\frac{7\zeta^{\prime}(-2)T^{2}}{8\pi^{2}}\left(\frac{m}{T}\right)^{4}+\mathcal{O}\left(\frac{m^{6}}{T^{6}}\right).$
(25)

The first term in mass was also obtained in Flachi:2017vlp  where the axial
vortical effect was evaluated with the statistical operator (6) but in curved
space-time. Low temperature behavior can also be extracted from (23), see
Buzzegoli:2017cqy . For a degenerate gas ($|\mu|>m$) at zero temperature we
obtain111Notice that $W^{\textup{A}}w^{\mu}\to(W^{\textup{A}}/T)\vec{\omega}$,
so there is no divergency for $T\to 0$.:

$\frac{W^{\textup{A}}}{T}=\frac{\mu^{2}}{2\pi^{2}}\frac{\sqrt{\mu^{2}-m^{2}}}{\mu}.$

Instead for a non degenerate gas ($|\mu|<m$) at low temperature $T<<m$ we have

$W^{\textup{A}}\approx\left(1+2\frac{T}{m}\right)\frac{(mT)^{3/2}}{\sqrt{2}\,\pi^{3/2}}\mathrm{e}^{|\beta|(\mu-m)}.$
(26)

Axial current corrections for rotating and accelerating fluids is also
discussed in Prokhorov:2017atp ; Prokhorov:2018qhq  for both massive and
massless fields using an ansatz for Wigner function with thermal vorticity.

##  3 Dirac Field in external electromagnetic field

Consider a Dirac field in external electromagnetic field. The Lagrangian of
the theory is given by

$\mathcal{L}=\frac{{\rm
i}}{2}\left[\bar{\Psi}\gamma^{\mu}\overrightarrow{\partial}_{\mu}\Psi-\bar{\Psi}\gamma^{\mu}\overleftarrow{\partial}_{\mu}\Psi\right]-m\bar{\Psi}\Psi-\widehat{j}^{\,\mu}A_{\mu},$

where $\widehat{j}^{\,\mu}=q\bar{\Psi}\gamma^{\mu}\Psi$, $q$ is the elementary
electric charge of the field and the gauge potential $A^{\mu}$ is an external
non dynamic field. This Lagrangian is obtained from the free Dirac one with
the minimal coupling substitution $\partial_{\mu}\to\partial_{\mu}+{\rm
i}qA_{\mu}$ which ensures gauge invariance to the theory. From Euler Eq.s we
obtain the equations of motion (EOM) for the Dirac field:

$\not{\partial}\Psi=-{\rm
i}(q\not{A}+m)\Psi,\quad\not{\partial}\bar{\Psi}=\bar{\Psi}{\rm
i}(q\not{A}+m).$

By applying Noether’s theorem to this Lagrangian we obtain the operators in
Eq. (4).

###  3.1 Symmetries in constant electromagnetic field

It is worth noticing that the symmetries of the theory of fermions in external
electromagnetic field are different from those of free fermions and from those
of quantum electrodynamics. While a system without external forces is
symmetric for the full Poincaré group, some of the symmetries are lost when
external fields are introduced. Indeed, external fields do not transform
together with the rest of the system. In this section we discuss the
symmetries of a system in the presence of an external constant homogeneous
electromagnetic field. We will examine the transformations that are still
symmetries of the theory, the consequent conserved quantities and the form of
the generators of such transformations.

If the Lagrangian of our theory is invariant under translations, from
Noether’s theorem, we can identify four operators. Those operators share three
properties: they are conserved quantities, they are the generators of
translations and they constitute the four-momentum of the system. However,
translation invariance by itself does not guarantees that the same quantity
must have all the three above properties altogether. Consider again a system
under an external electromagnetic field. In this situation, Poincaré symmetry
of space-time is broken. Only in the special case of a constant and
homogeneous electromagnetic field translation invariance is restored. However,
the Lagrangian is not invariant under space-time translation, but it acquires
a term that is a four-divergence. This term, under appropriate boundary
conditions, does not affect the action of the system and the overall
invariance is preserved. Nevertheless, the consequence of the additional term
is that we can distinguish between three different operators, each of them
having one of the three properties stated above. This is understood with the
Noether-Tassie-Buchdahl theorem MR0180181 ; MR0180182 ; Bacry:1970ye : given a
Lagrangian $\mathcal{L}(\Psi(x),\partial_{\mu}\Psi(x),x)$ and the
infinitesimal transformation

$x^{\prime\mu}=x^{\mu}+\delta x^{\mu},\quad\Psi^{\prime}=\Psi+\delta\Psi$

such that $\partial_{\mu}\delta x^{\mu}=0$, which transforms the Lagrangian in

$\mathcal{L}(\Psi^{\prime}(x^{\prime}),\partial^{\prime}_{\mu}\Psi^{\prime}(x^{\prime}),x^{\prime})=\mathcal{L}(\Psi(x),\partial_{\mu}\Psi(x),x)+\partial_{\mu}X^{\mu},$

where $X^{\mu}$ is a functional depending exclusively on $\Psi(x)$ and $x$,
the quantity

$\Gamma^{\mu}=\frac{\delta\mathcal{L}}{\delta\partial_{\mu}\Psi}\delta\Psi-\left(\frac{\delta\mathcal{L}}{\delta\partial_{\mu}\Psi}\partial_{\nu}\Psi-\mathcal{L}\,g^{\mu}_{\,\nu}\right)\delta
x^{\nu}-X^{\mu}$

is conserved, i.e. divergence-less.

Consider the Dirac Lagrangian in constant homogeneous electromagnetic field

$\mathcal{L}(\Psi(x),\partial_{\mu}\Psi(x),x)=\bar{\Psi}({\rm
i}\not{\partial}-m)\Psi-\widehat{j}^{\,\mu}A_{\mu}.$

The translation transformation ($\delta\Psi=0$, $\delta
x^{\mu}=\epsilon^{\mu}$) acts on the Dirac field but does not act directly on
the external gauge field. Therefore, a translation changes the Lagrangian by

$\begin{split}\delta\mathcal{L}=&\mathcal{L}(\Psi^{\prime}(x^{\prime}),\partial^{\prime}_{\mu}\Psi^{\prime}(x^{\prime}),x^{\prime})-\mathcal{L}(\Psi(x),\partial_{\mu}\Psi(x),x)\\\
=&\widehat{j}^{\,\mu}\partial_{\nu}A_{\mu}\epsilon^{\nu}=\widehat{j}^{\,\mu}(F_{\nu\mu}+\partial_{\mu}A_{\nu})\epsilon^{\nu}=-\widehat{j}^{\,\mu}(F_{\mu\nu}-\partial_{\mu}A_{\nu})\epsilon^{\nu}.\end{split}$

The quantity $X$ of Noether-Tassie-Buchdahl in this case is

$X^{\mu\nu}\equiv\widehat{j}^{\,\mu}(A^{\nu}+F^{\nu\lambda}x_{\lambda}).$

Indeed its divergence is the variation of the Lagrangian

$\epsilon_{\nu}\partial_{\mu}X^{\mu\nu}=(\partial_{\mu}\widehat{j}^{\,\mu})(A^{\nu}+F^{\nu\lambda}x_{\lambda})\epsilon_{\nu}+\widehat{j}^{\,\mu}(\partial_{\mu}A_{\nu}+F_{\nu\mu})\epsilon^{\nu}=-\widehat{j}^{\,\mu}(F_{\mu\nu}-\partial_{\mu}A_{\nu})\epsilon^{\nu}=\delta\mathcal{L}.$

Therefore the theorem implies that the system has a canonical conserved tensor
given by:

$\widehat{\pi}^{\,\mu\nu}_{\textup{can}}=\widehat{T}^{\mu\nu}_{0}-\widehat{j}^{\,\mu}A^{\nu}-\widehat{j}^{\,\mu}F^{\nu\lambda}x_{\lambda},$

where $\widehat{T}^{\mu\nu}_{0}$ is the free canonical Dirac stress-energy
tensor. Using Belinfante procedure we can transform
$\widehat{T}^{\mu\nu}_{0}-\widehat{j}^{\,\mu}A^{\nu}$ into the symmetric
stress-energy tensor of Dirac field in external magnetic field
$\widehat{T}^{\mu\nu}_{S}$ and the above conserved tensor can be written as:

$\widehat{\pi}^{\,\mu\nu}\equiv\widehat{T}^{\mu\nu}_{S}-\widehat{j}^{\,\mu}F^{\nu\lambda}x_{\lambda}.$
(27)

From the above equation, we can simply verify that
$\partial_{\mu}\widehat{\pi}^{\,\mu\nu}=0$ form Eq. (2). Note that
$\widehat{\pi}^{\,\mu\nu}$ is not symmetric and that it is gauge invariant.
The conserved quantities are obtained from the previous operators by

$\widehat{\pi}^{\,\mu}=\int\mathrm{d}^{3}x\,\widehat{\pi}^{\,0\mu}$

and we can show that this four-vector constitutes the generators of the
translation Bacry:1970ye . However, the momentum of the system is still given
by

$\widehat{P}^{\mu}=\int\mathrm{d}^{3}x\,\widehat{T}^{0\mu}$

but it is no longer a conserved quantity and it is no longer the generator of
translations. Another difference with the four-momenta is that different
components of this vectors do not commute, instead they satisfy the
commutation relation Bacry:1970ye

$[\widehat{\pi}^{\,\mu},\widehat{\pi}^{\,\nu}]={\rm i}\widehat{Q}F^{\mu\nu},$

where $\widehat{Q}$ is the electric charge operator.

As for Lorentz’s transformations we expect that the variation of the
Lagrangian is a full divergence only for specific forms of transformations.
For example, with a vanishing electric field and a constant magnetic field,
only the rotation along the direction of the magnetic field and the boost
along the magnetic field are symmetries of the theory. Therefore, only in
these cases, the Lagrangian variation could be vanishing or a full divergence.

Returning to the general case, by repeating the previous argument made for the
translations for the Lorentz transformation:

$\delta x^{\mu}=\omega^{\mu\nu}x_{\nu},\quad\delta\Psi=-\frac{{\rm
i}}{2}\omega_{\mu\nu}\sigma^{\mu\nu}\Psi,$

we find that the transformed Lagrangian is

$\begin{split}\delta\mathcal{L}=&\omega^{\mu\nu}\widehat{j}^{\lambda}\left[\partial_{\lambda}\left(x_{\mu}A_{\nu}-x_{\nu}A_{\mu}\right)+x_{\mu}F_{\nu\lambda}-x_{\nu}F_{\mu\lambda}\right].\end{split}$

We can show that the Lagrangian variation can also be written as

$\begin{split}\delta\mathcal{L}=&\frac{1}{2}\omega^{\mu\nu}\partial_{\lambda}\left[\widehat{j}^{\lambda}x_{\mu}\left(A_{\nu}-\frac{1}{2}x^{\sigma}F_{\sigma\nu}\right)-\widehat{j}^{\lambda}x_{\nu}\left(A_{\mu}-\frac{1}{2}x^{\sigma}F_{\sigma\mu}\right)\right]\\\
&-\frac{1}{2}x^{\rho}\widehat{j}^{\lambda}(\omega_{\lambda\sigma}F^{\sigma}_{\,\,\rho}-\omega_{\rho\sigma}F^{\sigma}_{\,\,\lambda}).\end{split}$

The first term of the r.h.s. is written as a four-divergence. The remaining
part can not be cast into a four-divergence but it is proportional to the
following product

$(\omega\wedge
F)_{\lambda\rho}=\omega_{\lambda\sigma}F^{\sigma}_{\,\,\rho}-\omega_{\rho\sigma}F^{\sigma}_{\,\,\lambda}.$

The product of two non-vanishing anti-symmetric tensor of rank two,
$\omega\wedge F$, is zero if and only if $\omega$ is a linear combination of
$F$ and its dual $F^{*}$ (or viceversa) Bacry:1970ye :

$(\omega\wedge
F)_{\lambda\rho}=0\quad\text{iff}\quad\omega_{\mu\nu}=k\,F_{\mu\nu}+k^{\prime}\,F^{*}_{\mu\nu},\quad
k,k^{\prime}\in\mathbb{R}.$ (28)

Therefore, the part of Lagrangian variation which is not a divergence is
vanishing when $\omega^{\mu\nu}$ is a linear combination of electromagnetic
stress-energy tensor and its dual:

$\omega^{\mu\nu}=a\,F^{\mu\nu}+\frac{b}{2}\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}.$
(29)

This means, as expected, that the theory is invariant only under certain type
of Lorentz transformations: the ones generated with parameters of the form
(29). For example, in the case of constant magnetic field, we recover that the
system is invariant only for rotation and boost along the magnetic field. Set
then $\omega$ either as $\omega_{\mu\nu}\propto F_{\mu\nu}$ or
$\omega_{\mu\nu}\propto F^{*}_{\mu\nu}$, so that the Lagrangian variation is a
four-divergence. In this case we can apply Noether-Tassie-Buchdahl theorem and
the two following quantities are divergence-less:

$\begin{split}\widehat{\Gamma}^{\lambda}=&\frac{F^{\mu\nu}}{2}\left[x_{\mu}\left(\widehat{T}^{\lambda}_{0\,\,\nu}-\widehat{j}^{\,\lambda}A_{\nu}+\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\nu}\right)\right.\\\
&\left.-x_{\nu}\left(\widehat{T}^{\lambda}_{0\,\,\mu}-\widehat{j}^{\,\lambda}A_{\mu}+\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\mu}\right)+\widehat{S}^{\lambda}_{\,\mu\nu}\right],\\\
\widehat{\Gamma}^{*\lambda}=&\frac{F^{*\mu\nu}}{2}\left[x_{\mu}\left(\widehat{T}^{\lambda}_{0\,\,\nu}-\widehat{j}^{\,\lambda}A_{\nu}+\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\nu}\right)\right.\\\
&\left.-x_{\nu}\left(\widehat{T}^{\lambda}_{0\,\,\mu}-\widehat{j}^{\,\lambda}A_{\mu}+\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\mu}\right)+\widehat{S}^{\lambda}_{\,\mu\nu}\right],\end{split}$

where $\widehat{S}^{\lambda}_{\mu\nu}$ is the canonical spin tensor of free
Dirac field. After Belinfante transformation, the quantities become

$\begin{split}\widehat{\Gamma}^{\lambda}=&\frac{1}{2}F^{\mu\nu}\widehat{M}^{\lambda}_{\mu\nu},\quad\widehat{\Gamma}^{*\lambda}=\frac{1}{2}F^{*\mu\nu}\widehat{M}^{\lambda}_{\mu\nu},\\\
\widehat{M}^{\lambda}_{\mu\nu}\equiv&x_{\mu}\left(\widehat{T}^{\lambda}_{S\,\nu}+\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\nu}\right)-x_{\nu}\left(\widehat{T}^{\lambda}_{S\,\mu}+\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\mu}\right)\\\
=&x_{\mu}\left(\widehat{\pi}^{\,\lambda}_{\,\,\,\nu}-\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\nu}\right)-x_{\nu}\left(\widehat{\pi}^{\,\lambda}_{\,\,\,\mu}-\frac{1}{2}\widehat{j}^{\,\lambda}x^{\rho}F_{\rho\mu}\right).\end{split}$
(30)

We can define the integrals:

$\widehat{M}_{\mu\nu}=\int\mathrm{d}^{3}x\,\widehat{M}^{0}_{\mu\nu},$

which are conserved quantities only if contracted with $F^{\mu\nu}$ or
$F^{*\mu\nu}$. The operators $\widehat{M}_{\mu\nu}$ are the generators of
Lorentz transformations if they are also a symmetry for the theory, otherwise
the Wigner’s theorem does not applies and we can not say that such
transformations admit an unitary and linear (or anti-unitary and anti-linear)
representation. For those operators, the following Algebra holds true
Bacry:1970ye :

$\begin{split}[\widehat{\pi}^{\,\mu},\widehat{\pi}^{\,\nu}]=&{\rm
i}F^{\mu\nu}\widehat{Q},\\\
\frac{1}{2}F_{\rho\sigma}[\widehat{\pi}{\,{}^{\mu}},\widehat{M}^{\rho\sigma}]=&\frac{{\rm
i}}{2}F_{\rho\sigma}\left(\eta^{\mu\rho}\widehat{\pi}^{\,\sigma}-\eta^{\mu\sigma}\widehat{\pi}^{\,\rho}\right),\\\
\frac{1}{2}F^{*}_{\rho\sigma}[\widehat{\pi}^{\,\mu},\widehat{M}^{\rho\sigma}]=&\frac{{\rm
i}}{2}F^{*}_{\rho\sigma}\left(\eta^{\mu\rho}\widehat{\pi}^{\,\sigma}-\eta^{\mu\sigma}\widehat{\pi}^{\,\rho}\right),\end{split}$
(31)

where $\widehat{Q}$ is the electric charge operator. In the particular case of
vanishing electric field and constant magnetic field along the $z$ axis, the
Algebra becomes:

$\displaystyle[\widehat{\pi}_{x},\,\widehat{\pi}_{y}]={\rm
i}|\vec{B}|\widehat{Q},$
$\displaystyle[\widehat{J}_{z},\widehat{\pi}_{x}]={\rm
i}\widehat{\pi}_{y},\quad[\widehat{J}_{z},\widehat{\pi}_{y}]=-{\rm
i}\widehat{\pi}_{x},$ $\displaystyle[\widehat{K}_{z},\widehat{\pi}_{t}]=-{\rm
i}\widehat{\pi}_{z},\quad[\widehat{K}_{z},\widehat{\pi}_{z}]=-{\rm
i}\widehat{\pi}_{t}.$

##  4 Chiral fermions in constant magnetic field

Consider now a system consisting of free chiral fermions in an external
homogeneous constant magnetic field $\vec{B}$ at global thermal equilibrium
with vanishing thermal vorticity. In this configuration, the chiral anomaly is
vanishing because there is no electric field ($B\cdot E=0$). It follows that
global equilibrium without vorticity is described by a constant inverse four-
temperature $\beta$ and a constant axial chemical potential
$\zeta_{\textup{A}}$ (see Sec. 2). Instead, the condition for the electric
chemical potential reads

$\partial^{\mu}\zeta(x)=F^{\nu\mu}\beta_{\nu}=\sqrt{\beta^{2}}F^{\nu\mu}u_{\nu}=\sqrt{\beta^{2}}E^{\mu},$

where $u$ is the fluid velocity directed along $\beta$ and $E$ is the comoving
electric field. Since we are considering the case without electric field the
global equilibrium condition is simply a constant $\zeta$. The global
equilibrium statistical operator then becomes

$\widehat{\rho}=\frac{1}{Z}\exp\left[-\widehat{P}^{\mu}\beta_{\mu}+\zeta\widehat{Q}+\zeta_{\textup{A}}\widehat{Q}_{\textup{A}}\right].$

Notice that the operators $\widehat{P}^{\mu}$ are not the generators of
translations, which are instead given by $\widehat{\pi}^{\,\mu}$ and are
obtained by integrating the conserved current in Eq. (27). However, in the
case of vanishing comoving electric field, the projection of the inverse
temperature along the four-momentum is equivalent to the projection along of
the generators of translations, that is:

$\widehat{\pi}^{\,\mu}\beta_{\mu}=\int_{\Sigma}\mathrm{d}\Sigma_{\lambda}\left(\widehat{T}^{\lambda\nu}-\widehat{j}^{\,\lambda}F^{\nu\sigma}x_{\sigma}\right)\beta_{\nu}=\widehat{P}^{\mu}\beta_{\mu}-\sqrt{\beta^{2}}E^{\sigma}\int_{\Sigma}\mathrm{d}\Sigma_{\lambda}\widehat{j}^{\,\lambda}x_{\sigma}=\widehat{P}^{\mu}\beta_{\mu}.$

The statistical operator can now be written as

$\widehat{\rho}=\frac{1}{Z}\exp\left[-\widehat{\pi}^{\,\mu}\beta_{\mu}+\zeta\widehat{Q}+\zeta_{\textup{A}}\widehat{Q}_{\textup{A}}\right].$

In this form, it is straightforward to use the algebra in Eq. (31) and
translate the statistical operator of a quantity $a^{\mu}$. We find

$\widehat{\sf T}(a)\,\widehat{\rho}\,\widehat{\sf T}^{-1}(a)=\mathrm{e}^{{\rm
i}a\cdot\widehat{\pi}}\widehat{\rho}\,\mathrm{e}^{-{\rm
i}a\cdot\widehat{\pi}}=\frac{1}{Z}\exp\left[-\widehat{\pi}^{\,\mu}\beta_{\mu}+\zeta\widehat{Q}+\zeta_{\textup{A}}\widehat{Q}_{\textup{A}}+a_{\mu}F^{\mu\nu}\beta_{\nu}\widehat{Q}\right].$

Since $F^{\mu\nu}\beta_{\nu}$ is the comoving electric field, which is
vanishing, the statistical operator is homogeneous:

$\widehat{\sf T}(a)\,\widehat{\rho}\,\widehat{\sf T}^{-1}(a)=\widehat{\rho}.$

###  4.1 Exact thermal solutions

Having established the basic quantities of thermal equilibrium with a constant
and homogeneous magnetic field, we now move on to introduce the techniques of
thermal field theory in order to find exact solutions for thermodynamic
equilibrium in a magnetic field. We start by giving a path integral
description of the partition function. Since the partition function is a
Lorentz invariant, we can choose the evaluate it in the local rest frame where
$u=(1,\vec{0})$. In this frame, without loss in generality, the magnetic field
is chosen along the $z$ axis and we adopt the Landau gauge
$A^{\mu}=(0,0,Bx_{1},0)$. The path integral formulation of the partition
function in local rest frame

$Z(T,\mu,\mu_{5})={\text{tr}}\left[\mathrm{e}^{-\beta(\widehat{H}-\mu\widehat{Q}-\mu_{\textup{A}}\widehat{Q}_{\textup{A}})}\right]$

is given by222We added a mass term for generalization, although with mass we
cannot have a conserved axial current.

$Z=C\int_{\Psi(\beta,\vec{x})=-\Psi(0,\vec{x})}\mathcal{D}\bar{\Psi}\,\mathcal{D}\Psi\,\,\exp\left(-S_{\textup{E}}(\Psi,\bar{\Psi},\mu_{\textup{A}})\right)$

where the Euclidean action of Dirac fermions in external electromagnetic field
is

$S_{\textup{E}}(\Psi,\bar{\Psi},\mu_{5})=\int_{0}^{\beta}\mathrm{d}\tau\int_{\vec{x}}\bar{\Psi}(X)\left[{\rm
i}(\gamma\cdot\pi^{+})+m-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right]\Psi(X)$

and $\pi^{+}_{\mu}\equiv P^{+}_{\mu}-qA_{\mu}$, which is not to be confused
with the generators of translations.

With regard to the exact solution, instead of solving the Dirac equation
directly, we use the Ritus method Ritus:1972ky  (see Leung:2005yq  for a brief
recap of the method). The core concept of Ritus method is that we can
construct a complete set of orthonormal function, called $E_{p}$ Ritus
functions, such that the Euclidean action is rendered formally identical to
the Euclidean action of a free Dirac field in absence of external fields. The
$E_{p}$ functions are constructed such that they are the matrix of the
contemporaneous eigenfunctions (eigenvectors) of the maximal set of mutually
commuting operators $\\{(\gamma\cdot\pi)^{2},{\rm
i}\gamma_{1}\gamma_{2},\gamma^{5}\\}$. From gamma algebra, it is
straightforward to check that

${\rm
i}\gamma_{1}\gamma_{2}\delta(\sigma)=\sigma\Delta(\sigma),\quad\frac{1+\chi\gamma^{5}}{2}\gamma^{5}=\chi\frac{1+\chi\gamma^{5}}{2},$

with $\sigma=\pm$ and $\chi=\pm$ and we defined

$\Delta(\sigma)\equiv\frac{1+{\rm i}\sigma\gamma_{1}\gamma_{2}}{2}.$

We can then show that Ritus:1972ky ; Leung:2005yq

$(\gamma\cdot\pi^{+})^{2}E_{\widehat{p}\sigma}(X)={\underline{P}^{+}}^{2}E_{\widehat{p}\sigma}(X)$

where $\widehat{p}$ is a label for the quantum numbers
$\\{l,\omega_{n},p_{2},p_{3}\\}$, the eigenvalues $\underline{P}^{+}$ are
given by

$\underline{P}^{+}\equiv(\omega_{n}+{\rm
i}\mu,0,-\bar{\sigma}\sqrt{2|qB|l},p_{3}),\quad\bar{\sigma}\equiv\text{sgn}(qB)$

and the form of eigenfunction is

$E_{\widehat{p}\sigma}(X)=N(n)\mathrm{e}^{{\rm
i}(P_{0}\tau+P_{2}X_{2}+P_{3}X_{3})}D_{n}(\rho)$ (32)

where $N(n)=(4\pi|qB|)^{1/4}/\sqrt{n!}$ is a normalization factor, and
$D_{n}(\rho)$ denotes the parabolic cylinder functions with argument
$\rho=\sqrt{2|qB|}(X_{1}-p_{2}/qB)$ and non-negative integer index
$n=0,1,2,\dots$ given by

$n=l+\frac{\sigma}{2}\text{sgn}(qB)-\frac{1}{2}.$

Note that the form of the functions (32) strongly depends on the gauge chosen,
in our case the Landau gauge. Since the eigenfunction
$E_{\widehat{p}\sigma}(X)$ does not depend on chirality, the maximal
eigenfunctions of the operators $\\{(\gamma\cdot\pi)^{2},{\rm
i}\gamma_{1}\gamma_{2},\gamma^{5}\\}$ are given by

$E_{\widehat{p}}(X)=\sideset{}{{}^{\prime}}{\sum}_{\sigma=\pm}E_{\widehat{p}\sigma}(X)\Delta(\sigma),\quad\bar{E}_{\widehat{p}}(X)=\gamma_{0}E^{\dagger}_{\widehat{p}}(X)\gamma_{0}=\sideset{}{{}^{\prime}}{\sum}_{\sigma^{\prime}=\pm}E^{*}_{\widehat{p}\sigma^{\prime}}(X)\Delta(\sigma^{\prime})$
(33)

where the prime on the summation symbol denotes that the sum is subject to the
constraint

$\sigma=\begin{cases}\text{sgn}(qB)&l=0\\\ \pm&l>0\end{cases}.$

Some important properties can be derived from these definitions. Firstly, that
the functions $E_{p}$ commute with $\gamma_{0}$ and with $\gamma^{5}$.
Secondly, they satisfy the orthogonality relation

$\int_{X}\bar{E}_{\widehat{q}}(X)E_{\widehat{p}}(X)=(2\pi)^{4}\widehat{\delta}^{(4)}(\widehat{p}-\widehat{q})\Pi(l)$

where we defined

$\displaystyle\delta^{(4)}(\widehat{p}-\widehat{p}^{\prime})\equiv$
$\displaystyle\delta_{l,l^{\prime}}\beta\delta_{\omega_{n},\omega_{n^{\prime}}}\delta(p_{2}-p_{2}^{\prime})\delta(p_{3}-p_{3}^{\prime})$
$\displaystyle\Pi(l)\equiv$ $\displaystyle\begin{cases}\frac{1+{\rm
i}\bar{\sigma}\gamma_{1}\gamma_{2}}{2}&l=0\\\ 1&l>0\end{cases}.$

And lastly, the action of the operator $(\gamma\cdot\pi^{+})$ on these
function is

$(\gamma\cdot\pi^{+})E_{\widehat{p}}(X)=E_{\widehat{p}}(X)\gamma\cdot\underline{P}.$

Since we showed that $E_{p}$ Ritus functions are a complete orthonormal
functions, we can expand the Dirac fields in these functions:

$\begin{split}\Psi(X)=&T\sum_{\\{\omega_{n}\\}}\sum_{l=0}^{\infty}\int\mathrm{d}p_{2}\int\frac{\mathrm{d}p_{3}}{(2\pi)^{3}}E_{\widehat{p}}(X)\Psi(\underline{P})\equiv\sumint_{\widehat{P}}E_{\widehat{p}}(X)\Psi(\underline{P}),\\\
\bar{\Psi}(X)=&\sumint_{\widehat{Q}}\bar{\Psi}(\underline{Q})\bar{E}_{\widehat{q}}(X).\end{split}$

Replacing this expansion on the Euclidean action, we find

$S_{\textup{E}}(\Psi,\bar{\Psi},\mu_{5})=\sumint_{\widehat{P}}\sumint_{\widehat{Q}}\int_{X}\bar{\Psi}(\underline{Q})\bar{E}_{\widehat{q}}(X)\left[{\rm
i}(\gamma\cdot\pi^{+})+m-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right]E_{\widehat{p}}(X)\Psi(\underline{P})$

and, using the above mentioned properties of $E_{p}$ functions, we obtain

$S_{\textup{E}}(\Psi,\bar{\Psi},\mu_{5})=\sumint_{\widehat{P}}\bar{\Psi}(\underline{P})\Pi(l)\left[{\rm
i}\gamma\cdot\underline{P}^{+}+m-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right]\Psi(\underline{P}).$

Notice that this is formally identical to the Euclidean action of free Dirac
field.

We can now proceed to evaluate the partition function. We first change the
integration variables in the partition function to the modes of $E_{p}$
functions $\bar{\Psi}(\underline{P})$ and $\Psi(\underline{P})$. The partition
function is then a Gaussian integral of Grassmann variables, whose result is
the exponent determinant. Hence the partition function becomes

$\begin{split}Z=&\tilde{C}\int_{\Psi(\beta,\vec{x})=-\Psi(0,\vec{x})}\mathcal{D}\bar{\Psi}(\underline{P})\,\mathcal{D}\Psi(\underline{P})\times\\\
&\times\exp\left\\{-\sumint_{\widehat{P}}\bar{\Psi}(\underline{P})\Pi(l)\left[{\rm
i}\gamma\cdot\underline{P}^{+}+m-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right]\Psi(\underline{P})\right\\}\\\
=&\tilde{C}\det\left[\Pi(l)\left({\rm
i}\gamma\cdot\underline{P}^{+}+m-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right)\right].\end{split}$
(34)

For the sake of clarity, for now on we will remove the factor $\Pi(l)$:

$\begin{split}Z=&\tilde{C}\det\left(\begin{array}[]{cc}m\mathbb{I}_{2\times
2}&[{\rm i}(\omega_{n}+{\rm i}\mu)-\mu_{\textup{A}}]\mathbb{I}_{2\times
2}+\sigma_{i}\underline{P}_{i}\\\ ({\rm i}(\omega_{n}+{\rm
i}\mu)+\mu_{\textup{A}})\mathbb{I}_{2\times
2}-\sigma_{i}\underline{P}_{i}&m\mathbb{I}_{2\times
2}\end{array}\right).\end{split}$

The determinant is evaluated using the standard formula for block matrices

$\det\left(\begin{array}[]{cc}A&B\\\ C&D\end{array}\right)=\det\left(AD-
BD^{-1}CD\right);$

replacing that into the partition function we have

$\begin{split}Z&=\tilde{C}\det\left[({\underline{P}^{+}}^{2}+m^{2}+\mu_{\textup{A}}^{2})\mathbb{I}_{2\times
2}-2\sigma_{i}\underline{P}_{i}\mu_{\textup{A}}\right]\\\
&=\tilde{C}\det\left[({\underline{P}^{+}}^{2}+m^{2}+\mu_{\textup{A}}^{2})^{2}-4|\vec{\underline{P}}|^{2}\mu_{\textup{A}}^{2}\right]\\\
&=\prod_{\omega_{n},l,p_{3}}\tilde{C}\left[({\underline{P}^{+}}^{2}+m^{2}+\mu_{\textup{A}}^{2})^{2}-4|\vec{\underline{P}}|^{2}\mu_{\textup{A}}^{2}\right],\end{split}$

where we evaluated the determinant as the product of the eigenvalues of the
matrix. To connect this quantity to the thermodynamics of the system, we are
actually interested in its logarithm:

$\log
Z=\sum_{\omega_{n},l,p_{3}}\log\left[({\underline{P}^{+}}^{2}+m^{2}+\mu_{\textup{A}}^{2})^{2}-4|\vec{\underline{P}}|^{2}\mu_{\textup{A}}^{2}\right]+\text{cnst}.$
(35)

In the next subsection we evaluate the thermodynamic potential of the system
from the logarithm of partition function. Then, by simple derivation, we could
obtain other thermodynamic properties. However, starting from partition
function in Eq. (34), we will obtain the thermal propagator of a chiral
fermion in a magnetic field. Once we have the propagator, we can use the
point-splitting procedure to evaluate other thermal properties that are not
related to the thermodynamic potential. We will use the thermal propagator to
evaluate the mean value of the electric current and of the axial current in
the following subsections.

###  4.2 Thermodynamic potential

The thermodynamic potential $\Omega$ is derived from the partition function as
following

$\Omega=\lim_{V\to\infty}-\frac{T}{V}\log Z,$

where the logarithm of partition function is given by Eq. (35). We can follow
the usual techniques used for free fermions to evaluate the thermodynamic
potential and to sum the Matsubara frequencies. However, we must first
consider that in this case the Landau levels generated by the magnetic field
have different degeneracy factors ad must be properly taken into account when
performing the infinite volume limit. Let be $S$ the area in the $x-y$ plane
and $p_{\perp 1}$ and $p_{\perp 2}$ the momenta in that plane. In the limit of
infinity area, the sum on modes becomes the following integrals:

$\lim_{S\to\infty}\frac{1}{S}\sum_{p_{\perp 1}}\sum_{p_{\perp
2}}=\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{\perp
1}}{2\pi}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{\perp 2}}{2\pi}.$

Each Landau level has a degeneracy associated with some quantum numbers; this
degeneracy is gauge independent and it is given by:

$d_{l}=\left\lfloor\frac{|qB|S}{2\pi}\right\rfloor.$

To obtain this degeneracy, we just have to evaluate the quantity

$\frac{\mathrm{d}p_{\perp 1}}{2\pi}\frac{\mathrm{d}p_{\perp 2}}{2\pi}$

between two consecutive energy levels:

$d_{l}=\left\lfloor\frac{|qB|S}{2\pi}\right\rfloor=\int_{l}^{l+1}\frac{\mathrm{d}p_{\perp
1}}{2\pi}\frac{\mathrm{d}p_{\perp 2}}{2\pi}.$

Therefore, removing the floor function, the infinite volume limit of the sum
on the states of the system gives:

$\lim_{V\to\infty}\frac{1}{V}\sum_{\omega_{n},l,p_{3}}=\frac{|qB|}{2\pi}\sum_{l=0}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{2\pi}\sum_{\\{\omega_{n}\\}}.$

Consequently the thermodynamic potential reads:

$\begin{split}\Omega=&\lim_{V\to\infty}-\frac{T}{V}\log Z\\\
=&-\frac{|qB|}{2\pi}\sum_{l=0}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{2\pi}T\sum_{\\{\omega_{n}\\}}\log\left[({\underline{P}^{+}}^{2}+m^{2}+\mu_{\textup{A}}^{2})^{2}-4|\vec{\underline{P}}|^{2}\mu_{\textup{A}}^{2}\right]+\text{cnst}.\end{split}$

The Matsubara sum can be performed as describe in Laine:2016hma . The final
result for thermodynamic potential is

$\Omega=-\frac{|qB|}{2\pi}\sum_{l=0}^{\infty}\sideset{}{{}^{\prime}}{\sum}_{s=\pm}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{2\pi}\left[E_{s}+T\sum_{\pm}\log\left(1+\mathrm{e}^{-\beta(E_{s}\pm\mu)}\right)\right]+\text{cnst},$

where $E_{s}^{2}=[(p_{3}^{2}+2qBl)^{1/2}+s\mu_{\textup{A}}]^{2}+m^{2}$ and the
constraint of $s=\bar{\sigma}$ for $l=0$ is caused by the projector $\Pi(l)$.
This same thermodynamic potential for chiral fermions in external magnetic
field was used in Fukushima:2008xe  to derive the chiral magnetic effect
(CME). This expression can be used to obtain the electric and axial charge
density, but instead we are using the point-splitting procedure because the
latter can also be used to evaluate other thermodynamic functions related to
currents. To do that, we first need the thermal propagator of chiral fermions.

###  4.3 Chiral fermion propagator in magnetic field

Since we have used the Ritus method, the form of Euclidean action is formally
identical to those of a free Dirac field. It is therefore not surprising that
the fermionic propagator is obtained in the same way as the free case. The
propagator in Fourier modes can be obtained in path integral formulation by
Laine:2016hma

$\langle\tilde{\Psi}_{a}(P)\bar{\tilde{\Psi}}_{b}(Q)\rangle_{T}=\frac{\int\mathcal{D}\tilde{\Psi}\,\mathcal{D}\tilde{\bar{\Psi}}\,\,\exp\left(-S_{\textup{E}}\right)\tilde{\Psi}_{a}(P)\bar{\tilde{\Psi}}_{b}(Q)}{\int\mathcal{D}\tilde{\Psi}\,\mathcal{D}\tilde{\bar{\Psi}}\,\,\exp\left(-S_{\textup{E}}\right)},$

where in our case the form of the partition function is given in the first
line of Eq. (34):

$\begin{split}Z=&\tilde{C}\int_{\Psi(\beta,\vec{x})=-\Psi(0,\vec{x})}\mathcal{D}\bar{\Psi}(\underline{P})\,\mathcal{D}\Psi(\underline{P})\,\,\times\\\
&\times\exp\left\\{-\sumint_{\widehat{P}}\bar{\Psi}(\underline{P})\Pi(l)\left[{\rm
i}\gamma\cdot\underline{P}^{+}+m-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right]\Psi(\underline{P})\right\\}.\end{split}$

The Grassmann integrals are straightforward and gives

$\langle\bar{\Psi}(\underline{Q})_{a}\Psi(\underline{P})_{b}\rangle=\delta^{(4)}(\underline{P}-\underline{Q})\mathcal{M}^{-1}_{ab}$

where $a,b$ denotes spinorial indices and

$\mathcal{M}=\left[{\rm
i}\gamma\cdot\underline{P}^{+}-\gamma_{0}\gamma^{5}\mu_{\textup{A}}\right].$

The inverse of $\mathcal{M}$ is easily written in terms of the projector into
right and left chirality states, which are defined by:

$\mathbb{P}_{\chi}=\frac{1+\chi\gamma_{5}}{2},\quad\text{i.e.}\quad\mathbb{P}_{\textup{R}}=\frac{1+\gamma_{5}}{2},\quad\mathbb{P}_{\textup{L}}=\frac{1-\gamma_{5}}{2}.$

We also introduce right and left chemical potential:

$\mu_{\textup{R}}\equiv\mu+\mu_{\textup{A}},\quad\mu_{\textup{L}}\equiv\mu-\mu_{\textup{A}},$

and we define right and left charged momenta by

$P_{\textup{R/L}}^{\pm}\equiv(\omega_{n}\pm{\rm
i}\mu_{\textup{R/L}},\vec{p}).$

With this notation, after inverting $\mathcal{M}$, the thermal propagator is

$\langle\bar{\Psi}(\underline{Q})_{a}\Psi(\underline{P})_{b}\rangle=\delta^{(4)}(\underline{P}-\underline{Q})\sum_{\chi}\left(\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right)_{ab}.$

This is the generalization in Euclidean space-time with chemical potentials of
the propagator in Ritus:1972ky ; Leung:2005yq . In the configuration space,
the two-point function is

$\begin{split}\langle\bar{\Psi}(X)_{a}\Psi(Y)_{b}\rangle&=\sumint_{\widehat{P}}\sumint_{\widehat{Q}}\bar{E}_{\widehat{p}}(X)_{a^{\prime}a}E_{\widehat{q}}(Y)_{bb^{\prime}}\langle\bar{\Psi}(\underline{P})_{a^{\prime}}\Psi(\underline{Q})_{b^{\prime}}\rangle\\\
&=-\sumint_{\widehat{P}}\sumint_{\widehat{Q}}\bar{E}_{\widehat{p}}(X)_{a^{\prime}a}E_{\widehat{q}}(Y)_{bb^{\prime}}\langle\Psi(\underline{Q})_{b^{\prime}}\bar{\Psi}(\underline{P})_{a^{\prime}}\rangle\\\
&=-\sumint_{\widehat{P}}\sumint_{\widehat{Q}}\bar{E}_{\widehat{p}}(X)_{a^{\prime}a}E_{\widehat{q}}(Y)_{bb^{\prime}}\delta^{(4)}(\underline{P}-\underline{Q})\sum_{\chi}\left(\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right)_{b^{\prime}a^{\prime}}\end{split}$

where to go to second line we used fermion anti-commutation. Finally,
integrating the delta we have

$\langle\bar{\Psi}(X)_{a}\Psi(Y)_{b}\rangle=-\sumint_{\widehat{P}}\sum_{\chi}E_{\widehat{p}}(Y)_{bb^{\prime}}\left(\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right)_{b^{\prime}a^{\prime}}\bar{E}_{\widehat{p}}(X)_{a^{\prime}a}.$
(36)

###  4.4 Electric current mean value

Having derived the propagator, we now proceed to evaluate the mean value of
electric current. The following method is similar to the one used in
Fukushima:2009ft , as we both use the Ritus method. We take advantage of the
point-splitting procedure to compute the thermal expectation value of electric
current. First, we write the current in Euclidean space-time and we split the
coordinate point in which the fields are evaluated as follows:

$\langle\widehat{j}_{\mu}(X)\rangle\mskip-3.0mu=\mskip-3.0mu(-{\rm
i})^{1-\delta_{\mu,0}}q\langle\widehat{\bar{\Psi}}(X)\gamma_{\mu}\widehat{\Psi}(X)\rangle\mskip-3.0mu=\mskip-3.0mu\lim_{X_{1},X_{2}\to
X}(-{\rm
i})^{1-\delta_{\mu,0}}q\left(\gamma_{\mu}\right)_{ab}\langle\widehat{\bar{\Psi}}_{a}(X_{1})\widehat{\Psi}_{b}(X_{2})\rangle.$

Then we plug the form of fermionic propagator (36) and we reconstruct the
trace on spinorial indices; eventually we obtain

$\langle\widehat{j}_{\mu}(X)\rangle=-(-{\rm
i})^{1-\delta_{\mu,0}}q\sumint_{\widehat{P}}\sum_{\chi}{\text{tr}}\left[\bar{E}_{\widehat{p}}(X)\,\gamma_{\mu}\,E_{\widehat{p}}(X)\,\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right].$

It is convenient to indicate the components parallel to the magnetic field,
which are the time component and the $z$ component, with the parallel symbol
“$\parallel$”. For those components the following commutator holds true:

$\left[\gamma^{\parallel}_{\mu}\,,E_{\widehat{p}}(X)\right]=0.$

Therefore, reminding the definitions of the Ritus $E_{\widehat{p}}(X)$
functions (33), for the parallel component of electric current we obtain

$\begin{split}\langle\widehat{j}^{\,\parallel}_{\mu}(X)\rangle=&-(-{\rm
i})^{1-\delta_{\mu,0}}q\sum_{l=0}^{\infty}\int\frac{\mathrm{d}p_{2}}{2\pi}\int\frac{\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\sideset{}{{}^{\prime}}{\sum}_{\sigma,\sigma^{\prime}=\pm}E^{*}_{\widehat{p}\sigma^{\prime}}(X)\,E_{\widehat{p}\sigma}(X)\times\\\
&\times{\text{tr}}\left[\Delta(\sigma^{\prime})\,\Delta(\sigma)\,\gamma_{\mu}^{\parallel}\,\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right].\end{split}$

We can simplify the previous expression by taking advantage of the following
identity

$\Delta(\sigma)\Delta(\sigma^{\prime})=\frac{1+\sigma\sigma^{\prime}+{\rm
i}(\sigma+\sigma^{\prime})\gamma_{1}\gamma_{2}}{4}=\delta_{\sigma,\sigma^{\prime}}\Delta(\sigma).$

Notice that the dependence on $p_{2}$ is only contained inside
$E^{*}_{\widehat{p}\sigma^{\prime}}(X)E_{\widehat{p}\sigma}(X)$. We can then
show that the integration on $p_{2}$ gives

$\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{2}}{2\pi}E^{*}_{\widehat{p}\sigma^{\prime}}(X)E_{\widehat{p}\sigma}(X)=|qB|\delta_{n,n^{\prime}}.$

Furthermore, it is convenient to split the sum on $l$ between the lowest
Landau level (LLL) $l=0$ and the higher Landau levels (HLL) $l>1$. For $l=0$
the sums on $\sigma$ are constrained to be equal to
$\sigma=\sigma^{\prime}=\bar{\sigma}=$sgn$(eB)$, and the momenta are given by
$\underline{P}^{+}_{\chi}=(\omega_{n}+{\rm i}\mu_{\chi},0,0,p_{3})$; then at
the lowest Landau level we have

$\langle\widehat{j}^{\,\parallel}_{\mu}(X)\rangle_{\textup{LLL}}=-(-{\rm
i})^{1-\delta_{\mu,0}}q|qB|\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}{\text{tr}}\left[\frac{1+{\rm
i}\bar{\sigma}\gamma_{1}\gamma_{2}}{2}\gamma_{\mu}^{\parallel}\,\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right].$

After computing the trace, we find that the zero component is

$\begin{split}\langle\widehat{j}_{0}(X)\rangle_{\textup{LLL}}=&-\int_{-\infty}^{\infty}\frac{q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\frac{-{\rm
i}\left[(\omega_{n}+{\rm i}\mu_{\chi})-{\rm
i}p_{3}\bar{\sigma}\chi\right]}{(\omega_{n}+{\rm
i}\mu_{\chi})^{2}+p_{3}^{2}}\\\
=&\int_{-\infty}^{\infty}\frac{q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\frac{{\rm
i}\left[(\omega_{n}+{\rm i}\mu_{\chi})\right]}{(\omega_{n}+{\rm
i}\mu_{\chi})^{2}+p_{3}^{2}},\end{split}$

while the $z$ component is

$\begin{split}\langle\widehat{j}_{3}(X)\rangle_{\textup{LLL}}=&{\rm
i}\int_{-\infty}^{\infty}\frac{q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\frac{-{\rm
i}\left[p_{3}+{\rm i}(\omega_{n}+{\rm
i}\mu_{\chi})\bar{\sigma}\chi\right]}{(\omega_{n}+{\rm
i}\mu_{\chi})^{2}+p_{3}^{2}}\\\
=&\int_{-\infty}^{\infty}\frac{\bar{\sigma}q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\frac{\chi{\rm
i}(\omega_{n}+{\rm i}\mu_{\chi})}{(\omega_{n}+{\rm
i}\mu_{\chi})^{2}+p_{3}^{2}},\end{split}$

where the linear terms in $p_{3}$ were dropped because they are odd on $p_{3}$
and as such they vanish when integrated. After the Matsubara sum we have

$\displaystyle\langle\widehat{j}_{0}(X)\rangle_{\textup{LLL}}=q|qB|\sum_{\chi}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{(2\pi)^{2}}\frac{1}{2}\left[n_{\textup{F}}(p_{3}-\mu_{\chi})-n_{\textup{F}}(p_{3}+\mu_{\chi})\right],$
$\displaystyle\langle\widehat{j}_{3}(X)\rangle_{\textup{LLL}}=q^{2}B\sum_{\chi}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{(2\pi)^{2}}\frac{\chi}{2}\left[n_{\textup{F}}(p_{3}-\mu_{\chi})-n_{\textup{F}}(p_{3}+\mu_{\chi})\right].$

Finally, taking advantage of the integral in $p_{3}$:

$\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{2}\left[n_{\textup{F}}(p_{3}-\mu_{\chi})-n_{\textup{F}}(p_{3}+\mu_{\chi})\right]=\mu_{\chi}$

and summing on chiralities, we obtain

$\displaystyle\langle\widehat{j}_{0}(X)\rangle_{\textup{LLL}}=\frac{q|qB|}{(2\pi)^{2}}(\mu_{\textup{R}}+\mu_{\textup{L}})=\frac{\mu
q|qB|}{2\pi^{2}},$
$\displaystyle\langle\widehat{j}_{3}(X)\rangle_{\textup{LLL}}=\frac{q^{2}B}{(2\pi)^{2}}(\mu_{\textup{R}}-\mu_{\textup{L}})=\frac{q^{2}\mu_{\textup{A}}}{2\pi^{2}}B.$

Moving on now to the higher Landau levels, consider

$\begin{split}\langle\widehat{j}^{\,\parallel}_{\mu}(X)\rangle_{\textup{HLL}}=&-(-{\rm
i})^{1-\delta_{\mu,0}}q|qB|\sum_{l=1}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\sum_{\sigma,\sigma^{\prime}=\pm}\delta_{n,n^{\prime}}\times\\\
&\times{\text{tr}}\left[\Delta(\sigma^{\prime})\,\Delta(\sigma)\,\gamma_{\mu}^{\parallel}\,\mathbb{P}_{\chi}\frac{-{\rm
i}\not{\underline{P}}_{\chi}^{+}}{{\underline{P}^{+}_{\chi}}^{2}}\right].\end{split}$

When $l$ is fixed we can replace the $\delta_{n,n^{\prime}}$ with the
$\delta_{\sigma,\sigma^{\prime}}$ and the sum on $\sigma^{\prime}$ becomes
straightforward. The expression is similar to the LLL case, we just have to
replace $\bar{\sigma}$ with $\sigma$ and sum over $\sigma=\pm$. Remind that
now $\underline{P}$ has also a $y$ component. After evaluating the trace and
removing $p_{3}$ odd terms, we obtain

$\begin{split}\langle\widehat{j}_{0}(X)\rangle_{\textup{HLL}}=&\sum_{l=1}^{\infty}\int_{-\infty}^{\infty}\frac{q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\frac{2{\rm
i}\left[(\omega_{n}+{\rm
i}\mu_{\chi})\right]}{{\underline{P}^{+}_{\chi}}^{2}},\\\
\langle\widehat{j}_{3}(X)\rangle_{\textup{HLL}}=&\sum_{l=1}^{\infty}\int_{-\infty}^{\infty}\frac{q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}T\sum_{\\{\omega_{n}\\}}\sum_{\sigma=\pm}\sigma\frac{\chi{\rm
i}(\omega_{n}+{\rm
i}\mu_{\chi})+p_{3}}{{\underline{P}^{+}_{\chi}}^{2}}=0.\end{split}$

We found that the third component does not get corrections from HLL. Instead
for the time component, after the frequency sum, we obtain

$\begin{split}\langle\widehat{j}_{0}(X)\rangle_{\textup{HLL}}&=\sum_{l=1}^{\infty}\int_{-\infty}^{\infty}\frac{q|qB|\mathrm{d}p_{3}}{(2\pi)^{2}}\sum_{\chi}\left[n_{\textup{F}}(E_{p_{3},l}-\mu_{\chi})-n_{\textup{F}}(E_{p_{3},l}+\mu_{\chi})\right],\end{split}$

where $E_{p_{3},l}\equiv\sqrt{p_{3}^{2}+2|qB|l}$. For the perpendicular
components ($\widehat{j}_{x}$ and $\widehat{j}_{y}$) we expect a vanishing
result because they are not allowed by the symmetries of the system. Indeed
the explicit calculations confirmed this expectation.

In summary, restoring covariant expression, we found that the electric current
has two thermodynamic function: the electric charge density $n_{\textup{c}}$
and the chiral magnetic effect (CME) conductivity $\sigma_{\textup{B}}$:

$\langle\widehat{j}_{\mu}(X)\rangle=n_{\textup{c}}\,u_{\mu}+\sigma_{\textup{B}}\,B_{\mu}.$

The electric charge density is given by the mean value
$\langle\widehat{j}_{0}(X)\rangle$ at the local rest frame:

$\begin{split}n_{\textup{c}}=\frac{q|qB|}{2\pi^{2}}\Bigg{\\{}\mu+\sum_{l=1}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{2}&\Big{[}n_{\textup{F}}(E_{p_{3},l}-\mu_{\textup{R}})-n_{\textup{F}}(E_{p_{3},l}+\mu_{\textup{R}})+\\\
&\hphantom{\Big{[}}+n_{\textup{F}}(E_{p_{3},l}-\mu_{\textup{L}})-n_{\textup{F}}(E_{p_{3},l}+\mu_{\textup{L}})\Big{]}\Bigg{\\}},\end{split}$

while the CME conductivity is given by $\langle\widehat{j}_{3}(X)\rangle/B$,
that is

$\sigma_{\textup{B}}=\frac{q^{2}\mu_{\textup{A}}}{2\pi^{2}}.$ (37)

To our knowledge the equation for the electric charge density of an electron
gas in a magnetic medium was first given in Canuto:1969cs  and coincides with
the expression above. The CME effect evaluated here coincides with the one
obtained with many other derivations Kharzeev:2015znc , however we want to
point out that this derivation is valid at thermal equilibrium, as the one in
Vilenkin:1980fu , and that is non-perturbative in magnetic field.

###  4.5 Axial current mean value

We can compute the axial current mean value exactly as described above for the
electric current. Because of that we omit all the calculations. The axial
current constitutive equation is written in terms of an axial charge density
$n_{\textup{A}}$ and a Chiral Separation Effect (CSE) conductivity
$\sigma_{\textup{s}}$:

$\langle\widehat{j}_{A\mu}\rangle=n_{\textup{A}}\,u_{\mu}+\sigma_{\textup{s}}B_{\mu}.$

In this case too we found that only the lowest Landau level contributes to CSE
and that the final result is

$\begin{split}n_{\textup{A}}=&\frac{|qB|}{2\pi^{2}}\Bigg{\\{}\mu_{\textup{A}}+\sum_{l=1}^{\infty}\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{2}\left[n_{\textup{F}}(E_{p_{3},l}-\mu_{\textup{R}})-n_{\textup{F}}(E_{p_{3},l}+\mu_{\textup{R}})+\right.\\\
&\left.-n_{\textup{F}}(E_{p_{3},l}-\mu_{\textup{L}})+n_{\textup{F}}(E_{p_{3},l}+\mu_{\textup{L}})\right]\Bigg{\\}},\\\
\sigma_{\textup{s}}=&\frac{q\mu}{2\pi^{2}},\end{split}$

with $E_{p_{3},l}=\sqrt{p_{3}^{2}+2|qB|l}$. The last thermal coefficient is
exactly the well known value of Chiral Separation Effect (CSE) conductivity
Kharzeev:2015znc .

The same procedure can be followed to evaluate the axial charge density and
the CSE conductivity of massive fermions with vanishing axial chemical
potential $\mu_{\textup{A}}=0$. In that case, as discussed previously, the
thermal equilibrium can be reached and all the quantity discussed in this
section are still well defined. The results for the thermodynamic functions
related to axial current are:

$\displaystyle n_{\textup{A}}=$ $\displaystyle 0,$
$\displaystyle\sigma_{\textup{s}}=$ $\displaystyle
qB\int_{-\infty}^{\infty}\frac{\mathrm{d}p_{3}}{(2\pi)^{2}}\left[n_{\textup{F}}(E_{p_{3}}-\mu)-n_{\textup{F}}(E_{p_{3}}+\mu)\right],$
(38)

with $E_{p_{3}}^{2}=p_{3}^{2}+m^{2}$. The CSE induces an axial current even if
the system is not chiral. It is apparent from the result above that CSE has an
explicit mass dependence, as it is known that it should have Metlitski:2005pr
.

##  5 Constant vorticity and electromagnetic field

So far this contribution has focused on the general properties of the
statistical operator of global thermodynamic equilibrium with both vorticity
and electromagnetic field (12). The following section will discuss the special
case of a constant homogeneous electromagnetic field ($F^{\mu\nu}=$constant)
for which we already studied the symmetries (Sec. 3.1). As discussed in Sec.
2, global equilibrium can only be reached if condition (9) is satisfied, which
in this case becomes

$\mathcal{L}_{\beta}(F^{\mu\nu})=\varpi^{\mu}_{\,\,\sigma}F^{\sigma\nu}-\varpi^{\nu}_{\,\,\sigma}F^{\sigma\mu}\equiv(\varpi\wedge
F)^{\mu\nu}=0.$ (39)

We already discussed this wedge product in Eq. (28). The Eq. (39) has two
independent solutions: $F=k\varpi$ and $F=k^{\prime}\varpi^{*}$, with $k$ and
$k^{\prime}$ real numbers. In terms of the gauge potential the condition (10)
must be satisfied. From Eq. (39), choosing the covariant gauge
$A^{\mu}=\frac{1}{2}F^{\rho\mu}x_{\rho}$, we find that condition (10) is
satisfied setting $\Phi=\frac{1}{2}b_{\sigma}F^{\sigma\lambda}x_{\lambda}$.
The equilibrium chemical potential (11) is then written as

$\zeta(x)=\zeta_{0}-\beta_{\sigma}(x)F^{\lambda\sigma}x_{\lambda}+\frac{1}{2}\varpi_{\sigma\rho}x^{\rho}F^{\lambda\sigma}x_{\lambda}.$
(40)

The same solution can also be obtained by directly solving Eq. (7) using Eq.
(39). This last method to obtain the solution explicitly shows that the
chemical potential in Eq. (40) is not gauge dependent. For constant magnetic
field and vanishing thermal vorticity the solution (40) reduces to
$\zeta=$constant, as it was correctly used in Sec. 4.

Plugging the form (40) inside the operator of Eq. (1) we find:

$\widehat{\rho}=\frac{1}{Z}\exp\left\\{-\int\mathrm{d}\Sigma_{\lambda}\left[\left(\widehat{T}^{\lambda\nu}-\widehat{j}^{\,\lambda}F^{\nu\rho}x_{\rho}\right)\beta_{\nu}-\frac{1}{2}\varpi_{\sigma\rho}\widehat{j}^{\,\lambda}x^{\rho}F^{\tau\sigma}x_{\tau}-\zeta_{0}\widehat{j}^{\,\lambda}\right]\right\\}.$

Inside the round bracket we recognize the divergence-less operator
$\widehat{\pi}^{\lambda\nu}$ of Eq. (27), whose integrals are the generators
of translations. Expressing the coordinate dependence of $\beta$, we can then
write

$\begin{split}\widehat{\rho}=&\frac{1}{Z}\exp\left\\{-\int\mathrm{d}\Sigma_{\lambda}\left[\widehat{\pi}^{\lambda\nu}b_{\nu}+\varpi_{\nu\tau}x^{\tau}\widehat{\pi}^{\lambda\nu}-\frac{1}{2}\varpi_{\mu\nu}\widehat{j}^{\,\lambda}x^{\nu}F^{\tau\mu}x_{\tau}-\zeta_{0}\widehat{j}^{\,\lambda}\right]\right\\}\\\
=&\frac{1}{Z}\exp\left\\{-\int\mathrm{d}\Sigma_{\lambda}\left[\widehat{\pi}^{\lambda\nu}b_{\nu}+\varpi_{\mu\nu}x^{\nu}\left(\widehat{\pi}^{\lambda\mu}-\frac{1}{2}\widehat{j}^{\,\lambda}F^{\rho\mu}x_{\rho}\right)-\zeta_{0}\widehat{j}^{\,\lambda}\right]\right\\}\\\
=&\frac{1}{Z}\exp\left\\{-\int\mathrm{d}\Sigma_{\lambda}\left[\widehat{\pi}^{\lambda\nu}b_{\nu}-\frac{1}{2}\varpi_{\mu\nu}\left[x^{\mu}\left(\widehat{\pi}^{\lambda\nu}-\frac{1}{2}\widehat{j}^{\,\lambda}F^{\rho\nu}x_{\rho}\right)\right.\right.\right.\\\
&\left.\left.\left.-x^{\nu}\left(\widehat{\pi}^{\lambda\mu}-\frac{1}{2}\widehat{j}^{\,\lambda}F^{\rho\mu}x_{\rho}\right)\right]-\zeta_{0}\widehat{j}^{\,\lambda}\right]\right\\};\end{split}$

this time we have recreated the divergence-less quantity
$\varpi_{\mu\nu}\widehat{M}^{\lambda,\mu\nu}$ of Eq. (30) that generates the
Lorentz transformations and that are symmetries of the system. We can then
integrate over the coordinate and we find:

$\widehat{\rho}=\frac{1}{Z}\exp\left\\{-b\cdot\widehat{\pi}+\frac{1}{2}\varpi:\widehat{M}+\zeta_{0}\widehat{Q}\right\\}.$

In the above form, the analogy with statistical operator without
electromagnetic field in Eq. (6) is evident. In both cases the statistical
operator is written with the sum of conserved operators, each one weighted
with a constant Lagrange multiplier. Moreover, starting from a fixed point $x$
we can write the constants thermal fields as

$b_{\mu}=\beta(x)_{\mu}-\varpi_{\mu\nu}x^{\nu},\quad\zeta_{0}=\zeta(x)+\beta_{\sigma}(x)F^{\lambda\sigma}x_{\lambda}-\frac{1}{2}\varpi_{\sigma\rho}x^{\rho}F^{\lambda\sigma}x_{\lambda},$

from which the statistical operator becomes

$\begin{split}\widehat{\rho}=&\frac{1}{Z}\exp\Big{\\{}-\beta(x)_{\mu}\left(\widehat{\pi}^{\mu}-F^{\lambda\mu}x_{\lambda}\widehat{Q}\right)+\\\
&+\frac{1}{2}\varpi_{\mu\nu}\left(\widehat{M}^{\mu\nu}+x^{\nu}\widehat{\pi}^{\mu}-x^{\mu}\widehat{\pi}^{\nu}-x^{\nu}F^{\lambda\mu}x_{\lambda}\widehat{Q}\right)+\zeta(x)\widehat{Q}\Big{\\}}.\end{split}$

It is important to point out that with an external magnetic field the Poincaré
algebra is modified and becomes the Algebra in Eq. (31), which we report here
for convenience:

$\begin{split}[\widehat{\pi}^{\mu},\widehat{\pi}^{\nu}]=&{\rm
i}F^{\mu\nu}\widehat{Q},\\\
\frac{1}{2}F_{\rho\sigma}[\widehat{\pi}^{\mu},\widehat{M}^{\rho\sigma}]=&\frac{{\rm
i}}{2}F_{\rho\sigma}\left(\eta^{\mu\rho}\widehat{\pi}^{\sigma}-\eta^{\mu\sigma}\widehat{\pi}^{\rho}\right),\\\
\frac{1}{2}F^{*}_{\rho\sigma}[\widehat{\pi}^{\mu},\widehat{M}^{\rho\sigma}]=&\frac{{\rm
i}}{2}F^{*}_{\rho\sigma}\left(\eta^{\mu\rho}\widehat{\pi}^{\sigma}-\eta^{\mu\sigma}\widehat{\pi}^{\rho}\right).\end{split}$

Notice that because $F$ is proportional to $\varpi$, if we replace $F$ with
$\varpi$ and $F^{*}$ with $\varpi^{*}$ the last two algebra identities still
hold true. Since the Algebra is known, we can translate the statistical
operator. Taking advantage of the unitary of the translation transformation,
the translated statistical operator is

$\begin{split}\widehat{\sf T}(x)\,\widehat{\rho}\,\widehat{\sf
T}^{-1}(x)=&\frac{1}{Z}\exp\Big{\\{}-\widehat{\sf
T}(x)\left(b\cdot\widehat{\pi}\right)\widehat{\sf T}^{-1}(x)+\\\
&+\widehat{\sf T}(x)\left(\frac{\varpi:\widehat{M}}{2}\right)\widehat{\sf
T}^{-1}(x)+\zeta_{0}\widehat{\sf T}(x)\,\widehat{Q}\,\widehat{\sf
T}^{-1}(x)\Big{\\}}.\end{split}$

Therefore we just need to evaluate how the operators
$\widehat{\pi},\,\widehat{M}$ and $\widehat{Q}$ transform under translations.
For a unitary transformation an operator $\widehat{K}$ transforms with

$\mathrm{e}^{{\rm i}\widehat{A}}\widehat{K}\mathrm{e}^{-{\rm
i}\widehat{A}}\simeq\widehat{K}-{\rm
i}\left[\widehat{K},\widehat{A}\right]-\frac{1}{2}\left[\left[\widehat{K},\widehat{A}\right],\widehat{A}\right]+\frac{{\rm
i}}{6}\left[\left[\left[\widehat{K},\widehat{A}\right],\widehat{A}\right],\widehat{A}\right]+\cdots.$

By applying this formula to our operators we obtain the complete
transformation because after a certain order all the commutators become
vanishing. In particular, for the Lorentz transformation generators we find

$\frac{1}{2}\varpi_{\mu\nu}\widehat{M}^{\mu\nu}_{x}\equiv\frac{1}{2}\varpi_{\mu\nu}\widehat{\sf
T}(x)\widehat{M}^{\mu\nu}\widehat{\sf
T}^{-1}(x)=\frac{1}{2}\varpi_{\mu\nu}\left(\widehat{M}^{\mu\nu}+x^{\nu}\widehat{\pi}^{\mu}-x^{\mu}\widehat{\pi}^{\nu}-x^{\nu}F^{\lambda\mu}x_{\lambda}\widehat{Q}\right).$

For the other operators instead we find:

$\widehat{\pi}^{\mu}_{x}\equiv\widehat{\sf
T}(x)\,\widehat{\pi}^{\mu}\,\widehat{\sf
T}^{-1}(x)=\widehat{\pi}^{\mu}-x_{\rho}F^{\rho\mu}\widehat{Q},\quad\widehat{\sf
T}(x)\,\widehat{Q}\,\widehat{\sf T}^{-1}(x)=\widehat{Q}.$

With these definitions a translation transformation on the statistical
operator act as following:

$\begin{split}\widehat{\sf T}(a)\,\widehat{\rho}\,\widehat{\sf
T}^{-1}(a)=&\frac{1}{Z}\exp\left\\{-\beta(x)\cdot\widehat{\pi}_{x+a}+\frac{1}{2}\varpi:\widehat{M}_{x+a}+\zeta(x)\widehat{Q}\right\\}\\\
=&\frac{1}{Z}\exp\left\\{-\beta(x-a)\cdot\widehat{\pi}_{x}+\frac{1}{2}\varpi:\widehat{M}_{x}+\zeta(x-a)\widehat{Q}\right\\}.\end{split}$

It follows that the statistical operator around a point $x$ can be written as:

$\widehat{\rho}=\frac{1}{Z}\exp\left\\{-\beta(x)\cdot\widehat{\pi}_{x}+\frac{1}{2}\varpi:\widehat{M}_{x}+\zeta(x)\widehat{Q}\right\\}.$
(41)

###  5.1 Expansion on thermal vorticity

Following ref. Buzzegoli:2018wpy  we use linear response theory to evaluate
thermal expectation values in the case of constant electromagnetic field. The
purpose of this section is to give the thermal expectation value of an
operator $\widehat{O}$ at the point $x$ as a thermal vorticity expansion.
Using the properties of the trace, we can transfer the $x$ dependence from the
operator $\widehat{O}$ to the statistical operator (41) written around the
same point $x$:

$\begin{split}\langle\widehat{O}(x)\rangle=&\frac{1}{Z}{\text{tr}}\left[\exp\left\\{-\beta(x)\cdot\widehat{\pi}_{x}+\frac{1}{2}\varpi:\widehat{M}_{x}+\zeta(x)\widehat{Q}\right\\}\widehat{O}(x)\right]\\\
=&\frac{1}{Z}{\text{tr}}\left[\widehat{\sf
T}(-x)\exp\left\\{-\beta(x)\cdot\widehat{\pi}_{x}+\frac{1}{2}\varpi:\widehat{M}_{x}+\zeta(x)\widehat{Q}\right\\}\widehat{\sf
T}^{-1}(-x)\widehat{O}(0)\right]\\\
=&\frac{1}{Z}{\text{tr}}\left[\exp\left\\{-\beta(x)\cdot\widehat{\pi}+\frac{1}{2}\varpi:\widehat{M}+\zeta(x)\widehat{Q}\right\\}\widehat{O}(0)\right].\end{split}$

To evaluate the mean value we expand the statistical operator of the last
equality around vanishing vorticity. First, we split the exponent of the
statistical operator in two parts as follows:

$\widehat{\rho}=\dfrac{1}{Z}\exp\left[\widehat{A}+\widehat{B}\right],\quad\widehat{A}\equiv-\beta_{\mu}(x)\widehat{\pi}^{\mu}+\zeta(x)\widehat{Q},\quad\widehat{B}\equiv\frac{1}{2}\varpi:\widehat{M},$

and then we expand on $\widehat{B}$, which is the part containing thermal
vorticity. Since $\widehat{B}$ and $\widehat{A}$ satisfy the same algebra of
the case discussed in Buzzegoli:2018wpy , the expansion will leads to the same
result, which is

$\begin{split}\langle\widehat{O}(x)\rangle=&\langle\widehat{O}(0)\rangle_{\beta(x)}-\alpha_{\rho}\langle\langle\,\widehat{K}^{\rho}\widehat{O}\,\rangle\rangle-
w_{\rho}\langle\langle\,\widehat{J}^{\rho}\widehat{O}\,\rangle\rangle+\frac{\alpha_{\rho}\alpha_{\sigma}}{2}\langle\langle\,\widehat{K}^{\rho}\widehat{K}^{\sigma}\widehat{O}\,\rangle\rangle\\\
&+\frac{w_{\rho}w_{\sigma}}{2}\langle\langle\,\widehat{J}^{\rho}\widehat{J}^{\sigma}\widehat{O}\,\rangle\rangle+\frac{\alpha_{\rho}w_{\sigma}}{2}\langle\langle\,\\{\widehat{K}^{\rho},\widehat{J}^{\sigma}\\}\widehat{O}\,\rangle\rangle+\mathcal{O}(\varpi^{3}),\end{split}$
(42)

where we defined

$\begin{split}\langle\langle\widehat{K}^{\rho_{1}}\cdots\widehat{K}^{\rho_{n}}\widehat{J}^{\sigma_{1}}\cdots\widehat{J}^{\sigma_{m}}&\widehat{O}\rangle\rangle\equiv\int_{0}^{|\beta|}\frac{\mathrm{d}\tau_{1}\cdots\mathrm{d}\tau_{n+m}}{|\beta|^{n+m}}\times\\\
&\times\langle{\rm T}_{\tau}\left(\widehat{K}^{\rho_{1}}_{-{\rm
i}\tau_{1}u}\cdots\widehat{K}^{\rho_{n}}_{-{\rm
i}\tau_{n}u}\widehat{J}^{\sigma_{1}}_{-{\rm
i}\tau_{n+1}u}\cdots\widehat{J}^{\sigma_{m}}_{-{\rm
i}\tau_{n+m}u}\widehat{O}(0)\right)\rangle_{\beta(x),c}.\end{split}$

As discussed in Sec. 3.1, the boost and rotation defined starting from
$\widehat{M}^{\mu\nu}$, i.e.

$\widehat{K}^{\mu}=u_{\lambda}\widehat{M}^{\lambda\mu},\quad\widehat{J}^{\mu}=\frac{1}{2}\epsilon^{\alpha\beta\gamma\mu}u_{\alpha}\widehat{M}_{\beta\gamma},$

are different from those of a system without external electromagnetic field.
The other difference with Buzzegoli:2018wpy  is that in this case the averages
$\langle\cdots\rangle_{\beta(x)}$ are made with the statistical operator

$\widehat{\rho}_{0}=\frac{1}{Z_{0}}\exp\left\\{-\beta(x)\cdot\widehat{\pi}+\zeta(x)\widehat{Q}\right\\}.$

###  5.2 Currents and chiral anomaly

In this section we determine the constitutive equations for the electric and
the axial current at first order in thermal vorticity and we investigate the
contributions from electric and magnetic fields to the thermal coefficients
related to vorticity. The constitutive equations are obtained from the
expansion on thermal vorticity given in the previous section. Instead of a
direct evaluation we use the conservation equations to show that indeed no
additional corrections from electric and magnetic field occur to first order
vorticous coefficients, such as the conductivity of the Chiral Vortical Effect
(CVE) and of the Axial Vortical Effect (AVE). In this way we obtain several
relations between those coefficients and their relation to the chiral anomaly.

Consider the case of global thermal equilibrium with constant vorticity
$\varpi_{\mu\nu}$ and an electromagnetic field with strength tensor
$F_{\mu\nu}=k\,\varpi_{\mu\nu}$, with $k$ a constant. It then follows that the
comoving magnetic and electric fields are parallel respectively to thermal
rotation and thermal acceleration:

$B^{\mu}(x)=-k\,w^{\mu}(x),\quad E^{\mu}(x)=k\,\alpha^{\mu}(x).$

For instance, in the case of a constant thermal vorticity caused by a rigid
rotation along the $z$ axis and a constant magnetic field along $z$, we have:

$\varpi_{\mu\nu}=\frac{\Omega}{T_{0}}\left(\eta_{\mu 1}\eta_{\nu 2}-\eta_{\nu
1}\eta_{\mu 2}\right),\quad F_{\mu\nu}=B\left(\eta_{\mu 1}\eta_{\nu
2}-\eta_{\nu 1}\eta_{\mu 2}\right),\quad k=\frac{BT_{0}}{\Omega},$

where $\Omega,\,T_{0}$ and $B$ are constants. In this example, electric and
magnetic fields are orthogonal and there is no chiral anomaly. However, in the
general case, the product $E\cdot B$ is non-vanishing. In that case, as we
showed in Sec. 2, we can still discuss global equilibrium with chiral
imbalance by defining a conserved Chern-Simons current.

By using the thermal vorticity expansion (42) we now proceed to write the
thermal expectation value of electric current at first order in thermal
vorticity. We want to stress that in the expansion (42) no approximations are
made on the effects of the external electric and magnetic fields; the
expansion (42) only approximates the effects of vorticity. At first order on
thermal vorticity the only quantities that can contribute to the mean value of
a current are the four-vectors $w^{\mu},\,\alpha^{\mu}$ and the scalars
$E\cdot\alpha$, $E\cdot w=-B\cdot\alpha$ and $B\cdot w$. We therefore write
the thermal expansion in terms of these quantities, which will define several
thermal coefficients. Taking into account the symmetries, the thermal
vorticity expansion of the electric current is

$\begin{split}\langle\,\widehat{j}^{\,\mu}(x)\rangle=&\left[n_{\textup{c}}^{0}+n_{\textup{c}}^{E\cdot\alpha}(E\cdot\alpha)+n_{\textup{c}}^{B\cdot
w}(B\cdot
w)\right]u^{\mu}+W^{\textup{V}}w^{\mu}+\sigma_{E}^{B\cdot\alpha}(B\cdot\alpha)E^{\mu}\\\
&+\left[\sigma_{B}^{0}+\sigma_{B}^{E\cdot\alpha}(E\cdot\alpha)+\sigma_{B}^{B\cdot
w}(B\cdot w)\right]B^{\mu}+\mathcal{O}\left(\varpi^{2}\right).\end{split}$
(43)

Since the thermal coefficients $n_{\textup{c}}^{0}$ and $\sigma_{B}^{0}$ must
be evaluated at vanishing thermal vorticity, they are exactly those computed
in Sec. 4.4 (for vanishing electric field). In particular, the Chiral Magnetic
Effect (CME) conductivity at vanishing vorticity, Eq. (37), is

$\sigma_{B}^{0}(x)=\frac{q\zeta_{A}}{2\pi^{2}|\beta(x)|}.$ (44)

All the other coefficients are related to thermal vorticity and they have the
following properties under parity, time-reversal and charge conjugation

$\begin{array}[]{lccc|cccc|ccc}&E\cdot\alpha&E\cdot w&B\cdot
w&n_{\textup{c}}^{0}&\sigma_{B}^{0}&W^{\textup{V}}&\sigma_{E}^{0}&\sigma_{B}^{E\cdot\alpha}&\sigma_{B}^{B\cdot
w}&\sigma_{E}^{B\cdot\alpha}\\\ \hline\cr\mathsf{P}&+&-&+&+&-&-&+&-&-&-\\\
\mathsf{T}&+&-&+&+&+&+&-&+&+&+\\\ \mathsf{C}&-&-&-&-&+&-&+&-&-&-\\\
\end{array}$ (45)

Similarly the axial current thermal expectation value is

$\begin{split}\langle\,\widehat{j}^{\,\mu}_{\textup{A}}(x)\rangle=&\left[n_{\textup{A}}^{0}+n_{\textup{A}}^{E\cdot\alpha}(E\cdot\alpha)+n_{\textup{A}}^{B\cdot
w}(B\cdot
w)\right]u^{\mu}+W^{\textup{A}}w^{\mu}+\sigma_{sE}^{B\cdot\alpha}(B\cdot\alpha)E^{\mu}\\\
&+\left[\sigma_{s}^{0}+\sigma_{s}^{E\cdot\alpha}(E\cdot\alpha)+\sigma_{s}^{B\cdot
w}(B\cdot w)\right]B^{\mu}+\mathcal{O}\left(\varpi^{2}\right).\end{split}$

Each thermal coefficients is a function depending only on

$|\beta|,\,\zeta,\,\zeta_{\textup{A}},\,B^{2},\,E^{2},\,E\cdot B.$ (46)

The coordinate dependence of any thermal coefficients is completely contained
inside the Lorentz scalars in (46).

With the constitutive equations written down, we are now looking for relations
and constraints between those thermodynamic coefficients. The conservation of
electric current implies that

$\partial_{\mu}\langle\,\widehat{j}^{\,\mu}(x)\rangle=\langle\,\partial_{\mu}\widehat{j}^{\,\mu}(x)\rangle=0.$

The coordinate derivative acts both on thermal coefficients and on
thermodynamic fields. We need to establish how the derivative acts on those
quantities. For thermodynamic fields, using the equilibrium conditions and the
identities in Appendix Appendix: Thermodynamic relations in beta frame, we
find

$\begin{split}\partial_{\mu}u^{\mu}=&0,\quad\partial_{\mu}w^{\mu}=-3\frac{w\cdot\alpha}{|\beta|},\quad\partial_{\mu}\alpha^{\mu}=\frac{2w^{2}-\alpha^{2}}{|\beta|},\\\
\partial_{\mu}B^{\mu}=&-3\frac{B\cdot\alpha}{|\beta|},\quad\partial_{\mu}E^{\mu}=-\frac{2(w\cdot
B)+(\alpha\cdot E)}{|\beta|},\quad\partial_{\mu}(B\cdot\alpha)=0,\\\
\partial_{\mu}(E\cdot\alpha)=&-\frac{2}{|\beta|}\left[(w\cdot
B)\alpha_{\mu}+(E\cdot w)w_{\mu}\right],\\\ \partial_{\mu}(B\cdot
w)=&\frac{2}{|\beta|}\left[(w\cdot B)\alpha_{\mu}+(E\cdot
w)w_{\mu}\right].\end{split}$

Moreover, we can also show that

$\begin{split}\partial_{\mu}|\beta|=-\alpha_{\mu},\quad\partial_{\mu}\zeta=\beta^{\nu}F_{\nu\mu}=-|\beta|E_{\mu},\quad\partial_{\mu}\zeta_{\textup{A}}=0,\quad\partial_{\mu}(E\cdot
B)=0,\\\ \partial_{\mu}|B|=-\frac{(E\cdot
B)w_{\mu}+B^{2}\alpha_{\mu}}{|\beta||B|},\quad\partial_{\mu}|E|=-\frac{(E\cdot
B)w_{\mu}+B^{2}\alpha_{\mu}}{|\beta||E|},\end{split}$

where

$|\beta|=\sqrt{\beta^{\sigma}\beta_{\sigma}},\quad|B|=\sqrt{-B^{\sigma}B_{\sigma}},\quad|E|=\sqrt{-E^{\sigma}E_{\sigma}}.$

The derivative respect to coordinates of a thermodynamic function is

$\begin{split}\partial_{\mu}f(|\beta|,\zeta,\zeta_{A},|B|,|E|,E\cdot
B)=&\left(-\partial_{\mu}|\beta|\frac{\partial}{\partial|\beta|}+\partial_{\mu}\zeta\frac{\partial}{\partial\zeta}+\partial_{\mu}\zeta_{\textup{A}}\frac{\partial}{\partial\zeta_{\textup{A}}}+\partial_{\mu}|B|\frac{\partial}{\partial|B|}\right.\\\
&\left.+\partial_{\mu}|E|\frac{\partial}{\partial|E|}+\partial_{\mu}(E\cdot
B)\frac{\partial}{\partial(E\cdot B)}\right)f.\end{split}$

Therefore, using the previous identities, the derivative of a thermodynamic
function becomes

$\begin{split}\partial_{\mu}f=&\left[-\alpha_{\mu}\left(\frac{\partial}{\partial|\beta|}-\frac{|B|^{2}}{|\beta|}\frac{1}{|B|}\frac{\partial}{\partial|B|}-\frac{|B|^{2}}{|\beta|}\frac{1}{|E|}\frac{\partial}{\partial|E|}\right)\right.\\\
&\left.-|\beta|E_{\mu}\frac{\partial}{\partial\zeta}-\frac{(E\cdot
B)w_{\mu}}{|\beta|}\left(\frac{1}{|B|}\frac{\partial}{\partial|B|}+\frac{1}{|E|}\frac{\partial}{\partial|E|}\right)\right]f.\end{split}$

We can also define the following short-hand notation:

$\partial_{\tilde{\beta}}\equiv\frac{\partial}{\partial|\beta|}-\frac{|B|^{2}}{|\beta|}\partial_{\tilde{B}},\quad\partial_{\tilde{B}}\equiv\frac{1}{|B|}\frac{\partial}{\partial|B|}+\frac{1}{|E|}\frac{\partial}{\partial|E|},$

from which the previous derivative is written as

$\partial_{\mu}f=\left[-\alpha_{\mu}\partial_{\tilde{\beta}}-|\beta|E_{\mu}\partial_{\zeta}-\frac{(E\cdot
B)w_{\mu}}{|\beta|}\partial_{\tilde{B}}\right]f.$

We can now use the previous relations to impose electric current conservation
by evaluating the divergence of the expansion in Eq. (43). For the terms
directed along the fluid velocity we find that no additional constraints are
required:

$\partial_{\mu}\left(n^{0}u^{\mu}\right)=\partial_{\mu}\left(n^{E\cdot\alpha}(E\cdot\alpha)u^{\mu}\right)=\partial_{\mu}\left(n^{B\cdot
w}(B\cdot w)u^{\mu}\right)=0.$

For the terms along the magnetic field we find:

$\begin{split}\partial_{\mu}\left(\sigma_{B}^{0}B^{\mu}\right)=&-(B\cdot\alpha)\left[\frac{3}{|\beta|}+\partial_{\tilde{\beta}}\right]\sigma_{B}^{0}-(E\cdot
B)|\beta|\partial_{\zeta}\sigma_{B}^{0}\\\ &-\frac{(E\cdot B)(B\cdot
w)}{|\beta|}\partial_{\tilde{B}}\sigma_{B}^{0},\\\
\partial_{\mu}\left(\sigma_{B}^{E\cdot\alpha}(E\cdot\alpha)B^{\mu}\right)=&-(B\cdot\alpha)(E\cdot\alpha)\left[\frac{3}{|\beta|}+\partial_{\tilde{\beta}}\right]\sigma_{B}^{E\cdot\alpha}-(E\cdot\alpha)(E\cdot
B)|\beta|\partial_{\zeta}\sigma_{B}^{E\cdot\alpha}\\\ &-\frac{(E\cdot
B)(B\cdot
w)(E\cdot\alpha)}{|\beta|}\partial_{\tilde{B}}\sigma_{B}^{E\cdot\alpha}\\\
&-\frac{2}{|\beta|}\left[(w\cdot B)(B\cdot\alpha)+(E\cdot w)(B\cdot
w)\right]\sigma_{B}^{E\cdot\alpha},\\\ \partial_{\mu}\left(\sigma_{B}^{B\cdot
w}(B\cdot w)B^{\mu}\right)=&-(B\cdot\alpha)(B\cdot
w)\left[\frac{3}{|\beta|}+\partial_{\tilde{\beta}}\right]\sigma_{B}^{B\cdot
w}-(B\cdot w)(E\cdot B)|\beta|\partial_{\zeta}\sigma_{B}^{B\cdot w}\\\
&-\frac{(E\cdot B)(B\cdot
w)^{2}}{|\beta|}\partial_{\tilde{B}}\sigma_{B}^{B\cdot w}\\\
&+\frac{2}{|\beta|}\left[(w\cdot B)(B\cdot\alpha)+(E\cdot w)(B\cdot
w)\right]\sigma_{B}^{B\cdot w}.\end{split}$

Along electric field we have

$\begin{split}\partial_{\mu}\left(\sigma_{E}(B\cdot\alpha)E^{\mu}\right)=&-(B\cdot\alpha)(\alpha\cdot
E)\left[\frac{1}{|\beta|}+\partial_{\tilde{\beta}}\right]\sigma_{E}-(B\cdot\alpha)E^{2}|\beta|\partial_{\zeta}\sigma_{E}\\\
&-\frac{(E\cdot B)(E\cdot
w)}{|\beta|}\partial_{\tilde{B}}\sigma_{E}.\end{split}$

Lastly, the divergence of the term along rotation is

$\partial_{\mu}\left(W^{\textup{V}}w^{\mu}\right)=-(w\cdot\alpha)\left[\frac{3}{|\beta|}+\partial_{\tilde{\beta}}\right]W^{\textup{V}}-(E\cdot
w)|\beta|\partial_{\zeta}W^{\textup{V}}-\frac{(E\cdot
B)w^{2}}{|\beta|}\partial_{\tilde{B}}W^{\textup{V}}.$

To impose that $\partial_{\mu}\langle\,\widehat{j}^{\,\mu}(x)\rangle=0$ we sum
all the previous pieces and we split between the linear independent terms.
Those terms must vanish independently of the values of the electromagnetic
field and of the thermal vorticity and several equalities are obtained. Among
those, we first consider the following identities:

$\begin{split}\partial_{\zeta}\sigma_{B}^{0}=&0,\quad\partial_{\tilde{B}}W^{\textup{V}}=0,\quad\partial_{\zeta}\sigma_{B}^{E\cdot\alpha}=0,\quad\partial_{\tilde{B}}\sigma_{B}^{E\cdot\alpha}=0,\quad\partial_{\tilde{B}}\sigma_{B}^{B\cdot
w}=0,\\\
\partial_{\zeta}\sigma_{E}^{B\cdot\alpha}=&0,\quad\partial_{\tilde{B}}\sigma_{E}^{B\cdot\alpha}=0.\end{split}$

Notice from the table in (45) that $\sigma_{B}^{E\cdot\alpha}$ and
$\sigma_{E}^{B\cdot\alpha}$ are related to $\mathsf{C}$-odd correlator.
Therefore they must be odd functions of the electric chemical potential
$\zeta$. But the previous constraints require that they do not depend on
$\zeta$, therefore they must be vanishing

$\sigma_{B}^{E\cdot\alpha}=0,\quad\sigma_{E}^{B\cdot\alpha}=0.$

The previous constraints also require that $W^{\textup{V}}$ and
$\sigma_{B}^{B\cdot w}$ do not depend on $|B|$ and $|E|$333Note that to reduce
the numbers of relations, we have indicated electric field and magnetic field
derivative together with one derivative $\partial_{\tilde{B}}$. However,
electric and magnetic fields are independent and each derivative must be
considered independently.. That is to say that the CVE conductivity
$W^{\textup{V}}$ is not affected by electric and magnetic field. Moreover,
electric current conservation also imposes that

$\begin{split}\left(3+|\beta|\partial_{\tilde{\beta}}\right)\sigma_{B}^{0}-|\beta|^{2}\partial_{\zeta}W^{\textup{V}}=&0,\\\
\partial_{\tilde{B}}\sigma_{B}^{0}+|\beta|^{2}\partial_{\zeta}\sigma_{B}^{B\cdot
w}=&0,\\\ \left(3+|\beta|\partial_{\tilde{\beta}}\right)\sigma_{B}^{B\cdot
w}=&0,\\\ \left(3+|\beta|\partial_{\beta}\right)W^{\textup{V}}=&0.\end{split}$

We can replace the known result for the CME conductivity $\sigma_{B}^{0}$ of
Eq. (44) in the previous constraints to find find that

$\begin{split}\partial_{\zeta}W^{\textup{V}}=&\frac{1}{|\beta|^{2}}\left(3+|\beta|\partial_{\tilde{\beta}}\right)\sigma_{B}^{0}=\frac{q\zeta_{\textup{A}}}{\pi^{2}|\beta|^{3}},\\\
\partial_{\zeta}\sigma_{B}^{B\cdot w}=&0,\\\
\left(3+|\beta|\frac{\partial}{\partial\beta}\right)W^{\textup{V}}=&0.\end{split}$

Again, since $\sigma_{B}^{B\cdot w}$ is $\mathsf{C}$-odd it follows form
second equation that it must be vanishing

$\sigma_{B}^{B\cdot w}=0.$

We want to empathize that we have found a relation between CVE and CME
conductivities:

$\partial_{\zeta}W^{\textup{V}}=\frac{1}{|\beta|^{2}}\left(3+|\beta|\partial_{\tilde{\beta}}\right)\sigma_{B}^{0}.$
(47)

The CVE conductivity $W^{\textup{V}}$ in Eq. (22) satisfies this relation. It
is important to notice that Eq. (47) completely determines the CVE
conductivity from the CME one. Indeed, since $W^{\textup{V}}$ is odd under
charge conjugation, by fixing the $\zeta$ part of $W^{\textup{V}}$ we obtain
the entire coefficient. This also implies that the $\zeta$ part of
$W^{\textup{V}}$ is dictated by the chiral anomaly as found in effective field
theories Sadofyev:2010is . Therefore, the CVE inherits all the properties
proved for the CME. For instance, it is known that the CME conductivity is
completely dictated by the chiral anomaly Kharzeev:2013ffa  and that it is
protected from corrections coming from interactions Zakharov:2012vv ;
Feng:2018tpb . Since the relation (47) holds not only for a free theory but
for any microscopic interactions, as long as global thermal equilibrium is
concerned, then also the CVE conductivity is dictated by the chiral anomaly
and it is universal. Despite the CVE can be related to the vector-axial
anomalous term of vector current anomaly Landsteiner:2012kd , Eq. (47) shows
that it can be explained with just the electric charge conservation and the
chiral anomaly.

Let us now move to the axial current. Similar steps can be followed to derive
the constraint equations between the thermal coefficients of axial current. In
this case we must impose the following identities between thermal expectation
values:

$\partial_{\mu}\langle\,\widehat{j}^{\,\mu}_{\textup{A}}(x)\rangle=2m\langle\,{\rm
i}\bar{\Psi}\gamma^{5}\Psi\rangle-\frac{q^{2}(E\cdot B)}{2\pi^{2}},$

where we also added the naive divergence term $2m\langle\,{\rm
i}\bar{\Psi}\gamma^{5}\Psi\rangle$ which is due to the mass of the field. From
symmetries the constitutive equation for the pseudo-scalar is

$\begin{split}\langle\,{\rm i}\bar{\Psi}\gamma^{5}\Psi\rangle=&L^{E\cdot
B}(E\cdot B)+L^{\alpha\cdot w}(\alpha\cdot w)+L^{E\cdot w}(E\cdot w)\\\
&+L^{(E\cdot B)\alpha^{2}}(E\cdot B)\alpha^{2}+L^{(E\cdot B)w^{2}}(E\cdot
B)w^{2}\\\ &+L^{(E\cdot w)(E\cdot\alpha)}(E\cdot w)(E\cdot\alpha)+L^{(E\cdot
w)(B\cdot w)}(E\cdot w)(B\cdot
w)+\mathcal{O}\left(\varpi^{3}\right).\end{split}$

The value of $L^{\alpha\cdot w}$ for free Dirac field has been reported in
(20), and the other coefficients related to the electromagnetic field can be
computed with Ritus method.

Because the axial current is not conserved we find different identities
compared to the previous case of electric current. For instance, the
identities related to the Chiral Separation Effect (CSE) conductivity
$\sigma_{s}^{0}$ and to the AVE conductivity $W^{\textup{A}}$ are

$\begin{split}\left(3+|\beta|\frac{\partial}{\partial\beta}\right)W^{\textup{A}}=&-2mL^{\alpha\cdot
w},\\\
\partial_{\zeta}\sigma_{s}^{0}=&\frac{q^{2}}{2\pi^{2}|\beta|}-2mL^{E\cdot
B},\\\
\partial_{\tilde{B}}\sigma^{0}_{s}=&-|\beta|^{2}\partial_{\zeta}\sigma^{B\cdot
w}_{s},\\\ \partial_{\tilde{B}}W^{\textup{A}}=&-2m|\beta|L^{(E\cdot
B)w^{2}},\\\
\partial_{\zeta}W^{\textup{A}}=&\frac{1}{|\beta|^{2}}\left(3+|\beta|\partial_{\tilde{\beta}}\right)\sigma_{s}^{0}-\frac{2m}{|\beta|}L^{E\cdot
w}.\end{split}$

The first equation has been discussed in Buzzegoli:2018wpy  and in Sec. 2.1.
The second equation is similar to the first: the first term on the r.h.s. is
coming from the chiral anomaly and the second from the naive anomaly.
Therefore for a massive field, as discussed for the AVE, the CSE is not
entirely dictated by the anomaly. It is then not surprising to find a
corrections to the CSE Gorbar:2013upa  and that it is affected by the mass,
see Eq. (38). In the massive case we also expect corrections from the external
electromagnetic field both in the AVE and in the CSE conductivities.

On the other hand, for massless field those constraints becomes

$\begin{split}\partial_{\zeta}\sigma_{s}^{0}=&\frac{q^{2}}{2\pi^{2}|\beta|},\\\
\left(3+|\beta|\frac{\partial}{\partial\beta}\right)W^{\textup{A}}=&0,\quad\partial_{\tilde{B}}W^{\textup{A}}=0,\\\
\partial_{\zeta}W^{\textup{A}}=&\frac{1}{|\beta|^{2}}\left(3+|\beta|\partial_{\beta}\right)\sigma_{s}^{0}.\end{split}$

In this case the CSE conductivity is completely fixed by the chiral anomaly as
it is clear by the first equation and by the fact that $\sigma_{s}^{0}$ must
be an odd function of $\zeta$. For the symmetries of axial current,
$W^{\textup{A}}$ has both terms which depend on $\zeta$ and terms which depend
only on $\zeta_{\textup{A}}$ and $\beta$. All these terms must satisfy the
equations in the second line. In particular we conclude that the AVE is not
affected by the electromagnetic field. Moreover, from the third line we see
that the terms related to $\zeta$ are fixed by the CSE conductivity and
consequently they are dictated by the chiral anomaly. As it is evident from
previous discussion, this only occurs for the massless field.

In summary, by imposing the conservation equation we conclude that, at global
thermodynamic equilibrium with thermal vorticity and constant homogeneous
electromagnetic field, the chiral vortical effect is dictated by the chiral
magnetic effect. Then it is not affected by the mass of the particle, by the
external electromagnetic field or by radiative corrections. For the axial
current this analysis has showed that we need to distinguish between the
massive and the massless case. In the latter case we found that the chiral
anomaly completely fixes the whole chiral separation effect (CSE) but fixes
only the part of axial vortical effect (AVE) conductivity which depends on the
electric chemical potential. We also found that the AVE is not affected by the
external electromagnetic field. For the massive case, despite it exists a
relation between the CSE and the AVE, both of them are affected by the mass of
the field, the external electromagnetic field and by radiative corrections.

###### Acknowledgements.

I carried out part of this work while visiting Stony Brook University (New
York, U.S.A.). I would like to thank F. Becattini, E. Grossi and D. Kharzeev
for stimulating discussions on the subject matter. This research was supported
in part by the Florence University with the fellowship “Polarizzazione nei
fluidi relativistici”.

## Appendix: Thermodynamic relations in beta frame

At global thermal equilibrium with thermal vorticity, thermodynamic fields
satisfy several equilibrium relations which constraints their coordinate
dependence. In the $\beta$-frame we can build several quantities from the
four-vector $\beta$ and thermal vorticity $\varpi$:

$\displaystyle
u_{\mu}=\frac{\beta_{\mu}}{\sqrt{\beta^{2}}};\quad\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu};\quad\varpi_{\mu\nu}=\partial_{\nu}\beta_{\mu}=\epsilon_{\mu\nu\rho\sigma}w^{\rho}u^{\sigma}+\alpha_{\mu}u_{\nu}-\alpha_{\nu}u_{\mu};$
$\displaystyle\alpha_{\mu}=\varpi_{\mu\nu}u^{\nu};\quad
w_{\mu}=-\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}\varpi^{\nu\rho}u^{\sigma};\quad\gamma_{\mu}=(\alpha\cdot\varpi)^{\lambda}\Delta_{\lambda\mu}=\epsilon_{\mu\nu\rho\sigma}w^{\nu}\alpha^{\rho}u^{\sigma}.$

Most of these quantities depend on coordinates and their derivatives are
Buzzegoli:2017cqy :

$\displaystyle\partial_{\nu}\beta_{\mu}=\varpi_{\mu\nu};\quad\partial_{\nu}=-\alpha_{\nu}\frac{\partial}{\partial\sqrt{\beta^{2}}};\quad\varpi:\varpi=2\left(\alpha^{2}-w^{2}\right)$
$\displaystyle\partial_{\nu}u_{\mu}=\frac{1}{\sqrt{\beta^{2}}}\big{(}\varpi_{\mu\nu}+\alpha_{\nu}u_{\mu}\big{)};\quad\partial^{\alpha}u_{\alpha}=0;\quad
u_{\alpha}\partial^{\alpha}u_{\mu}=\frac{\alpha_{\mu}}{\sqrt{\beta^{2}}};$
$\displaystyle\partial_{\mu}\alpha_{\nu}=\frac{1}{\sqrt{\beta^{2}}}\big{(}\varpi_{\nu\rho}\varpi^{\rho}_{\,\,\mu}+\alpha_{\mu}\alpha_{\nu}\big{)};\quad\partial^{\alpha}\alpha_{\alpha}=\frac{1}{\sqrt{\beta^{2}}}\big{(}2w^{2}-\alpha^{2}\big{)};\quad
u_{\alpha}\partial^{\alpha}\alpha^{2}=0;$
$\displaystyle\partial_{\mu}w_{\nu}=\frac{1}{\sqrt{\beta^{2}}}\big{(}\alpha_{\mu}w_{\nu}-\frac{1}{2}\epsilon_{\nu\rho\sigma\lambda}\varpi^{\rho\sigma}\varpi^{\lambda}_{\,\,\mu}\big{)};\quad\partial^{\alpha}w_{\alpha}=-3\frac{w\cdot\alpha}{\sqrt{\beta^{2}}};\quad
u_{\alpha}\partial^{\alpha}w^{2}=0;$
$\displaystyle\alpha^{\sigma}\partial_{\mu}\alpha_{\sigma}=w^{\sigma}\partial_{\mu}w_{\sigma}=\frac{1}{\sqrt{\beta^{2}}}\big{(}w^{2}\alpha_{\mu}-(\alpha\cdot
w)w_{\mu}\big{)};\quad\partial_{\mu}(\alpha\cdot w)=0;$
$\displaystyle\partial_{\alpha}\gamma^{\alpha}=0;\quad\partial^{\alpha}\Delta_{\alpha\beta}=-\frac{\alpha_{\beta}}{\sqrt{\beta^{2}}}.$

## References

  * (1) L. Adamczyk, et al., Nature 548, 62 (2017).  DOI 10.1038/nature23004
  * (2) D.E. Kharzeev, J. Liao, S.A. Voloshin, G. Wang, Prog. Part. Nucl. Phys. 88, 1 (2016).  DOI 10.1016/j.ppnp.2016.01.001
  * (3) D.E. Kharzeev, L.D. McLerran, H.J. Warringa, Nucl. Phys. A 803, 227 (2008).  DOI 10.1016/j.nuclphysa.2008.02.298
  * (4) K. Fukushima, D.E. Kharzeev, H.J. Warringa, Phys. Rev. D78, 074033 (2008).  DOI 10.1103/PhysRevD.78.074033
  * (5) U. Heinz, R. Snellings, Ann. Rev. Nucl. Part. Sci. 63, 123 (2013).  DOI 10.1146/annurev-nucl-102212-170540
  * (6) R. Hakim, _Introduction to Relativistic Statistical Mechanics: Classical and Quantum_ (World Scientific Publishing Co Inc, 2011)
  * (7) W. Israel, Gen. Rel. Grav. 9, 451 (1978).  DOI 10.1007/BF00759845
  * (8) V. Canuto, H.Y. Chiu, Phys. Rev. 173, 1220 (1968).  DOI 10.1103/PhysRev.173.1220
  * (9) X.G. Huang, A. Sedrakian, D.H. Rischke, Annals Phys. 326, 3075 (2011).  DOI 10.1016/j.aop.2011.08.001
  * (10) P. Kovtun, JHEP 07, 028 (2016).  DOI 10.1007/JHEP07(2016)028
  * (11) N. Weickgenannt, X.L. Sheng, E. Speranza, Q. Wang, D.H. Rischke, Phys. Rev. D100, 056018 (2019).  DOI 10.1103/PhysRevD.100.056018
  * (12) J.H. Gao, Z.T. Liang, Phys. Rev. D100, 056021 (2019).  DOI 10.1103/PhysRevD.100.056021
  * (13) J.W. Chen, S. Pu, Q. Wang, X.N. Wang, Phys. Rev. Lett. 110(26), 262301 (2013).  DOI 10.1103/PhysRevLett.110.262301
  * (14) D.N. Zubarev, A.V. Prozorkevich, S.A. Smolyanskii, Theoretical and Mathematical Physics 40(3), 821 (1979).  DOI 10.1007/BF01032069
  * (15) C.G. van Weert, Annals of Physics 140, 133 (1982).  DOI 10.1016/0003-4916(82)90338-4
  * (16) F. Becattini, L. Bucciantini, E. Grossi, L. Tinti, Eur. Phys. J. C75(5), 191 (2015).  DOI 10.1140/epjc/s10052-015-3384-y
  * (17) T. Hayata, Y. Hidaka, T. Noumi, M. Hongo, Phys. Rev. D92(6), 065008 (2015).  DOI 10.1103/PhysRevD.92.065008
  * (18) M. Hongo, Annals Phys. 383, 1 (2017).  DOI 10.1016/j.aop.2017.04.004
  * (19) F. Becattini, M. Buzzegoli, E. Grossi, Particles 2(2), 197 (2019).  DOI 10.3390/particles2020014
  * (20) F. Becattini, W. Florkowski, E. Speranza, Phys. Lett. B789, 419 (2019).  DOI 10.1016/j.physletb.2018.12.016
  * (21) F. Becattini, Phys. Rev. Lett. 108, 244502 (2012).  DOI 10.1103/PhysRevLett.108.244502
  * (22) F. Becattini, E. Grossi, Phys. Rev. D92, 045037 (2015).  DOI 10.1103/PhysRevD.92.045037
  * (23) M. Buzzegoli, F. Becattini, JHEP 12, 002 (2018).  DOI 10.1007/JHEP12(2018)002
  * (24) K. Jensen, R. Loganayagam, A. Yarom, JHEP 05, 134 (2014).  DOI 10.1007/JHEP05(2014)134
  * (25) V.E. Ambruş, E. Winstanley, Phys. Lett. B734, 296 (2014).  DOI 10.1016/j.physletb.2014.05.031
  * (26) V.E. Ambrus, E. Winstanley, arXiv: 1908.10244 (2019)
  * (27) F. Becattini, M. Buzzegoli, A. Palermo, arXiv:2007.08249 (2020)
  * (28) F. Becattini, Phys. Rev. D97(8), 085013 (2018).  DOI 10.1103/PhysRevD.97.085013
  * (29) F. Becattini, D. Rindori, Phys. Rev. D99(12), 125011 (2019).  DOI 10.1103/PhysRevD.99.125011
  * (30) G.Y. Prokhorov, O.V. Teryaev, V.I. Zakharov, Particles 3(1), 1 (2020).  DOI 10.3390/particles3010001
  * (31) G.Y. Prokhorov, O.V. Teryaev, V.I. Zakharov, Phys. Rev. D99(7), 071901 (2019).  DOI 10.1103/PhysRevD.99.071901
  * (32) G.Y. Prokhorov, O.V. Teryaev, V.I. Zakharov, JHEP 03, 137 (2020).  DOI 10.1007/JHEP03(2020)137
  * (33) M. Buzzegoli, E. Grossi, F. Becattini, JHEP 10, 091 (2017).  DOI 10.1007/JHEP07(2018)119,10.1007/JHEP10(2017)091.  [Erratum: JHEP07,119(2018)]
  * (34) M. Hongo, Y. Hidaka, Particles 2(2), 261 (2019).  DOI 10.3390/particles2020018
  * (35) N. Landsman, C. van Weert, Physics Reports 145(3), 141 (1987).  DOI http://dx.doi.org/10.1016/0370-1573(87)90121-9
  * (36) A. Flachi, K. Fukushima, Phys. Rev. D98(9), 096011 (2018).  DOI 10.1103/PhysRevD.98.096011
  * (37) G. Prokhorov, O. Teryaev, Phys. Rev. D97(7), 076013 (2018).  DOI 10.1103/PhysRevD.97.076013
  * (38) G. Prokhorov, O. Teryaev, V. Zakharov, Phys. Rev. D98(7), 071901 (2018).  DOI 10.1103/PhysRevD.98.071901
  * (39) L.J. Tassie, H.A. Buchdahl, Austral. J. Phys. 17, 431 (1964).  DOI 10.1071/PH640431
  * (40) H.A. Buchdahl, L.J. Tassie, Austral. J. Phys. 18, 109 (1965).  DOI 10.1071/PH650109
  * (41) H. Bacry, P. Combe, J.L. Richard, Nuovo Cim. A67, 267 (1970).  DOI 10.1007/BF02725178
  * (42) V.I. Ritus, Annals Phys. 69, 555 (1972).  DOI 10.1016/0003-4916(72)90191-1
  * (43) C.N. Leung, S.Y. Wang, Nucl. Phys. B747, 266 (2006).  DOI 10.1016/j.nuclphysb.2006.04.028
  * (44) M. Laine, A. Vuorinen, Lect. Notes Phys. 925, pp.1 (2016).  DOI 10.1007/978-3-319-31933-9
  * (45) K. Fukushima, D.E. Kharzeev, H.J. Warringa, Nucl. Phys. A836, 311 (2010).  DOI 10.1016/j.nuclphysa.2010.02.003
  * (46) A. Vilenkin, Phys. Rev. D22, 3080 (1980).  DOI 10.1103/PhysRevD.22.3080
  * (47) M.A. Metlitski, A.R. Zhitnitsky, Phys. Rev. D72, 045011 (2005).  DOI 10.1103/PhysRevD.72.045011
  * (48) A. Sadofyev, V. Shevchenko, V. Zakharov, Phys. Rev. D 83, 105025 (2011).  DOI 10.1103/PhysRevD.83.105025
  * (49) D.E. Kharzeev, Prog. Part. Nucl. Phys. 75, 133 (2014).  DOI 10.1016/j.ppnp.2014.01.002
  * (50) V.I. Zakharov, Lect. Notes Phys. 871, 295 (2013).  DOI 10.1007/978-3-642-37305-3˙11
  * (51) B. Feng, D.F. Hou, H.C. Ren, Phys. Rev. D99(3), 036010 (2019).  DOI 10.1103/PhysRevD.99.036010
  * (52) K. Landsteiner, E. Megias, F. Pena-Benitez, Lect. Notes Phys. 871, 433 (2013).  DOI 10.1007/978-3-642-37305-3˙17
  * (53) E.V. Gorbar, V.A. Miransky, I.A. Shovkovy, X. Wang, Phys. Rev. D88(2), 025025 (2013).  DOI 10.1103/PhysRevD.88.025025