††thanks: These authors contribute equally to this work.††thanks: These
authors contribute equally to this work.

# Abnormal Phonon Angular Momentum due to Off-diagonal Elements in Density
Matrix induced by Temperature Gradient

Jinxin Zhong  School of Physics Science and Engineering, Tongji University,
Shanghai 200092, China  Hong Sun  Phonon Engineering Research Center of
Jiangsu Province, Center for Quantum Transport and Thermal Energy Science,
Institute of Physics and Interdisciplinary Science, School of Physics and
Technology, Nanjing Normal University, Nanjing 210023, China  Yang Pan  Phonon
Engineering Research Center of Jiangsu Province, Center for Quantum Transport
and Thermal Energy Science, Institute of Physics and Interdisciplinary
Science, School of Physics and Technology, Nanjing Normal University, Nanjing
210023, China  Zhiguo Wang  School of Physics Science and Engineering, Tongji
University, Shanghai 200092, China  Xiangfan Xu  School of Physics Science and
Engineering, Tongji University, Shanghai 200092, China  Lifa Zhang
<EMAIL_ADDRESS>Phonon Engineering Research Center of Jiangsu Province,
Center for Quantum Transport and Thermal Energy Science, Institute of Physics
and Interdisciplinary Science, School of Physics and Technology, Nanjing
Normal University, Nanjing 210023, China  Jun Zhou<EMAIL_ADDRESS>Phonon Engineering Research Center of Jiangsu Province, Center for Quantum
Transport and Thermal Energy Science, Institute of Physics and
Interdisciplinary Science, School of Physics and Technology, Nanjing Normal
University, Nanjing 210023, China

###### Abstract

Nonzero mean value of phonon angular momentum (PAM) in chiral materials can be
generated when a temperature gradient is applied. We find that both diagonal
and off-diagonal terms of PAM contribute to mean PAM by using the Kubo formula
where both diagonal and off-diagonal elements of the heat current operator are
considered. The calculation results show that the off-diagonal term is
dominant when the phonon scattering is strong enough. This finding reveals
that the quantum transition between different phonon modes induced by
temperature gradient strongly affects the local atomic rotation. Our discovery
provides an explanation of the recently observed chiral phonon activated spin
Seebeck effect.

Phonons are the quanta of atomic displacement field in solids Ashcroft and
Mermin (1976). For a given monatomic lattice, the displacement of each atom
can be obtained through solving the equations of motion once the dynamical
matrix ${\bf D}({\bf k})$ is known where ${\bf k}$ is the phonon wave vector.
When the lattice has inversion symmetry, ${\bf D}({\bf k})$ is a real matrix
and ${\bf D}({\bf k})={\bf D}(-{\bf k})$. Consequently, the polarization
vectors ${\bm{\epsilon}}_{{\bf k}\sigma}$ must be real, where $\sigma$ is the
branch index, and the system contains only pure vibration. In contrast, when
the inversion symmetry is broken, ${\bf D}({\bf k})$ is a complex matrix and
${\bm{\epsilon}}_{{\bf k}\sigma}$ could be complex. The system does not only
contain vibration but also rotation.

Classically, the displacement of atom in unit cell $l$ at ${\bf R}_{l}$ is
${\bf u}_{l}\sim{\text{Re}}[{\bm{\epsilon}}_{{\bf k}\sigma}e^{i({\bf
k}\cdot{\bf R}_{l}-\omega_{{\bf k}\sigma}t)}]$ for given state $(\sigma,{\bf
k})$ Ashcroft and Mermin (1976) where $\omega_{{\bf k}\sigma}$ is the phonon
frequency. For real polarization vector, ${\bf u}_{l}\sim{\bm{\epsilon}}_{{\bf
k}\sigma}{\text{cos}}({\bf k}\cdot{\bf R}_{l}-\omega_{{\bf k}\sigma}t)$. The
motion of each atom can be regarded as three “in-phase” harmonic oscillators
which gives rise to zero angular momentum because ${\bf u}_{l}$ is always
parallel to $\dot{\bf u}_{l}$. The case is different for complex polarization
vector, the displacement along $\mu$-direction ($\mu=x,y,z$) is

$\displaystyle u_{l}^{\mu}\sim$ $\displaystyle{\text{Re}}(\epsilon_{{\bf
k}\sigma}^{\mu}){\text{cos}}({\bf k}\cdot{\bf R}_{l}-\omega_{{\bf
k}\sigma}t)-{\text{Im}}(\epsilon_{{\bf k}\sigma}^{\mu}){\text{sin}}({\bf
k}\cdot{\bf R}_{l}-\omega_{{\bf k}\sigma}t)$ $\displaystyle=$
$\displaystyle{\text{cos}}({\bf k}\cdot{\bf R}_{l}-\omega_{{\bf
k}\sigma}t-\phi^{\mu}_{{\bf k}\sigma}),$ (1)

where the phase shift is determined by $\tan(\phi^{\mu}_{{\bf
k}\sigma})={\text{Im}}(\epsilon_{{\bf
k}\sigma}^{\mu})/{\text{Re}}(\epsilon_{{\bf k}\sigma}^{\mu})$. Eq. (1) shows
that the motion of atom can be regarded as three “out-of-phase” harmonic
oscillators, because $\phi^{\mu}_{{\bf k}\sigma}\neq\phi^{\nu}_{{\bf
k}\sigma}$ when $\mu\neq\nu$. In other words, atoms rotate around their
equilibrium positions circularly or elliptically. Such microscopic local
rotation gives rise to nonzero angular momentum, which has been experimentally
observed in WSe2 Zhu  _et al._ (2018), is termed as phonon angular momentum
(PAM) Zhang and Niu (2014, 2015).

A quantum mechanical theory of PAM was first given by McLellan McLellan
(1988). The overall PAM of a lattice with $n$ atoms in each unit cell can be
written as McLellan (1988):

${\bm{L}}=\sum_{l\kappa}{\bf u}_{l\kappa}\times(m_{\kappa}\dot{\bf
u}_{l\kappa}).$ (2)

$m_{\kappa}$ is the mass of $\kappa$-th atom with $\kappa=1,2,...n$. In 2014,
Zhang and Niu Zhang and Niu (2014) presented a comprehensive second
quantization form of PAM when the inversion symmetry is absent. They found
that, when the system is in equilibrium and has time-reversal symmetry, mean
PAM vanishes. Nonzero mean PAM can be obtained by two possible ways: (1)
breaking the time-reversal symmetry Zhang and Niu (2014); (2) driving the
system into non-equilibrium Hamada  _et al._ (2018); Park and Yang (2020);
Huang  _et al._ (2020). Later on Hamada et al. Hamada  _et al._ (2018) found a
nonzero PAM by using the Boltzmann transport equation under relaxation time
approximation when a temperature gradient was applied. The $\mu$-component of
mean PAM is calculated as $\langle L_{\mu}\rangle=\Lambda_{\mu\nu}\partial
T/\partial x_{\nu}$ where $\Lambda_{\mu\nu}$ is a response tensor. As a
result, a nonzero phonon magnetic moment due to PAM was calculated accordingly
Hamada  _et al._ (2018); Park and Yang (2020); Ren  _et al._ (2021).

However, the aforementioned calculation results are too small to explain the
recent observed chiral phonon activated Seebeck effect Kim  _et al._ (2022). A
significant magnetic effect was measured by the time-resolved magneto-optical
Kerr effect (TR MOKE) in chiral organic-inorganic hybrid perovskite with laser
pulse heating. Moreover, the current theory is not able to explain the
observed large magnetic moment in Dirac semimetal Cd3As2 Cheng  _et al._
(2020) and ErFeO3 Nova  _et al._ (2017). Therefore, there must be other
mechanism beyond Hamada’s calculation Hamada  _et al._ (2018) which can leads
to larger mean PAM and phonon magnetic moment.

In this Letter, we revisit the derivation from Eq. (2) and propose a new
mechanism to generate nonzero mean PAM by keeping both diagonal and off-
diagonal terms of density matrix in Kubo formula. We find that the off-
diagonal terms of PAM, which describe the quantum transition between different
phonon states, result in notable mean PAM only when the system is not in
equilibrium.

We start our study from the $\mu$-component of Eq. (2) which can be written as
Hamada (2021):

$\displaystyle L_{\mu}\approx\dfrac{\hbar}{2N}\sum\limits_{l}\sum\limits_{{\bf
k},{\bf k}^{\prime}}\sum\limits_{\sigma,\sigma^{\prime}}\left(\epsilon_{{\bf
k}\sigma}^{\dagger}M_{\mu}\epsilon_{{\bf
k}^{\prime}\sigma^{\prime}}\sqrt{\dfrac{\omega_{{\bf
k}^{\prime}\sigma^{\prime}}}{\omega_{{\bf k}\sigma}}}a_{{\bf
k}\sigma}^{\dagger}a_{{\bf k}^{\prime}\sigma^{\prime}}-\epsilon_{{\bf
k}^{\prime}\sigma^{\prime}}^{T}M_{\mu}\epsilon_{{\bf
k}\sigma}^{*}\sqrt{\frac{\omega_{{\bf
k}\sigma}}{\omega_{{\bf{k}^{\prime}\sigma^{\prime}}}}}a_{{\bf
k}^{\prime}\sigma^{\prime}}a_{{\bf k}\sigma}^{\dagger}\right)e^{i{\bf
R}_{l}\cdot({\bf k}^{\prime}-{\bf k})}.$ (3)

$a_{{\bf k}\sigma}^{\dagger}$ and $a_{{\bf k}\sigma}$ are the creation and
annihilation operators of phonons, respectively. Both $aa$ and
$a^{\dagger}a^{\dagger}$ terms are neglected since they vary rapidly with time
and have marginal contribution. When the system is a crystal, Eq. (3) can be
further simplified by using $\frac{1}{N}\sum_{l}e^{i{\bf R}_{l}\cdot({\bf
k}^{\prime}-{\bf k})}=\delta_{{\bf k},{\bf k}^{\prime}}$ and $-\epsilon_{{\bf
k}^{\prime}\sigma^{\prime}}^{T}M_{\mu}\epsilon_{{\bf
k}\sigma}^{*}=\epsilon_{{\bf k}\sigma}^{\dagger}M_{\mu}\epsilon_{{\bf
k}^{\prime}\sigma^{\prime}}$. Then $L_{\mu}$ can be divided into two parts:
the diagonal term ($L_{\mu}^{\text{D}}$) and off-diagonal term
($L_{\mu}^{\text{OD}}$) which can be written as

$L_{\mu}^{\text{D}}=\sum\limits_{{\bf k}\sigma}l_{{\bf
k}\sigma}^{\mu}\left[a_{{\bf k}\sigma}^{\dagger}a_{{\bf
k}\sigma}+\frac{1}{2}\right],$ (4a) $L_{\mu}^{\text{OD}}=\sum\limits_{\bf
k}\sum^{\sigma\neq\sigma^{\prime}}\limits_{\sigma\sigma^{\prime}}l^{\mu}_{{\bf
k}\sigma\sigma^{\prime}}a_{{\bf k}\sigma}^{\dagger}a_{{\bf
k}\sigma^{\prime}}.$ (4b)

We note the matrix elements in Eqs. (4a) and (4b) as follows:

$l_{{\bf k}\sigma}^{\mu}=\hbar(\epsilon_{{\bf
k}\sigma}^{\dagger}M_{\mu}\epsilon_{{\bf k}\sigma}),$ (5a) $l^{\mu}_{{\bf
k}\sigma\sigma^{\prime}}=\frac{\hbar}{2}\epsilon_{{\bf
k}\sigma}^{\dagger}M_{\mu}\epsilon_{{\bf
k}\sigma^{\prime}}\left(\sqrt{\frac{\omega_{{\bf
k}\sigma^{\prime}}}{\omega_{{\bf k}\sigma}}}+\sqrt{\frac{\omega_{{\bf
k}\sigma}}{\omega_{{\bf
k}\sigma^{\prime}}}}\right),(\sigma\neq\sigma^{\prime}).$ (5b)

where $M_{\mu}=I_{n\times n}\bigotimes(-i)\varepsilon_{\mu\nu\gamma}$
($\varepsilon$ is Levi-Civita tensor) and $\hbar$ is the Planck constant. The
diagonal term in Eq. (4a) and its matrix elements in Eq. (5a) have been
studied in our previous work Zhang and Niu (2014). The off-diagonal term in
Eq.(4b) and its matrix elements in Eq. (5b) show that the off-diagonal PAM
describes the quantum transition between state ($\sigma,{\bf k}$) and state
($\sigma^{\prime},{\bf k}$).

We modify the Kubo formula of thermal conductivity Kubo  _et al._ (1957);
Zwanzig (1965); Mori  _et al._ (1962), which is similar to the electrical Kubo
formula Kubo (1957), to calculate frequency-dependent mean PAM as:

$\displaystyle\langle L_{\mu}(\omega)\rangle$ $\displaystyle=$
$\displaystyle-V\sum\limits_{\nu}\dfrac{(\nabla
T)_{\nu}}{T}\sum\limits_{n}\dfrac{e^{-\beta
E_{n}}}{Z}\int_{0}^{\beta}d\lambda\int_{0}^{\infty}dte^{i(\omega+i\eta)t}$ (6)
$\displaystyle\times\langle S_{\nu}(-i\hbar\lambda)L_{\mu}(t)\rangle.$

Here $T$ is temperature, $\nabla T$ is temperature gradient, $V$ is volume,
$Z$ is the partition function, $\beta=1/k_{B}T$ with $k_{B}$ the Boltzmann
constant. $\lambda$ is a parameter. Detailed derivation of Eq. (6) is shown in
Appendix A. The heat current density operator (${\bf S}$) in Eq.(6) is also
written into diagonal part (${\bf S}^{\text{D}}$) and off-diagonal part (${\bf
S}^{\text{OD}}$) as follows Allen and Feldman (1993); Hardy (1963):

${\bf S}^{\text{D}}=\sum_{{\bf k}\sigma}{\bf s}_{\bf k\sigma}\left[a_{{\bf
k}\sigma}^{\dagger}a_{{\bf k}\sigma}+\frac{1}{2}\right],$ (7a) ${\bf
S}^{\text{OD}}=\sum_{{\bf
k}}\sum^{\sigma\neq\sigma^{\prime}}_{\sigma\sigma^{\prime}}{\bf s}_{\bf
k\sigma\sigma^{\prime}}a_{{\bf k}\sigma}^{\dagger}a_{{\bf k}\sigma^{\prime}}.$
(7b)

The matrix elements are

${\bf s}_{{\bf k}\sigma}=\dfrac{\hbar}{2V}\widetilde{\epsilon}_{{\bf
k}\sigma}^{\hskip 2.84544pt\dagger}\dfrac{\partial{\widetilde{\text{D}}(\bf
k)}}{\partial{\bf k}}\widetilde{\epsilon}_{{\bf k}{\sigma}},$ (8a) ${\bf
s}_{{\bf k}\sigma\sigma^{\prime}}=\dfrac{\hbar}{4V}\widetilde{\epsilon}_{{\bf
k}\sigma}^{\hskip 2.84544pt\dagger}\dfrac{\partial{\widetilde{\text{D}}(\bf
k)}}{\partial{\bf k}}\widetilde{\epsilon}_{{\bf
k}\sigma^{\prime}}\left(\sqrt{\dfrac{\omega_{{\bf
k}{\sigma^{\prime}}}}{\omega_{{\bf k}{\sigma}}}}+\sqrt{\dfrac{\omega_{{\bf
k}{\sigma}}}{\omega_{{\bf k}{\sigma^{\prime}}}}}\right),$ (8b)

where $\widetilde{\text{D}}(\bf k)$ is the redefined dynamical matrix and
$\widetilde{\epsilon}_{{\bf k}\sigma}$ is the corresponding eigenstate. Detail
derivations of above terms are given in Appendix B. It has been proved that
${\bf s}_{\bf k\sigma}=\frac{1}{V}\hbar\omega_{{\bf k}\sigma}{\bf v}_{{\bf
k}\sigma}$, which recover the conventional heat current operators Hardy
(1963), where ${\bf v}_{{\bf k}\sigma}$ is the phonon group velocity.

Substituting Eqs. (7a) and (7b) into Eq. (6), one can obtain the mean PAM per
unit volume as
$\langle{\mathcal{L}}_{\mu}(\omega)\rangle=\langle{\mathcal{L}}^{\text{D}}_{\mu}(\omega)\rangle+\langle{\mathcal{L}}^{\text{OD}}_{\mu}(\omega)\rangle$.
The diagonal part
$\langle{\mathcal{L}}_{\mu}^{\text{D}}(\omega)\rangle=\langle
L_{\mu}^{\text{D}}(\omega)\rangle/V$ and the off-diagonal part
$\langle{\mathcal{L}}_{\mu}^{\text{OD}}(\omega)\rangle=\langle
L_{\mu}^{\text{OD}}(\omega)\rangle/V$ are written as follows, respectively,

$\langle{\mathcal{L}}_{\mu}^{\text{D}}(\omega)\rangle=\sum_{\nu}\frac{i(\nabla
T)_{\nu}}{T}\frac{1}{\omega+i\eta}\sum\limits_{{\bf k}\sigma}\frac{\partial
f^{0}_{{\bf k}\sigma}}{\partial(\hbar\omega_{{\bf k}\sigma})}({\bf s}_{{\bf
k}\sigma})_{\nu}l^{\mu}_{{\bf
k\sigma}}=\sum_{\nu}\Lambda_{\mu\nu}^{\text{D}}(\omega)\frac{\partial
T}{\partial x_{\nu}},$ (9a)

$\langle{\mathcal{L}}_{\mu}^{\rm OD}(\omega)\rangle=-\sum_{\nu}\frac{i(\nabla
T)_{\nu}}{T}\sum\limits_{{\bf k}\sigma\sigma^{\prime}}\frac{f^{0}_{{\bf
k}\sigma}-f^{0}_{{\bf k}\sigma^{\prime}}}{\hbar(\omega_{{\bf
k}\sigma}-\omega_{{\bf k}\sigma^{\prime}})}\frac{({\bf s}_{{\bf
k}\sigma\sigma^{\prime}})_{\nu}l^{\mu}_{{\bf
k}\sigma^{\prime}\sigma}}{\omega_{\bf k\sigma}-\omega_{\bf
k\sigma^{\prime}}-\omega-i\eta}=\sum_{\nu}\Lambda_{\mu\nu}^{\text{OD}}(\omega)\frac{\partial
T}{\partial x_{\nu}},$ (9b)

The detailed derivations are shown in Appendix A.

We now numerically calculate the response tensor $\Lambda_{\mu\nu}$ in a dice
lattice model Wang and Ran (2011) as an example. We are interested in the dc
current in the $\omega\rightarrow 0$ limit and $\eta$ is simply approximate as
the phonon relaxation rate $1/\tau$ Allen and Feldman (1993). The detailed
parameters of material in calculation are given in Supporting Material.

The symmetry of the lattice we studied results in nonzero $xy$ and $yx$
components of the response tensor. We first calculate the $\Lambda_{yx}$ at
$T$ = 300 K with phonon relaxation time $\tau$ ranging from 0.1 ps to 10 ps as
shown in Fig. 1. As for the diagonal part, $\Lambda_{yx}^{\text{D}}$ varies
linearly with $\tau$, which is the same as the Boltzmann transport equation
calculations Hamada  _et al._ (2018). In our calculations, we find that
$\Lambda_{yx}^{\text{D}}\sim$ $-10^{-8}$ $\times$ [$\tau$/(1 s)] J s
${\text{m}}^{-2}$ ${\text{K}}^{-1}$. More interestingly, we find that the
contribution from the off-diagonal part $\Lambda_{yx}^{\text{OD}}$ is also
important, especially when $\tau<5$ ps. One can see that
$\Lambda^{\text{OD}}_{yx}$ decreases first and then increases with increasing
$\tau$. The minimum value of $\Lambda^{\text{OD}}_{yx}$ is determined by
$\eta\sim\omega_{{\bf k}\sigma}-\omega_{{\bf k}\sigma^{\prime}}$.

Figure 1: (color online) The calculated response tensor $\Lambda_{yx}$ (blue
solid line), $\Lambda_{yx}^{\text{D}}$ (green dashed line), and
$\Lambda_{yx}^{\text{O}\text{D}}$ (red dotted line) as functions of $\tau$
when $T$ = 300 K. Figure 2: (color online) (a) Contribution to
$\Lambda_{yx}^{\text{D}}$ from different phonon branches. (b) Contribution to
$\Lambda_{yx}^{\text{OD}}$ from transition between different phonon branches.

Figure 2 (a) shows the contributions from different phonon branches to
$\Lambda_{yx}^{\text{D}}$. The calculation results show that the contribution
from acoustic phonon is more significant than the contribution from optical
phonon. The reason is that the acoustic phonon has larger phonon group
velocity which is included in $s_{{\bf k}\sigma}$ in Eq. (9a). In order to
distinguish the contribution from the quantum transitions between different
phonon branches to off-diagonal PAM, we plot the calculated acoustic-acoustic
transition, acoustic-optical transition, and optical-optical transition
components of $\Lambda_{yx}^{\text{OD}}$ in Fig. 2 (b). The calculation
results show that the quantum transition between two optical branches is
dominant. The contributions from quantum transition between two acoustic
branches and that between one acoustic branch and one optical branch are
relatively small. This finding implies that more optical phonon branches are
helpful to enlarge the mean PAM. Therefore, chiral materials with many atoms
in unit cell, such as the chiral hybrid perovskites, is a good candidate to
achieve large PAM Kim  _et al._ (2022).

Figure 3: (color online). $\Lambda_{yx}^{\text{D}}$ and
$\Lambda_{yx}^{\text{OD}}$ vs temperature $T$ when $\tau=1$ ps.

Figure 3 shows the temperature dependence of response tensors. One can see
that when $T\rightarrow 0$, both $\Lambda_{yx}^{\text{D}}$ and
$\Lambda_{yx}^{\text{OD}}$ vanish. $\Lambda_{yx}^{\text{D}}$ at low
temperature is positive and increases slowly with increasing $T$. Then
$\Lambda_{yx}^{\text{D}}$ begins a gradual decline and becomes negative. When
$T$ is higher than the Debye temperature, $\Lambda_{yx}^{\text{D}}$ approaches
the high temperature limit. The calculation results show that
$\Lambda_{yx}^{\text{OD}}$ is more sensitive to temperature.
$\Lambda_{yx}^{\text{OD}}$ decreases rapidly with increasing temperature and
reaches a high temperature limit. The analytical expressions of response
tensor at low-temperature limit and at high-temperature limit as given in
Supporting Materials.

In summary, we calculate the phonon angular momentum induced by temperature
gradient using the Kubo formula. The off-diagonal elements of density matrix
due to the quantum transition between different phonon branches are found to
be important to generate the nonzero PAM. Both diagonal and off-diagonal
response tensors for certain dice lattice are calculated numerically as an
example. Chiral materials with many atoms in unit cell is find to be helpful
to observe large PAM and large phonon magnetic moment. In addition, we point
out that the off-diagonal PAM could be further enhanced in non-perfect
crystals in which the momentum conservation is absent, in other words, the
transition between $({\bf k},\sigma)$ and $({\bf k}^{\prime},\sigma^{\prime})$
should be also considered. It will be considered in our future work. To our
best knowledge, chiral organic-inorganic hybrid perovskites is a possible
candidate to observe large phonon PAM and phonon magnetic moment.

##  I Acknowledgments

This work is supported by National Natural Science Foundation of China (No.
11890703). JZ is also supported by the “Shuangchuang” Doctor program of
Jiangsu Province (JSSCBS20210341).

##  II Appendix A: Kubo formula for PAM

The thermal Kubo formula Kubo  _et al._ (1957) is different from the
electrical Kubo formulaKubo (1957). The reason is that there is a well-defined
external force which drives the electrical current. However, there is no
similar term involving the temperature gradient in the Hamiltonian to drive a
heat current. Thus the thermal Kubo formula requires an additional statistical
hypothesis Zwanzig (1965); Mori  _et al._ (1962), which assumes a local space-
dependent temperature $T(x)=[k_{\rm B}\beta(x)]^{-1}$. Then the local density
matrix is

$\rho=\dfrac{e^{-\int d^{3}x\beta(x)h(x)}}{Z},$ (A1)

where $h(x)$ is the Hamiltonian density operator, $Z$ is the partition
function, and the Hamiltonian $H=\int d^{3}xh(x)$. A heat current density
operator ${\bf S}(x)$ is now defined by the condition of local energy
conservation

$\dfrac{\partial h(x)}{\partial t}+\nabla\cdot{\bf S}(x)=0.$ (A2)

The total heat current operator is ${\bf S}=\frac{1}{V}\int d^{3}x{\bf S}(x)$.
If the temperature variation $\delta T(x)$ is weak, $\beta(x)$ can be written
as $\beta[1-\delta T(x)]/T$, where $(k_{\rm B}\beta)^{-1}$ is the average
temperature $T$. Then Eq. (A1) becomes $\rho=e^{-\beta(H+H^{\prime})}/Z$ where
the operator $H^{\prime}=-\frac{1}{T}\int d^{3}x\delta T(x)h(x)$ formally
plays the role of a perturbation. Using the integrated form of Eq. (A2), the
perturbation due to the temperature gradient is

$H^{\prime}=-\sum\limits_{\nu}\dfrac{(\nabla
T)_{\nu}}{T}\int_{-\infty}^{0}dt\int d^{3}xS_{\nu}(x,t).$ (A3)

The density matrix can be expended in powers of the perturbation by
$e^{-\beta(H+H^{\prime})}=e^{-\beta H}+e^{-\beta H}\int^{\beta}_{0}d\lambda
e^{\lambda H}H^{\prime}e^{-\lambda H}+\cdots=\rho_{0}+\rho_{1}$ with parameter
$\lambda$. It is obvious that ${\text{tr}}\rho_{0}{L_{\mu}}\equiv 0$, and the
nonzero PAM is determined by ${\text{tr}}\rho_{1}L_{\mu}$ which is

$\displaystyle\langle\ L_{\mu}\rangle$ $\displaystyle=$
$\displaystyle-V\sum\limits_{\nu}\dfrac{(\nabla
T)_{\nu}}{T}\sum\limits_{n}\dfrac{e^{-\beta
E_{n}}}{Z}\int_{0}^{\beta}d\lambda\int_{-\infty}^{0}dt$ (A4)
$\displaystyle\times\langle e^{\lambda H}S_{\nu}(t)e^{-\lambda
H}L_{\mu}(0)\rangle.$

$e^{\lambda H}Se^{-\lambda H}=S(-i\hbar\lambda)$, $E_{n}$ is the energy of
state $({\bf k},\sigma)$ which is noted as $n$ for simplicity. Then the Kubo
formula for PAM in Eq. (6) can be obtained Allen and Feldman (1989, 1993).
$\omega$ represents the frequency of an ac applied temperature gradient. This
Kubo expression can be written in the Lehman representation

$\displaystyle\langle L_{\mu}(\omega)\rangle=$
$\displaystyle-V\sum\limits_{\nu}\dfrac{(\nabla
T)_{\nu}}{T}\sum\limits_{n,m}\dfrac{e^{-\beta
E_{n}}}{Z}\int_{0}^{\beta}d\lambda\int_{0}^{\infty}dte^{i(\omega+i\eta)t}\langle
n|e^{\lambda H}S_{\nu}e^{-\lambda H}|m\rangle\langle
m|e^{\frac{iHt}{\hbar}}L_{\mu}e^{\frac{-iHt}{\hbar}}|n\rangle,$ (A5)
$\displaystyle=$ $\displaystyle V\sum\limits_{\nu}\dfrac{(\nabla
T)_{\nu}}{T}\sum\limits_{n,m}\dfrac{e^{-\beta
E_{n}}}{Z}\left[\dfrac{e^{\beta(E_{n}-E_{m})}-1}{E_{n}-E_{m}}\right]\left(\dfrac{i\hbar}{E_{n}-E_{m}-\hbar\omega-i\hbar\eta}\right)\langle
n|S_{\nu}|m\rangle\langle m|L_{\mu}|n\rangle.$

This expression can be divided into two parts: the contribution from diagonal
terms $\langle L_{\mu}^{\text{D}}(\omega)\rangle$ when $n=m$ and the
contribution from offdiagonal term $\langle
L_{\mu}^{\text{OD}}(\omega)\rangle$ when $n\neq m$.

$\displaystyle\langle L_{\mu}^{\rm OD}(\omega)\rangle=$ $\displaystyle
V\sum\limits_{\nu}\dfrac{(\nabla T)_{\nu}}{T}\sum\limits_{nm}\dfrac{e^{-\beta
E_{n}}}{Z}\left[\dfrac{e^{\beta(E_{n}-E_{m})}-1}{E_{n}-E_{m}}\right]\left(\dfrac{i\hbar}{E_{n}-E_{m}-\hbar\omega-i\hbar\eta}\right)$
(A6) $\displaystyle\times\sum\limits_{{\bf k}{\bf
k}^{\prime}\alpha\alpha^{\prime}\sigma\sigma^{\prime}}\langle n|a_{{\bf
k}^{\prime}\alpha}^{\dagger}a_{{\bf
k}^{\prime}\alpha^{\prime}}|m\rangle\langle m|a_{{\bf
k}\sigma}^{\dagger}a_{{\bf k}\sigma^{\prime}}|n\rangle({\bf s}_{{\bf
k}^{\prime}\alpha\alpha^{\prime}})_{\nu}{l}^{\mu}_{{\bf
k}\sigma\sigma^{\prime}}.$

The exact state $|m\rangle$ is the unique state which couples to $a_{{\bf
k}\sigma}^{\dagger}a_{{\bf k}\sigma^{\prime}}|n\rangle$, denoted as $|n_{{\bf
k}\sigma^{\prime}\sigma}\rangle$. We can use the standard expressions $\langle
m|a_{{\bf k}\sigma}^{\dagger}a_{{\bf
k}\sigma^{\prime}}|n\rangle=\sqrt{(n_{{\bf k}\sigma}+1)n_{{\bf
k}\sigma^{\prime}}}\delta_{mn_{{\bf k}\sigma^{\prime}\sigma}}$ and $\langle
n|a_{{\bf k}^{\prime}\alpha}^{\dagger}a_{{\bf
k}^{\prime}\alpha^{\prime}}|n_{{\bf
k}\sigma^{\prime}\sigma}\rangle=\sqrt{(n_{{\bf k}\sigma}+1)n_{{\bf
k}\sigma^{\prime}}}\delta_{\alpha^{\prime}\sigma}\delta_{\alpha\sigma^{\prime}}\delta_{{\bf
k}{\bf k}^{\prime}}$. In terms of the normal mode frequencies, $E_{n}-E_{m}$
becomes $\hbar(\omega_{{\bf k}\sigma^{\prime}}-\omega_{{\bf k}\sigma})$,
$\alpha=\sigma^{\prime}$, and $\alpha^{\prime}=\sigma$. Now since
$\sigma\neq\sigma^{\prime}$, we can write $\sum\limits_{n}\frac{e^{-\beta
E_{n}}}{Z}n_{{\bf k}\sigma^{\prime}}(n_{{\bf k}\sigma}+1)=\langle(n_{{\bf
k}\sigma}+1)n_{{\bf k}\sigma^{\prime}}\rangle=\langle n_{{\bf
k}\sigma^{\prime}}\rangle(\langle n_{{\bf k}\sigma}\rangle+1)$. Using the
identity $(e^{\beta\hbar(\omega_{{\bf k}\sigma^{\prime}}-\omega_{{\bf
k}\sigma})}-1)\langle n_{{\bf k}\sigma^{\prime}}\rangle(\langle n_{{\bf
k}\sigma}\rangle+1)=\langle n_{{\bf k}\sigma}\rangle-\langle n_{{\bf
k}\sigma^{\prime}}\rangle=f^{0}_{{\bf k}\sigma}-f^{0}_{{\bf
k}\sigma^{\prime}}$, $\langle{\mathcal{L}}_{\mu}^{\text{OD}}(\omega)\rangle$
in Eq. (9b) can be obtained by $\sigma\leftrightarrow\sigma^{\prime}$.

Considering $\langle
m|a^{\dagger}_{\textbf{{k}}\sigma}a_{\textbf{{k}}\sigma}|n\rangle=n_{\textbf{{k}}\sigma}\delta_{mn}$,
the contribution from diagonal term can be written as

$\displaystyle\langle L_{\mu}^{\rm
D}(\omega)\rangle=-V\beta\sum\limits_{\nu}\dfrac{(\nabla
T)_{\nu}}{T}\left(\dfrac{i}{\omega+i\eta}\right)\sum\limits_{{\bf k}{\bf
k}^{\prime}\sigma\alpha}$ $\displaystyle\left[\langle n_{{\bf
k}^{\prime}\alpha}n_{{\bf k}\sigma}\rangle(\textbf{{s}}_{{\bf
k}^{\prime}\alpha})_{\nu}l^{\mu}_{{\bf k}\sigma}+\frac{1}{2}\langle n_{{\bf
k}^{\prime}\alpha}\rangle({\bf s}_{{\bf k}^{\prime}\alpha})_{\nu}l^{\mu}_{{\bf
k}\sigma}+\frac{1}{2}({\bf s}_{{\bf k}^{\prime}\alpha})_{\nu}\langle n_{{\bf
k}\sigma}\rangle l^{\mu}_{{\bf k}\sigma}+\frac{1}{4}({\bf s}_{{\bf
k}^{\prime}\alpha})_{\nu}l^{\mu}_{{\bf k}\sigma}\right].$ (A7)

According the Wick’s theorem, the factor $\langle n_{{\bf
k}^{\prime}\alpha}n_{{\bf k}\sigma}\rangle$ can be written as $\langle n_{{\bf
k}^{\prime}\alpha}\rangle\langle n_{{\bf k}\sigma}\rangle$ plus a correction
$\langle n_{{\bf k}\sigma}\rangle(\langle n_{{\bf k}\sigma}\rangle+1)$. Using
the condition of zero average energy current and PAM in equilibrium and the
identity $\langle n_{{\bf k}\sigma}\rangle(\langle n_{{\bf
k}\sigma}\rangle+1)=-k_{\rm B}T\frac{\partial\langle n_{{\bf
k}\sigma}\rangle}{\partial(\hbar\omega_{{\bf k}\sigma})}$,
$\langle{\mathcal{L}}_{\mu}^{\text{D}}(\omega)\rangle$ in Eq. (9a) can be
obtained by dropping the last three terms in the bracket in Eq. (A7).

##  III Appendix B: Derivation of off-diagonal heat current density operator

Hardy Hardy (1963) has shown that the quadratic terms of the heat current
operator is

$\displaystyle{\bf S}$ $\displaystyle=$
$\displaystyle\dfrac{1}{2V}\sum\limits_{ll^{\prime}}\sum\limits_{\kappa\kappa^{\prime}}({\bf
R}_{l\kappa}-{\bf
R}_{l^{\prime}\kappa^{\prime}})\sum\limits_{\alpha}\left\\{\dfrac{p_{l\kappa\alpha}}{m_{\kappa}}\dfrac{1}{i\hbar}[p_{l\kappa\alpha},V({\bf
R}_{l^{\prime}\kappa^{\prime}})]\right.$ (B1)
$\displaystyle\left.+\dfrac{1}{i\hbar}[p_{l\kappa\alpha},V({\bf
R}_{l^{\prime}\kappa^{\prime}})]\dfrac{p_{l\kappa\alpha}}{m_{\kappa}}\right\\}.$

where $V({\bf R}_{l^{\prime}\kappa^{\prime}})$ is harmonic potential energy

$V({\bf
R}_{l^{\prime}\kappa^{\prime}})=\dfrac{1}{2}\sum\limits_{l}\sum\limits_{\kappa}\sum\limits_{\alpha\gamma}\Phi_{\alpha\gamma}(l\kappa,l^{\prime}\kappa^{\prime})u_{l\kappa\alpha}u_{l^{\prime}\kappa\prime\gamma}.$
(B2)

Here $\alpha,\gamma=x,y,z$. $\Phi$ is the force constant matrix. It is easy to
verify that $\frac{1}{i\hbar}[p_{l\kappa\alpha},V({\bf
R}_{l^{\prime}\kappa^{\prime}})]=-\frac{1}{2}\sum_{\beta}\Phi_{\alpha\gamma}(l\kappa,l^{\prime}\kappa^{\prime})u_{l^{\prime}\kappa^{\prime}\gamma}$,
then

${\bf
S}=\dfrac{1}{2V}\sum\limits_{ll^{\prime}}\sum\limits_{\kappa\kappa^{\prime}}\sum\limits_{\alpha\gamma}({\bf
R}_{l^{\prime}\kappa^{\prime}}-{\bf
R}_{l\kappa})\Phi_{\alpha\gamma}({l\kappa},{l^{\prime}\kappa^{\prime}})\dfrac{p_{l\kappa\alpha}}{m_{\kappa}}u_{l^{\prime}\kappa^{\prime}\gamma}.$
(B3)

By using the second quantization form of $u_{l^{\prime}\kappa^{\prime}\gamma}$
and $p_{l\kappa\alpha}$, we obtain

$\displaystyle{\bf S}\approx$
$\displaystyle\dfrac{1}{2V}\dfrac{i\hbar}{2N}\sum\limits_{ll^{\prime}}\sum\limits_{\kappa\kappa^{\prime}}\sum\limits_{{\bf
k}{\bf k}^{\prime}}\sum\limits_{\sigma\sigma^{\prime}}({\bf
R}_{l^{\prime}\kappa^{\prime}}-{\bf
R}_{l\kappa})\sum\limits_{\alpha\gamma}\dfrac{\Phi_{\alpha\gamma}({l\kappa},{l^{\prime}\kappa^{\prime}})}{\sqrt{m_{\kappa}m_{\kappa^{\prime}}}}$
(B4) $\displaystyle\times\epsilon_{{\bf
k}\sigma\kappa\alpha}^{*}\epsilon_{{{\bf
k}^{\prime}}\sigma^{\prime}\kappa^{\prime}\gamma}\left(\sqrt{\dfrac{\omega_{{\bf
k}^{\prime}\sigma^{\prime}}}{\omega_{{\bf k}\sigma}}}a_{{\bf
k}^{\prime}\sigma^{\prime}}a_{{\bf
k}\sigma}^{\dagger}+\sqrt{\dfrac{\omega_{{\bf k}\sigma}}{\omega_{{\bf
k}^{\prime}\sigma^{\prime}}}}a_{{\bf k}\sigma}^{\dagger}a_{{\bf
k}^{\prime}\sigma^{\prime}}\right)$ $\displaystyle\times e^{i({\bf
k}^{\prime}\cdot{\bf R}_{l^{\prime}\kappa^{\prime}}-{\bf k}\cdot{\bf
R}_{l\kappa})}\times e^{i({\bf k}\cdot{\bf d}_{k}-{\bf k^{\prime}}\cdot{\bf
d}_{k^{\prime}})}.$

Here ${\bf d}_{i}$ is the equilibrium position of the $i$th atom relative to
the center of the unit cell and we have ignored the terms like $aa$,
$a^{\dagger}a^{\dagger}$. Using the definition
$\widetilde{\text{D}}_{\alpha\gamma}({\bf
k^{\prime}},\kappa\kappa^{\prime})$=$\sum_{l^{\prime}}\Phi_{\alpha\gamma}(l\kappa,l^{\prime}\kappa^{\prime})e^{i{\bf
k^{\prime}}\cdot({\bf R}_{l^{\prime}\kappa^{\prime}}-{\bf
R}_{l\kappa})}/\sqrt{m_{\kappa}m_{\kappa^{\prime}}}$,
$\widetilde{\epsilon}_{{{\bf
k}^{\prime}}\sigma^{\prime}\kappa^{\prime}\gamma}=\epsilon_{{{\bf
k}^{\prime}}\sigma^{\prime}\kappa^{\prime}\gamma}e^{-i{\bf
k}^{\prime}\cdot{\bf d}_{k^{\prime}}}$ and $\frac{1}{N}\sum_{l}e^{i{\bf
R}_{l\kappa}\cdot(\textbf{{k}}^{\prime}-\textbf{{k}})}=\delta_{\textbf{{k}},\textbf{{k}}^{\prime}}$,
we rewrite Eq. (B4) as

$\displaystyle{\bf S}\approx$ $\displaystyle\dfrac{\hbar}{4V}\sum\limits_{\bf
k}\sum\limits_{\kappa\kappa^{\prime}}\sum\limits_{\sigma\sigma^{\prime}}\sum\limits_{\alpha\gamma}\widetilde{\epsilon}^{\hskip
2.84544pt*}_{{\bf
k}\sigma\kappa\alpha}\frac{\partial{\widetilde{\text{D}}_{\alpha\gamma}({\bf
k},\kappa\kappa^{\prime})}}{\partial\bf k}\widetilde{\epsilon}_{{\bf
k}\sigma^{\prime}\kappa^{\prime}\gamma}$ (B5)
$\displaystyle\times\left(\sqrt{\dfrac{\omega_{{\bf
k}\sigma^{\prime}}}{\omega_{{\bf k}\sigma}}}a_{{\bf k}\sigma^{\prime}}a_{{\bf
k}\sigma}^{\dagger}+\sqrt{\dfrac{\omega_{{\bf k}\sigma}}{\omega_{{\bf
k}\sigma^{\prime}}}}a_{{\bf k}\sigma}^{\dagger}a_{{\bf
k}\sigma^{\prime}}\right).$

Due to the communication relation
$[a_{\textbf{{k}}\sigma^{\prime}},a_{\textbf{{k}}\sigma}^{\dagger}]=\delta_{\sigma\sigma^{\prime}}$,
we finally get the diagonal term in Eq. (7a) and off-diagonal term in Eq. (7b)

## References

  * Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Thomson Learning, Inc., 1976) p. 453.
  * Zhu  _et al._ (2018) H. Zhu, J. Yi, M. Li, J. Xiao, L. Zhang, C. Yang, R. A. Kaindl, L. Li, Y. Wang, and X. Zhang, Science 359, 579 (2018).
  * Zhang and Niu (2014) L. Zhang and Q. Niu, Phys. Rev. Lett. 112, 085503 (2014).
  * Zhang and Niu (2015) L. Zhang and Q. Niu, Phys. Rev. Lett. 115, 115502 (2015).
  * McLellan (1988) A. G. McLellan, J. Phys. C: Solid State Phys. 21, 1177 (1988).
  * Hamada  _et al._ (2018) M. Hamada, E. Minamitani, M. Hirayama, and S. Murakami, Phys. Rev. Lett. 121, 175301 (2018).
  * Park and Yang (2020) S. Park and B.-J. Yang, Nano Lett. 20, 7694 (2020).
  * Huang  _et al._ (2020) Z. Huang, B. P. Bloom, X. Ni, Z. N. Georgieva, M. Marciesky, E. Vetter, F. Liu, D. H. Waldeck, and D. Sun, ACS Nano 14, 10370 (2020).
  * Ren  _et al._ (2021) Y. Ren, C. Xiao, D. Sapaov, and Q. Niu, Phys. Rev. Lett. 127, 186403 (2021).
  * Kim  _et al._ (2022) K. Kim, E. Vetter, L. Yan, C. Yang, Z. Wang, R. Sun, Y. Yang, A. Comstock, X. Li, J. Zhou, L. Zhang, W. You, D. Sun, and J. Liu, Nat. Mater., in revision (2022).
  * Cheng  _et al._ (2020) B. Cheng, T. Schumann, Y. Wang, X. Zhang, D. Barbalas, S. Stemmer, and N. P. Armitage, Nano Lett. 20, 5991 (2020).
  * Nova  _et al._ (2017) T. F. Nova, A. Cartella, A. Cantaluppi, M. Först, D. Bossini, R. V. Mikhaylovskiy, A. Kimel, R. Merlin, and A. Cavalleri, Nat. Phys. 13, 132 (2017).
  * Hamada (2021) M. Hamada,  _Ph.D. Thesis_ (Tokyo Institute of Technology, 2021).
  * Kubo  _et al._ (1957) R. Kubo, M. Yokota, and S. Nakajima, J. Phys. Soc. Jpn. 12, 1203 (1957).
  * Zwanzig (1965) R. Zwanzig, Ann. Rev. Phys. Chem. 16, 67 (1965).
  * Mori  _et al._ (1962) H. Mori, I. Oppenheim, and J. Ross,  _Studies in Statistical Mechanics_ , Vol. 1 (North Holland, Amsterdam, 1962) p. 213.
  * Kubo (1957) R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957).
  * Allen and Feldman (1993) P. B. Allen and J. L. Feldman, Phys. Rev. B 48, 12581 (1993).
  * Hardy (1963) R. J. Hardy, Phys. Rev. 132, 168 (1963).
  * Wang and Ran (2011) F. Wang and Y. Ran, Phys.rev.b 84, 277 (2011).
  * Allen and Feldman (1989) P. B. Allen and J. L. Feldman, Phys. Rev. Lett. 62, 645 (1989).


#  S-wave pairing in a two-orbital t-J model on triangular lattice: possible
application to Pb10-xCux(PO4)6O

Hanbit Oh1 Ya-Hui Zhang2<EMAIL_ADDRESS>1Department of Physics, KAIST,
Daejeon, 34126, Republic of Korea  2William H. Miller III Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, Maryland, 21218,
USA

###### Abstract

Recently room temperature superconductor was claimed in Pb10-xCux(PO4)6O (also
known as LK-99) with $x\in(0.9,1.1)$. Density functional theory (DFT)
calculations suggest that the conduction electrons are from the doped Cu atoms
with valence close to $d^{9}$. Motivated by this picture, we build a two-
orbital Hubbard model on triangular lattice formed by the $d_{xz}$ and
$d_{yz}$ orbitals with total hole density (summed over spin and orbital)
$n=1-p$. When $p=0$, the system is in a Mott insulator within this model. When
$p>0$, we derive a $t-J$ model and perform a self-consistent slave boson mean
field calculation. Interestingly we find a s wave pairing in contrast to one-
orbital t-J model which favors $d+id$ pairing. S wave pairing should be more
robust to disorder and may lead to high Tc superconductor with sufficiently
large value of $t$ and $J$. However, the DFT calculations predict a very small
value of $t$ and then the $T_{c}$ is expected to be small. If LK99 is really a
high Tc superconductor, ingredients beyond the current model is needed. We
conjecture that the doped Cu atoms may distort the original lattice and form
local clusters with smaller Cu -Cu distance and thus larger value of $t$ and
$J$. Within these clusters we may locally apply our t-J model calculation and
expect high Tc s-wave superconductor. Then the superconducting islands couple
together, which may eventually become a global superconductor, an insulator or
even an anomalous metal depending on sample details.

Introduction Recently there is report of room temperature superconductivity in
Pb10-xCux(PO4)6O, also called LK-99Lee  _et al._ (2023a, b). The experimental
reproduction of the exciting discovery is still ongoingLiu  _et al._ (2023);
Kumar  _et al._ (2023); Hou  _et al._ (2023); Wu  _et al._ (2023); Abramian
_et al._ (2023). In the theoretical sideGriffin (2023); Kurleto  _et al._
(2023); Si and Held (2023); Tavakol and Scaffidi (2023); Baskaran (2023),
density functional theory (DFT) calculations suggest a picture of narrow bands
formed by Cu 3d orbitals with the valence close to d9Griffin (2023); Kurleto
_et al._ (2023); Si and Held (2023), similar to the high Tc cuprates.

In this paper, we build a two-orbital model on triangular lattice based on the
Cu 3d $d_{xz}$ and $d_{yz}$ orbitals. Based on symmetry analysis, we find that
three tight binding parameters are allowed with only the nearest neighbor
hopping. We suggest one specific choice to fit the DFT band. To capture the
strong correlations from the Cu 3d orbitals, we propose a two-orbital Hubbard
model at total hole density $n=1-p$. $p=0$ is a Mott insulator and we derive a
two-orbital t-J model for the small finite $p$ regime. Then we apply the salve
boson mean field theory to analyze the t-J model. Slave boson theory is known
to be able to reproduce several essential properties of the superconductor in
cuprates, including the pairing symmetry and the doping dependence of the
pairing strengthLee  _et al._ (2006). Based on slave boson mean field
calculation, we identify a s′-wave spin-singlet pairing with both intra-
orbital and inter-orbital components, which is unusual given that usually t-J
model calculation predicts higher angular momentum pairing. A s-wave pairing
is more robust to disorder and thus may support higher Tc at the same level of
pairing strength, compared to $d$ wave pairing in cuprates.

While high Tc superconductor exists in this model with sufficiently large
value of $t$ and $J$, we need a much higher value of $t$ than the DFT
prediction to obtain a superconductor at the order of $100$ K. If LK-99 is
indeed a high temperature (or even room temperature) superconductor, new
ingredients must be included beyond the DFT calculations to greatly enhance
the mobility of the electrons. One conjecture we have is that the doped Cu
atoms may distort the original lattice and form clusters with much smaller
Cu-O or Cu-Cu distance. In this picture, we have regions which is locally
described by our t-J model with large value of t and J. Then we can have
s-wave high Tc superconductor islands according to our calculation. These
islands then need to couple together through the Josephson tunneling. Whether
the system can be a global superconductor depends on details such as disorder
strength and may vary a lot in different samples.

Figure 1: (a) Illustrations of atomic structure with Cu (red)and Pb (blue)
atoms. The Cu itself forms a hexagonal structure with primitive vectors
$\vec{a}_{i}$. Two d-orbitals ($d_{xz},d_{yz}$) are living at each Cu site.
(b) Energy band structure of a tight binding model (in hole picture). We set
$t_{\sigma}=-1,t_{\pi}=0,t_{2}=-\sqrt{3}/4$. The Fermi surfaces at total hole
filling $n=1$ are illustrated in the inset. There is only a single $C_{3}$
symmetric electron-pocket near the $\Gamma$ point.

Model We consider a two-orbital based on the $d_{xz}$ and $d_{yz}$ orbitals,
living on layered triangular lattice. For simplicity, we will ignore the
inter-layer tunneling $t_{z}$ in this paper and focus on a two-dimensional
(2D) model on a triangular lattice, as illustrated in Fig1 (a).

We label $d_{i;1\sigma}$ and $d_{i;2\sigma}$ as the annihilation operator of a
hole in $d_{xz}$ and $d_{yz}$ orbital respectively at site $i$ for spin
$\sigma=\uparrow,\downarrow$. For convenience we use t he hole picture with
the hole density
$n_{i}=\sum_{a=1,2}\sum_{\sigma=\uparrow,\downarrow}d^{\dagger}_{i;a\sigma}d_{i;a\sigma}$
has an average value of $n=1-p$ per site. We will mainly consider the small
positive $p$ regime. We define
$\psi_{i;\sigma}=(d_{i;1\sigma},d_{i;2\sigma})^{T}$. Our system have a $C_{3}$
rotation symmetry around the site $i$. There is also a mirror reflection
symmetry $\mathcal{M}_{y}$: $x\rightarrow-x,y\rightarrow y$. Under $C_{3}$ and
$\mathcal{M}_{y}$, we have $\psi_{i}\rightarrow U_{C_{3}}\psi_{C_{3}i}$ and
$\psi_{i}\rightarrow U_{M_{y}}\psi_{M_{y}i}$. We have
$U_{C_{3}}=\begin{pmatrix}-\frac{1}{2}&\frac{\sqrt{3}}{2}\\\
-\frac{\sqrt{3}}{2}&-\frac{1}{2}\end{pmatrix}$ and
$U_{M_{y}}=\begin{pmatrix}-1&0\\\ 0&1\end{pmatrix}$.

We consider the generic form of a tight-binding model with nearest-neighbor
hopping,

$\displaystyle H_{K}$ $\displaystyle=$
$\displaystyle-\sum_{\sigma=\uparrow,\downarrow}\sum_{i;l=1,2,3,4,5,6}\Psi^{\dagger}_{i;\sigma}\bm{T}_{l}\Psi_{i+\vec{a}_{l};\sigma}$
(1)

with $2\times 2$ matrix of hopping matrix $\bm{T}_{l}=\bm{T}_{-l}^{\dagger}$
as required by the Hermitian condition. Here, we used the lattice vectors
which connect the nearest neighbor,
$\vec{a}_{1}=(1,0),\vec{a}_{2}=(1,\sqrt{3})/2,\vec{a}_{3}=(-1,\sqrt{3})/2$.

The hopping matrix is derived by the central symmetries of the system, $C_{3}$
rotation, and mirror symmetry $\mathcal{M}_{y}$. This constrains the hopping
matrix to be in the form,

$T_{1}=\begin{pmatrix}t_{\sigma}&-t_{2}\\\ t_{2}&t_{\pi}\end{pmatrix},$ (2)
$T_{2}=\begin{pmatrix}\frac{1}{4}(t_{\sigma}+3t_{\pi})&t_{2}+\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})\\\
-t_{2}+\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})&\frac{1}{4}(3t_{\sigma}+t_{\pi})\end{pmatrix}$
(3)

and

$T_{3}=\begin{pmatrix}\frac{1}{4}(t_{\sigma}+3t_{\pi})&-t_{2}-\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})\\\
t_{2}-\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})&\frac{1}{4}(3t_{\sigma}+t_{\pi})\end{pmatrix}$
(4)

In the above $t_{\sigma}$ and $t_{\pi}$ are from the $\sigma$ and $\pi$ bond.
$t_{2}$ arises from breaking the $C_{6}$ rotation symmetry. We note that
$t_{2}$ is necessary to split the two-fold degeneracy at K and K′ point in the
Brillouin zone (BZ). While we focus on the simple 2D model, extension of it to
3D is straightforward by simply adding a term
$-t_{z}\sum_{i}\Psi^{\dagger}_{i;\sigma}\Psi_{i+{\widehat{z}};\sigma}$.

In the Supplemental Material (SM), we added more discussion to emphasize the
symmetry action on the two orbitals. In Figure 1(b), we illustrate the band
structure model, Eq. 1 for a specific choice of hopping
parameters:$t_{\sigma}=-1$, $t_{\pi}=0$ and $t_{2}=-\frac{\sqrt{3}}{4}$. Here
we are using the hole picture and add a negative sign to the hopping. In
particular, at total hole filling $n=1$, there is a single electron
pocket111The pocket looks like a hole pocket in the plot, but because we are
using the hole picture, this is an electron pocket. near $\Gamma$ pocket whose
Fermi-surface shape is shown in the inset. We also note that a two-orbital
model was recently already proposedTavakol and Scaffidi (2023). But the model
there seems to be not equivalent to our model at any parameter.

To also incorporate the strong on-site repulsion, we consider a Hubbard model:

$H_{Hubbard}=H_{K}+\frac{U}{2}\sum_{i}n_{i}(n_{i}-1)$ (5)

Note that here we ignore the difference between the intra-orbital and inter-
orbital repulsion and also the Hund’s coupling.

t-J model When $n=1$, at large U/t regime the system is in a Mott insulator
and described by a spin-orbital model at low energy. The model can be obtained
from the standard t/U expansion. To make the calculation in the next part
convenient, here we represent the spin-orbital model using the Abrikosov-
fermion representation. We assume $f^{\dagger}_{i;a\sigma}$ creates one hole
with orbital $a$ and spin $\sigma$ at the site $i$. With the constraint
$n_{i;f}=1$, we can recover the spin-orbital model Hilbert space. Then the
spin-orbital model is written as:

$H_{J}=\frac{J}{2}\sum_{i}\sum_{l=1,2,3}T_{l;aa^{\prime}}T^{\dagger}_{l;bb^{\prime}}(f_{i;a\sigma}^{\dagger}f_{i;b^{\prime}\sigma^{\prime}})(f^{\dagger}_{i+a_{l};b\sigma^{\prime}}f_{i+a_{l};a^{\prime}\sigma})$
(6)

where $J=\frac{4}{U}$ assuming $t_{\sigma}=1$. We have assumed Einstein
summation convention.

At filling $n=1-p$, the low energy physics is described by a $t-J$ model:

$H_{t-J}=PH_{K}P+H_{J}$ (7)

where $P$ is the projection operator to remove the double occupancy.

We will focus on the electron doped side with $p>0$. The other side with hole
density $n>1$ may be different in the sense that the additional holes may
enter oxygen or Pb atoms, similar to the hole doped cuprates. We leave it to
future to model the hole doped side.

Slave boson mean field theory We use the standard slave boson constructionLee
_et al._ (2006): $d_{i;a\sigma}=b_{i}^{\dagger}f_{i;a\sigma}$ with the
constraint $n_{i;b}+n_{i;f}=1$. On average $n_{i;f}=1-p$ and $n_{i;b}=p$.
Assuming that the slave boson condenses with $\langle b_{i}\rangle=\sqrt{p}$,
we can get the mean field equation for $f$. In the following for convenience
we still use the convention
$\psi_{i;\sigma}=(f_{i;1\sigma},f_{i;2\sigma})^{T}$.

We have the mean field ansatz:

$H_{M}=H_{0}+H_{D}$ (8)

where

$H_{0}=-\sum_{\sigma}\sum_{i}\sum_{l=1,2,3,4,5,6}\psi_{i;\sigma}^{\dagger}(pT_{l}+C_{l})\psi_{i+a_{l};\sigma}$
(9)

The Hermitian condition again constrains that $T_{-l}=T_{l}^{\dagger}$ and
$C_{-l}=C_{l}^{\dagger}$. $C_{l}$ is a $2\times 2$ complex matrix from
decoupling of the super-exchange J term. Note also $-l$ means the opposite
direction of $a_{l}$.

Meanwhile, we have the pairing term:

$H_{D}=-\sum_{i}\sum_{l=1,2,3,4,5,6}\psi^{\dagger}_{i;\uparrow}D_{l}\psi^{\dagger}_{i+a_{l};\downarrow}+H.C.$
(10)

We will restrict to spin-singlet pairing, which gives us the constraint that
$D_{-l}=D_{l}^{T}$.

We have the self consistent equations:

$C_{l}=\frac{J}{2}T_{l}\text{Tr}(\chi_{l}^{\dagger}T_{l}^{*})$ (11)

and

$D_{l}=\frac{J}{4}T_{l}\Delta_{l}^{T}T^{*}_{l}$ (12)

with the definition:

$\displaystyle\chi_{l;ab}$
$\displaystyle=2\langle\psi^{\dagger}_{i;a\uparrow}\psi_{i+a_{l};b\uparrow}\rangle$
$\displaystyle\Delta_{l;ab}$
$\displaystyle=2\langle\psi_{i;a\downarrow}\psi_{i+a_{l};b\uparrow}\rangle$
(13)

s wave pairing We perform the self-consistent calculation at $J/t=0.5$. We
start from initial ansatz with $C_{l}$ and $D_{l}$ as random complex matrix.
Nevertheless, we always reach a time reversal invariant ansatz:
$D_{1}=\begin{pmatrix}0.0121&-0.0044\\\ 0.0044&-0.0016\end{pmatrix}$,
$D_{2}=\begin{pmatrix}0.0018&0.0104\\\ 0.0015&0.0087\end{pmatrix}$ and
$D_{3}=\begin{pmatrix}0.0018&-0.0104\\\ -0.0015&0.0087\end{pmatrix}$ for
$p=0.1$.

Under $C_{3}$, we know $\psi\rightarrow U\psi$ with
$U=\begin{pmatrix}-\frac{1}{2}&\frac{\sqrt{3}}{2}\\\
-\frac{\sqrt{3}}{2}&-\frac{1}{2}\end{pmatrix}$ and then $D_{l}\rightarrow
U^{\dagger}D_{l}U^{*}$. We have checked that our ansatz is $C_{3}$ invariant,
like a s-wave pairing. Similarly it is invariant under $M_{y}$. We can project
the pairing term $D_{ij}$ to the lower band and get a scalar
$\Delta(\mathbf{k})$ in momentum space, so the projected Hamiltonian is
$\sum_{k}\Delta(\mathbf{k})c^{\dagger}_{\uparrow}(\mathbf{k})c^{\dagger}_{\downarrow}(-\mathbf{k})$,
where $c_{\sigma}(\mathbf{k})$ is the projected operator to the lower band.
$\Delta(\mathbf{k})$ is real as guaranteed by the time reversal symmetry:
$\psi_{i;a\uparrow}\rightarrow\psi_{i;a\downarrow},\psi_{i;a\downarrow}\rightarrow-\psi_{i;a\uparrow}$,
which must act as $c_{\uparrow}(\mathbf{k})\rightarrow
c_{\downarrow}(-\mathbf{k}),c_{\downarrow}(\mathbf{k})\rightarrow-
c_{\uparrow}(-\mathbf{k})$ projected to the lower band. Meanwhile in the lower
band $C_{3}$ acts trivially and we simply have a requirement
$\Delta(C_{3}\mathbf{k})=\Delta(\mathbf{k})$, in agreement with a s-wave
pairing. In Fig. 2 we show that the minimal gap in momentum space is non-zero,
suggesting that there is no node.

Figure 2: (a) Superconductor gap dependence on doping $p$. $\Delta_{m}$ is the
minimal gap in momentum space. $\Delta_{m}>0$ suggests that there is no node.
(b) At each $\theta$, the radius $r(\theta)$ of the blue line indicates the
minimal gap $\Delta_{min}(\theta)$ along this direction. We can see that the
gap is non-zero and almost uniform along all directions. Given that the ansatz
is time reversal invariant, this is a s-wave pairing.

The s-wave pairing ansatz is robust to variation of $t_{2}$ and $t_{\pi}$. We
also tried to relax the constraint $D_{-l}=D_{l}^{T}$, but the convergent
solution always satisfy this equation, suggesting that the s-wave spin-singlet
pairing is the appropriate solution. We also show the dependence of the
pairing gap on doping $p$ in Fig. 2(a).

Discussion In the slave boson treatment of the two-orbital t-J model, we find
a s-wave pairing and its pairing strength decreases with the doping level $p$,
similar to the solution found in the simple t-J model on square lattice. It is
known that the real $T_{c}$ is decided by the phase stiffness at small doping,
thus we expect a superconductor dome. At $p=0.1$, the superconductor gap is
$\Delta\approx 0.06t_{\sigma}$ at $J/t=0.5$. The DFT calculations suggest that
the bandwidth is $W\approx 130$ meV, which suggests that $t_{\sigma}\approx
20$ meV. Then the gap is only $1.2$ meV and $T_{c}$ may be only at order of
$10$ K. Besides, if $t$ is too small, then $J/t\sim\frac{t}{U}$ can not reach
$0.5$ and $T_{c}$ must be even much smaller. We note that the $p=0$ insulator
may be a charge transfer insulator and there may be other paths to generate a
larger $J$, but $J/t$ seems to be already quite significant. Even assuming
$J/t=0.5$, to get a $T_{c}$ around 100 K, we need to increase the hopping $t$
by at least one order of magnitude.

Therefore it appears that the model with the value predicted by DFT can not
explain the potential high Tc superconductor. New ingredients are needed to
enhance the mobility of electron to reach a high temperature superconductor in
our model. One may wonder whether a flat band is beneficial in the phonon
driven mechanism. We note that if the hopping is too small, the system is deep
inside a strong Mott insulator at $p=0$. Then at finite $p$, the doped
electrons will have even reduced hopping due to Mott physics and may tend to
localize. It is likely that the Bardeen–Cooper–Schrieffer theory (BCS) theory
can not be applied in this case. Due to strong repulsion, a phonon-driven
strong superconductor at the small hopping regime does not seem very likely. A
larger hopping is probably also needed to reach a high Tc superconductor in
the phonon scenario.

So how can we get a higher mobility? We note that the small hopping in this
system is largely due to the large Cu-Cu distance and Cu-oxygen distance.
LK-99 system differs from other superconductor materials in the following
aspects. In conventional material like cuprates, the doped atoms only provide
additional charges and usually do not play any essential roles. The electrons
are still moving in the crystal formed by the original atoms. In contrast, in
the LK-99 system, one Pb atom per unit cell is replaced with Cu atom. The
doped electrons are moving along the array formed by the Cu atoms. The DFT
calculations assume that the Cu atoms just form a periodic crystal, but this
assumption is highly questionable. The real chemistry is likely much more
complicated. We conjecture a scenario that the Cu atoms may greatly distort
the original crystal and form clusters where the lattice is compressed in a
local region and Cu-O and Cu-Cu distances are smaller. The hopping in these
clusters can thus be large and locally form s wave pairing according to our
calculation. Then these superconductor islands need to couple together through
Josephson tunneling. In the ideal case, there is a global condensation.
However, this is not guaranteed. Similar to previous studies of disorder
driven superconductor to insulator transition, the system may also be
insulating or in the more exotic scenario forms an anomalous metalKapitulnik
_et al._ (2019). Such an anomalous metal has preformed pairs and may explain
certain experimental result with diamagnetism, but finite resistivityLee  _et
al._ (2023b).

Summary In conclusion, we propose a two-orbital t-J model to describe the
essential physics of LK-99, the candidate material with possible room
temperature superconductivity. Our theoretical calculation predicts a s wave
pairing. Given that t-J model usually predicts a higher angular momentum
pairing, the model is certainly conceptual interesting and worth future
theoretical and numerical studies. It is also interesting to explore the
possibility that phonon cooperates with the super-exchange $J$ to further
enhance the pairing strength. The relevance of the model to LK-99 remains to
be investigated in future experimental studies. To obtain a superconductor
with $T_{c}\sim 100$K, we need the value of $t$ to be at least one magnitude
larger than the predicted values from DFT calculations. The exact mechanism to
generate such a large hopping is not clear now, but may be related to the
distortion caused by the Cu atoms. We conjecture that there may be local
region which is compressed and thus has a large hopping and a strong s wave
pairing according to our model. This picture suggest that LK-99 may be in the
category of granular superconductor. If true, there is likely strong sample
dependence with possible superconductor, insulator and even anomalous metal
depending on sample details. We also propose to realize our model and the
possible high Tc s-wave superconductor in other materials with active $d_{xz}$
and $d_{yz}$ orbitals.

Acknowledgement YHZ was supported by the National Science Foundation under
Grant No.DMR-2237031.

## References

  * Lee  _et al._ (2023a) S. Lee, J.-H. Kim, and Y.-W. Kwon, “The first room-temperature ambient-pressure superconductor,” (2023a), arXiv:2307.12008 [cond-mat.supr-con] .
  * Lee  _et al._ (2023b) S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, “Superconductor pb10-xcux(po4)6o showing levitation at room temperature and atmospheric pressure and mechanism,” (2023b), arXiv:2307.12037 [cond-mat.supr-con] .
  * Liu  _et al._ (2023) L. Liu, Z. Meng, X. Wang, H. Chen, Z. Duan, X. Zhou, H. Yan, P. Qin, and Z. Liu, arXiv preprint arXiv:2307.16802 (2023).
  * Kumar  _et al._ (2023) K. Kumar, N. Karn, and V. Awana, arXiv preprint arXiv:2307.16402 (2023).
  * Hou  _et al._ (2023) Q. Hou, W. Wei, X. Zhou, Y. Sun, and Z. Shi, arXiv preprint arXiv:2308.01192 (2023).
  * Wu  _et al._ (2023) H. Wu, L. Yang, B. Xiao, and H. Chang, arXiv preprint arXiv:2308.01516 (2023).
  * Abramian  _et al._ (2023) P. Abramian, A. Kuzanyan, V. Nikoghosyan, S. Teknowijoyo, and A. Gulian, arXiv preprint arXiv:2308.01723 (2023).
  * Griffin (2023) S. M. Griffin, arXiv preprint arXiv:2307.16892 (2023).
  * Kurleto  _et al._ (2023) R. Kurleto, S. Lany, D. Pashov, S. Acharya, M. van Schilfgaarde, and D. S. Dessau, arXiv preprint arXiv:2308.00698 (2023).
  * Si and Held (2023) L. Si and K. Held, arXiv preprint arXiv:2308.00676 (2023).
  * Tavakol and Scaffidi (2023) O. Tavakol and T. Scaffidi, arXiv preprint arXiv:2308.01315 (2023).
  * Baskaran (2023) G. Baskaran, arXiv preprint arXiv:2308.01307 (2023).
  * Lee  _et al._ (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Reviews of modern physics 78, 17 (2006).
  * Note (1) The pocket looks like a hole pocket in the plot, but because we are using the hole picture, this is an electron pocket.
  * Kapitulnik  _et al._ (2019) A. Kapitulnik, S. A. Kivelson, and B. Spivak, Reviews of Modern Physics 91, 011002 (2019).

##  Appendix A Details on Symmetry analysis

We consider the generic form of tight-binding model with an nearest-neighbor
hopping,

$\displaystyle H_{K}$ $\displaystyle=$
$\displaystyle-\sum_{\sigma=\uparrow,\downarrow}\sum_{i,l}\Psi^{\dagger}_{i;\sigma}\bm{T}_{l}\Psi_{i+\vec{a}_{l};\sigma}+H.c,$
(14)

with $2\times 2$ matrix of hopping matrix $\bm{T}_{l}=\bm{T}_{-l}^{\dagger}$.
We start with a general form of $T_{1}$,

$T_{1}=\frac{1}{2}\begin{pmatrix}t_{\sigma}&t_{1}-t_{2}\\\
t_{1}+t_{2}&t_{\pi}\end{pmatrix}$ (15)

One can find that the mirror symmetry ($\mathcal{M}_{y}$) imposes $t_{1}=0$,
since $T_{-1}=U_{M_{y}}^{\dagger}T_{1}U_{M_{y}}$ while
$T_{-1}=T_{1}^{\dagger}$. We have $U_{M_{y}}=\begin{pmatrix}-1&0\\\
0&1\end{pmatrix}$. It is easy to see that we need $t_{1}=0$.

Hence we introduce a general tight-binding model with three free variables
$t_{\sigma},t_{\pi},t_{2}$. From $T_{1}$ we can generate $T_{2}$ and $T_{3}$
by applying the $C_{3}$ transformation:
$T_{3}=U^{\dagger}_{C_{3}}T_{1}U_{C_{3}}$ and
$T_{-2}=T_{2}^{\dagger}=U^{\dagger}_{C_{3}}T_{3}U_{C_{3}}$. The final tight
binding model is constrained to be:

$\displaystyle H_{K}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{\sigma}\sum_{i,l}\Psi_{i}^{\dagger}\left(\begin{array}[]{cc}t_{\sigma}&-t_{2}\\\
t_{2}&t_{\pi}\end{array}\right)\Psi_{i+\vec{a}_{1}}+H.c.$ (24)
$\displaystyle+\Psi_{i}^{\dagger}\left(\begin{array}[]{cc}\frac{1}{4}(t_{\sigma}+3t_{\pi})&t_{2}+\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})\\\
-t_{2}+\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})&\frac{1}{4}(3t_{\sigma}+t_{\pi})\end{array}\right)\Psi_{i+\vec{a}_{2}}+H.c.$
$\displaystyle+\Psi_{i}^{\dagger}\left(\begin{array}[]{cc}\frac{1}{4}(t_{\sigma}+3t_{\pi})&-t_{2}-\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})\\\
t_{2}-\frac{\sqrt{3}}{4}(t_{\sigma}-t_{\pi})&\frac{1}{4}(3t_{\sigma}+t_{\pi})\end{array}\right)\Psi_{i+\vec{a}_{3}}+H.c.$
$\displaystyle=$
$\displaystyle\sum_{\sigma}\Psi^{\dagger}_{k;\sigma}\left(\begin{array}[]{cc}h_{11}(k)&h_{12}(k)\\\
h_{21}(k)&h_{22}(k)\end{array}\right)\Psi_{k;\sigma}$ (27)

with

$\displaystyle h_{ab}(k)$ $\displaystyle=$
$\displaystyle-\sum_{l}\left[(t_{\sigma}\phi^{a}_{\sigma,l}\phi^{b}_{\sigma,l}+t_{\pi}\sum_{l}\phi^{a}_{\pi,l}\phi^{b}_{\pi,l})\cos(k\cdot\vec{a}_{l})+t_{2}\sigma^{ab}_{y}s^{l}\sin(k\cdot\vec{a}_{l})\right]$
(28)

with $\phi_{\sigma,l}=\vec{a}_{l}/a$, and

$\displaystyle\phi_{\sigma,1}=(1,0),\quad\phi_{\sigma,2}=(\frac{1}{2},\frac{\sqrt{3}}{2}),\quad\phi_{\sigma,3}=(-\frac{1}{2},\frac{\sqrt{3}}{2}),$
(29)
$\displaystyle\phi_{\pi,1}=(0,1),\quad\phi_{\pi,2}=(-\frac{\sqrt{3}}{2},\frac{1}{2}),\quad\phi_{\pi,3}=(-\frac{\sqrt{3}}{2},-\frac{1}{2}).$
(30)

with $s_{1}=-s_{2}=s_{3}=1$.

##  Appendix B Derivation of the self consistent equations.

###  B.1 mean field Hamiltonian

We first write down the most general mean field ansatz. We assume translation
symmetry and spin rotation symmetry( only spin-singlet pairing). But allow
general other ansatz.

We introduce the notation that
$\psi_{i;\sigma}=(f_{i;a\sigma},f_{i;b\sigma})^{T}$. Ourm ean field can be
written as:

$H_{M}=H_{0}+H_{D}$ (31)

where

$H_{0}=-\sum_{i,\sigma}\sum_{l=1,2,3,4,5,6}\psi_{i;\sigma}^{\dagger}(pT_{l}+C_{l})\psi_{i+a_{l};\sigma}$
(32)

We always have constraint that $T_{-l}=T_{l}^{\dagger}$ and
$C_{-l}=C_{l}^{\dagger}$. $C_{l}$ is a $2\times 2$ complex matrix from
decoupling of the super-exchange J term.

Meanwhile, we have the pairing term:

$H_{D}=-\sum_{i}\sum_{l=1,2,3,4,5,6}\psi^{\dagger}_{i;\uparrow}D_{l}\psi^{*}_{i+a_{l};\downarrow}+H.C.$
(33)

We have the constraint that $D_{-l}=D_{l}^{T}$ for spin-singlet pairing.

Here $D_{l}$ is a $2\times 2$ matrix decoupled from the J term.

###  B.2 Self consistent equations

Before we continue, we define the following two $2\times 2$ complex matrix:

$\displaystyle\chi_{l;ab}$
$\displaystyle=2\langle\psi^{\dagger}_{i;a\uparrow}\psi_{i+a_{l};b\uparrow}\rangle$
$\displaystyle\Delta_{l;ab}$
$\displaystyle=2\langle\psi_{i;a\downarrow}\psi_{i+a_{l};b\uparrow}\rangle$
(34)

We fix the overall hopping term energy scale to be $1$. Then we define
$J=\frac{4}{U}$ to match the usual definition. From the second order
perturbation, we can write the $J$ term as:

$H_{J}=-\frac{2}{U}\sum_{i}\sum_{l=1,2,3}\psi_{i;a\sigma}^{\dagger}T_{l;aa^{\prime}}\psi_{i+a_{l};a^{\prime}\sigma}\psi^{\dagger}_{i+a_{l};b\sigma^{\prime}}T^{\dagger}_{l;bb^{\prime}}\psi_{i;b^{\prime}\sigma^{\prime}}$
(35)

Then we can get the following self consistent equation for the hopping term:

$C_{l;aa^{\prime}}=T_{l;aa^{\prime}}\frac{J}{2}2\sum_{b,b^{\prime}}\langle\psi^{\dagger}_{i+a_{l};b\uparrow}T^{\dagger}_{l;bb^{\prime}}\psi_{i;b^{\prime}\uparrow}\rangle$
(36)

In matrix notation, it is basically:

$C_{l}=\frac{J}{2}T_{l}\sum_{b,b^{\prime}}\chi^{\dagger}_{l;bb^{\prime}}T^{*}_{l;b^{\prime}b}=\frac{J}{2}T_{l}\mathrm{Tr}\left(\chi^{\dagger}_{l}T^{*}_{l}\right)$
(37)

For pairing, we can get:

$D_{l;ab}=\frac{J}{4}\sum_{a^{\prime}b^{\prime}}T_{l;aa^{\prime}}T^{\dagger}_{l;bb^{\prime}}\Delta_{l;b^{\prime}a^{\prime}}$
(38)

which is equivalent to:

$D_{l}=\frac{J}{4}T_{l}\Delta_{l}^{T}T^{*}_{l}$ (39)


[a]Matteo Buzzegoli

# Spin polarization in heavy-ion collisions induced by thermal vorticity and
thermal shear

###### Abstract

The vorticity is a quantity defined in a relativistic fluid that describes how
much a fluid element is rotating and is accelerating. By measuring the spin
polarization of hadrons, it was found that the quark gluon plasma produced in
heavy-ion collisions is the most “vorticous” fluid ever observed. More
generally, this opens the possibility to study the physics of QCD matter using
spin. Here I use statistical quantum field theory applied to a fluid in local
thermodynamic equilibrium to show how to connect the average spin of a fermion
with hydrodynamic quantities, and in particular with the thermal vorticity and
the thermal shear. I show that the spin polarization of a Dirac particle
induced by thermal vorticity is related to the gravitational (in medium) form
factor related to spin-rotation coupling. For these reasons, as we are
understanding the role of spin in hydrodynamics and in heavy-ion collisions,
spin is becoming a promising tool to investigate the properties of QCD and
whose applications are just begun to be explored.

##  1 Introduction

One of the goals of relativistic heavy-ion collisions is to study the physics
of finite temperature QCD. By measuring the momentum spectrum of particles and
comparing with the theoretical predictions of momentum distributions and
correlations, we can infer what happens during a collision. Most notably, in
this way it was found that the nuclear system goes through a phase of matter
named quark gluon plasma (QGP) that behaves like a nearly perfect fluid.
Before the measurement of spin polarization of $\Lambda$ particles in 2017
[1], all the studies were based on the observation of momentum.

This paradigm started to change when it was realized [2, 3] that off central
collisions should generate a large angular momentum and that particles, due to
the spin-orbit coupling, should have their spin aligned along the direction of
the total orbital angular momentum. Later, describing the spin degrees of
freedom as close to the local thermal equilibrium (in the same fashion as
successfully done for the momenta), the spin polarization was connected to the
collective motion of the plasma, and precisely to the thermal vorticity [4, 5,
6]. The observation of global spin polarization of $\Lambda$s [1] is in good
agreement with the hydrodynamic predictions, for more detail see for instance
the review [7].

More recently, the measurement of the global polarization of $\Xi$ hyperons
[8] further confirmed the hydrodynamic picture since it showed that spin
polarization in heavy-ion collisions, unlike $pp$ collisions, is not sensitive
to specific hadron properties. Over the past decade there has been a
significant development of spin physics for heavy-ion collisions, and spin
polarization has also been measured at very high [9, 10] and at very low [11,
12] energy.

The first measurement of spin polarization as a function of momentum, referred
to as local spin polarization, was reported in [13] and is in disagreement
with the predictions obtained from thermal vorticity. As I will show below,
the agreement with the hydrodynamic picture is restored once it is realized
that also other collective motions, in particular the shear flow, can induce a
spin polarization. I will show how to use quantum statistical field theory to
derive the connections between the spin polarization and the hydrodynamic
quantities, and I will discuss the features of the shear induced polarization.
I will then connect the vorticity induced polarization with the spin-orbit
coupling and the related form factor.

##  2 Spin polarization in quantum statistical field theory

The spin of a massive relativistic particle in quantum mechanics is defined
through the Pauli-Lubanski operator

$\widehat{S}^{\mu}=-\frac{1}{2m}\epsilon^{\mu\nu\rho\sigma}{\widehat{J}}_{\nu\rho}\widehat{P}_{\sigma},$
(1)

where ${\widehat{J}}$ is the angular momentum-boost operator and
${\widehat{P}}$ the energy-momentum operator. The spin polarization measured
in heavy-ion collisions is the average of the Pauli-Lubanski operator for a
particle of given momentum $p$, denoted by $S^{\mu}(p)$, and it is obtained as
the mean value of the operator (1) with the density matrix operator
${\widehat{\rho}}$, that is

$S^{\mu}(p)={\rm Tr}(\,{\widehat{\rho}}\,\widehat{S}^{\mu}(p)\,).$ (2)

It is convenient to relate the average spin vector defined in Eq. (2) with the
covariant Wigner functions, which, for a the Dirac field $\Psi$, are defined
as

$\displaystyle W_{\alpha\beta}(x,p)$ $\displaystyle=$
$\displaystyle\int\frac{{\rm d}^{4}y}{(2\pi)^{4}}\;{\rm e}^{-ip\cdot
y}\langle:{\overline{\Psi}}_{\beta}\left(x+\frac{y}{2}\right)\Psi_{\alpha}\left(x-\frac{y}{2}\right):\rangle\,,$
(3)

where $\alpha$ and $\beta$ are spinorial indices, $:\;:$ denotes the normal
ordering of creation and destruction operators, and angle bracket
$\langle\>\rangle$ denotes the trace with the density matrix operator. When
the momentum $p$ is future time-like, the Wigner functions contains only the
contribution form the particle part of the Dirac field, and the following
traces

$\mathcal{F}_{+}(x,p)={\rm
tr}\left[\,W(x,p)\theta(p^{2})\theta(p^{0})\right],\quad\mathcal{A}_{+}^{\mu}(x,p)={\rm
tr}\left[\gamma^{\mu}\gamma^{5}W(x,p)\theta(p^{2})\theta(p^{0})\right],$ (4)

are the particle part of the scalar and axial Wigner functions. For a weakly
interacting field, it can be shown [14] that the average spin vector can be
evaluated from the integrals of $\mathcal{F}_{+}$ and $\mathcal{A}_{+}^{\mu}$
over a space-like hyper-surface $\Sigma$ as follows

$S^{\mu}(p)=\frac{1}{2}\frac{\int_{\Sigma}{\rm d}\Sigma\cdot
p\;\mathcal{A}_{+}^{\mu}(x,p)}{\int_{\Sigma}{\rm d}\Sigma\cdot
p\;\mathcal{F}_{+}(x,p)}.$ (5)

Predictions for the spin polarization can be obtained from this formula if the
statistical ensemble, described by a density matrix ${\widehat{\rho}}$, is
known.

The covariant density matrix at the hypersurface $\Sigma$ is derived assuming
that the system reached the local thermodynamic equilibrium at an initial
stage and maximizing the entropy of the system. Neglecting the dissipative
effects, this procedure yields [15]

${\widehat{\rho}}_{\rm LE}=\frac{1}{Z}\exp\left[-\int_{\Sigma}{\rm
d}\Sigma_{\mu}(y)\left({\widehat{T}}^{\mu\nu}(y)\beta_{\nu}(y)-\zeta(y){\widehat{j}}^{\mu}(y)\right)\right],$
(6)

with ${\widehat{T}}$ the energy momentum tensor, ${\widehat{j}}$ a conserved
current, $\beta=u/T$ is the inverse four temperature vector and $\zeta=\mu/T$
where $\mu$ is the chemical potential. This operator is needed to obtain the
Wigner functions (4) at the point $x$. Since the interaction scales are
shorter than the macroscopic ones, the Wigner functions (4) can be obtained
expanding the hydrodynamic fields in the operator (6) around the same point
$x$ where they are to be evaluated, that is, setting $\zeta=0$ for simplicity,

$\beta_{\nu}(y)\simeq\beta_{\nu}(x)+\partial_{\lambda}\beta_{\nu}(x)(y-x)^{\lambda}.$
(7)

The resulting density matrix is

$\widehat{\rho}_{\rm
LE}\simeq\frac{1}{Z}\exp\left[-\beta_{\nu}(x)\widehat{P}^{\nu}+\frac{1}{2}\varpi_{\mu\nu}(x)\widehat{J}^{\mu\nu}_{x}-\frac{1}{2}\xi_{\mu\nu}(x)\widehat{Q}^{\mu\nu}_{x}+\cdots\right]\,,$
(8)

where the antisymmetric derivative of $\beta$ is the thermal vorticity

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

and the symmetric one is the thermal shear

$\xi_{\mu\nu}=\frac{1}{2}\left(\partial_{\mu}\beta_{\nu}+\partial_{\nu}\beta_{\mu}\right).$
(10)

Thermal vorticity couples with the conserved angular momentum operator

$\widehat{J}^{\mu\nu}_{x}=\\!\int\\!{\rm
d}\Sigma_{\lambda}\left[(y-x)^{\mu}\widehat{T}^{\lambda\nu}(y)-(y-x)^{\nu}\widehat{T}^{\lambda\mu}(y)\right],$
(11)

while the thermal shear couples with a non-conserved symmetric quadrupole like
operator

$\widehat{Q}^{\mu\nu}_{x}=\\!\int_{\Sigma_{FO}}\\!{\rm
d}\Sigma_{\lambda}\left[(y-x)^{\mu}\widehat{T}^{\lambda\nu}(y)+(y-x)^{\nu}\widehat{T}^{\lambda\mu}(y)\right].$
(12)

The definition of thermal vorticity (9) implies that $\varpi$ is related to
the acceleration of the fluid $a^{\mu}=u\cdot\partial u^{\mu}$, the
relativistic angular velocity
$\omega^{\mu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}u_{\sigma}\partial_{\nu}u_{\rho}$
and to the temperature $T=1/\sqrt{\beta^{2}}$ as follows

$\varpi^{\mu\nu}=-\frac{1}{2}\left[\partial^{\mu}\left(\frac{1}{T}\right)u^{\nu}-\partial^{\nu}\left(\frac{1}{T}\right)u^{\mu}\right]+\epsilon^{\mu\nu\rho\sigma}\frac{\omega_{\rho}}{T}u_{\sigma}+\frac{1}{2T}\left(a^{\mu}u^{\nu}-a^{\nu}u^{\mu}\right)\,.$
(13)

While the thermal shear is expressed as

$\xi^{\mu\nu}=\frac{1}{2}\left[\partial^{\mu}\left(\frac{1}{T}\right)u^{\nu}+\partial^{\nu}\left(\frac{1}{T}\right)u^{\mu}\right]+\frac{1}{2T}\left(a^{\mu}u^{\nu}+a^{\nu}u^{\mu}\right)+\frac{1}{T}\sigma^{\mu\nu}+\frac{1}{3T}\theta\Delta^{\mu\nu}\,,$
(14)

where $\nabla_{\mu}=\partial_{\mu}-u_{\mu}(u\cdot\partial)$,
$\Delta^{\mu\nu}=\eta^{\mu\nu}-u^{\mu}u^{\nu}$, $\theta=\nabla\cdot u$ and
$\sigma$ is the shear tensor

$\sigma^{\mu\nu}=\frac{1}{2}\left(\nabla^{\mu}u^{\nu}+\nabla^{\nu}u^{\mu}\right)-\frac{1}{3}\Delta^{\mu\nu}\theta\,.$
(15)

Using the statistical operator in Eq. (8) and linear response theory, the
average of the spin vector (5) for a free Dirac field result in [16]

$S^{\mu}(p)=-\frac{\epsilon^{\mu\rho\sigma\tau}p_{\tau}}{8m\int_{\Sigma}{\rm
d}\Sigma\cdot p\;n_{F}}\int_{\Sigma}{\rm d}\Sigma\cdot
p\;n_{F}(1-n_{F})\left[\varpi_{\rho\sigma}+2\hat{t}_{\rho}\frac{p^{\lambda}}{E_{p}}\xi_{\lambda\sigma}\right]\,,$
(16)

where $\hat{t}$ is the time direction in the laboratory frame and $n_{F}$ is
the Fermi-Dirac phase-space distribution function:

$n_{F}=\frac{1}{\exp[\beta\cdot p-\mu q]+1},$ (17)

where $q$ is the charge of the particle.

This recently found [17, 16] contribution of spin polarization from thermal
shear does not affect the predictions of global spin polarization that are
well reproduced by thermal vorticity alone, but it is the dominant
contribution in local spin polarization [18, 19]. Quantitative agreement with
the data of local spin polarization has been found if hadronization occurs at
nearly constant temperature [19], as expected at the energies where data has
been taken. Further analysis of the impact of thermal shear on spin
polarization [20, 21, 22] revealed a significant dependence on the initial
conditions of the hydrodynamic model. As the theory for spin polarization does
not require any additional parameter compared to the standard model of heavy-
ion collisions, the measurements of local spin polarization have the potential
to further constraint the initial conditions of the hydrodynamic equations and
to probe the properties of the QGP, such as the shear and bulk viscosities.

##  3 In medium form factors

As mentioned, the spin polarization can be understood as a result of spin-
rotation coupling. Stating that a medium is rotating is equivalent of saying
that there is an effective force causing an acceleration resulting in the
circular motion. The interaction of a single constituent of the medium with
the overall rotation can then be studied as an inertial effect caused by a
rotating observer. The analysis of the Dirac equation in rotating coordinates
revealed [23, 24] that a Dirac particle in a medium rotating with angular
velocity $\vec{\Omega}$ has the Hamiltonian

$H=H_{0}-\vec{J}\cdot\vec{\Omega},$ (18)

where $\vec{J}=\vec{L}+\vec{S}$ is the total angular momentum of the fermion
and $H_{0}$ is the Hamiltonian in absence of rotation. This is analogous to
the interaction of a fermion with an external magnetic field, resulting in the
Zeeman effect. It follows that the energy is lower when the spin is aligned
along the rotation, hence spin is polarized along the rotation when the system
is close to equilibrium.

More generally, we can describe the interaction with weak gravitational field
through the Lagrangian

$\mathcal{L}_{int}=\frac{1}{2}(g_{\mu\nu}-\eta_{\mu\nu})T^{\mu\nu},$ (19)

which reproduces the Hamiltonian (18) if we choose $g$ as the metric of
rotating coordinates. Using scattering theory, the coupling between spin and
rotation can be studied looking at the matrix elements of the energy momentum
tensor, that can be decomposed in terms of the form factors $f_{1}$ and
$g_{\Omega}$ as follows

$\langle
p^{\prime},\,s^{\prime}|\widehat{T}^{\mu\nu}(0)|p,\,s\rangle=\bar{u}(p^{\prime},s^{\prime})\left[f_{1}(q^{2})\frac{P^{\mu}P^{\nu}}{m}+g_{\Omega}(q^{2})\frac{\sigma^{(\mu\alpha}q_{\alpha}P^{\nu)}}{2m}+O(q^{2})\right]u(p,s),$
(20)

where $|p,\,s\rangle$ is a one particle state with momentum $p$ and spin $s$,
$u(p,s)$ and $\bar{u}(p^{\prime},s^{\prime})$ are spinors, $q=p^{\prime}-p$
and $P=(p+p^{\prime})/2$. The form factor $g_{\Omega}$ is also called the
gravitomagnetic moment and it can be shown [25] that describes exactly the
spin rotation coupling in (18). It has long been known that, differently form
other form factors like the magnetic moment, $g_{\Omega}$ can not receive
radiative corrections and should always be $g_{\Omega}=1$ otherwise the theory
is not Lorentz covariant and would violate the Einstein Equivalence Principle
(EEP) [26, 27, 28].

This restriction is lifted when considering a system in thermal equilibrium.
Indeed, the presence of a thermal bath breaks the Lorentz invariance and as a
result finite temperature field theory violates the EEP [29, 30]. As a
consequence, the matrix elements of an operator in a medium can have
additional structure and additional form factors [25, 31]. In addition to
four-vectors $P$ and $q$, matrix elements can also be given in terms of $u$,
that denotes the fluid velocity, and the vector
$l^{\mu}=\epsilon^{\mu\nu\rho\sigma}u_{\nu}P_{\rho}q_{\sigma}$. The energy-
momentum tensor is then decomposed as

$\begin{split}\langle
p^{\prime},s^{\prime}|\widehat{T}_{\mu\nu}(0)&|p,s\rangle=\bar{u}(p^{\prime},s^{\prime})\Big{\\{}I_{P\gamma}(P,q)\left(P_{\mu}\gamma_{\nu}+P_{\nu}\gamma_{\mu}\right)+I_{u\gamma}(P,q)\left(u_{\mu}\gamma_{\nu}+u_{\nu}\gamma_{\mu}\right)+\\\
&+I_{Pl}(P,q)\not{\hat{l}}\left(P_{\mu}\hat{l}_{\nu}+P_{\nu}\hat{l}_{\mu}\right)+I_{ul}(P,q)\not{\hat{l}}\left(u_{\mu}\hat{l}_{\nu}+u_{\nu}\hat{l}_{\mu}\right)\Big{\\}}u(p,s)+\cdots,\end{split}$
(21)

where $\hat{l}^{\mu}=l^{\mu}/\sqrt{-l^{2}}$ and the in-medium spin-rotation
coupling is [25]

$g_{\Omega}(P,q)=4\left(I_{P\gamma}(P,q)+\frac{I_{u\gamma}(P,q)}{(P\cdot
u)}-I_{Pl}(P,q)-\frac{I_{ul}(P,q)}{(P\cdot u)}\right).$ (22)

Explicit calculations in finite temperature QED [25], revealed that indeed the
gravitomagnetic moment receives radiative corrections, for instance, in the
low temperature limit $T\ll m$, the 1-loop gravitomagnetic moment is

$\lim_{\begin{subarray}{c}q\to 0\\\ P^{2}\to(P\cdot
u)^{2}\end{subarray}}g_{\Omega}=1-\frac{1}{6}\frac{e^{2}T^{2}}{m^{2}}.$ (23)

In [25] it was shown that $g_{\Omega}$ is directly connected to the axial
vortical effect, that is the generation of an axial current induced by the
rotation of the medium. Here I use the same method to show how spin-rotation
coupling is related to the spin polarization (2). In order to apply the
formula (5), I first need to compute the Wigner function in a system in
thermal equilibrium with thermal vorticity. Such a system is described by the
density matrix in Eq. (6). The contribution of thermal vorticity $\varpi$ to
the Wigner functions (3), obtained using linear response theory to the density
matrix in Eq. (6), is

$\Delta_{\varpi}W_{ab}^{+}(x,k)=-\int_{0}^{1}{\rm d}z\int_{\Sigma}{\rm
d}\Sigma_{\lambda}(y)\varpi_{\rho\kappa}(y-x)^{\kappa}\langle{\widehat{W}}_{ab}^{+}(x,k){\widehat{T}}^{\lambda\rho}(y+iz\beta(x))\rangle_{c,\beta},$
(24)

where the bracket $\langle\;\rangle_{c,\beta}$ denotes the connected
correlator traced with the density matrix in absence of thermal vorticity

${\widehat{\rho}}_{\beta}=\frac{1}{Z_{\beta}}\exp\left[-{\widehat{P}}\cdot\beta(x)-\zeta(x){\widehat{Q}}\right].$
(25)

To obtain the spin vector (5) we need the axial part of (24), which, according
to the definition (4), is

$\begin{split}\Delta_{\varpi}\mathcal{A}_{+}^{\mu}(x,k)=&-\int\frac{{\rm
d}^{4}s}{(2\pi)^{4}}\int_{0}^{1}{\rm d}z\int_{\Sigma}{\rm
d}\Sigma_{\lambda}(y)\varpi_{\rho\kappa}(y-x)^{\kappa}{\rm e}^{-ik\cdot s}\\\
&\times\langle\bar{\psi}^{+}(x+s/2)\gamma^{\mu}\gamma^{5}\psi^{+}(x-s/2){\widehat{T}}^{\lambda\rho}(y+iz\beta(x))\rangle_{c,\beta}.\end{split}$
(26)

Expanding the operators inside the correlators in terms of multi-particle
states, see [25], the leading contribution of the Wigner function is given in
terms of the matrix elements of the axial current

$\langle
q^{\prime},\tau|{\widehat{j}}_{A}^{\mu}(0)|q,\tau\rangle=\bar{u}_{\tau}(q^{\prime})A^{\mu}(q,q^{\prime})u_{\tau}(q)=\bar{u}_{\tau}(q^{\prime})\left[F_{A1}\gamma^{\mu}\gamma^{5}+F_{A2}\frac{(q^{\prime}-q)^{\mu}}{2m}\gamma^{5}\right]u_{\tau}(q),$
(27)

and of the energy momentum tensor:

$\langle
q^{\prime},\tau|{\widehat{T}}^{\mu\nu}(0)|q,\tau\rangle=\bar{u}_{\tau}(q^{\prime})M^{\mu\nu}\left((q^{\prime}+q)/2,q^{\prime}-q\right)u_{\tau}(q).$
(28)

As the system is supposed to be in the hydrodynamic regime, we have a
separation of scales between microscopic and macroscopic interactions. In this
approximation, the short-distance interactions are taken into account in the
matrix elements (27) and (28) which take place at finite temperature. The
long-distance correlations are taken into account in the thermal expectation
values of the creation and annihilation operators
$\langle\widehat{a}_{\tau^{\prime}}^{\dagger}(q^{\prime})\widehat{a}_{\tau}(q)\rangle=\delta_{\tau,\tau^{\prime}}\delta^{3}(\vec{q}-\vec{q}^{\prime})n_{F}(\beta\cdot
q-\zeta)$. One eventually obtains that

$\begin{split}\Delta_{\varpi}\mathcal{A}_{+}^{\mu}(x,k)=&F_{A1}(k,0)\varpi_{\kappa\rho}u_{\lambda}\Delta^{\kappa}_{\;\kappa^{\prime}}\frac{\theta(k_{0})}{2\varepsilon_{k}}\frac{\delta(k^{2}-m^{2})}{(2\pi)^{3}}n_{F}(k)(1-n_{F}(k))\\\
&\times\left\\{i\frac{\partial}{\partial p^{\prime}_{\kappa^{\prime}}}{\rm
tr}\left[M^{\lambda\rho}\left((p^{\prime}+k)/2,p^{\prime}-k\right)(\not{k}+m)\gamma^{\mu}\gamma^{5}(\not{p}^{\prime}+m)\right]\right\\}_{p^{\prime}=k},\end{split}$
(29)

where
$\Delta^{\kappa}_{\kappa^{\prime}}=\eta^{\kappa}_{\kappa^{\prime}}-u^{\kappa}u_{\kappa^{\prime}}$.
Finally, plugging the Eq. (21) in Eq. (29) and tracing gives the axial Wigner
function induced by thermal vorticity:

$\Delta_{\varpi}\mathcal{A}_{+}^{\mu}(x,k)=-g_{\Omega}(k,0)F_{A1}(k,0)\frac{\theta(k_{0})\delta(k^{2}-m^{2})}{(2\pi)^{3}}n_{F}(k)(1-n_{F}(k))\epsilon^{\mu\nu\sigma\tau}\varpi_{\nu\sigma}k_{\tau}.$
(30)

The spin polarization of a free field, using $g_{\Omega}(k,0)=1$,
$F_{A1}(k,0)=1$, and applying the formula (5), is the same as Eq. (16). For an
interacting field, the form factor $g_{\Omega}$ can receive radiative
corrections, for instance like in Eq. (23), and the spin polarization induced
by thermal vorticity can also receive radiative corrections.

##  4 Summary and outlook

I showed how spin polarization in a medium arises from the flow of the plasma
and I discussed on how it depends on the properties of such medium. I also
showed how the spin polarization contains information about the gravitational
form factor related to the spin-rotation coupling. Therefore, measurements of
the spin polarization in heavy-ion collisions provide a new tool to probe the
properties of QCD matter at finite temperature.

In addition to probe the shear and vortical flow of the fluid and to test the
dynamics of spin, several ideas on how to use spin polarization for
investigating the QGP have already been proposed. For instance, spin
polarization can reveal local parity violations in QCD [32, 33], it can probe
the presence of the critical point [34], and it can be used to estimate the
energy loss of a jet traversing the medium [35, 36].

Acknowledgments. M.B. thanks the U.S. Department of Energy for the support
through the Grant No. DE-SC0023692.

## References

  * [1] STAR collaboration, _Global $\Lambda$ hyperon polarization in nuclear collisions: evidence for the most vortical fluid_, _Nature_ 548 (2017) 62 [1701.06657].
  * [2] Z.-T. Liang and X.-N. Wang, _Globally polarized quark-gluon plasma in non-central A+A collisions_ , _Phys. Rev. Lett._ 94 (2005) 102301 [nucl-th/0410079].
  * [3] S.A. Voloshin, _Polarized secondary particles in unpolarized high energy hadron-hadron collisions?_ , nucl-th/0410089.
  * [4] F. Becattini, F. Piccinini and J. Rizzo, _Angular momentum conservation in heavy ion collisions at very high energy_ , _Phys. Rev. C_ 77 (2008) 024906 [0711.1253].
  * [5] F. Becattini and F. Piccinini, _The Ideal relativistic spinning gas: Polarization and spectra_ , _Annals Phys._ 323 (2008) 2452 [0710.5694].
  * [6] F. Becattini, V. Chandra, L. Del Zanna and E. Grossi, _Relativistic distribution function for particles with spin at local thermodynamical equilibrium_ , _Annals Phys._ 338 (2013) 32 [1303.3431].
  * [7] F. Becattini and M.A. Lisa, _Polarization and Vorticity in the Quark–Gluon Plasma_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 395 [2003.03640].
  * [8] STAR collaboration, _Global Polarization of $\Xi$ and $\Omega$ Hyperons in Au+Au Collisions at $\sqrt{s_{NN}}$ = 200 GeV_, _Phys. Rev. Lett._ 126 (2021) 162301 [2012.13601].
  * [9] ALICE collaboration, _Global polarization of $\Lambda\bar{\Lambda}$ hyperons in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 and 5.02 TeV_, _Phys. Rev. C_ 101 (2020) 044611 [1909.01281].
  * [10] ALICE collaboration, _Polarization of $\Lambda$ and $\bar{\Lambda}$ Hyperons along the Beam Direction in Pb-Pb Collisions at $\sqrt{s_{NN}}$=5.02 TeV_, _Phys. Rev. Lett._ 128 (2022) 172005 [2107.11183].
  * [11] STAR collaboration, _Global $\Lambda$-hyperon polarization in Au+Au collisions at $\sqrt{s_{NN}}$=3 GeV_, _Phys. Rev. C_ 104 L061901 [2108.00044].
  * [12] HADES collaboration, _$\Lambda$ Polarization in Au+Au Collisions at $\sqrt{s}_{NN}=2.4\,{\text{GeV}}$ Measured with HADES_, _Springer Proc. Phys._ 250 (2020) 435.
  * [13] STAR collaboration, _Polarization of $\Lambda$ ($\bar{\Lambda}$) hyperons along the beam direction in Au+Au collisions at $\sqrt{s_{{}_{NN}}}$ = 200 GeV_, _Phys. Rev. Lett._ 123 (2019) 132301 [1905.11917].
  * [14] F. Becattini, _Polarization in relativistic fluids: a quantum field theoretical derivation_ , _Lect. Notes Phys._ 987 (2021) 15 [2004.04050].
  * [15] F. Becattini, M. Buzzegoli and E. Grossi, _Reworking the Zubarev’s approach to non-equilibrium quantum statistical mechanics_ , _Particles_ 2 (2019) 197 [1902.01089].
  * [16] F. Becattini, M. Buzzegoli and A. Palermo, _Spin-thermal shear coupling in a relativistic fluid_ , _Phys. Lett. B_ 820 (2021) 136519 [2103.10917].
  * [17] S.Y.F. Liu and Y. Yin, _Spin polarization induced by the hydrodynamic gradients_ , _JHEP_ 07 (2021) 188 [2103.09200].
  * [18] B. Fu, S.Y.F. Liu, L. Pang, H. Song and Y. Yin, _Shear-Induced Spin Polarization in Heavy-Ion Collisions_ , _Phys. Rev. Lett._ 127 (2021) 142301 [2103.10403].
  * [19] F. Becattini, M. Buzzegoli, G. Inghirami, I. Karpenko and A. Palermo, _Local Polarization and Isothermal Local Equilibrium in Relativistic Heavy Ion Collisions_ , _Phys. Rev. Lett._ 127 (2021) 272302 [2103.14621].
  * [20] S. Ryu, V. Jupic and C. Shen, _Probing early-time longitudinal dynamics with the $\Lambda$ hyperon’s spin polarization in relativistic heavy-ion collisions_, _Phys. Rev. C_ 104 (2021) 054908 [2106.08125].
  * [21] S. Alzhrani, S. Ryu and C. Shen, _$\Lambda$ spin polarization in event-by-event relativistic heavy-ion collisions_, _Phys. Rev. C_ 106 (2022) 014905 [2203.15718].
  * [22] X.-Y. Wu, C. Yi, G.-Y. Qin and S. Pu, _Local and global polarization of $\Lambda$ hyperons across RHIC-BES energies: The roles of spin hall effect, initial condition, and baryon diffusion_, _Phys. Rev. C_ 105 (2022) 064909 [2204.02218].
  * [23] C.G. de Oliveira and J. Tiomno, _Representations of Dirac equation in general relativity_ , _Nuovo Cim._ 24 (1962) 672.
  * [24] F.W. Hehl and W.-T. Ni, _Inertial effects of a Dirac particle_ , _Phys. Rev. D_ 42 (1990) 2045.
  * [25] M. Buzzegoli and D.E. Kharzeev, _Anomalous gravitomagnetic moment and nonuniversality of the axial vortical effect at finite temperature_ , _Phys. Rev. D_ 103 (2021) 116005 [2102.01676].
  * [26] I.Y. Kobzarev and L.B. Okun, _Gravitational interaction of fermions_ , _Zh. Eksp. Teor. Fiz._ 43 (1962) 1904.
  * [27] C.F. Cho and N.D. Hari Dass, _Equivalence Principle, Stress Tensor and Long Range Behavior of Gravitational Interactions_ , _Phys. Rev. D_ 14 (1976) 2511.
  * [28] O.V. Teryaev, _Gravitational form factors and nucleon spin structure_ , _Front. Phys. (Beijing)_ 11 (2016) 111207.
  * [29] J.F. Donoghue, B.R. Holstein and R.W. Robinett, _Renormalization of the Energy Momentum Tensor and the Validity of the Equivalence Principle at Finite Temperature_ , _Phys. Rev. D_ 30 (1984) 2561.
  * [30] J.F. Donoghue, B.R. Holstein and R.W. Robinett, _The Principle of Equivalence at Finite Temperature_ , _Gen. Rel. Grav._ 17 (1985) 207.
  * [31] S. Lin and T. Jia-Yuan, _Medium correction to gravitational form factors_ , _Acta Phys. Sin._ 72 (2023) 071201 [2302.12450].
  * [32] F. Du, L.E. Finch and J. Sandweiss, _Observing spontaneous strong CP violation through hyperon helicity correlations_ , _Phys. Rev. C_ 78 (2008) 044908.
  * [33] F. Becattini, M. Buzzegoli, A. Palermo and G. Prokhorov, _Polarization as a signature of local parity violation in hot QCD matter_ , _Phys. Lett. B_ 822 (2021) 136706 [2009.13449].
  * [34] S.K. Singh and J.-e. Alam, _Suppression of spin polarization as an indicator of QCD critical point_ , _Eur. Phys. J. C_ 83 (2023) 585 [2110.15604].
  * [35] W.M. Serenone, J.a.G.P. Barbon, D.D. Chinellato, M.A. Lisa, C. Shen, J. Takahashi et al., _$\Lambda$ polarization from thermalized jet energy_, _Phys. Lett. B_ 820 (2021) 136500 [2102.11919].
  * [36] V.H. Ribeiro, D. Dobrigkeit Chinellato, M.A. Lisa, W. Matioli Serenone, C. Shen, J. Takahashi et al., _$\Lambda$ polarization from vortex ring as medium response for jet thermalization_, 2305.02428.