Revisiting the spin effects induced by thermal vorticity^*

Spin polarization effects have been observed in heavy-ion collisions at RHIC [1−4] and LHC [5, 6] since their pioneering theoretical prediction [7−9]. Relevant reviews on spin effects in heavy ion collisions can be found in Refs. [10−19]. However, some recent measurement results are contradictory to the theoretical calculations. This is the case of the local longitudinal polarization of hyperons [20] and spin alignment of vector mesons [8]. These spin puzzles have prompted extensive theoretical research [21−43] on these topics. As a result of this research, new physical mechanisms have been proposed to explain the aforementioned contradictions, such as the shear contribution for the local longitudinal polarization of hyperons [25−28] and strong force fields for the spin-alignment of vector mesons [32−35]. Although new physical mechanisms are being employed to interpret these unexpected results, it is also necessary to revisit the original theoretical methods which produced such discrepancies to determine whether they could be modified and improved. Such revisiting process is indispensable to quantitatively elucidate the real physical mechanisms underlying the spin polarization effects in heavy-ion collisions.

For the spin polarization of the hyperon, the numerical prediction in the formalism of relativistic hydrodynamics is based on the spin Cooper-Frye formula [44], which measures the mean spin vector $S^\mu(k)$ by integrating the local spin vector $S^\mu(x,k)$ over the freeze-out surface $\Sigma_\alpha$ in heavy-ion collisions:

$$S^\mu(k) = \frac{\int {\rm d} \Sigma_\alpha k^\alpha S^\mu(x,k) f(x,k)} {\int {\rm d}\Sigma_\alpha k^\alpha f(x,k) }, \tag{ 1 }$$

where $f(x,k)$ is a single-particle distribution function at the space-time point $x^\mu=(t,{\bf{x}})$ with four momentum $k^\mu= (E_{\bf{k}},{\bf{k}})$. For the particle with mass m, the energy is given by $E_{\bf{k}}=\sqrt{m^2+{\bf{k}}^2}$. At global thermodynamical equilibrium with small thermal vorticity $\varpi_{\mu\nu}$, the first order contribution for the spin vector $S^\mu(x,k)$ is given by

$$S^\mu(x,k) =\frac{f'_{\rm F}}{8mf_{\rm F}} \epsilon^{\mu\nu\rho\sigma} k_\nu \varpi_{\rho\sigma}\, , \tag{ 2 }$$

where $f_{\rm F}$ is the Fermi-Dirac distribution function:

$$f_{\rm F}=\frac{1}{{\rm e}^{\beta^\mu k_\mu-\alpha}+1} \tag{ 3 }$$

and $f'_{\rm F} = \partial f_{\rm F}/\partial(\beta\cdot k)$. The four-temperature vector $\beta^\mu = (\beta^0, {\boldsymbol\beta})$ is related to the fluid velocity $u^\mu$, with $u^2=1$ and temperature T, by $\beta^\mu = u^\mu/T$, and $\alpha=\mu/T$ denotes the chemical potential μ scaled by temperature T. The thermal vorticity is defined as

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

with components

$$\begin{aligned}[b] \varepsilon^{i}=\; & \varpi^{i 0} = -\frac{1}{2}\left(\partial^i \beta^0 - \partial_t \beta^i\right),\\ \omega^i =\; & \frac{1}{2}\epsilon^{0ijk}\varpi_{jk} = -\frac{1}{2} \epsilon^{0ijk}\partial_j\beta_k,\ \ \end{aligned} \tag{ 5 }$$

or in 3-vector form

$${\boldsymbol\varepsilon} = \frac{1}{2}\left({\pmb{\nabla}} \beta^0 + \partial_t \boldsymbol{\beta}\right), \ \ \ \boldsymbol{\omega} = \frac{1}{2}\pmb{\nabla}\times \boldsymbol{\beta}. \tag{ 6 }$$

After transforming mean spin vector $S^\mu=(S^0,{\bf{S}})$ in Eq. (1) into the rest frame of the particle with momentum ${\bf{k}}$

$$S^{* \mu} = \left(0 , {\bf{S}}^* \right),\ \ \ {\bf{S}}^* = {\bf{S }} - \frac{({\bf{S}}\cdot {\bf{k}}) {\bf{k}}}{E_{\bf{k}}(E_{\bf{k}}+m)}, \tag{ 7 }$$

the final polarization P along some quantization direction ${\bf{n}}_3$ is given by

$$P={\bf{n}}_3\cdot {\bf{S}}^* = {\bf{n}}_3\cdot {\bf{S }} - \frac{({\bf{S}}\cdot {\bf{k}})({\bf{n}}_3\cdot {\bf{k}})}{E_{\bf{k}}(E_{\bf{k}}+m)}. \tag{ 8 }$$

According to the specific expression in Eq. (2) for a local spin vector, the 3-vector form reads

$$S^0(x,k)=\frac{f'_{\rm F}}{4m f_{\rm F}} {\boldsymbol\omega}\cdot {\bf{k}},\ \ \ {\bf{S}}(x,k)=-\frac{f'_{\rm F}}{4m f_{\rm F}}\left( E_{\bf{k}}{\boldsymbol\omega} - {\boldsymbol\varepsilon}\times {\bf{k}}\right). \tag{ 9 }$$

Then, the local polarization is given by

$$P(x,k)= -\frac{ f'_{\rm F}}{2f_{\rm F}}\cdot\frac{E_{\bf{k}}}{2m} \left[{\boldsymbol\omega}-\frac{({\boldsymbol\omega}\cdot {\bf{k}}){\bf{k}}}{E_{\bf{k}} (E_{\bf{k}} +m)} - \frac{{\boldsymbol\varepsilon}\times {\bf{k}}}{E_{\bf{k}}}\right]\cdot{\bf{n}}_3 . \tag{ 10 }$$

It is well known that the polarization P can be obtained directly from the particle distribution $f_{rs}(x,k)$ with indexes $r,s=\pm 1$ (sometimes expressed as $r,s=\pm$ for brevity) corresponding to the spin $\pm 1/2$ along the spin quantization direction ${\bf{n}}_3$

$$P(x,k)=\frac{f_{+,+}(x,k)-f_{-,-}(x,k)}{f_{+,+}(x,k)+f_{-,-}(x,k)}. \tag{ 11 }$$

Hence, we can also calculate the final polarization in heavy-ion collisions with

$$P(k) = \frac{\int {\rm d}\Sigma_\alpha k^\alpha P(x,k) f(x,k)} {\int {\rm d}\Sigma_\alpha k^\alpha f(x,k) }, \tag{ 12 }$$

where $f(x,k)\equiv f_{++}(x,k)+f_{-}(x,k)$ expresses the sum of spin up and spin down along the direction ${\bf{n}}_3$. We will demonstrate that the polarization expressed by Eq. (12) is different from that expressed by Eq. (1).

For the spin polarization of the vector meson, theoretical predictions focus on the spin alignment, and there is no similar formula to that of Eq. (1) for the vector meson yet. Most predictions rely on the quark coalescence model. In the formalism of relativistic hydrodynamics, the measured spin density matrix can be calculated from the particle distribution function $f_{rs}(x,k)$ with spin indexes $r,s = 0,\pm 1$

$$\rho_{rs}(k) = \frac{\int {\rm d}\Sigma_\alpha k^\alpha f_{rs}(x,k) } {\int {\rm d}\Sigma_\alpha k^\alpha f(x,k) } , \tag{ 13 }$$

where $f(x,k)\equiv f_{11}(x,k)+ f_{00}(x,k) + f_{-1-1}(x,k)$ expresses the sum of all the diagonal components along the direction ${\bf{n}}_3$. We will derive a specific expression for this density matrix by calculating the particle distribution $f_{rs}(x,k)$. We will show that the spin alignment receives only a second-order contribution from acceleration or vorticity.

Note that the particle distributions $f_{rs}(x,k)$ with spin will be the crucial elements in Eqs. (12) and (13) rather than the spin vector in Eqs. (1) and (2). Given that we only revisit the spin polarization by thermal vorticity in this study, we calculate the particle distributions with spin for a free particle in global equilibrium with thermal vorticity. We assume that these results still dominate in local equilibrium. The exact equilibrium distributions with thermal vorticity have been recently obtained by analytical continuation in Refs. [45−47]. In this study, we calculate these distribution functions in a more direct and usual manner, and expand them specifically in terms of vorticity and acceleration in first or second order.

In Sec. II, we calculate the particle distributions for the scalar field in global equilibrium with thermal vorticity. We obtain the particle distributions for the Dirac field in Sec. III and for the vector field in Sec. IV. A summary of our results is presented in Sec. V. In this paper, we use the metric $g^{\mu\nu}=\mathrm{diag}(1,-1,-1,-1)$ and Levi-Civita tensor $\epsilon^{0123}=1$.

Let us first review well-known results for the scalar field. The Lagrange density of charged scalar field reads

$$\mathcal{L}=\partial^\mu\phi^\dagger \partial_\mu\phi -m^2\phi^\dagger\phi. \tag{ 14 }$$

The Euler-Lagrangian equation is just the Klein-Gordon equation:

$$\partial^\mu \partial_\mu\phi + m^2\phi = 0, \quad \partial^\mu \partial_\mu\phi^\dagger + m^2\phi^\dagger = 0. \tag{ 15 }$$

From the Noether' theorem, we have the conserved charge current,

$$j^{\mu} = i\left( \partial^\mu\phi^\dagger\phi - \phi^\dagger \partial^\mu\phi\right), \tag{ 16 }$$

the canonical energy-momentum tensor,

$$T^{\mu\nu} =\partial^\mu\phi^\dagger \partial^\nu\phi + \partial^\nu\phi^\dagger \partial^\mu\phi - g^{\mu\nu}\left(\partial_\alpha\phi^\dagger \partial^\alpha\phi -m^2\phi^\dagger\phi\right), \tag{ 17 }$$

and the angular momentum density,

$$\mathcal{M}^{\mu\alpha\beta} = x^\alpha T^{\mu\beta} - x^\beta T^{\mu\alpha} . \tag{ 18 }$$

We can expand the Klein-Gordon field in terms of the annihilation and creation operators:

$$\begin{aligned}[b] \phi(x) =\; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\frac{1}{\sqrt{2E_{\bf{p}}}} \left[a_{\bf{p}} {\rm e}^{-{\rm i}p\cdot x} + b_{\bf{p}}^\dagger {\rm e}^{{\rm i}p\cdot x} \right],\\ \phi^\dagger(x) =\; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\frac{1}{\sqrt{2E_{\bf{p}}}} \left[b_{\bf{p}} {\rm e}^{-{\rm i}p\cdot x} + a_{\bf{p}}^\dagger {\rm e}^{{\rm i}p\cdot x} \right], \end{aligned} \tag{ 19 }$$

where $p^\mu = (p^0,{\bf{p}})$ with $p^0=E_{\bf{p}}$ and the creation and annihilation operators obey the commutation rules:

$$\left[a_{\bf{p}} , a_{\bar{\bf{p}}}^\dagger\right]= \left[b_{\bf{p}} , b_{\bar{\bf{p}}}^\dagger\right] =(2\pi)^3 \delta^{(3)}({\bf{p}}-\bar{\bf{p}}). \tag{ 20 }$$

Inserting Eq. (19) into Eqs. (16) and (17) and integrating over the whole space, we obtain the conserved charge and energy-momentum:

$$Q = \int {\rm d}^3 {\bf{x}} \,j^{0} = \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3} \left( a_{\bf{p}}^{\dagger} a_{\bf{p}} - b_{\bf{p}} b_{\bf{p}}^{\dagger} \right), \tag{ 21 }$$

$$P^\mu = \int {\rm d}^3 {\bf{x}} \, T^{0\mu} = \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3} p^\nu \left( a_{\bf{p}}^{\dagger} a_{\bf{p}} + b_{\bf{p}} b_{\bf{p}}^{\dagger} \right). \tag{ 22 }$$

The angular-momentum tensor is defined by

$$J^{\mu\nu}= \int {\rm d}^3 {\bf{x}} \, \mathcal{M}^{0\mu\nu}, \tag{ 23 }$$

with the components

$$\begin{aligned}[b] K^{i}=\; & J^{0i} = \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3} \Big[\sqrt{E_{\bf{p}}}a_{\bf{p}}^{\dagger} {\rm i}\partial^i_p \left( \sqrt{E_{\bf{p}}} a_{\bf{p}}\right) \\ & - \sqrt{E_{\bf{p}}} b_{\bf{p}} {\rm i} \partial^i_p \left( \sqrt{E_{\bf{p}}} b_{\bf{p}}^{\dagger}\right)\Big], \end{aligned} \tag{ 24 }$$

$$J^{i}=\frac{1}{2} \epsilon^{0ijk} J_{jk} = \epsilon^{0ijk} \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3} p^k\left[a_{\bf{p}}^{\dagger}{\rm i}\partial_j^p a_{\bf{p}} - b_{\bf{p}} {\rm i} \partial_j^p b_{\bf{p}}^{\dagger}\right], \tag{ 25 }$$

or in 3-vector form,

$$\begin{aligned}[b] {\bf{K}}=\; & - \int \frac{{\rm d}^3{\bf{p}}}{(2\pi)^3} \Big[ a_{\bf{p}}^{\dagger} \sqrt{E_{\bf{p}}}\pmb{\nabla}_p \left( \sqrt{E_{\bf{p}}}a_{\bf{p}}\right)\\ & - b_{\bf{p}}\sqrt{E_{\bf{p}}}\pmb{\nabla}_p\left( \sqrt{E_{\bf{p}}} b_{\bf{p}}^{\dagger}\right)\Big], \end{aligned} \tag{ 26 }$$

$${\bf{J}} =-\int \frac{{\rm d}^3{\bf{p}}}{(2\pi)^3} \left[a_{\bf{p}}^{\dagger} \left({\bf{p}}\times {\rm i}\,\pmb{\nabla}_p\right) a_{\bf{p}} - b_{\bf{p}} \left({\bf{p}}\times {\rm i}\,\pmb{\nabla}_p\right) b_{\bf{p}}^{\dagger}\frac{}{}\right], \tag{ 27 }$$

where

$$\partial_p^i = \frac{\partial}{\partial p_i} = - \frac{\partial}{\partial p^i}, \ \ \ \ \ \pmb{\nabla}_p = \left(\frac{\partial}{\partial p^1},\frac{\partial}{\partial p^2}, \frac{\partial}{\partial p^3} \right). \tag{ 28 }$$

Now, let us analyze the particle distribution function, which is defined by [48]

$$f(x,k)\equiv \frac{1}{(2\pi)^3}\int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t + {\rm i}{\bf{q}}\cdot{\bf{x}}} \langle a^{\dagger}_{{\bf{k}}-{\bf{q}}/{2}} a_{ {\bf{k}}+{\bf{q}}/{2}}\rangle, \tag{ 29 }$$

where $\langle Q \rangle = \rm{Tr}\left(\rho Q\right)$ denotes the ensemble average of some operator Q with density matrix ρ. Here, we choose the density operator in global equilibrium [49]:

$$\rho=\frac{1}{Z} {\rm exp}\left(-b_\mu P^\mu +\alpha Q + \frac{1}{2}\omega_{\mu\nu}J^{\mu\nu} \right), \tag{ 30 }$$

where Z is the partition function, $b_\mu$ is a constant timelike vector, α is a constant scalar, and $\omega_{\mu\nu}$ is a constant antisymmetric tensor. The operators $P^\mu$, Q, and $J^{\mu\nu}$ denote the energy-momentum, charge, and angular momentum tensors, respectively. The contribution from the angular momentum tensor can be rewritten as

$$\frac{1}{2}\omega_{\mu\nu}J^{\mu\nu} = {{\boldsymbol\varepsilon}}\cdot{\bf{K}} + {{\boldsymbol\omega}}\cdot{\bf{J}}. \tag{ 31 }$$

The density operator expressed by Eq. (30) is a general result, valid for any field. For the free charged scalar field, $P^\mu$, Q, and $J^{\mu\nu}$ will take the specific forms in Eqs. (21), (22), (24), and (25).

To obtain the distribution function $f(x,k)$, we first calculate

$$\langle a^{\dagger}_{\bf{p}} a_{ \bar{\bf{p}}}\rangle =\frac{1}{Z} {\rm Tr}\left[ {\rm exp}\left(-b_\mu P^\mu +\alpha Q + \frac{1}{2}\omega_{\mu\nu}J^{\mu\nu} \right) a^{\dagger}_{\bf{p}} a_{ \bar{\bf{p}}} \right]. \tag{ 32 }$$

From Eqs. (21), (22), (24), and (25), we obtain the following commutation relations:

$$\left[Q, a_{\bf{p}}^{\dagger}\right] = a_{{\bf{p}}}^{\dagger},\ \ \ \left[b_\mu P^\mu, a_{\bf{p}}^{\dagger}\right] = b_\mu p^\mu a_{{\bf{p}}}^{\dagger},\ \ \ \left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}} + {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{\dagger}\right] = \Lambda_{\bf{p}}{a_{\bf{p}}}^\dagger. \tag{ 33 }$$

In the last equation above, $\Lambda_{\bf{p}}$ denotes an operator defined by

$$\Lambda_{\bf{p}} =E_{\bf{p}} \pmb{\mathcal{ E}}_{\bf{p}} \cdot {\rm i}\,{\pmb\nabla}_{p} + \frac{\rm i}{2E_{\bf{p}}}{\boldsymbol\varepsilon}\cdot {\bf{p}}, \tag{ 34 }$$

where we introduced an effective acceleration vector $\pmb{\mathcal{E}}_{\bf{p}}$

$${\bf\mathcal{E}}_{\bf{p}} = {\boldsymbol\varepsilon} +\frac{1}{E_{\bf{p}}} {\boldsymbol\omega}\times {\bf{p}}. \tag{ 35 }$$

Using the Baker-Hausdorff formula

$$\begin{aligned}[b] {\rm e}^{- A} B {\rm e}^{ A} =\; & B +\frac{(-1)^1}{1!} \left[A,B\right] + \frac{(-1)^2}{2!}\left[A,\left[A,B\right]\right] \\ & +\frac{(-1)^3}{3!}\left[A,\left[A, \left[A,B\right]\right]\right] +\cdots \end{aligned}$$

we have the identity

$$\begin{aligned}[b] & {\rm exp}\left(-b_\mu P^\mu +\alpha Q + \frac{1}{2}\omega_{\mu\nu}J^{\mu\nu} \right) a_{\bf{p}}^{\dagger} {\rm exp}\Big(b_\mu P^\mu -\alpha Q \\ & - \frac{1}{2}\omega_{\mu\nu}J^{\mu\nu} \Big) = {\rm e}^{-b_\mu p^\mu + \alpha + \Lambda_{\bf{p}}}a_{\bf{p}}^{\dagger}. \end{aligned}$$

which leads to

$$\langle a^{\dagger}_{\bf{p}} a_{ \bar{\bf{p}}}\rangle = {\rm e}^{-b_\mu p^\mu + \alpha + \Lambda_{\bf{p}}}\langle a_{ \bar{\bf{p}}} a^{\dagger}_{\bf{p}}\rangle . \tag{ 36 }$$

From the commutation relation for the Klein-Gorden field, we also have

$$\langle a_{ \bar{\bf{p}}} a^{\dagger}_{\bf{p}} \rangle =(2\pi)^3\delta\left({\bf{p}}-\bar{\bf{p}} \right) + \langle a^{\dagger}_{\bf{p}} a_{ \bar{\bf{p}}}\rangle. \tag{ 37 }$$

Substituting Eq. (37) into Eq. (36), we obtain

$$\langle a^{\dagger}_{\bf{p}} a_{ \bar{\bf{p}}}\rangle =(2\pi)^3\left[\frac{1}{{\rm e}^{b_\mu p^\mu - \alpha - \Lambda_{\bf{p}}}-1}\right]\delta\left({\bf{p}}-\bar{\bf{p}} \right). \tag{ 38 }$$

Replacing ${\bf{p}}$ and $\bar{\bf{p}}$ with ${{\bf{k}}-{\bf{q}}/{2}}$ and ${{\bf{k}}-{\bf{q}}/{2}}$, respectively, and inserting Eq. (38) into Eq. (29), we obtain the distribution function

$$\begin{aligned}[b] f(x,k) =\; & \int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+{\rm i}{\bf{q}}\cdot{\bf{x}}} \\ & \times\frac{1}{{\rm e}^{b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha - \Lambda_{{\bf{k}}-{\bf{q}}/2}}-1} \delta\left({\bf{q}}\right), \end{aligned} \tag{ 39 }$$

where the operator $\Lambda_{{\bf{k}}-{\bf{q}}/2}$ reads

$$\Lambda_{{\bf{k}}-{\bf{q}}/2} =E_{{\bf{k}}-{\bf{q}}/2} \pmb{\mathcal{E}}_{{\bf{k}}-{\bf{q}}/2} \cdot {\rm i}\left(\frac{1}{2}\pmb{\nabla}_{k} - \pmb{\nabla}_{q}\right) + \frac{\rm i}{2E_{{\bf{k}}-{\bf{q}}/2}}{\boldsymbol\varepsilon}\cdot \left( {\bf{k}}-\frac{\bf{q}}{2}\right). \tag{ 40 }$$

We assume that ${{\boldsymbol\varepsilon}}$ and ${{\boldsymbol\omega}}$ are both small variables and expand the distribution function as the Taylor series of these variables. In particular, we employ the following expansion:

$$\begin{aligned}[b] \frac{1}{{\rm e}^{X+Y} - 1} =\; & f_{\rm B}(X)+f'_{\rm B}(X)Y - \frac{1}{2}f''_{\rm B}(X)C +\frac{1}{2}f''_{\rm B}(X)Y^2 \\ & - \frac{1}{6}f'''_{\rm B} (X)YC - \frac{1}{3}f'''_{\rm B}(X)CY +\frac{1}{8}f''''_{\rm B}(X)C^2+\cdots \end{aligned} \tag{ 41 }$$

where $C=[X,Y]$ and $f_{\rm B}(X)$ denotes the Bose-Einstein distribution function $1/({\rm e}^{X}-1)$. This expansion is valid up to the second order of Y. In our case, we can identify X and Y as

$$X=b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha,\ \ \ \ Y=-\Lambda_{{\bf{k}}-{\bf{q}}/2} , \tag{ 42 }$$

$$C=[X,Y]={\rm i}\left[\left(b_0 {\boldsymbol\varepsilon}+ {\boldsymbol\omega}\times {\bf{b}}\right) \cdot \left( {\bf{k}}-\frac{\bf{q}}{2}\right) - ({\bf{b}}\cdot{\boldsymbol\varepsilon})E_{{\bf{k}}-{\bf{q}}/2}\right]. \tag{ 43 }$$

Now, we can calculate the distribution functions order by order after integrating over the momentum ${\bf{q}}$. The zeroth-order result is trivial; it is given by the Bose-Einstein distribution function:

$$f^{(0)}(x,k) =f_{\rm B}(b\cdot k -\alpha) , \tag{ 44 }$$

The first-order result is also simple:

$$f^{(1)}(x,k) = f'_{\rm B}(b\cdot k -\alpha) \left[E_{\bf{k}}{\boldsymbol\varepsilon}\cdot{\bf{x}} - {\bf{k}}\cdot ({\boldsymbol\varepsilon}t + {\boldsymbol\omega}\times{\bf{x}})\right] . \tag{ 45 }$$

The second-order result is more involved:

$$\begin{aligned}[b] f^{(2)}(x,k) =\; & \frac{1}{2} f''_{\rm B}(b\cdot k -\alpha) \left[E_{\bf{k}}{\boldsymbol\varepsilon}\cdot{\bf{x}} - {\bf{k}}\cdot ({\boldsymbol\varepsilon}t + {\boldsymbol\omega}\times{\bf{x}})\right]^2\\ & + \frac{1}{8}f''_{\rm B} (b\cdot k -\alpha) \left[\frac{({\bf{k}}\cdot {\boldsymbol\varepsilon})^2}{E_{\bf{k}}^2} + \frac{2{\bf{k}}\cdot ({\boldsymbol\omega}\times{\boldsymbol\varepsilon} )}{E_{\bf{k}}} -2{\boldsymbol\omega}^2 \right]\\ & + \frac{1}{12}f'''_{\rm B} (b\cdot k -\alpha) \left(E_{\bf{k}}{\boldsymbol\varepsilon}+{\boldsymbol\omega}\times{\bf{k}} \right) \cdot \left(b_0 {\boldsymbol\varepsilon}+ {\boldsymbol\omega}\times {\bf{b}}\right). \end{aligned} \tag{ 46 }$$

Note that there exist some terms which depend on the time t and space coordinates ${\bf{x}}$. When t or ${\bf{x}}$ are large, our expansion will not hold. We can absorb these terms into the vector $b^\mu$ and obtain a new vector $\beta^\mu$, which can be regarded as the zeroth order contribution:

$$\beta^0 = b^0 + {\boldsymbol\varepsilon}\cdot{\bf{x}},\ \ \ \ {\boldsymbol\beta} = {\bf{b}} + {\boldsymbol\varepsilon} t + {\boldsymbol\omega}\times {\bf{x}} . \tag{ 47 }$$

Summing over the particle distributions up to the second order, we finally obtain

$$\begin{aligned}[b] f(x,k) =\; & f_{\rm B}(\beta\cdot k -\alpha) + \frac{1}{8}f''_{\rm B} (\beta\cdot k -\alpha)\\ & \times \left[\frac{({\bf{k}}\cdot {\boldsymbol\varepsilon})^2}{E_{\bf{k}}^2} + \frac{2{\bf{k}}\cdot ({\boldsymbol\omega}\times{\boldsymbol\varepsilon} )}{E_{\bf{k}}} -2{\boldsymbol\omega}^2 \right]\\ & + \frac{1}{12}f'''_{\rm B}(\beta\cdot k -\alpha) \left(E_{\bf{k}}{\boldsymbol\varepsilon}+{\boldsymbol\omega}\times{\bf{k}} \right) \cdot \left(b^0 {\boldsymbol\varepsilon}+ {\boldsymbol\omega}\times {\bf{b}}\right). \end{aligned} \tag{ 48 }$$

If we identify $\beta^\mu$ as the inverse temperature vector $u^\mu/T$, we obtain the well-established conclusion that the spin chemical potential $\omega^{\mu\nu}$ is equal to the thermal vorticity $\varpi^{\mu\nu}$ at global equilibrium.

Let us now consider the Dirac fermions with spin-1/2. The Lagrangian for the free Dirac field is given by

$$\mathcal{L}=\bar\psi({\rm i}\gamma^\mu \partial_\mu -m)\psi \tag{ 49 }$$

from which we can obtain the Dirac equations,

$${\rm i}\gamma^\mu \partial_\mu \psi(x)-m \psi(x) =0,\ \ \ {\rm i}\partial_\mu \bar\psi(x)\gamma^\mu +m\bar\psi(x) =0, \tag{ 50 }$$

the electric currents,

$$j^{\mu} = \bar \psi \gamma^\mu \psi , \tag{ 51 }$$

the canonical energy-momentum tensor,

$$T^{\mu\nu} = \frac{\rm i}{2}\bar\psi \gamma^\mu \left[ \overrightarrow{\partial}^\nu - \overleftarrow{\partial}^\nu \right]\psi, \tag{ 52 }$$

and the angular momentum tensor density,

$$\begin{aligned}[b] \mathcal{M}^{\mu,\alpha\beta} = \; & x^\alpha T^{\mu\beta} - x^\beta T^{\mu\alpha} + \frac{1}{2}\bar\psi \left\{\gamma^\mu, \frac{\sigma^{\alpha\beta}}{2}\right\}\psi\\ =\; & x^\alpha T^{\mu\beta} - x^\beta T^{\mu\alpha} - \frac{1}{2} \epsilon^{\mu\alpha\beta\nu} \bar\psi \gamma_\nu \gamma_5 \psi. \end{aligned} \tag{ 53 }$$

Then, the charge, energy-momentum, and angular momentum tensor read, respectively,

$$Q= \int {\rm d}^3 {\bf{x}} \, j^{0}, \ P^\mu = \int {\rm d}^3 {\bf{x}} \, T^{0\mu}, J^{\mu\nu} = \int {\rm d}^3 {\bf{x}} \, \mathcal{M}^{0\mu\nu}, \tag{ 54 }$$

Expanding the free Dirac field in terms of the annihilation and creation operators, we obtain

$$\begin{aligned}[b] \psi(x) =\; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\frac{1}{\sqrt{2E_{\bf{p}}}}\sum\limits_s \left[a^s_{\bf{p}} u^s({\bf{p}}) {\rm e}^{-{\rm i}p\cdot x} + {b^s_{\bf{p}}}^\dagger v^s({\bf{p}}) {\rm e}^{{\rm i}p\cdot x} \right];\\ \bar\psi(x) =\; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\frac{1}{\sqrt{2E_{\bf{p}}}}\sum\limits_s \left[b^s_{\bf{p}} \bar v^s({\bf{p}}) {\rm e}^{-{\rm i}p\cdot x} + {a^s_{\bf{p}}}^\dagger \bar u^s({\bf{p}}) {\rm e}^{{\rm i}p\cdot x} \right], \end{aligned} \tag{ 55 }$$

where the creation and annihilation operators obey the anticommutation rules

$$\left\{a^s_{\bf{p}} , {a^r_{\bar{\bf{p}}}}^\dagger\right\}= \left\{b^s_{\bf{p}} , {b^r_{\bar{\bf{p}}}}^\dagger\right\} =(2\pi)^3 \delta^{(3)}({\bf{p}}-\bar{\bf{p}})\delta^{sr} \tag{ 56 }$$

and the indexes $s,r=\pm 1$ denote the spin $\pm 1/2$. Substituting the expansion described by Eq. (55) into the conserved charges described by Eq. (54), we have

$$Q=\int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\sum\limits_s \left[ {a^s_{\bf{p}}}^\dagger a^s_{\bf{p}}+ b^s_{\bf{p}} {b^s_{\bf{p}}}^\dagger \right], \tag{ 57 }$$

$$P^\mu = \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\sum\limits_s \left[p^\mu {a^s_{\bf{p}}}^\dagger a^s_{\bf{p}}- p^\mu b^s_{\bf{p}} {b^s_{\bf{p}}}^\dagger \right] \tag{ 58 }$$

and the angular momentum tensor with the components as defined in Eqs. (24) and (25) is expressed as

$$\begin{aligned}[b] K^{i}=\; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\frac{1}{2} \sum\limits_s \sum\limits_r \left\{{ v^{s\dagger}({\bf{p}})} b^s_{\bf{p}} ( {\rm i} \partial^i_p ) \left[ {b^r_{\bf{p}}}^\dagger v^r({\bf{p}}) \right]\right.\\ & \left. +{u^{s\dagger}({\bf{p}})} {a^s_{\bf{p}}}^\dagger ( {\rm i} \partial^i_p) \left[ a^r_{\bf{p}} u^r({\bf{p}})\right] \frac{}{} \right\}, \end{aligned} \tag{ 59 }$$

$$\begin{aligned}[b] J^i = \; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3} \frac{1}{ 2 E_{\bf{p}} }\sum\limits_s \sum\limits_r \left\{ v^{s\dagger}({\bf{p}}) b^s_{\bf{p}} \left[\epsilon^{ijk} p^j ({\rm i} \partial^k_p) +\frac{1}{2}\Sigma^i \right]\right.\\ & \times \left[{b^{r\dagger}_{\bf{p}}} v^r({\bf{p}})\right] + u^{s\dagger}({\bf{p}}) {a^{s\dagger}_{\bf{p}}} \end{aligned}$$

$$\begin{aligned}[b] \left. \times \left[\epsilon^{ijk} p^j ( {\rm i}\partial^k_p) +\frac{1}{2}\Sigma^{i} \right] \left[ a^r_{\bf{p}} u^r({\bf{p}})\right]\right\} \end{aligned} \tag{ 60 }$$

or in 3-vector form

$$\begin{aligned}[b] {\bf{K}}=\; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3}\frac{1}{2} \sum\limits_s \sum\limits_r \left\{{ v^{s\dagger}({\bf{p}})} b^s_{\bf{p}} (- {\rm i}\,\pmb{ \nabla}_p ) \left[ {b^r_{\bf{p}}}^\dagger v^r({\bf{p}}) \right]\right.\\ & \left. +{u^{s\dagger}({\bf{p}})} {a^s_{\bf{p}}}^\dagger ( - {\rm i}\,\pmb{ \nabla}_p) \left[ a^r_{\bf{p}} u^r({\bf{p}})\right] \right\}, \end{aligned} \tag{ 61 }$$

$$\begin{aligned}[b] {\bf{J}} = \; & \int \frac{{\rm d}^3 {\bf{p}}}{(2\pi)^3} \frac{1}{ 2 E_{\bf{p}} }\sum\limits_s \sum\limits_r \left\{ v^{s\dagger}({\bf{p}}) b^s_{\bf{p}} \left[- {\bf{p}}\times ({\rm i}\,\pmb{\nabla}_p) +\frac{1}{2}\pmb{\Sigma} \right]\right.\\ & \times \left[{b^{r\dagger}_{\bf{p}}} v^r({\bf{p}})\right] + u^{s\dagger}({\bf{p}}) {a^{s\dagger}_{\bf{p}}}\\ & \times\left. \left[- {\bf{p}}\times ({\rm i}\,\pmb{\nabla}_p) +\frac{1}{2}\pmb{\Sigma} \right] \left[ a^r_{\bf{p}} u^r({\bf{p}})\right]\right\} , \end{aligned} \tag{ 62 }$$

where ${\pmb \Sigma}$ is defined from the Pauli matrix $\pmb{ \sigma}$

$${\boldsymbol \Sigma} = \left( \begin{array}{cc} {\boldsymbol \sigma} & 0 \\ 0 & {\boldsymbol \sigma} \\ \end{array} \right) . \tag{ 63 }$$

Then, we can obtain the commutation relations

$$\begin{aligned}[b] & \left[Q, a_{\bf{p}}^{s\dagger}\right] = a_{{\bf{p}}}^{s\dagger},\ \ \ \left[b_\mu P^\mu, a_{\bf{p}}^{s\dagger}\right] = b_\mu p^\mu a_{{\bf{p}}}^{s\dagger},\\ & \left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}} + {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{s\dagger}\right] =\sum\limits_r \Lambda_{\bf{p}}^{sr}{a_{\bf{p}}}^{r\dagger} \, , \end{aligned} \tag{ 64 }$$

where the operator $\Lambda_{\bf{p}}^{sr}$ with spin index is given by

$$\begin{aligned}[b] \Lambda_{\bf{p}}^{sr} =\; & \frac{1}{2E_{\bf{p}}}{u^{r\dagger}({\bf{p}})} u^s({\bf{p}}) (E_{\bf{p}}{\boldsymbol\varepsilon} - {\bf{p}}\times {\boldsymbol\omega})\cdot {\rm i}\,{\pmb{ \nabla}}_p \\ & + \frac{1}{2E_{\bf{p}}}{u^{r\dagger}({\bf{p}})} \left(\frac{1}{2}{\boldsymbol\omega}\cdot{\pmb{\Sigma}}\right) u^s({\bf{p}}) \\ & +\frac{1}{2E_{\bf{p}}}(E_{\bf{p}}{\boldsymbol\varepsilon} - {\bf{p}}\times {\boldsymbol\omega}) \cdot \left[{\rm i}\, {\pmb{ \nabla}}_p {u^{r\dagger}({\bf{p}})} \right] u^s({\bf{p}}). \end{aligned} \tag{ 65 }$$

Following the same route to obtain Eq. (36) in the scalar field, we have

$$\langle a^{s\dagger}_{\bf{p}} a^{r}_{ \bar{\bf{p}}}\rangle =\left({\rm e}^{-b_\mu p^\mu + \alpha + \Lambda}\right)^{s s'}\langle a^{r}_{ \bar{\bf{p}}} a^{s'\dagger}_{\bf{p}}\rangle. \tag{ 66 }$$

According to the anticommutation relation for the fermion field, we also have

$$\langle a^r_{ \bar{\bf{p}}} a^{s\dagger}_{\bf{p}} \rangle =(2\pi)^3\delta\left({\bf{p}}-\bar{\bf{p}} \right)\delta^{sr} - \langle a^{s\dagger}_{\bf{p}} a^r_{ \bar{\bf{p}}}\rangle , \tag{ 67 }$$

which leads to

$$\langle a^{s\dagger}_{\bf{p}} a^{r}_{ \bar{\bf{p}}}\rangle =(2\pi)^3\left[\frac{1}{{\rm e}^{b_\mu p^\mu - \alpha - \Lambda}+1}\right]^{sr}\delta\left({\bf{p}}-\bar{\bf{p}} \right). \tag{ 68 }$$

It follows that the distribution function is given by

$$\begin{aligned}[b] f_{rs}(x,k)\equiv\; & \frac{1}{(2\pi)^3}\int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+{\rm i}{\bf{q}}\cdot{\bf{x}}} \langle a^{s\dagger}_{{\bf{k}}-{\bf{q}}/{2}} a^r_{ {\bf{k}}+{\bf{q}}/{2}}\rangle_0\\ =\; & \int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+{\rm i}{\bf{q}}\cdot{\bf{x}}}\\ & \times \left[\frac{1}{{\rm e}^{b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha - \Lambda_{{\bf{k}}-{\bf{q}}/2}}+1}\right]^{sr} \delta\left({\bf{q}}\right). \end{aligned} \tag{ 69 }$$

Further calculations on $\Lambda_{\bf{p}}$ or $\Lambda_{{\bf{k}}-{\bf{q}}/2}$ need specific expressions for the four-component Dirac spinors $u^s({\bf{p}})$

$$u^s({\bf{p}})= \left( \begin{array}{c} \sqrt{p\cdot\sigma} \xi^s\\ \sqrt{p\cdot \bar\sigma}\xi^s \\ \end{array} \right), \tag{ 70 }$$

where $\sigma^\mu = (1, \pmb{\sigma})$ and $\bar\sigma^\mu = (1, -\pmb{\sigma})$, and the two-component spinor $\xi^s$ is chosen as

$$\xi^{+1}=\left(\begin{array}{c} \cos\dfrac{\vartheta}{2}{\rm e}^{-{\rm i}\varphi/2} \\ \sin\dfrac{\vartheta}{2}{\rm e}^{{\rm i}\varphi/2} \end{array} \right),\ \ \ \ \xi^{-1}=\left(\begin{array}{c} -\sin\dfrac{\vartheta}{2}{\rm e}^{-{\rm i}\varphi/2} \\ \cos\dfrac{\vartheta}{2}{\rm e}^{{\rm i}\varphi/2} \end{array} \right) \tag{ 71 }$$

or in a unified form

$$\xi^{s}=(- {\rm i})^{\frac{1-s}{2}}\left(\begin{array}{c} \cos\dfrac{s\vartheta+(1-s)\pi/2}{2}{\rm e}^{-{\rm i}\frac{\varphi+(1-s)\pi/2}{2}} \\ \sin\dfrac{s\vartheta+(1-s)\pi/2}{2} {\rm e}^{{\rm i}\frac{\varphi+(1-s)\pi/2}{2}} \end{array} \right) . \tag{ 72 }$$

They are the eigenstates of the spin operator along the direction ${\bf{n}}_3=(\sin\vartheta \cos\varphi,\; \sin\vartheta \sin \varphi,\; \cos\vartheta)$:

$$\pmb{\sigma}\cdot{\bf{n}}_3 \xi^{s}= s \xi^{s} . \tag{ 73 }$$

It is easy to verify the following relations:

$$\xi^{s}\xi^{r\dagger} =\frac{1}{2}\lambda^{sr}\cdot \bar\sigma,\ \ \ \sqrt{p\cdot\sigma} =\varrho \cdot \sigma,\ \ \ \sqrt{p\cdot\bar\sigma} =\varrho \cdot \bar\sigma , \tag{ 74 }$$

where we have defined the 4-vector $\lambda^{sr,\mu}$ with spin index and 4-vector $\varrho^\mu$ as

$$\begin{aligned}[b] & \lambda^{sr,\mu}=(\lambda^{sr,0},\pmb{\lambda}^{sr})=\left(\delta^{s,r},s \delta^{s,r} {\bf{n}}_3 + {\rm i} s \delta^{-s,r} {\bf{n}}_2 + \delta^{-s,r} {\bf{n}}_1 \right),\\ & \varrho^\mu =\frac{1}{2}\left( \sqrt{E+p} + \sqrt{E-p}, (\sqrt{E+p} - \sqrt{E-p})\hat{\bf{p}}\right) , \end{aligned} \tag{ 75 }$$

where $p=|{\bf{p}}|$. We can rewrite $\lambda^{sr,\mu}$ in matrix form as

$$\lambda^{\mu}=\left(1,{\bf{n}}_3 \sigma_3^T + {\bf{n}}_2\sigma_2^T + {\bf{n}}_1 \sigma_1^T \right), \tag{ 76 }$$

where the superscript T denotes the transpose of a matrix. Here, we have introduced two transverse unit 3-vector, ${\bf{n}}_1$ and ${\bf{n}}_2$, orthogonal to ${\bf{n}}_3$

$${\bf{n}}_2 = \frac{\hat{ \bf z}\times {\bf{n}}_3}{|\hat{ \bf z}\times {\bf{n}}_3|} =\left(-\sin\varphi, \cos\varphi, 0 \right), \tag{ 77 }$$

$${\bf{n}}_1 ={\bf{n}}_2\times {\bf{n}}_3 = \left( \cos\vartheta \cos\varphi, \cos\vartheta \sin\varphi, -\sin\vartheta\right) \tag{ 78 }$$

where $\hat{\bf{z}}$ is the unit vector along the z-axis. When we set ${\bf{n}}_3=\hat{\bf{z}}$, we obtain that ${\bf{n}}_2=\hat{\bf{y}}$ and ${\bf{n}}_1=\hat{\bf{x}}$. The first and second terms in Eq. (65) can be further developed by using the following identities:

$${u^{r\dagger}({\bf{p}})} u^s({\bf{p}})=2E \lambda^{sr,0},\ \ \ {u^{r\dagger}({\bf{p}})}{\pmb{\Sigma}} u^s({\bf{p}}) = 2 m {\boldsymbol\lambda}^{sr} + 4 ({\boldsymbol\varrho}\cdot{\boldsymbol\lambda}^{rs}){\boldsymbol\varrho} . \tag{ 79 }$$

However, to deal with the last term including the derivative with the momentum on the Dirac spinors in Eq. (65), we need to know whether the two-component spinor $\xi^s$ depends on the momentum or not.

A. Polarization along a fixed direction

If the spin quantization direction ${\bf{n}}_3$ does not depend on the momentum ${\bf{p}}$, then the derivative does not act on the two-component spinor $\xi^s$. Using the identities

$${\pmb{ \nabla}}_p\sqrt{p\cdot\sigma} =\frac{1}{2E}\pmb{\varrho} - \frac{1}{2E} \varrho_0\hat{\bf{p}}(\hat{\bf{p}} \cdot\pmb{\sigma}) -\frac{1}{p} (\pmb{\varrho}\cdot \hat{\bf{p}})\left[\pmb{\sigma}-\hat{\bf{p}}(\hat{\bf{p}} \cdot\pmb{\sigma}) \right], \tag{ 80 }$$

$${\pmb{ \nabla}}_p\sqrt{p\cdot\bar\sigma} =\frac{1}{2E}{\boldsymbol\varrho} + \frac{1}{2E} \varrho_0\hat{\bf{p}}(\hat{\bf{p}} \cdot\pmb{\sigma}) +\frac{1}{p} (\pmb{\varrho}\cdot \hat{\bf{p}})\left[\pmb{\sigma}-\hat{\bf{p}}(\hat{\bf{p}} \cdot\pmb{\sigma}) \right], \tag{ 81 }$$

we have

$$\left[ {\pmb{ \nabla}}_p {u^{r\dagger}({\bf{p}})} \right] u^s({\bf{p}}) = \frac{2}{E}\varrho_0\lambda_0^{sr}{\boldsymbol\varrho} +\frac{2}{p} {\rm i} (\hat{\bf{p}}\cdot {\boldsymbol\varrho}){\boldsymbol\varrho} \times {\boldsymbol\lambda}^{sr} . \tag{ 82 }$$

Substituting Eqs. (79) and Eq. (82) into the operator $\Lambda^{sr}_{\bf{p}}$ in Eq. (65), we obtain

$$\Lambda^{sr}_{\bf{p}}= \Lambda_{\bf{p}}\delta^{sr} + \frac{1}{2}{\boldsymbol\Omega}_{\bf{p}} \cdot{\boldsymbol\lambda}^{sr} \tag{ 83 }$$

where $\Lambda_{\bf{p}}$ is given in Eq. (34) and ${\boldsymbol\Omega}_{\bf{p}}$ is an effective vorticity vector defined by

$${\bf\Omega}_{\bf{p}}= \left({\boldsymbol\omega}-\frac{{\boldsymbol\varepsilon}\times {\bf{p}}} {E_{{\bf{p}}}+m}\right) . \tag{ 84 }$$

Then, $\Lambda^{sr}_{{\bf{k}}-{\bf{q}}/2}$ in the distribution function expressed by Eq. (69) is given by

$$\Lambda^{sr}_{{\bf{k}}-{\bf{q}}/2}= \Lambda_{{\bf{k}}-{\bf{q}}/2}\delta^{sr} +\frac{1}{2} {\boldsymbol\Omega}_{{\bf{k}}-{\bf{q}}/2}^{ \rm{D}} \cdot{\boldsymbol\lambda}^{sr} . \tag{ 85 }$$

Using a similar expansion to that of Eq. (41) for fermions up to the first order, that is,

$$\frac{1}{{\rm e}^{X+Y} + 1} =f_{\rm F}(X)+f'_{\rm F}(X)Y - \frac{1}{2}f''_{\rm F}(X)C +\cdots \quad \tag{ 86 }$$

where

$$X=b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha,\ \ \ Y = -\Lambda_{{\bf{k}}-{\bf{q}}/2} , \tag{ 87 }$$

$$C=[X,Y]= {\rm i} \left[\left(b_0 {\boldsymbol\varepsilon}+ {\boldsymbol\omega}\times {\bf{b}}\right) \cdot \left( {\bf{k}}-\frac{\bf{q}}{2}\right) - ({\bf{b}}\cdot{\boldsymbol\varepsilon})E_{{\bf{k}}-{\bf{q}}/2}\right], \tag{ 88 }$$

we can expand the distribution function up to the first order. The zeroth-order result is just the Fermi-Dirac distribution,

$$f^{(0)}_{rs}(x,k) =f_{\rm F}(b\cdot p-\alpha)\delta^{sr}, \tag{ 89 }$$

and the first-order result is given by

$$\begin{aligned}[b] f^{(1)}_{rs}(x,k) =\; & f'_{\rm F}(b\cdot p-\alpha)\left[E_{\bf{k}}{\boldsymbol\varepsilon}\cdot{\bf{x}} - {\bf{k}}\cdot ({\boldsymbol\varepsilon}t + {\boldsymbol\omega}\times{\bf{x}})\right]\delta^{sr} \\ & -\frac{1}{2}f'_{\rm F}(b\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}}\cdot{\boldsymbol\lambda}^{sr}. \end{aligned} \tag{ 90 }$$

Similar to the distribution function for the scalar field, the first term in $f^{(1)}_{rs}(x,k)$ above can be absorbed into $f^{(0)}_{rs}(x,k)$ with $b^\mu$ replaced by $\beta^\mu$. After this rearrangement of the contribution, we obtain the final distribution function up to the first order:

$$f_{rs}(x,k) = f_{\rm F}(\beta\cdot p-\alpha)\delta^{sr} -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \cdot{\boldsymbol\lambda}^{sr} \tag{ 91 }$$

or in matrix form:

$$\begin{aligned}[b] f(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\cdot 1 -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \\ & \cdot\left({\bf{n}}_3 \sigma_3 +{\bf{n}}_2 \sigma_2 + {\bf{n}}_1 \sigma_1 \right). \end{aligned} \tag{ 92 }$$

Then, the local polarization along the fixed direction ${\bf{n}}_3$ is expressed as

$$P(x,k) \equiv\frac{f_{+,+}(x,k) - f_{-,-}(x,k)} {f_{+,+}(x,k) + f_{-,-}(x,k)} =-\frac{ f'_{\rm F}}{2f_{\rm F} }{\boldsymbol\Omega}_{\bf{k}}\cdot{\bf{n}}_3 , \tag{ 93 }$$

where we have suppressed the argument $\beta\cdot p-\alpha$ in the last expression for brevity.

Let us compare our results with that of Eq. (10). Three differences are evident. First, we have no $1/m$ term, which becomes singular when the mass approaches zero. Second, we have no term proportional to $({\boldsymbol\omega}\cdot {\bf{k}}) {\bf{k}}$ in the middle term of Eq. (10). Third, the contribution from the acceleration term ${{\boldsymbol\varepsilon}\times{\bf{k}}}$ is suppressed by $E_{\bf{k}}/(E_{\bf{k}}+m)$ relative to the term proportional to ${\boldsymbol\omega}$. These differences originate from the fact that the definition of the polarization is different from the conventional Wigner functions. In the appendix, we provide a specific relation between the distribution function with spin and the Wigner function.

B. Polarization along the momentum direction

For the helicity polarization, the spin quantization direction ${\bf{n}}_3$ is that along the particle's momentum ${\bf{p}}= p (\sin\theta \cos\phi,~ \sin\theta \sin\phi,~ \cos\theta)$, and we have $\vartheta=\theta, \varphi=\phi$

$${\bf{n}}_3 = \hat{\bf{p}} = {\bf{e}}_{ p},\ \ \ {\bf{n}}_2 = \frac{\hat{ \bf z}\times \hat{\bf{p}}}{|\hat{ \bf z}\times \hat{\bf{p}}|} = {\bf{e}}_{ \phi},\ \ \ {\bf{n}}_1 = {\bf{n}}_2\times \hat{\bf{p}}={\bf{e}}_{ \theta} \tag{ 94 }$$

with the helicity spinor expressed as

$$\xi^{s}=(- {\rm i})^{\frac{1-s}{2}}\left(\begin{array}{c} \cos\dfrac{s\theta+(1-s)\pi/2}{2}{\rm e}^{-{\rm i}\frac{\phi+(1-s)\pi/2}{2}} \\ \sin\dfrac{s\theta+(1-s)\pi/2}{2} {\rm e}^{{\rm i}\frac{\phi+(1-s)\pi/2}{2}} \end{array} \right) . \tag{ 95 }$$

In this case, the derivative in the last term of Eq. (65) does act on the two-component spinor $\xi^r$. It is easy to verify that

$$\pmb{\nabla}_p\xi^{r}= -\frac{{\rm i}{\bf{e}}_\phi}{2p\sin\theta} \sigma_3\xi^r + \frac{ {\bf{e}}_\theta}{2p} r \xi^{-r}. \tag{ 96 }$$

With this contribution, we obtain an extra term $\delta\Lambda^{sr}_{{\bf{p}}\parallel}$ compared with $\Lambda^{sr}_{\bf{p}}$ for the fixed spin direction:

$$\left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}}+ {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{s\dagger}\right] = \left(\Lambda^{sr}_{\bf{p}}+\delta\Lambda^{sr}_{{\bf{p}}\parallel} \right){a^r_{\bf{p}}}^\dagger, \tag{ 97 }$$

where $\Lambda^{sr}_{\bf{p}}$ is given by Eq. (83) with ${\bf{n}}_3, {\bf{n}}_2, {\bf{n}}_1$ designated in Eq. (94) and the extra term given by

$$\delta\Lambda^{sr}_{{\bf{p}}\parallel} \equiv \frac{E_{\bf{p}}}{2p}\pmb{\mathcal{E}}_{\bf{p}}\cdot \left(-{{\bf{e}}_\phi}\cot \theta s\delta^{rs} + {\bf{e}}_\phi \delta^{-r,s} - {\rm i}{\bf{e}}_\theta s\delta^{-r,s}\right) . \tag{ 98 }$$

Then, the distribution function with helicity index is expressed as

$$\begin{aligned}[b] f_{rs}(x,k)=\; & \int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+{\rm i}{\bf{q}}\cdot{\bf{x}}} \\ & \times\left[\frac{1}{{\rm e}^{b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha - \Lambda_{{\bf{k}}-{\bf{q}}/2} -\delta\Lambda_{{\bf{k}}-{\bf{q}}/2\parallel}}+1}\right]^{sr} \delta\left({\bf{q}}\right). \end{aligned} \tag{ 99 }$$

The extra term contributes to the first-order distribution function:

$$\begin{aligned}[b] f^{(1)}_{rs}(x,k) = \; & f'_{\rm F}(b\cdot k-\alpha)\left[E_{\bf{k}}{\boldsymbol\varepsilon}\cdot{\bf{x}} - {\bf{k}}\cdot ({\boldsymbol\varepsilon}t + {\boldsymbol\omega}\times{\bf{x}})\right]\delta^{sr} \\ & -\frac{1}{2}f'_{\rm F}(b\cdot k-\alpha) {\boldsymbol\Omega}_{\bf{k}}\cdot{\boldsymbol\lambda}^{sr}+\frac{1}{2}f'_{\rm F}(b\cdot k-\alpha)\frac{E_{\bf{k}}}{k}\pmb{\mathcal{E}}_{\bf{k}}\\ & \cdot \left({{\bf{e}}_\phi}\cot \theta s\delta^{rs} -{\bf{e}}_\phi \delta^{-r,s} + {\rm i} {\bf{e}}_\theta s\delta^{-r,s}\right), \end{aligned} \tag{ 100 }$$

where the last term is an additional contribution with respect to the fixed spin direction. After the rearrangement from $b^\mu$ to $\beta^\mu$, we obtain the final distribution function up to the first order:

$$\begin{aligned}[b] f_{rs}(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\delta^{sr} -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}}\cdot{\boldsymbol\lambda}^{sr}\\ & +\frac{1}{2}f'_{\rm F}(\beta\cdot k-\alpha)\frac{E_{\bf{k}}}{k}\pmb{\mathcal{E}}_{\bf{k}}\\ & \cdot \left({{\bf{e}}_\phi}\cot \theta s\delta^{rs} -{\bf{e}}_\phi \delta^{-r,s} + {\rm i} {\bf{e}}_\theta s\delta^{-r,s}\right) \end{aligned} \tag{ 101 }$$

or in matrix form:

$$\begin{aligned}[b] f(x,k) = \; & f_{\rm F}(\beta\cdot p-\alpha)\cdot 1 -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \\ & \cdot \left({\bf{n}}_3 \sigma_3 +{\bf{n}}_2 \sigma_2 + {\bf{n}}_1 \sigma_1 \right)\\ & +\frac{1}{2}f'_{\rm F}(\beta\cdot k-\alpha)\frac{E_{\bf{k}}}{k}\pmb{\mathcal{E}}_{\bf{k}} \cdot \left({{\bf{e}}_\phi}\cot \theta \sigma_3 -{\bf{e}}_\phi \sigma_1 + {\bf{e}}_\theta\sigma_2\right). \end{aligned} \tag{ 102 }$$

The helicity polarization with ${\bf{n}}_3={\bf{e}}_p$ is given by

$$P(x,k) =-\frac{ f'_{\rm F}}{2f_{\rm F} } {\boldsymbol\Omega}_{\bf{k}}\cdot{\bf{n}}_3 +\frac{f'_{\rm F}}{2f_{\rm F}} \frac{E_{\bf{k}}}{k}\cot\theta\, \pmb{\mathcal{E}}_{\bf{k}} \cdot {{\bf{e}}_\phi}. \tag{ 103 }$$

The last term will contribute to additional helicity polarization, which is absent in the formalism expressed in Eq. (10).

C. Polarization perpendicular to the momentum

It is also interesting to consider the polarization perpendicular to the momentum. We have two independent directions perpendicular to the momentum. One choice of the spin quantization ${\bf{n}}_3$ is along the direction ${\bf{e}}_\phi$. This choice can be fulfilled by $\vartheta=\pi/2, \varphi=\phi+\pi/2$, yielding

$${\bf{n}}_3= {\bf{e}}_\phi,\ \ \ {\bf{n}}_2 = \hat{ \bf k}\times {\bf{n}}_3 ,\ \ \ {\bf{n}}_1 = {\bf{n}}_2\times {\bf{n}}_3= - \hat{\bf{k}}. \tag{ 104 }$$

Substituting Eq. (104) into Eq. (65) and using the following relation for this specific case:

$$\pmb{\nabla}_p\xi^{r}= -\frac{{\rm i} {\bf{e}}_\phi}{2p\sin\theta} \sigma_3\xi^r , \tag{ 105 }$$

we obtain an extra term $\delta\Lambda^{sr}_{{\bf{p}}\phi}$ to be added to $\Lambda^{sr}_{\bf{p}}$ for the fixed spin direction given in Eq. (83):

$$\left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}}+ {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{s\dagger}\right] = \left(\Lambda^{sr}_{\bf{p}}+\delta\Lambda^{sr}_{{\bf{p}}\phi} \right){a^r_{\bf{p}}}^\dagger , \tag{ 106 }$$

where

$$\delta\Lambda^{sr}_{\bf{p}\phi} = \frac{E_{\bf{p}}}{2p} \pmb{\mathcal{E}}_{\bf{p}}\cdot \left(- {{\rm i} {\bf{e}}_\phi }\cot\theta s\delta^{r,-s} + {{\bf{e}}_\phi} \delta^{-r,s}\right) , \tag{ 107 }$$

Substituting this contribution into the distribution function

$$\begin{aligned}[b] f_{rs}(x,k) =\; & \int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t} {\rm e}^{{\rm i}{\bf{q}}\cdot{\bf{x}}}\\ & \times \left[\frac{1}{{\rm e}^{b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha - \Lambda_{\bf{p}} -\delta\Lambda_{{\bf{p}}\phi}}+1}\right]^{sr} \delta\left({\bf{q}}\right),\end{aligned} \tag{ 108 }$$

we obtain the final distribution function up to the first order:

$$\begin{aligned}[b] f_{rs}(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\delta^{sr} -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \cdot {\pmb{\lambda}}^{sr}\\ & +\frac{1}{2}f'_{\rm F}(\beta\cdot k-\alpha)\frac{E_{\bf{k}}}{k}\pmb{\mathcal{E}}_{\bf{k}} \cdot \left( {\rm i} {{\bf{e}}_\phi}\cot \theta s\delta^{r,-s} -{\bf{e}}_\phi \delta^{-r,s}\right) \end{aligned} \tag{ 109 }$$

or in matrix form:

$$\begin{aligned}[b] f(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\cdot 1 -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \\ & \cdot \left({\bf{n}}_3 \sigma_3 +{\bf{n}}_2 \sigma_2 + {\bf{n}}_1 \sigma_1 \right)\\ & +\frac{1}{2}f'_{\rm F}(\beta\cdot k-\alpha)\frac{E_{\bf{k}}}{k}\pmb{\mathcal{E}}_{\bf{k}} \cdot \left({{\bf{e}}_\phi}\cot \theta \sigma_3 -{\bf{e}}_\phi \sigma_1 \right). \end{aligned} \tag{ 110 }$$

The polarization can be expressed as

$$P(x,k) =-\frac{ f'_0}{2f_0 }{\boldsymbol\Omega}_{\bf{k}}\cdot{\bf{n}}_3. \tag{ 111 }$$

Note that the polarization receives no extra contribution except for the designation of ${\bf{n}}_3$ as ${\bf{e}}_\phi$.

Another choice of the independent transverse direction is along the direction ${\bf{e}}_\theta$, which can be obtained by $\vartheta =\theta+\pi/2, \varphi=\phi$. In such a case, we have

$${\bf{n}}_3= {\bf{e}}_\theta,\ \ \ {\bf{n}}_2 = {\bf{e}}_\phi ,\ \ \ {\bf{n}}_1 = -{\bf{e}}_p \tag{ 112 }$$

and the derivative on the spinor satisfies the same relation as the helicity polarization expressed by Eq. (96). Following the same procedure as in the cases ${\bf{n}}_3={\bf{e}}_p$ and ${\bf{n}}_3={\bf{e}}_\phi$ described above, we obtain the relation

$$\left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}}+ {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{s\dagger}\right] = \left(\Lambda^{sr}+\delta\Lambda^{sr}_{{\bf{p}}\theta} \right){a^r_{\bf{p}}}^\dagger , \tag{ 113 }$$

where

$$\delta\Lambda^{sr}_{{\bf{p}}\theta} \equiv \frac{E_{\bf{p}}}{2p} \pmb{\mathcal{E}}_{\bf{k}} \cdot \left( { {\bf{e}}_\phi}\cot\theta \delta^{r,-s} - {\bf{e}}_\phi s \delta^{r,s} -{ {\rm i}{\bf{e}}_\theta} s\delta^{-r,s} \right), \tag{ 114 }$$

which leads to the final first-order result for the distribution function:

$$\begin{aligned}[b] f_{sr}(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\delta^{sr} -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \cdot\pmb{\lambda}^{sr}\\ & +\frac{1}{2}f'_{\rm F}(\beta\cdot k-\alpha)\frac{E_{\bf{k}}}{k} \pmb{\mathcal{E}}_{\bf{k}} \\ & \times \left(- {{\bf{e}}_\phi}\cot \theta \delta^{r,-s} +{\bf{e}}_\phi s \delta^{r,s} + {\rm i} {\bf{e}}_\theta s\delta^{-r,s}\right), \end{aligned} \tag{ 115 }$$

or in matrix form:

$$\begin{aligned}[b] f(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\cdot 1 -\frac{1}{2}f'_{\rm F}(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \\ & \cdot \left({\bf{n}}_3 \sigma_3 +{\bf{n}}_2 \sigma_2 + {\bf{n}}_1 \sigma_1 \right)\\ & +\frac{1}{2}f'_{\rm F}(\beta\cdot k-\alpha)\frac{E_{\bf{k}}}{k} \pmb{\mathcal{E}}_{\bf{k}} \cdot \left(-{{\bf{e}}_\phi}\cot \theta \sigma_1 +{\bf{e}}_\phi \sigma_3 + {\bf{e}}_\theta\sigma_2\right). \end{aligned} \tag{ 116 }$$

The polarization can be expressed as

$$P =-\frac{ f'_{\rm F}}{2f_{\rm F} } {\boldsymbol\Omega}_{\bf{k}}\cdot{\bf{n}}_3 +\frac{f'_{\rm F}}{2f_{\rm F} }\frac{E_{\bf{k}}}{k} \pmb{\mathcal{E}}_{\bf{k}} \cdot {{\bf{e}}_\phi} . \tag{ 117 }$$

Note that the last term is an extra contribution to the polarization with respect to the fixed spin direction.

Next, we proceed with the charge vector field with Lagrangian density

$$\mathcal{L} = -\frac{1}{2} F^\dagger_{\mu\nu}F^{\mu\nu} + m^2 A^\dagger_\mu A^\mu , \tag{ 118 }$$

where the field tensor $F^{\mu\nu}$ is defined as

$$F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu . \tag{ 119 }$$

The Euler-Lagrangian equation leads to the Proca equation:

$$\partial_\mu \partial^\mu A_\nu + m^2 A_\nu = 0 ,\ \ \ \partial_\mu \partial^\mu A_\nu^\dagger + m^2 A_\nu^\dagger = 0 \tag{ 120 }$$

with the following constraint condition:

$$\partial^\nu A_\nu = 0. \tag{ 121 }$$

to remove the spin-0 contribution. The general solution can be expressed as a Fourier transform:

$$A_\mu (x) = \int \frac{{\rm d}^3{\bf{p}}}{(2\pi)^3}\frac{1}{\sqrt{2E_{\bf{p}}}} \sum\limits_{r=1}^3\left( a_{\bf{p}}^{r}\epsilon_\mu^{r}({\bf{p}}) {\rm e}^{-{\rm i}p\cdot x} + b_{\bf{p}}^{r\dagger}\eta_\mu^{r*}({\bf{p}}) {\rm e}^{{\rm i}p\cdot x}\right), \tag{ 122 }$$

with the following constraint from Eq. (121)

$$p^\mu \epsilon_\mu^{r}({\bf{p}}) = 0,\ \ \ p^\mu \eta_\mu^{r}({\bf{p}}) = 0 \tag{ 123 }$$

and the expression

$$\epsilon^{r\mu}({\bf{p}})=\left(\frac{{\bf{p}}\cdot{\bf{n}}_r}{m}, {\bf{n}}_r +\frac{{\bf{p}}\cdot{\bf{n}}_r}{m(E+m)}{\bf{p}}\right). \tag{ 124 }$$

Here, ${\bf{n}}_r$ ($r=1,2,3$) is a real orthogonal unit vector satisfying ${\bf{n}}_3 = {\bf{n}}_1\times {\bf{n}}_2$. For our choice of $\epsilon_\mu^{r}({\bf{p}})$, we have $\epsilon^{r*}_\mu({\bf{p}}) = \epsilon^{r}_\mu({\bf{p}})$. The antiparticle part $\eta_\mu^{r}({\bf{p}})$ can be obtained in the same manner as for $\epsilon_\mu^{r}({\bf{p}})$. As usual, we choose ${\bf{n}}_3$ as the spin quantization direction. The creation and annihilation operators obey the commutation rules:

$$\begin{aligned}[b] \left[ a^r_{\bf{p}}, a_{\bar{\bf{p}}}^{s\dagger}\right]=(2\pi)^3 \delta^3({\bf{p}}-\bar{\bf{p}})\delta^{rs}, \;\left[ b^r_{\bf{p}}, b_{\bar{\bf{p}}}^{s\dagger}\right]=(2\pi)^3 \delta^3({\bf{p}}-\bar{\bf{p}})\delta^{rs}. \end{aligned} \tag{ 125 }$$

The conserved charge current of the Proca theory is given by

$$j^{\mu} = -i\left( \partial^\mu A_\nu^\dagger A^\nu - A^\dagger_\nu \partial^\mu A^\nu\right), \tag{ 126 }$$

The canonical energy-momentum tensor of the Proca theory is

$$T^{\mu\nu} = \frac{1}{2}g^{\mu\nu}F^\dagger_{\alpha\beta}F^{\alpha\beta} - F^{\mu\alpha\dagger}\partial^\nu A_\alpha - \partial^\nu A^\dagger_\alpha F^{\mu\alpha} -m^2 g^{\mu\nu}A^\dagger_\alpha A^\alpha . \tag{ 127 }$$

The angular momentum density is then obtained:

$$\begin{aligned}[b] \mathcal{M}^{\mu\alpha\beta} =\; & x^\alpha T^{\mu\beta} - x^\beta T^{\mu\alpha} +\left( A^{\alpha\dagger} F^{\mu\beta} - A^{\beta\dagger} F^{\mu\alpha}\right)\\ & +\left( F^{\mu\beta\dagger} A^{\alpha} - F^{\mu\alpha\dagger} A^{\beta}\right). \end{aligned} \tag{ 128 }$$

The conserved charge and energy-momentum can be expressed as

$$\begin{aligned}[b] & Q = \int\frac{{\rm d}^3{\bf{p}}}{(2\pi)^3} \sum\limits_{s=1}^3 \left( a_{\bf{p}}^{s\dagger} a_{\bf{p}}^{s} -b_{\bf{p}}^{s} b_{\bf{p}}^{s\dagger} \right),\\ & P^\mu = \int\frac{{\rm d}^3{\bf{p}}}{(2\pi)^3} \sum\limits_{s=1}^3 p^{\mu}\left( a_{\bf{p}}^{s\dagger} a_{\bf{p}}^{s} +b_{\bf{p}}^{s} b_{\bf{p}}^{s\dagger} \right). \end{aligned} \tag{ 129 }$$

The angular momentum tensor can be expressed as

$$\begin{aligned}[b] K^{i} =\; & \int \frac{{\rm d}^3{\bf{p}}}{(2\pi)^3} \sum\limits_{s=1}^3\sum\limits_{r=1}^3\\ & \left\{ \sqrt{E_{\bf{p}}} \epsilon_\alpha^{s*}({\bf{p}})a_{\bf{p}}^{s\dagger} {\rm i} \partial_i^p \left[ \sqrt{E_{\bf{p}}} \epsilon^{r\alpha }({\bf{p}})a_{\bf{p}}^{r}\right] \right.\\ & - \sqrt{E_{\bf{p}}} \eta_\alpha^{s}({\bf{p}}) b_{\bf{p}}^{s} {\rm i} \partial_i^p \left[\sqrt{E_{\bf{p}}} \eta^{r \alpha *}({\bf{p}})b_{\bf{p}}^{r\dagger}\right]\\ & \left. -{\rm i} a_{\bf{p}}^{s\dagger}a_{\bf{p}}^{r} \left[ \epsilon^{s0*}({\bf{p}}) \epsilon^{r i }({\bf{p}}) - \epsilon^{s i*}({\bf{p}}) \epsilon^{r 0}({\bf{p}}) \right]\right.\\ & \left. + {\rm i} b_{\bf{p}}^{s} b_{\bf{p}}^{r\dagger} \left[ \eta^{s 0}({\bf{p}}) \eta^{r i *}({\bf{p}}) - \eta^{s i}({\bf{p}}) \eta^{r 0 *}({\bf{q}}) \right]\frac{}{}\right\}, \end{aligned} \tag{ 130 }$$

$$\begin{aligned}[b] J^{i} =\; & \epsilon^{0ijk} \int\frac{{\rm d}^3{\bf{p}}}{(2\pi)^3} \sum\limits_{s=1}^3\sum\limits_{r=1}^3\\ & \left\{\frac{}{} - a_{\bf{p}}^{s\dagger} \epsilon^{s \alpha *}({\bf{p}}) p^k {\rm i} \partial^p_j\left[ \epsilon^{r }_\alpha({\bf{p}}) a_{\bf{p}}^{r}\right] \right.\\ & - b_{\bf{p}}^{s} \eta^{s\alpha }({\bf{p}}) p^k {\rm i} \partial_j^p \left[\eta^{r*}_\alpha({\bf{p}})b_{\bf{p}}^{r\dagger}\right] \\ & \left.-{\rm i} a_{\bf{p}}^{s\dagger} a_{\bf{p}}^{r} \epsilon^{s j *}({\bf{p}}) \epsilon^{r k}({\bf{p}}) + {\rm i} b_{\bf{p}}^{s} b_{\bf{p}}^{r\dagger} \eta^{s j}({\bf{p}})\eta^{r k* }({\bf{p}}) \frac{}{}\right\} . \end{aligned} \tag{ 131 }$$

It is straightforward to obtain the commutation relations:

$$\begin{aligned}[b] & \left[Q, a_{\bf{p}}^{s\dagger}\right] = a_{{\bf{p}}}^{s\dagger},\ \ \ \left[b_\mu P^\mu, a_{\bf{p}}^{s\dagger}\right] = b_\mu p^\mu a_{{\bf{p}}}^{s\dagger},\ \ \\ & \left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}} + {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{s\dagger}\right] =\sum\limits_r \Lambda_{\bf{p}}^{sr}{a_{\bf{p}}}^{r\dagger} , \end{aligned} \tag{ 132 }$$

where the operator $\Lambda_{\bf{p}}^{sr}$ for the vector field is defined by

$$\begin{aligned}[b] \Lambda_{\bf{p}}^{sr}=\; & - \,\epsilon^{s}_\alpha({\bf{p}}) \epsilon^{r \alpha*}({\bf{p}}) {E_{\bf{p}}} \pmb{\mathcal{E}}_{\bf{p}} \cdot {\rm i}\,\pmb{\nabla}_{p} -\frac{\rm i}{2E_{\bf{p}}} \,\epsilon^{s}_\alpha({\bf{p}})\epsilon^{r \alpha*}({\bf{p}}) {\boldsymbol\varepsilon}\cdot {\bf{p}} \\ & - \,\epsilon^{s}_\alpha({\bf{p}})\left[{\rm i}\,\pmb{\nabla}_{p} \epsilon^{r \alpha*}({\bf{p}})a_{\bf{p}}^{r\dagger} \right] \cdot {E_{\bf{p}}} \pmb{\mathcal{E}}_{\bf{p}}\\ & - {\rm i}\, {\boldsymbol\varepsilon}\cdot \left[ \epsilon^{r0*}({\bf{p}}) {\boldsymbol \epsilon}^{s }({\bf{p}}) - \epsilon^{s 0}({\bf{p}}){ \boldsymbol\epsilon}^{r*}({\bf{p}}) \right]\\ & -{\rm i}\,{\boldsymbol\omega}\cdot \left[ {\boldsymbol\epsilon}^{r*}({\bf{p}})\times {\boldsymbol \epsilon}^{s }({\bf{p}}) \right]. \end{aligned} \tag{ 133 }$$

It follows that the distribution function is given by

$$\begin{aligned}[b]f_{rs}(x,k)\equiv\; & \frac{1}{(2\pi)^3}\int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+{\rm i} {\bf{q}}\cdot{\bf{x}}} \langle a^{s\dagger}_{{\bf{k}}-{\bf{q}}/{2}} a^r_{ {\bf{k}}+{\bf{q}}/{2}}\rangle_0\\ =\; & \int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+ {\rm i} {\bf{q}}\cdot{\bf{x}}}\\ & \times \left[\frac{1}{{\rm e}^{b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha - \Lambda_{{\bf{k}}-{\bf{q}}/2}}-1}\right]^{sr} \delta\left({\bf{q}}\right). \end{aligned} \tag{ 134 }$$

In the following, we further analyze this distribution function by using the following identities:

$$\begin{aligned}[b] \epsilon_\mu^{r}({\bf{p}})\epsilon^{ s \mu*}({\bf{p}}) =\; & -\delta^{r s} , \\ \epsilon^{s }_\alpha({\bf{p}})\partial^i_p \left[ \epsilon^{r\alpha *}({\bf{p}})\right] =\; & - \frac{\epsilon^{ri*}({\bf{p}})}{E_{\bf{p}}+m}\epsilon^{s 0 }({\bf{p}}) \\ & +\frac{\epsilon^{r0*}({\bf{p}})}{E_{\bf{p}}+m}\epsilon^{s i }({\bf{p}}) {-{\bf{n}}_s\cdot \partial^i_p {\bf{n}}^*_r }. \end{aligned} \tag{ 135 }$$

A. Polarization along a fixed direction

If the unit vectors ${\bf{n}}_1$, ${\bf{n}}_2$, ${\bf{n}}_3$ are independent on the momentum, the last term in Eq. (135) will vanish and we obtain

$$\Lambda^{sr}_{\bf{p}} = \Lambda_{\bf{p}}\delta^{sr} + {\rm i}\, {\boldsymbol\Omega}_{\bf{p}} \cdot \left( {\bf{n}}_s \times {\bf{n}}_r\right). \tag{ 136 }$$

Following the same procedure applied for the scalar and Dirac fields, and using the expansion in Eq. (41), we obtain the zeroth-order result for the vector particle:

$$f^{(0)}_{rs}(x,k)=f_{\rm B}(b\cdot k-\alpha)\delta^{sr} . \tag{ 137 }$$

The first-order result is given by

$$\begin{aligned}[b] f^{(1)}_{rs}(x,k) =\; & f'_B(b\cdot k-\alpha)\left[E_{\bf{k}}{\boldsymbol\varepsilon}\cdot{\bf{x}} - {\bf{k}}\cdot ({\boldsymbol\varepsilon}t + \pmb{\omega}\times{\bf{x}})\right] \delta^{sr}\\ & - {\rm i} f'_B(b\cdot k-\alpha) {\boldsymbol\Omega}_{{\bf{k}}} \cdot \left( {\bf{n}}_s \times {\bf{n}}_r\right) .\end{aligned} \tag{ 138 }$$

As we all know, the vector $\epsilon_\mu^{r}({\bf{p}})$ with $r=1,2,3$ denotes the linear polarization vector and does not correspond to the spin eigenstate, which can be achieved by introducing circular polarization operators:

$$a^+_{\bf{p}} =- \frac{1}{\sqrt{2}}\left(a^{1}_{\bf{p}}-ia^{2}_{\bf{p}} \right),\ \ a^-_{\bf{p}} = \frac{1}{\sqrt{2}}\left(a^{1}_{\bf{p}}+ia^{2}_{\bf{p}} \right),\ \ a^{0}_{\bf{p}} = a^{3}_{\bf{p}}, \tag{ 139 }$$

where $0$ and $\pm$ denote the spin components with $0$, $\pm 1$, respectively. With circular polarization indices, the vector field can be expressed as

$$A_\mu (x) = \int \frac{{\rm d}^3{\bf{p}}}{(2\pi)^3}\frac{1}{\sqrt{2E_{\bf{p}}}} \sum\limits_{r=0,\pm}\left( a_{\bf{p}}^{r}\epsilon_\mu^{r}({\bf{p}}) {\rm e}^{-{\rm i}p\cdot x} + b_{\bf{p}}^{r\dagger}\eta_\mu^{r*}({\bf{p}}) {\rm e}^{{\rm i}p\cdot x}\right), \tag{ 140 }$$

where the polarization four-vectors with circular polarization index is defined by

$$\begin{aligned}[b] & \epsilon_\mu^+({\bf{p}})=-\frac{1}{\sqrt{2}}\left[\epsilon_\mu^{1}({\bf{p}})+ {\rm i}\epsilon_\mu^{2}({\bf{p}}) \right], \\ & \epsilon_\mu^-({\bf{p}})=\frac{1}{\sqrt{2}}\left[\epsilon_\mu^{1}({\bf{p}})- {\rm i}\epsilon_\mu^{2}({\bf{p}}) \right],\\ & \epsilon_\mu^{0}({\bf{p}})=\epsilon_\mu^{3}({\bf{p}}) . \end{aligned} \tag{ 141 }$$

We use the same indices $r,s$ to denote linear or circular polarization. With circular polarization indices $0,\pm$, we have

$$\begin{aligned}[b] f^{(1)}_{rs}(x,p) =\; & f'_{\rm B}(b\cdot k-\alpha)\left[E_{\bf{k}}{\boldsymbol\varepsilon}\cdot{\bf{x}} - {\bf{k}}\cdot ({\boldsymbol\varepsilon}t + \pmb{\omega}\times{\bf{x}})\right] \delta^{sr}\\ & - {\rm i} f'_{\rm B}(b\cdot k-\alpha) {\boldsymbol\Omega}_{{\bf{k}}} \cdot \left( {\bf{n}}_s \times {\bf{n}}_r^*\right), \end{aligned} \tag{ 142 }$$

where the polarization three-vector with circular indices is given by

$${\bf{n}}_+=-\frac{1}{\sqrt{2}}\left({\bf{n}}_1 + {\rm i} {\bf{n}}_2 \right), \ \ {\bf{n}}_-=\frac{1}{\sqrt{2}}\left({\bf{n}}_1 - i {\bf{n}}_2 \right),\ \ {\bf{n}}_0={\bf{n}}_3 . \tag{ 143 }$$

To obtain the non-trivial contribution for the spin alignment, we need the second-order result for the diagonal components with linear polarization indices:

$$\begin{aligned}[b] f_{rr}^{(2)}(x,k) =\; & f^{(2)}(x,k)+\frac{1}{2}f''_{\rm B}(b\cdot k-\alpha)\left( {\boldsymbol\Omega}_{{\bf{k}}}\right) ^2 \\ & -\frac{1}{2}f''_{\rm B}(b\cdot k-\alpha)\left( {\boldsymbol\Omega}_{{\bf{k}}} \cdot {\bf{n}}_r\right)^2 ,\end{aligned} \tag{ 144 }$$

where $f^{(2)}$ denotes the second-order contribution to the distribution function for the scalar field as expressed in Eq. (46). Summing $f^{(0)}_{rs}$ and $f^{(1)}_{rs}$ for non-diagonal components and replacing $b^\mu$ with $\beta^\mu$, we obtain the spin distribution function up to the first order:

$$f_{rs}(x,k) = f_{\rm B}(\beta\cdot k-\alpha)\delta^{sr}- i f'_{\rm B}(\beta\cdot k-\alpha) {\boldsymbol\Omega}_{{\bf{k}}} \cdot \left( {\bf{n}}_s \times {\bf{n}}_r^*\right). \tag{ 145 }$$

This expression is valid for both linear polarization indices $r,s=1,2,3$ and circular polarization indices $r,s=0,\pm$. For the linear polarization $r,s=1,2,3$, we have ${\bf{n}}_r^*={\bf{n}}_r$. Similar to the result for the Dirac fermion expressed in Eq. (92), we can express the distribution function with circular polarization indices in the conventional matrix form:

$$\begin{aligned}[b] f(x,k) =\; & f_{\rm F}(\beta\cdot p-\alpha)\cdot 1 - f'_F(\beta\cdot p-\alpha) {\boldsymbol\Omega}_{\bf{k}} \\ & \cdot\left({\bf{n}}_3 S_3 +{\bf{n}}_2 S_2 + {\bf{n}}_1 S_1 \right), \end{aligned} \tag{ 146 }$$

where $S_i$ denotes the $3\times 3$ spin matrices for spin-1

$$\begin{aligned}[b] & S_1=\frac{1}{\sqrt{2}}\left( \begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \\ \end{array} \right),\\ & S_2=\frac{\rm i}{\sqrt{2}}\left( \begin{array}{ccc} 0 & -1 & 0 \\ 1 & 0 & -1 \\ 0 & 1 & 0 \\ \end{array} \right),\\ & S_3=\left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \\ \end{array} \right) \end{aligned} \tag{ 147 }$$

Summing $f^{(0)}_{rr}$, $f^{(1)}_{rr}$, and $f^{(2)}_{rr}$ for the diagonal components with linear polarization indices and replacing $b^\mu$ with $\beta^\mu$, we obtain the diagonal distribution function up to the second order:

$$\begin{aligned}[b] f_{rr}(x,k) =\; & f(x,k) +\frac{1}{2}f''_{\rm B}(\beta\cdot k-\alpha) \left( {\boldsymbol\Omega}_{{\bf{k}}} \right)^2 \\ & -\frac{1}{2}f''_{\rm B}(\beta\cdot k-\alpha)\left({\boldsymbol\Omega}_{{\bf{k}}} \cdot {\bf{n}}_r\right)^2 , \end{aligned} \tag{ 148 }$$

where $f(x,k)$ is the distribution function for scalar field as expressed in Eq. (48). It should be noted that this second-order expression holds only for linear polarization indices $r,s=1,2,3$.

As we all know, the spin density matrix for vector particles such as ϕ and $K^{*0}$ mesons can be measured by their two-body decay channels $\phi\rightarrow KK$ and $K^{*0}\rightarrow K\pi$, in which the distribution of the decay products is related to the elements of spin density matrix by

$$\begin{aligned}[b]\frac{{\rm d} N}{{\rm d}\cos\theta^* {\rm d}\phi^*} =\; & \frac{3}{8\pi} \Big[ 1 - \rho_{00}+(3\rho_{00}-1)\cos^2\theta^* \\ & -\sqrt{2} \rm{Re}\left( \rho_{+0} - \rho_{0-} \right)\sin\theta^*\cos\phi^*\\ & +\sqrt{2} \rm{Im} \left(\rho_{+0} - \rho_{0-} \right)\sin\theta^*\sin\phi^*\\ & -2 \rm{Re} \rho_{+-} \sin^2\theta^*\cos(2\phi^*) \\ & -2 \rm{Im} \rho_{+-} \sin^2\theta^*\sin(2\phi^*)\Big] , \end{aligned} \tag{ 149 }$$

where $\theta^*$ and $\phi^*$ are the polar and azimuthal angles of the momentum of one final meson in the rest frame of the initial vector mesons. When we approximate the spin distribution function up to the first order as given in Eq. (145), all the non-diagonal elements vanish and the diagonal element $\rho_{00}$ becomes $1/3$, which means that there is no spin alignment. At the second order, the spin alignment receives a nonzero contribution because the local spin alignment can be expressed as

$$\rho_{00}(x,k) = \frac{f_{00}}{f_{++}+f_{-}+f_{00}} =\frac{1}{3} -\frac{1}{6}\frac{f_{\rm B}''}{f_{\rm B}} \left[ \left( {\boldsymbol\Omega}_{{\bf{k}}} \cdot {\bf{n}}_3\right)^2 -\frac{1}{3}\left({\boldsymbol\Omega}_{{\bf{k}}}\right)^2 \right]. \tag{ 150 }$$

Note that whether the spin alignment is less or greater than $1/3$ depends on the balance between the contribution $\left({\boldsymbol\Omega}_{{\bf{k}}} \cdot {\bf{n}}_3\right)^2$ and $\left({\boldsymbol\Omega}_{{\bf{k}}}\right)^2/3$.

B. Polarization along the momentum direction

Next, let us consider the helicity polarization with $\hat{\bf{p}}$ as the spin quantization direction:

$${\bf{n}}_3 = {\bf{e}}_p,\ \ {\bf{n}}_2={\bf{e}}_\phi, \ \ {\bf{n}}_1={\bf{e}}_\theta . \tag{ 151 }$$

To calculate the last term in the second identity in Eq. (135), we need the following relations:

$$\frac{\partial}{\partial p}{\bf{n}}_s=0,\ \ \ \frac{\partial}{\partial\phi}{\bf{n}}_s=\hat{\bf{z}}\times {\bf{n}}_s,\ \ \ \frac{\partial}{\partial\theta}{\bf{n}}_s = {\bf{e}}_\phi\times {\bf{n}}_s . \tag{ 152 }$$

With these relations, some additional terms contribute to the commutation relation:

$$\left[{{\boldsymbol\varepsilon}}\cdot{\bf{K}} + {{\boldsymbol\omega}}\cdot{\bf{J}}, a_{\bf{p}}^{s\dagger}\right] \equiv\left(\Lambda^{sr}_{\bf{p}}{ +\delta \Lambda^{sr}_{{\bf{p}}}} \right){a^r_{\bf{p}}}^\dagger , \tag{ 153 }$$

where $\Lambda^{sr}_{\bf{p}}$ is given by Eq. (136) and the additional term $\delta\Lambda^{sr}_{{\bf{p}}\parallel}$ is given by

$$\delta\Lambda^{sr}_{{\bf{p}}} =- \frac{{\rm i} E_{\bf{p}} }{p}\left( \pmb{\mathcal{E}}_{\bf{p}}\cdot {\bf{e}}_\theta\right) {\bf{e}}_\phi \cdot\left({\bf{n}}_s\times {\bf{n}}_r\right) - \frac{{\rm i} E_{\bf{p}} }{p\sin\theta} \left( \pmb{\mathcal{E}}_{\bf{p}}\cdot {\bf{e}}_\phi\right) \hat{\bf{z}} \cdot \left({\bf{n}}_s\times {\bf{n}}_r\right). \tag{ 154 }$$

From the definition

$$\begin{aligned}[b] f_{rs}(x,k) =\; & \int {\rm d}^3 {\bf{q}} {\rm e}^{-{\rm i}(E_{{\bf{k}}+{\bf{q}}/2}-E_{{\bf{k}}-{\bf{q}}/2})t+ {\rm i}{\bf{q}}\cdot{\bf{x}}} \\ & \times\left[\frac{1}{{\rm e}^{b_0 E_{{\bf{k}}-{\bf{q}}/2} - {\bf{b}}\cdot({\bf{k}}-{\bf{q}}/2) - \alpha - \Lambda_{{\bf{k}}-{\bf{q}}/2} -\delta\Lambda_{{\bf{k}}-{\bf{q}}/2,\parallel}}-1}\right]^{sr} \delta\left({\bf{q}}\right), \end{aligned} \tag{ 155 }$$

we obtain the distribution function up to the first order:

$$\begin{aligned}[b] f_{rs} =\; & f_{\rm B} - i f'_{\rm B} \left[{\boldsymbol\Omega}_{\bf{k}} -\frac{E_{\bf{k}} }{k}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\theta\right) {\bf{e}}_\phi -\frac{E_{\bf{k}}}{k\sin\theta}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\phi \right) \hat{\bf{z}}\right]\\ & \cdot \left( {\bf{n}}_s \times {\bf{n}}_r^*\right).\end{aligned} \tag{ 156 }$$

from which we find that all the non-diagonal elements of the spin density matrix vanish and there is no spin alignment. To obtain nonzero spin alignment, we need the second-order result for the diagonal elements:

$$\begin{aligned}[b] f_{rr} =\; & f_{\rm B}+\frac{1}{2}f''_{\rm B} \left[{\boldsymbol\Omega}_{\bf{k}} -\frac{E_{\bf{k}} }{k}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\theta\right) {\bf{e}}_\phi -\frac{E_{\bf{k}}}{k\sin\theta}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\phi \right) \hat{\bf{z}}\right]^2\\ & -\frac{1}{2}f''_{\rm B}\left\{\left[{\boldsymbol\Omega}_{\bf{k}} -\frac{E_{\bf{k}} }{k}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\theta\right) {\bf{e}}_\phi -\frac{E_{\bf{k}}}{k\sin\theta}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\phi \right) \hat{\bf{z}}\right]\cdot {\bf{n}}_r\right\}^2 . \end{aligned} \tag{ 157 }$$

It follows that the spin alignment is given by

$$\begin{aligned}[b] \rho_{00} = \; & \frac{1}{3} - \frac{1}{6} \cdot \frac{f_{\rm B}''}{f_{\rm B}} \left\{ \left[ \left( {\boldsymbol\Omega}_{\bf{k}} -\frac{E_{\bf{k}} }{k}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\theta\right) {\bf{e}}_\phi -\frac{E_{\bf{k}}}{k\sin\theta}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\phi \right) \hat{\bf{z}} \right) \cdot {\bf{n}}_3 \right]^2\right.\\ & \left. -\frac{1}{3}\left[{\boldsymbol\Omega}_{\bf{k}} -\frac{E_{\bf{k}} }{k}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\theta\right) {\bf{e}}_\phi -\frac{E_{\bf{k}}}{k\sin\theta}\left( \pmb{\mathcal{E}}_{\bf{k}} \cdot {\bf{e}}_\phi \right) \hat{\bf{z}}\right]^2\right\}. \end{aligned} \tag{ 158 }$$

C. Polarization perpendicular to the momentum

Finally, let us consider the transverse polarization, which is orthogonal to the momentum. Similar to the Dirac particle, we have two independent basis vectors. We can choose one group of basis vectors as

$${\bf{n}}_3 ={\bf{e}}_\phi ,\ \ {\bf{n}}_2={\bf{e}}_\theta, \ \ {\bf{n}}_1={\bf{e}}_p \tag{ 159 }$$

and the other group as

$${\bf{n}}_3 = {\bf{e}}_\theta,\ \ {\bf{n}}_2={\bf{e}}_p, \ \ {\bf{n}}_1={\bf{e}}_\phi . \tag{ 160 }$$

For both groups, we follow the same procedure used for the helicity polarization and find that all the results coincide with those expressed in subsection IV.B by Eqs. (152) to (158). The only difference is that we need to replace ${\bf{n}}_1, {\bf{n}}_2, {\bf{n}}_3$ in Eq. (151) by either Eq. (159) or Eq. (160).

We have revisited the spin polarization by thermal vorticity and proposed another formalism to calculate them directly from the spin-dependent distribution function. We have calculated these spin-dependent distribution functions for spin-1/2 and spin-1 in global equilibrium with thermal vorticity. For the Dirac field with spin 1/2, the local polarization along the fixed direction ${\bf{n}}_3$ is given by

$$P(x,k)=-\frac{ f'_{\rm F}}{2f_{\rm F} }\left({\boldsymbol\omega}-\frac{{\boldsymbol\varepsilon}\times {\bf{k}}}{E_{\bf{k}}+m}\right)\cdot{\bf{n}}_3 . \tag{ 161 }$$

For a vector field with spin 1, the spin alignment and all the diagonal elements of the spin density matrix are absent up to the first order. At the second order, the local spin alignment along the fixed direction ${\bf{n}}_3$ is given by

$$\begin{aligned}[b] \rho_{00}(x,k) =\; & \frac{1}{3} -\frac{1}{6}\frac{f_{\rm B}''}{f_{\rm B}} \Bigg\{ \left[\left({\boldsymbol\omega}-\frac{{\boldsymbol\varepsilon}\times {\bf{k}}}{E_{\bf{k}}+m}\right) \cdot {\bf{n}}_3\right]^2 \\ & -\frac{1}{3}\left({\boldsymbol\omega}-\frac{{\boldsymbol\varepsilon}\times {\bf{k}}}{E_{\bf{k}}+m}\right)^2 \Bigg\}. \end{aligned} \tag{ 162 }$$

Besides, we have found that when the spin quantization direction is dependent on the momentum, an additional contribution emerges, departing from an earlier prediction. It will be valuable to conduct numerical simulations with these results and quantitatively study the differences in the future.

To finish this study, we would like to discuss on the source of these differences. When we predict the spin polarization in heavy-ion collisions by the relativistic hydrodynamics, we need the local spin distribution function in phase space, which depends on both the coordinates and momentum. However, as we all know, such distribution functions in quantum mechanics are not unique. Different distribution functions might lead to different predictions for some physical results. We do not know which distribution function is more appropriate to describe the realistic physics in heavy-ion collisions. More research on this topic is needed in the future.

In this appendix, we relate the distribution functions with spin used in this work to the Wigner functions. We take the Dirac field as an example. From the free Dirac field given in Eq. (55) with the following normalization relations:

$${ u^{s\dagger}({\bf{p}})} u^r({\bf{p}})={ v^{s\dagger}({\bf{p}})} v^r({\bf{p}})=2E_p \delta^{sr}, \tag{A1}$$

$${u^{s\dagger}({\bf{p}})} v^r({-\bf p}) ={ v^{s\dagger}({\bf{p}})} u^r({\bf{-p}}) =0, \tag{A2}$$

we have

$$a^s_{\bf{p}} =\frac{1 }{\sqrt{2E_{p}}}\int {\rm d}^3{\bf{x}} {\rm e}^{{\rm i}p\cdot x} u^{s\dagger}({\bf{p}}) \psi(x), \tag{A3}$$

$$a^{s\dagger}_{\bf{p}} = \frac{1 }{\sqrt{2E_{p}}}\int {\rm d}^3{\bf{x}} {\rm e}^{-{\rm i}p\cdot x} \psi^\dagger(x)u^{s}({\bf{p}}) . \tag{A4}$$

Substituting these expressions into the spin distribution function in Eq. (69), we obtain

$$\begin{aligned}[b] f_{rs}(x,k) =\; & \frac{1}{(2\pi)^3}\int{{\rm d}^3 {\bf{q}}} {\rm e}^{-{\rm i}( E_{{\bf{k}}+{\bf{q}}/2} -E_{{\bf{k}}-{\bf{q}}/2}) t} {\rm e}^{{\rm i}{\bf{q}}\cdot{\bf{x}}} \langle a^{s\dagger}_{\bf{k}-\bf{q}/{2}} a^r_{ \bf{k}+\bf{q}/{2}}\rangle\\ =\; & \frac{1}{(2\pi)^3}\int{{\rm d}^3 {\bf{q}}} {\rm d}^3 {\bf{y}}_1 {\rm d}^3 {\bf{y}}_2 {\rm e}^{{\rm i}{\bf{q}}\cdot{\bf{x}}} {\rm e}^{{\rm i}({\bf{k}}-{\bf{q}}/2)\cdot {\bf{y}}_1} {\rm e}^{-{\rm i}({\bf{k}}+{\bf{q}}/2)\cdot {\bf{y}}_2}\\ & \times \frac{1 }{\sqrt{2E_{{\bf{k}}+{\bf{q}}/2}}}\frac{1 }{\sqrt{2E_{{\bf{k}}-{\bf{q}}/2}}} \langle \psi^\dagger(t,{\bf{y}}_1)u^{s}({\bf{k}}-{\bf{q}}/2)\\ & \times u^{r\dagger}({\bf{k}}+{\bf{q}}/2) \psi(t,{\bf{y}}_2) \rangle . \end{aligned}\tag{A5}$$

With the definitions ${\bf{y}}=({\bf{y}}_1+{\bf{y}}_2)/2$ and ${\bf{z}}={\bf{y}}_1-{\bf{y}}_2$, the distribution function becomes

$$\begin{aligned}[b] f_{rs}(x,k) =\; & \frac{1}{(2\pi)^3}\int{{\rm d}^3 {\bf{q}}} {\rm e}^{{\rm i}{\bf{q}}\cdot{\bf{x}}} {\rm d}^3 {\bf{y}} {\rm d}^3 {\bf{z}} {\rm e}^{{\rm i} {\bf{k}}\cdot {\bf{z}}} e^{-i {\bf{q}}\cdot {\bf{y}}}\\ & \times \frac{1 }{\sqrt{2E_{{\bf{k}}+{\bf{q}}/2}}}\frac{1 }{\sqrt{2E_{{\bf{k}}-{\bf{q}}/2}}}\\ & \times \rm{Tr}\Big[ \langle \psi^\dagger(t,{\bf{y}}+ \frac{{\bf{z}}}{2}) \psi(t,{\bf{y}}- \frac{{\bf{z}}}{2}) \rangle \\ & \times u^{s}({\bf{k}}-{\bf{q}}/2)u^{r\dagger}({\bf{k}}+{\bf{q}}/2)\Big], \end{aligned} \tag{A6}$$

where we have expressed the spinor matrix into the trace form. Recall that the equal-time Wigner function is defined as

$$W_{ab}(t,{\bf{y}},{\bf{k}}) =\int \frac{{\rm d}^3 {\bf{z}}}{(2\pi)^3} {\rm e}^{{\rm i} {\bf{k}}\cdot {\bf{z}}} \langle \bar\psi_b(t,{\bf{y}}+ \frac{ {\bf{z}}}{2}) \psi_a(t,{\bf{y}}- \frac{{\bf{z}}}{2}) \rangle. \tag{A7}$$

We can relate the distribution functions with spin to this Wigner function by

$$\begin{aligned}[b] f_{rs}(x,k) =\; & \int{{\rm d}^3 {\bf{q}}}\int {\rm d}^3 {\bf{y}} {\rm e}^{{\rm i}{\bf{q}}\cdot{\bf{x}}} {\rm e}^{-{\rm i} {\bf{q}}\cdot {\bf{y}}} \frac{1 }{\sqrt{2E_{{\bf{k}}+{\bf{q}}/2}}}\frac{1 }{\sqrt{2E_{{\bf{k}}-{\bf{q}}/2}}}\\ & \times \rm{Tr}\left[W(t,{\bf{y}},{\bf{k}}) \gamma^0 u^{s}({\bf{k}}-{\bf{q}}/2)u^{r\dagger}({\bf{k}}+{\bf{q}}/2)\right]. \end{aligned}\tag{A8}$$

It is evident that the distribution function with spin and the usual Wigner function are related to each other in a nontrivial manner.

It is convenient to compare this result with the usual definition given in [48], which leads to the relation between the Wigner and distribution functions:

$$f_{rs}(x, {\bf{k}}) =\frac{(2\pi)^3}{2E_{\bf{k}}} {\rm Tr}\left[ W(t,{\bf{x}}, {\bf{k}})\gamma^0 u^s({\bf{k}}) u^{r\dagger}({\bf{k}})\right]. \tag{A9}$$

Evidently, these two definitions only coincide with each other after integration over the whole coordinate space ${\rm d}^3 {\bf{x}}$. It is easy to show that this definition can lead to conventional results. To extract the polarization, we only need to consider the diagonal elements, which can be recast into

$$f_{ss}(x, {\bf{k}}) =\frac{(2\pi)^3}{2E_{\bf{k}}} \rm{Tr}\left[\gamma^0 W(t,{\bf{x}}, {\bf{k}})\gamma^0 u^s({\bf{k}}) \bar u^{s}({\bf{k}})\right]. \tag{A10}$$

Using the identity in [50] for the spinor $u^s({\bf{k}})$

$$u^s({\bf{k}}) \bar u^s ({\bf{k}}) =\frac{1}{2}\left[ (1+\gamma^5 {\not s})( {\not k} +m )\right] \tag{A11}$$

with the polarization vector $s^\mu$ defined by

$$s^\mu \equiv \left( \frac{{\bf{k}}\cdot {\bf{n}}_3}{m}, {\bf{n}}_3+ \frac{({\bf{k}}\cdot {\bf{n}}){\bf{k}}}{m(m+E_{\bf{k}})}\right), \tag{A12}$$

and the Wigner function's decomposition in Clifford algebra

$$W(t,{\bf{x}}, {\bf{k}})= \frac{1}{4}\left[\mathcal{F}+ {\rm i}\gamma^5 \mathcal{P} +\gamma^\mu \mathcal{V}_\mu +\gamma^5 \gamma^\mu \mathcal{A}_\mu +\frac{1}{2}\sigma^{\mu\nu} \mathcal{S}_{\mu\nu}\right] , \tag{A13}$$

we obtain

$$\begin{aligned}[b] f_{ss}(x, {\bf{k}}) =\; & \frac{(2\pi)^3}{8E_{\bf{k}}}\left(\frac{}{} m \rm{Tr}\left[ 1 \right] \mathcal{F} + \rm{Tr}\left[\gamma^0 \gamma^{\mu}\gamma^ 0 {\not k} \right] \mathcal{V}_{\mu} \right.\\ & \left. + m \rm{Tr}\left[ \gamma^0 \gamma^\mu \gamma^0 {\not s} \right] \mathcal{A}_\mu +\frac{1}{2} \rm{Tr}\left[ \gamma^5\gamma^0 \sigma^{\mu\nu}\gamma^0 {\not s} {\not k} \right] \mathcal{S}_{\mu\nu} \right)\\ =\; & \frac{(2\pi)^3}{2E_{\bf{k}}}\left(\frac{}{} m \mathcal{F} - k^\mu \mathcal{V}_{\mu} + 2 E_{\bf{k}} \mathcal{V}_0 - m s^\mu \mathcal{A}_\mu + 2m s^0 \mathcal{A}_0\right.\\ & \left. -\frac{1}{2}\epsilon^{\mu\nu\alpha\beta} s_\alpha k_\beta \mathcal{S}_{\mu\nu} + s^0 \epsilon^{0\mu\nu\alpha} k_\alpha \mathcal{S}_{\mu\nu} - E_{\bf{k}} \epsilon^{0\mu\nu\alpha} s_\alpha \mathcal{S}_{\mu\nu} \right). \end{aligned}\tag{A14}$$

Using the Wigner equations [51, 52] up to the semiclassical limit, we finally obtain

$$f_{ss}(x, {\bf{k}}) =(2\pi)^3\frac{E_{\bf{k}}}{m}\left(\frac{}{} \mathcal{F} - s^\mu \mathcal{A}_\mu\right). \tag{A15}$$

The conventional result for the polarization is as follows:

$$P(x,{\bf{k}}) = - \frac{ s^\mu \mathcal{A}_\mu }{ \mathcal{F} }. \tag{A16}$$

It should be noted that only positive energy parts are included in the Wigner functions.

We thank F. Becattini and S. Pu for helpful discussions.