# Phonon-mediated spin transport in quantum paraelectric metals

Kyoung-Min Kim<EMAIL_ADDRESS>Center for Theoretical Physics of Complex
Systems, Institute for Basic Science, Daejeon 34126, Republic of Korea  Suk
Bum Chung<EMAIL_ADDRESS>Department of Physics and Natural Science
Research Institute, University of Seoul, Seoul 02504, Republic of Korea
School of Physics, Korea Institute for Advanced Study, Seoul 02455, Republic
of Korea

###### Abstract

The concept of ferroelectricity is now often extended to include continuous
inversion symmetry-breaking transitions in various metals and doped
semiconductors. Paraelectric metals near ferroelectric quantum criticality,
which we term ‘quantum paraelectric metals,’ typically possess soft transverse
optical phonons that have Rashba-type coupling to itinerant electrons in the
presence of spin-orbit coupling. We find through the Kubo formula calculation
that such Rashba electron-phonon coupling has a profound impact on electron
spin transport. While the spin Hall effect arising from non-trivial electronic
band structures has been studied extensively, we find here the presence of the
Rashba electron-phonon coupling can give rise to spin current, including spin
Hall current, in response to an inhomogeneous electric field even with a
completely trivial band structure. Furthermore, this spin conductivity
displays unconventional characteristics, such as quadrupolar symmetry
associated with the wave vector of the electric field and a thermal activation
behavior characterized by scaling laws dependent on the phonon frequency to
temperature ratio. These findings shed light on exotic electronic transport
phenomena originating from ferroelectric quantum criticality, highlighting the
intricate interplay of charge and spin degrees of freedom.

spin conductivity, ferroelectric quantum criticality, Rashba spin-orbit
coupling, transverse optical phonons, Kubo formula

###### pacs:

72.10.Di, 72.25.Bq, 77.90.+k

##  I Introduction

Figure 1: Spin conductivity in quantum paraelectric metals. Schematic phase
diagram of quantum paraelectric metals with spin-orbit coupling near the
ferroelectric quantum critical point $(\delta,T)=(0,0)$ ($\delta$: tuning
parameter, $T$: temperature) and distinct scaling laws of spin conductivity
($\sigma_{ij}^{\alpha}$) in each phase, where $\alpha,i,j$ denote the
direction of spin quantization axis, spin current and electric fields,
respectively. In these quantum paraelectric metals, an inhomogeneous electric
field induces a spin current. In the paraelectric phase (yellow area), the
phonon-mediated spin conductivity possesses a thermal-activated form:
$\sigma_{ij}^{\alpha}\sim\exp{(-E_{g}/k_{B}T)}$, where $E_{g}\sim\delta$ is
the energy gap of soft transverse optical (TO) phonons associated with the
phase transition. In contrast, within the gapless-phonon region (red area),
the spin conductivity adopts a power-law behavior: $\sigma_{ij}^{\alpha}\sim
T^{2}$. This power law of spin conductivity may undergo modifications in the
quantum critical (QC) region (blue area) below a specific crossover
temperature scale $T_{*}$ (indicated by the dashed line), wherein electrons
exhibit non-Fermi liquid characteristics.

The search for materials combining ferroelectricity/polarity with metallicity
has been longstanding in condensed matter physics, dating back to the first
proposal by Anderson and Blount over 50 years ago [1]. This endeavor has made
significant progress, especially in the last decade, resulting in the
accumulation of numerous experimentally confirmed examples [2, 3, 4], starting
with LiOsO3 [5]. Other noteworthy examples include doped quantum paraelectrics
such as SrTiO3 [6, 7, 8], IV-VI compounds [9, 10], and specific bilayer
transition metal dichalcogenides [11, 12, 13]. These so-called ferroelectric
(or polar) metals, typically doped ferroelectrics in semimetals and
semiconductors, present the intriguing coexistence of ferroelectricity and
metallicity, contrary to their apparent mutual exclusivity. In addition, the
possibility of various correlated electronic phenomena arising from
ferroelectric quantum fluctuations near a ferroelectric quantum critical
point, including the augmentation of the critical temperature for
superconductivity, has attracted strong interest [14, 15, 6, 16, 17, 18, 19,
20].

For the displacive ferroelectrics under consideration, the continuous
ferroelectric phase transition involves the softening of transverse optical
(TO) phonon modes associated with the displacement in proximity to the
critical point [21, 22], as this transition is characterized by a collective
displacement of ions from their centrosymmetric positions [7, 19, 23]. Given
that the TO mode displacement breaks the inversion symmetry while preserving
the time-reversal symmetry, the interactions between the TO phonons and
itinerant electrons in the presence of any finite atomic spin-orbit coupling
takes the unconventional form of a Rashba-type spin-orbit coupling, which
couples the momentum and spin of itinerant electrons [24, 25, 26, 27]. We
refer to these distinctive interactions as “phonon-mediated spin-orbit
coupling” (PM-SOC). Previous theoretical studies explored the impacts of the
PM-SOC on correlated electronic phenomena in the quantum critical region, such
as non-Fermi liquid behavior [27] as well as augmented superconducting
instability [27, 20]. However, the effect of the PM-SOC on electronic
transport remains unexplored so far, remaining a missing piece of the physics
near the ferroelectric quantum critical region.

In this study, we investigate the influence of the PM-SOC on the electronic
transport properties of a centrosymmetric metal (i.e. possessing finite
carrier concentration) near the ferroelectric quantum critical point, which
may be termed as a ‘quantum paraelectric metal.’ From the Kubo formula, we
obtain a nonzero spin conductivity, even in the centrosymmetric paraelectric
phase (Fig. 1), from a single orbital. This phenomenon may seem counter-
intuitive at first glance since not only is the Rashba spin-orbit coupling in
the electronic band structure, which results in the finite spin Hall
conductivity [28, 29], symmetry-forbidden in such a phase, but any orbital
Hall effect [30, 31, 32, 33] is also absent. However, the Rashba-type spin-
orbit coupling to TO phonons [34], i.e., the PM-SOC, in conjunction with
inhomogeneous external electric fields, gives rise to an unconventional type
of spin conductivity. Notably, we shall show in Sec. IV that this phonon-
mediated spin conductivity exhibits a unique directional dependence on the
wave vector of external electric fields, displaying a quadrupolar symmetry
with respect to the wave vector that is, however, distinct from the
quadrupolar symmetry predicted for electrical Hall resistivity in quantum Hall
states [35, 36] or spin Hall conductivity in Rashba metals [37]. Furthermore,
we demonstrate in Sec. V that our phonon-mediated spin conductivity also
exhibits peculiar scaling laws as a function of a tuning parameter and
temperature (Fig. 1). Whereas most theoretical research on spin transport has
been based on band structure considerations, our findings point to new
possibilities in interaction-induced spin transport. Moreover, our results
highlight the intriguing aspect of the emergent exotic transport phenomena
arising from the intricate interplay of charge and spin degrees of freedom in
itinerant electrons in the realm of ferroelectric quantum criticality.

##  II Model

Quantum paraelectric metals are characterized by the emergence of soft TO
phonon modes and their distinctive electron-phonon interactions, which, in
combination with atomic spin-orbit coupling, exhibit a Rashba-type spin-orbit
coupling for itinerant electrons [32, 18, 25]. The minimal model for a quantum
paraelectric metal is given by the following effective Hamiltonian [26, 27]:

$\displaystyle\hat{H}=\sum_{s,\bm{p}}\xi_{\bm{p}}c_{s,\bm{p}}^{\dagger}c_{s,\bm{p}}+\frac{g}{\sqrt{V}}\sum_{\bm{k}}\sum_{s,s^{\prime},\bm{p}}\bm{\phi}_{\bm{k}}\cdot(\bm{\sigma}_{ss^{\prime}}\times\bm{p})c_{s,\bm{p}+\bm{k}/2}^{\dagger}c_{s^{\prime},\bm{p}-\bm{k}/2}$
(1)

for electrons, where $c_{s,\bm{p}}^{\dagger}$ is the electron creation
operator ($s$, $\bm{p}$ denoting the spin and wave vector of the electron,
respectively), the electron energy dispersion is isotropic, i.e.
$\xi_{\bm{p}}=\frac{\hbar^{2}|\bm{p}|^{2}}{2m}-\mu$ ($m$ and $\mu$ denoting
the electron effective mass and the chemical potential, respectively), $g$ is
the coupling constant for the electron-phonon interaction, and
$\bm{\phi}_{\bm{k}}$ is the transverse phonon displacement field, whose
dynamics, in the free limit, is given by the action

$S_{ph}=\frac{M_{ph}}{\hbar^{2}}\sum_{i\nu_{n}}\sum_{\bm{k}}\sum_{ij}\phi_{i}(-i\nu_{m},-\bm{k})(\nu_{n}^{2}+\omega^{2}_{\bm{k}})\mathcal{P}_{ij}(\bm{k})\phi_{j}(i\nu_{m},\bm{k}),$
(2)

where
$\mathcal{P}_{ij}(\bm{k})\equiv\delta_{ij}-(\hat{\bm{e}}_{i}\cdot\hat{\bm{k}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{k}})$
is the transverse projection operator, $\nu_{n}$ is the boson Matsubara
frequency, $M_{ph}$ is the phonon effective mass, and the phonon energy
dispersion is given by $\omega_{\bm{k}}^{2}=(\hbar c|\bm{k}|)^{2}+E_{g}^{2}$,
where $c$ and $E_{g}$ denote the phonon velocity and energy gap, respectively.

Figure 2: Relevant Feynman diagrams for the computation of phonon-mediated
spin conductivity. (a) The correlation function involving charge and spin
currents at the tree level. (b) The correlation function at the one-loop
order. Here, the black solid and blue wavy lines represent electron and phonon
propagators, respectively. The black vertex points denote the electron-phonon
interactions that result in Rashba-type spin-orbit coupling for electrons. The
red lines represent external electric fields. In each diagram, $\bm{p}$,
$\bm{k}$, and $\bm{q}$ represent the wave vectors for electrons, phonons, and
external electric fields, respectively.

Whereas previous studies of the Eq. (1) quantum paraelectric metal model
focused on its instability to superconductivity [26, 27], we calculate in this
work its spin conductivity, denoted as $\sigma_{ij}^{\alpha}(\omega,\bm{q})$,
using the Kubo formula:

$\displaystyle\sigma_{ij}^{\alpha}(\omega,\bm{q})=\frac{i}{\omega}\pi_{ij}^{\alpha}(i\omega_{n},\bm{q})\big{|}_{i\omega_{n}\rightarrow\omega+i\delta}.$
(3)

Here, the indices $i$ and $j$ denote the directions of the spin and charge
currents, respectively, while $\omega$ and $\mathbf{q}$ denote the frequency
and wave vector of the external electric field.
$\pi_{ij}^{\alpha}(i\omega_{n},\bm{q})$ represents the current-current
correlation function, defined as

$\displaystyle\pi_{ij}^{\alpha}(i\omega_{n},\bm{q})=$
$\displaystyle\;-\frac{1}{V}\int^{\beta}_{0}d\tau
e^{i\omega_{n}\tau}\big{\langle}j_{\textrm{spin},i}^{\alpha}(\tau,\bm{q})j_{j}(0,-\bm{q})\big{\rangle},$
(4)

where $j_{\textrm{spin},i}^{\alpha}(\bm{q})$ and $j_{j}(\bm{q})$ denote the
charge and spin current operators, and $\omega_{n}$ the fermion Matsubara
frequency. The average $\langle\cdots\rangle$ denotes the ensemble average
over the quantum partition function. The charge and spin current operators are
explicitly given by

$\displaystyle j_{\textrm{spin},i}^{\alpha}(\bm{q})=$
$\displaystyle\;\frac{\hbar^{2}}{m}\sum_{s,s^{\prime},\bm{p}}\bigg{[}\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\frac{1}{2}\bm{q}\bigg{)}\bigg{]}(\sigma_{\alpha})_{ss^{\prime}}c_{s,\bm{p}+\bm{q}}^{\dagger}c_{s^{\prime},\bm{p}},$
(5) $\displaystyle j_{j}(\bm{q})=$
$\displaystyle\;\frac{e\hbar}{m}\sum_{s,\bm{p}}\bigg{[}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}+\frac{1}{2}\bm{q}\bigg{)}\bigg{]}c_{s,\bm{p}+\bm{q}}^{\dagger}c_{s,\bm{p}}.$
(6)

Here, the unit vectors $\hat{\bm{e}}_{i}$ and $\hat{\bm{e}}_{j}$ denote the
directions of the spin and charge currents, respectively. We compute
$\sigma_{ij}^{\alpha}(\omega,\bm{q})$ through a diagrammatic expansion in $g$.
As a result, we obtain the following expression for
$\sigma^{\alpha}_{ij}(\omega,\bm{q})$ [38]:

$\displaystyle\sigma^{\alpha}_{ij}(\omega,\bm{q})=$
$\displaystyle\;\frac{e\hbar^{4}}{m^{2}}\int\frac{d^{3}p}{(2\pi)^{3}}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}\int^{\infty}_{-\infty}\frac{d\epsilon}{2\pi}\frac{1}{\omega}\bigg{[}n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\bigg{]}$
$\displaystyle\times
G_{\textrm{ret}}(\epsilon+\omega,\bm{p}+\bm{q})G_{\textrm{adv}}(\epsilon,\bm{p})\gamma_{i}^{\alpha}(\epsilon+\omega+i\delta,\epsilon-i\delta;\bm{p}+\bm{q},\bm{p}).$
(7)

Here, $n_{F}(\epsilon)$ denotes the Fermi-Dirac distribution function and
$G_{\textrm{ret}}(\epsilon+\omega,\bm{p}+\bm{q})$,
$G_{\textrm{adv}}(\epsilon,\bm{p})$ the full retarded and advanced propagators
for the electron and hole states, respectively. We first focus on the
paraelectric phase near the critical point, specifically, away from the
quantum critical region. In this case, we may posit that the electron
propagators possess well-defined quasiparticle peaks:
$G_{\textrm{ret}}(\epsilon+\omega,\bm{p}+\bm{q})=(\epsilon+\omega-\xi_{\bm{p}+\bm{q}}+i\hbar/2\tau)^{-1}$
and
$G_{\textrm{adv}}(\epsilon,\bm{p})=(\epsilon-\xi_{\bm{p}}-i\hbar/2\tau)^{-1}$,
where $\tau$ is the electron scattering time.

##  III Phonon-mediated spin conductivity

Figure 3: Phonon-mediated spin polarization. (a) Schematic illustration of the
vector calculations involved in computing phonon-mediated spin polarization.
The outcome, $\hat{\bm{e}}_{\alpha}\cdot\bm{P}_{\perp}$, represents the value
of the desired spin polarization. Here, the unit vectors
$\hat{\bm{e}}_{\alpha}$ and $\hat{\bm{k}}$ denote the directions of the spin
polarization and the phonon wave vector, respectively.
$\bm{P}=\bm{q}\times(\bm{p}+\bm{k}/2)$ denotes the vector product of $\bm{q}$
and $\bm{p}+\bm{k}/2$. (b-d) Plots of
$\hat{\bm{e}}_{\alpha}\cdot\bm{P}_{\perp}$ as a function of the azimuthal
($\phi_{k}$) and polar ($\theta_{k}$) angles of
$\hat{\bm{k}}=(\sin{\theta_{k}}\cos{\phi_{k}},\sin{\theta_{k}}\sin{\phi_{k}},\cos{\theta_{k}})$.
The color scale denotes the value of
$\hat{\bm{e}}_{\alpha}\cdot\bm{P}_{\perp}$. Specific relative orientations
between $\hat{\bm{e}}_{\alpha}$ and $\bm{P}$ are utilized in each plot. (b)
Plot for $\bm{P}\parallel\hat{\bm{e}}_{\alpha}$, where
$\hat{\bm{e}}_{\alpha}=\hat{z}$ and $\bm{P}=|\bm{P}|\hat{z}$. (c) Plot for
$\hat{\bm{e}}_{\alpha}=\hat{z}$ and
$\bm{P}=|\bm{P}|(\sqrt{2}/2,0,\sqrt{2}/2)$. (d) Plot for
$\bm{P}\perp\hat{\bm{e}}_{\alpha}$, where $\hat{\bm{e}}_{\alpha}=\hat{z}$ and
$\bm{P}=|\bm{P}|\hat{x}$. In (b-d), the cyan and magenta lines indicate the
conditions under which $\hat{\bm{e}}_{\alpha}\cdot\bm{P}_{\perp}$ vanishes,
corresponding to $\hat{\bm{k}}||\pm\hat{\bm{e}}_{\alpha}$ and
$\hat{\bm{e}}_{\alpha}\perp\bm{P}_{\perp}$, respectively.

From the yet unelucidated vertex term $\gamma^{\alpha}_{i}$, we can show that
the spin conductivity of Eq. (II) vanishes in lieu of the virtual TO phonon
exchange. At the tree level as shown in Fig. 2(a), where such exchange is
absent, it arises entirely from $j_{\textrm{spin},i}^{\alpha}(\bm{q})$ and
takes the form of
$\gamma_{i}^{\alpha}(\epsilon-i\delta,\epsilon+\omega+i\delta;\bm{p},\bm{p}+\bm{q})=\hat{\bm{e}}_{i}\cdot(\bm{p}+\frac{1}{2}\bm{q})\textrm{tr}(\sigma_{\alpha})$,
where “tr” represents a trace over spin matrices, and hence vanishes. However,
at the one-loop order, as shown in Fig. 2(b), the PM-SOC within this system
allows for a non-vanishing form of the vertex function:

$\displaystyle\gamma_{i}^{\alpha}(\epsilon-i\delta,\epsilon+\omega+i\delta;\bm{p},\bm{p}+\bm{q})=$
$\displaystyle\;\frac{g^{2}\hbar^{2}}{M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})\mathcal{S}(\epsilon,\bm{p},\bm{k};\omega,\bm{q}).$
(8)

Here, $\mathcal{S}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$ accounts for the
creation of a virtual phonon with a wave vector $\bm{k}$, leading to the
transition of the electron and hole states into $(\bm{p}+\bm{q}+\bm{k})$ and
$(\bm{p}+\bm{k})$, respectively. Such a factor is a common occurrence in
conventional electron-phonon systems [38].
$\mathcal{S}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$ is explicitly given by

$\displaystyle\mathcal{S}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})=$
$\displaystyle\;\bigg{[}n_{B}(\omega_{\bm{k}})+\frac{1}{2}n_{F}(\epsilon+\omega_{\bm{k}})+\frac{1}{2}n_{F}(\epsilon+\omega+\omega_{\bm{k}})\bigg{]}G_{\textrm{adv}}(\epsilon+\omega_{\bm{k}},\bm{p}+\bm{k})G_{\textrm{ret}}(\epsilon+\omega+\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})$
$\displaystyle+\bigg{[}n_{B}(\omega_{\bm{k}})+1-\frac{1}{2}n_{F}(\epsilon-\omega_{\bm{k}})-\frac{1}{2}n_{F}(\epsilon+\omega-\omega_{\bm{k}})\bigg{]}G_{\textrm{adv}}(\epsilon-\omega_{\bm{k}},\bm{p}+\bm{k})G_{\textrm{ret}}(\epsilon+\omega-\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q}).$
(9)

$\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$ now includes the effects of
the PM-SOC vertex pair as shown in Fig. 2(b):

$\displaystyle\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})=\sum_{k,l}\mathcal{P}_{kl}(\bm{k})\textrm{tr}\Big{[}\sigma_{\alpha}\hat{\bm{e}}_{k}\cdot(\bm{\sigma}\times(\bm{p}+\bm{q}+\bm{k}/2))\hat{\bm{e}}_{l}\cdot(\bm{\sigma}\times(\bm{p}+\bm{k}/2))\Big{]}.$
(10)

Parts of this trace arise from the annihilation and the creation of a virtual
TO phonon through the Rashba-type spin-orbit interactions, which result in
$\hat{\bm{e}}_{k}\cdot(\bm{\sigma}\times(\bm{p}+\bm{q}+\bm{k}/2))$ and
$\hat{\bm{e}}_{l}\cdot(\bm{\sigma}\times(\bm{p}+\bm{k}/2))$ respectively.
Their product can be written as:

$\displaystyle\hat{\bm{e}}_{k}\cdot(\bm{\sigma}\times(\bm{p}+\bm{q}+\bm{k}/2))\hat{\bm{e}}_{l}\cdot(\bm{\sigma}\times(\bm{p}+\bm{k}/2))$
$\displaystyle=$
$\displaystyle\;((\bm{p}+\bm{q}+\bm{k}/2)\cdot(\bm{p}+\bm{k}/2))(\hat{\bm{e}}_{k}\cdot\hat{\bm{e}}_{l})-((\bm{p}+\bm{q}+\bm{k}/2)\cdot\hat{\bm{e}}_{l})((\bm{p}+\bm{k}/2)\cdot\hat{\bm{e}}_{k})$
$\displaystyle+i(\bm{p}+\bm{q}+\bm{k}/2)\cdot(\hat{\bm{e}}_{k}\times\hat{\bm{e}}_{l})((\bm{p}+\bm{k}/2)\cdot\bm{\sigma})+i\hat{\bm{e}}_{k}\cdot(\bm{q}\times(\bm{p}+\bm{k}/2))(\hat{\bm{e}}_{l}\cdot\bm{\sigma}).$
(11)

In this expression, the first two terms vanish after taking the trace over
spin due to $\sigma_{\alpha}$. The third term vanishes after averaging over
all phonon polarizations ($k$ and $l$) due to antisymmetry in
$k\leftrightarrow l$, as the transverse projection operator
$\mathcal{P}_{kl}(\bm{k})$ is symmetric. However, the last term may survive in
both tracing over spin and summing over the phonon polarizations, leaving a
nonzero spin polarization,

$\displaystyle\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})=2i\Big{\\{}\hat{\bm{e}}_{\alpha}\cdot\big{[}\bm{q}\times(\bm{p}+\bm{k}/2)\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\hat{\bm{k}}\cdot\big{[}\bm{q}\times(\bm{p}+\bm{k}/2)\big{]}\Big{\\}},$
(12)

in contrast to the tree level, where $\hat{\bm{e}}_{\alpha}$ denotes the unit
vector in the direction of $\alpha$. This quantity is termed “phonon-mediated
spin polarization” (PMSP).

A simple picture of the constraint on PMSP can be obtained through a simple
re-writing of our PMSP. The first step is to take the vector product of the
two momenta associated with the electron and hole states, namely
$(\bm{p}+\bm{q}+\bm{k}/2)$ and $(\bm{p}+\bm{k}/2)$:
$\bm{P}=\bm{q}\times(\bm{p}+\bm{k}/2)$. The next step is to subtract from
$\bm{P}$ its $\hat{\bm{k}}$ projection, resulting in the vector
$\bm{P}_{\perp}=\bm{P}-\hat{\bm{k}}(\hat{\bm{k}}\cdot\bm{P})$, as illustrated
in Fig. 3(a). The final step is to take the scalar product of $\bm{P}_{\perp}$
and $\hat{\bm{e}}_{\alpha}$, resulting in
$\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})=2i\hat{\bm{e}}_{\alpha}\cdot\bm{P}_{\perp}$.
We can now easily see how $\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$
vanishes for specific orientations of $\hat{\bm{k}}$:

$\displaystyle\hat{\bm{k}}\parallel\pm\hat{\bm{e}}_{\alpha}\Rightarrow\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})=0,$
(13)

as depicted by the cyan lines in Fig. 3(b–d). In other words, TO phonons with
momenta parallel to $\hat{\bm{e}}_{\alpha}$ do not contribute to the PMSP and,
ultimately, to spin conductivity. Additionally, there exists a more general
vanishing condition:

$\displaystyle\hat{\bm{e}}_{\alpha}\perp\bm{P}_{\perp}\Rightarrow\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})=0,$
(14)

as illustrated by the magenta lines in Fig. 3(b–d). Consequently, TO phonons
that result in $\bm{P}_{\perp}$ being perpendicular to $\hat{\bm{e}}_{\alpha}$
also do not contribute to the generation of the PMSP.

From the PMSP we have obtained, we find that the lowest order term in ${\bf
q}$ of the DC spin conductivity, i.e.

$\sigma^{\alpha}_{ij}(\bm{q})=\lim_{\omega\rightarrow
0}\sigma^{\alpha}_{ij}(\omega,\bm{q}),$

is quadratic, which means that the spin current will arise in response to an
inhomogeneous electric field. By substituting Eqs. (8), (III) and (10) into
Eq. (II) and taking the limit $\omega\rightarrow 0$, we obtain

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}$
$\displaystyle\times\hat{\bm{e}}_{\alpha}\cdot\Big{\\{}[\bm{q}\times(\bm{p}+\bm{k}/2)]-\hat{\bm{k}}\big{(}\hat{\bm{k}}\cdot[\bm{q}\times(\bm{p}+\bm{k}/2)]\big{)}\Big{\\}}$
$\displaystyle\times\Big{[}G_{\textrm{ret}}(0,\bm{p}+\bm{q})G_{\textrm{adv}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})G_{\textrm{adv}}(\omega_{\bm{k}},\bm{p}+\bm{k})+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}})\Big{]}.$
(15)

Given that our PMSP factor is proportional to $|\bm{q}|$,
$\sigma^{\alpha}_{ij}(\bm{q})$ obviously vanishes as $|\bm{q}|\rightarrow 0$.
However, $\sigma^{\alpha}_{ij}(\bm{q})$ turns out to be even, not odd, in
$\bm{q}$. Given that the electron dispersion is even in momentum, i.e.
$\xi_{-\bm{p}}=\xi_{\bm{p}}$, we expect the electron and the hole propagators
to be even in their momenta. Since the phonon dispersion is likewise even in
momentum, i.e. $\omega_{-\bm{k}}=\omega_{\bm{k}}$, it is straightforward to
show that Eq. (III) is unchanged by reversing the sign of all its momenta,
i.e. $\bm{p}\leftrightarrow-\bm{p}$, $\bm{k}\leftrightarrow-\bm{k}$ and
$\bm{q}\leftrightarrow-\bm{q}$. Consequently, if we expand
$\sigma^{\alpha}_{ij}(\bm{q})$ in the powers of $\bm{q}$, the first non-
vanishing term is quadratic, hence the predominant $\bm{q}$-dependence of
$\sigma^{\alpha}_{ij}(\bm{q})$. Our numerical integration on Eq. (III)
confirms this quadratic behavior within a low-$|\bm{q}|$ regime (Fig. 4(b)).
Assuming $E_{F}=1$ eV and $\tau=10^{-14}$ s, we determined that the quartic
term is significantly smaller than the quadratic term for
$|\bm{q}|\ll\frac{\hbar}{E_{F}\tau}k_{F}\approx 0.0212k_{F}$, where $k_{F}$
denotes the Fermi wave vector. Further details can be found in the Appendix C.

$\sigma_{yy}^{x}(\bm{q})=\chi_{0}(-2q_{z}q_{y})$
---
$\sigma_{zz}^{x}(\bm{q})=\chi_{0}(2q_{y}q_{z})$
$\sigma_{xx}^{y}(\bm{q})=\chi_{0}(2q_{z}q_{x})$
$\sigma_{zz}^{y}(\bm{q})=\chi_{0}(-2q_{x}q_{z})$
$\sigma_{xx}^{z}(\bm{q})=\chi_{0}(-2q_{y}q_{x})$
$\sigma_{yy}^{z}(\bm{q})=\chi_{0}(2q_{x}q_{y})$
Table 1: Nonvanishing components of longitudinal conductivity.
$\sigma_{ii}^{\alpha}(\bm{q})$ for $i,\alpha=\\{x,y,z\\}$ denotes the
longitudinal conductivity for spin currents, as defined in Eq. (17).
$\bm{q}=(q_{x},q_{y},q_{z})$ denotes a wave vector of external electric
fields. $\chi_{0}$ is a $\bm{q}$-independent constant, as defined in Eq. (19).
$\sigma_{xy}^{x}(\bm{q})=[\sigma_{yx}^{x}(\bm{q})]^{*}=\chi_{0}(-q_{x}q_{z})+i\kappa_{0}(-q_{x}q_{z})$
---
$\sigma_{yz}^{x}(\bm{q})=[\sigma_{zy}^{x}(\bm{q})]^{*}=\chi_{0}(q_{y}^{2}-q_{z}^{2})+i\kappa_{0}(q_{y}^{2}+q_{z}^{2})$
$\sigma_{zx}^{x}(\bm{q})=[\sigma_{xz}^{x}(\bm{q})]^{*}=\chi_{0}(q_{x}q_{y})+i\kappa_{0}(-q_{x}q_{y})$
$\sigma_{xy}^{y}(\bm{q})=[\sigma_{yx}^{y}(\bm{q})]^{*}=\chi_{0}(q_{y}q_{z})+i\kappa_{0}(-q_{y}q_{z})$
$\sigma_{yz}^{y}(\bm{q})=[\sigma_{zy}^{y}(\bm{q})]^{*}=\chi_{0}(-q_{x}q_{y})+i\kappa_{0}(-q_{x}q_{y})$
$\sigma_{zx}^{y}(\bm{q})=[\sigma_{xz}^{y}(\bm{q})]^{*}=\chi_{0}(q_{z}^{2}-q_{x}^{2})+i\kappa_{0}(q_{z}^{2}+q_{x}^{2})$
$\sigma_{xy}^{z}(\bm{q})=[\sigma_{yx}^{z}(\bm{q})]^{*}=\chi_{0}(q_{x}^{2}-q_{y}^{2})+i\kappa_{0}(q_{x}^{2}+q_{y}^{2})$
$\sigma_{yz}^{z}(\bm{q})=[\sigma_{zy}^{z}(\bm{q})]^{*}=\chi_{0}(q_{x}q_{z})+i\kappa_{0}(-q_{x}q_{z})$
$\sigma_{zx}^{z}(\bm{q})=[\sigma_{xz}^{z}(\bm{q})]^{*}=\chi_{0}(-q_{y}q_{z})+i\kappa_{0}(-q_{y}q_{z})$
Table 2: Nonvanishing components of transverse conductivity.
$\sigma_{ij}^{\alpha}(\bm{q})$ for $i,j,\alpha=\\{x,y,z\\}$ denotes the spin
conductivity, as defined in Eq. (17). $\bm{q}=(q_{x},q_{y},q_{z})$ denotes a
wave vector of external electric fields. $\chi_{0}$ and $\kappa_{0}$ are
$\bm{q}$-independent constants, as defined in Eqs. (19) and (20),
respectively. $*$ denotes complex conjugation.

To obtain this first nonvanishing term of $\sigma^{\alpha}_{ij}(\bm{q})$,

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})\simeq$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}$
$\displaystyle\times
p_{i}^{2}\big{[}-\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}(1-\hat{k}_{\alpha}^{2})|G_{\textrm{ret}}(0,\bm{p})|^{4}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}$
$\displaystyle\times
p_{i}p_{j}\big{[}\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})\big{]}(1-\hat{k}_{\alpha}^{2})\frac{4\hbar^{2}}{m}(\bm{p}\cdot\bm{q})|G_{\textrm{ret}}(0,\bm{p})|^{4}G_{\textrm{ret}}(0,\bm{p})+\mathcal{O}(q^{4}),$
(16)

we need to consider how the electron propagator should be expanded in the
powers of $\bm{q}$ and how the dependence on $\bm{k}$ should be approximated.
In this equation, the first term is derived from Eq. (III) by substituting
$\bm{q}=0$ into the electron propagators and collecting nonvanishing even
terms in $\bm{q}$. This term contributes to the imaginary part of the spin
conductivity as the integrand is purely real. But there is also the second
term arising from expanding the electron propagators in $\bm{q}$. This
expansion leads to the additional factor
$\frac{4\hbar^{2}\bm{p}\cdot\bm{q}}{m}G_{\textrm{ret}}(0,\bm{p})$.
Conventionally, this term is disregarded because it is of higher order in
$\bm{q}$. However, in our case, both this term and the first term are
quadratic in $\bm{q}$, and hence, both must be retained. The second term
contributes to the real part of the spin conductivity as the additional
propagator $G_{\textrm{ret}}(0,\bm{p})$ in the integrand introduces a purely
imaginary factor $\big{(}-i\frac{E_{F}\tau}{\hbar}\big{)}$. It is noteworthy
that this real part predominates in typical metals as
$\frac{E_{F}\tau}{\hbar}\gg 1$ is required for having well-defined
quasiparticles; conversely, close to the ferroelectric quantum criticality,
the predominance of the real part cannot be taken for granted. It needs to be
noted here that, in deriving Eq. (III), the phonon momentum $\bm{k}$ is set to
be zero except for
$[n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})]/\omega_{\bm{k}}$ that
effectively constrains $|\bm{k}|$ at small values $|\bm{k}|\ll|\bm{p}|\sim
k_{F}$ in the $\bm{k}$-integral [38]. Further details can be found in Methods.

The lowest-order $\bm{q}$-dependence and the temperature dependence of
$\sigma^{\alpha}_{ij}(\bm{q})$ can be obtained analytically by performing the
Eq. (III) integration, details of which can be found in the Methods. The
result can be written in the form

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=\chi_{0}\big{(}\nu_{ij}^{\alpha}(\bm{q})+\nu_{ji}^{\alpha}(\bm{q})\big{)}+i\kappa_{0}\big{(}-\nu_{ij}^{\alpha}(\bm{q})+\nu_{ji}^{\alpha}(\bm{q})\big{)},$
(17)

where $\nu_{ij}^{\alpha}(\bm{q})$ is a quadratic function of $\bm{q}$, defined
as:

$\displaystyle\nu_{ij}^{\alpha}(\bm{q})=\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j}).$
(18)

$\chi_{0}$ and $\kappa_{0}$ are $\bm{q}$-independent constants, which are
defined as

$\displaystyle\chi_{0}$
$\displaystyle\equiv\bigg{(}\frac{ne\hbar\tau}{m}\bigg{)}\bigg{(}\frac{g^{2}}{\hbar
M_{ph}c^{3}}\bigg{)}\bigg{(}\frac{k_{B}T\tau}{\hbar}\bigg{)}^{2}\bigg{(}\frac{E_{F}\tau}{\hbar}\bigg{)}\frac{8}{5\pi^{2}}I\bigg{(}\frac{E_{g}}{k_{B}T}\bigg{)},$
(19) $\displaystyle\kappa_{0}$
$\displaystyle\equiv\bigg{(}\frac{ne\hbar\tau}{m}\bigg{)}\bigg{(}\frac{g^{2}}{\hbar
M_{ph}c^{3}}\bigg{)}\bigg{(}\frac{k_{B}T\tau}{\hbar}\bigg{)}^{2}\frac{2}{3\pi^{2}}I\bigg{(}\frac{E_{g}}{k_{B}T}\bigg{)};$
(20)

which confirms $\chi_{0}/\kappa_{0}=4E_{F}\tau/15\hbar\gg 1$ when the
quasiparticles are well-defined. In this expression,
$n=\frac{k_{F}^{3}}{3\pi^{2}}$ and $E_{F}=\frac{\hbar^{2}k_{F}^{2}}{2m}$
denote the electron density and the Fermi energy, respectively.
$I\big{(}\frac{E_{g}}{k_{B}T}\big{)}$ denotes an integral function
representing the integral over $|\bm{k}|$, which is given by:

$\displaystyle
I(y)=\int^{\infty}_{0}dx\frac{x^{2}}{2\sqrt{x^{2}+y^{2}}\sinh{\big{[}\sqrt{x^{2}+y^{2}}\;\big{]}}}\approx\frac{\pi^{2}}{8}\exp\big{(}-ay^{b}\big{)},$
(21)

where $a\approx 0.75$ and $b\approx 1.05$.

##  IV Quadrupolar symmetry

Figure 4: Quadrupole symmetry of the spin conductivity. (a) Real part of the
spin conductivity ($\textrm{Re}[\sigma_{xy}^{z}(\bm{q})]$) as a function of
$\bm{q}=(q_{x},q_{y},0)$ (in units of $q_{x}/k_{F}$ and $q_{y}/k_{F}$)
obtained from a numerical integration using the formula in Eq. (III). The
color scale indicates the magnitude of the spin conductivity divided by the
overall factor $\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}$. (b)
$\textrm{Re}[\sigma_{xy}^{z}(\bm{q})]$ along the line $q_{y}=0$, showcasing
the quadratic dependence on $q_{x}^{2}$. The parameter values used are as
follows: $E_{F}=3\textrm{eV}$, $v_{F}=10^{6}\textrm{m/s}$,
$\tau=10^{-14}\textrm{s}$, $c=10^{4}\textrm{m/s}$, $E_{g}=1\textrm{meV}$, and
$T=1\textrm{K}$.

The nonvanishing components of rank-3 tensor $\sigma_{ij}^{\alpha}(\bm{q})$
exhibit the unconventional quadrupolar symmetry in $\bm{q}$. In principle,
$\sigma_{ij}^{\alpha}(\bm{q})$ could have twenty-seven distinct components,
among various combinations of $(i,j,\alpha)$ where $i,j,\alpha=\\{x,y,z\\}$.
However, three longitudinal components with $i=j=\alpha$ vanish, remaining
only six components with $i=j\neq\alpha$. Moreover, we can see from Eq. (17)
that each pair of transverse components with $i\neq j$ is related by a
symmetric relation in the permutation $i\leftrightarrow j$:

$\displaystyle\sigma_{ij}^{\alpha}(\bm{q})=[\sigma_{ji}^{\alpha}(\bm{q})]^{*},$
(22)

where ∗ denotes a complex conjugation. The symmetry of Eq. (22) also reduces
the number of independent components to nine for the transverse components.
The six longitudinal and nine transverse components are summarized in Tabs. 1
and 2, respectively. These response functions collectively characterize the
spin conductivity in the quantum paraelectric metal. Among these nonvanishing
components, of particular interest are $\sigma_{xy}^{z}(\bm{q})$ and
$\sigma_{yx}^{z}(\bm{q})$. These transverse spin conductivity components are
explicitly given, in accordance with Eq. (22) as:

$\displaystyle\sigma_{xy}^{z}(\bm{q})$
$\displaystyle=\chi_{0}(q_{x}^{2}-q_{y}^{2})+i\kappa_{0}(q_{x}^{2}+q_{y}^{2}),$
(23) $\displaystyle\sigma_{yx}^{z}(\bm{q})$
$\displaystyle=\chi_{0}(q_{x}^{2}-q_{y}^{2})-i\kappa_{0}(q_{x}^{2}+q_{y}^{2}).$
(24)

Note here the quadrupolar dependence on $\bm{q}$ characterized by
$q_{x}^{2}-q_{y}^{2}$ for the real parts of these components. Our numerical
computation on Eq. (III) confirms this quadrupolar behavior within a
low-$|\bm{q}|$ regime (Fig. 4(a)).

Figure 5: Schematic illustration of spin currents generated by external
electric fields. (a) The longitudinal electric field
($\bm{E}_{\parallel}(\bm{r})$) and its wave vector ($\bm{q}$) align in the
same $x$ direction, as indicated by the black and green arrows, respectively.
A sinusoidal form
$\bm{E}_{\parallel}(\bm{r})=E_{0}\hat{x}\sin{(\bm{q}\cdot\bm{r})}$ is assumed
where translation symmetry is assumed in the other directions. The spin
current ($\bm{J}^{\alpha}_{\textrm{spin}}(\bm{r})$) and the spin quantization
axis ($\hat{\bm{e}}_{\alpha}$) point in the $y$ and $z$ directions,
respectively, as shown by the red and blue arrows. (b) The transverse electric
field ($\bm{E}_{\perp}(\bm{r})$) and its wave vector ($\bm{q}$) point in the
different $z$ and $x$ directions, as denoted by the black and green arrows,
respectively. A sinusoidal form
$\bm{E}_{\perp}(\bm{r})=E_{0}\hat{z}\sin{(\bm{q}\cdot\bm{r})}$ is assumed. The
spin current ($\bm{J}^{\alpha}_{\textrm{spin}}(\bm{r})$) and spin quantization
axis ($\hat{\bm{e}}_{\alpha}$) point in the $y$ and $z$ directions, as denoted
by the red and blue arrows, respectively. In each plot, the thin black line
guides the real-space modulation of $\bm{E}_{\parallel}(\bm{r})$ and
$\bm{E}_{\perp}(\bm{r})$ along the $x$ direction. The sizes of the black and
red arrows indicate the magnitudes of the electric field and the spin current,
respectively.

The phonon-mediated spin conductivity can also be written as the sum of very
distinct responses to the longitudinal and the transverse electric field,
respectively, allowing us to clarify its quadrupolar symmetry further. This
can be derived by using $\sigma_{ij}^{\alpha}(\bm{q})$ in Eq. (17) to
calculate spin currents generated by external electric fields, which are
expressed as:

$\displaystyle\bm{J}_{\textrm{spin}}^{\alpha}(\bm{q})=\sum_{i,j}\hat{\bm{e}}_{i}\sigma_{ij}^{\alpha}(\bm{q})E_{j}(\bm{q}).$
(25)

Here, $E_{j}(\bm{q})$ denotes an electric field applied in the $j$ direction.
Substituting Eq. (17) into Eq. (25), we obtain the first of our main results,
an explicit expression of $\bm{J}_{\textrm{spin}}^{\alpha}(\bm{q})$:

$\displaystyle\bm{J}_{\textrm{spin}}^{\alpha}(\bm{q})=(\chi_{0}-i\kappa_{0})(\hat{\bm{e}}_{\alpha}\times\bm{q})[\bm{q}\cdot\bm{E}_{\parallel}(\bm{q})]+(\chi_{0}+i\kappa_{0})\bm{q}[\hat{\bm{e}}_{\alpha}\cdot(\bm{q}\times\bm{E}_{\perp}(\bm{q}))].$
(26)

Here, $\bm{E}_{\perp}(\bm{q})$ and $\bm{E}_{\parallel}(\bm{q})$ represent
longitudinal and transverse components of $\bm{E}(\bm{q})$ with respect to
$\bm{q}$, respectively, which are defined as

$\displaystyle\bm{E}_{\parallel}(\bm{q})$
$\displaystyle=\frac{\bm{q}(\bm{q}\cdot\bm{E}(\bm{q}))}{|\bm{q}|^{2}},$ (27)
$\displaystyle\bm{E}_{\perp}(\bm{q})$
$\displaystyle=\bm{E}(\bm{q})-\frac{\bm{q}(\bm{q}\cdot\bm{E}(\bm{q}))}{|\bm{q}|^{2}}.$
(28)

We note that the spin current of Eq. (26) will be non-uniform in general, e.g.
when $\kappa_{0}/\chi_{0}$ can be ignored, it can be written in the real space
approximately as

$\bm{J}_{\textrm{spin}}^{\alpha}=-\chi_{0}\bigg{[}(\hat{\bm{e}}_{\alpha}\times\bm{\nabla})(\bm{\nabla}\cdot\bm{E})+\bm{\nabla}\\{\hat{\bm{e}}\cdot(\bm{\nabla}\times\bm{E})\\}\bigg{]}.$

Eq. (26) tells us that $\bm{E}_{\parallel}$ generates the spin Hall current
perpendicular to the spin polarization direction $\hat{\bm{e}}_{\alpha}$
whereas $\bm{E}_{\perp}$ generates the $\bm{q}$-parallel spin Hall current
whose magnitude is maximized for $\hat{\bm{e}}_{\alpha}\perp\bm{q}$. For an
example of response when $\bm{E}=\bm{E}_{\parallel}$, we can consider
$\bm{q}=q\hat{\bm{x}}$, $\bm{E}(\bm{q})=E_{0}\hat{\bm{x}}$, and
$\hat{\bm{e}}_{\alpha}=\hat{\bm{z}}$, for which the resulting spin Hall
current flows in the $y$ direction (Fig. 5(a)):

$\displaystyle\bm{J}_{\textrm{spin}}^{z}(\bm{q})=(\chi_{0}-i\kappa_{0})\hat{\bm{y}}q^{2}E_{0}.$
(29)

In fact, this spin current is of the form analogous to the charge current in
quantum Hall states induced by inhomogeneous electric field [35, 36] due to
the Hall viscosity [39, 40]; the spin analogue in the quantum spin Hall state
[41] and the Rashba metal [37] has also been discussed. For an example of
response when $\bm{E}=\bm{E}_{\perp}$, we can consider $\bm{q}=q\hat{\bm{x}}$,
$\bm{E}(\bm{q})=E_{0}\hat{\bm{z}}$, and $\hat{\bm{e}}_{\alpha}=\hat{\bm{y}}$,
for which the resulting spin current flows in the $x$ direction (Fig. 5(b)):

$\displaystyle\bm{J}_{\textrm{spin}}^{y}(\bm{q})=(\chi_{0}+i\kappa_{0})(-\hat{\bm{x}})q^{2}E_{0}.$
(30)

This $\bm{q}$-parallel Hall spin current induced by $\bm{E}_{\perp}$ has, to
the best of our knowledge, no analogue reported in the transport of either the
Rashba metal or the quantum Hall state and can be considered as the most
distinct feature of the phonon-mediated spin conductivity.

The above separation of the two components of the spin conductivity is
physically relevant as a metal exhibits different responses to
$\bm{E}_{\parallel}$ and $\bm{E}_{\perp}$. Briefly, $\bm{E}_{\parallel}$ is
screened below the plasma frequency, whereas, as dictated by the Faraday
effect, $\bm{E}_{\perp}$ is necessarily dynamic 111The coupling of electric
field to the TO phonon, omitted in our analysis, can be adequately treated
with the static dielectric constant for frequency far below $E_{g}/\hbar$ ..
Indeed, the response to $\bm{E}_{\perp}$ can be interpreted as the generation
of spin current in response to an electromagnetic wave if we additionally take
into account the frequency dependence of the spin conductivity. For instance,
we have derived in the Appendix E the approximate relation, valid for the low
frequency $\omega\tau\ll 1$, between the static and dynamic conductivity in
the first order in $\omega$:

$\displaystyle\sigma_{ij}^{\alpha}(\omega,\bm{q})=\bigg{(}1+\frac{3}{2}i\omega\tau\bigg{)}\sigma_{ij}^{\alpha}(\bm{q}).$
(31)

Physically, large differences between $|\bm{E}_{\parallel}|$ and
$|\bm{E}_{\perp}|$ is required for obtaining a nearly pure spin Hall current;
$|\bm{E}_{\parallel}(\bm{q})|\neq 0$ and $|\bm{E}_{\perp}(\bm{q})|\neq 0$, on
the other hands, will give us a nonzero longitudinal spin current. For
instance, with
$\bm{E}(\bm{q}=q\hat{\bm{x}})=E_{0}\frac{\hat{\bm{x}}+\hat{\bm{z}}}{\sqrt{2}}$
and $\hat{\bm{e}}_{\alpha}=\hat{\bm{y}}$, the spin Hall and longitudinal
currents are given by

$\displaystyle\bm{J}_{\textrm{sH}}^{y}(\bm{q}=q\hat{\bm{x}})$
$\displaystyle=-i\kappa_{0}\frac{\hat{\bm{x}}-\hat{\bm{z}}}{\sqrt{2}}q^{2}E_{0},$
(32) $\displaystyle\bm{J}_{\textrm{lon}}^{y}(\bm{q}=q\hat{\bm{x}})$
$\displaystyle=-\chi_{0}\frac{\hat{\bm{x}}+\hat{\bm{z}}}{\sqrt{2}}q^{2}E_{0}.$
(33)

Likewise, for any of the longitudinal spin conductivity components of Table 1
to be nonzero, both $|\bm{E}_{\parallel}|\neq 0$ and $|\bm{E}_{\perp}|\neq 0$
need to hold.

##  V Thermal activation behavior

Figure 6: Scaling laws in spin conductivity. (a) Real part of the spin
conductivity ($\textrm{Re}[\sigma_{xy}^{z}(\bm{q})]$) as a function of
temperature ($T$). (b) Imaginary part of the spin conductivity
($\textrm{Im}[\sigma_{xy}^{z}(\bm{q})]$) as a function of $T$. Different
curves represent distinct values of the energy gap ($E_{g}$). (c) Real part of
the spin conductivity ($\textrm{Re}[\sigma_{xy}^{z}(\bm{q})]$) as a function
of $E_{g}$. (d) Imaginary part of the spin conductivity
($\textrm{Im}[\sigma_{xy}^{z}(\bm{q})]$) as a function of $E_{g}$. Different
curves represent distinct values of $T$. In each plot, the dashed line
represents the theoretical curve for the spin conductivity based on the
analytical formula in Eq. (III). Additionally, the magnitude of the spin
conductivity is normalized by the overall factor
$\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}$. The parameter values used are as
follows: $E_{F}=3\textrm{eV}$, $v_{F}=10^{6}\textrm{m/s}$,
$\tau=10^{-14}\textrm{s}$, $c=10^{4}\textrm{m/s}$, and
$\bm{q}=(0.001k_{f},0,0)$.

The second of our main results is the scaling laws with respect to temperature
governing the phonon-mediated spin conductivity, as represented by:

$\displaystyle\sigma_{ij}^{\alpha}(\bm{q})\propto|\bm{q}|^{2}\bigg{(}\frac{k_{B}T\tau}{\hbar}\bigg{)}^{2}\exp\Bigg{[}-a\bigg{(}\frac{E_{g}}{k_{B}T}\bigg{)}^{b}\Bigg{]},$
(34)

where $a\approx 0.75$ and $b\approx 1.05$, as obtained in Eq. (21).
Importantly, $\sigma_{ij}^{\alpha}(\bm{q})$ becomes zero at absolute zero
temperature ($T=0$). Conversely, the spin conductivity undergoes enhancement
as the temperature increases. This behavior stems from the mediation of
thermally excited TO phonons. This thermal activation behavior stands in stark
contrast to the previous intrinsic spin Hall conductivity, which maintains a
nonzero value at the zero temperature [28]. It also represents a significant
departure from the impact of acoustic phonons, which lead to the degradation
of electrical conductivity [43]. It is worth noting that the quadratic
temperature dependence in the scaling formula originates from the phase volume
factor $k_{\textrm{cutoff}}^{3}=\big{(}\frac{k_{B}T}{\hbar c}\big{)}^{3}$
divided by an additional factor of $k_{B}T$ stemming from the phonon energy
$\omega_{\bm{k}}$. Roughly speaking, this factor accounts for the density of
thermally excited TO phonons that give rise to the generation of the PMSP.
Additionally, the presence of four electron propagators results in a cubic
dependence on the electron’s scattering time $\tau$, another crucial departure
from the usual linear dependence observed in conductivity.

When $E_{g}\ll k_{B}T$, as indicated by the red region in Fig. 1, TO phonons
exhibit gapless behavior, and thermal effects dominate. In this gapless-phonon
regime, the exponential factor can be disregarded, simplifying the scaling
formula in Eq. (34) as follows:

$\displaystyle\sigma_{ij}^{\alpha}(\bm{q})\propto|\bm{q}|^{2}\left(\frac{k_{B}T\tau}{\hbar}\right)^{2},\quad(E_{g}\ll
k_{B}T).$ (35)

Our numerical results confirm this power-law relationship, as depicted in Fig.
6(a,b). This quadratic temperature scaling is a distinctive feature of the
phonon-mediated spin conductivity, specifically in the gapless-phonon regime.

However, it needs to be recognized that as temperatures decrease, the system
eventually enters the quantum critical regime $(E_{g}\ll k_{B}T\lesssim
k_{B}T_{*})$, as indicated by the blue region in Fig. 1. In this regime, our
results may not be applicable due to the significant impact of self-energy
corrections, leading to the breakdown of the quasiparticle concept within the
Fermi liquid theory [27], which forms the basis of our calculations. Our
analysis of self-energy corrections provides an estimate of $T_{*}\sim
T_{0}\exp(-R)$, where
$R=16\Big{(}\frac{E_{F}^{2}}{g^{2}n}\Big{)}\Big{(}\frac{M_{ph}c^{2}}{E_{F}}\Big{)}$
and $T_{0}=T_{F}\frac{R}{6\pi^{2}}\Big{(}\frac{\Lambda}{k_{F}}\Big{)}^{3}$,
with $\Lambda$ representing the phonon cutoff wave vector (see the Appendix G
for the derivation of $T_{*}$). Consequently, the quadratic scaling law for
phonon-mediated spin conductivity remains valid within a limited temperature
range $T>T_{*}$.

Conversely, when $E_{g}\gg k_{B}T$, the exponential factor dominates, and the
power-law component can be neglected. In such cases, the scaling formula in
Eq. (34) can be expressed simply as:

$\displaystyle\sigma_{ij}^{\alpha}(\bm{q})\propto|\bm{q}|^{2}\exp\left[-a\left(\frac{E_{g}}{k_{B}T}\right)^{b}\right],\quad(E_{g}\gg
k_{B}T).$ (36)

Our numerical results confirm this exponential behavior in the gapped-phonon
regime, as shown in Fig. 6(c,d). This thermal activation behavior emphasizes
that the spin Hall effect is facilitated by thermally excited TO phonons, and
the presence of TO phonons becomes increasingly favorable as their energy gap
diminishes, enhancing the spin conductivity. We suggest that this scaling law
is applicable in the “classical” paraelectric phase near the quantum critical
point where $E_{g}\gg k_{B}T$, as depicted by the yellow region in Fig. 1.

In principle, our calculations can be extended to the quantum critical regime
by accounting for self-energy effects. In this scenario, the detailed
characteristics of the spin conductivity, such as the quadratic scaling law,
may undergo potential modifications. Nevertheless, the presence of the spin
conductivity remains robust since the mediation of the spin conductivity by TO
phonons, as described by the PMSP in Eq. (10), sorely relies on Rashba-type
spin-orbit coupling, which is a fundamental aspect of the system, irrespective
of the specific phonon or electron propagator. Consequently, we argue that the
spin conductivity persists even in the quantum critical regime, though it may
not necessarily adhere to the specific form, as presented in Eq. (34), and
could manifest with a different temperature dependence characterized by a
different exponent and numerical coefficients.

##  VI Anisotropy effects

Many metals near the ferroelectric phase transition exhibit quasi-two-
dimensional behavior or belong to the category of two-dimensional (2D)
materials [2, 3, 4]. The crystalline anisotropy in such systems significantly
impacts the energy dispersion of TO phonons, leading to either an easy-plane
or easy-axis behavior. Additionally, the anisotropic character can be
experimentally manipulated using external stimuli, such as strain [44].
Furthermore, the electron-phonon coupling effect on the electrons’ dynamics is
quite different in 2D anisotropic systems compared to 3D isotropic systems
[27]. Hence, it’s crucial to consider anisotropy when applying our theoretical
model to the candidate materials for the quantum paraelectric metal.

Under these considerations, we extend our previous three-dimensional (3D)
isotropic results, as presented in Eqs. (III) and (17), to encompass
anisotropic 2D scenarios, including easy-plane and easy-axis cases. The
primary modification arises in the factor
$\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$ in Eq. (10), which is adjusted
as follows:

$\displaystyle\tilde{\gamma}^{\alpha,\textrm{EP}}(\bm{p},\bm{k};\bm{q})$
$\displaystyle=(1-\delta_{\alpha
z})\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q}),$ (37)
$\displaystyle\tilde{\gamma}^{\alpha,\textrm{EA}}(\bm{p},\bm{,}k;\bm{q})$
$\displaystyle=\delta_{\alpha
z}\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q}).$ (38)

Here, the superscripts “EP” and “EA” denote “easy-plane” and “easy-axis,”
respectively. In comparison to the isotropic case,
$\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$ now exhibits directional
preferences denoted by projection factors such as $(1-\delta_{\alpha z})$ and
$\delta_{\alpha z}$. These projection factors constrain the direction of the
spin polarization to the easy-plane $\alpha=\\{x,y\\}$ or to the easy-axis
$\alpha=z$ in each case. This spin polarization constraint clearly
differentiates the spin conductivity in anisotropic systems from that of
isotropic systems.

To precisely quantify the deviations from the isotropic scenario, we calculate
the phonon-mediated for both easy-plane and easy-axis cases, which are denoted
as $\sigma^{\alpha,\textrm{EP}}_{ij}(\bm{q})$ and
$\sigma^{\alpha,\textrm{EA}}_{ij}(\bm{q})$, respectively, and find the same
quadrupolar symmetry with the different response for the transverse and the
longitudinal field, albeit with the spin polarization constraint of Eqs. (37)
and (38). The calculations closely mirror those of the isotropic case; for
further details, refer to the Appendix F. As a result, we obtain
$\sigma^{\alpha,\textrm{EP}}_{ij}(\bm{q})$ and
$\sigma^{\alpha,\textrm{EA}}_{ij}(\bm{q})$ up to quadratic order in $\bm{q}$
as:

$\displaystyle\sigma^{\alpha,\textrm{EP}}_{ij}(\bm{q})=$
$\displaystyle\;(1-\delta_{\alpha
z})\Big{[}\chi_{0}^{\textrm{2D}}\big{(}\nu_{ij}^{\alpha}(\bm{q})+\nu_{ji}^{\alpha}(\bm{q})\big{)}+i\kappa_{0}^{\textrm{2D}}\big{(}-\nu_{ij}^{\alpha}(\bm{q})+\nu_{ji}^{\alpha}(\bm{q})\big{)}\Big{]},$
(39) $\displaystyle\sigma^{\alpha,\textrm{EA}}_{ij}(\bm{q})=$
$\displaystyle\;\delta_{\alpha
z}\Big{[}\chi_{0}^{\textrm{2D}}\big{(}\nu_{ij}^{\alpha}(\bm{q})+\nu_{ji}^{\alpha}(\bm{q})\big{)}+i\kappa_{0}^{\textrm{2D}}\big{(}-\nu_{ij}^{\alpha}(\bm{q})+\nu_{ji}^{\alpha}(\bm{q})\big{)}\Big{]}.$
(40)

In contrast to the isotropic case, the spin quantization direction $\alpha$
cannot be freely adjusted any longer; for the easy-plane case, $\alpha$ is
constrained to the easy-plane directions, while for the easy-axis case, it
aligns with the easy-axis direction as indicated by the factors
$1-\delta_{\alpha z}$ and $\delta_{\alpha z}$, respectively. For the easy axis
case, the spin current response to the longitudinal electric field is now
exactly the same form as that of the Rashba metal [37].

The qualitative change due to the 2D nature does occur in the temperature
dependence, as can be seen from the coefficients of the spin conductivity,

$\displaystyle\chi_{0}^{\textrm{2D}}$
$\displaystyle=\bigg{(}\frac{n_{\textrm{2D}}e\hbar\tau}{m}\bigg{)}\bigg{(}\frac{g^{2}\tau}{\hbar
M_{ph}c^{2}}\bigg{)}\bigg{(}\frac{k_{B}T\tau}{\hbar}\bigg{)}\bigg{(}\frac{E_{F}\tau}{\hbar}\bigg{)}\frac{3}{2\pi}I_{\textrm{2D}}\bigg{(}\frac{E_{g}}{k_{B}T}\bigg{)},$
(41) $\displaystyle\kappa_{0}^{\textrm{2D}}$
$\displaystyle=\bigg{(}\frac{n_{\textrm{2D}}e\hbar\tau}{m}\bigg{)}\bigg{(}\frac{g^{2}\tau}{\hbar
M_{ph}c^{2}}\bigg{)}\bigg{(}\frac{k_{B}T\tau}{\hbar}\bigg{)}\frac{1}{2\pi}I_{\textrm{2D}}\bigg{(}\frac{E_{g}}{k_{B}T}\bigg{)},$
(42)

where $n_{\textrm{2D}}=\frac{k_{F}^{2}}{2\pi}$ represents the electron density
in 2D. It is $I_{\textrm{2D}}\big{(}\frac{E_{g}}{k_{B}T}\big{)}$, an integral
function that encapsulates the integration over $|\bm{k}|$,

$\displaystyle
I_{\textrm{2D}}(y)=\int^{\infty}_{0}dx\frac{x}{2\sqrt{x^{2}+y^{2}}\sinh{\big{[}\sqrt{x^{2}+y^{2}}\;\big{]}}}=\frac{1}{2}\ln\bigg{(}\coth\frac{y}{2}\bigg{)},$
(43)

that gives rise to the most important qualitative change in the temperature
and the phonon energy gap dependence. In the gapless-phonon regime, $E_{g}\ll
k_{B}T$ or $y\ll 1$, (the red region in Fig. 1), $I_{\textrm{2D}}(y)$ can be
well-approximated as
$I_{\textrm{2D}}(y)\approx\frac{1}{2}\ln\big{(}\frac{2}{y}\big{)}$,
simplifying the scaling formula in Eq. (39) as follows:

$\displaystyle\sigma_{ij}^{\alpha,\textrm{EP(EA)}}(\bm{q})\propto|\bm{q}|^{2}\left(\frac{k_{B}T\tau}{\hbar}\right)\ln\bigg{(}\frac{2k_{B}T}{E_{g}}\bigg{)},\quad(E_{g}\ll
k_{B}T).$ (44)

Conversely, in the gapped-phonon regime, $E_{g}\gg k_{B}T$ or $y\gg 1$, (the
yellow region in Fig. 1), $I_{\textrm{2D}}(y)$ can be approximated as
$I_{\textrm{2D}}(y)\approx\exp(-y)$, simplifying the scaling formula in Eq.
(39) as follows:

$\displaystyle\sigma_{ij}^{\alpha,\textrm{EP(EA)}}(\bm{q})\propto|\bm{q}|^{2}\exp\left[-\left(\frac{E_{g}}{k_{B}T}\right)\right],\quad(E_{g}\gg
k_{B}T).$ (45)

It is possible that in the easy-axis case, the temperature dependence in Eqs.
(44) and (45) may undergo substantial modifications at the QCP where the
effect of self-corrections becomes significant. However, in the easy-plane
case, such effect is weaker due to the phase volume restriction allowed for
electron scattering, and the dynamics of electrons maintains a Fermi liquid
character even in the quantum critical regime [27]. Consequently, the above
scaling laws may remain valid even in the quantum critical regime within our
leading order-perturbation theory in $g$.

##  VII Conclusion

In this study, we have demonstrated that quantum paraelectric metals near
ferroelectric quantum criticality can exhibit an unconventional phonon-
mediated spin transport in response to an inhomogeneous electric field, an
example of spin transport arising from interaction rather than band structure.
Our rigorous calculations, employing the Kubo formula and a perturbative
expansion in electron-phonon interaction, have unveiled that soft transverse
optical phonons, with their intrinsic Rashba-type spin-orbit coupling to
electrons, can serve as unconventional contributors to the spin conductivity
in response to the inhomogeneous electric field. Furthermore, the resulting
spin conductivity displays a couple of intriguing and unique characteristics.
One is that it exhibits unconventional quadrupolar symmetry associated with
the $\bm{q}$-vector, leading to a possible nonzero response, in contrast to
the theoretical prediction for quantum Hall states and Rashba metals, even
when $\bm{\nabla}\cdot\bm{E}=0$. The resulting spin current is, therefore,
non-uniform, and its observation may require a local spin probe such as the
X-ray magnetic circular dichroism [45] or the spin-torque transfer
ferromagnetic resonance [46]. The other is that it follows distinctive scaling
laws in temperature and phonon energy gap; the conductivity increases as
temperature rises or the energy gap diminishes, as it is mediated by thermally
excited phonons. Consequently, the proposed spin transport may be best
observed in the vicinity of a ferroelectric quantum criticality, where the
phonon energy gap diminishes while thermal effects amplify, even if how much
of Fermi liquid theory-based temperature dependence would hold in the quantum
critical regime remains to be examined in future research. Therefore, the
recent report of the ferroelectric quantum critical point in the n-type SrTiO3
may provide one example [47, 48]; additionally, it would also be interesting
to look further and search for other and different instances of interaction-
induced spin transport. Lastly, we note that the impurity effect on the spin
conductivity is an interesting issue to be addressed, e.g. whether it would be
analogous to that of the Rashba metal spin Hall conductivity [49, 50].

##  VIII Methods

###  VIII.1 Expansion of spin conductivity to the wave vector of electric
fields

To obtain the quadratic term of $\sigma^{\alpha}_{ij}(\bm{q})$ presented in
Eq. (III), we expanded the full expression of Eq. (III) to the wave vector
$\bm{q}$ of electric fields. Utilizing the expansion of
$G_{\textrm{ret}}(0,\bm{k}+\bm{q})$ to $\bm{q}$:

$\displaystyle G_{\textrm{ret}}(0,\bm{k}+\bm{q})$
$\displaystyle=\;G_{\textrm{ret}}(0,\bm{k})\bigg{(}1-G_{\textrm{ret}}(0,\bm{k})\frac{\hbar^{2}\bm{k}\cdot\bm{q}}{2m}-G_{\textrm{ret}}(0,\bm{k})\frac{\hbar^{2}\bm{q}^{2}}{2m}\bigg{)}^{-1}$
$\displaystyle\approx
G_{\textrm{ret}}(0,\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{k}\cdot\bm{q}}{2m}+\frac{\hbar^{2}\bm{q}^{2}}{2m}\bigg{)}G_{\textrm{ret}}(0,\bm{k})^{2}+\bigg{(}\frac{\hbar^{2}\bm{k}\cdot\bm{q}}{2m}\bigg{)}^{2}G_{\textrm{ret}}(0,\bm{k})^{3}+\mathcal{O}(q^{3}),$
(46)

and a similar expression for
$G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})$, we expanded the
product of the electron propagators in Eq. (III) as:

$\displaystyle
G_{\textrm{ret}}(0,\bm{p}+\bm{q})G_{\textrm{adv}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})G_{\textrm{adv}}(\omega_{\bm{k}},\bm{p}+\bm{k})+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}})$
$\displaystyle=$
$\displaystyle\;|G_{\textrm{ret}}(0,\bm{p})|^{2}|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})|^{2}\bigg{\\{}1+\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})$
$\displaystyle+\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{2}[G_{\textrm{ret}}(0,\bm{p})]^{2}+\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{2}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}$
$\displaystyle+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\mathcal{O}(q^{3})\bigg{\\}}+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}}).$
(47)

Substituting Eq. (47) into Eq. (III) and gathering terms of quadratic order in
$\bm{q}$, we obtain

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}$
$\displaystyle\times\Bigg{[}\Bigg{\\{}|G_{\textrm{ret}}(0,\bm{p})|^{2}|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})|^{2}\Bigg{(}\frac{1}{2}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{p})-\frac{1}{2}\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{q})$
$\displaystyle+\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{p})\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})\bigg{)}\Bigg{)}\Bigg{\\}}+\Bigg{\\{}\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}}\Bigg{\\}}\Bigg{]}.$
(48)

This expression is further simplified by approximating
$(\bm{p}+\bm{k}/2)\approx\bm{p}$, $(\bm{p}+\bm{k})\approx\bm{p}$, and
$G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})\approx
G_{\textrm{ret}}(0,\bm{p})$ [38]. This approximation is ensured by the
distribution factor
$\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{\omega_{\bm{k}}}$ that
exponentially decays for a large $|\bm{k}|$, effectively constraining
$|\bm{k}|\ll\frac{Ck_{B}T}{\hbar c}$, where $C$ is an arbitrary constant. In a
low-temperature regime $T\ll\frac{\hbar ck_{F}}{Ck_{B}}$ of our interest, this
constraint indicates $|\bm{k}|\ll k_{F}\sim|\bm{p}|$, thus validating the
stated approximation. Within this framework, we approximate Eq. (VIII.1) as:

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}|G_{\textrm{ret}}(0,\bm{p})|^{4}$
$\displaystyle\times\Big{\\{}\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot(\bm{p}\times\hat{\bm{k}})\Big{\\}}(q_{i}p_{j}-p_{i}q_{j})$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}|G_{\textrm{ret}}(0,\bm{p})|^{4}G_{\textrm{ret}}(0,\bm{p})$
$\displaystyle\times\Big{\\{}\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot(\bm{p}\times\hat{\bm{k}})\Big{\\}}p_{i}p_{j}\frac{4\hbar^{2}\bm{p}\cdot\bm{q}}{m}.$
(49)

The requirement of the integrand being even in $\bm{k}$ and $\bm{p}$ leads to:

$\displaystyle\Big{\\{}\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot(\bm{p}\times\hat{\bm{k}})\Big{\\}}(q_{i}p_{j}-p_{i}q_{j})\rightarrow\Big{\\{}\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})q_{i}p_{j}^{2}-\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})q_{j}p_{i}^{2}\Big{\\}}(1-\hat{k}_{\alpha}^{2}),$
$\displaystyle\Big{\\{}\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot(\bm{p}\times\hat{\bm{k}})\Big{\\}}p_{i}p_{j}\rightarrow\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})(1-\hat{k}_{\alpha}^{2})p_{i}p_{j}.$
(50)

Substituting Eq. (VIII.1) into Eq. (VIII.1), we obtain Eq. (III).

###  VIII.2 Evaluation of the spin conductivity

To obtain the explicit expression for $\sigma^{\alpha}_{ij}(\bm{q})$ presented
in Eq. (17), we performed an analytic integration of Eq. (III) by
approximating the integrated as follows: (i) the electron propagator is
represented as $G_{\textrm{ret}}(0,\bm{p})\approx(-\hbar
v_{F}p+i\hbar/2\tau)^{-1}\equiv G_{\textrm{ret}}(0,p)$, where
$p\equiv|\bm{p}|-k_{F}$, (ii) the factors of $\bm{p}$ are replaced by
$k_{F}\hat{\bm{p}}$, and (iii) the integral measure is represented as $\int
d^{3}p\approx\int^{\infty}_{-\infty}dpk_{F}^{2}\int d\Omega_{p}$, where $\int
d\Omega_{p}\equiv\int^{1}_{-1}d\cos{\theta_{p}}\int^{2\pi}_{0}d\phi_{p}$.
These approximations are justified under the condition the condition
$E_{F}\tau\gg\hbar$, ensuring that the integral over $\bm{p}$ sharply peeks
around $|\bm{p}|\approx k_{F}$. Additionally, we represent the integral over
$\bm{k}$ as $\int d^{3}k=\int^{\infty}_{0}dkk^{2}\int d\Omega_{k}$. Within
this framework, we approximate $\sigma^{\alpha}_{ij}(\bm{q})$ in Eq. (III) as:

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\frac{k_{F}^{4}q^{2}}{(2\pi)^{6}}\int^{\infty}_{0}dkk^{2}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\int
d\Omega_{k}(1-\hat{k}_{\alpha}^{2})$
$\displaystyle\times\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}\int
d\Omega_{p}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{\\}}-\big{\\{}i\leftrightarrow
j\big{\\}}\Big{]},$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\frac{4\hbar^{2}k_{F}^{6}q^{2}}{(2\pi)^{6}m}\int^{\infty}_{0}dkk^{2}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\int
d\Omega_{k}(1-\hat{k}_{\alpha}^{2})$
$\displaystyle\times\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)\int
d\Omega_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{p}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}.$
(51)

The integrals for $\Omega_{k}$, $p$, and $\Omega_{p}$ can be conducted
analytically, yielding

$\displaystyle\int d\Omega_{k}(1-\hat{k}_{\alpha}^{2})$
$\displaystyle=\frac{8\pi}{3},$
$\displaystyle\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}$
$\displaystyle=\frac{4\pi}{\hbar v_{F}(\hbar/\tau)^{3}},$
$\displaystyle\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)$
$\displaystyle=\frac{-6\pi i}{\hbar v_{F}(\hbar/\tau)^{4}},$ (52)

and

$\displaystyle\int
d\Omega_{p}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{\\}}-\big{\\{}i\leftrightarrow
j\big{\\}}\Big{]}$
$\displaystyle=\frac{4\pi}{3}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{\\}}-\big{\\{}i\leftrightarrow
j\big{\\}}\Big{]},$ $\displaystyle\int
d\Omega_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{p}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}(\hat{\bm{p}}\cdot\hat{\bm{q}})$
$\displaystyle=\frac{4\pi}{15}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{]}\big{\\}}+\\{i\leftrightarrow
j\\}\Big{]}.$ (53)

The remaining integral for $k$ is represented as:

$\displaystyle\int^{\infty}_{0}dkk^{2}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}=\frac{1}{k_{B}T}\bigg{(}\frac{k_{B}T}{\hbar
c}\bigg{)}^{3}I\bigg{(}\frac{E_{g}}{k_{B}T}\bigg{)},$ (54)

where $I(y)$ is defined in Eq. (21). Substituting Eqs. (VIII.2) to (54) into
Eq. (VIII.2), we obtain Eq. (17).

###  VIII.3 Numerical integration method

To produce the results depicted in Figs. 4 and 6, we performed numerical
integration of the $\sigma_{ij}^{\alpha}(\bm{q})$ formula given in Eq. (III)
using a Riemann sum approach. The discretization of the electron wave vector
$\bm{p}=(p_{x},p_{y},p_{z})$ spanned the Fermi surface region, defined as
$\sqrt{k_{F}^{2}-9.95\cdot(\hbar/2\tau)}<|\bm{p}|<\sqrt{k_{F}^{2}+9.95\cdot(\hbar/2\tau)}$.
At the boundaries of this region, we ensured that the electron spectral
function, $-\frac{1}{\pi}\textrm{Im}[G_{\textrm{ret}}]$, descended below 1% of
its maximum value occurring at $|\bm{p}|=k_{F}$. Consequently, contributions
from outside this region were deemed negligible. Additionally, we established
a mesh for the phonon wave vector $\bm{k}=(k_{x},k_{y},k_{z})$ within the
range $|\bm{k}|\leq k_{\textrm{cutoff}}=10\frac{k_{B}T}{\hbar c}$. At the
boundaries of this region, we ensured that the distribution factor
${n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}$ fell below a few percent of
its maximum value occurring at $|\bm{k}|=0$ (the ratio was maintained at 1%
for $E_{g}=0$ and less than 10% for $E_{g}\approx 10k_{B}T$). Consequently,
contributions from outside this region were also considered negligible. For
the computation, we used a grid size of $(51\times 51\times 51)$ for both
$(p_{x},p_{y},p_{z})$ and $(k_{x},k_{y},k_{z})$ within the specified regions
for numerical integration. Our analyses confirmed that the numerical
integration results were well-converged and reliable for our intended
purposes.

## Acknowledgements

We would like to thank Jung Hoon Han, Beom Hyun Kim, Changyoung Kim, Ki-Seok
Kim, Se Kwon Kim, Hyun-Woo Lee, and S. Raghu for sharing their insights.
K.-M.K. was supported by the Institute for Basic Science in the Republic of
Korea through the project IBS-R024-D1. S.B.C. was supported by the National
Research Foundation of Korea (NRF) grants funded by the Korea government
(MSIT) (NRF-2023R1A2C1006144, NRF-2020R1A2C1007554, and
NRF-2018R1A6A1A06024977).

## References

  * Anderson and Blount [1965] P. W. Anderson and E. I. Blount, Symmetry considerations on martensitic transformations: “ferroelectric” metals?, Phys. Rev. Lett. 14, 217 (1965).
  * Zhou and Ariando [2020] W. X. Zhou and A. Ariando, Review on ferroelectric/polar metals, Jpn. J. Appl. Phys. 59, SI0802 (2020).
  * Hickox-Young  _et al._ [2023] D. Hickox-Young, D. Puggioni, and J. M. Rondinelli, Polar metals taxonomy for materials classification and discovery, Phys. Rev. Mater. 7, 010301 (2023).
  * Bhowal and Spaldin [2023] S. Bhowal and N. A. Spaldin, Polar metals: Principles and prospects, Annu. Rev. Mater. Res. 53, 53 (2023).
  * Shi  _et al._ [2013] Y. Shi, Y. Guo, X. Wang, A. J. Princep, D. Khalyavin, P. Manuel, Y. Michiue, A. Sato, K. Tsuda, S. Yu, M. Arai, Y. Shirako, M. Akaogi, N. Wang, K. Yamaura, and A. T. Boothroyd, A ferroelectric-like structural transition in a metal, Nat. Mater. 12, 1024 (2013).
  * Rischau  _et al._ [2017] C. W. Rischau, X. Lin, C. P. Grams, D. Finck, S. Harms, J. Engelmayer, T. Lorenz, Y. Gallais, B. Fauqué, J. Hemberger, and K. Behnia, A ferroelectric quantum phase transition inside the superconducting dome of Sr1-xCaxTiO3-δ, Nat. Phys. 13, 643 (2017).
  * Salmani-Rezaie  _et al._ [2020a] S. Salmani-Rezaie, K. Ahadi, W. M. Strickland, and S. Stemmer, Order-disorder ferroelectric transition of strained ${\mathrm{srtio}}_{3}$, Phys. Rev. Lett. 125, 087601 (2020a).
  * Salmani-Rezaie  _et al._ [2020b] S. Salmani-Rezaie, K. Ahadi, and S. Stemmer, Polar nanodomains in a ferroelectric superconductor, Nano Lett. 20, 6542 (2020b).
  * Yu  _et al._ [2018] H. Yu, D. Gao, X. Wang, X. Du, X. Lin, W. Guo, R. Zou, C. Jin, K. Li, and Y. Chen, Unraveling a novel ferroelectric gese phase and its transformation into a topological crystalline insulator under high pressure, NPG Asia Mater. 10, 882 (2018).
  * Bilz  _et al._ [1983] H. Bilz, A. Bussmann-Holder, W. Jantsch, and P. Vogl,  _Dynamical Properties of IV-VI Compounds_ (Springer Berlin Heidelberg, 1983).
  * Fei  _et al._ [2018] Z. Fei, W. Zhao, T. A. Palomaki, B. Sun, M. K. Miller, Z. Zhao, J. Yan, X. Xu, and D. H. Cobden, Ferroelectric switching of a two-dimensional metal, Nature 560, 336 (2018).
  * Sharma  _et al._ [2019] P. Sharma, F.-X. Xiang, D.-F. Shao, D. Zhang, E. Y. Tsymbal, A. R. Hamilton, and J. Seidel, A room-temperature ferroelectric semimetal, Sci. Adv. 5, eaax5080 (2019).
  * Sakai  _et al._ [2016] H. Sakai, K. Ikeura, M. S. Bahramy, N. Ogawa, D. Hashizume, J. Fujioka, Y. Tokura, and S. Ishiwata, Critical enhancement of thermopower in a chemically tuned polar semimetal MoTe2, Sci. Adv. 2, e1601378 (2016).
  * Rowley  _et al._ [2014] S. E. Rowley, L. J. Spalek, R. P. Smith, M. P. M. Dean, M. Itoh, J. F. Scott, G. G. Lonzarich, and S. S. Saxena, Ferroelectric quantum criticality, Nat. Phys. 10, 367 (2014).
  * Edge  _et al._ [2015] J. M. Edge, Y. Kedem, U. Aschauer, N. A. Spaldin, and A. V. Balatsky, Quantum critical origin of the superconducting dome in SrTiO3, Phys. Rev. Lett. 115, 247002 (2015).
  * Kozii  _et al._ [2019] V. Kozii, Z. Bi, and J. Ruhman, Superconductivity near a ferroelectric quantum critical point in ultralow-density dirac materials, Phys. Rev. X 9, 031046 (2019).
  * Volkov  _et al._ [2022] P. A. Volkov, P. Chandra, and P. Coleman, Superconductivity from energy fluctuations in dilute quantum critical polar metals, Nat. Commun. 13, 4599 (2022).
  * Gastiasoro  _et al._ [2020a] M. N. Gastiasoro, J. Ruhman, and R. M. Fernandes, Superconductivity in dilute SrTiO3: A review, Ann. Phys. 417, 168107 (2020a), eliashberg theory at 60: Strong-coupling superconductivity and beyond.
  * Chandra  _et al._ [2017] P. Chandra, G. G. Lonzarich, S. E. Rowley, and J. F. Scott, Prospects and applications near ferroelectric quantum phase transitions: a key issues review, Rep. Prog. Phys. 80, 112502 (2017).
  * Yu  _et al._ [2022] Y. Yu, H. Y. Hwang, S. Raghu, and S. B. Chung, Theory of superconductivity in doped quantum paraelectrics, npj Quantum Mater. 7, 63 (2022).
  * Cohen [1992] R. E. Cohen, Origin of ferroelectricity in perovskite oxides, Nature 358, 136 (1992).
  * Wang  _et al._ [2012] Y. Wang, X. Liu, J. D. Burton, S. S. Jaswal, and E. Y. Tsymbal, Ferroelectric instability under screened coulomb interactions, Phys. Rev. Lett. 109, 247601 (2012).
  * Lines and Glass [2001] M. E. Lines and A. M. Glass,  _Principles and Applications of Ferroelectrics and Related Materials_ (Oxford University Press, 2001).
  * Gastiasoro  _et al._ [2019] M. N. Gastiasoro, A. V. Chubukov, and R. M. Fernandes, Phonon-mediated superconductivity in low carrier-density systems, Phys. Rev. B 99, 094524 (2019).
  * Lee  _et al._ [2020] M. Lee, H.-J. Lee, J. H. Lee, and S. B. Chung, Topological superconductivity from transverse optical phonons in oxide heterostructures, Phys. Rev. Mater. 4, 034202 (2020).
  * Gastiasoro  _et al._ [2020b] M. N. Gastiasoro, T. V. Trevisan, and R. M. Fernandes, Anisotropic superconductivity mediated by ferroelectric fluctuations in cubic systems with spin-orbit coupling, Phys. Rev. B 101, 174501 (2020b).
  * Klein  _et al._ [2023] A. Klein, V. Kozii, J. Ruhman, and R. M. Fernandes, Theory of criticality for quantum ferroelectric metals, Phys. Rev. B 107, 165110 (2023).
  * Sinova  _et al._ [2004] J. Sinova, D. Culcer, Q. Niu, N. A. Sinitsyn, T. Jungwirth, and A. H. MacDonald, Universal intrinsic spin hall effect, Phys. Rev. Lett. 92, 126603 (2004).
  * Shen [2004] S.-Q. Shen, Spin hall effect and berry phase in two-dimensional electron gas, Phys. Rev. B 70, 081311 (2004).
  * Tanaka  _et al._ [2008] T. Tanaka, H. Kontani, M. Naito, T. Naito, D. S. Hirashima, K. Yamada, and J. Inoue, Intrinsic spin hall effect and orbital hall effect in $4d$ and $5d$ transition metals, Phys. Rev. B 77, 165117 (2008).
  * Kontani  _et al._ [2009] H. Kontani, T. Tanaka, D. S. Hirashima, K. Yamada, and J. Inoue, Giant orbital hall effect in transition metals: Origin of large spin and anomalous hall effects, Phys. Rev. Lett. 102, 016601 (2009).
  * Park  _et al._ [2011] S. R. Park, C. H. Kim, J. Yu, J. H. Han, and C. Kim, Orbital-angular-momentum based origin of rashba-type surface band splitting, Phys. Rev. Lett. 107, 156803 (2011).
  * Go  _et al._ [2018] D. Go, D. Jo, C. Kim, and H.-W. Lee, Intrinsic spin and orbital hall effects from orbital texture, Phys. Rev. Lett. 121, 086602 (2018).
  * Schlipf and Giustino [2021] M. Schlipf and F. Giustino, Dynamic rashba-dresselhaus effect, Phys. Rev. Lett. 127, 237601 (2021).
  * Haldane [2009] F. D. M. Haldane, “hall viscosity” and intrinsic metric of incompressible fractional hall fluids (2009), arXiv:0906.1854 [cond-mat.str-el] .
  * Hoyos and Son [2012] C. Hoyos and D. T. Son, Hall viscosity and electromagnetic response, Phys. Rev. Lett. 108, 066805 (2012).
  * Zhang and Rhim [2022] A. Zhang and J.-W. Rhim, Geometric origin of intrinsic spin hall effect in an inhomogeneous electric field, Commun. Phys. 5, 195 (2022).
  * Mahan [2000] G. D. Mahan,  _Many Particle Physics, Third Edition_ (Plenum, New York, 2000).
  * Avron  _et al._ [1995] J. E. Avron, R. Seiler, and P. G. Zograf, Viscosity of quantum hall fluids, Phys. Rev. Lett. 75, 697 (1995).
  * Read [2009] N. Read, Non-abelian adiabatic statistics and hall viscosity in quantum hall states and ${p}_{x}+i{p}_{y}$ paired superfluids, Phys. Rev. B 79, 045308 (2009).
  * Kimura [2021] T. Kimura, Hall and spin hall viscosity in 2d topological systems, J. Phys. Soc. Jpn. 90, 064705 (2021).
  * Note [1] The coupling of electric field to the TO phonon, omitted in our analysis, can be adequately treated with the static dielectric constant for frequency far below $E_{g}/\hbar$ .
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Holt-Saunders, 1976).
  * Haeni  _et al._ [2004] J. H. Haeni, P. Irvin, W. Chang, R. Uecker, P. Reiche, Y. L. Li, S. Choudhury, W. Tian, M. E. Hawley, B. Craigo, A. K. Tagantsev, X. Q. Pan, S. K. Streiffer, L. Q. Chen, S. W. Kirchoefer, J. Levy, and D. G. Schlom, Room-temperature ferroelectricity in strained srtio3, Nature 430, 758 (2004).
  * Stöhr [1999] J. Stöhr, Exploring the microscopic origin of magnetic anisotropies with x-ray magnetic circular dichroism (xmcd) spectroscopy, J. Magn. Magn. Mater. 200, 470 (1999).
  * Wang  _et al._ [2018] Y. Wang, R. Ramaswamy, and H. Yang, Fmr-related phenomena in spintronic devices, J. Phys. D: Appl. Phys. 51, 273002 (2018).
  * Rischau  _et al._ [2022] C. W. Rischau, D. Pulmannová, G. W. Scheerer, A. Stucky, E. Giannini, and D. van der Marel, Isotope tuning of the superconducting dome of strontium titanate, Phys. Rev. Res. 4, 013019 (2022).
  * Tomioka  _et al._ [2022] Y. Tomioka, N. Shirakawa, and I. H. Inoue, Superconductivity enhancement in polar metal regions of sr 0.95 ba 0.05 tio 3 and sr 0.985 ca 0.015 tio 3 revealed by systematic nb doping, npj Quantum Mater. 7, 111 (2022).
  * Mishchenko  _et al._ [2004] E. G. Mishchenko, A. V. Shytov, and B. I. Halperin, Spin current and polarization in impure two-dimensional electron systems with spin-orbit coupling, Phys. Rev. Lett. 93, 226602 (2004).
  * Nomura  _et al._ [2005] K. Nomura, J. Sinova, N. A. Sinitsyn, and A. H. MacDonald, Dependence of the intrinsic spin-hall effect on spin-orbit interaction character, Phys. Rev. B 72, 165316 (2005).

##  Appendix A Formal expression of the conductivity

We derive the formula of the spin conductivity presented in Eq. (II) of the
main text. Using the diagrammatic expansion (Fig. 2), we obtain the following
expression for the correlation function:

$\displaystyle\pi^{\alpha}_{ij}(i\omega_{n},\bm{q})=\frac{e\hbar^{3}}{m^{2}}\int\frac{d^{3}p}{(2\pi)^{3}}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}S_{i}^{\alpha}(i\omega_{n};\bm{p}+\bm{q},\bm{p}),$
(55)

where $S_{i}^{\alpha}(i\omega_{n};\bm{p}+\bm{q},\bm{p})$ is given by

$\displaystyle S^{\alpha}_{i}(i\omega_{n};\bm{p}+\bm{q},\bm{p})=$
$\displaystyle\;\frac{1}{\beta}\sum_{ip_{m}}\mathcal{G}(ip_{m}+i\omega_{n},\bm{p}+\bm{q})\mathcal{G}(ip_{m},\bm{p})\gamma_{i}^{\alpha}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{q},\bm{p}).$
(56)

We can further calculate $S_{i}^{\alpha}(i\omega_{n};\bm{p}+\bm{q},\bm{p})$,
without assuming a specific form of
$\gamma_{i}^{\alpha}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{q},\bm{p})$. For
brevity, we omit the notations for the variables $\bm{p}+\bm{q}$ and $\bm{p}$
in $S^{\alpha}_{i}(i\omega_{n})$, $\mathcal{G}(ip_{m}+i\omega_{n})$,
$\mathcal{G}(ip_{m})$, and $\gamma_{i}^{\alpha}(ip_{m}+i\omega_{n},ip_{m})$,
and rewrite $S^{\alpha}_{i}(i\omega_{n})$ as

$\displaystyle
S^{\alpha}_{i}(i\omega_{n})=\frac{1}{\beta}\sum_{ip_{m}}\mathcal{G}(ip_{m}+i\omega_{n})\mathcal{G}(ip_{m})\gamma_{i}^{\alpha}(ip_{m}+i\omega_{n},ip_{m}).$
(57)

We transform the Matsubara frequency summation over $ip_{n}$ into a contour
integral as

$\displaystyle S^{\alpha}_{i}(i\omega_{n})=$
$\displaystyle\;-\oint_{\gamma}\frac{dz}{2\pi
i}n_{F}(z)\mathcal{G}(z+i\omega_{n})\mathcal{G}(z)\gamma_{i}^{\alpha}(z+i\omega_{n},z),$
(58)

where the contour $\gamma$ expands to infinity in the complex plane of $z$
while avoiding two branch cuts along $z=\epsilon$ and $z=\epsilon-i\omega_{n}$
with $\epsilon\in(-\infty,\infty)$. The integral is reduced along the branch
cuts as

$\displaystyle S^{\alpha}_{i}(i\omega_{n})=$
$\displaystyle\;-\int^{\infty}_{-\infty}\frac{d\epsilon}{2\pi
i}n_{F}(\epsilon)\Big{\\{}\mathcal{G}(\epsilon+i\omega_{n})G_{\textrm{ret}}(\epsilon)\gamma_{i}^{\alpha}(\epsilon+i\omega_{n},\epsilon+i\delta)-\mathcal{G}(\epsilon+i\omega_{n})G_{\textrm{adv}}(\epsilon)\gamma_{i}^{\alpha}(\epsilon+i\omega_{n},\epsilon-i\delta)$
$\displaystyle+G_{\textrm{ret}}(\epsilon)\mathcal{G}(\epsilon-i\omega_{n})\gamma_{i}^{\alpha}(\epsilon+i\delta,\epsilon-i\omega_{n})-G_{\textrm{adv}}(\epsilon)\mathcal{G}(\epsilon-i\omega_{n})\gamma_{i}^{\alpha}(\epsilon-i\delta,\epsilon-i\omega_{n})\Big{\\}},$
(59)

where $\delta$ is an infinitesimal. We perform an analytic continuation
$i\omega_{n}\rightarrow\omega+i\delta$ and a variable change
$\epsilon\rightarrow\epsilon+\omega$ in the last two terms to obtain

$\displaystyle S^{\alpha}_{i}(\omega+i\delta)=$
$\displaystyle\;\int^{\infty}_{-\infty}\frac{d\epsilon}{2\pi
i}\Big{\\{}\big{(}n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\big{)}G_{\textrm{ret}}(\epsilon+\omega)G_{\textrm{adv}}(\epsilon)\gamma_{i}^{\alpha}(\epsilon+\omega+i\delta,\epsilon-i\delta)$
$\displaystyle-
n_{F}(\epsilon)G_{\textrm{ret}}(\epsilon+\omega)G_{\textrm{ret}}(\epsilon)\gamma_{i}^{\alpha}(\epsilon+\omega+i\delta,\epsilon+i\delta)+n_{F}(\epsilon+\omega)G_{\textrm{adv}}(\epsilon+\omega)G_{\textrm{adv}}(\epsilon)\gamma_{i}^{\alpha}(\epsilon+\omega-i\delta,\epsilon-i\delta)\Big{\\}}.$
(60)

We drop the terms in the second line, which are given as the products of two
retarded propagators or two advanced propagators, as their contribution is
negligible when $\omega\ll E_{F}$ [38]. Under such approximation, we arrive at
the result:

$\displaystyle S^{\alpha}_{i}(\omega+i\delta)\approx$
$\displaystyle\;\int^{\infty}_{-\infty}\frac{d\epsilon}{2\pi
i}\big{[}n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\big{]}G_{\textrm{ret}}(\epsilon+\omega)G_{\textrm{adv}}(\epsilon)\gamma_{i}^{\alpha}(\epsilon+\omega+i\delta,\epsilon-i\delta).$
(61)

Substituting Eq. (61) into Eq. (55), we obtain
$\pi^{\alpha}_{ij}(\omega+i\delta,\bm{q})$ as

$\displaystyle\pi^{\alpha}_{ij}(\omega+i\delta,\bm{q})=$
$\displaystyle\;\frac{e\hbar^{3}}{m^{2}}\int\frac{d^{3}p}{(2\pi)^{3}}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}\int^{\infty}_{-\infty}\frac{d\epsilon}{2\pi
i}\big{[}n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\big{]}$ $\displaystyle\times
G_{\textrm{ret}}(\epsilon+\omega,\bm{p}+\bm{q})G_{\textrm{adv}}(\epsilon,\bm{p})\gamma_{i}^{\alpha}(\epsilon+\omega+i\delta,\epsilon-i\delta;\bm{p}+\bm{q},\bm{p}).$
(62)

Lastly, substituting this expression into the Kubo formula in Eq. (3), we
finally obtain Eq. (II) of the main text.

##  Appendix B Vertex function for spin current

We next compute the vertex function
$\gamma_{i}^{\alpha}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{q},\bm{p})$. In the
one-loop order, it is given by

$\displaystyle\gamma^{\alpha}_{i}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{q},\bm{p})=\frac{g^{2}\hbar^{2}}{M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})\mathcal{S}(ip_{m}+\omega_{n},ip_{m};\bm{p}+\bm{k}+\bm{q},\bm{p}+\bm{k}),$
(63)

where $\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$ and
$\mathcal{S}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{k}+\bm{q},\bm{p}+\bm{k})$
are defined as

$\displaystyle\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$
$\displaystyle=\sum_{a,b,c,d,e,f}\Big{[}\mathcal{P}_{ad}(\bm{k})(p+q+k/2)_{c}(p+k/2)_{f}\varepsilon_{abc}\varepsilon_{def}\textrm{tr}\big{(}\sigma_{\alpha}\sigma_{b}\sigma_{e}\big{)}\Big{]},$
(64)
$\displaystyle\mathcal{S}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{k}+\bm{q},\bm{p}+\bm{k})$
$\displaystyle=\frac{1}{\beta}\sum_{ik_{l}}\frac{2\omega_{\bm{k}}}{(ik_{l})^{2}-\omega_{\bm{k}}^{2}}\mathcal{G}(ip_{m}+i\omega_{n}+ik_{l},\bm{p}+\bm{k}+\bm{q})\mathcal{G}(ip_{m}+ik_{l},\bm{p}+\bm{k}).$
(65)

We first calculate $\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$ as

$\displaystyle\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})=$
$\displaystyle\;\sum_{a,b,c,d,e,f}\Big{[}\mathcal{P}_{ad}(\bm{k})(p+q+k/2)_{c}(p+k/2)_{f}\varepsilon_{abc}\varepsilon_{def}\textrm{tr}\big{(}\sigma_{\alpha}\sigma_{b}\sigma_{e}\big{)}\Big{]}$
$\displaystyle=$
$\displaystyle\;\sum_{a,b,c,d,e,f}\Big{[}\mathcal{P}_{ad}(\bm{k})(p+q+k/2)_{c}(p+k/2)_{f}\varepsilon_{abc}\varepsilon_{def}2i\varepsilon_{\alpha
be}\Big{]}$ $\displaystyle=$
$\displaystyle\;2i\sum_{a,c,d,f}\Big{[}\mathcal{P}_{ad}(\bm{k})(p+q+k/2)_{c}(p+k/2)_{f}(\delta_{a\alpha}\varepsilon_{dcf}-\delta_{c\alpha}\varepsilon_{daf})\Big{]}$
$\displaystyle=$ $\displaystyle\;2i\sum_{d}\mathcal{P}_{\alpha
d}(\bm{k})\big{[}(\bm{p}+\bm{q}+\bm{k}/2)\times(\bm{p}+\bm{k}/2)\big{]}_{d}.$
(66)

Summing over $d$, we obtain $\tilde{\gamma}^{\alpha}(\bm{p},\bm{p}+\bm{k})$ as

$\displaystyle\tilde{\gamma}^{\alpha}(\bm{p},\bm{p}+\bm{k})=2i\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}.$
(67)

We next turn to
$\mathcal{S}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{k}+\bm{q},\bm{p}+\bm{k})$ in
Eq. (65). For brevity, we omit the notations for the variables,
$\bm{p}+\bm{k}+\bm{q}$ and $\bm{p}+\bm{k}$, in
$\mathcal{S}(ip_{m}+i\omega_{n},ip_{m};\bm{p}+\bm{k}+\bm{q},\bm{p}+\bm{k})$,
$\mathcal{G}(ip_{m}+i\omega_{n}+ik_{l},\bm{p}+\bm{k}+\bm{q})$ and
$\mathcal{G}(ip_{m}+ik_{l},\bm{p}+\bm{k})$. In addition, we switch the
notation as $\omega_{\bm{k}}\rightarrow u$. Then, we convert the Matsubara
frequency summation over $ik_{n}$ into a contour integral as

$\displaystyle\mathcal{S}(ip_{m}+i\omega_{n},ip_{m})=\oint_{\gamma}\frac{dz}{2\pi
i}n_{B}(z)\frac{2u}{z^{2}-u^{2}}\mathcal{G}(ip_{m}+i\omega_{n}+z)\mathcal{G}(ip_{m}+z).$
(68)

The contour $\gamma$ expands to infinity in the complex plane of $z$ while
avoiding two isolated poles at $z=\pm u$ as well as two branch cuts along
$z=\epsilon^{\prime}-ip_{m}$ and $z=\epsilon^{\prime}-ip_{m}-i\omega_{n}$
where $\epsilon^{\prime}\in(-\infty,\infty)$. The contribution from the poles
is readily evaluated using the residue theorem, and the remaining integral
coming from the branch cuts is expressed as

$\displaystyle\mathcal{S}(ip_{m}+i\omega_{n},ip_{m})=$
$\displaystyle\;n_{B}(u)\mathcal{G}(ip_{m}+i\omega_{n}+u)\mathcal{G}(ip_{m}+u)+[n_{B}(u)+1]\mathcal{G}(ip_{m}+i\omega_{n}-u)\mathcal{G}(ip_{m}-u)$
$\displaystyle+\int^{\infty}_{-\infty}\frac{d\epsilon^{\prime}}{2\pi
i}n_{F}(\epsilon^{\prime})\Bigg{[}\frac{2u}{(\epsilon^{\prime}-ip_{m})^{2}-u^{2}}\mathcal{G}(\epsilon^{\prime}+i\omega_{n})\big{\\{}G_{\textrm{ret}}(\epsilon^{\prime})-G_{\textrm{adv}}(\epsilon^{\prime})\big{\\}}$
$\displaystyle+\frac{2u}{(\epsilon^{\prime}-ip_{m}-i\omega_{n})^{2}-u^{2}}\big{\\{}G_{\textrm{ret}}(\epsilon^{\prime})-G_{\textrm{adv}}(\epsilon^{\prime})\big{\\}}\mathcal{G}(\epsilon^{\prime}-i\omega_{n})\Bigg{]},$
(69)

where $n_{B}(\epsilon)$ is the Bose-Einstein distribution function. To obtain
$\gamma_{y}^{\alpha}(\epsilon+\omega+i\delta,\epsilon-i\delta)$, we perform
the analytic continuations: $i\omega_{n}\rightarrow\omega+i\delta$,
$ip_{m}\rightarrow\epsilon-i\delta$, and
$ip_{m}+i\omega_{n}\rightarrow\epsilon+\omega+i\delta$. As a result, we obtain

$\displaystyle\mathcal{S}(\epsilon+\omega+i\delta,\epsilon-i\delta)=$
$\displaystyle\;n_{B}(u)G_{\textrm{ret}}(\epsilon+\omega+u)G_{\textrm{adv}}(\epsilon+u)+[n_{B}(u)+1]G_{\textrm{ret}}(\epsilon+\omega-u)G_{\textrm{adv}}(\epsilon-u)$
$\displaystyle+\int^{\infty}_{-\infty}\frac{d\epsilon^{\prime}}{2\pi
i}n_{F}(\epsilon^{\prime})\Bigg{[}\frac{2u}{(\epsilon^{\prime}-\epsilon+i\delta)^{2}-u^{2}}G_{\textrm{ret}}(\epsilon^{\prime}+\omega)\big{\\{}G_{\textrm{ret}}(\epsilon^{\prime})-G_{\textrm{adv}}(\epsilon^{\prime})\big{\\}}$
$\displaystyle+\frac{2u}{(\epsilon^{\prime}-\epsilon-\omega-i\delta)^{2}-u^{2}}\big{\\{}G_{\textrm{ret}}(\epsilon^{\prime})-G_{\textrm{adv}}(\epsilon^{\prime})\big{\\}}G_{\textrm{adv}}(\epsilon^{\prime}-\omega)\Bigg{]}.$
(70)

We perform the integral over $\epsilon^{\prime}$ in the above equation
utilizing the following identity:

$\displaystyle\lim_{\delta\rightarrow
0}\frac{2u}{(\epsilon^{\prime}-\epsilon+i\delta)^{2}-u^{2}}$
$\displaystyle=\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon-u}\bigg{)}-\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon+u}\bigg{)}-i\pi\delta(\epsilon^{\prime}-\epsilon-u)+i\pi\delta(\epsilon^{\prime}-\epsilon+u),$
(71) $\displaystyle\lim_{\delta\rightarrow
0}\frac{2u}{(\epsilon^{\prime}-\epsilon-\omega-i\delta)^{2}-u^{2}}$
$\displaystyle=\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon-\omega-u}\bigg{)}-\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon-\omega+u}\bigg{)}+i\pi\delta(\epsilon^{\prime}-\epsilon-\omega-u)-i\pi\delta(\epsilon^{\prime}-\epsilon-\omega+u),$
(72)

where $\textrm{P}(x)$ denotes the principal part of $x$. As a result, we
obtain

$\displaystyle\mathcal{S}(\epsilon+\omega+i\delta,\epsilon-i\delta)=$
$\displaystyle\;n_{B}(u)G_{\textrm{ret}}(\epsilon+\omega+u)G_{\textrm{adv}}(\epsilon+u)+\big{\\{}n_{B}(u)+1\big{\\}}G_{\textrm{ret}}(\epsilon+\omega-u)G_{\textrm{adv}}(\epsilon-u)$
$\displaystyle-\frac{1}{2}n_{F}(\epsilon+u)G_{\textrm{ret}}(\epsilon+\omega+u)\big{\\{}G_{\textrm{ret}}(\epsilon+u)-G_{\textrm{adv}}(\epsilon+u)\big{\\}}$
$\displaystyle+\frac{1}{2}n_{F}(\epsilon-u)G_{\textrm{ret}}(\epsilon+\omega-u)\big{\\{}G_{\textrm{ret}}(\epsilon-u)-G_{\textrm{adv}}(\epsilon-u)\big{\\}}$
$\displaystyle+\frac{1}{2}n_{F}(\epsilon+\omega+u)\big{\\{}G_{\textrm{ret}}(\epsilon+\omega+u)-G_{\textrm{adv}}(\epsilon+\omega+u)\big{\\}}G_{\textrm{adv}}(\epsilon+u)$
$\displaystyle-\frac{1}{2}n_{F}(\epsilon+\omega-u)\big{\\{}G_{\textrm{ret}}(\epsilon+\omega-u)-G_{\textrm{adv}}(\epsilon+\omega-u)\big{\\}}G_{\textrm{adv}}(\epsilon-u)$
$\displaystyle+\tilde{P}(\epsilon+\omega,\epsilon),$ (73)

where $\tilde{P}(\epsilon+\omega,\epsilon)$ represents the contribution from
the principal part, given by

$\displaystyle P(\epsilon+\omega,\epsilon)=$
$\displaystyle\;\int^{\infty}_{-\infty}\frac{d\epsilon^{\prime}}{2\pi
i}n_{F}(\epsilon^{\prime})\Bigg{[}\bigg{\\{}\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon-u}\bigg{)}-\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon+u}\bigg{)}\bigg{\\}}G_{\textrm{ret}}(\epsilon^{\prime}+\omega)\big{\\{}G_{\textrm{ret}}(\epsilon^{\prime})-G_{\textrm{adv}}(\epsilon^{\prime})\big{\\}}$
$\displaystyle+\bigg{\\{}\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon-\omega-u}\bigg{)}-\textrm{P}\bigg{(}\frac{1}{\epsilon^{\prime}-\epsilon-\omega+u}\bigg{)}\bigg{\\}}\big{\\{}G_{\textrm{ret}}(\epsilon^{\prime})-G_{\textrm{adv}}(\epsilon^{\prime})\big{\\}}G_{\textrm{adv}}(\epsilon^{\prime}-\omega)\Bigg{]}.$
(74)

In $\mathcal{S}(\epsilon+\omega+i\delta,\epsilon-i\delta)$, we drop the terms
that are given as the products of two retarded propagators or two advanced
propagators as well as the contribution from the principal part since their
contribution is negligible when $\omega\ll E_{F}$ [38]. As a result, we obtain

$\displaystyle\mathcal{S}(\epsilon+\omega+i\delta,\epsilon-i\delta)\approx$
$\displaystyle\;\bigg{\\{}n_{B}(u)+\frac{1}{2}n_{F}(\epsilon+u)+\frac{1}{2}n_{F}(\epsilon+\omega+u)\bigg{\\}}G_{\textrm{ret}}(\epsilon+\omega+u)G_{\textrm{adv}}(\epsilon+u)$
$\displaystyle+\bigg{\\{}n_{B}(u)+1-\frac{1}{2}n_{F}(\epsilon-u)-\frac{1}{2}n_{F}(\epsilon+\omega-u)\bigg{\\}}G_{\textrm{ret}}(\epsilon+\omega-u)G_{\textrm{adv}}(\epsilon-u).$
(75)

Substituting Eqs. (67) and (B) into Eq. (63), we finally obtain

$\displaystyle\gamma^{\alpha}_{i}(\epsilon+\omega+i\delta,\epsilon-i\delta)\approx$
$\displaystyle\;\frac{2ig^{2}\hbar^{2}}{M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}$
$\displaystyle\times\Bigg{[}\bigg{\\{}n_{B}(u)+\frac{1}{2}n_{F}(\epsilon+u)+\frac{1}{2}n_{F}(\epsilon+\omega+u)\bigg{\\}}G_{\textrm{ret}}(\epsilon+\omega+u)G_{\textrm{adv}}(\epsilon+u)$
$\displaystyle+\bigg{\\{}n_{B}(u)+1-\frac{1}{2}n_{F}(\epsilon-u)-\frac{1}{2}n_{F}(\epsilon+\omega-u)\bigg{\\}}G_{\textrm{ret}}(\epsilon+\omega-u)G_{\textrm{adv}}(\epsilon-u)\bigg{]}.$
(76)

##  Appendix C Conductivity in the one-loop order

We now compute the spin conductivity using the formulae derived in the
preceding sections. By substituting the vertex function from Eq. (B) into Eq.
(A), we obtain the correlation function as

$\displaystyle\pi_{ij,\textrm{ret}}^{\alpha}(\omega,\bm{q})=$
$\displaystyle\;\frac{1}{\pi}\frac{g^{2}e\hbar^{5}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\int^{\infty}_{-\infty}d\epsilon\\{n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\\}G_{\textrm{ret}}(\epsilon+\omega,\bm{p}+\bm{q})G_{\textrm{adv}}(\epsilon,\bm{p})$
$\displaystyle\times\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}$
$\displaystyle\times\Bigg{[}\bigg{\\{}n_{B}(\omega_{\bm{k}})+\frac{1}{2}n_{F}(\epsilon+\omega_{\bm{k}})+\frac{1}{2}n_{F}(\epsilon+\omega+\omega_{\bm{k}})\bigg{\\}}G_{\textrm{ret}}(\epsilon+\omega+\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})G_{\textrm{adv}}(\epsilon+\omega_{\bm{k}},\bm{p}+\bm{k})$
$\displaystyle+\bigg{\\{}n_{B}(\omega_{\bm{k}})+1-\frac{1}{2}n_{F}(\epsilon-\omega_{\bm{k}})-\frac{1}{2}n_{F}(\epsilon+\omega-\omega_{\bm{k}})\bigg{\\}}G_{\textrm{ret}}(\epsilon+\omega-\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})G_{\textrm{adv}}(\epsilon-\omega_{\bm{k}},\bm{p}+\bm{k}).$
(77)

Substituting Eq. (C) into the Kobu formula in Eq. (3), we finally obtain the
spin conductivity in the one-loop order as

$\displaystyle\sigma_{ij}^{\alpha}(\omega,\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\int^{\infty}_{-\infty}d\epsilon\frac{1}{\omega}\\{n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\\}K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q}),$
(78)

where $K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$ is defined as

$\displaystyle\;K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$
$\displaystyle=$
$\displaystyle\;G_{\textrm{ret}}(\epsilon+\omega,\bm{p}+\bm{q})G_{\textrm{adv}}(\epsilon,\bm{p})\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}$
$\displaystyle\times\Bigg{[}\bigg{\\{}n_{B}(\omega_{\bm{k}})+\frac{1}{2}n_{F}(\epsilon+\omega_{\bm{k}})+\frac{1}{2}n_{F}(\epsilon+\omega+\omega_{\bm{k}})\bigg{\\}}G_{\textrm{adv}}(\epsilon+\omega_{\bm{k}},\bm{p}+\bm{k})G_{\textrm{ret}}(\epsilon+\omega+\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})$
$\displaystyle+\bigg{\\{}n_{B}(\omega_{\bm{k}})+1-\frac{1}{2}n_{F}(\epsilon-\omega_{\bm{k}})-\frac{1}{2}n_{F}(\epsilon+\omega-\omega_{\bm{k}})\bigg{\\}}G_{\textrm{adv}}(\epsilon-\omega_{\bm{k}},\bm{p}+\bm{k})G_{\textrm{ret}}(\epsilon+\omega-\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})\Bigg{]}.$
(79)

To determine the dominant dynamic behavior and temperature dependence, we
utilize a double expansion in the integral over $\epsilon$: a Taylor expansion
in $\omega$ around $\omega=0$ and a Sommerfeld expansion in $k_{B}T$ at $T=0$.
Carrying out this expansion on
$\int^{\infty}_{-\infty}d\epsilon\frac{1}{\omega}\big{\\{}n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\big{\\}}K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$,
we obtain:

$\displaystyle\int^{\infty}_{-\infty}d\epsilon\frac{1}{\omega}\big{\\{}n_{F}(\epsilon)-n_{F}(\epsilon+\omega)\big{\\}}K_{ij}^{\alpha}(\epsilon;\omega)$
(80) $\displaystyle\approx$
$\displaystyle\frac{1}{\omega}\int^{0}_{-\omega}d\epsilon
K_{ij}^{\alpha}(\epsilon;\omega)+\frac{\pi^{2}}{6}(k_{B}T)^{2}\frac{1}{\omega}\left(\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}-\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=-\omega}\right)K_{ij}^{\alpha}(\epsilon;\omega)$
$\displaystyle\approx$ $\displaystyle
K_{ij}^{\alpha}(0;\omega)-\frac{1}{2}\omega\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}K_{ij}^{\alpha}(\epsilon;\omega)+\frac{\pi^{2}}{6}(k_{B}T)^{2}\left.\frac{\partial^{2}}{\partial\epsilon^{2}}\right|_{\epsilon=0}K_{ij}^{\alpha}(\epsilon;\omega)$
$\displaystyle\approx$ $\displaystyle
K_{ij}^{\alpha}(0;0)+\omega\left[\left.\frac{\partial}{\partial\omega}\right|_{\omega=0}K_{ij}^{\alpha}(0;\omega)-\frac{1}{2}\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}K_{ij}^{\alpha}(\epsilon;0)\right]+\frac{\pi^{2}}{6}(k_{B}T)^{2}\left.\frac{\partial^{2}}{\partial\epsilon^{2}}\right|_{\epsilon=0}K_{ij}^{\alpha}(\epsilon;0),$

where $K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$ has been written
as $K_{ij}^{\alpha}(\epsilon;\omega)$ for brevity. Substituting this
expression into Eq. (78), we obtain an expanded form of
$\sigma^{\alpha}_{ij}(\omega,\bm{q})$:

$\displaystyle\sigma^{\alpha}_{ij}(\omega,\bm{q})=\sigma^{\alpha}_{ij}(\bm{q})+\omega\sigma^{\alpha,(1)}_{ij}(\bm{q})+(k_{B}T)^{2}\sigma^{\alpha,(2)}_{ij}(\bm{q})+\mathcal{O}(\omega^{2}),$
(81)

where $\sigma^{\alpha}_{ij}(\bm{q})$, $\sigma^{\alpha,(1)}_{ij}(\bm{q})$, and
$\sigma^{\alpha,(2)}_{ij}(\bm{q})$ are given by

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=$
$\displaystyle\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}K_{ij}^{\alpha}(0,\bm{p},\bm{k};0,\bm{q}),$
(82) $\displaystyle\sigma^{\alpha,(1)}_{ij}(\bm{q})=$
$\displaystyle\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\left[\frac{\partial}{\partial\omega}K_{ij}^{\alpha}(0,\bm{p},\bm{k};0,\bm{q})-\frac{1}{2}\frac{\partial}{\partial\epsilon}K_{ij}^{\alpha}(0,\bm{p},\bm{k};0,\bm{q})\right],$
(83) $\displaystyle\sigma^{\alpha,(2)}_{ij}(\bm{q})=$
$\displaystyle\frac{\pi^{2}}{6}\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{2\omega_{\bm{k}}}\frac{\partial^{2}}{\partial\epsilon^{2}}K_{ij}^{\alpha}(0,\bm{p},\bm{k};0,\bm{q}).$
(84)

##  Appendix D DC conductivity

In this section, we derive the DC conductivity in Eq. (17) in the main text.
It is defined in Eq. (82), and is expressed here as

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}$
$\displaystyle\times\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{\alpha}\cdot\Big{\\{}[\bm{q}\times(\bm{p}+\bm{k}/2)]-\hat{\bm{k}}\big{(}\hat{\bm{k}}\cdot[\bm{q}\times(\bm{p}+\bm{k}/2)]\big{)}\Big{\\}}$
$\displaystyle\times\Big{[}G_{\textrm{ret}}(0,\bm{p}+\bm{q})G_{\textrm{adv}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})G_{\textrm{adv}}(\omega_{\bm{k}},\bm{p}+\bm{k})+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}})\Big{]}.$
(85)

The second line can be expanded as:

$\displaystyle\hat{\bm{e}}_{i}\cdot\bigg{(}\bm{p}+\bm{k}+\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{j}\cdot\bigg{(}\bm{p}-\frac{1}{2}\bm{q}\bigg{)}\hat{\bm{e}}_{\alpha}\cdot\Big{\\{}[\bm{q}\times(\bm{p}+\bm{k}/2)]-\hat{\bm{k}}\big{(}\hat{\bm{k}}\cdot[\bm{q}\times(\bm{p}+\bm{k}/2)]\big{)}\Big{\\}}$
$\displaystyle=$
$\displaystyle\bigg{(}\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{p})+\frac{1}{2}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{p})-\frac{1}{2}\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{q})-\frac{1}{4}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{q})\bigg{)}$
$\displaystyle\times\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}.$
(86)

Additionally, we expand the product of the electron propagators in the third
line in $\bm{q}$ as:

$\displaystyle
G_{\textrm{ret}}(0,\bm{p}+\bm{q})G_{\textrm{adv}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k}+\bm{q})G_{\textrm{adv}}(\omega_{\bm{k}},\bm{p}+\bm{k})+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}})$
$\displaystyle=$
$\displaystyle|G_{\textrm{ret}}(0,\bm{p})|^{2}|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})|^{2}\bigg{\\{}1+\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})$
$\displaystyle+\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{2}[G_{\textrm{ret}}(0,\bm{p})]^{2}+\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{2}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}$
$\displaystyle+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})$
$\displaystyle+2\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}[G_{\textrm{ret}}(0,\bm{p})]^{2}+2\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}+\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})$
$\displaystyle+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{3}[G_{\textrm{ret}}(0,\bm{p})]^{3}+\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{3}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{3}$
$\displaystyle+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{2}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}[G_{\textrm{ret}}(0,\bm{p})]^{2}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{2}G_{\textrm{ret}}(0,\bm{p})[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}$
$\displaystyle+\mathcal{O}(q^{4})\bigg{\\}}+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}}).$
(87)

Utilizing Eqs. (86) and (87), we expand $\sigma^{\alpha}_{ij}(\bm{q})$ in
$\bm{q}$ as:

$\displaystyle\sigma^{\alpha}_{ij}(\bm{q})=[\sigma^{\alpha}_{ij}(\bm{q})]^{(1)}+[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}+[\sigma^{\alpha}_{ij}(\bm{q})]^{(3)}+[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}+\mathcal{O}(q^{5}),$
(88)

where $[\sigma^{\alpha}_{ij}(\bm{q})]^{(o)}$ denotes the $o$-th order term in
the expansion. The odd terms are zero since $\sigma^{\alpha}_{ij}(\bm{q})$ is
symmetric in $\bm{q}\rightarrow-\bm{q}$. The quadratic term is given by

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}$
$\displaystyle\times\Bigg{[}\Bigg{\\{}|G_{\textrm{ret}}(0,\bm{p})|^{2}|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})|^{2}\Bigg{(}\frac{1}{2}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{p})-\frac{1}{2}\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{q})$
$\displaystyle+\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{p})\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})\bigg{)}\Bigg{)}\Bigg{\\}}+\Bigg{\\{}\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}}\Bigg{\\}}\Bigg{]}.$
(89)

The quartic order term is given by

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}$ $\displaystyle=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\Big{\\{}\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{e}}_{\alpha}\big{]}-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot\big{[}(\bm{p}+\bm{k}/2)\times\hat{\bm{k}}\big{]}\Big{\\}}$
$\displaystyle\times|G_{\textrm{ret}}(0,\bm{p})|^{2}|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})|^{2}\Bigg{[}\Bigg{\\{}\bigg{(}-\frac{1}{4}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{q})\bigg{)}\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}G_{\textrm{ret}}(0,\bm{p})+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})\bigg{)}\Bigg{\\}}$
$\displaystyle+\Bigg{\\{}\bigg{(}\frac{1}{2}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{p})-\frac{1}{2}\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{q})\bigg{)}\bigg{(}\frac{\hbar^{2}\bm{q}^{2}}{2m}[G_{\textrm{ret}}(0,\bm{p})+G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{2}[G_{\textrm{ret}}(0,\bm{p})]^{2}$
$\displaystyle+\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{2}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})\bigg{)}\Bigg{\\}}$
$\displaystyle+\Bigg{\\{}\hat{\bm{e}}_{i}\cdot(\bm{p}+\bm{k})(\hat{\bm{e}}_{j}\cdot\bm{p})\Bigg{(}2\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}[G_{\textrm{ret}}(0,\bm{p})]^{2}+2\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}+\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})$
$\displaystyle+\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\frac{\hbar^{2}\bm{q}^{2}}{2m}G_{\textrm{ret}}(0,\bm{p})G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{3}[G_{\textrm{ret}}(0,\bm{p})]^{3}+\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{3}[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{3}$
$\displaystyle+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}^{2}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}[G_{\textrm{ret}}(0,\bm{p})]^{2}G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})+\bigg{(}\frac{\hbar^{2}\bm{p}\cdot\bm{q}}{m}\bigg{)}\bigg{(}\frac{\hbar^{2}(\bm{p}+\bm{k})\cdot\bm{q}}{m}\bigg{)}^{2}G_{\textrm{ret}}(0,\bm{p})[G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p}+\bm{k})]^{2}\Bigg{)}\Bigg{\\}}\Bigg{]}$
$\displaystyle+\Bigg{[}\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}}\Bigg{]}.$
(90)

We first calculate $[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$. In the integral,
the factor $n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})$, as it decays
exponentially for sufficiently large $|\bm{k}|$, effectively constrains
$|\bm{k}|\lesssim k_{\textrm{cutoff}}$. In the regime of our interest, where
$T\ll T_{\textrm{k}}\equiv\frac{\hbar ck_{F}}{k_{B}}$, we have
$|\bm{k}|\lesssim k_{\textrm{cutoff}}\ll k_{F}$. This motivates us to neglect
the $\bm{k}$-dependence in the integrand except for
$\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{\omega_{\bm{k}}}$, the
denominator necessary for obtaining the correct behavior in the limit where
the phonon energy gap vanishes, i.e. $E_{g}\to 0$. Furthermore, we observe
that the integral over $\bm{p}$ sharply peaks around $|\bm{p}|\approx k_{F}$
when $E_{F}\tau\gg\hbar$. Consequently, we can approximate the propagator as
$G_{\textrm{ret}}(0,\bm{p})\approx G_{\textrm{ret}}(0,p)=(-\hbar
v_{F}p+i\hbar/2\tau)^{-1}$. Here, we introduce the shorthand notation
$p\equiv|\bm{p}|-k_{F}$. For the same reason, we replace $\bm{p}$ with
$k_{F}\hat{\bm{p}}$ in the other factors involving $\bm{p}$ in Eq. (D). As a
result, the integral over $\bm{p}$ can be approximated as: $\int
d^{3}p\approx\int^{\infty}_{-\infty}dpk_{F}^{2}\int d\Omega_{p}$. Within this
approximation, we obtain a more concise expression for
$[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$ as:

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{k_{F}^{4}q^{2}}{(2\pi)^{3}}[\zeta_{ij}^{\alpha}(\bm{q})-\zeta_{ji}^{\alpha}(\bm{q})]\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{4\hbar^{2}k_{F}^{6}q^{2}}{(2\pi)^{3}m}\eta_{ij}^{\alpha}(\bm{q})\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p).$
(91)

Here, $\zeta_{ij}^{\alpha}(\bm{q})$ and $\eta_{ij}^{\alpha}(\bm{q})$ are
defined as

$\displaystyle\zeta_{ij}^{\alpha}(\bm{q})=$ $\displaystyle\int
d\Omega_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}=\frac{4\pi}{3}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{]},$
(92) $\displaystyle\eta_{ij}^{\alpha}(\bm{q})=$ $\displaystyle\int
d\Omega_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{p}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}(\hat{\bm{p}}\cdot\hat{\bm{q}})=\frac{4\pi}{15}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{]}\big{\\}}+\\{i\leftrightarrow
j\\}\Big{]},$ (93)

where we have denoted $\hat{\bm{p}}=\bm{p}/|\bm{p}|$ and
$\hat{\bm{q}}=\bm{q}/|\bm{q}|$. The remaining integrals for $p$ are evaluated
as follows:

$\displaystyle\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}$
$\displaystyle=\;\int^{\infty}_{-\infty}dp\frac{1}{[(\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}]^{2}}=\frac{4\pi}{\hbar v_{F}(\hbar/\tau)^{3}},$
(94)
$\displaystyle\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)$
$\displaystyle=\;\int^{\infty}_{-\infty}dp\frac{1}{[(\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}]^{2}}\frac{1}{-\hbar
v_{F}p+i(\hbar/2\tau)}=\frac{-6\pi i}{\hbar v_{F}(\hbar/\tau)^{4}}.$ (95)

Additionally, we represent the integral over $\bm{k}$ as $\int
d^{3}k=\int^{\infty}_{0}dkk^{2}\int d\Omega_{k}$. The integral over the
angular variables, denoted as $\Omega_{k}$, is evaluated as

$\displaystyle\int d\Omega_{k}(1-\hat{k}_{\alpha}^{2})=\frac{8\pi}{3}.$ (96)

Substituting all the obtained integral expressions, from Eq. (92) to Eq. (96),
into Eq. (D), we derive the following expression for
$[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$:

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\big{[}-\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}$
$\displaystyle\times\int^{\infty}_{0}dk\Bigg{[}\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{k_{F}^{4}}{(2\pi)^{6}}\frac{32\pi^{2}}{9}\frac{4\pi}{\hbar
v_{F}(\hbar/\tau)^{3}}\Bigg{]}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\big{[}\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}$
$\displaystyle\times\int^{\infty}_{0}dk\Bigg{[}\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{k_{F}^{6}}{(2\pi)^{6}}\frac{4\hbar^{2}}{m}\frac{32\pi^{2}}{45}\frac{-6\pi
i}{\hbar v_{F}(\hbar/\tau)^{4}}\Bigg{]}.$ (97)

This expression can be represented as Eq. (17) in the main text, utilizing the
scaled expression for the integral over $k$ provided in Eq. (21) in the main
text, along with the notations $n=\frac{k_{F}^{3}}{3\pi^{2}}$ and
$E_{F}=\frac{\hbar^{2}k_{F}^{2}}{2m}$.

We next calculate $[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}$. The calculation is
similar with $[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$. Given $|\bm{k}|\lesssim
k_{\textrm{cutoff}}\ll k_{F}$, $[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}$ is
represented as:

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\Big{\\{}\bm{q}\cdot(\bm{p}\times\hat{\bm{e}}_{\alpha})-(\hat{\bm{e}}_{\alpha}\cdot\hat{\bm{k}})\bm{q}\cdot(\bm{p}\times\hat{\bm{k}})\Big{\\}}$
$\displaystyle\times\Bigg{[}\bigg{(}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{p})-(\hat{\bm{e}}_{i}\cdot\bm{p})(\hat{\bm{e}}_{j}\cdot\bm{q})\bigg{)}\bigg{(}\frac{\hbar^{2}\bm{q}^{2}}{m}\bigg{)}|G_{\textrm{ret}}(0,\bm{p})|^{4}G_{\textrm{ret}}(0,\bm{p})$
$\displaystyle+\bigg{(}(\hat{\bm{e}}_{i}\cdot\bm{q})(\hat{\bm{e}}_{j}\cdot\bm{p})-(\hat{\bm{e}}_{i}\cdot\bm{p})(\hat{\bm{e}}_{j}\cdot\bm{q})\bigg{)}\bigg{(}3\frac{\hbar^{4}(\bm{p}\cdot\bm{q})^{2}}{m^{2}}\bigg{)}|G_{\textrm{ret}}(0,\bm{p})|^{4}G_{\textrm{ret}}(0,\bm{p})^{2}$
$\displaystyle+(\hat{\bm{e}}_{i}\cdot\bm{p})(\hat{\bm{e}}_{j}\cdot\bm{p})\bigg{(}6\frac{\hbar^{4}\bm{q}^{2}(\bm{p}\cdot\bm{q})}{m^{2}}\bigg{)}\bigg{\\}}|G_{\textrm{ret}}(0,\bm{p})|^{4}G_{\textrm{ret}}(0,\bm{p})^{2}$
$\displaystyle+(\hat{\bm{e}}_{i}\cdot\bm{p})(\hat{\bm{e}}_{j}\cdot\bm{p})\bigg{(}8\frac{\hbar^{6}(\bm{p}\cdot\bm{q})^{3}}{m^{3}}\bigg{)}|G_{\textrm{ret}}(0,\bm{p})|^{4}G_{\textrm{ret}}(0,\bm{p})^{3}\Bigg{]}.$
(98)

Given $E_{F}\tau\gg\hbar$, this can be further simplified as:

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{k_{F}^{4}}{(2\pi)^{3}}\frac{\hbar^{2}q^{4}}{m}[\zeta_{ij}^{\alpha}(\bm{q})-\zeta_{ji}^{\alpha}(\bm{q})]\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{k_{F}^{6}}{(2\pi)^{3}}\frac{3\hbar^{4}q^{4}}{m^{2}}[\tilde{\zeta}_{ij}^{\alpha}(\bm{q})-\tilde{\zeta}_{ji}^{\alpha}(\bm{q})]\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)^{2}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{k_{F}^{6}}{(2\pi)^{3}}\frac{6\hbar^{4}q^{4}}{m^{2}}\eta_{ij}^{\alpha}(\bm{q})\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)^{2}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{k_{F}^{8}}{(2\pi)^{3}}\frac{8\hbar^{6}q^{4}}{m^{3}}\tilde{\eta}_{ij}^{\alpha}(\bm{q})\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)^{3}.$
(99)

Here, $\tilde{\zeta}_{ij}^{\alpha}(\bm{q})$ and
$\tilde{\eta}_{ij}^{\alpha}(\bm{q})$ are defined as

$\displaystyle\tilde{\zeta}_{ij}^{\alpha}(\bm{q})=$ $\displaystyle\int
d\Omega_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}(\hat{\bm{p}}\cdot\hat{\bm{q}})^{2}=\frac{4\pi}{15}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})=\frac{1}{5}\zeta_{ij}^{\alpha}(\bm{q}),$
(100) $\displaystyle\tilde{\eta}_{ij}^{\alpha}(\bm{q})=$ $\displaystyle\int
d\Omega_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{p}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}(\hat{\bm{p}}\cdot\hat{\bm{q}})^{3}=\frac{4\pi}{35}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{]}\big{\\}}+\\{i\leftrightarrow
j\\}\Big{]}=\frac{3}{7}\eta_{ij}^{\alpha}(\bm{q}).$ (101)

The remaining integrals for $p$ are evaluated as follows:

$\displaystyle\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)^{2}$
$\displaystyle=\;\int^{\infty}_{-\infty}dp\frac{1}{[(\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}]^{2}}\frac{1}{[-\hbar
v_{F}p+i(\hbar/2\tau)]^{2}}=\frac{-8\pi}{\hbar v_{F}(\hbar/\tau)^{5}},$ (102)
$\displaystyle\int^{\infty}_{-\infty}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p)^{2}$
$\displaystyle=\;\int^{\infty}_{-\infty}dp\frac{1}{[(\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}]^{2}}\frac{1}{[-\hbar
v_{F}p+i(\hbar/2\tau)]^{3}}=\frac{-10\pi i}{\hbar v_{F}(\hbar/\tau)^{6}}.$
(103)

Substituting all the obtained integral expressions, Eqs. from (92) to (96) and
Eqs. from (100) to (103), into Eq. (D), we derive the following expression for
$[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}$:

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\big{[}-\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}$
$\displaystyle\times\int^{\infty}_{0}dk\Bigg{[}\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{k_{F}^{4}}{(2\pi)^{6}}\frac{\hbar^{2}q^{2}}{m}\frac{32\pi^{2}}{9}\frac{-6\pi
i}{\hbar v_{F}(\hbar/\tau)^{4}}\Bigg{]}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\big{[}-\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}$
$\displaystyle\times\int^{\infty}_{0}dk\Bigg{[}\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{k_{F}^{6}}{(2\pi)^{6}}\frac{3\hbar^{4}q^{2}}{m^{2}}\frac{32\pi^{2}}{45}\frac{-8\pi}{\hbar
v_{F}(\hbar/\tau)^{5}}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\big{[}\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}$
$\displaystyle\times\int^{\infty}_{0}dk\Bigg{[}\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{k_{F}^{6}}{(2\pi)^{6}}\frac{6\hbar^{4}q^{2}}{m^{2}}\frac{32\pi^{2}}{45}\frac{-8\pi}{\hbar
v_{F}(\hbar/\tau)^{5}}\Bigg{]}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\big{[}\bm{q}\cdot(\hat{\bm{e}}_{i}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{j})+\bm{q}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})(\bm{q}\cdot\hat{\bm{e}}_{i})\big{]}$
$\displaystyle\times\int^{\infty}_{0}dk\Bigg{[}\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{k_{F}^{8}}{(2\pi)^{6}}\frac{8\hbar^{6}q^{2}}{m^{3}}\frac{32\pi^{2}}{105}\frac{-10\pi
i}{\hbar v_{F}(\hbar/\tau)^{6}}\Bigg{]}.$ (104)

This is rewritten as:

$\displaystyle[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}=$
$\displaystyle\;\chi_{0}\frac{5q^{2}}{4k_{F}^{2}}\big{[}-\nu_{ij}^{\alpha}(\bm{q})+\nu_{ij}^{\alpha}(\bm{q})\big{]}$
$\displaystyle-i\kappa_{0}\frac{24q^{2}}{5k_{F}^{2}}\bigg{(}\frac{E_{F}\tau}{\hbar}\bigg{)}^{2}\big{[}-\nu_{ij}^{\alpha}(\bm{q})+\nu_{ij}^{\alpha}(\bm{q})\big{]}$
$\displaystyle-i\kappa_{0}\frac{48q^{2}}{5k_{F}^{2}}\bigg{(}\frac{E_{F}\tau}{\hbar}\bigg{)}^{2}\big{[}\nu_{ij}^{\alpha}(\bm{q})+\nu_{ij}^{\alpha}(\bm{q})\big{]}$
$\displaystyle+\chi_{0}\frac{40q^{2}}{7k_{F}^{2}}\bigg{(}\frac{E_{F}\tau}{\hbar}\bigg{)}^{2}\big{[}\nu_{ij}^{\alpha}(\bm{q})+\nu_{ij}^{\alpha}(\bm{q})\big{]}.$
(105)

When $q\ll k_{F}$, the first term is negligible compared to the real part of
$[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$. Moreover, if
$\frac{48q^{2}}{5k_{F}^{2}}\big{(}\frac{E_{F}\tau}{\hbar}\big{)}^{2}\ll 1$ or
$q\ll\sqrt{5/48}\frac{\hbar}{E_{F}\tau}k_{F}\approx
0.32\frac{\hbar}{E_{F}\tau}k_{F}$, the last term is negligible compared to the
real part of $[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$ and the second and third
terms are negligible compared to the imaginary part of
$[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$. Since typical metals manifest
$\hbar/(E_{F}\tau)\ll 1$, the latter condition implies the former condition.
Consequently, we conclude that $[\sigma^{\alpha}_{ij}(\bm{q})]^{(4)}$ is
negligible compared to $[\sigma^{\alpha}_{ij}(\bm{q})]^{(2)}$ when $q\ll
0.32\frac{\hbar}{E_{F}\tau}k_{F}$.

##  Appendix E Dynamic conductivity

The dynamic spin conductivity up to the first order in frequency derived from
Eq. (80) gives us the Drude form for the spin conductivity only in the limit
where the ratio of the phonon frequency to the scattering rate vanishes. The
starting point would be to carry out the first-order differentiation in Eq.
(80) where we drop any dependence on ${\bf q}$ and ${\bm{k}}$ except for
$W^{\alpha}_{ij}({\bf p},{\bm{k}};{\bf q})$ and $\omega_{\bm{k}}$ for
simplicity:

$\displaystyle\left.\frac{\partial}{\partial\omega}\right|_{\omega=0}K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$
$\displaystyle=W^{\alpha}_{ij}({\bf p},{\bm{k}};{\bf
q})G_{\textrm{adv}}(0,\bm{p})$
$\displaystyle\times\Bigg{[}\\{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})\\}\bigg{\\{}G_{\textrm{adv}}(\omega_{\bm{k}},\bm{p})\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}\\{G_{\textrm{ret}}(\epsilon,\bm{p})G_{\textrm{ret}}(\epsilon+\omega_{\bm{k}},\bm{p})\\}+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}})\bigg{\\}}$
$\displaystyle+\frac{1}{2}G_{\textrm{ret}}(0,\bm{p})\\{|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p})|^{2}-|G_{\textrm{ret}}(-\omega_{\bm{k}},\bm{p})|^{2}\\}\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}n_{F}(\epsilon+\omega_{\bm{k}})\Bigg{]},$
$\displaystyle\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}K_{ij}^{\alpha}(\epsilon,\bm{p},\bm{k};\omega,\bm{q})$
$\displaystyle=W^{\alpha}_{ij}({\bf p},{\bm{k}};{\bf
q})\\{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})\\}$
$\displaystyle\times\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}\Bigg{[}|G_{\textrm{ret}}(\epsilon,\bm{p})|^{2}\\{|G_{\textrm{ret}}(\epsilon+\omega_{\bm{k}},\bm{p})|^{2}+|G_{\textrm{ret}}(\epsilon-\omega_{\bm{k}},\bm{p})|^{2}\\}\Bigg{]}$
$\displaystyle+W^{\alpha}_{ij}({\bf p},{\bm{k}};{\bf
q})|G_{\textrm{ret}}(0,\bm{p})|^{2}\\{|G_{\textrm{ret}}(\omega_{\bm{k}},\bm{p})|^{2}-|G_{\textrm{ret}}(-\omega_{\bm{k}},\bm{p})|^{2}\\}\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}n_{F}(\epsilon+\omega_{\bm{k}}).$
(106)

At low temperatures, we can clearly ignore the terms involving the derivative
of the Fermi-Dirac distribution function. In addition, it is fairly
straightforward to see that the integral $\int
p^{2}dp|G_{\textrm{ret}}(\epsilon,\bm{p})|^{2}|G_{\textrm{ret}}(\epsilon\pm\omega_{\bm{k}},\bm{p})|^{2}$
does not depend on $\epsilon$. We therefore obtain

$\displaystyle\sigma^{\alpha,(1)}_{ij}(\bm{q})=$
$\displaystyle\frac{ig^{2}}{\pi}\frac{e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}k}{(2\pi)^{3}}W^{\alpha}_{ij}({\bm{p}},{\bm{k}};{\bm{q}})G_{\textrm{adv}}(0,\bm{p})$
$\displaystyle\times\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}\bigg{\\{}G_{\textrm{adv}}(\omega_{\bm{k}},\bm{p})\left.\frac{\partial}{\partial\epsilon}\right|_{\epsilon=0}\\{G_{\textrm{ret}}(\epsilon,\bm{p})G_{\textrm{ret}}(\epsilon+\omega_{\bm{k}},\bm{p})\\}+(\omega_{\bm{k}}\longleftrightarrow-\omega_{\bm{k}})\bigg{\\}}$
$\displaystyle\simeq\frac{ig^{2}}{\pi}\frac{e\hbar^{6}}{m^{2}M_{ph}}\frac{k_{F}^{4}}{(2\pi)^{6}}4(\hbar/2\tau)\int
d\Omega_{p}\int
d\Omega_{k}W^{\alpha}_{ij}({\hat{\bm{p}}},{\hat{\bm{k}}};{\bm{q}})$
$\displaystyle\times\int^{\infty}_{0}dk\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\int^{\infty}_{-\infty}dp\frac{1}{\big{[}(\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}\big{]}^{2}}\frac{1}{(\omega_{k}-\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}}$
$\displaystyle=\frac{ig^{2}}{\pi}\frac{e\hbar^{6}}{m^{2}M_{ph}}\frac{k_{F}^{4}}{(2\pi)^{6}}\frac{3\pi}{2v_{F}(\hbar/2\tau)^{4}}\int
d\Omega_{p}\int
d\Omega_{k}W^{\alpha}_{ij}({\hat{\bm{p}}},{\hat{\bm{k}}};{\bm{q}})\int^{\infty}_{0}dk\frac{k^{2}}{2\omega_{\bm{k}}\sinh{(\beta\omega_{k})}}\frac{1+(\tau\omega_{k}/\hbar)^{2}/3}{(1+(\tau\omega_{k}/\hbar)^{2})^{2}},$
(107)

where the integral over $p$ is done as

$\int^{\infty}_{-\infty}dp\frac{1}{\big{[}(\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}\big{]}^{2}\big{[}(\omega_{k}-\hbar
v_{F}p)^{2}+(\hbar/2\tau)^{2}\big{]}}=\frac{3\pi}{8v_{F}(\hbar/2\tau)^{4}}\frac{1+(\tau\omega_{k}/\hbar)^{2}/3}{(1+(\tau\omega_{k}/\hbar)^{2})^{2}}.$

Only for the limit $\tau\omega_{k}/\hbar\to 0$, we can write the resulting
dynamic spin conductivity up to the first order in frequency as

$\sigma^{\alpha}_{ij}(\omega,{\bm{q}})=\left(1+i\frac{3}{2}\omega\tau\right)\sigma^{\alpha}_{ij}({\bm{q}}).$
(108)

##  Appendix F Anisotropic cases

###  F.1 Easy-plane case

In this section, we extend the spin conductivity from the isotropic case in
Eq. (82) to anisotropic cases. Our initial focus is on the easy-plane case,
wherein the unit polarization vector $\bm{e}_{\lambda}$ of TO phonons is
constrained in the $xy$-plane, i.e., $\bm{e}_{\lambda}=\\{\hat{x},\hat{y}\\}$.
The derivation closely follows the approach outlined in the isotropic case
described in Sec. B, except that
$\tilde{\gamma}^{\alpha}(\bm{p},\bm{k};\bm{q})$ in Eq. (67) should be modified
as:

$\displaystyle\tilde{\gamma}^{\alpha,\textrm{EP}}(\bm{p},\bm{k};\bm{q})=$
$\displaystyle\sum_{a,d\in
x,y}\mathcal{P}_{ad}(\bm{k})\textrm{tr}\Big{[}\sigma_{\alpha}\hat{\bm{e}}_{a}\cdot(\bm{\sigma}\times(\bm{p}+\bm{q}+\bm{k}/2))\hat{\bm{e}}_{d}\cdot(\bm{\sigma}\times(\bm{p}+\bm{k}/2))\Big{]}$
$\displaystyle=$ $\displaystyle\;2i\sum_{a,d\in
x,y}\sum_{c,f}\Big{[}\mathcal{P}_{ad}(\bm{k})(p+q+k/2)_{c}(p+k/2)_{f}(\delta_{a\alpha}\varepsilon_{dcf}-\delta_{c\alpha}\varepsilon_{daf})\Big{]}$
$\displaystyle=$ $\displaystyle\;2i(1-\delta_{\alpha
z})\hat{\bm{e}}_{\alpha}\cdot\Big{\\{}[\bm{q}\times(\bm{p}+\bm{k}/2)]-\hat{\bm{k}}\big{(}\hat{\bm{k}}\cdot[\bm{q}\times(\bm{p}+\bm{k}/2)]\big{)}\Big{\\}}.$
(109)

Here, the factor $(1-\delta_{\alpha z})$ indicates the constraint
$\alpha=\\{x,y\\}$. If we adhere to the isotropic 3D dispersion for electrons
and phonons, i.e.,
$\xi_{\bm{p}}=\frac{\hbar^{2}(p_{x}^{2}+p_{y}^{2}+p_{z}^{2})}{2m}-\mu$ and
$\omega_{\bm{k}}=\sqrt{\hbar^{2}c^{2}(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})+E_{g}^{2}}$,
the remaining calculations are identical with those of the isotropic case, as
presented in Secs. C and D. As a result, we obtain the spin conductivity in
the easy-plane case, denoted as $\sigma_{ij}^{\alpha,\textrm{EP}}(\bm{q})$, up
to the quadratic order as follows:

$\displaystyle\sigma_{ij}^{\alpha,\textrm{EP}}(\bm{q})=(1-\delta_{\alpha
z})\sigma_{ij}^{\alpha}(\bm{q}),$ (110)

where $\sigma_{ij}^{\alpha}(\bm{q})$ in the right-hand side represents the
spin conductivity formulae in the isotropic case, as presented in Eqs. (III)
or (17).

Alternatively, in a (quasi-) 2D scenario where both electron and phonon
motions are constrained to the $xy$-plane, the electron and phonon dispersions
undergo modification, given by
$\omega_{\bm{k}}=\sqrt{\hbar^{2}c^{2}(k_{x}^{2}+k_{y}^{2})+E_{g}^{2}}$.
Consequently, the integral presented in Eq. (D) must be adapted into a 2D
integral form as follows:

$\displaystyle\sigma^{\alpha,\textrm{EP}}_{ij}(\bm{q})=$
$\displaystyle\;\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{k_{F}^{3}q^{2}}{(2\pi)^{2}}[\zeta_{ij}^{\alpha,}(\bm{q})-\zeta_{ji}^{\alpha}(\bm{q})]\int^{\infty}_{0}dp|G_{\textrm{ret}}(0,p)|^{4}$
$\displaystyle+\frac{i}{\pi}\frac{g^{2}e\hbar^{6}}{m^{2}M_{ph}}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{n_{B}(\omega_{\bm{k}})+n_{F}(\omega_{\bm{k}})}{2\omega_{\bm{k}}}(1-\hat{k}_{\alpha}^{2})\frac{4\hbar^{2}k_{F}^{5}q^{2}}{(2\pi)^{2}m}\eta_{ij}^{\alpha}(\bm{q})\int^{\infty}_{0}dp|G_{\textrm{ret}}(0,p)|^{4}G_{\textrm{ret}}(0,p).$
(111)

Here, $\zeta_{ij}^{\alpha}(\bm{q})$ and $\eta_{ij}^{\alpha}(\bm{q})$ are
defined as

$\displaystyle\zeta_{ij}^{\alpha}(\bm{q})=$
$\displaystyle\int^{2\pi}_{0}d\phi_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}=\pi(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{]},$
(112) $\displaystyle\eta_{ij}^{\alpha}(\bm{q})=$
$\displaystyle\int^{2\pi}_{0}d\phi_{p}(\hat{\bm{e}}_{i}\cdot\hat{\bm{p}})(\hat{\bm{e}}_{j}\cdot\hat{\bm{p}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{p}}\times\hat{\bm{e}}_{\alpha})\big{]}(\hat{\bm{p}}\cdot\hat{\bm{q}})=\frac{\pi}{4}\Big{[}\big{\\{}(\hat{\bm{e}}_{i}\cdot\hat{\bm{q}})\big{[}\hat{\bm{q}}\cdot(\hat{\bm{e}}_{j}\times\hat{\bm{e}}_{\alpha})\big{]}\big{\\}}+\\{i\leftrightarrow
j\\}\Big{]},$ (113)

To obtain these expressions, we performed the following approximation as in
the 3D isotropic case: (i) The electron propagator is approximated as
$G_{\textrm{ret}}(0,p)\approx(-\hbar v_{F}p+i\hbar/2\tau)^{-1}$. (ii) We
replace $\bm{p}$ with $k_{F}\hat{\bm{p}}$ in other factors involving $\bm{p}$.
(iii) The integral over $\bm{p}$ can be approximated as follows: $\int
d^{2}p\approx\int^{\infty}_{-\infty}dpk_{F}\int^{2\pi}_{0}d\phi_{p}$. (iv) We
represent the integral over $\bm{k}$ as $\int
d^{2}k=\int^{\infty}_{0}dkk\int^{2\pi}_{0}d\phi_{k}$. The integral over
$\phi_{k}$ is computed as:

$\displaystyle\int^{2\pi}_{0}d\phi_{k}(1-\hat{k}_{\alpha}^{2})=\pi.$ (114)

The remaining integrals for $p$ are provided in Eqs. (94) and (95).
Substituting all the derived integral expressions, from Eq. (112) to (114) and
Eqs. (94) and (95), into Eq. (F.1), we obtain the following expression for
$\sigma^{\alpha,\textrm{EP}}_{ij}(\bm{q})$:

$\displaystyle\sigma^{\alpha,\textrm{EP}}_{ij}(\bm{q})=$
$\displaystyle\;(1-\delta_{\alpha