††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

#  Spin-momentum Locking and Topological Vector Charge Response with Conserved
Spin

Yoonseok Hwang  Penghao Zhu  Taylor L. Hughes  Department of Physics and
Institute for Condensed Matter Theory, University of Illinois at Urbana-
Champaign, Illinois 61801, USA

###### Abstract

Spin-momentum locking plays a fundamental role in spintronics and, more
broadly, is an important concept in condensed matter physics. In 2D and 3D,
spin-momentum locking typically does not allow spin-conservation because the
spin-1/2 operators of electrons anticommute. Instead, here we study spin-
momentum locking terms with conserved, commuting pseudospins built from a
combination of spin and orbitals. We find that 2D spin-momentum locking terms
with conserved pseudospins generally lead to linearly dispersing modes at low-
energy with anomalous charge and pseudospin currents. To cure the anomaly we
show that such anomalous modes can be realized on the surface of a 3D Weyl
semimetal (or an associated weak topological insulator) with a nonzero mixed
spin-momentum quadrupole moment, which is determined by the momentum location
and pseudospin eigenvalues of Weyl points at the Fermi level. Crucially, this
mixed quadrupole moment captures a mixed pseudospin-charge bulk response that
cancels the anomaly of surface modes, and can generate a giant 3D spin Hall
effect, among other phenomena.

Spintronics, and in general the electrical control of spin, is built upon the
relativistic phenomenon of spin-orbit coupling (SOC) Dresselhaus (1955);
Bychkov and Rashba (1984a, b); Winkler (2003); Manchon  _et al._ (2015);
Soumyanarayanan  _et al._ (2016). In solids, one remarkable, low-energy
manifestation of atomic SOC is the concept of spin-momentum locking, where the
orientation of an electron’s spin is inherently tied to its direction of
motion. Spin-momentum locking is the foundation for many of the most exciting
spintronic phenomena Žutić  _et al._ (2004) in two dimensional (2D) electron
systems with Rashba SOC Bychkov and Rashba (1984b, a), the 2D surface states
of topological insulators (TIs) Liu  _et al._ (2010); Qi and Zhang (2011), the
1D helical edge states of quantum spin Hall insulators Kane and Mele (2005);
Bernevig and Zhang (2006); Wu  _et al._ (2006); Xu and Moore (2006); Roushan
_et al._ (2009).

Previous studies of spin-momentum locking have focused on couplings between
momenta and the intrinsic spin-1/2 of the electrons. Because of the
anticommuting nature of the spin-1/2 Pauli matrices, spin is not conserved in
2D and 3D systems exhibiting spin-momentum locking as more than one Pauli
matrix generally appears in the low-energy Bloch Hamiltonian. This is not
ideal for spintronics applications as the nonconserved spins have a finite
lifetime and any stored information will decay as the spins undergo
unconstrained equilibration. With this problem in mind, it is important to ask
if we can have both spin-momentum locking and conserved spins.

Below we will find that the answer to this question is yes, and furthermore,
that in some cases these systems are anomalous in the sense that they have a
mixed spin-charge ‘t-Hooft anomaly, which can be remedied when attached to the
boundary of a higher dimensional bulk. In addition, we find that in some
crystal symmetry classes, the higher dimensional bulk response that
compensates for the surface anomaly is captured by an inherently non-
relativistic, 3D rank-2 symmetric or antisymmetric axial vector charge Chern-
Simons term, where the axial vector charge is the conserved spin.
Interestingly, we show that the bulk response is naturally generated by a
topological Weyl semimetal with a nodal configuration that generates a mixed
spin-momentum Weyl quadrupole. This is in analogy to the recently proposed
rank-2 chiral fermions which carry a crystal momentum vector charge and a
momentum-momentum Weyl quadrupole Dubinkin  _et al._ (2021); Zhu  _et al._
(2023).

Explicitly, here we study 2D spin-momentum locking terms in nonmagnetic
crystalline systems with rotation and mirror symmetries. We focus on effective
$k\cdot p$ Hamiltonians having the form

$\displaystyle h_{\mathbf{k}}=\sum_{i=x,y}\sum_{A=X,Y}g_{iA}k_{i}S_{A},$ (1)

where $(S_{X},S_{Y})$ are spin operators satisfying $S_{A}^{\dagger}=S_{A}$,
and ${\rm Tr}[S_{A}]=0$, and $g_{iA}$ is $2\times 2$ matrix with real entries
that will be constrained by the imposed symmetries. Note that we work in units
$\hbar=c=1$ and the electron charge is $-e$. We expect that effective
Hamiltonians of this form may appear when the full Bloch Hamiltonian is
expanded around time-reversal invariant momenta (TRIM) in the Brillouin zone
(BZ). In the neighborhood of a TRIM, time-reversal symmetry simply flips both
spin and momentum in Eq. (1), and hence the effective Hamiltonian is
manifestly time-reversal invariant.

The simplest case of spin-momentum locking with conserved spin in Eq. (1) is
when only one spin operator appears. This happens in the 1D helical edge
states of the quantum spin Hall insulator, but in 2D, this case is allowed
only in lower-symmetry systems having 2D point groups $C_{1,2}$ or $C_{1v,2v}$
(and even then would require extreme anisotropy or fine-tuning) Winkler
(2003); Samokhin (2022). We note that for $C_{3}$ and $C_{3v}$ symmetries,
there can also be spin-momentum locking terms with one conserved spin operator
$S_{Z}$, which are cubic in momentum Fu (2009); Samokhin (2022), but we leave
further analysis of higher-order terms as a future work.

Our primary focus will instead be on systems with $C_{3,4,6}$ or
$C_{3v,4v,6v}$, where, as shown in the Supplementary Material (SM) Sec. S1 sup
, the spin-momentum locking terms that are _linear_ in momentum always involve
both $S_{X}$ and $S_{Y}$. In these cases the spin is conserved only when
$S_{X}$ and $S_{Y}$ commute with each other, i.e., $[S_{X},S_{Y}]=0$. Although
the intrinsic spin-1/2 of the electrons cannot satisfy this condition on their
own, we will show that commuting operators consisting of both orbital- and
spin-degrees of freedom of electrons can be constructed. Hence, we study
generalized spin-momentum locking terms in Eq. (1) with commuting pseudospin
operators $S_{X}$ and $S_{Y}$ that are both conserved.

Rashba term with commuting pseudospins.— For our primary example we focus on
the point group $C_{4v}$, which is promoted to the $P4mm$ space group for
lattice Hamiltonians. We further limit our study to nonmagnetic systems having
time reversal symmetry (TRS) $\mathcal{T}$ with $\mathcal{T}^{2}=-1$. The
relevant symmetry elements in $C_{4v}$ act on the position-space coordinates
as $M_{x}:(x,y,z)\to(-x,y,z)$, $M_{y}:(x,y,z)\to(x,-y,z)$, and
$C_{4z}:(x,y,z)\to(-y,x,z)$. We choose the following representation of the
symmetry elements:

$\displaystyle
M_{x,y}=-i\tau_{0}\sigma_{x,y},\,C_{4z}=\tau_{z}e^{-i\frac{\pi}{4}\sigma_{z}},\,\mathcal{T}=-i\tau_{0}\sigma_{y}\mathcal{K},$
(2)

where $\mathcal{K}$ is complex conjugation, $\sigma_{a}$ are the Pauli
matrices for the spin-1/2 degree of freedom, and the $\tau_{a}$ Pauli matrices
represent a pair of orbitals transforming under the trivial (e.g.,
$s$-orbital) and $B_{1}$ (e.g., $d_{x^{2}-y^{2}}$-orbital) representations of
$C_{4v}$. From these degrees of freedom we can construct commuting pseudospin
operators $S_{X}$ and $S_{Y}$:

$\displaystyle S_{X}=\frac{1}{\sqrt{2}}(\tau_{z}+\tau_{x})\sigma_{x},\quad
S_{Y}=\frac{1}{\sqrt{2}}(\tau_{z}-\tau_{x})\sigma_{y}.$ (3)

We see that $S_{X}$ and $S_{Y}$ commute with each other and transform under
$C_{4v}$ and $\mathcal{T}$ as a axial vector: $M_{x}$ ($M_{y}$) flips $S_{Y}$
($S_{X}$), $\mathcal{T}$ flips both $S_{X}$ and $S_{Y}$, and $C_{4z}$ rotates
$(S_{X},S_{Y})$ to $(S_{Y},-S_{X})$.

Using this construction we can now consider spin-momentum locking phenomena
with commuting spins. Independent of the spin commutation relations, $C_{4v}$
and $\mathcal{T}$ symmetries restrict the coefficients $g_{iA}$ in Eq. 1 to
$g_{xY}=-g_{yX}=v$ and $g_{xX}=g_{yY}=0$ (see the SM Sec. S1 sup  for
details). Without loss of generality, we set $v=1$ in the following unless
otherwise specified. Thus, Eq. (1) is reduced to

$\displaystyle h_{{\mathbf{k}}}=k_{x}S_{Y}-k_{y}S_{X},$ (4)

which is an analog of the Rashba term Bychkov and Rashba (1984b) but with
commuting, instead of anticommuting, spin operators. For simplicity, we call
the Hamiltonian described by Eq. (4) a commuting-spin Rashba (CSR) term. The
spectrum of the CSR term has four bands $E_{{\mathbf{k}}}=\pm k_{x}\mp k_{y}$
and they, and the Fermi-surface contours of $h_{{\mathbf{k}}}$, are shown in
Figs. 1(a) and (b). We see that $h_{{\mathbf{k}}}$ has independent symmetries
$U(1)_{A}$ generated by $S_{A}$ ($A=X,Y$), and each eigenstate of
$h_{{\mathbf{k}}}$ can be labeled by the quantum numbers $(s_{X},s_{Y})$.
Hence, each of the four open Fermi-surfaces has a fixed spin polarization that
determines the Fermi-velocity $v_{{\mathbf{k}}}=(s_{Y},-s_{X})$. In
comparison, instead of having four separate Fermi-surfaces that transform into
each other under $C_{4},$ the conventional Rashba Fermi-surface is a single
closed circle where the spin-polarization and group velocity continuously wind
around the Fermi-surface.

Figure 1:  Energy spectrum and anomalous responses of $h_{{\mathbf{k}}}$. (a)
A portion of the energy spectrum for eigenstates of $h_{{\mathbf{k}}}$ with
different pseudospins $(s_{X},s_{Y})$. (b) Fermi-surface contours of
$h_{{\mathbf{k}}}$. (c)-(d) Low-energy spectra and fillings as a function of
$k_{x}$ at a fixed $k_{y}>0$ when (c) the gauge fields are zero, (d) the
charge $U(1)_{\rm EM}$ gauge field is adiabatically shifted by
$2\pi/(eL_{x})$, and (e) the spin $U(1)_{Y}$ gauge field is adiabatically
shifted by $2\pi/L_{x}$. Filled (empty) circles represent occupied
(unoccupied) states. Throughout this figure, pseudospin sectors $(-,+)$,
$(+,+)$, $(+,-)$, and $(-,-)$ are color-coded by blue, yellow, red, and gray,
respectively.

Anomalous charge and pseudospin currents.— We will now show that the CSR Fermi
surfaces have a mixed anomaly between the electric charge $U(1)_{\rm EM}$ and
the spin $U(1)_{A}$ symmetries. To demonstrate this, we view the CSR mode as a
family of 1D chiral modes propagating in the $x$-direction ($y$-direction)
with chirality $\chi_{x}(s_{Y})={\rm sgn}(s_{Y})$ [$\chi_{y}(s_{X})=-{\rm
sgn}(s_{X})$] parameterized by each $k_{y}$ ($k_{x}$). For example, at a fixed
value of $k_{y}$ there are four linearly-dispersing modes propagating in the
$x$-direction with chirality ${\rm sgn}(s_{Y}),$ each with an effective
chemical potential proportional to the fixed $k_{y}.$ Hence, the charge and
pseudospin anomalies of the CSR modes can be straightforwardly determined from
the anomalies of 1D chiral modes.

Hence, to determine the anomalies of the CSR modes, let us first recall the
anomaly of a 1D chiral mode with dispersion $E_{k_{x}}=\chi_{x}k_{x}.$ Let us
turn on an electric field by adiabatically shifting the vector potential,
$A_{x}$, by $2\pi/(eL_{x})$ over a time period $T$. Explicitly, we have
$A_{x}(t)=2\pi t/(eL_{x}T)$ and electric field $E_{x}=-2\pi/(eL_{x}T)$ for
$t\in[0,T]$. As $t$ goes from $0$ to $T$, the number of occupied right-moving
modes ($\chi_{x}=+1$) increases by one, and the number of occupied left-moving
modes ($\chi_{x}=-1$) decreases by one. To interpret this as an anomalous
conservation law we define the charge density $J^{0}=-\frac{e}{V}\sum_{n\in
occ.}1$ where $V$ is the volume (area in 2D, length in 1D) of the system.
Thus, during this process $\Delta J^{0}=-e\chi_{x}/V,$ and the corresponding
anomalous conservation law is

$\displaystyle\partial_{\mu}J^{\mu}=\chi_{x}\frac{e^{2}}{2\pi}E_{x}.$ (5)

Equipped with Eq. (5), we first consider the anomalies of the CSR mode induced
by the $U(1)_{\rm EM}$ gauge field $A_{i}$ ($i=x,y$). The gauge field $A_{i}$
couples with the CSR mode as:

$\displaystyle h_{\mathbf{k}}[A]=(k_{x}+eA_{x})S_{Y}-(k_{y}+eA_{y})S_{X}.$ (6)

Considering the CSR mode as a family of 1D chiral modes, we see that the net
chiralities along both the $x$\- and $y$-directions are zero because chiral
modes with opposite pseudospin eigenvalues have opposite chiralities as
required by $\mathcal{T}$. Therefore, the electric field will not induce an
anomaly for the electric charge. However, we can define the pseudospin density
$\mathcal{J}_{A}^{0}=\frac{1}{V}\sum_{n\in{\rm occ}}s_{n,A}$ where $A=X,Y$ and
$s_{n,A}$ is the pseudospin eigenvalue of the $n$-th eigenstate.
Interestingly, we find that the pseudospin density changes by a nonzero amount
during an adiabatic shift of $A_{i}$, because the modes with opposite
chiralities also have opposite pseudospin eigenvalues. As illustrated in Figs.
1(c) and (d), for a shift $\Delta A_{x}=2\pi/(eL_{x})$ we find
$\Delta\mathcal{J}_{Y}^{0}=\frac{1}{V}\sum_{|k_{y}|<\Lambda_{y}}4\,=\frac{4\Lambda_{y}}{\pi
L_{x}},$ and similarly, for $\Delta A_{y}=2\pi/(eL_{y})$ we find
$\Delta\mathcal{J}_{X}^{0}=\frac{1}{L_{x}L_{y}}\sum_{|k_{x}|<\Lambda_{x}}-4=-\frac{4\Lambda_{x}}{\pi
L_{y}}$. Here, $L_{i}$ and $\Lambda_{i}$ are the linear system sizes and
momentum cutoffs along the $i=x,y$ directions. These responses can be
expressed by the anomalous conservation laws:

$\displaystyle\partial_{\mu}\mathcal{J}^{\mu}_{X}=\frac{2e\Lambda_{x}}{\pi^{2}}E_{y},\quad\partial_{\mu}\mathcal{J}^{\mu}_{Y}=-\frac{2e\Lambda_{y}}{\pi^{2}}E_{x},$
(7)

with the electric field $E_{i}=-\partial_{t}A_{i}+\partial_{i}A_{0}$
($i=x,y$).111We use the mostly-plus metric such that $A^{0}=-A_{0}$ and
$A^{i}=A_{i}$. Furthermore, because of the $C_{4}$ symmetry it is natural to
choose $\Lambda_{x}=\Lambda_{y}=\Lambda.$

As befitting a mixed anomaly, we now discuss the inverse effect where a
pseudospin electric field leads to an electric charge anomaly. We denote the
pseudospin $U(1)_{A}$ gauge fields as $a^{A}_{i}$ ($A=X,Y$),222One can
understand $a^{A}_{i}$ as the pseudospin connection. and couple them to
$h_{\mathbf{k}}$ via:

$\displaystyle
h_{\mathbf{k}}[a]=(k_{x}-S_{A}a^{A}_{x})S_{Y}-(k_{y}-S_{A}a^{A}_{y})S_{X}.$
(8)

Note that $a^{A}_{i}$ always appears as $S_{A}a^{A}_{i}$, which indicates that
the modes with opposite pseudospins feel opposite _effective_ electric fields.
If we adiabatically shift $a_{x}^{X}$ by $2\pi/L_{x}$, then each $k_{y}$
contributes, $\Delta J^{0}=-e$ $[+e]$ in sectors with pseudospin eigenvalues
$(+1,-1)$ and $(-1,+1)$ [$(+1,+1)$ and $(-1,-1)$]. Therefore, the shift of
$a_{x}^{X}$ does not lead to charge anomaly, and similarly for $a_{y}^{Y}$.
However, as illustrated in Figs. 1 (c) and (e), shifting $a^{Y}_{x}$ or
$a^{X}_{y}$, generates nonzero charge density. These responses are encoded in
the anomalous charge conservation law:

$\displaystyle\partial_{\mu}J^{\mu}=\frac{2e\Lambda}{\pi^{2}}(\mathcal{E}^{X}_{y}-\mathcal{E}^{Y}_{x}),$
(9)

where the spin electric field is
$\mathcal{E}^{A}_{i}=-\partial_{t}a^{A}_{i}+\partial_{i}a^{A}_{0}$.

While Eq. (9) describes a charge anomaly in response to an antisymmetric
combination of pseudospin electric fields, we can also find cases generated by
a symmetric combination. For example, if we choose $C_{4z}$ rotation such that
$C_{4z}:(x,y,z)\to(-y,x,z)$ and $(S_{X},S_{Y})\to(-S_{Y},S_{X})$, such that
the pseudospin rotates opposite to our previous convention, then we find
$g_{xY}=g_{yX}=v$ and $g_{xX}=g_{yY}=0$. This leads to a CSR term
$h_{{\mathbf{k}}}=v(k_{x}S_{Y}+k_{y}S_{X})$, which has its charge anomaly
induced by the symmetric combination of pseudospin electric fields, i.e.,
$\mathcal{E}^{X}_{y}+\mathcal{E}^{Y}_{x}$. Similar symmetric vector charge
anomalies have been discussed in the recently proposed rank-2 chiral modes
with dispersion $E=k_{x}k_{y}$, where the charge anomaly is induced by a
symmetric strain field, and the vector charge is momentum instead of
pseudospin Dubinkin  _et al._ (2021); Zhu  _et al._ (2023).

Realization of the CSR on the surface of a 3D semimetal.— The anomalies
discussed above imply that the CSR modes cannot exist in isolation in a 2D
lattice system. There are two natural ways to treat this issue: (i) introduce
symmetry preserving terms that generate closed, non-anomalous Fermi surfaces
intrinsically in 2D, or (ii) realize the CSR on the surface of a 3D system
such that the bulk cancels the anomalous current via the Callan-Harvey anomaly
inflow Callan Jr and Harvey (1985). For (i) we can consider a Hamiltonian

$\displaystyle h_{{\mathbf{k}}}=\alpha(k_{x}^{2}+k_{y}^{2})\mathds{1}_{4\times
4}+vk_{x}S_{Y}-vk_{y}S_{X},$ (10)

where we restored $v$. The $\alpha$-dependent term preserves all of the
symmetries and modifies the energy bands as
$E_{{\mathbf{k}}}=\alpha(k_{x}^{2}+k_{y}^{2})\pm vk_{x}\pm vk_{y}.$ These
bands generate four closed, circular Fermi surfaces centered around $(\pm
v/2\alpha,\pm v/2\alpha).$ Furthermore these Fermi-surfaces will not be
anomalous if
$|\Lambda|>|\tfrac{v}{2\alpha}|+\sqrt{\tfrac{\mu}{\alpha}+\tfrac{v^{2}}{2\alpha^{2}}}$,
i.e., if the region inside the momentum cutoff includes all four closed
circles generated at the chemical potential $\mu$. This is what we would
expect for a 2D lattice model Bloch Hamiltonian that realizes Eq. (4) when
expanded around a TRIM. Note that the $\alpha$-term is not the only possible
symmetry preserving perturbation, and discussions of another allowed term can
be found in SM Sec. S2 sup .

For a more interesting resolution to the anomaly let us consider mechanism
(ii). To realize the anomalous CSR mode as a boundary mode we consider a 3D
tight-binding (TB) model,

$\displaystyle H_{3D}({\mathbf{k}})$ $\displaystyle=\sin
k_{x}\mu_{x}S_{Y}-\sin k_{y}\mu_{x}S_{X}+\sin k_{z}\mu_{y}$ (11)
$\displaystyle+(m-1+\cos k_{x}+\cos k_{y}+\cos k_{z})\mu_{z},$

which is an eight-band model expressed using three sets of Pauli matrices
$\mu$, $\tau$, and $\sigma$; the latter two have the same interpretation as in
Eq. (3). This model has space group $P4mm$ (No. 99) whose symmetry generators
are $C_{4v}+\mathcal{T}$ and the 3D translations. In addition,
$H_{3D}({\mathbf{k}})$ has two conserved spins $S_{X,Y}$ such that
$[H_{3D}({\mathbf{k}}),S_{X,Y}]=0$.

Figure 2:  Tight-binding model $H_{3D}({\mathbf{k}})$ for 3D semimetal with
nonzero $Q_{iA}$. (a) Energy spectra in the $k_{z}=0$ plane for $m=-1$. Eight
Weyl points are located in pairs at $(\pm\pi/3,\pm\pi/3,0)$. [See also Fig.
3(a).] The representative point $(\pi/3,\pi/3,0)$ is indicated by the filled
black dots. (b) [001] surface energy spectrum for a 20-layer slab. The CSR
modes (red color) are realized as the surface states of 3D semimetal, [c.f.
Fig. 1(a)]. (c) Spin density induced by magnetic flux $\Phi=\pi/2$. Note that
$2\pi$ flux is equivalent to a magnetic flux quantum. (d) Charge density
induced by spin magnetic flux $\varphi_{X}=\pi/2$. Note that
$\widetilde{\mathcal{J}}^{0}_{X}=\frac{4\pi^{2}}{Q_{0}}\mathcal{J}^{0}_{X}$
and $\widetilde{J}^{0}=\frac{4\pi^{2}}{Q_{0}}J^{0}$. A pair of (spin) magnetic
flux tubes with opposite strength is inserted along the periodic
$x$-direction. The size of the lattice along the $y$\- and $z$-directions is
$(N_{y},N_{z})=(32,32)$.

Let us focus on the regime $-2<m<0$, which is a semimetal where 8 Weyl points
are located at ${\mathbf{K}}=(\pm K_{*},\pm K_{*},0)$ with
$K_{*}=\cos^{-1}(-m/2)$ [see Fig. 2(a) and Fig. 3(a)]. At each ${\mathbf{K}}$,
there are two Weyl points with opposite chiralities, and they cannot be gapped
because they are in different pseudospin sectors. For example, the two Weyl
points with chirality $\chi=1$ and $-1$ at $(K_{*},K_{*},0)$ have pseudospin
eigenvalues $(1,1)$ and $(-1,-1)$, respectively [see Fig. 2(a)].

If we have an open boundary along the $z$-direction, we can analytically
derive the Hamiltonian $H_{\rm surf}(k_{x},k_{y})$ of the low-energy surface
localized modes following Ref. Dwivedi and Ramamurthy, 2016. We find that
$H_{\rm surf}(k_{x},k_{y})=k_{x}S_{Y}-k_{y}S_{X}$, which is exactly the CSR
term. As shown in Fig. 2(b), the surface CSR modes exist only in a restricted
region of the surface BZ, $-K_{*}\leq k_{x,y}\leq K_{*}$, defined by the
location of the Weyl points. Thus, the cutoff $\Lambda$ in the anomalous
conservation laws for the CSR modes is now naturally identified as
$\Lambda=K_{*}$.

Bulk response and Weyl quadrupole.— Because of the anomalous surface response
we expect the Weyl semimetal to produce a compensating bulk response. We
calculated the quasi-topological bulk charge and pseudo-spin responses to
external fields and found (detailed Kubo formula derivations can be found in
SM Sec. S3 sup ):

$\displaystyle J^{\mu}_{3D}=-\frac{e}{4\pi^{2}}\epsilon^{\mu
i\nu\rho}Q_{iA}\partial_{\nu}a^{A}_{\rho},$
$\displaystyle\mathcal{J}^{\mu}_{A,3D}=-\frac{e}{4\pi^{2}}\epsilon^{\mu
i\nu\rho}Q_{iA}\partial_{\nu}A_{\rho},$ (12)

where the response coefficient $Q_{iA}$ is the Weyl spin-momentum quadrupole
moment defined as

$\displaystyle Q_{iA}=\sum_{\alpha}\chi_{\alpha}K_{\alpha,i}s_{\alpha,A},$
(13)

where the $\alpha$ in Eq. (13) is summed over all Weyl points sitting on the
bulk FS. The matrix $Q_{iA}$ is only well-defined when the Weyl momentum and
spin dipole moments vanish. Because of time-reversal the former vanishes,
while a spin-resolved version of the Nielsen-Ninomiya no-go theorem Nielsen
and Ninomiya (1981) requires the latter to vanish. For $C_{nv}$ with $n>2$,
$Q_{iA}$ is proportional to the identity matrix, i.e.,
$Q_{iA}=Q_{0}\delta_{iA}$ (see the SM Sec. S3). Indeed, direct calculations of
$Q_{iA}$ for our model Eq. (11) show that $Q_{xY}=Q_{yX}=0$ and
$Q_{xX}=Q_{yY}=8K_{*}\equiv Q_{0}$. These non-zero coefficients yield a mixed
Chern-Simons term effective action

$\displaystyle S_{CS}$ $\displaystyle=\frac{e}{4\pi^{2}}\int
d^{4}x\epsilon^{\mu\nu\rho\sigma}Q_{\mu
A}a^{A}_{\nu}\partial_{\rho}A_{\sigma}$ (14) $\displaystyle=\frac{Q_{0}e}{2\pi
h}\int
d^{4}x[A_{z}(\mathcal{E}^{Y}_{x}-\mathcal{E}^{X}_{y})-A_{x}\mathcal{E}^{Y}_{z}$
$\displaystyle+A_{y}\mathcal{E}^{X}_{z}+A_{0}(\mathcal{B}_{x}^{X}+\mathcal{B}_{y}^{Y})],$
(15)

where $\mathcal{B}^{A}_{i}=\epsilon_{ijk}\partial_{j}a^{A}_{k}$ is the
pseudospin magnetic field.

We can now study some of the response phenomena from Eqs. (12) or (14). First,
we have a 3D pseudospin Hall effect given by
$\mathcal{J}^{i}_{A}=-\frac{e}{4\pi^{2}}Q_{kA}\epsilon^{ijk}E_{j}=\sigma_{S}\epsilon^{Aij}E_{j}$
where the spin Hall conductivity $\sigma_{S}=-\frac{e}{4\pi^{2}}Q_{0}$. For a
general symmetry class, the pseudospin Hall conductivity tensor
$\sigma_{A}^{ij}=\frac{e}{4\pi^{2}}Q_{kA}\epsilon^{ijk}$. Additionally, Eq.
(14) implies that magnetic field can induce pseudospin density in the bulk

$\displaystyle\mathcal{J}^{0}_{X,3D}=-\frac{e}{4\pi^{2}}Q_{0}B_{x},\
\mathcal{J}^{0}_{Y,3D}=-\frac{e}{4\pi^{2}}Q_{0}B_{y},$ (16)

and, inversely, that the pseudospin magnetic flux can bind electric charge
density

$\displaystyle
J^{0}_{3D}=-\frac{e}{4\pi^{2}}Q_{0}(\mathcal{B}^{X}_{x}+\mathcal{B}^{Y}_{y}).$
(17)

We confirm these bulk predictions with numerical calculations, and in Figs.
2(c) and (d) we show the results for magnetic field (pseudospin magnetic
field) induced pseudospin (charge) density, which match Eqs. (16) and (17)
(see the SM Sec. S5 sup  for details).

The bulk responses in Eq. (12) compensate the anomalous CSR modes that appear
on the surface of $H_{3D}({\mathbf{k}})$. Without a boundary, it is
straightforward to see $\partial_{\mu}J^{\mu}_{3D}=0$. Now consider a boundary
at $z=0$ such that there is a non-vanishing (vanishing) response in the system
(vacuum) for $z<0$ ($z>0$). Then we can model $Q_{iA}(z)=Q_{iA}\Theta(-z)$ as
a $z$-dependent step function, and hence there will be anomaly ascribed to the
boundary:

$\displaystyle\partial_{\mu}J^{\mu}_{3D}=\partial_{z}J^{z}_{3D}=-\delta(z)\frac{Q_{0}e}{4\pi^{2}}(\mathcal{E}^{X}_{y}-\mathcal{E}^{Y}_{x}).$
(18)

By recalling $Q_{0}=8K_{*}=8\Lambda$, it is clear that $\int
dz\partial_{\mu}J^{\mu}_{3D}$ exactly cancels the $\partial_{\mu}J^{\mu}$ from
the surface CSR mode [see Eq. (9)], which is nothing but the Callan-Harvey
mechanism Callan Jr and Harvey (1985). In a similar way, the surface
pseudospin anomaly is compensated by the bulk inflow.

Let us make three further remarks about the bulk response. (i) We note that
symmetry-preserving perturbations can change the locations of Weyl points in
momentum space, and thus deform the surface Hamiltonian away from the exact
form in Eq. (4). As an illustration, we show the splitting of the Weyl point
pairs, and the deformation of the surface mode, induced by a nontrivial
symmetry preserving perturbation term, $\beta\sin k_{x}\sin
k_{y}\mu_{x}S_{X}S_{Y}$ in Fig. 3 (see SM sup  Sec. S2 for more detailed
discussion). However, this deformation does not qualitatively affect the bulk
response and the anomalies of the surface mode, which are still determined by
$Q_{iA}$. (ii) We note that if we take the alternative $C_{4z}$ mentioned
below Eq. (9), then in Eq. (15)
$\mathcal{E}^{Y}_{i}\rightarrow-\mathcal{E}^{Y}_{i}$ and
$\mathcal{B}^{Y}_{i}\rightarrow-\mathcal{B}^{Y}_{i}$ like the anomalous
conservation equation Eq. (9). With this change, the bulk response becomes a
pseudospin version of the symmetric rank-2 vector charge Chern-Simons term
discussed in Refs. Dubinkin  _et al._ , 2021; Zhu  _et al._ , 2023, which is
the bulk response term for a 3D system with surface rank-2 chiral mode. (iii)
If we annihilate the Weyl points with opposite chiralities in the model Eq.
(11), we can derive a weak TI with $Q_{xX}=Q_{yY}=8\pi/a$ where $a$ is the
lattice constant. This weak TI would exhibit a giant spin-Hall conductivity
$\sigma_{S}=2e/(a\pi)$. We could furthermore introduce the translation-
symmetry gauge field $\mathfrak{e}^{l}_{\mu}$ as discussed in Refs. Nissinen
and Volovik, 2018, 2019; Gioia  _et al._ , 2021; Wang  _et al._ , 2021, to
find the effective action

$\displaystyle S_{CS}$ $\displaystyle=\frac{e}{4\pi^{2}}\int
d^{4}x\epsilon^{\mu\nu\rho\sigma}Q_{lA}\mathfrak{e}^{l}_{\mu}a^{A}_{\nu}\partial_{\rho}A_{\sigma}.$
(19)

In addition to the responses above, this term captures responses of the weak
TI to dislocations Ran  _et al._ (2009). Interested readers are referred to
the SM Sec. S6 sup  for more details.

Figure 3:  Effect of symmetry-preserving perturbation. (a) Splitting of Weyl
points in the $k_{z}=0$ plane by the perturbation $\beta\sin k_{x}\sin
k_{y}\mu_{x}S_{X}S_{Y}$. Before adding the perturbation ($\beta=0$), two Weyl
points are degenerate at the gray dots. (b) Fermi-energy contours for
$\beta=0.8,m=-1$. (c) [001] surface energy spectrum for 20 layers when
$\beta=0.8,m=-1$. The deformed CSR surface modes on the top surface are
indicated by the red color.

Conclusions.— Here we studied $C_{4v}$-symmetric spin-momentum locking terms
with conserved spin, and discussed the dynamics and anomalies of a commuting
spin Rashba term. We showed that the CSR is anomalous by itself and can appear
on the surface of a 3D topological semimetal. Such a semimetal will have a
vector charge Chern-Simons theory bulk response with a response coefficient
determined by the Weyl spin-momentum quadrupole moment of the Fermi surface.
Beyond the specific crystalline symmetry group studied in this Letter, our
approach can be applied to general commuting spin-momentum locking terms in
Eq. (1). Interested readers can find relevant discussions in SM Sec. S4 sup .
This work opens an important new direction in the quasi-topological semimetal
responses characterized by mixed multipole moments of the nodal Fermi
surfaces. Furthermore, we expect this work to be relevant for both quantum
materials and metamaterials. Indeed, a pre-cursor to this work Dubinkin  _et
al._ (2021), where the vector charge was realized by momentum instead of spin,
has recently been explored experimentally in a topo-electric circuit material
Zhu  _et al._ (2023). Additionally it would be interesting to investigate
connections between our work and the recently proposed spin-resolved topology
Lin  _et al._ (2022).

###### Acknowledgements.

Y.H. thanks the Air Force Office of Scientific Research under award number
FA9550-21-1-0131 for support. Y.H. and P.Z. thank the US Office of Naval
Research (ONR) Multidisciplinary University Research Initiative (MURI) grant
N00014-20-1-2325 on Robust Photonic Materials with High-Order Topological
Protection. T.L.H. thanks the ARO under the MURI Grant W911NF2020166 for
support.

## References

  * Dresselhaus (1955) Gene Dresselhaus, “Spin-orbit coupling effects in zinc blende structures,” Physical Review 100, 580 (1955).
  * Bychkov and Rashba (1984a) Yua A Bychkov and E I Rashba, “Properties of a 2D electron gas with lifted spectral degeneracy,” JETP lett 39, 78 (1984a).
  * Bychkov and Rashba (1984b) Yu A Bychkov and Emmanuel I Rashba, “Oscillatory effects and the magnetic susceptibility of carriers in inversion layers,” Journal of physics C: Solid state physics 17, 6039 (1984b).
  * Winkler (2003) Roland Winkler, “Spin-orbit coupling effects in two-dimensional electron and hole systems,” (2003).
  * Manchon  _et al._ (2015) Aurelien Manchon, Hyun Cheol Koo, Junsaku Nitta, Sergey M Frolov, and Rembert A Duine, “New perspectives for Rashba spin-orbit coupling,” Nature materials 14, 871–882 (2015).
  * Soumyanarayanan  _et al._ (2016) Anjan Soumyanarayanan, Nicolas Reyren, Albert Fert, and Christos Panagopoulos, “Emergent phenomena induced by spin-orbit coupling at surfaces and interfaces,” Nature 539, 509–517 (2016).
  * Žutić  _et al._ (2004) Igor Žutić, Jaroslav Fabian, and S Das Sarma, “Spintronics: Fundamentals and applications,” Reviews of Modern Physics 76, 323 (2004).
  * Liu  _et al._ (2010) Chao-Xing Liu, Xiao-Liang Qi, HaiJun Zhang, Xi Dai, Zhong Fang, and Shou-Cheng Zhang, “Model hamiltonian for topological insulators,” Physical Review B 82, 045122 (2010).
  * Qi and Zhang (2011) Xiao-Liang Qi and Shou-Cheng Zhang, “Topological insulators and superconductors,” Review of Modern Physics 83, 1057–1110 (2011).
  * Kane and Mele (2005) Charles L Kane and Eugene J Mele, “Quantum spin Hall effect in graphene,” Physical Review Letters 95, 226801 (2005).
  * Bernevig and Zhang (2006) B Andrei Bernevig and Shou-Cheng Zhang, “Quantum spin Hall effect,” Physical Review Letters 96, 106802 (2006).
  * Wu  _et al._ (2006) Congjun Wu, B Andrei Bernevig, and Shou-Cheng Zhang, “Helical liquid and the edge of quantum spin Hall systems,” Physical Review Letters 96, 106401 (2006).
  * Xu and Moore (2006) Cenke Xu and Joel E Moore, “Stability of the quantum spin Hall effect: Effects of interactions, disorder, and $\mathbb{Z}_{2}$ topology,” Physical Review B 73, 045322 (2006).
  * Roushan  _et al._ (2009) Pedram Roushan, Jungpil Seo, Colin V Parker, Yew San Hor, David Hsieh, Dong Qian, Anthony Richardella, M Zahid Hasan, Robert Joseph Cava, and Ali Yazdani, “Topological surface states protected from backscattering by chiral spin texture,” Nature 460, 1106–1109 (2009).
  * Dubinkin  _et al._ (2021) Oleg Dubinkin, FJ Burnell, and Taylor L Hughes, “Higher rank chiral fermions in 3d Weyl semimetals,” arXiv preprint arXiv:2102.08959 (2021).
  * Zhu  _et al._ (2023) Penghao Zhu, Xiao-Qi Sun, Taylor L Hughes, and Gaurav Bahl, “Higher rank chirality and non-Hermitian skin effect in a topolectrical circuit,” Nature Communications 14, 720 (2023).
  * Samokhin (2022) KV Samokhin, “On the effective models of spin-orbit coupling in a two-dimensional electron gas,” Annals of Physics 437, 168710 (2022).
  * Fu (2009) Liang Fu, “Hexagonal warping effects in the surface states of the topological insulator Bi2Te3,” Physical Review Letters 103, 266801 (2009).
  * (19) See the Supplemental Material for more details about the symmetry analysis for two-dimensional spin-momentum locking terms, the effect of symmetry-preserving perturbation, derivation of linear responses, and the discussion on the responses to dislocation.
  * Callan Jr and Harvey (1985) Curtis G Callan Jr and Jeffrey A Harvey, “Anomalies and fermion zero modes on strings and domain walls,” Nuclear Physics B 250, 427–436 (1985).
  * Dwivedi and Ramamurthy (2016) Vatsal Dwivedi and Srinidhi T Ramamurthy, “Connecting the dots: Time-reversal symmetric Weyl semimetals with tunable Fermi arcs,” Physical Review B 94, 245143 (2016).
  * Nielsen and Ninomiya (1981) H.B. Nielsen and M. Ninomiya, “A no-go theorem for regularizing chiral fermions,” Physics Letters B 105, 219–223 (1981).
  * Nissinen and Volovik (2018) J Nissinen and GE Volovik, “Tetrads in solids: from elasticity theory to topological quantum Hall systems and Weyl fermions,” Journal of Experimental and Theoretical Physics 127, 948–957 (2018).
  * Nissinen and Volovik (2019) J Nissinen and GE Volovik, “Elasticity tetrads, mixed axial-gravitational anomalies, and (3+1)-d quantum Hall effect,” Physical Review Research 1, 023007 (2019).
  * Gioia  _et al._ (2021) L. Gioia, Chong Wang, and A. A. Burkov, “Unquantized anomalies in topological semimetals,” Phys. Rev. Res. 3, 043067 (2021).
  * Wang  _et al._ (2021) Chong Wang, Alexander Hickey, Xuzhe Ying, and A. A. Burkov, “Emergent anomalies and generalized Luttinger theorems in metals and semimetals,” Phys. Rev. B 104, 235113 (2021).
  * Ran  _et al._ (2009) Ying Ran, Yi Zhang, and Ashvin Vishwanath, “One-dimensional topologically protected modes in topological insulators with lattice dislocations,” Nature Physics 5, 298–303 (2009).
  * Lin  _et al._ (2022) Kuan-Sen Lin, Giandomenico Palumbo, Zhaopeng Guo, Jeremy Blackburn, Daniel P Shoemaker, Fahad Mahmood, Zhijun Wang, Gregory A Fiete, Benjamin J Wieder, and Barry Bradlyn, “Spin-resolved topology and partial axion angles in three-dimensional insulators,” arXiv preprint arXiv:2207.10099 (2022).

Supplemental Material for “ Spin-momentum Locking and Topological Vector
Charge Response with Conserved Spin”

###### Contents

  1. S1 $C_{n}$ and $C_{nv}$ symmetric spin-momentum locking terms in 2D
    1. S1.1 General formula for leading terms
    2. S1.2 Symmetry constraints of $g_{iA}$
  2. S2 Perturbations effect on the commuting-spin Rashba mode
  3. S3 Bulk responses
    1. S3.1 Derivation of linear responses
    2. S3.2 Symmetry constraints on $Q_{iA}$
  4. S4 Anomalies of the generalized spin-momentum locked mode
    1. S4.1 Current anomaly
    2. S4.2 Anomaly cancellation and spin-momentum quadrupole
  5. S5 Numerical scheme for linear response
  6. S6 Response to dislocations

##  S1 $C_{n}$ and $C_{nv}$ symmetric spin-momentum locking terms in 2D

In this section, we study how the point groups in two dimensions (2D)
constrain the spin-momentum locking term. In 2D, point groups are classified
into $C_{n}$ and $C_{nv}$ with $n=1,2,3,4,6$. $C_{n}$ is the $n$-fold rotation
symmetry around the axis perpendicular to the physical 2D plane, and the
$C_{nv}$ symmetry group is generated by $C_{n}$ rotation and mirror symmetry
along the $x$-direction, $M_{x}$.

###  S1.1 General formula for leading terms

We start by considering the most general spin-momentum locking term in 2D
which takes the form

$\displaystyle
d_{x}(k_{x},k_{y})S_{X}+d_{y}(k_{x},k_{y})S_{Y}+d_{z}(k_{x},k_{y})S_{Z}.$ (S1)

When time-reversal symmetry (TRS) $\mathcal{T}$ exists, $d_{x,y,z}$ must be
odd in $k_{x,y}$. In the presence of $C_{2,4,6}$ and $C_{2v,4v,6v}$
symmetries, $d_{z}$ should be even in $k_{x,y}$. This implies that $d_{z}=0$
if we have both $C_{2,4,6}$ ($C_{2v,4v,6v}$) and TRS. For $C_{3}$ and $C_{3v}$
symmetric systems, $d_{z}$ can be nonzero but the leading term is cubic in
momentum: In the $C_{3}$ symmetric case, $d_{z}=k_{+}^{3}+k_{-}^{3}$ and/or
$d_{z}=i(k_{+}^{3}-k_{-}^{3})$ where $k_{\pm}=k_{x}\pm ik_{y}$. In the case of
$C_{3v}$, because $M_{x}$ flips $S_{Z}$, $d_{z}$ can only be
$k_{+}^{3}+k_{-}^{3}$ that is odd in $k_{x}$ Fu (2009); Samokhin (2022). On
the other hand, in $C_{1}$ and $C_{1v}$ symmetric systems, $d_{z}$ can be
nonzero and linear in $k_{x,y}$. In summary, for $C_{n}$ and $C_{nv}$ groups
with $n\geq 2$, the most general spin-momentum locking terms that are linear
in $k_{x,y}$ (i.e., at lowest order) Samokhin (2022) take the form

$\displaystyle h({\mathbf{k}},{\mathbf{S}})=\sum_{i,A}g_{iA}k_{i}S_{A},$ (S2)

which is Eq. (1) of the main text. Here, $i=x,y$ and $A=X,Y$,
${\mathbf{S}}=(S_{X},S_{Y})$. Under different symmetry groups, the allowed
$g_{iA}$’s are different. We will discuss the symmetry constraints for
$g_{iA}$ in the next subsection.

###  S1.2 Symmetry constraints of $g_{iA}$

Consider a symmetry operation $\sigma$, which acts on real-space coordinate
${\mathbf{r}}$: ${\mathbf{r}}\to
O_{\sigma}{\mathbf{r}}+{\mathbf{\delta}}_{\sigma}$. Here, $O_{\sigma}$ is an
element in $O(3)$. Hence, $O_{\sigma}=l_{\sigma}R_{\sigma}$. Note
$l_{\sigma}={\rm Det}O_{\sigma}\in\\{+1,-1\\}$. $R_{\sigma}$, which is $SO(3)$
part of $O_{\sigma}$, can be characterized by $({\mathbf{n}},\phi)$, rotation
axis and angle. Thus, momentum and spin transform under $\sigma$ as
${\mathbf{k}}\to O_{\sigma}{\mathbf{k}}$ and ${\mathbf{S}}\to
R_{\sigma}{\mathbf{S}}$. Thus, to be
$h({\mathbf{k}},{\mathbf{S}})=h(O_{\sigma}{\mathbf{k}},R_{\sigma}{\mathbf{S}})$,
$g_{iA}$ must satisfy

$\displaystyle
g_{iA}=l_{\sigma}[O_{\sigma}]_{ij}g_{jB}[O_{\sigma}^{-1}]_{BA}.$ (S3)

From Eq. (S3), let us consider the constraints by $M_{x}$ and $C_{4}$. Note
that $M_{x}:(x,y)\to(-x,y)$ with $O_{M_{x}}={\rm Diag}(-1,1)$ and
$C_{4}:(x,y)\to(-y,x)$ with $O_{C_{4}}=\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}$. Thus, $g_{iA}$ takes the form
$g_{iA}=\begin{pmatrix}0&a_{1}\\\ a_{2}&0\end{pmatrix}_{iA}$ under $M_{x}$ and
$g_{iA}=\begin{pmatrix}a_{1}&a_{2}\\\ -a_{2}&a_{1}\end{pmatrix}_{iA}$ under
$C_{4}$. More explicitly, $h({\mathbf{k}},{\mathbf{S}})$ are constrained by
$M_{x}$ and $C_{4}$ as

$\displaystyle M_{x}:a_{1}k_{x}S_{Y}+a_{2}k_{y}S_{X}$ (S4) $\displaystyle
C_{4}:a_{1}(k_{x}S_{X}+k_{y}S_{Y})+a_{2}(k_{x}S_{Y}-k_{y}S_{X}).$ (S5)

When $C_{4v}+\mathcal{T}$ exists, then

$\displaystyle h({\mathbf{k}},{\mathbf{S}})=v(k_{x}S_{Y}-k_{y}S_{X}),$ (S6)

which was our primary focus in the main text.

##  S2 Perturbations effect on the commuting-spin Rashba mode

In this section, we discuss the symmetric perturbation terms for the CSR mode.
Up to the quadratic order in ${\mathbf{k}}$, there exist only two perturbation
terms that preserve both $C_{4v}+\mathcal{T}$ and spin-$U(1)$ symmetries, and
they can be added into the CSR Hamiltonian:

$\displaystyle h_{\mathbf{k}}=$ $\displaystyle vk_{x}S_{Y}-vk_{y}S_{X}$
$\displaystyle+\alpha(k_{x}^{2}+k_{y}^{2})\mathds{1}_{4\times 4}+\beta
k_{x}k_{y}S_{X}S_{Y}.$ (S7)

The role of the first term in the second line closes the Fermi surface, as
discussed in the main text. $h_{\mathbf{k}}$ has four energy bands
$E_{{\mathbf{k}}}^{1,2}=\alpha(k_{x}^{2}+k_{y}^{2})\pm v(k_{x}-k_{y})+\beta
k_{x}k_{y}$ and $E_{{\mathbf{k}}}^{3,4}=\alpha(k_{x}^{2}+k_{y}^{2})\pm
v(k_{x}+k_{y})-\beta k_{x}k_{y}$. The Fermi surfaces of these four bands are
generally conic sections that are closed when $|\alpha|>|\beta|/2$ if
${\mathbf{k}}$ can be arbitrarily large. More precisely, in order for theses
Fermi surfaces to be closed, the momentum cutoff $\Lambda$ should be larger
than a bound,
$2(|v(\alpha+\beta/2)|+\sqrt{v^{2}(\alpha+\beta/2)^{2}+(\alpha^{2}-\beta^{2}/4)(v^{2}+4\mu\alpha)})/(4\alpha^{2}-\beta^{2})$,
where $\mu$ is the chemical potential.

For the surface CSR mode of the 3D topological semimetal [c.f. Eq. (11) in the
main text], the two allowed perturbation terms can be added through bulk terms
$(2-\cos k_{x}-\cos k_{y})\mu_{i}\mathbbm{1}_{4\times 4}$ and $\sin k_{x}\sin
k_{y}\mu_{i}S_{X}S_{Y}$, where $i=0,x,y,z$ and $\mu_{0}$ is the $2\times 2$
identity matrix. The terms with $\mu_{0}$ only shift the energy of the Weyl
points, and the terms with $\mu_{x,y,z}$ can move and/or split the Weyl points
in the momentum space.

##  S3 Bulk responses

###  S3.1 Derivation of linear responses

To study the mixed spin-charge linear response, we start from the Peierls
substitution to derive the current operator. With charge and spin $U(1)$ gauge
fields, the Hamiltonian can be expressed as

$\displaystyle H=$
$\displaystyle\sum_{\mathbf{R}\mathbf{R}^{\prime}\alpha\beta}\,h_{\mathbf{R}\mathbf{R}^{\prime}}^{\alpha\beta}\,c_{\mathbf{R}\alpha}^{{\dagger}}c_{\mathbf{R}^{\prime}\beta}\,e^{i\phi_{\mathbf{R}\alpha,\mathbf{R^{\prime}}\beta}},$
(S8)

where $\mathbf{R},\mathbf{R}^{\prime}$ are coordinates of lattice sites, and
$\mathbf{x}_{\alpha}$ denotes the position of $\alpha$-th orbital within unit
cell. The quantity $\phi_{\mathbf{R}\alpha,\mathbf{R^{\prime}}\beta}$ is the
Peierls phase,

$\displaystyle\phi_{\mathbf{R}\alpha,\mathbf{R^{\prime}}\beta}=\int_{\mathbf{R}^{\prime}+\mathbf{x}_{\beta}}^{\mathbf{R}+\mathbf{x}_{\alpha}}(-e\mathbf{A}+S^{\alpha\beta}_{A}\mathbf{a}^{A})\cdot
d\mathbf{l}.$ (S9)

Note that we set $\hbar=c=1$ and the charge of electron to be $-e$.

With spin $U(1)$ symmetries, $h_{\mathbf{R}\mathbf{R}^{\prime}}^{\alpha\beta}$
is nonzero only when $\alpha$ and $\beta$ are in the same spin sector. With
translation symmetry, $h^{\alpha\beta}_{\mathbf{R},\mathbf{R}^{\prime}}$
depends on only the distance between two sites, i.e.
$h^{\alpha\beta}_{\mathbf{R}\mathbf{R}^{\prime}}=h^{\alpha\beta}_{\mathbf{R}-\mathbf{R}^{\prime}}\equiv
h^{\alpha\beta}_{\bm{\Delta\mathbf{R}}}$. We Fourier transform and define

$\displaystyle
c_{\mathbf{R}\alpha}^{{\dagger}}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}e^{-i{\mathbf{k}}\cdot(\mathbf{R}+\mathbf{x}_{\alpha})}c^{{\dagger}}_{{\mathbf{k}}\alpha},$
(S10) $\displaystyle
h^{\alpha\beta}_{{\mathbf{k}}}=e^{-i\bf{k}\cdot(\mathbf{x}_{\alpha}-\mathbf{x}_{\beta})}\sum_{\Delta\mathbf{R}}e^{-i{\mathbf{k}}\cdot\Delta\mathbf{R}}h^{\alpha\beta}_{\Delta\mathbf{R}},$
(S11)

where $N$ is the number of lattice sites. Then, up to second order in gauge
fields, Eq. (S8) can be written as

$\displaystyle H\simeq$
$\displaystyle\frac{1}{N}\sum_{\Delta\mathbf{R}\mathbf{R}^{\prime}\alpha\beta}\sum_{{\mathbf{k}}_{1}{\mathbf{k}}_{2}}e^{-i{\mathbf{k}}_{1}\cdot\Delta\mathbf{R}}e^{i({\mathbf{k}}_{2}-{\mathbf{k}}_{1})\cdot\mathbf{R}^{\prime}}h_{\Delta\mathbf{R}}^{\alpha\beta}c^{{\dagger}}_{{\mathbf{k}}_{1}\alpha}c_{{\mathbf{k}}_{2}\beta}$
$\displaystyle\times
e^{-i{\mathbf{k}}_{1}\cdot{\mathbf{x}}_{\alpha}+i{\mathbf{k}}_{2}\cdot{\mathbf{x}}_{\beta}}\sum_{n=0}^{2}\frac{1}{n!}(i\phi_{(\mathbf{R}^{\prime}+\Delta\mathbf{R})\alpha,\mathbf{R}^{\prime}\beta})^{n}$
(S12) $\displaystyle\equiv H_{0}+H_{\rm ext}.$ (S13)

We use an approximation

$\displaystyle\phi_{{\mathbf{R}}\alpha,{\mathbf{R}}^{\prime}\beta}\simeq\Delta{\mathbf{r}}\cdot(-e{\mathbf{A}}({\mathbf{r}})+S^{\alpha\beta}_{A}{\mathbf{a}}^{A}({\mathbf{r}}))_{{\mathbf{r}}={\mathbf{r}}_{\rm
mid}}$ (S14)

where $\Delta{\mathbf{r}}={\mathbf{r}}_{1}-{\mathbf{r}}_{2}$ and
${\mathbf{r}}_{\rm mid}=\frac{1}{2}({\mathbf{r}}_{1}+{\mathbf{r}}_{2})$ for
${\mathbf{r}}_{1}={\mathbf{R}}+{\mathbf{x}}_{\alpha}$,
${\mathbf{r}}_{2}={\mathbf{R}}^{\prime}+{\mathbf{x}}_{\beta}$. We also define
the Fourier transform of gauge fields,

$\displaystyle\mathbf{A}(\mathbf{r})=\frac{1}{N}\sum_{{\mathbf{k}}}e^{-i{\mathbf{k}}\cdot\mathbf{r}}\mathbf{A}_{\mathbf{k}},$
(S15)
$\displaystyle\mathbf{a}^{A}(\mathbf{r})=\frac{1}{N}\sum_{{\mathbf{k}}}e^{-i{\mathbf{k}}\cdot\mathbf{r}}\mathbf{a}^{A}_{\mathbf{k}}.$
(S16)

By substituting these equations in Eq. (S3.1), the term linear in gauge fields
becomes

$\displaystyle H_{\rm
ext}^{(1)}=\frac{i}{N}\sum_{\Delta{\mathbf{R}}{\mathbf{k}}{\mathbf{q}}}\sum_{\alpha\beta}e^{-i({\mathbf{k}}+\frac{{\mathbf{q}}}{2})\cdot(\Delta{\mathbf{R}}+{\mathbf{x}}_{\alpha}-{\mathbf{x}}_{\beta})}h^{\alpha\beta}_{\Delta{\mathbf{R}}}c_{{\mathbf{k}}\alpha}^{{\dagger}}c_{{\mathbf{k}}+{\mathbf{q}}\beta}$
$\displaystyle\times(\Delta{\mathbf{R}}+{\mathbf{x}}_{\alpha}-{\mathbf{x}}_{\beta})\cdot(-e\mathbf{A}_{\mathbf{q}}+S^{\alpha\beta}_{A}\mathbf{a}^{A}_{\mathbf{q}})$
$\displaystyle=-\frac{1}{N}\sum_{{\mathbf{k}}{\mathbf{q}}}\sum_{\alpha\beta}c^{{\dagger}}_{{\mathbf{k}}-\frac{{\mathbf{q}}}{2}\alpha}c_{{\mathbf{k}}+\frac{{\mathbf{q}}}{2}\beta}\frac{\partial
h^{\alpha\beta}_{{\mathbf{k}}}}{\partial{\mathbf{k}}}\cdot(-e\mathbf{A}_{{\mathbf{q}}}+S^{\alpha\beta}_{A}\mathbf{a}^{A}_{{\mathbf{q}}}),$
(S17)

and the term quadratic in gauge fields becomes

$\displaystyle H_{\rm
ext}^{(2)}=\frac{1}{2N^{2}}\sum_{{\mathbf{k}}{\mathbf{q}}{\mathbf{q}}^{\prime}}\sum_{\alpha\beta}c^{{\dagger}}_{{\mathbf{k}}-\frac{{\mathbf{q}}+{\mathbf{q}}^{\prime}}{2}\alpha}c_{{\mathbf{k}}+\frac{{\mathbf{q}}+{\mathbf{q}}^{\prime}}{2}\beta}\frac{\partial^{2}h^{\alpha\beta}_{{\mathbf{k}}}}{\partial
k_{i}\partial k_{j}}$
$\displaystyle\times(-e{\mathbf{A}}_{{\mathbf{q}}}+S^{\alpha\beta}_{A}a^{A}_{{\mathbf{q}}})_{i}(-e{\mathbf{A}}_{{\mathbf{q}}^{\prime}}+S^{\alpha\beta}_{A}a^{A}_{{\mathbf{q}}^{\prime}})_{j}.$
(S18)

We define the total charge and spin current at momentum $\mathbf{q}$ through

$\displaystyle\mathbf{J}_{\text{tot}}(\mathbf{r})=\frac{1}{N}\sum_{\mathbf{q}}e^{-i\mathbf{q}\cdot\mathbf{r}}\mathbf{J}_{\text{tot}}(\mathbf{q}),$
(S19)
$\displaystyle\bm{\mathcal{J}}_{A,\text{tot}}(\mathbf{r})=\frac{1}{N}\sum_{\mathbf{q}}e^{-i\mathbf{q}\cdot\mathbf{r}}\bm{\mathcal{J}}_{A,\text{tot}}(\mathbf{q}).$
(S20)

The current operators, ${\mathbf{J}}_{\rm tot}({\mathbf{q}})$ and
${\mathbf{\mathcal{missing}}}{J}_{A,\rm tot}({\mathbf{q}})$, can be derived
from $H_{\rm ext}=H^{(1)}_{\rm ext}+H^{(2)}_{\rm ext}$ through

$\displaystyle J^{i}_{\rm tot}({\mathbf{q}})=-\frac{N}{V_{\rm
cell}}\frac{\delta H_{\rm ext}}{\delta A_{i,-{\mathbf{q}}}},$ (S21)
$\displaystyle\mathcal{J}^{i}_{A,\rm tot}({\mathbf{q}})=-\frac{N}{V_{\rm
cell}}\frac{\delta H_{\rm ext}}{\delta a^{A}_{i,-{\mathbf{q}}}},$ (S22)

where $V_{\rm cell}$ is the volume of the unit cell, and $\delta
A_{i,{\mathbf{q}}}/\delta
A_{j,{\mathbf{q}}^{\prime}}=\delta_{ij}\delta_{{\mathbf{q}}{\mathbf{q}}^{\prime}}$.
We now consider uniform gauge fields such that
$A_{i,{\mathbf{q}}\neq{\mathbf{0}}}=0$ and
$a^{A}_{i,{\mathbf{q}}\neq{\mathbf{0}}}=0$, and apply Eqs. (S3.1), (S3.1) and
(S21), we have

$\displaystyle J^{i}_{\rm tot}=\frac{1}{V_{\rm
cell}}\sum_{{\mathbf{k}}\alpha\beta}c_{{\mathbf{k}}\alpha}^{{\dagger}}(J^{i}({\mathbf{k}})+J^{i}_{\rm
diam}({\mathbf{k}}))_{\alpha\beta}c_{{\mathbf{k}}\beta},$ (S23)

where $J^{i}_{\rm tot}$ is short for $J^{i}_{\rm
tot}({\mathbf{q}}={\mathbf{0}})$, and

$\displaystyle J^{i}({\mathbf{k}})=-e\partial_{i}h_{{\mathbf{k}}},$ (S24)
$\displaystyle J^{i}_{\rm
diam}({\mathbf{k}})=\frac{e}{N}\partial_{i}\partial_{j}h_{{\mathbf{k}}}(-eA_{j}+S_{A}a^{A}_{j}),$
(S25)

are the paramagnetic current and diamagnetic current, respectively. Note that
$\partial_{i}$ is short for $\partial_{k_{i}}$, and $A_{i}$ ($a^{A}_{i}$) is
short for $A_{i,{\mathbf{q}}={\mathbf{0}}}$
($a_{i,{\mathbf{q}}={\mathbf{0}}}^{A}$). Similarly, we have

$\displaystyle\mathcal{J}^{i}_{A,\rm tot}=\frac{1}{V_{\rm
cell}}\sum_{{\mathbf{k}}\alpha\beta}c_{{\mathbf{k}}\alpha}^{{\dagger}}(\mathcal{J}^{i}_{A}({\mathbf{k}})+\mathcal{J}^{i}_{A,\rm
diam}({\mathbf{k}}))_{\alpha\beta}c_{{\mathbf{k}}\beta},$ (S26)

where $\mathcal{J}^{i}_{A,\rm tot}$ is short for $\mathcal{J}^{i}_{A,\rm
tot}({\mathbf{q}}={\mathbf{0}})$, and

$\displaystyle\mathcal{J}^{i}_{A}({\mathbf{k}})=S_{A}\partial_{i}h_{{\mathbf{k}}},$
(S27) $\displaystyle\mathcal{J}^{i}_{A,\rm
diam}({\mathbf{k}})=-\frac{1}{N}S_{A}\partial_{i}\partial_{j}h_{{\mathbf{k}}}(-eA_{j}+S_{A}a^{A}_{j}).$
(S28)

Then the conductivity for mutual charge-spin response at zero frequency is
given by

$\displaystyle\sigma_{A}^{ij}=\lim_{\nu_{m}\rightarrow
0}\frac{\Pi^{ij}_{A}(i\nu_{m})}{\nu_{m}}+\sigma_{A,\rm diam}^{ij},$ (S29)

where

$\displaystyle\Pi^{ij}_{A}(i\nu_{m})=$ $\displaystyle\frac{1}{\beta
V}\sum_{{\mathbf{k}},n}\,{\rm
Tr}[\mathcal{J}^{i}_{A}({\mathbf{k}})G({\mathbf{k}},i\omega_{n}+i\nu_{m})$
$\displaystyle\times J^{j}({\mathbf{k}})G({\mathbf{k}},i\omega_{n})]$ (S30)

is the Kubo formula. Here, $\beta=(k_{B}T)^{-1}$. $V=NV_{\rm cell}$ is the
volume of the whole system. $\omega_{n}=(2n+1)\pi/\beta$ is the Matsubara
frequency for fermion. Also,
$G({\mathbf{k}},i\omega_{n})=(i\omega_{n}-h_{{\mathbf{k}}})^{-1}$ is the Green
function. Meanwhile, the diamagnetic part $\sigma_{A,\rm diam}^{ij}$ takes the
form

$\displaystyle\sigma_{A,\rm diam}^{ij}=\lim_{\omega\rightarrow
0}\frac{1}{V}\sum_{{\mathbf{k}},a\in\rm occ}\frac{e}{i\omega}\langle
a,{\mathbf{k}}|S_{A}\partial_{i}\partial_{j}h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle$
$\displaystyle=\lim_{\omega\rightarrow 0}\sum_{a\in\rm
occ}\frac{e}{i\omega}\int\frac{d^{d}{\mathbf{k}}}{(2\pi)^{d}}\langle
a,{\mathbf{k}}|S_{A}\partial_{i}\partial_{j}h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle,$
(S31)

which is divergent. Importantly, we will show later that this divergence in
$\sigma_{A,\rm diam}^{ij}$ cancels the divergence in the symmetric
contribution to the off-diagonal $\sigma_{xy}$ response in the Kubo formula
[c.f. Eq. (S3.1)].

To calculate the ${\rm Tr}[\dots]$ in Eq. (S3.1), we introduce energy
eigenstates $|a,{\mathbf{k}}\rangle$ such that
$h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle=\epsilon_{a}({\mathbf{k}})|a,{\mathbf{k}}\rangle$.
Then,

$\displaystyle G({\mathbf{k}},i\omega_{n})=\sum_{a\in{\rm
all}}|a,{\mathbf{k}}\rangle\langle
a,{\mathbf{k}}|\,[i\omega_{n}-\epsilon_{a}({\mathbf{k}})]^{-1}$ (S32)
$\displaystyle{\rm Tr}[\dots]=-e\sum_{a,b\in{\rm all}}\langle
a,{\mathbf{k}}|S_{A}\partial_{i}h_{{\mathbf{k}}}|b,{\mathbf{k}}\rangle\langle
b,{\mathbf{k}}|\partial_{j}h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle$
$\displaystyle\times[i\omega_{n}-\epsilon_{a}({\mathbf{k}})]^{-1}[i\omega_{n}+i\nu_{m}-\epsilon_{b}({\mathbf{k}})]^{-1}$
(S33)

Note that

$\displaystyle\frac{1}{\beta}\sum_{n}\,[i\omega_{n}-\epsilon_{a}({\mathbf{k}})]^{-1}[i\omega_{n}+i\nu_{m}-\epsilon_{b}({\mathbf{k}})]^{-1}$
(S34)
$\displaystyle=\frac{\theta(E_{F}-\epsilon_{a}({\mathbf{k}}))-\theta(E_{F}-\epsilon_{b}({\mathbf{k}}))}{i\nu_{m}-\epsilon_{b}({\mathbf{k}})+\epsilon_{a}({\mathbf{k}})}$
(S35)

where the Fermi energy $E_{F}$. Since $[h_{{\mathbf{k}}},S_{A}]=0$ and
$h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle=\epsilon_{a}({\mathbf{k}})|a,{\mathbf{k}}\rangle$,
$\langle
a,{\mathbf{k}}|S_{A}\partial_{i}h_{{\mathbf{k}}}|b,{\mathbf{k}}\rangle=0$
unless the spin eigenvalues of states $a$ and $b$ are the same, i.e.
$S_{a,A}=S_{b,A}$. We define $v_{iab}({\mathbf{k}})=\langle
a,{\mathbf{k}}|\partial_{i}h_{{\mathbf{k}}}|b,{\mathbf{k}}\rangle$, and then
$\Pi_{A}^{ij}(i\nu_{m})$ can be expressed as

$\displaystyle\frac{e}{V}\sum_{{\mathbf{k}},a,\bar{b}}\left[\frac{v_{ja\bar{b}}({\mathbf{k}})S_{a,A}v_{i\bar{b}a}({\mathbf{k}})}{i\nu_{m}+\Delta_{\bar{b}a}({\mathbf{k}})}-\frac{S_{a,A}v_{ia\bar{b}}({\mathbf{k}})v_{j\bar{b}a}({\mathbf{k}})}{i\nu_{m}-\Delta_{\bar{b}a}({\mathbf{k}})}\right],$
(S36)

where $a$ labels eigenstates in the occupied subspace, $\bar{b}$ labels
eigenstates in unoccupied subspace, and
$\Delta_{\bar{b}a}({\mathbf{k}})=\epsilon_{\bar{b}}({\mathbf{k}})-\epsilon_{a}({\mathbf{k}})$.
The symmetric part of Eq. (S36) is then

$\displaystyle\frac{e}{V}\sum_{{\mathbf{k}},a,\bar{b}}\frac{\Delta_{\bar{b}a}({\mathbf{k}})}{\nu_{m}^{2}+\Delta^{2}_{\bar{b}a}}S_{a,A}[v_{ia\bar{b}}({\mathbf{k}})v_{j\bar{b}a}({\mathbf{k}})+(i\leftrightarrow
j)].$ (S37)

Replacing $\nu_{m}$ by $-i\omega$ and replacing the sum over ${\mathbf{k}}$ by
an integral, the symmetric part of $\Pi_{A}^{ij}$ contributes to the charge-
spin mixed conductivity as

$\displaystyle\lim_{\omega\rightarrow 0}\sum_{a,\bar{b}}$
$\displaystyle-\frac{e}{i\omega}\int\frac{d{\mathbf{k}}}{(2\pi)^{d}}\Delta_{\bar{b}a}({\mathbf{k}})S_{a,A}$
$\displaystyle\times[A_{ia\overline{b}}({\mathbf{k}})A_{j\overline{b}a}({\mathbf{k}})+(i\leftrightarrow
j)].$ (S38)

where we define the non-Abelian Berry connection

$\displaystyle A_{iab}({\mathbf{k}})=-i\langle
a,{\mathbf{k}}|\partial_{i}|b,{\mathbf{k}}\rangle$ (S39)

and
$v_{ia\bar{b}}({\mathbf{k}})=i\Delta_{\bar{b}a}({\mathbf{k}})A_{ia\overline{b}}({\mathbf{k}})$
is used. Next, we show explicitly that Eq. (S3.1) cancels Eq. (S38). First, we
use $\sum_{a\in\rm occ}\langle
a,{\mathbf{k}}|S_{A}\partial_{i}\partial_{j}h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle$
in Eq. (S3.1) is equal to

$\displaystyle\sum_{a\in\rm
occ}S_{a,A}[\partial_{i}v_{jaa}({\mathbf{k}})-\langle\partial_{i}a,{\mathbf{k}}|\partial_{j}h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle$
$\displaystyle-\langle
a,{\mathbf{k}}|\partial_{j}h_{{\mathbf{k}}}|\partial_{i}a,{\mathbf{k}}\rangle],$
(S40)

where the first term vanishes after integrating over ${\mathbf{k}}$ due to the
periodicity. Ignoring the first term, the above equation becomes

$\displaystyle i\sum_{a,a^{\prime}\in\rm
occ}S_{a,A}[A_{iaa^{\prime}}({\mathbf{k}})v_{ja^{\prime}a}({\mathbf{k}})-v_{jaa^{\prime}}({\mathbf{k}})A_{ia^{\prime}a}({\mathbf{k}})]$
$\displaystyle+i\sum_{a,\overline{b}}S_{a,A}[A_{ia\overline{b}}({\mathbf{k}})v_{j\overline{b}a}({\mathbf{k}})-v_{ja\overline{b}}({\mathbf{k}})A_{i\overline{b}a}({\mathbf{k}})],$
(S41)

where the first two lines vanish because $a$ and $a^{\prime}$ are symmetric.
Using
$v_{ja\bar{b}}({\mathbf{k}})=i\Delta_{\bar{b}a}({\mathbf{k}})A_{ja\overline{b}}({\mathbf{k}})$
again, we eventually have

$\displaystyle\sum_{a\in\rm occ}\langle
a,{\mathbf{k}}|S_{A}\partial_{i}\partial_{j}h_{{\mathbf{k}}}|a,{\mathbf{k}}\rangle$
$\displaystyle=\sum_{a,\bar{b}}\Delta_{\bar{b}a}({\mathbf{k}})S_{a,A}[A_{ia\overline{b}}({\mathbf{k}})A_{j\overline{b}a}({\mathbf{k}})+(i\leftrightarrow
j)].$ (S42)

Substituting the above equation into Eq. (S3.1), it is straightforward to see
that $\sigma_{A,\rm diag}^{ij}$ cancels the symmetric part from the Kubo
formula.

Therefore, only the antisymmetric part of the Kubo formula is important for
the linear response, which takes the form

$\displaystyle\Pi^{ij}_{A}(i\nu_{m})=\frac{e}{V}\sum_{{\mathbf{k}},a,\overline{b}}\,\frac{i\nu_{m}}{[\epsilon_{a}({\mathbf{k}})-\epsilon_{\overline{b}}({\mathbf{k}})]^{2}}S_{a,A}$
$\displaystyle\times[v_{ia\overline{b}}({\mathbf{k}})v_{j\overline{b}a}({\mathbf{k}})-(i\leftrightarrow
j)].$ (S43)

Note that the Abelian Berry connection and curvature are defined by

$\displaystyle A_{i,a}({\mathbf{k}})=A_{iaa}({\mathbf{k}}),$ (S44)
$\displaystyle
F^{a}_{ij}({\mathbf{k}})=\partial_{i}A_{j,a}({\mathbf{k}})-\partial_{j}A_{i,a}({\mathbf{k}})$
$\displaystyle=-i\sum_{\overline{b}\notin\rm
occ}\left[A_{ia\overline{b}}({\mathbf{k}})A_{j\overline{b}a}({\mathbf{k}})-(i\leftrightarrow
j)\right].$ (S45)

Then, Eq. (S43) can be simplified to

$\displaystyle\Pi^{ij}_{A}(i\nu_{m})=-\frac{e}{V}\nu_{m}\sum_{{\mathbf{k}},a\in\rm
occ}S_{a,A}F_{ij}^{a}({\mathbf{k}}).$ (S46)

Thus, the conductivity $\sigma_{A}^{ij}$ becomes

$\displaystyle\sigma^{ij}_{A}=$ $\displaystyle\lim_{\nu_{m}\to
0}\frac{1}{\nu_{m}}\Pi^{ij}_{A}({\mathbf{k}})$ (S47) $\displaystyle=$
$\displaystyle-\frac{e}{(2\pi)^{3}}\sum_{a\in\text{occ}}\int
d{\mathbf{k}}S_{a,A}F_{ij}^{a}({\mathbf{k}}).$ (S48)

From the Stokes theorem,

$\displaystyle\int
d{\mathbf{k}}\sum_{a\in\text{occ}}S_{a,A}F_{ij}^{a}({\mathbf{k}})$
$\displaystyle=\sum_{\gamma}\int_{FS_{\gamma}}\epsilon_{ijk}k_{k}\sum_{a\in\text{occ}}S_{a,A}\left(\frac{1}{2}F_{\mu\nu}^{a}({\mathbf{k}})dS_{\mu}\wedge
dS_{\nu}\right)$
$\displaystyle=2\pi\epsilon_{ijk}\sum_{\alpha}\chi_{\alpha}K_{\alpha,k}S_{\alpha,A}=2\pi\epsilon_{ijk}Q_{kA}.$
(S49)

Here $FS_{\gamma}$ denotes a set of disjoint Fermi surfaces. The index
$\alpha=(\gamma,a)$ in the last line labels both the pseudospin sector and the
location in momentum space of the Fermi surfaces. Note that we set the Fermi
level slightly above/below the Weyl points. Hence, we arrive at

$\displaystyle\sigma^{ij}_{A}=-\frac{e}{4\pi^{2}}\epsilon_{ijk}Q_{kA}.$ (S50)

Then, from $J^{i}_{\rm bulk}=\sigma^{ij}_{A}\mathcal{E}^{j}_{\alpha}$ and
$\mathcal{J}^{i}_{A}=\sigma^{ij}_{A}E_{j}$, it is straightforward to derive
the bulk responses shown in the main text.

###  S3.2 Symmetry constraints on $Q_{iA}$

In the main text, we showed that the Weyl spin-momentum quadrupole $Q_{iA}$
plays the role of a response coefficient. In this subsection we discuss the
symmetry constraints on $Q_{iA}$. We first note that the chiral charge
$\chi_{\alpha}$, momentum ${\mathbf{K}}_{\alpha}$, and pseudospin
${\mathbf{S}}_{\alpha}$ of Weyl point $\alpha$ transform like a pseudoscalar,
a vector, and an axial vector, respectively. Thus, for a symmetry $\sigma$
with $O(3)$ element $O_{\sigma}$ (see Sec. S1 for the definition of
$O_{\sigma}$), we derive

$\displaystyle Q_{iA}=[O_{\sigma}]_{ij}Q_{jB}[O_{\sigma}^{-1}]_{BA}.$ (S51)

For example, mirror $M_{x}$ and $C_{n}$ rotation impose the following
constraints on $Q_{iA}$:

$\displaystyle M_{x}:Q_{iA}=\begin{pmatrix}a_{1}&0\\\
0&a_{2}\end{pmatrix}_{iA},$ (S52) $\displaystyle
C_{2}:Q_{iA}=\begin{pmatrix}a_{1}&a_{2}\\\ a_{3}&a_{4}\end{pmatrix}_{iA},$
(S53) $\displaystyle C_{n>2}:Q_{iA}=\begin{pmatrix}a_{1}&a_{2}\\\
-a_{2}&a_{1}\end{pmatrix}_{iA},$ (S54)

where $a_{1,2,3,4}\in\mathbb{R}$. This implies that $Q_{iA}$ is proportional
to the identity matrix when $C_{nv}$ ($n>2$) exists, in a suitable coordinate
system.

##  S4 Anomalies of the generalized spin-momentum locked mode

###  S4.1 Current anomaly

In this section, we explicitly write down the anomalous conservation laws for
the generalized Hamiltonian in Eq. (S2). To simplify the equations, we define
$g_{i\pm}=g_{iX}\pm g_{iY}$, and $\sigma_{i\pm}={\rm sgn}(g_{i\pm})$. The
anomalous conservation laws for charge and spin currents can be obtained by
inserting gauge fluxes in the same way described in the main text. Then, the
divergence of charge and spin currents are given by

$\displaystyle\partial_{\mu}\mathcal{J}^{\mu}_{X}=-\frac{e\Lambda_{y}}{\pi^{2}}(\sigma_{x+}+\sigma_{x-})E_{x}-\frac{e\Lambda_{x}}{\pi^{2}}(\sigma_{y+}+\sigma_{y-})E_{y},$
$\displaystyle\partial_{\mu}\mathcal{J}^{\mu}_{Y}=-\frac{e\Lambda_{y}}{\pi^{2}}(\sigma_{x+}-\sigma_{x-})E_{x}-\frac{e\Lambda_{x}}{\pi^{2}}(\sigma_{y+}-\sigma_{y-})E_{y},$
$\displaystyle\partial_{\mu}J^{\mu}=-\frac{e\Lambda_{x}}{\pi^{2}}\large[\sigma_{y+}(\mathcal{E}^{X}_{y}+\mathcal{E}^{Y}_{y})+\sigma_{y-}(\mathcal{E}^{X}_{y}-\mathcal{E}^{Y}_{y})\large]$
$\displaystyle-\frac{e\Lambda_{y}}{\pi^{2}}\large[\sigma_{x+}(\mathcal{E}^{X}_{x}+\mathcal{E}^{Y}_{x})+\sigma_{x-}(\mathcal{E}^{X}_{x}-\mathcal{E}^{Y}_{x})\large].$
(S55)

The last equation can be further simplified if we introduce
$\mathcal{E}^{\pm}_{i}=\mathcal{E}^{X}_{i}\pm\mathcal{E}^{Y}_{i}$:

$\displaystyle\partial_{\mu}J^{\mu}=-\frac{e\Lambda_{x}}{\pi^{2}}(\sigma_{y+}\mathcal{E}^{+}_{y}+\sigma_{y-}\mathcal{E}^{-}_{y})$
$\displaystyle-\frac{e\Lambda_{y}}{\pi^{2}}(\sigma_{x+}\mathcal{E}^{+}_{x}+\sigma_{x-}\mathcal{E}^{-}_{x}).$
(S56)

This reproduce the anomalous conservation laws discussed in the main text. By
substituting $g_{xY}=1$, $g_{yX}=-1$, and $g_{xX}=g_{yY}=0$,

$\displaystyle\partial_{\mu}\mathcal{J}^{\mu}_{X}=\frac{2e\Lambda_{x}}{\pi^{2}}E_{y},\quad\partial_{\mu}\mathcal{J}^{\mu}_{Y}=-\frac{2e\Lambda_{y}}{\pi^{2}}E_{x},$
$\displaystyle\partial_{\mu}J^{\mu}=\frac{2e\Lambda_{x}}{\pi^{2}}\mathcal{E}^{X}_{y}-\frac{2e\Lambda_{y}}{\pi^{2}}\mathcal{E}^{Y}_{x}.$
(S57)

Note that $\Lambda_{x}=\Lambda_{y}=\Lambda$ is set in the main text because of
$C_{4v}$ symmetry.

###  S4.2 Anomaly cancellation and spin-momentum quadrupole

Generalized spin-momentum locked modes cannot exist solely because of the
current anomalies. As discussed in the main text, they can live on the surface
of 3D semimetals with proper spin-momentum quadrupole $Q_{iA}$. We now discuss
what the word _proper_ means in terms of the Callan-Harvey anomaly
cancellation mechanism Callan Jr and Harvey (1985). The bulk response term
determined by the spin-momentum quadrupole generally takes the form

$\displaystyle S_{CS}[A,a]=\int d^{4}x\alpha\epsilon^{\mu\nu\rho\sigma}Q_{\mu
A}\partial_{\nu}a^{A}_{\rho}A_{\sigma},$ (S58) $\displaystyle
J^{\mu}=\frac{\delta S}{\delta A_{\mu}}=-\alpha\epsilon^{\mu
i\nu\rho}Q_{iA}\partial_{\nu}a^{A}_{\rho},$ (S59)
$\displaystyle\mathcal{J}^{\mu}_{A}=\frac{\delta S}{\delta
a^{A}_{\mu}}=-\alpha\epsilon^{\mu i\nu\rho}Q_{iA}\partial_{\nu}A_{\rho}.$
(S60)

When the bulk has a constant spin-momentum quadrupole that is denoted as
$Q_{0,iA}$, then the presence of boundary with $Q_{iA}(z)=Q_{0,iA}\theta(z)$,
this contributes
$\partial_{\mu}J_{b}^{\mu}=\int_{0_{-}}^{0^{+}}dz\partial_{\mu}J^{\mu}$ and
$\partial_{\mu}\mathcal{J}_{b,A}^{\mu}=\int_{0_{-}}^{0^{+}}dz\partial_{\mu}\mathcal{J}_{b,A}^{\mu}$.
Hence, considering the contributions from both surface and bulk, we have

$\displaystyle\partial_{\mu}(J_{b}^{\mu}+J_{s}^{\mu})=-f_{iA}\mathcal{E}^{A}_{i},$
(S61) $\displaystyle\partial_{\mu}(J_{b,A}^{\mu}+J_{s,A}^{\mu})=-f_{iA}E_{i},$
(S62)

where the coefficients $f_{iA}$ are given by

$\displaystyle f_{xX}=\alpha
Q_{0,yX}+\frac{e\Lambda_{y}}{\pi^{2}}(\sigma_{x+}+\sigma_{x-}),$ (S63)
$\displaystyle f_{xY}=\alpha
Q_{0,yY}+\frac{e\Lambda_{y}}{\pi^{2}}(\sigma_{x+}-\sigma_{x-}),$ (S64)
$\displaystyle f_{yX}=-\alpha
Q_{0,xX}+\frac{e\Lambda_{x}}{\pi^{2}}(\sigma_{y+}+\sigma_{y-}),$ (S65)
$\displaystyle f_{yY}=-\alpha
Q_{0,xY}+\frac{e\Lambda_{x}}{\pi^{2}}(\sigma_{y+}-\sigma_{y-}).$ (S66)

All of the four $f_{iA}$ must be zero to have the anomaly cancellation between
bulk and surface. This is equivalent to say that generalized spin-momentum
locked modes can only live on the surface of 3D systems with bulk spin-
momentum quadrupole satisfying $f_{iA}=0$.

Finally, we comment that, in the presence of symmetries, $Q_{iA}$ and $g_{iA}$
are further constrained by Eqs. (S3) and (S51).

##  S5 Numerical scheme for linear response

In the main text, we studied the bulk responses under the electromagnetic (EM)
$U(1)_{\rm EM}$ and pseudospin $U(1)_{A=X,Y}$ magnetic fields. Here, we
explain how (pseudospin) magnetic fluxes are implemented in our tight-binding
model.

Figure S1:  A pair of flux tubes. Flux tubes with magnitudes $-\Phi$ and
$\Phi$ are inserted at $(y,z)={\mathbf{r}}_{-}$ and ${\mathbf{r}}_{+}$,
respectively. These flux tubes are stretched along $x$ direction. The red line
$L$ connects two flux tubes. The Peierls substitution is implemented such that
a hopping amplitude $t$ is substituted by $te^{\mp i\Phi}$ when it passes
through $L$ depending on the direction of hopping process $\hat{t}$;
$te^{-i\Phi}$ for $\hat{t}\cdot\hat{y}>0$ and $te^{i\Phi}$ for
$\hat{t}\cdot\hat{y}<0$. Pseudospin magnetic flux $\mathcal{B}^{A}_{i}$ can be
implemented in a similar way. For example, a pair of flux tubes for
$\mathcal{B}^{X}_{x}$ can be described by the following Peierls substitution.
For hopping amplitudes pass through $L$, $t$ is substituted by
$te^{iS_{X}\varphi_{X}}$ if $\hat{t}\cdot\hat{y}>0$ and
$te^{-iS_{X}\varphi_{X}}$ if $\hat{t}\cdot\hat{y}<0$. Here, $\varphi_{X}$ is
the flux of $\mathcal{B}^{X}_{x}$.

First, let us explain the case of $U(1)_{\rm EM}$ magnetic field. The magnetic
field can be defined by spatial components of gauge field $A_{\mu}$. In the
presence of $A_{\mu}$, the hopping amplitudes
$h^{\alpha\beta}_{{\mathbf{R}},{\mathbf{R}}^{\prime}}$ are changed by

$\displaystyle
h^{\alpha\beta}_{{\mathbf{R}},{\mathbf{R}}^{\prime}}e^{-ie\int_{\mathbf{R}^{\prime}+\mathbf{x}_{\beta}}^{\mathbf{R}+\mathbf{x}_{\alpha}}\mathbf{A}\cdot
d\mathbf{l}}$ (S67)

according to the Peierls substitution. Note that we use the same tight-binding
notation introduced in Sec. S3.

Let us consider the magnetic field $B_{x}$ along the $x$ direction. For this,
we introduce a pair of flux tubes with magnitudes $\Phi$ and $-\Phi$. In this
case, we can impose the periodic boundary condition along all three
directions. Moreover, $k_{x}$ remains as well-defined quantum number because
the flux tubes keep the translation symmetry along the $x$ direction. A
schematic description of flux tubes is shown in Fig. S1. The configuration of
gauge field $A_{\mu}$, which is compatible with the flux tubes at
$(y,z)={\mathbf{r}}_{-}$ and ${\mathbf{r}}_{+}$, is described by a string
connecting two flux tubes. (In fact, it is a strip not string because flux
tubes are stretched along the $x$ direction.) For a given string $L$, the
Peierls substitution can be easily implemented. When a hopping amplitude $t$
passes through $L$, then it is substituted by $te^{\mp i\Phi}$ depending on
the direction of hopping process (see Fig. S1 for the details). Note that we
normalize the magnetic flux $\Phi$ such that $\Phi=2\pi$ is equivalent to a
magnetic flux quantum $\phi_{0}=2\pi/e$. Thus, $\Phi=2\pi\phi/\phi_{0}$ where
$\phi$ is the magnetic flux per unit cell. The same normalization is also
applied to the pseudospin magnetic flux.

The pseudospin magnetic flux can be implemented in a similar way. For the
pseudospin gauge field $a^{A}_{i}$, the relevant Peierls substitution is
defined by

$\displaystyle
h^{\alpha\beta}_{{\mathbf{R}},{\mathbf{R}}^{\prime}}e^{ie\int_{{\mathbf{R}}^{\prime}+{\mathbf{x}}_{\beta}}^{{\mathbf{R}}+{\mathbf{x}}_{\alpha}}\sum_{A}S_{A}^{\alpha\beta}{\mathbf{a}}^{A}\cdot
d\mathbf{l}}.$ (S68)

Thus, for a hopping amplitude $t$ passes through $L$, the Peierls substitution
leads $t\to te^{\pm iS_{A}\varphi_{A}}$, where $\varphi_{A}$ is the spin
magnetic flux corresponding to $U(1)_{A}$ symmetry.

As representative cases, we considered (1) a pair of magnetic flux tubes with
strength $\Phi=\pm\pi/2$ and (2) a pair of pseudospin magnetic flux tubes with
strength $\varphi=\pm\pi/2$, in the model Hamiltonian in Eq. (11) presented
the main text. Note that flux tubes stretches along the $x$ direction and are
located at $(y,z)=(N_{y}/2+1/2,N_{z}/4+1/2)$ and $(N_{y}/2+1/2,3N_{z}/4+1/2)$
in both calculations. Here, $N_{y}=32$ and $N_{z}=32$ are the dimension of the
lattice along $y$ and $z$ directions. To compute the charge and spin
densities, $J^{0}(y,z)$ and $\mathcal{J}^{0}_{A}(y,z)$, we take a half filling
and sum the square of eigenstates weighted by charge and spin over filled
states. Then, these values are averaged for $N_{k}=21$ different $k_{x}$
values, which corresponds to $N_{x}=21$ system size along the $x$ direction.
[See Fig. 2 in the main text for the distribution of $J^{0}(y,z)$ and
$\mathcal{J}^{0}_{A}(y,z)$.]

Our theoretical analysis predicts $\int_{\rm
flux}dydz\mathcal{J}^{0}_{X}(y,z)=-\frac{Q_{0}}{4\pi^{2}}\Phi$ in the case
(1), and $\int_{\rm flux}dydzJ^{0}(y,z)=-\frac{Q_{0}}{4\pi^{2}}\varphi$ in the
case (2), respectively. Here, $Q_{0}=Q_{xX}=8\pi/3$ is the nonzero component
in $Q_{iA}$, and $\int_{\rm flux}$ means that the charge/pseudospin density is
summed over a restricted region close to one of flux tubes. When 16 sites near
a flux tube is considered, the numerical result is consistent with the
predicted value within 8.08% error. The deviation of numerical result from the
theoretical prediction becomes reduced if more sites are included. For
example, the deviation is 4.01% if 64 sites considered.

##  S6 Response to dislocations

In this section, we discuss the response to dislocations of the weak TI with
nonzero pseudospin-momentum quadrupole $Q_{iA}$.

Figure S2:  Chiral modes localized on the screw dislocations. Energy spectrum
of pseudospin sectors (a) $(+,+)$ and $(+,-)$ and (b) $(-,+)$ and $(-,-)$ of
the weak TI model $H_{3D}({\mathbf{k}},t=1)$. Note that in our numerical
calculation, we have a lattice that has $24\times 24$ sites in the $yz$-plane,
and has two screw dislocations with Burgers vectors $\mathbf{b}=\pm
a\mathbf{e}_{x}$ along the $x$-direction. The center of the screw dislocations
are located at $(y,z)=(12,12)$ and $(12,24)$. Modes highlighted by red (blue)
color are chiral modes localized on the dislocation with Burgers vector
$b=a\mathbf{e}_{x}$ ($b=-a\mathbf{e}_{x}$).

First, we consider a continuous deformation from the Weyl semimetal, which is
described by Eq. (11) in the main text, to a weak TI. To this end, we define a
one-parameter family of Hamiltonian,

$\displaystyle H_{3D}({\mathbf{k}},t)=(1-t)(\sin k_{x}\mu_{x}S_{Y}-\sin
k_{y}\mu_{x}S_{X})$ $\displaystyle+[m-1+(1-t)(\cos k_{x}+\cos k_{y})+\cos
k_{z}]\mu_{z}$ $\displaystyle+\sin k_{z}\mu_{y}+t\mu_{x}U_{1}+t\mu_{z}U_{2},$
(S69)

which connects the Weyl semimetallic phase and weak TI phase adiabatically.
The diagonal matrices $U_{1}$ and $U_{2}$ takes the form $U_{1}={\rm
Diag}(\sin(k_{x}-k_{y}),\,\sin(-k_{x}+k_{y}),\,\sin(-k_{x}-k_{y}),\,\sin(k_{x}+k_{y}))$
and $U_{2}={\rm
Diag}(\cos(k_{x}-k_{y}),\cos(-k_{x}+k_{y}),\cos(-k_{x}-k_{y}),\cos(k_{x}+k_{y}))$
under the basis where $S_{X}={\rm Diag}(1,-1,1,-1)$ and $S_{Y}={\rm
Diag}(1,-1,-1,1)$. Also, we set $m=-0.5$ throughout this section. By tuning
$t$ from $0$ to $1$, one can continuously move the Weyl points with opposite
chiralities toward ${\mathbf{k}}=(\pi,\pi,0)$, and then gap them out to derive
a weak TI.

Then, let us consider screw dislocations along the $x$-direction in the weak
TI model $H_{3D}({\mathbf{k}},t=1)$. We require periodic boundary conditions
along all three directions, and thus must have two screw dislocations that
have opposite Burgers vectors $\mathbf{b}=\pm a\mathbf{e}_{x}$, where $a$ and
$\mathbf{e}_{x}$ are the lattice constant and unit vector along the
$x$-direction respectively. From numerical calculations, we observe chiral and
anti-chiral modes trapped on the dislocations in different pseudospin sectors,
which is consistent with the cut and glue picture discussed in Ref. Ran  _et
al._ , 2009. As shown in Fig. S2(a), in pseudospin sector $(+1,+1)$ and
$(+1,-1)$, there is a chiral mode with positive (negative) chirality on the
screw dislocation with Burgers vector $\mathbf{b}=-a\mathbf{e}_{x}$
($\mathbf{b}=a\mathbf{e}_{x}$). In contrast, as shown in Fig. S2(b), in
pseudospin sector $(-1,+1)$ and $(-1,-1)$, there is a chiral mode with
negative (positive) chirality on the screw dislocation with Burgers vector
$\mathbf{b}=-a\mathbf{e}_{x}$ ($\mathbf{b}=a\mathbf{e}_{x}$). Therefore, a
single screw dislocation does not produce a charge anomaly, but can produce a
pseudospin anomaly, which is captured by the action in Eq. (19) of the main
text:

$\displaystyle S_{CS}$ $\displaystyle=\frac{e}{4\pi^{2}}\int
d^{4}x\epsilon^{\mu\nu\rho\sigma}Q_{lA}\mathfrak{e}^{l}_{\mu}a^{A}_{\nu}\partial_{\rho}A_{\sigma}.$
(S70)

By variations with respect to charge or pseudospin gauge fields, we can derive
the charge or pseudospin current as

$\displaystyle
J^{\mu}=\frac{e}{4\pi^{2}}\epsilon^{\mu\nu\rho\sigma}Q_{lA}\partial_{\rho}(\mathfrak{e}^{l}_{\sigma}a^{A}_{\nu}),$
(S71)
$\displaystyle\mathcal{J}^{\mu}_{A}=-\frac{e}{4\pi^{2}}\epsilon^{\mu\nu\rho\sigma}Q_{lA}\mathfrak{e}^{l}_{\nu}\partial_{\rho}A_{\sigma}.$
(S72)

It is straightforward to observe that $J^{\mu}$ is conserved and satisfies
$\partial_{\mu}J^{\mu}=0$. However, $\mathcal{J}^{\mu}_{A}$ has an anomalous
conservation law:

$\displaystyle\partial_{\mu}\mathcal{J}^{\mu}_{A}=-\frac{e}{4\pi^{2}}\epsilon^{\mu\nu\rho\sigma}Q_{lA}\partial_{\mu}\mathfrak{e}^{l}_{\nu}\partial_{\rho}A_{\sigma},$
(S73)

which says that if we have both nonzero torsional magnetic (electric) field
and electric (magnetic) field, then the pseudospin current is anomalous. In
lattice systems, a dislocation can provide nonzero torsional magnetic flux
(i.e., the flux of translation gauge field $\mathfrak{e}^{l}_{\nu}$). For
example, a screw dislocation with Burgers vector $\mathbf{b}=-a\mathbf{e}_{x}$
will lead to a nonzero torsional flux $\int
dydz(\partial_{y}\mathfrak{e}^{x}_{z}-\partial_{z}\mathfrak{e}^{x}_{y})=-a$
localized on the screw dislocation. With such a dislocation, if we furthermore
have a nonzero electric field along the $x$-direction, then the pseudospin
current will be anomalous and satisfies $\int
dydz\partial_{\mu}\mathcal{J}^{\mu}_{X}=-\frac{2e}{\pi}E_{x}$ (note that
$Q_{xX}=8\pi/a$). This is consistent with the chiral modes in different
pseudospin sectors that are bound to the screw dislocation.