Axial chiral vortical effect in a sphere with finite size effect

Relativistic heavy ion collisions (RHICs) are utilized to produce quark-gluon plasmas (QGPs) at high temperature and nonzero baryon density. A typical (off-central) collision exposes the QGP thus generated under an ultra-strong magnetic field and endows it with a high angular momentum. A number of novel transport phenomena [1–15]have been proposed, including the axial-chiral-vortical-effect (ACVE). The ACVE refers to the axial vector current, i.e., the spin density of fermions in response to the global angular momentum, and it is expected to be detected via the polarization of Lambda post hadronization. The ACVE is also expected inside the core of a fast spinning neutron star [16–18]. In this work, we shall focus on the theoretical aspect of the ACVE.

In a thermal equilibrium ensemble, the ACVE is represented in terms of the global angular velocity ω by the formula

where the coefficient σ is referred to as the axial vortical conductivity and the ellipsis represents higher power in ω. Based on a pioneer work by Son and Surowka [7, 19] and a supplemental work by Neiman and Oz [8], the axial vortical conductivity is restricted by thermodynamic laws to the following general form in the chiral limit,

where $ \mu_V $ and $ \mu_A $ are the vector and axial vector chemical potentials, respectively, and the coefficient c in front of the temperature square has to be determined by other means.

Besides the hydrodynamic approach, Eq. (2) with $ c=1/6 $ was first derived by Vilenkin via the solution of a free Dirac equation in a rotating cylinder [20–22], and the axial vortical conductivity for non-interacting fermions reads

The same expression was obtained by Landsteiner et al. via the Kubo formula to one-loop order [23]. There is also a large body of literature on the derivation of Eq. (3) from kinetic theory [24, 25] or holography [26, 27]. Beyond Eq. (3), the authors of [10, 14] discovered higher order corrections to the coefficient c in QED or QCD coupling. The authors of Ref. [28] determined the higher order terms in ω, i.e., the ellipsis in Eq. (1) for massless fermions, and ended up with a closed form of the axial-vector current, and the authors of Ref. [9] derived the leading order correction of the fermion mass. A recent calculation [29] of axial current for massless fermions in a general thermodynamic equilibrium with rotation and acceleration (within a formalism “far from the boundary,” that is, without enforcing boundary conditions) reproduces the known results for rotating equilibrium, such as those in Ref. [28], but it extends them to systems including acceleration.

In this study, we explore the axial vortical effect in a finite sphere of radius R subject to the MIT boundary condition. Our motivation is twofold. First, a system rotating with constant angular velocity has to be finite in the direction transverse to the rotation axis, as restricted by the subluminal linear speed on the boundary. Second, a finite sphere serves as a better approximation to the shape of the QGP fireball in heavy ion collisions and the quark matter core of a neutron star than the infinitely long cylinder considered in literature. The MIT boundary condition effectively separates the deconfinement phase of the interior and the confinement phase outside. However, we could neither include the strong coupling underlying the near-perfect fluidity inside an actual QGP fireball nor describe its rapid expansion, especially in the early stage of its evolution. Far from the boundary, where the finite size effect can be ignored, we reproduce in spherical coordinates exactly the same form of the axial-vector current in the chiral limit derived in cylindrical coordinates [28]. We also carry out the fermion mass correction to all orders with the leading order matching the result in Ref. [9], which was derived with the Kubo formulation. The infinite series in powers of the mass correction indicates that the leading order correction for the mass of an s quark at the RIHC temperature is quite accurate. More importantly, we obtain an analytic approximation of the axial vector current on the spherical boundary with the aid of the asymptotic formula of the Bessel function of large argument and large order. For $ {\boldsymbol{\omega}}=\omega\hat{{\boldsymbol{z}}} $, we find that

where $ {\boldsymbol{e}}_\rho $ is the unit radial vector of the cylindrical coordinate systems $ (\rho, \phi, z) $. For $ T\gg1/R $ and the fermion mass $ M\ll T $, the axial vortical conductivity parallel to $ {\boldsymbol{\omega}} $ is

and that perpendicular to $ {\boldsymbol{\omega}} $ is

where θ is the polar angle with respect to the direction of the angular velocity. Note that we have to set $ \mu_A=0 $ and $ \mu_V=\mu $ because the MIT boundary condition breaks the chiral symmetry even for massless fermions. To the best of our knowledge, the perpendicular component has never been reported in literature, and its existence may shed some light on the longitudinal (with respect to the beam direction) polarization in heavy ion collisions.

The organization of the paper is as follows. In Sec. II, general properties of the axial vortical effect are discussed from symmetry perspectives. In Sec. III, we lay out the general formulation of the chiral magnetic effect in spherical coordinates with the MIT boundary condition. The axial vortical effects of massless and massive fermions are calculated in Sec. IV and Sec. V, respectively. Sec. VI concludes the paper with a qualitative speculation on the impact of the finite size effect for heavy ion collisions. Some technical details are deferred to Appendices. We also include two additional Appendices for self-containment, one for an alternative derivation of the closed end formula of the axial-current in cylindrical coordinates and the other one for the mass correction via the Kubo formula under dimensional regularization. Throughout the paper, we shall stay with the notation of Eqs. (5) and (6) by setting $ \mu_A=0 $ and $ \mu\equiv\mu_V $. Furthermore, the size of the sphere is assumed to be sufficiently large in comparison with the length scale corresponding to the temperature or chemical potential for the boundary condition to be analytically soluble.

In this section, we explore the axial vortical effect from symmetry perspectives. The validity of the conclusion reached here is not limited to a free Dirac as considered in literature and the subsequent sections of this work.

The axial vortical effect refers to the thermal average of the spatial component of the axial vector current density $ {\boldsymbol{\mathcal{J}}}_A $ in the presence of a nonzero angular momentum. Taking the direction of the angular momentum as z-axis, we have

In terms of the field theoretic Hamiltonian $ \mathcal{H} $, conserved charge $ \mathcal{Q} $, and z-component of the angular momentum $ \mathcal{J}_{z} $, the density matrix at thermal equilibrium reads

where T is the temperature, μ is the chemical potential, ω is the angular velocity, and Z is the normalization constant such that $ {\rm Tr}\varrho=1 $.

Introducing the basic vector of cylindrical coordinates $ \hat{{\boldsymbol{z}}} $ and

the ensemble average (7) can be decomposed into its longitudinal component

and its transverse components

with

where the dependence on the cylindrical coordinates, chemical potential, and angular velocity is explicitly indicated and will be suppressed in subsequent sections. We have

and consequently

Assuming that the Hamiltonian and the boundary condition are invariant under spatial rotation, spatial inversion, time reversal, and charge conjugation, we have

where $ \mathcal{R}(\alpha) $ is a Hilbert space operator of a rotation about the z-axis by an angle α, and $ \mathcal{P} $, $ \mathcal{T} $, and $ \mathcal{C} $are Hilbert space operators for spatial inversion, time reversal, and charge conjugation, respectively. Together with the transformation laws of the axial vector current $ {\boldsymbol{\mathcal{J}}}_{A}({\boldsymbol{r}}) $

it follows that

where $ \overleftrightarrow{D}(\alpha) $ is the dyadic notation of the $ 3\times 3 $ rotation matrix

Because of the relations

and $ {\boldsymbol{e}}_s(\phi+\pi)=-{\boldsymbol{e}}_s(\phi) $s, Eqs. (23) and (24) imply that the longitudinal and transverse components of the axial current defined in Eq. (10) and (11) are independent of the azimuthal angle as expected, and the transverse component is odd in z, i.e.,

Consequently, there cannot be a transverse axial vortical effect for an infinitely long cylinder because the axial vector current is independent of z. This, however, is not the case with a sphere as the z dependence cannot be ignored. The oddness with respect to z implies only zero transverse axial vector current on the equatorial plane of the sphere. Indeed, the subsequent sections show that the transverse component of the axial-vector current does exist on the spherical boundary for a free Dirac field and does vanish on the equatorial plane. The other two equations, Eqs. (25) and (26), imply that the thermal average of the axial-vector current is always odd with respect to the angular velocity and even with the chemical potential.

Before concluding this section, we remark that some of the relations above can be readily generalized to a non-equilibrium density matrix with its time development dictated by the Liouville theorem. For instance^①, for a homogeneous and expanding system, as long as relations (15), (16), and (18) hold initially, they will hold always. Then, relation (29) and its implications discussed above remain valid always.

The Hamiltonian for a Dirac fermion in a uniformly rotating system with angular velocity $ {\boldsymbol{\omega}}=\omega{\boldsymbol{e}}_{z} $ can be written as [28, 31]

where $H_{0}=-{\rm i}{\boldsymbol{\alpha}}\cdot{\boldsymbol{\nabla}}+\beta M$ is the free Hamiltonian, $J_{z}= \dfrac{1}{2}\Sigma_{3}-{\rm i}(x\partial_{y}-y\partial_{x})$ is the z-component of the total angular momentum, μ is the the chemical potential of the system, M is the mass of the Dirac fermion, and $ {\boldsymbol{\alpha}}=\gamma^{0}{\boldsymbol{\gamma}} $, $ \beta=\gamma^{0} $. We work in the Dirac representation for gamma matrices $ \gamma^{\mu} $ as follows:

The last two terms of Eq. (30) are included in the single particle Dirac Hamiltonian because it is Eq. (30), when being sandwiched between Dirac fields $ \psi, \psi^\dagger $

to define the density operator ${\rm e}^{(-\frac{\mathcal{H}}{T})}$ for the thermal average.

In this section, we consider the eigenfunctions of the Hamiltonian in spherical coordinates. The eigenfunctions of the Hamiltonian satisfy

where E is the eigen-energy of $ H_{0} $. The solutions of Eq. (33) can be chosen as the common eigenfunctions of these four commutative Hermitian operators: Hamiltonian H, square of total angular momentum $ {\boldsymbol{J}}^{2} $, z-component of total angular momentum $ J_{z} $, and parity operator P. We list the eigenfunctions in spherical coordinates as follows,

where $j,~ l,~ m$ denote the eigenvalues $ j(j+1), (-1)^l, m $ of $ {\boldsymbol{J}}^{2}, P, J_{z} $, respectively; the spinor spherical harmonics $ Z_{j,j\pm\frac{1}{2},m}(\theta,\phi) $ are defined as

and $ f(r) $, $ g(r) $ are the radial wave functions. Making use of the following relations,

one can obtain the following differential equation satisfied by $ f(r) $ in Eq. (34),

which is the l-th order spherical Bessel equation. The radial function $ g(r) $ in Eq. (34) can be expressed by $ f(r) $,

We list the solutions of $ f(r) $ and $ g(r) $ in Table 1, with $ k>0 $, $ E_{k}=\sqrt{k^{2}+M^{2}} $, and C the normalization factor.

$ l=j\pm\dfrac{1}{2} $	$ E=E_{k} $	$ E=-E_{k} $
$ f(r) $	$ Cj_{l}(kr) $	$ -C\sqrt{\dfrac{E_{k}-M}{E_{k}+M}}j_{l}(kr) $
$ g(r) $	$ C\sqrt{\dfrac{E_{k}-M}{E_{k}+M}}j_{l\mp1}(kr) $	$ Cj_{l\mp1}(kr) $

Table 1.  Solutions of $ f(r) $ and $ g(r) $.

For a spherical volume of radius R, the quantization of the radial momentum k depends on the boundary condition. An approximate boundary condition for a fireball of QGP is derived from the MIT bag model [30] and reads

which requires that the solution of the Dirac equation on the boundary implements the eigenfunction of $-{\rm i}{\boldsymbol{\gamma}}\cdot{\boldsymbol{\hat r}}$ of eigenvalue one. As $ \gamma_5{\boldsymbol{\gamma}}\cdot{\boldsymbol{\hat r}}=-{\boldsymbol{\gamma}}\cdot{\boldsymbol{\hat r}}\gamma_5 $, the MIT boundary condition breaks the chiral symmetry even for massless fermions. In accordance with Eqs. (34) and (36), the radial wave function satisfies

for $ l=j\mp 1/2 $. For the solutions of the free Dirac equation in Table 1^②, the MIT boundary condition reads

for the positive energy state of $ l=j-1/2 $ and

for the positive energy state of $ l=j+1/2 $, where

The boundary conditions for the negative energy states are based on the charge conjugation, i.e.,

Employing the integration formula

and the formulas of the derivative $ J_\nu^\prime(z) $ in terms of $ J_\nu(z) $ and $ J_{\nu\pm1}(z) $, the normalization constant in Table 1 is readily determined

The boundary conditions Eqs. (41) and (42) can be solved approximately for $ kR\gg1 $ and $ j\ll kR $ with the aid of the asymptotic formula of the spherical Bessel function

We find

for $ l=j\pm1/2 $. The summation of k can be converted to an integral

and the normalization constant in Eq. (46) under both conditions of Eqs. (41) and (42) is simplified to

The quantized Dirac field can be expressed by the eigenfunctions of the Hamiltonian H as follows:

where $ a_{kjlm}^{\dagger} $ and $ a_{kjlm} $ are the creation and annihilation operators of particles, whereas $ b_{kjlm}^{\dagger} $ and $ b_{kjlm} $ are those of anti-particles. The explicit forms of $ u_{kjlm}({\boldsymbol{r}}) $ and $ v_{kjlm}({\boldsymbol{r}}) $ are

We have

where $ E_{k}=\sqrt{k^{2}+M^{2}} $. The ensemble average (7) of $ a_{kjlm}^{\dagger}a_{kjlm} $ and $ b_{kjlm}^{\dagger}b_{kjlm} $ with the density operator (8) gives rise to the Fermi-Dirac distribution functions,

and the thermal expectation values of $ a_{kjlm}^{\dagger}b_{kjlm}^{\dagger} $, $ b_{kjlm}^{\dagger}a_{kjlm}^{\dagger} $, $ a_{kjlm}b_{kjlm} $, and $ b_{kjlm}a_{kjlm} $ are all zero.

In the following, we calculate the axial vector current of the uniformly rotating system of Dirac fermions. The axial vector current is the ensemble average of the corresponding operator, i.e.,

where $ {\boldsymbol{J}}_{A}^{\mathrm{vac}}=-\displaystyle\sum_{kjlm}v_{kjlm}^{\dagger}{\bf{\Sigma}}v_{kjlm}=0 $ is the vacuum term, and the charge conjugation relation in Eq. (39) has been employed in the last step. It follows from the relation

that

and Eq. (55) becomes

Introducing the following ϕ independent functions

and

we have

and

with $ |C_\pm|^2 $ given by the upper (lower) line of Eq. (46) In particular, the expression of $ \eta_{jlm}(\theta) $ can be reduced to

and

with

It follows from the property $Y_{lm'}(\pi-\theta,\phi+\pi)= $ $ (-)^lY_{lm'}(\theta,\phi)$ that $ \eta_{jlm}(\pi-\theta)=-\eta_{jlm}(\theta) $ and thereby $ \eta_{jlm}\left(\dfrac{\pi}{2}\right)=0 $.

Before concluding this section, we point out an interesting property of the MIT boundary condition, which is not dictated by symmetries: the axial vector current vanishes along the equator of the fireball. Indeed, Eq. (40) implies that

with

Writing $ \Theta_{jm}(\theta) $ in terms of associated Legendre functions $ P_l^\mu(\cos\theta) $ and using the explicit form of $ P_l^\mu(0) $, we find that

See Appendix A for details of the proof.

For massless fermions, $ M=0 $ and $ E_{k}=k $ in Table 1. Far from the boundary, the main support of the axial vector current comes from the spherical Bessel function with $ l=O(1) $. Together with the condition $ T\gg 1/R $ and $ k\sim T $, we have $ kR\gg 1 $ for typical radial momentum and the approximation in the last paragraph of Sec. III.A becomes useful. Using the relations

the z-component of Eq. (55) reads

where we have turned the summation over k to integral according to Eq. (49) and extended the integration domain to $ (-\infty,\infty) $ via Eq. (69).

The Taylor expansion of the axial current in Eq. (70) reads

where the coefficient

with $f(x)=1/({\rm e}^x+1)$ the Fermi-Dirac distribution function and $ f^{(n)}(x) $ its n-th derivative. In the second step of Eq. (72), we have substituted the explicit form of the wave function in Eq. (51) together with Eq. (35) for the spinor spherical harmonics. Applying the addition formula

for $ r^\prime=r $, $ \theta^\prime=\theta $, and $ \phi^\prime=\phi+\epsilon $, we find

where $ \xi\equiv 2r\sin\theta\sin\dfrac{\epsilon}{2} $. After $ 2n $ times of integration by part with respect to k, we obtain that

Only the $ (2n+1) $-th power of ξ inside the curly brackets contributes to $ C_n $. Together with the integrals

we have

with $ \rho=r\sin\theta $. Substituting into Eq. (71) and summing up the series, we end up with

which is in agreement with the closed-end formula derived in cylindrical coordinates in Ref. [28]. An alternative derivation in cylindrical coordinates is presented in Appendix B. To the cubic order in ω, Eq. (78) yields the formula derived in Ref. [28]. As is shown in the step from Eq. (74) to Eq. (75), the key reason for having the closed form of the axial current Eq. (78) is that the density of states for massless fermions is proportional to an integer power of the energy $ E_k=k $ so that the integration by part terminates with a finite number of terms for arbitrary n. This is no longer the case for massive fermions.

Equation (78) is plotted in Fig. 1, where we set $ \omega=0.01T $, a rough estimate of the vorticity of the QGP fireball created in RHICs. The pole at $ \omega r\sin\theta=1 $ occurs where the linear speed of rotation reaches the speed of light and the linear speed beyond the pole becomes superluminal, which is not admissible. Therefore, the Hamiltonian in Eq. (30) applies only to a finite volume, which in the case of the sphere discussed in this section requires its radius below $ 1/\omega $. The axial vector current in Eq. (78) is thereby free from the pole within the sphere, but the finite size effect becomes significant. Unless the finite size effect falls to zero faster than a power series in $ r/R>\omega r $, its contribution will be of the same order of importance as that of the higher order terms of Eq. (78).

Figure 1.  (color online) Ratio of axial vector current $ J_A^z $ over $ T^3 $ of massless fermions far from the boundary in Eq. (78) with the angular velocity $ \omega=0.01\,T $ as a function of $ \rho T $, where ρ and T are the radius and temperature, respectively. The black, blue, and red lines indicate $ \mu/T $ = 0.5, 1.0, and 1.5, respectively. The inner panel is for $ \rho T $= 0–1.0.

Regarding the transverse component far from the boundary, the typical contribution to the thermal average comes from $ j\ll kR $, and the sum over k in Eq. (55) and wave function normalization can be approximated by Eqs. (49) and (50), respectively. Following Eqs. (60), (63), (64), and (65), we obtain

where

The absence of the transverse components is expected because the finite size effect can be neglected in the bulk, and the spherical and cylindrical shapes of the volume make no difference there.

Regarding the boundary of a QGP fireball, we have to distinguish the radial momentum k of the wave function $ u_{kjlm} $ for $ l=j-1/2 $ and $ l=j+1/2 $ because of the different quantization conditions in Eqs. (41) and (42). Based on Table 1 and Eq. (50), for $ M=0 $,

The radial momentum k of Eq. (81) follows from Eqs. (41) and (42). The axial vector current on the boundary is obtained upon substitution of Eq. (81) into Eq. (55). An analytic expression of the boundary axial vector current can be derived for the linear order term of the Taylor expansion in ω, i.e., the chiral conductivity, at high temperature, i.e. $ T\gg1/R $. We have

where $ k_{nj}^\mp $ stands for the solutions of Eq. (41) ("$ - $" sign) or Eq. (42) ("$ + $" sign). According to the definition in Eq. (67) and the explicit form of the spinor spherical harmonics in Eq. (35),

where

and the addition formula of the spherical harmonics has been employed. Combining Eqs. (82), (83), and (84), we arrive at

To evaluate the summation over k and j under the condition $ T\gg 1/R $ or $ \mu\gg 1/R $, we notice that $ kR\gg 1 $and the wave functions of the large j become important because of the centrifugal force. The asymptotic formula Eq. (47) is no longer sufficient to serve the purpose, and one has to switch to the Debye formula [32] for the Bessel function of a large argument and large order^③,

which implies

for a spherical Bessel function. Then, the MIT boundary conditions in Eqs. (41) and (42) become

with

The large j serves as the guideline to sort the order of approximation. Eq. (89) gives rise to the leading order relation between β and $ \beta^\prime $

Substituting Eq. (90) to the RHS of Eq. (88) and dropping the terms beyond the order of $ 1/j $, the boundary condition is reduced to

with the solutions

for the upper sign and

for the lower sign, where n is a positive integer. Together with the relation between β and the radial momentum k in Eq. (89), we have [33]

to the leading order of a large j for both signs in Eq. (88). Converting the summation over k and j in Eq. (85) to integrals, we obtain the leading order axial vector current on the boundary

Therefore, the longitudinal axial vortical conductivity vanishes along the equator, which is consistent with the general statement according to Eq. (68) and matches the axial vortical conductivity far from the boundary at the poles ($ \theta=0,\pi $).

To the linear order in ω, the transverse component of the axial vector current $ {\boldsymbol{J}}_A $ is obtained by replacing $ \Theta_{jm}(\theta) $ of the formula for the longitudinal component Eq. (82) by

i.e.

The summation over m can be carried out in a similar manner to that in Eq. (84), and we find, with $ m=m'+\dfrac{1}{2} $ and $ l=j\pm\dfrac{1}{2} $, that

where the derivative formula of the Legendre polynomial

is employed. Approximating the sum over n and j by integrals of Eqs. (97) and (98), we obtain the transverse component of the axial vector current

which is of the same order of magnitude as the longitudinal component. Restoring the cylindrical coordinates via

we have

which is independent of the azimuthal angle and odd in z and is consistent with the symmetry argument in Sec. II.

For massive fermions, the same approximation of the MIT boundary condition applied to massless fermions reduces the axial vector current $ {\boldsymbol{J}}_{A} $ in Eq. (55) far from the boundary to

with $ E_k=\sqrt{k^2+M^2} $. As ${\rm d}kk^2={\rm d}E_kE_k\sqrt{E_k^2-M^2}$, the density of states is no longer an integer power of the energy $ E_k $, and a closed-end formula such as Eq. (78) does not exist. We shall stay with the linear response of $ J_A^z $ to ω in what follows and calculate the axial vertical conductivity. It is straightforward to verify that the combination

with the radial wave functions in Table 1 and the normalization constant Eq. (50) at a given k, is independent of the mass M and thereby takes the same massless form. For the longitudinal component, the spinor spherical harmonics part can be reduced the same way as in Sec. IV.A, and Eq. (103) becomes, to the order ω,

Using the relation of Eq. (73), Eq. (105) becomes

where we have transformed the integration variable from k to $ x=E/T $ with $ E=\sqrt{k^2+M^2} $. The integral Eq. (106) can be converted to a contour integral by the observation that

where the first two terms of the Taylor expansion of $ \sqrt{z^2-\dfrac{M^2}{T^2}} $ in the powers of M are included in $ I^\prime $, i.e.,

with $ a=\dfrac{M}{T} $. Then, the integrand of I vanishes sufficiently fast at infinity so that the integration path can be closed from infinity on the upper or lower z-plane and the integral equals the sum of residues at the poles of the distribution function within the contour. Closing the path from the upper plane, we have the poles

within the contour, i.e., $ n=0,1,2,... $. Consequently,

Combining Eqs. (106), (107), and (108), we have

with the axial vertical conductivity of massive fermions

The binomial expansions of the square roots in Eq. (112) enable us to write

where $ \zeta(...) $ denotes the Hurwitz zeta function defined by

Away from the branch points of the square roots in the summands, the infinite series Eq. (110) converges uniformly with respect to a, and thereby, the radius of convergence of the power series Eq. (113) corresponds to the absolute value of the closest branch point to the origin of the complex a-plane, i.e.,$ \sqrt{\pi^2+(\mu/T)^2} $. This can also be inferred from the asymptotic behavior of the expansion coefficients of Eq. (113). We also obtain Eq. (113) in Appendix B by using a cylindrical coordinate system and in Appendix C by the Kubo formula via a thermal diagram, which shows that this result, derived by different methods, is robust. In particular, the thermal diagram requires UV regularization, but the result is independent of regularization schemes.

At zero temperature, the summation over n in Eq. (112) can be converted to an integral, and we obtain that

The zero σ for $ \mu<M $ is obvious from Eq. (105), where the derivative of the distribution function vanishes exponentially in the limit $ T\to 0 $ for all k. The case with $ \mu>M $ returns the massless result derived in Sec. IV.A for $ M=0 $.

The axial vector current with mass correction is plotted in Fig. 2, where the solid line is $ J_{A}^{z(n)} $, and the dashed line is $ J_{A}^{z(1)} $. Here, $ J_{A}^{z(0)} $ is the axial vector current at $ M=0 $, $ J_{A}^{z(1)} $ is the axial vector current with only $ M^{2} $ correction (the result of $ J_{A}^{z(1)} $ is also obtained in Ref. [9]), and $ J_{A}^{z(n)} $ is the result including mass correction up to $ M^{2n} $. Their concrete expressions are

Figure 2.  (color online) Ratio of axial vector current $ J_{A}^{z(n)} $ including mass correction up to $ M^{2n} $ over massless current $ J_{A}^{z(0)} $ ($ J_{A}^{z(n)}/J_{A}^{z(0)} $) as a function of the product of M over T ($ M/T $), where the black, blue, green, and red lines indicate $ \mu/T $=0, 0.5, 1.0, and 1.5, respectively. The dashed lines are $ J_{A}^{z(1)}/J_{A}^{z(0)} $, and the solid lines are $ J_{A}^{z(n)}/J_{A}^{z(0)} $.

We can see clearly that $ J_{A}^{z(n)}/J_{A}^{z(0)} $ decreases with $ M/T $. This is because the presence of mass generally inhibits the fluidity, thus suppressing the vortical conductivity. While the presence of chemical potential slows down this inhibition, when we fix $ M/T $, $ J_{A}^{z(n)}/J_{A}^{z(0)} $ and $ J_{A}^{z(1)}/J_{A}^{z(0)} $ increase with increasing $ \mu/T $.

An s quark is taken as an example. We set $ M=150 $ MeV, $ \mu/T=1.0 $, and $ n=2000 $ and list the numerical values of the mass correction in Table 2.

$T/\mathrm{MeV}$	$ J_{A}^{z(1)} / J_{A}^{z(0)} $	$ J_{A}^{z(n)} / J_{A}^{z(0)} $	$ \left(T^2\omega\displaystyle\sum\limits_{r=2}^{n}A_{r}M^{2r} \right) /J_{A}^{z(0)} $
$ 100 $	$ 0.737754 $	$ 0.741961 $	$ 4.20713\times10^{-3} $
$ 150 $	$ 0.883446 $	$ 0.884285 $	$ 8.38337\times10^{-4} $
$ 200 $	$ 0.934439 $	$ 0.934704 $	$ 2.65816\times10^{-4} $
$ 250 $	$ 0.958041 $	$ 0.95815 $	$ 1.08964\times10^{-4} $

Table 2.  Mass correction of axial current when $ M=150 $ MeV, $ \mu/T=1.0 $, $ n=2000 $.

Far from the boundary, the mass correction for the s quark is modest for the selected temperature and chemical potential and is dominated by the leading order $ O(M^2) $ correction. On the boundary, the leading order mass correction is $ O(M) $, as shown below. The mass suppression for the s quark is thereby much stronger there.

For the transverse component of the axial vector current of massive fermions far from the boundary, all we need is to replace k in $ g_m(k) $ of Eq. (79) with $ E_k $, and the result remains zero, the same as in the massless case.

An analytical result can also be obtained for the leading order mass correction on the spherical boundary under the same approximation of Sec. IV.B, i.e., $ T\gg 1/R $. For massive fermions, it follows from Eq. (46) that Eq. (81) is replaced by

with $ b=E_k/\sqrt{E_{k}^2-M^2} $, where we have substituted Eq. (43) for the trigonometric functions in the normalization constant Eq. (46) and made the approximation $ 2j+2\simeq 2j $in the last term inside the parentheses for a large j. The conversion from the sum of the radial momentum into an integral proceeds in the same way as for the massless case in Sec. IV.B, and we obtain the following form of the axial vector current to the order$ O(\omega) $

The integration over u can be carried out readily

Consequently, the leading order mass correction is $ O(M) $, stronger than $ O(M^2) $ for the mass correction far from the boundary. Substituting Eq. (119) into Eq. (118) and setting $ E_k=k $, we find

where the first term, $ J_{A}^{z(0)} $, is the axial-vector current of massless fermions given by Eq. (95), and the leading order mass correction reads

which is an even function of μ. Adding Eqs. (95) and (121), we have the longitudinal axial vector current on the boundary up to the leading order mass correction.

where $ J_{A}^{z(B)} $is the axial vector current with only leading order mass correction on the boundary. We can clearly observe that the mass correction is stronger on the boundary than that far from the boundary. The coefficient of ω of Eq. (116) gives rise to the axial vortical conductivity on the boundary, Eq. (5), presented in the introduction. As$ T\to 0 $,

With the aid of the integral Eq. (119), together with the definition of b, we obtain a closed-end formula of the axial vortical conductivivity to all orders of mass on the boundary

in parallel to Eq. (115) in the bulk.

It is straightforward to extend the above analysis to the transverse component. Starting with Eq. (62) and going through the gymnastics from Eq. (118) to Eq. (116) with $ \cos^2\theta $ replaced by $ \dfrac{1}{8}\sin2\theta $, we find the transverse axial vector current on the boundary up to the leading order of mass correction, i.e.,

At zero temperature, we have

This is valid up to all orders of the mass M.

Our study can be summarized as follows. We started with a general discussion of the axial vortical effect from symmetry perspectives and investigated the axial vortical effect of a free Dirac field in a finite sphere rotating with a given angular velocity ω. For massless fermions far from the boundary, we were able to reproduce the closed-end formula derived within a cylinder in literature. On the boundary, the axial vector current displays both longitudinal and transverse components with respect to the rotation axis, and the magnitude of each component depends on the colatitude angle of the spherical coordinates. For massive fermions, we obtained the mass correction of the chiral conductivity far from and on the boundary. In the former case, we expanded the chiral conductivity to all orders of mass with the leading order correction in agreement with what was reported in literature. In the latter case, we found that the leading order mass correction was stronger than that of the former, $ O(M) $ versus $ O(M^2) $. To the best of our knowledge, the axial vortical effect on the boundary, especially the emergence of the transverse component, has not been explored in literature.

While the values of the above results are mainly theoretical and cannot describe quantitatively the ACVE of a strongly interacting and expanding fireball of QGP, some qualitative speculations on the finite size effect in heavy ion collisions remain instructive. The quadrupole factor $ \cos^2\theta $ in Eq. (5) would suppress the global polarization (z-component of Eq. (4)) and the perpendicular component in Eq. (4) would contribute to the polarization in the reaction plane shown in Fig. 3, e.g., the longitudinal polarization (the polarization along the beam).

Figure 3.  $ x' $ and the beamline define the reaction plane.

To observe the latter effect clearly, we assume that the beam is along $ \hat{{\boldsymbol{y}}} $ and rotate the coordinate system by $ 90^\circ $ around the x-axis, i.e., $ y=-z^\prime=-r\cos\theta^\prime $, $ z=y^\prime= r\sin\theta^\prime\sin\phi^\prime $, and $ x=x^\prime=r\sin\theta^\prime\cos\phi^\prime $, with r being the radial coordinate. In terms of the polar angle, $ \theta^\prime $, and azimuthal angle $ \phi^\prime $ associated with the primed coordinates, the longitudinal component in Eqs. (4) and (5) takes the form

with $b=\left\{\dfrac{\mu^2}{16\pi^2} + \dfrac{1}{48}T^2-\dfrac{M}{32\pi}\left[\mu+2T\ln \left(1+{\rm e}^{-\frac{\mu}{T}}\right)\right]\right\}\omega$. As the fragmenthadrons, e.g., Λ hyperons originating from the boundary layer, are more likely flying in the radial direction, Eq. (127) maps out the longitudinal polarization profile of these hadrons, with$ \phi^\prime $being the angle of the transverse momentum with respect to the reaction plan and $ \theta^\prime $being related to the pseudorapidity via $ \eta=-\ln\tan\dfrac{\theta^\prime}{2} $.

More investigations are required for the finite size effect discovered in this work to be practical with respect to the phenomenology of heavy ion collisions. These include exploring the ACVE with the solution of the Dirac equation in an expanding sphere and/or incorporating the anisotropic ACVE conductivity in Eq. (4) in hydrodynamic models. We hope to report the progress along this line in the near future.

We thank Ren-Da Dong and Xin-Li Sheng for fruitful discussions.

To prove Eq. (68), we substitute the explicit form of $ Z_{j,l,m}^{\dagger}(\theta,\phi) $ into Eq. (67), i.e.,

As $ \Theta_{jm}(\theta) $ is odd in m, we only need to consider the case with $ m>0 $. Setting $ j=l+1/2 $ and $ m=\mu+1/2 $ and using the expression of spherical harmonics in terms of the associated Legendre function, we have

with $ \mu\geqslant 0 $. It follows from the generating function of Legendre polynomials

and the definition

that

Setting $ z=0 $ and comparing the coefficients of $ t^l $ on both sides, we obtain that [34]

It is straightforward to verify that

and

Equation (68) is thereby proved.

In this appendix, we first solve the free Dirac equation in a cylindrical coordinate system and then calculate the axial vector current of the system of massive Dirac fermions, which uniformly rotates with angular velocity $ {\boldsymbol{\omega}}=\omega{\boldsymbol{e}}_{z} $ along the z-axis. We consider only the axial vector current far from the boundary and thereby ignore the finite size effect.

We work in the chiral representation of gamma matrices as adopted in Ref. [35],

with $\sigma_{i} (i=x,~y,~z)$ being the three Pauli matrices. The equation of motion for the free Dirac field $\Psi(t,~{\boldsymbol{r}})$ can be written as

with the Hamiltonian $\hat{H}=-{\rm i}\gamma^{0}{\boldsymbol{\gamma}}\cdot\nabla+\gamma^{0}M$ and the Dirac fermion mass M. Suppose that $ \Psi(t,{\boldsymbol{r}}) $ is an energy eigenstate with eigenvalue E, i.e., $\Psi(t,{\boldsymbol{r}})={\rm e}^{-{\rm i}Et}\psi({\boldsymbol{r}})$; then, Eq. (B2) becomes

which is the energy eigenvalue equation of the Hamiltonian. It can be proved that these four Hermitian operators, $ \hat{H},\hat{p}_{z},\hat{J}_{z},{\bf{\Sigma}}\cdot\hat{{\boldsymbol{p}}} $, are commutative with each other, where $ {\bf{\Sigma}}=\mathrm{diag}\,({\boldsymbol{\sigma}},{\boldsymbol{\sigma}}) $, $ \hat{{\boldsymbol{p}}}=-i\nabla $, $ \hat{{\boldsymbol{J}}}={\boldsymbol{r}}\times\hat{{\boldsymbol{p}}}+\dfrac{1}{2}{\bf{\Sigma}} $, and $ \hat{p}_{z},\hat{J}_{z} $ are the z-components of $ \hat{{\boldsymbol{p}}} $ and $ \hat{{\boldsymbol{J}}} $ , respectively. In the following, we will calculate the common eigenstates of these four operators in a cylindrical coordinate system. We set $ \psi=(\psi_{1},\psi_{2})^{T} $, where $ \psi_{1},\psi_{2} $ are both two-component spinors; then, Eq. (B3) can be replaced by the following two equations,

In a cylindrical coordinate system, the form of $ \nabla^{2} $ is

Now, we solve $ \psi_{1} $ from Eq. (B4). $ \psi_{1} $ can be chosen as

which is the common eigenstate of $ \hat{p}_{z} $ and $-{\rm i}\partial_{\phi}+\dfrac{1}{2}\sigma_{z}$ with eigenvalues $ p_{z} $ and j. Plugging Eq. (B7) into Eq. (B4) gives

which are the Bessel equations of order $ \left(j\mp\dfrac{1}{2}\right) $. The boundary conditions of $ \psi_{1} $ at $ r=0 $ and $ r=\infty $ require that $ E^{2}>M^{2}+p_{z}^{2} $. We can introduce a transverse momentum $ \alpha=\sqrt{E^{2}-M^{2}-p_{z}^{2}} $; then, the eigen-energy becomes $ E= \lambda\sqrt{M^{2}+p_{z}^{2}+\alpha^{2}} $, with $ \lambda=\pm1 $ corresponding to the positive and negative modes. Now, one can obtain $ \psi_{1} $ as

where A is a constant to be determined. As ψ is also an eigenstate of $-{\rm i}{\bf{\Sigma}}\cdot\nabla$, then

where $ \epsilon=\sqrt{\alpha^{2}+p_{z}^{2}} $ is the magnitude of the total momentum and $ s=\pm1 $ correspond to the two opposite helicities. From Eq. (B11), one can obtain $-{\rm i}{\boldsymbol{\sigma}}\cdot\nabla\psi_{1}=s\epsilon\psi_{1}$, which leads to $A=\dfrac{\rm i}{\alpha}(s\epsilon-p_{z})$ and

Finally, we obtain the eigenfunctions and corresponding eigen-energy as follows,

where $ X=\sqrt{M^{2}+\epsilon^{2}} $, and $ \lambda=\pm1 $ correspond to the positive and negative modes. All eigenfunctions are orthonormal,

The Dirac equation in a uniformly rotating system with angular velocity $ {\boldsymbol{\omega}}=\omega{\boldsymbol{e}}_{z} $ can be written as [28, 31]

Compared with the free case in Sec. B.1, it can be observed that the eigenfunctions of Eq. (B16) are the same as the free case but with an energy shift $ \Delta E=-j\omega $. Now, we consider a uniformly rotating system of massive Dirac fermions with angular velocity $ {\boldsymbol{\omega}}=\omega{\boldsymbol{e}}_{z} $, where the interaction among fermions is ignored. This system is in equilibrium with a reservoir, which keeps a constant temperature T and constant chemical potential μ. In the following, we will calculate the axial vector current $ J_{A}^{\mu} $ of this system. According to the rotational symmetry along the z-axis of the system, we can obtain $ J_{A}^{x}=J_{A}^{y}=0 $. Due to the absence of axial chemical potential $ \mu_{5} $ in our formalism, $ J_{A}^{0} $ vanishes [24]. The unique non-zero component is $ J_{A}^{z} $. From the approach of statistical mechanics used in Refs. [20, 22], one can obtain

where the Fermi-Dirac distribution has been inserted, and $ \beta=1/T $. Making use of the following series for Bessel function $ J_{n}(x) $ with $ n\in\mathbb{N} $,

Equation (B17) becomes

where we have defined four dimensionless quantities, $ \rho=rT $, $ \Omega=\omega/T $, $ \alpha=\mu/T $, and $ c=M/T $. $ C_{N,n} $, $ \mathrm{I}_{N}(\alpha,c) $ are defined as

The coefficient $ C_{N,n} $ can also be expressed as follows:

where we have used the variable transformation $x={\rm e}^{t}$ in the second line. According to the Taylor expansion of $ \sinh t $, one can readily show that $ C_{N,n}=0 $ for $ n<N $. In principle, one can calculate $ C_{N,n} $ for any $ n\geqslant N $ from Eq. (B22). For example, for $ n=N,N+1 $, one can obtain

According to the calculation method in the appendixes of the recent articles by some of us [6, 36], the integral $ \mathrm{I}_{N}(\alpha,c) $ in Eq. (B21) can be expanded at $ c=0 $ as follows,

with $ D_{N,l}(\alpha) $ expanded at $ \alpha=0 $ as

Plugging Eqs. (B24) and (B25) into Eq. (B19), one can obtain the series expansion of $ J_{A}^{z} $ at $ \rho=0 $, $ \Omega=0 $, $ \alpha=0 $, $ c=0 $ or $ r=0 $, $ \omega=0 $, $ M=0 $, $ \mu=0 $ as follows:

If we only keep the linear term of Ω and set $ \alpha=0 $ in Eq. (B26), then

For the massless fermion case, we can obtain an analytic expression for $ J_{A}^{z} $,

which is divergent as the speed-of-light surface is approached [28].

The Kubo formula relates the axial vortical conductivity to the static Fourier component of the correlation between the axial vector current $ J_{A}^{i} $ and the stress tensor $ T^{0j} $ via

in the limit $ {\boldsymbol{k}}\to0 $. Ignoring the interactions, the LHS is represented by the one-loop thermal diagram in Fig. 4 [9, 14, 23, 37]. Calculating the thermal diagram with the Matsubara formulation, we have

Figure 4.  One-loop correction to the vortical conductivity [23].

where the Matsubara frequency $v_{n}=(2n+1)\pi T-{\rm i}\mu$. The summation and integral in Eq. (C2) appear UV-divergent, and we apply the dimensional regularization by extending the spatial components of the loop momentum from 3-dimensional to D-dimensional, i.e.,

It is straightforward to obtain that

where the D-dimensional solid angle $\omega_{D}={2\pi^{D}}/{\Gamma\left(\dfrac{D}{2}\right)}$, and $ B\left(x,y\right) $ is the beta function. The last line of Eq. (C4) separates σ into two terms, where $ \sigma|_{M=0} $ is the vortical conductivity of massless fermions, and $ \Delta\sigma $ is the mass correction of the vortical conductivity. Here,

and

With the dimensionality $ D=3-\epsilon $ and taking the limit $ \epsilon\to 0 $, we find

Upon expanding $ \Delta\sigma $ in the power of $ M^2 $, the leading term, the $ M^2 $ term, is of the form $ 0\times\infty $ at $ D=3 $, and the limit has to be taken carefully. Let $ c_{D} $ be the coefficient of $ M^{2} $; we have

For a higher power of $ M^2 $, however, the naive limit works. Taken together, we obtain the limit

Adding Eq. (C7) and Eq. (C9), we replicate Eq. (112). We have also verified that the same result emerges with the Pauli-Villars regularization.